Almkvist-Giullera formula for pi: Difference between revisions
Alextretyak (talk | contribs) (Added 11l) |
Alextretyak (talk | contribs) m (→{{header|11l}}) |
||
Line 32: | Line 32: | ||
=={{header|11l}}== |
=={{header|11l}}== |
||
{{trans| |
{{trans|C#}} |
||
<lang 11l>F isqrt(BigInt =x) |
<lang 11l>F isqrt(BigInt =x) |
Revision as of 23:56, 1 August 2021
You are encouraged to solve this task according to the task description, using any language you may know.
The Almkvist-Giullera formula for calculating 1/π2 is based on the Calabi-Yau differential equations of order 4 and 5, which were originally used to describe certain manifolds in string theory.
The formula is:
- 1/π2 = (25/3) ∑0∞ ((6n)! / (n!6))(532n2 + 126n + 9) / 10002n+1
This formula can be used to calculate the constant π-2, and thus to calculate π.
Note that, because the product of all terms but the power of 1000 can be calculated as an integer, the terms in the series can be separated into a large integer term:
- (25) (6n)! (532n2 + 126n + 9) / (3(n!)6) (***)
multiplied by a negative integer power of 10:
- 10-(6n + 3)
- Task
-
- Print the integer portions (the starred formula, which is without the power of 1000 divisor) of the first 10 terms of the series.
- Use the complete formula to calculate and print π to 70 decimal digits of precision.
11l
<lang 11l>F isqrt(BigInt =x)
BigInt q = 1 BigInt r = 0 BigInt t L q <= x q *= 4 L q > 1 q I/= 4 t = x - r - q r I/= 2 I t >= 0 x = t r += q R r
F dump(=digs, show)
V gb = 1 digs++ V dg = digs + gb BigInt t1 = 1 BigInt t2 = 9 BigInt t3 = 1 BigInt te BigInt su = 0 V t = BigInt(10) ^ (I dg <= 60 {0} E dg - 60) BigInt d = -1 BigInt _fn_ = 1
V n = 0 L n < dg I n > 0 t3 *= BigInt(n) ^ 6 te = t1 * t2 I/ t3 V z = dg - 1 - n * 6 I z > 0 te *= BigInt(10) ^ z E te I/= BigInt(10) ^ -z I show & n < 10 print(‘#2 #62’.format(n, te * 32 I/ 3 I/ t)) su += te
I te < 10 I show digs-- print("\n#. iterations required for #. digits after the decimal point.\n".format(n, digs)) L.break
L(j) n * 6 + 1 .. n * 6 + 6 t1 *= j d += 2 t2 += 126 + 532 * d
n++
V s = String(isqrt(BigInt(10) ^ (dg * 2 + 3) I/ su I/ 32 * 3 * BigInt(10) ^ (dg + 5))) R s[0]‘.’s[1 .+ digs]
print(dump(70, 1B))</lang>
- Output:
0 9600000000000000000000000000000000000000000000000000000000000 1 512256000000000000000000000000000000000000000000000000000000 2 19072247040000000000000000000000000000000000000000000000000 3 757482485760000000000000000000000000000000000000000000000 4 31254615037245600000000000000000000000000000000000000000 5 1320787470322549142065152000000000000000000000000000000 6 56727391979308908329225994240000000000000000000000000 7 2465060024817298714011276371558400000000000000000000 8 108065785435463945367040747443956640000000000000000 9 4770177939159496628747057049083997888000000000000 53 iterations required for 70 digits after the decimal point. 3.1415926535897932384626433832795028841971693993751058209749445923078164
C#
A little challenging due to lack of BigFloat or BigRational. Note the extended precision integers displayed for each term, not extended precision floats. Also features the next term based on the last term, rather than computing each term from scratch. And the multiply by 32, divide by 3 is reserved for final sum, instead of each term (except for the 0..9th displayed terms). <lang csharp>using System; using BI = System.Numerics.BigInteger; using static System.Console;
class Program {
static BI isqrt(BI x) { BI q = 1, r = 0, t; while (q <= x) q <<= 2; while (q > 1) { q >>= 2; t = x - r - q; r >>= 1; if (t >= 0) { x = t; r += q; } } return r; }
static string dump(int digs, bool show = false) {
int gb = 1, dg = ++digs + gb, z; BI t1 = 1, t2 = 9, t3 = 1, te, su = 0, t = BI.Pow(10, dg <= 60 ? 0 : dg - 60), d = -1, fn = 1; for (BI n = 0; n < dg; n++) { if (n > 0) t3 *= BI.Pow(n, 6); te = t1 * t2 / t3; if ((z = dg - 1 - (int)n * 6) > 0) te *= BI.Pow (10, z); else te /= BI.Pow (10, -z); if (show && n < 10) WriteLine("{0,2} {1,62}", n, te * 32 / 3 / t); su += te; if (te < 10) { if (show) WriteLine("\n{0} iterations required for {1} digits " + "after the decimal point.\n", n, --digs); break; } for (BI j = n * 6 + 1; j <= n * 6 + 6; j++) t1 *= j; t2 += 126 + 532 * (d += 2); } string s = string.Format("{0}", isqrt(BI.Pow(10, dg * 2 + 3) / su / 32 * 3 * BI.Pow((BI)10, dg + 5))); return s[0] + "." + s.Substring(1, digs); }
static void Main(string[] args) { WriteLine(dump(70, true)); }
}</lang>
- Output:
0 9600000000000000000000000000000000000000000000000000000000000 1 512256000000000000000000000000000000000000000000000000000000 2 19072247040000000000000000000000000000000000000000000000000 3 757482485760000000000000000000000000000000000000000000000 4 31254615037245600000000000000000000000000000000000000000 5 1320787470322549142065152000000000000000000000000000000 6 56727391979308908329225994240000000000000000000000000 7 2465060024817298714011276371558400000000000000000000 8 108065785435463945367040747443956640000000000000000 9 4770177939159496628747057049083997888000000000000 53 iterations required for 70 digits after the decimal point. 3.1415926535897932384626433832795028841971693993751058209749445923078164
Factor
<lang factor>USING: continuations formatting io kernel locals math math.factorials math.functions sequences ;
- integer-term ( n -- m )
32 6 n * factorial * 532 n sq * 126 n * + 9 + * n factorial 6 ^ 3 * / ;
- exponent-term ( n -- m ) 6 * 3 + neg ;
- nth-term ( n -- x )
[ integer-term ] [ exponent-term 10^ * ] bi ;
! Factor doesn't have an arbitrary-precision square root afaik, ! so make one using Heron's method.
- sqrt-approx ( r x -- r' x ) [ over / + 2 / ] keep ;
- almkvist-guillera ( precision -- x )
0 0 :> ( summed! next-add! ) [ 100,000,000 <iota> [| n | summed n nth-term + next-add! next-add summed - abs precision neg 10^ < [ return ] when next-add summed! ] each ] with-return next-add ;
CONSTANT: 1/pi 113/355 ! Use as initial guess for square root approximation
- pi ( -- )
1/pi 70 almkvist-guillera 5 [ sqrt-approx ] times drop recip "%.70f\n" printf ;
! Task "N Integer Portion Pow Nth Term (33 dp)" print 89 CHAR: - <repetition> print 10 [
dup [ integer-term ] [ exponent-term ] [ nth-term ] tri "%d %44d %3d %.33f\n" printf
] each-integer nl "Pi to 70 decimal places:" print pi</lang>
- Output:
N Integer Portion Pow Nth Term (33 dp) ----------------------------------------------------------------------------------------- 0 96 -3 0.096000000000000000000000000000000 1 5122560 -9 0.005122560000000000000000000000000 2 190722470400 -15 0.000190722470400000000000000000000 3 7574824857600000 -21 0.000007574824857600000000000000000 4 312546150372456000000 -27 0.000000312546150372456000000000000 5 13207874703225491420651520 -33 0.000000013207874703225491420651520 6 567273919793089083292259942400 -39 0.000000000567273919793089083292260 7 24650600248172987140112763715584000 -45 0.000000000024650600248172987140113 8 1080657854354639453670407474439566400000 -51 0.000000000001080657854354639453670 9 47701779391594966287470570490839978880000000 -57 0.000000000000047701779391594966287 Pi to 70 decimal places: 3.1415926535897932384626433832795028841971693993751058209749445923078164
Go
<lang go>package main
import (
"fmt" "math/big" "strings"
)
func factorial(n int64) *big.Int {
var z big.Int return z.MulRange(1, n)
}
var one = big.NewInt(1) var three = big.NewInt(3) var six = big.NewInt(6) var ten = big.NewInt(10) var seventy = big.NewInt(70)
func almkvistGiullera(n int64, print bool) *big.Rat {
t1 := big.NewInt(32) t1.Mul(factorial(6*n), t1) t2 := big.NewInt(532*n*n + 126*n + 9) t3 := new(big.Int) t3.Exp(factorial(n), six, nil) t3.Mul(t3, three) ip := new(big.Int) ip.Mul(t1, t2) ip.Quo(ip, t3) pw := 6*n + 3 t1.SetInt64(pw) tm := new(big.Rat).SetFrac(ip, t1.Exp(ten, t1, nil)) if print { fmt.Printf("%d %44d %3d %-35s\n", n, ip, -pw, tm.FloatString(33)) } return tm
}
func main() {
fmt.Println("N Integer Portion Pow Nth Term (33 dp)") fmt.Println(strings.Repeat("-", 89)) for n := int64(0); n < 10; n++ { almkvistGiullera(n, true) }
sum := new(big.Rat) prev := new(big.Rat) pow70 := new(big.Int).Exp(ten, seventy, nil) prec := new(big.Rat).SetFrac(one, pow70) n := int64(0) for { term := almkvistGiullera(n, false) sum.Add(sum, term) z := new(big.Rat).Sub(sum, prev) z.Abs(z) if z.Cmp(prec) < 0 { break } prev.Set(sum) n++ } sum.Inv(sum) pi := new(big.Float).SetPrec(256).SetRat(sum) pi.Sqrt(pi) fmt.Println("\nPi to 70 decimal places is:") fmt.Println(pi.Text('f', 70))
}</lang>
- Output:
N Integer Portion Pow Nth Term (33 dp) ----------------------------------------------------------------------------------------- 0 96 -3 0.096000000000000000000000000000000 1 5122560 -9 0.005122560000000000000000000000000 2 190722470400 -15 0.000190722470400000000000000000000 3 7574824857600000 -21 0.000007574824857600000000000000000 4 312546150372456000000 -27 0.000000312546150372456000000000000 5 13207874703225491420651520 -33 0.000000013207874703225491420651520 6 567273919793089083292259942400 -39 0.000000000567273919793089083292260 7 24650600248172987140112763715584000 -45 0.000000000024650600248172987140113 8 1080657854354639453670407474439566400000 -51 0.000000000001080657854354639453670 9 47701779391594966287470570490839978880000000 -57 0.000000000000047701779391594966287 Pi to 70 decimal places is: 3.1415926535897932384626433832795028841971693993751058209749445923078164
Julia
<lang julia>using Formatting
setprecision(BigFloat, 300)
function integerterm(n)
p = BigInt(532) * n * n + BigInt(126) * n + 9 return (p * BigInt(2)^5 * factorial(BigInt(6) * n)) ÷ (3 * factorial(BigInt(n))^6)
end
exponentterm(n) = -(6n + 3)
nthterm(n) = integerterm(n) * big"10.0"^exponentterm(n)
println(" N Integer Term Power of 10 Nth Term") println("-"^90) for n in 0:9
println(lpad(n, 3), lpad(integerterm(n), 48), lpad(exponentterm(n), 4), lpad(format("{1:22.19e}", nthterm(n)), 35))
end
function AlmkvistGuillera(floatprecision)
summed = nthterm(0) for n in 1:10000000 next = summed + nthterm(n) if abs(next - summed) < big"10.0"^(-floatprecision) return next end summed = next end
end
println("\nπ to 70 digits is ", format(big"1.0" / sqrt(AlmkvistGuillera(70)), precision=70))
println("Computer π is ", format(π + big"0.0", precision=70))
</lang>
- Output:
N Integer Term Power of 10 Nth Term ------------------------------------------------------------------------------------------ 0 96 -3 9.6000000000000000000e-02 1 5122560 -9 5.1225600000000000000e-03 2 190722470400 -15 1.9072247040000000000e-04 3 7574824857600000 -21 7.5748248576000000000e-06 4 312546150372456000000 -27 3.1254615037245600000e-07 5 13207874703225491420651520 -33 1.3207874703225491421e-08 6 567273919793089083292259942400 -39 5.6727391979308908329e-10 7 24650600248172987140112763715584000 -45 2.4650600248172987140e-11 8 1080657854354639453670407474439566400000 -51 1.0806578543546394537e-12 9 47701779391594966287470570490839978880000000 -57 4.7701779391594966287e-14 π to 70 digits is 3.1415926535897932384626433832795028841971693993751058209749445923078164 Computer π is 3.1415926535897932384626433832795028841971693993751058209749445923078164
Mathematica/Wolfram Language
<lang Mathematica>ClearAll[numerator, denominator] numerator[n_] := (2^5) ((6 n)!) (532 n^2 + 126 n + 9)/(3 (n!)^6) denominator[n_] := 10^(6 n + 3) numerator /@ Range[0, 9] val = 1/Sqrt[Total[numerator[#]/denominator[#] & /@ Range[0, 100]]]; N[val, 70]</lang>
- Output:
{96,5122560,190722470400,7574824857600000,312546150372456000000,13207874703225491420651520,567273919793089083292259942400,24650600248172987140112763715584000,1080657854354639453670407474439566400000,47701779391594966287470570490839978880000000} 3.141592653589793238462643383279502884197169399375105820974944592307816
Nim
Derived from Wren version with some simplifications. <lang Nim>import strformat, strutils import decimal
proc fact(n: int): DecimalType =
result = newDecimal(1) if n < 2: return for i in 2..n: result *= i
proc almkvistGiullera(n: int): DecimalType =
## Return the integer portion of the nth term of Almkvist-Giullera sequence. let t1 = fact(6 * n) * 32 let t2 = 532 * n * n + 126 * n + 9 let t3 = fact(n) ^ 6 * 3 result = t1 * t2 / t3
let One = newDecimal(1)
setPrec(78) echo "N Integer portion" echo repeat('-', 47) for n in 0..9:
echo &"{n} {almkvistGiullera(n):>44}"
echo()
echo "Pi to 70 decimal places:" var
sum = newDecimal(0) prev = newDecimal(0) prec = One.scaleb(newDecimal(-70)) n = 0
while true:
sum += almkvistGiullera(n) / One.scaleb(newDecimal(6 * n + 3)) if abs(sum - prev) < prec: break prev = sum.clone inc n
let pi = 1 / sqrt(sum) echo ($pi)[0..71]</lang>
- Output:
N Integer portion ----------------------------------------------- 0 96 1 5122560 2 190722470400 3 7574824857600000 4 312546150372456000000 5 13207874703225491420651520 6 567273919793089083292259942400 7 24650600248172987140112763715584000 8 1080657854354639453670407474439566400000 9 47701779391594966287470570490839978880000000 Pi to 70 decimal places: 3.1415926535897932384626433832795028841971693993751058209749445923078164
Perl
<lang perl>use strict; use warnings; use feature <say state>; use Math::AnyNum <:overload factorial>;
my $Precision = 70;
sub Integral {
my($n) = @_; int (2**5 * factorial(6*$n) * (532*$n**2 + 126*$n + 9)) / (3*factorial($n)**6)
}
sub A_G {
my($n) = @_; Integral($n) / (10**(6*$n + 3));
}
sub Pi {
my($n) = @_; local $Math::AnyNum::PREC = 5 * $Precision; state $AGcache = 0; $AGcache += A_G($n); Math::AnyNum::pow($AGcache,-0.5)->round(-$Precision)->stringify;
}
say 'First 10 integer portions: '; say "$_ " . Integral($_) for 0..9;
my $next = ; my $target = Pi(my $Nth = 0); while () {
last if ($next = Pi(++$Nth)) eq $target; $target = $next;
}
printf "π to $Precision decimal places is:\n$target";</lang>
- Output:
First 10 integer portions: 0 96 1 5122560 2 190722470400 3 7574824857600000 4 312546150372456000000 5 13207874703225491420651520 6 567273919793089083292259942400 7 24650600248172987140112763715584000 8 1080657854354639453670407474439566400000 9 47701779391594966287470570490839978880000000 π to 70 decimal places is: 3.1415926535897932384626433832795028841971693993751058209749445923078164
Phix
with javascript_semantics requires("1.0.0") include mpfr.e mpfr_set_default_precision(-70) function almkvistGiullera(integer n, bool bPrint) mpz {t1,t2,ip} = mpz_inits(3) mpz_fac_ui(t1,6*n) mpz_mul_si(t1,t1,32) -- t1:=2^5*(6n)! mpz_fac_ui(t2,n) mpz_pow_ui(t2,t2,6) mpz_mul_si(t2,t2,3) -- t2:=3*(n!)^6 mpz_mul_si(ip,t1,532*n*n+126*n+9) -- ip:=t1*(532n^2+126n+9) mpz_fdiv_q(ip,ip,t2) -- ip:=ip/t2 integer pw := 6*n+3 mpz_ui_pow_ui(t1,10,pw) -- t1 := 10^(6n+3) mpq tm = mpq_init_set_z(ip,t1) -- tm := rat(ip/t1) if bPrint then string ips = mpz_get_str(ip), tms = mpfr_get_fixed(mpfr_init_set_q(tm),50) tms = trim_tail(tms,"0") printf(1,"%d %44s %3d %s\n", {n, ips, -pw, tms}) end if return tm end function constant hdr = "N --------------- Integer portion ------------- Pow ----------------- Nth term (50 dp) -----------------" printf(1,"%s\n%s\n",{hdr,repeat('-',length(hdr))}) for n=0 to 9 do {} = almkvistGiullera(n, true) end for mpq {res,prev,z} = mpq_inits(3), prec = mpq_init_set_str(sprintf("1/1%s",repeat('0',70))) integer n = 0 while true do mpq term := almkvistGiullera(n, false) mpq_add(res,res,term) mpq_sub(z,res,prev) mpq_abs(z,z) if mpq_cmp(z,prec) < 0 then exit end if mpq_set(prev,res) n += 1 end while mpq_inv(res,res) mpfr pi = mpfr_init_set_q(res) mpfr_sqrt(pi,pi) printf(1,"\nCalculation of pi took %d iterations using the Almkvist-Giullera formula.\n\n",n) printf(1,"Pi to 70 d.p.: %s\n",mpfr_get_fixed(pi,70)) mpfr_const_pi(pi) printf(1,"Pi (builtin) : %s\n",mpfr_get_fixed(pi,70))
- Output:
N --------------- Integer portion ------------- Pow ----------------- Nth term (50 dp) ----------------- ---------------------------------------------------------------------------------------------------------- 0 96 -3 0.096 1 5122560 -9 0.00512256 2 190722470400 -15 0.0001907224704 3 7574824857600000 -21 0.0000075748248576 4 312546150372456000000 -27 0.000000312546150372456 5 13207874703225491420651520 -33 0.00000001320787470322549142065152 6 567273919793089083292259942400 -39 0.0000000005672739197930890832922599424 7 24650600248172987140112763715584000 -45 0.000000000024650600248172987140112763715584 8 1080657854354639453670407474439566400000 -51 0.0000000000010806578543546394536704074744395664 9 47701779391594966287470570490839978880000000 -57 0.00000000000004770177939159496628747057049083997888 Calculation of pi took 52 iterations using the Almkvist-Giullera formula. Pi to 70 d.p.: 3.1415926535897932384626433832795028841971693993751058209749445923078164 Pi (builtin) : 3.1415926535897932384626433832795028841971693993751058209749445923078164
Python
<lang python>import mpmath as mp
with mp.workdps(72):
def integer_term(n): p = 532 * n * n + 126 * n + 9 return (p * 2**5 * mp.factorial(6 * n)) / (3 * mp.factorial(n)**6)
def exponent_term(n): return -(mp.mpf("6.0") * n + 3)
def nthterm(n): return integer_term(n) * mp.mpf("10.0")**exponent_term(n)
for n in range(10): print("Term ", n, ' ', int(integer_term(n)))
def almkvist_guillera(floatprecision): summed, nextadd = mp.mpf('0.0'), mp.mpf('0.0') for n in range(100000000): nextadd = summed + nthterm(n) if abs(nextadd - summed) < 10.0**(-floatprecision): break
summed = nextadd
return nextadd
print('\nπ to 70 digits is ', end=) mp.nprint(mp.mpf(1.0 / mp.sqrt(almkvist_guillera(70))), 71) print('mpmath π is ', end=) mp.nprint(mp.pi, 71)
</lang>
- Output:
Term 0 96 Term 1 5122560 Term 2 190722470400 Term 3 7574824857600000 Term 4 312546150372456000000 Term 5 13207874703225491420651520 Term 6 567273919793089083292259942400 Term 7 24650600248172987140112763715584000 Term 8 1080657854354639453670407474439566400000 Term 9 47701779391594966287470570490839978880000000 π to 70 digits is 3.1415926535897932384626433832795028841971693993751058209749445923078164 mpmath π is 3.1415926535897932384626433832795028841971693993751058209749445923078164
Quackery
<lang Quackery> [ $ "bigrat.qky" loadfile ] now!
[ 1 swap times [ i^ 1+ * ] ] is ! ( n --> n )
[ dup dup 2 ** 532 * over 126 * + 9 + swap 6 * ! * 32 * swap ! 6 ** 3 * / ] is intterm ( n --> n )
[ dup intterm 10 rot 6 * 3 + ** reduce ] is vterm ( n --> n/d )
10 times [ i^ intterm echo cr ] cr 0 n->v 53 times [ i^ vterm v+ ] 1/v 70 vsqrt drop 70 point$ echo$ cr</lang>
- Output:
96 5122560 190722470400 7574824857600000 312546150372456000000 13207874703225491420651520 567273919793089083292259942400 24650600248172987140112763715584000 1080657854354639453670407474439566400000 47701779391594966287470570490839978880000000 3.1415926535897932384626433832795028841971693993751058209749445923078164
Raku
<lang perl6># 20201013 Raku programming solution
use BigRoot; use Rat::Precise; use experimental :cached;
BigRoot.precision = 75 ; my $Precision = 70 ; my $AGcache = 0 ;
sub postfix:<!>(Int $n --> Int) is cached { [*] 1 .. $n }
sub Integral(Int $n --> Int) is cached {
(2⁵*(6*$n)! * (532*$n² + 126*$n + 9)) div (3*($n!)⁶)
}
sub A-G(Int $n --> FatRat) is cached { # Almkvist-Giullera
Integral($n).FatRat / (10**(6*$n + 3)).FatRat
}
sub Pi(Int $n --> Str) {
(1/(BigRoot.newton's-sqrt: $AGcache += A-G $n)).precise($Precision)
}
say "First 10 integer portions : "; say $_, "\t", Integral $_ for ^10;
my $target = Pi my $Nth = 0;
loop { $target eq ( my $next = Pi ++$Nth ) ?? ( last ) !! $target = $next }
say "π to $Precision decimal places is :\n$target"</lang>
- Output:
First 10 integer portions : 0 96 1 5122560 2 190722470400 3 7574824857600000 4 312546150372456000000 5 13207874703225491420651520 6 567273919793089083292259942400 7 24650600248172987140112763715584000 8 1080657854354639453670407474439566400000 9 47701779391594966287470570490839978880000000 π to 70 decimal places is : 3.1415926535897932384626433832795028841971693993751058209749445923078164
REXX
<lang rexx>/*REXX program uses the Almkvist─Giullera formula for 1 / (pi**2) [or pi ** -2]. */ numeric digits length( pi() ) + length(.); w= 102 say $( , 3) $( , w%2) $('power', 5) $( , w) say $('N', 3) $('integer term', w%2) $('of 10', 5) $('Nth term', w) say $( , 3, "─") $( , w%2, "─") $( , 5, "─") $( , w, "─")
s= 0 /*initialize S (the sum) to zero. */ do n=0 until old=s; old= s /*use the "older" value of S for OLD.*/ a= 2**5 * !(6*n) * (532 * n**2 + 126*n + 9) / (3 * !(n)**6) z= 10 ** (- (6*n + 3) ) s= s + a * z if n>10 then do; do 3*(n==11); say ' .'; end; iterate; end say $(n, 3) right(a, w%2) $(powX(z), 5) right( lowE( format(a*z, 1, w-6, 2, 0)), w) end /*n*/
say say 'The calculation of pi took ' n " iterations with " digits() ,
" decimal digits precision using" subword( sourceLine(1), 4, 3).
say numeric digits length( pi() ) - length(.); d= digits() - length(.); @= ' ↓↓↓ ' say center(@ 'calculated pi to ' d " fractional decimal digits (below) is "@, d+4, '─') say ' 'sqrt(1/s); say say ' 'pi(); @= ' ↑↑↑ ' say center(@ 'the true pi to ' d " fractional decimal digits (above) is" @, d+4, '─') exit 0 /*stick a fork in it, we're all done. */ /*──────────────────────────────────────────────────────────────────────────────────────*/ $: procedure; parse arg text,width,fill; return center(text, width, left(fill, 1) ) !: procedure; parse arg x; !=1;; do j=2 to x; != !*j; end; return ! lowE: procedure; parse arg x; return translate(x, 'e', "E") powX: procedure; parse arg p; return right( format( p, 1, 3, 2, 0), 3) + 0 /*──────────────────────────────────────────────────────────────────────────────────────*/ pi: pi=3.141592653589793238462643383279502884197169399375105820974944592307816406286208,
||9986280348253421170679821480865132823066470938446095505822317253594081284811174503 return pi
/*──────────────────────────────────────────────────────────────────────────────────────*/ sqrt: procedure; parse arg x; if x=0 then return 0; d=digits(); numeric digits; h=d+6
m.=9; numeric form; parse value format(x,2,1,,0) 'E0' with g 'E' _ .; g=g*.5'e'_ % 2 do j=0 while h>9; m.j= h; h= h % 2 + 1; end /*j*/ do k=j+5 to 0 by -1; numeric digits m.k; g= (g + x/g) * .5; end /*k*/ numeric digits d; return g/1</lang>
- output when using the internal default input:
(Shown at two─thirds size.)
power N integer term of 10 Nth term ─── ─────────────────────────────────────────────────── ───── ────────────────────────────────────────────────────────────────────────────────────────────────────── 0 96 -3 9.600000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000e-02 1 5122560 -9 5.122560000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000e-03 2 190722470400 -15 1.907224704000000000000000000000000000000000000000000000000000000000000000000000000000000000000000e-04 3 7574824857600000 -21 7.574824857600000000000000000000000000000000000000000000000000000000000000000000000000000000000000e-06 4 312546150372456000000 -27 3.125461503724560000000000000000000000000000000000000000000000000000000000000000000000000000000000e-07 5 13207874703225491420651520 -33 1.320787470322549142065152000000000000000000000000000000000000000000000000000000000000000000000000e-08 6 567273919793089083292259942400 -39 5.672739197930890832922599424000000000000000000000000000000000000000000000000000000000000000000000e-10 7 24650600248172987140112763715584000 -45 2.465060024817298714011276371558400000000000000000000000000000000000000000000000000000000000000000e-11 8 1080657854354639453670407474439566400000 -51 1.080657854354639453670407474439566400000000000000000000000000000000000000000000000000000000000000e-12 9 47701779391594966287470570490839978880000000 -57 4.770177939159496628747057049083997888000000000000000000000000000000000000000000000000000000000000e-14 10 2117262852373157207626265529989139651218848358400 -63 2.117262852373157207626265529989139651218848358400000000000000000000000000000000000000000000000000e-15 . . . The calculation of pi took 122 iterations with 163 decimal digits precision using the Almkvist─Giullera formula. ────────────────────────────────────────────── ↓↓↓ calculated pi to 160 fractional decimal digits (below) is ↓↓↓ ─────────────────────────────────────────────── 3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821480865132823066470938446095505822317253594081284811174503 3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821480865132823066470938446095505822317253594081284811174503 ────────────────────────────────────────────── ↑↑↑ the true pi to 160 fractional decimal digits (above) is ↑↑↑ ───────────────────────────────────────────────
Visual Basic .NET
<lang vbnet>Imports System, BI = System.Numerics.BigInteger, System.Console
Module Module1
Function isqrt(ByVal x As BI) As BI Dim t As BI, q As BI = 1, r As BI = 0 While q <= x : q <<= 2 : End While While q > 1 : q >>= 2 : t = x - r - q : r >>= 1 If t >= 0 Then x = t : r += q End While : Return r End Function
Function dump(ByVal digs As Integer, ByVal Optional show As Boolean = False) As String digs += 1 Dim z As Integer, gb As Integer = 1, dg As Integer = digs + gb Dim te As BI, t1 As BI = 1, t2 As BI = 9, t3 As BI = 1, su As BI = 0, t As BI = BI.Pow(10, If(dg <= 60, 0, dg - 60)), d As BI = -1, fn As BI = 1 For n As BI = 0 To dg - 1 If n > 0 Then t3 = t3 * BI.Pow(n, 6) te = t1 * t2 / t3 : z = dg - 1 - CInt(n) * 6 If z > 0 Then te = te * BI.Pow(10, z) Else te = te / BI.Pow(10, -z) If show AndAlso n < 10 Then WriteLine("{0,2} {1,62}", n, te * 32 / 3 / t) su += te : If te < 10 Then digs -= 1 If show Then WriteLine(vbLf & "{0} iterations required for {1} digits " & _ "after the decimal point." & vbLf, n, digs) Exit For End If For j As BI = n * 6 + 1 To n * 6 + 6 t1 = t1 * j : Next d += 2 : t2 += 126 + 532 * d Next Dim s As String = String.Format("{0}", isqrt(BI.Pow(10, dg * 2 + 3) _ / su / 32 * 3 * BI.Pow(CType(10, BI), dg + 5))) Return s(0) & "." & s.Substring(1, digs) End Function
Sub Main(ByVal args As String()) WriteLine(dump(70, true)) End Sub
End Module</lang>
- Output:
0 9600000000000000000000000000000000000000000000000000000000000 1 512256000000000000000000000000000000000000000000000000000000 2 19072247040000000000000000000000000000000000000000000000000 3 757482485760000000000000000000000000000000000000000000000 4 31254615037245600000000000000000000000000000000000000000 5 1320787470322549142065152000000000000000000000000000000 6 56727391979308908329225994240000000000000000000000000 7 2465060024817298714011276371558400000000000000000000 8 108065785435463945367040747443956640000000000000000 9 4770177939159496628747057049083997888000000000000 53 iterations required for 70 digits after the decimal point. 3.1415926535897932384626433832795028841971693993751058209749445923078164
Wren
<lang ecmascript>import "/big" for BigInt, BigRat import "/fmt" for Fmt
var factorial = Fn.new { |n|
if (n < 2) return BigInt.one var fact = BigInt.one for (i in 2..n) fact = fact * i return fact
}
var almkvistGiullera = Fn.new { |n, print|
var t1 = factorial.call(6*n) * 32 var t2 = 532*n*n + 126*n + 9 var t3 = factorial.call(n).pow(6)*3 var ip = t1 * t2 / t3 var pw = 6*n + 3 var tm = BigRat.new(ip, BigInt.ten.pow(pw)) if (print) { Fmt.print("$d $44i $3d $-35s", n, ip, -pw, tm.toDecimal(33)) } else { return tm }
}
System.print("N Integer Portion Pow Nth Term (33 dp)") System.print("-" * 89) for (n in 0..9) {
almkvistGiullera.call(n, true)
}
var sum = BigRat.zero var prev = BigRat.zero var prec = BigRat.new(BigInt.one, BigInt.ten.pow(70)) var n = 0 while(true) {
var term = almkvistGiullera.call(n, false) sum = sum + term if ((sum-prev).abs < prec) break prev = sum n = n + 1
} var pi = BigRat.one/sum.sqrt(70) System.print("\nPi to 70 decimal places is:") System.print(pi.toDecimal(70))</lang>
- Output:
N Integer Portion Pow Nth Term (33 dp) ----------------------------------------------------------------------------------------- 0 96 -3 0.096 1 5122560 -9 0.00512256 2 190722470400 -15 0.0001907224704 3 7574824857600000 -21 0.0000075748248576 4 312546150372456000000 -27 0.000000312546150372456 5 13207874703225491420651520 -33 0.00000001320787470322549142065152 6 567273919793089083292259942400 -39 0.000000000567273919793089083292260 7 24650600248172987140112763715584000 -45 0.000000000024650600248172987140113 8 1080657854354639453670407474439566400000 -51 0.000000000001080657854354639453670 9 47701779391594966287470570490839978880000000 -57 0.000000000000047701779391594966287 Pi to 70 decimal places is: 3.1415926535897932384626433832795028841971693993751058209749445923078164