Almkvist-Giullera formula for pi
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
AArch64 Assembly
<lang AArch64 Assembly> /* ARM assembly AARCH64 Raspberry PI 3B */ /* program calculPi64.s */ /* this program use gmp library */ /* link with gcc option -lgmp */
/*******************************************/ /* Constantes file */ /*******************************************/ /* for this file see task include a file in language AArch64 assembly*/ .include "../includeConstantesARM64.inc"
.equ MAXI, 10 .equ SIZEBIG, 100
/*********************************/ /* Initialized data */ /*********************************/ .data szMessDebutPgm: .asciz "Program 64 bits start. \n" szMessPi: .asciz "\nPI = \n" szCarriageReturn: .asciz "\n"
szFormat: .asciz " %Zd\n" szFormatFloat: .asciz " %.*Ff\n" /*********************************/ /* UnInitialized data */ /*********************************/ .bss Result1: .skip SIZEBIG Result2: .skip SIZEBIG Result3: .skip SIZEBIG Result4: .skip SIZEBIG fIntex5: .skip SIZEBIG fIntex6: .skip SIZEBIG fIntex7: .skip SIZEBIG fSum: .skip SIZEBIG fSum1: .skip SIZEBIG sBuffer: .skip SIZEBIG fEpsilon: .skip SIZEBIG fPrec: .skip SIZEBIG fPI: .skip SIZEBIG fTEN: .skip SIZEBIG fONE: .skip SIZEBIG /*********************************/ /* code section */ /*********************************/ .text .global main main: // entry of program
ldr x0,qAdrszMessDebutPgm bl affichageMess mov x20,#0 // loop indice
1:
mov x0,x20 bl computeAlmkvist // compute mov x1,x0 ldr x0,qAdrszFormat // print big integer bl __gmp_printf add x20,x20,#1 cmp x20,#MAXI blt 1b // and loop mov x0,#560 // float précision in bits bl __gmpf_set_default_prec mov x19,#0 // compute indice ldr x0,qAdrfSum // init to zéro bl __gmpf_init ldr x0,qAdrfSum1 // init to zéro bl __gmpf_init ldr x0,qAdrfONE // result address mov x1,#1 // init à 1 bl __gmpf_init_set_ui ldr x0,qAdrfIntex5 // init to zéro bl __gmpf_init ldr x0,qAdrfIntex6 // init to zéro bl __gmpf_init ldr x0,qAdrfIntex7 // init to zéro bl __gmpf_init ldr x0,qAdrfEpsilon // init to zéro bl __gmpf_init ldr x0,qAdrfPrec // init to zéro bl __gmpf_init ldr x0,qAdrfPI // init to zéro bl __gmpf_init ldr x0,qAdrfTEN mov x1,#10 // init to 10 bl __gmpf_init_set_ui ldr x0,qAdrfIntex6 // compute 10 pow 70 ldr x1,qAdrfTEN mov x2,#70 bl __gmpf_pow_ui ldr x0,qAdrfEpsilon // divide 1 by 10 pow 70 ldr x1,qAdrfONE // dividende ldr x2,qAdrfIntex6 // divisor bl __gmpf_div
2: // PI compute loop
mov x0,x19
bl computeAlmkvist
mov x20,x0
mov x1,#6
mul x2,x1,x19
add x6,x2,#3 // compute 6n + 3
ldr x0,qAdrfIntex6 // compute 10 pow (6n+3)
ldr x1,qAdrfTEN
mov x2,x6
bl __gmpf_pow_ui
ldr x0,qAdrfIntex7 // compute 1 / 10 pow (6n+3)
ldr x1,qAdrfONE // dividende
ldr x2,qAdrfIntex6 // divisor
bl __gmpf_div
ldr x0,qAdrfIntex6 // result big float
mov x1,x20 // big integer Almkvist
bl __gmpf_set_z // conversion in big float
ldr x0,qAdrfIntex5 // result Almkvist * 1 / 10 pow (6n+3)
ldr x1,qAdrfIntex7 // operator 1
ldr x2,qAdrfIntex6 // operator 2
bl __gmpf_mul
ldr x0,qAdrfSum1 // terms addition
ldr x1,qAdrfSum
ldr x2,qAdrfIntex5
bl __gmpf_add
ldr x0,qAdrfSum // copy terms
ldr x1,qAdrfSum1
bl __gmpf_set
ldr x0,qAdrfIntex7 // compute 1 / sum
ldr x1,qAdrfONE // dividende
ldr x2,qAdrfSum // divisor
bl __gmpf_div
ldr x0,qAdrfPI // compute square root (1 / sum )
ldr x1,qAdrfIntex7
bl __gmpf_sqrt
ldr x0,qAdrfIntex6 // compute variation PI
ldr x1,qAdrfPrec
ldr x2,qAdrfPI
bl __gmpf_sub
ldr x0,qAdrfIntex6 // absolue value
ldr x1,qAdrfIntex5
bl __gmpf_abs
add x19,x19,#1 // increment indice
ldr x0,qAdrfPrec // copy PI -> prévious
ldr x1,qAdrfPI
bl __gmpf_set
ldr x0,qAdrfIntex6 // compare gap and epsilon
ldr x1,qAdrfEpsilon
bl __gmpf_cmp
cmp w0,#0 // !!! cmp return result on 32 bits
bgt 2b // if gap is highter -> loop
ldr x0,qAdrszMessPi // title display
bl affichageMess
ldr x2,qAdrfPI // PI display
ldr x0,qAdrszFormatFloat
mov x1,#70
bl __gmp_printf
100: // standard end of the program
mov x0, #0 // return code mov x8, #EXIT // request to exit program svc #0 // perform the system call
qAdrszMessDebutPgm: .quad szMessDebutPgm qAdrszCarriageReturn: .quad szCarriageReturn qAdrfIntex5: .quad fIntex5 qAdrfIntex6: .quad fIntex6 qAdrfIntex7: .quad fIntex7 qAdrfSum: .quad fSum qAdrfSum1: .quad fSum1 qAdrszFormatFloat: .quad szFormatFloat qAdrszMessPi: .quad szMessPi qAdrfEpsilon: .quad fEpsilon qAdrfPrec: .quad fPrec qAdrfPI: .quad fPI qAdrfTEN: .quad fTEN qAdrfONE: .quad fONE /***************************************************/ /* compute almkvist_giullera formula */ /***************************************************/ /* x0 contains the number */ computeAlmkvist:
stp x19,lr,[sp,-16]! // save registers mov x19,x0 mov x1,#6 mul x0,x1,x0 ldr x1,qAdrResult1 // result address bl computeFactorielle // compute (n*6)! mov x1,#532 mul x2,x19,x19 mul x2,x1,x2 mov x1,#126 mul x3,x19,x1 add x2,x2,x3 add x2,x2,#9 lsl x2,x2,#5 // * 32 ldr x0,qAdrResult2 // result ldr x1,qAdrResult1 // operator bl __gmpz_mul_ui mov x0,x19 ldr x1,qAdrResult1 bl computeFactorielle ldr x0,qAdrResult3 bl __gmpz_init // init to 0 ldr x0,qAdrResult3 // result ldr x1,qAdrResult1 // operator mov x2,#6 bl __gmpz_pow_ui ldr x0,qAdrResult1 // result ldr x1,qAdrResult3 // operator mov x2,#3 bl __gmpz_mul_ui ldr x0,qAdrResult3 // result ldr x1,qAdrResult2 // operator ldr x2,qAdrResult1 // operator bl __gmpz_cdiv_q ldr x0,qAdrResult3 // return result address
100:
ldp x19,lr,[sp],16 // restaur 2 registers ret // return to address lr x30
qAdrszFormat: .quad szFormat qAdrResult1: .quad Result1 qAdrResult2: .quad Result2 qAdrResult3: .quad Result3 /***************************************************/ /* compute factorielle N */ /***************************************************/ /* x0 contains the number */ /* x1 contains big number result address */ computeFactorielle:
stp x19,lr,[sp,-16]! // save registers stp x20,x21,[sp,-16]! // save registers mov x19,x0 // save N mov x20,x1 // save result address mov x0,x1 // result address mov x1,#1 // init to 1 bl __gmpz_init_set_ui ldr x0,qAdrResult4 bl __gmpz_init // init to 0 mov x21,#1
1: // loop
ldr x0,qAdrResult4 // result mov x1,x20 // operator 1 mov x2,x21 // operator 2 bl __gmpz_mul_ui mov x0,x20 // copy result4 -> result ldr x1,qAdrResult4 bl __gmpz_set add x21,x21,#1 // increment indice cmp x21,x19 // N ? ble 1b // no -> loop ldr x0,qAdrResult4 ldp x20,x21,[sp],16 // restaur 2 registers ldp x19,lr,[sp],16 // restaur 2 registers ret // return to address lr x30
qAdrResult4: .quad Result4 /********************************************************/ /* File Include fonctions */ /********************************************************/ /* for this file see task include a file in language AArch64 assembly */ .include "../includeARM64.inc" </lang>
Program 64 bits start. 96 5122560 190722470400 7574824857600000 312546150372456000000 13207874703225491420651520 567273919793089083292259942400 24650600248172987140112763715584000 1080657854354639453670407474439566400000 47701779391594966287470570490839978880000000 PI = 3.1415926535897932384626433832795028841971693993751058209749445923078164
ARM Assembly
<lang ARM Assembly> /* ARM assembly Raspberry PI */ /* program calculPi.s */ /* this program use gmp library package : libgmp3-dev */ /* link with gcc option -lgmp */
/* REMARK 1 : this program use routines in a include file see task Include a file language arm assembly for the routine affichageMess conversion10 see at end of this program the instruction include */
/* for constantes see task include a file in arm assembly */ /************************************/ /* Constantes */ /************************************/ .include "../constantes.inc"
.equ MAXI, 10 .equ SIZEBIG, 100
/*********************************/ /* Initialized data */ /*********************************/ .data szMessPi: .asciz "\nPI = \n" szCarriageReturn: .asciz "\n"
szFormat: .asciz " %Zd\n" szFormatFloat: .asciz " %.*Ff\n" /*********************************/ /* UnInitialized data */ /*********************************/ .bss Result1: .skip SIZEBIG Result2: .skip SIZEBIG Result3: .skip SIZEBIG Result4: .skip SIZEBIG fInter5: .skip SIZEBIG fInter6: .skip SIZEBIG fInter7: .skip SIZEBIG fSum: .skip SIZEBIG fSum1: .skip SIZEBIG sBuffer: .skip SIZEBIG fEpsilon: .skip SIZEBIG fPrec: .skip SIZEBIG fPI: .skip SIZEBIG fTEN: .skip SIZEBIG fONE: .skip SIZEBIG /*********************************/ /* code section */ /*********************************/ .text .global main main: @ entry of program
mov r4,#0 @ loop indice
1:
mov r0,r4 bl computeAlmkvist @ compute mov r1,r0 ldr r0,iAdrszFormat @ print big integer bl __gmp_printf add r4,r4,#1 cmp r4,#MAXI blt 1b @ and loop mov r0,#560 @ float précision in bits bl __gmpf_set_default_prec mov r4,#0 @ compute indice ldr r0,iAdrfSum @ init to zéro bl __gmpf_init ldr r0,iAdrfSum1 @ init to zéro bl __gmpf_init ldr r0,iAdrfONE @ result address mov r1,#1 @ init à 1 bl __gmpf_init_set_ui ldr r0,iAdrfInter5 @ init to zéro bl __gmpf_init ldr r0,iAdrfInter6 @ init to zéro bl __gmpf_init ldr r0,iAdrfInter7 @ init to zéro bl __gmpf_init ldr r0,iAdrfEpsilon @ init to zéro bl __gmpf_init ldr r0,iAdrfPrec @ init to zéro bl __gmpf_init ldr r0,iAdrfPI @ init to zéro bl __gmpf_init ldr r0,iAdrfTEN mov r1,#10 @ init to 10 bl __gmpf_init_set_ui ldr r0,iAdrfInter6 @ compute 10 pow 70 ldr r1,iAdrfTEN mov r2,#70 bl __gmpf_pow_ui ldr r0,iAdrfEpsilon @ divide 1 by 10 pow 70 ldr r1,iAdrfONE @ dividende ldr r2,iAdrfInter6 @ divisor bl __gmpf_div
2: @ PI compute loop
mov r0,r4
bl computeAlmkvist
mov r5,r0
mov r1,#6
mul r2,r1,r4
add r6,r2,#3 @ compute 6n + 3
ldr r0,iAdrfInter6 @ compute 10 pow (6n+3)
ldr r1,iAdrfTEN
mov r2,r6
bl __gmpf_pow_ui
ldr r0,iAdrfInter7 @ compute 1 / 10 pow (6n+3)
ldr r1,iAdrfONE @ dividende
ldr r2,iAdrfInter6 @ divisor
bl __gmpf_div
ldr r0,iAdrfInter6 @ result big float
mov r1,r5 @ big integer Almkvist
bl __gmpf_set_z @ conversion in big float
ldr r0,iAdrfInter5 @ result Almkvist * 1 / 10 pow (6n+3)
ldr r1,iAdrfInter7 @ operator 1
ldr r2,iAdrfInter6 @ operator 2
bl __gmpf_mul
ldr r0,iAdrfSum1 @ terms addition
ldr r1,iAdrfSum
ldr r2,iAdrfInter5
bl __gmpf_add
ldr r0,iAdrfSum @ copy terms
ldr r1,iAdrfSum1
bl __gmpf_set
ldr r0,iAdrfInter7 @ compute 1 / sum
ldr r1,iAdrfONE @ dividende
ldr r2,iAdrfSum @ divisor
bl __gmpf_div
ldr r0,iAdrfPI @ compute square root (1 / sum )
ldr r1,iAdrfInter7
bl __gmpf_sqrt
ldr r0,iAdrfInter6 @ compute variation PI
ldr r1,iAdrfPrec
ldr r2,iAdrfPI
bl __gmpf_sub
ldr r0,iAdrfInter6 @ absolue value
ldr r1,iAdrfInter5
bl __gmpf_abs
add r4,r4,#1 @ increment indice
ldr r0,iAdrfPrec @ copy PI -> prévious
ldr r1,iAdrfPI
bl __gmpf_set
ldr r0,iAdrfInter6 @ compare gap and epsilon
ldr r1,iAdrfEpsilon
bl __gmpf_cmp
cmp r0,#0
bgt 2b @ if gap is highter -> loop
ldr r0,iAdrszMessPi @ title display
bl affichageMess
ldr r2,iAdrfPI @ PI display
ldr r0,iAdrszFormatFloat
mov r1,#70
bl __gmp_printf
100: @ standard end of the program
mov r0, #0 @ return code mov r7, #EXIT @ request to exit program svc #0 @ perform the system call
iAdrszCarriageReturn: .int szCarriageReturn iAdrfInter5: .int fInter5 iAdrfInter6: .int fInter6 iAdrfInter7: .int fInter7 iAdrfSum: .int fSum iAdrfSum1: .int fSum1 iAdrszFormatFloat: .int szFormatFloat iAdrszMessPi: .int szMessPi iAdrfEpsilon: .int fEpsilon iAdrfPrec: .int fPrec iAdrfPI: .int fPI iAdrfTEN: .int fTEN iAdrfONE: .int fONE /***************************************************/ /* compute almkvist_giullera formula */ /***************************************************/ /* r0 contains the number */ computeAlmkvist:
push {r1-r4,lr} @ save registers
mov r4,r0
mov r1,#6
mul r0,r1,r0
ldr r1,iAdrResult1 @ result address
bl computeFactorielle @ compute (n*6)!
mov r1,#532
mul r2,r4,r4
mul r2,r1,r2
mov r1,#126
mul r3,r4,r1
add r2,r2,r3
add r2,#9
lsl r2,r2,#5 @ * 32
ldr r0,iAdrResult2 @ result
ldr r1,iAdrResult1 @ operator
bl __gmpz_mul_ui
mov r0,r4
ldr r1,iAdrResult1
bl computeFactorielle
ldr r0,iAdrResult3
bl __gmpz_init @ init to 0
ldr r0,iAdrResult3 @ result
ldr r1,iAdrResult1 @ operator
mov r2,#6
bl __gmpz_pow_ui
ldr r0,iAdrResult1 @ result
ldr r1,iAdrResult3 @ operator
mov r2,#3
bl __gmpz_mul_ui
ldr r0,iAdrResult3 @ result
ldr r1,iAdrResult2 @ operator
ldr r2,iAdrResult1 @ operator
bl __gmpz_cdiv_q
ldr r0,iAdrResult3 @ return result address
pop {r1-r4,pc} @ restaur des registres
iAdrszFormat: .int szFormat iAdrResult1: .int Result1 iAdrResult2: .int Result2 iAdrResult3: .int Result3 /***************************************************/ /* compute factorielle N */ /***************************************************/ /* r0 contains the number */ /* r1 contains big number result address */ computeFactorielle:
push {r1-r6,lr} @ save registers
mov r5,r0 @ save N
mov r6,r1 @ save result address
mov r0,r1 @ result address
mov r1,#1 @ init to 1
bl __gmpz_init_set_ui
ldr r0,iAdrResult4
bl __gmpz_init @ init to 0
mov r4,#1
1: @ loop
ldr r0,iAdrResult4 @ result
mov r1,r6 @ operator 1
mov r2,r4 @ operator 2
bl __gmpz_mul_ui
mov r0,r6 @ copy result4 -> result
ldr r1,iAdrResult4
bl __gmpz_set
add r4,r4,#1 @ increment indice
cmp r4,r5 @ N ?
ble 1b @ no -> loop
mov r0,r1
pop {r1-r6,pc} @ restaur des registres
iAdrResult4: .int Result4
/***************************************************/ /* ROUTINES INCLUDE */ /***************************************************/ .include "../affichage.inc" </lang>
96 5122560 190722470400 7574824857600000 312546150372456000000 13207874703225491420651520 567273919793089083292259942400 24650600248172987140112763715584000 1080657854354639453670407474439566400000 47701779391594966287470570490839978880000000 PI = 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
C++
<lang cpp>#include <boost/multiprecision/cpp_dec_float.hpp>
- include <boost/multiprecision/gmp.hpp>
- include <iomanip>
- include <iostream>
namespace mp = boost::multiprecision; using big_int = mp::mpz_int; using big_float = mp::cpp_dec_float_100; using rational = mp::mpq_rational;
big_int factorial(int n) {
big_int result = 1;
for (int i = 2; i <= n; ++i)
result *= i;
return result;
}
// Return the integer portion of the nth term of Almkvist-Giullera sequence. big_int almkvist_giullera(int n) {
return factorial(6 * n) * 32 * (532 * n * n + 126 * n + 9) /
(pow(factorial(n), 6) * 3);
}
int main() {
std::cout << "n | Integer portion of nth term\n"
<< "------------------------------------------------\n";
for (int n = 0; n < 10; ++n)
std::cout << n << " | " << std::setw(44) << almkvist_giullera(n)
<< '\n';
big_float epsilon(pow(big_float(10), -70));
big_float prev = 0, pi = 0;
rational sum = 0;
for (int n = 0;; ++n) {
rational term(almkvist_giullera(n), pow(big_int(10), 6 * n + 3));
sum += term;
pi = sqrt(big_float(1 / sum));
if (abs(pi - prev) < epsilon)
break;
prev = pi;
}
std::cout << "\nPi to 70 decimal places is:\n"
<< std::fixed << std::setprecision(70) << pi << '\n';
}</lang>
- Output:
n | Integer portion of nth term ------------------------------------------------ 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 is: 3.1415926535897932384626433832795028841971693993751058209749445923078164
Common Lisp
<lang lisp>(ql:quickload :computable-reals :silent t) (use-package :computable-reals) (setq *print-prec* 70) (defparameter *iterations* 52)
- factorial using computable-reals multiplication op to keep precision
(defun !r (n)
(let ((p 1)) (loop for i from 2 to n doing (setq p (*r p i))) p))
- the nth integer term
(defun integral (n)
(let* ((polynomial (+r (*r 532 n n) (*r 126 n) 9))
(numer (*r 32 (!r (*r 6 n)) polynomial))
(denom (*r 3 (expt-r (!r n) 6))))
(/r numer denom)))
- the exponent for 10 in the nth term of the series
(defun power (n) (- 3 (* 6 (1+ n))))
- the nth term of the series
(defun almkvist-giullera (n)
(/r (integral n) (expt-r 10 (abs (power n)))))
- the sum of the first n terms
(defun almkvist-giullera-sigma (n)
(let ((s 0)) (loop for i from 0 to n doing (setq s (+r s (almkvist-giullera i)))) s))
- the approximation to pi after n terms
(defun almkvist-giullera-pi (n)
(sqrt-r (/r 1 (almkvist-giullera-sigma n))))
(format t "~A. ~44A~4A ~A~%" "N" "Integral part of Nth term" "×10^" "=Actual value of Nth term") (loop for i from 0 to 9 doing
(format t "~&~a. ~44d ~3d " i (integral i) (power i)) (finish-output *standard-output*) (print-r (almkvist-giullera i) 50 nil))
(format t "~%~%Pi after ~a iterations: " *iterations*) (print-r (almkvist-giullera-pi *iterations*) *print-prec*) </lang>
- Output:
N. Integral part of Nth term ×10^ =Actual value of Nth term 0. 96 -3 +0.09600000000000000000000000000000000000000000000000... 1. 5122560 -9 +0.00512256000000000000000000000000000000000000000000... 2. 190722470400 -15 +0.00019072247040000000000000000000000000000000000000... 3. 7574824857600000 -21 +0.00000757482485760000000000000000000000000000000000... 4. 312546150372456000000 -27 +0.00000031254615037245600000000000000000000000000000... 5. 13207874703225491420651520 -33 +0.00000001320787470322549142065152000000000000000000... 6. 567273919793089083292259942400 -39 +0.00000000056727391979308908329225994240000000000000... 7. 24650600248172987140112763715584000 -45 +0.00000000002465060024817298714011276371558400000000... 8. 1080657854354639453670407474439566400000 -51 +0.00000000000108065785435463945367040747443956640000... 9. 47701779391594966287470570490839978880000000 -57 +0.00000000000004770177939159496628747057049083997888... Pi after 52 iterations: +3.1415926535897932384626433832795028841971693993751058209749445923078164...
dc
<lang dc>[* factorial *]sz [ 1 Sp [ d lp * sp 1 - d 1 <f ]Sf d 1 <f Lfsz sz Lp ]sF
[* nth integral term *]sz [ sn 32 6 ln * lFx 532 ln * ln * 126 ln * + 9 + * * 3 ln lFx 6 ^ * / ]sI
[* nth exponent of 10 *]sz [ 1 + 6 * 3 r - ]sE
[* nth term in series *]sz [ d lIx r 10 r lEx _1 * ^ / ]sA
[* sum of the first n terms *]sz [ [li lAx ls + ss li 1 - d si 0 r !<L]sL si 0ss lLx ls]sS
[* approximation of pi after n terms *]sz [ lSx 1 r / v ]sP
[* count digits in a number *]sz [sn 0 sd lCx ld]sD [ld 1 + sd ln 10 0k / d sn 0 !=C]sC
[* print a number in a given column width *]sz [sw d lDx si lw li <T n]sW [[ ]n li 1 + si lw li <T]sT
[* main loop: print values for first 10 terms *]sz [N. Integral part of Nth term .................. × 10^ =Actual value of Nth term]p 0 sj [
lj n [. ]n lj lIx 0k 1 / 44 lWx [ ]n lj lEx 4 lWx [ ]n lj 99k lAx 50k 1 / p lj 1 + d sj 10 >M
] sM lMx
[]p
[* print resulting value of pi to 70 places *]sz [Pi after ]n 52n [ iterations:]p 99k 52 lPx 70k 1 / p</lang>
- Output:
N. Integral part of Nth term .................. × 10^ =Actual value of Nth term 0. 96 -3 .09600000000000000000000000000000000000000000000000 1. 5122560 -9 .00512256000000000000000000000000000000000000000000 2. 190722470400 -15 .00019072247040000000000000000000000000000000000000 3. 7574824857600000 -21 .00000757482485760000000000000000000000000000000000 4. 312546150372456000000 -27 .00000031254615037245600000000000000000000000000000 5. 13207874703225491420651520 -33 .00000001320787470322549142065152000000000000000000 6. 567273919793089083292259942400 -39 .00000000056727391979308908329225994240000000000000 7. 24650600248172987140112763715584000 -45 .00000000002465060024817298714011276371558400000000 8. 1080657854354639453670407474439566400000 -51 .00000000000108065785435463945367040747443956640000 9. 47701779391594966287470570490839978880000000 -57 .00000000000004770177939159496628747057049083997888 3.1415926535897932384626433832795028841971693993751058209749445923078\ 164
Erlang
This version uses integer math only (does not resort to a rational number package) Since the denominator is always a power of 10, it's possible to just keep track of the log of the denominator and scale the numerator accordingly; to keep track of the accuracy we get the order of magnitude of the difference between terms by subtracting the log of the numerator from the log of the denominator, so again, no rational arithmetic is needed.
However, Erlang does not have much in the way of calculating with large integers beyond basic arithmetic, so this version implements integer powers, logs, square roots, as well as the factorial function. <lang Erlang> -mode(compile).
% Integer math routines: factorial, power, square root, integer logarithm. % fac(N) -> fac(N, 1). fac(N, A) when N < 2 -> A; fac(N, A) -> fac(N - 1, N*A).
pow(_, N) when N < 0 -> pow_domain_error;
pow(2, N) -> 1 bsl N;
pow(A, N) -> ipow(A, N).
ipow(_, 0) -> 1; ipow(A, 1) -> A; ipow(A, 2) -> A*A; ipow(A, N) ->
case N band 1 of
0 -> X = ipow(A, N bsr 1), X*X;
1 -> A * ipow(A, N - 1)
end.
% integer logarithm, based on Zeckendorf representations of integers. % https://www.keithschwarz.com/interesting/code/?dir=zeckendorf-logarithm % we need this, since the exponents get larger than IEEE-754 double can handle.
lognext({A, B, S, T}) -> {B, A+B, T, S*T}. logprev({A, B, S, T}) -> {B-A, A, T div S, S}.
ilog(A, B) when (A =< 0) or (B < 2) -> ilog_domain_error; ilog(A, B) ->
UBound = bracket(A, {0, 1, 1, B}),
backlog(A, UBound, 0).
bracket(A, State = {_, _, _, T}) when T > A -> State; bracket(A, State) -> bracket(A, lognext(State)).
backlog(_, {0, _, 1, _}, Log) -> Log; backlog(N, State = {A, _, S, _}, Log) when S =< N ->
backlog(N div S, logprev(State), Log + A);
backlog(N, State, Log) -> backlog(N, logprev(State), Log).
isqrt(N) when N < 0 -> isqrt_domain_error;
isqrt(N) ->
X0 = pow(2, ilog(N, 2) div 2), iterate(N, newton(X0, N), N).
iterate(A, B, _) when A =< B -> A; iterate(_, B, N) -> iterate(B, newton(B, N), N).
newton(X, N) -> (X + N div X) div 2.
% With this out of the way, we can get down to some serious calculation.
%
term(N) -> { % returns numerator and log10 of the denominator.
(fac(6*N)*(N*(532*N + 126) + 9) bsl 5) div (3*pow(fac(N), 6)), 6*N + 3 }.
neg_term({N, D}) -> {-N, D}. abs_term({N, D}) -> {abs(N), D}.
add_term(T1 = {_, D1}, T2 = {_, D2}) when D1 > D2 -> add_term(T2, T1); add_term({N1, D1}, {N2, D2}) ->
Scale = pow(10, D2 - D1),
{N1*Scale + N2, D2}.
calculate(Prec) -> calculate(Prec, {0, 0}, 0). calculate(Prec, T0, K) ->
T1 = add_term(T0, term(K)),
{N, D} = abs_term(add_term(neg_term(T1), T0)),
Accuracy = D - ilog(N, 10),
if
Accuracy < Prec -> calculate(Prec, T1, K + 1);
true -> T1
end.
get_pi(Prec) ->
{N0, D0} = calculate(Prec),
% from the term, t = n0/10^d0, calculate 1/√t
% if the denominator is an odd power of 10, add 1 to the denominator and multiply the numerator by 10.
{N, D} = case D0 band 1 of
0 -> {N0, D0};
1 -> {10*N0, D0 + 1}
end,
[Three|Rest] = lists:sublist(
integer_to_list(pow(10, D) div isqrt(N)), Prec),
[Three, $. | Rest].
show_term({A, Decimals}) ->
Str = integer_to_list(A), [$0, $.] ++ lists:duplicate(Decimals - length(Str), $0) ++ Str.
main(_) ->
Terms = [term(N) || N <- lists:seq(0, 9)],
io:format("The first 10 terms as scaled decimals are:~n"),
[io:format(" ~s~n", [show_term(T)]) || T <- Terms],
io:format("~nThe sum of these terms (pi^-2) is ~s~n",
[show_term(lists:foldl(fun add_term/2, {0, 0}, Terms))]),
Pi70 = get_pi(71),
io:format("~npi to 70 decimal places:~n"),
io:format("~s~n", [Pi70]).
</lang>
- Output:
The first 10 terms as scaled decimals are:
0.096
0.005122560
0.000190722470400
0.000007574824857600000
0.000000312546150372456000000
0.000000013207874703225491420651520
0.000000000567273919793089083292259942400
0.000000000024650600248172987140112763715584000
0.000000000001080657854354639453670407474439566400000
0.000000000000047701779391594966287470570490839978880000000
The sum of these terms (pi^-2) is 0.101321183642335555356499725503850584160514406378880000000
pi to 70 decimal places:
3.1415926535897932384626433832795028841971693993751058209749445923078164
F#
This task uses Isqrt_(integer_square_root)_of_X#F.23 <lang fsharp> // Almkvist-Giullera formula for pi. Nigel Galloway: August 17th., 2021 let factorial(n:bigint)=MathNet.Numerics.SpecialFunctions.Factorial n let fN g=(532I*g*g+126I*g+9I)*(factorial(6I*g))/(3I*(factorial g)**6) [0..9]|>Seq.iter(bigint>>fN>>(*)32I>>printfn "%A\n") let _,n=Seq.unfold(fun(n,g)->let n=n*(10I**6)+fN g in Some(Isqrt((10I**(145+6*(int g)))/(32I*n)),(n,g+1I)))(0I,0I)|>Seq.pairwise|>Seq.find(fun(n,g)->n=g) printfn $"""pi to 70 decimal places is %s{(n.ToString()).Insert(1,".")}""" </lang>
- Output:
96 5122560 190722470400 7574824857600000 312546150372456000000 13207874703225491420651520 567273919793089083292259942400 24650600248172987140112763715584000 1080657854354639453670407474439566400000 47701779391594966287470570490839978880000000 pi to 70 decimal places is 3.14159265358979323846264338327950288419716939937510582097494459230781640
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
Haskell
<lang haskell>import Control.Monad import Data.Number.CReal import GHC.Integer import Text.Printf
iterations = 52 main = do
printf "N. %44s %4s %s\n"
"Integral part of Nth term" "×10^" "=Actual value of Nth term"
forM_ [0..9] $ \n ->
printf "%d. %44d %4d %s\n" n
(almkvistGiulleraIntegral n)
(tenExponent n)
(showCReal 50 (almkvistGiullera n))
printf "\nPi after %d iterations:\n" iterations putStrLn $ showCReal 70 $ almkvistGiulleraPi iterations
-- The integral part of the Nth term in the Almkvist-Giullera series almkvistGiulleraIntegral n =
let polynomial = (532 `timesInteger` n `timesInteger` n) `plusInteger` (126 `timesInteger` n) `plusInteger` 9
numerator = 32 `timesInteger` (facInteger (6 `timesInteger` n)) `timesInteger` polynomial
denominator = 3 `timesInteger` (powInteger (facInteger n) 6)
in numerator `divInteger` denominator
-- The exponent for 10 in the Nth term of the series tenExponent n = 3 `minusInteger` (6 `timesInteger` (1 `plusInteger` n))
-- The Nth term in the series (integral * 10^tenExponent) almkvistGiullera n = fromInteger (almkvistGiulleraIntegral n) / fromInteger (powInteger 10 (abs (tenExponent n)))
-- The sum of the first N terms almkvistGiulleraSum n = sum $ map almkvistGiullera [0 .. n]
-- The approximation of pi from the first N terms almkvistGiulleraPi n = sqrt $ 1 / almkvistGiulleraSum n
-- Utility: factorial for arbitrary-precision integers facInteger n = if n `leInteger` 1 then 1 else n `timesInteger` facInteger (n `minusInteger` 1)
-- Utility: exponentiation for arbitrary-precision integers powInteger 1 _ = 1 powInteger _ 0 = 1 powInteger b 1 = b powInteger b e = b `timesInteger` powInteger b (e `minusInteger` 1) </lang>
- Output:
N. Integral part of Nth term ×10^ =Actual value of Nth term 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 Pi after 52 iterations: 3.1415926535897932384626433832795028841971693993751058209749445923078164
J
This solution just has it hard-coded that 53 iterations is necessary for 70 decimals. It would be possible to write a loop with a test, though in practice it would also be acceptable to just experiment to find the number of iterations.
sqrt is noticeably slow, bringing execution time to over 1 second. I'm not sure if it's because it's coded imperatively using traditional loops vs. J point-free style, or if it's due to the fact that the numbers are very large. I suspect the latter since it only takes 4 iterations of Heron's method to get the square root. <lang J> numerator =: monad define "0
(3 * (! x: y)^6) %~ 32 * (!6x*y) * (y*(126 + 532*y)) + 9x
)
term =: numerator % 10x ^ 3 + 6&*
echo 'The first 10 numerators are:' echo ,. numerator i.10
echo echo 'The sum of the first 10 terms (pi^-2) is ', 0j15 ": +/ term i.10
heron =: [: -: ] + %
sqrt =: dyad define NB. usage: x0 tolerance sqrt x
NB. e.g.: (1, %10^100x) sqrt 2 -> √2 to 100 decimals as a ratio p/q
x0 =. }: x
eps =. }. x
x1 =. y heron x0
while. (| x1 - x0) > eps do.
x2 =. y heron x1
x0 =. x1
x1 =. x2
end.
x1
)
pi70 =. (355r113, %10^70x) sqrt % +/ term i.53 echo echo 'pi to 70 decimals: ', 0j70 ": pi70 exit </lang>
- Output:
The first 10 numerators are:
96
5122560
190722470400
7574824857600000
312546150372456000000
13207874703225491420651520
567273919793089083292259942400
24650600248172987140112763715584000
1080657854354639453670407474439566400000
47701779391594966287470570490839978880000000
The sum of the first 10 terms (pi^-2) is 0.101321183642336
pi to 70 decimals: 3.1415926535897932384626433832795028841971693993751058209749445923078164
JavaScript
to support use of module as main code
<lang javascript>import esMain from 'es-main'; import { BigFloat, set_precision as SetPrecision } from 'bigfloat-esnext';
const Iterations = 52;
export const demo = function() {
SetPrecision(-75);
console.log("N." + "Integral part of Nth term".padStart(45) + " ×10^ =Actual value of Nth term");
for (let i=0; i<10; i++) {
let line = `${i}. `;
line += `${integral(i)} `.padStart(45);
line += `${tenExponent(i)} `.padStart(5);
line += nthTerm(i);
console.log(line);
}
let pi = approximatePi(Iterations);
SetPrecision(-70);
pi = pi.dividedBy(100000).times(100000);
console.log(`\nPi after ${Iterations} iterations: ${pi}`)
}
export const bigFactorial = n => n <= 1n ? 1n : n * bigFactorial(n-1n);
// the nth integer term export const integral = function(i) {
let n = BigInt(i); const polynomial = 532n * n * n + 126n * n + 9n; const numerator = 32n * bigFactorial(6n * n) * polynomial; const denominator = 3n * bigFactorial(n) ** 6n; return numerator / denominator;
}
// the exponent for 10 in the nth term of the series export const tenExponent = n => 3n - 6n * (BigInt(n) + 1n);
// the nth term of the series export const nthTerm = n =>
new BigFloat(integral(n)).dividedBy(new BigFloat(10n ** -tenExponent(n)))
// the sum of the first n terms export const sumThrough = function(n) {
let sum = new BigFloat(0);
for (let i=0; i<=n; ++i) {
sum = sum.plus(nthTerm(i));
}
return sum;
}
// the approximation to pi after n terms export const approximatePi = n =>
new BigFloat(1).dividedBy(sumThrough(n)).sqrt();
if (esMain(import.meta))
demo();
</lang>
- Output:
N. Integral part of Nth term ×10^ =Actual value of Nth term 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 Pi after 52 iterations: 3.1415926535897932384626433832795028841971693993751058209749445923078164
jq
Adapted from Wren
Works with gojq, the Go implementation of jq
This entry uses the "rational" module, which can be found at Arithmetic/Rational#jq.
Preliminaries <lang jq># A reminder to include the "rational" module:
- include "rational";
def lpad($len): tostring | ($len - length) as $l | (" " * $l)[:$l] + .;
- To take advantage of gojq's arbitrary-precision integer arithmetic:
def power($b): . as $in | reduce range(0;$b) as $i (1; . * $in);
def factorial:
if . < 2 then 1 else reduce range(2;.+1) as $i (1; .*$i) end; </lang>
Almkvist-Giullera Formula <lang jq> def almkvistGiullera(print):
. as $n | ((6*$n) | factorial * 32) as $t1 | (532*$n*$n + 126*$n + 9) as $t2 | (($n | factorial | power(6))*3) as $t3 | ($t1 * $t2 / $t3) as $ip | ( 6*$n + 3) as $pw | r($ip; 10 | power($pw)) as $tm | if print then "\($n|lpad(2)) \($ip|lpad(44)) \(-$pw|lpad(3)), \($tm|r_to_decimal(100))" else $tm end; </lang>
The Tasks <lang jq> def task1:
"N Integer Portion Pow Nth Term",
("-" * 89),
(range(0;10) | almkvistGiullera(true)) ;
def task2($precision):
r(1; 10 | power($precision)) as $p
| {sum: r(0;1), prev: r(0;1), n: 0 }
| until(.stop;
.sum = radd(.sum; .n | almkvistGiullera(false))
| if rminus(.sum; .prev) | rabs | rlessthan($p)
then .stop = true
else .prev = .sum
| .n += 1
end)
| .sum | rinv
| rsqrt($precision)
| "\nPi to \($precision) decimal places is:",
"\(r_to_decimal($precision))" ;
task1, "" task2(70)</lang>
- Output:
N Integer Portion Pow Nth Term ----------------------------------------------------------------------------------------- 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 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 qw(say); use Math::AnyNum qw(:overload factorial);
sub almkvist_giullera_integral {
my($n) = @_; (32 * (14*$n * (38*$n + 9) + 9) * factorial(6*$n)) / (3*factorial($n)**6);
}
sub almkvist_giullera {
my($n) = @_; almkvist_giullera_integral($n) / (10**(6*$n + 3));
}
sub almkvist_giullera_pi {
my ($prec) = @_;
local $Math::AnyNum::PREC = 4*($prec+1);
my $sum = 0; my $target = ;
for (my $n = 0; ; ++$n) {
$sum += almkvist_giullera($n);
my $curr = ($sum**-.5)->as_dec;
return $target if ($curr eq $target);
$target = $curr;
}
}
say 'First 10 integer portions: '; say "$_ " . almkvist_giullera_integral($_) for 0..9;
my $precision = 70;
printf("π to %s decimal places is:\n%s\n",
$precision, almkvist_giullera_pi($precision));</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 ↑↑↑ ───────────────────────────────────────────────
Sidef
<lang ruby>func almkvist_giullera(n) {
(32 * (14*n * (38*n + 9) + 9) * (6*n)!) / (3 * n!**6)
}
func almkvist_giullera_pi(prec = 70) {
local Num!PREC = (4*(prec+1)).numify
var sum = 0 var target = -1
for n in (0..Inf) {
sum += (almkvist_giullera(n) / (10**(6*n + 3)))
var curr = (sum**-.5).as_dec
return target if (target == curr)
target = curr
}
}
say 'First 10 integer portions: '
10.of {|n|
say "#{n} #{almkvist_giullera(n)}"
}
with(70) {|n|
say "π to #{n} decimal places is:"
say almkvist_giullera_pi(n)
}</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
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
- Programming Tasks
- Solutions by Programming Task
- 11l
- AArch64 Assembly
- ARM Assembly
- C sharp
- System.Numerics
- C++
- Boost
- GMP
- Common Lisp
- Computable-reals
- Dc
- Erlang
- F Sharp
- Factor
- Go
- Haskell
- Numbers
- J
- JavaScript
- Bigfloat-esnext
- Es-main
- Jq
- Julia
- Mathematica
- Wolfram Language
- Nim
- Nim-decimal
- Perl
- Phix
- Python
- Quackery
- Raku
- REXX
- Sidef
- Visual Basic .NET
- Wren
- Wren-big
- Wren-fmt