Almkvist-Giullera formula for pi: Difference between revisions
Add Kotlin implementation
m (→{{header|Haskell}}: reorganize) |
(Add Kotlin implementation) |
||
| (26 intermediate revisions by 13 users not shown) | |||
Line 34:
{{trans|C#}}
<
BigInt q = 1
BigInt r = 0
Line 92:
R s[0]‘.’s[1 .+ digs]
print(dump(70, 1B))</
{{out}}
Line 111:
3.1415926535897932384626433832795028841971693993751058209749445923078164
</pre>
=={{header|AArch64 Assembly}}==
{{works with|as|Raspberry Pi 3B version Buster 64 bits <br> or android 64 bits with application Termux }}
<syntaxhighlight 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"
</syntaxhighlight>
<pre>
Program 64 bits start.
96
5122560
190722470400
7574824857600000
312546150372456000000
13207874703225491420651520
567273919793089083292259942400
24650600248172987140112763715584000
1080657854354639453670407474439566400000
47701779391594966287470570490839978880000000
PI =
3.1415926535897932384626433832795028841971693993751058209749445923078164
</pre>
=={{header|ALGOL 68}}==
{{works with|ALGOL 68G|Any - tested with release 2.8.3.win32}}
{{Trans|C++}}
Uses code from the [[Arithmetic/Rational#ALGOL_68|Arithmetic/Rational]] task.
<syntaxhighlight lang="algol68">
# Almkvist-Giullera formula for pi - translated from the C++ sample #
PR precision 1024 PR # set precision for LONG LONG modes #
MODE INTEGER = LONG LONG INT;
MODE FLOAT = LONG LONG REAL;
INTEGER zero = 0, one = 1, ten = 10;
# iterative Greatest Common Divisor routine, returns the gcd of m and n #
PROC gcd = ( INTEGER m, n )INTEGER:
BEGIN
INTEGER a := ABS m, b := ABS n;
WHILE b /= 0 DO
INTEGER new a = b;
b := a MOD b;
a := new a
OD;
a
END # gcd # ;
# code from the Arithmetic/Rational task modified to use LONG LONG INT #
MODE RATIONAL = STRUCT( INTEGER num #erator#, den #ominator# );
PROC lcm = ( INTEGER a, b )INTEGER: # least common multiple #
a OVER gcd(a, b) * b;
PRIO // = 9; # higher then the ** operator #
OP // = ( INTEGER num, den )RATIONAL: ( # initialise and normalise #
INTEGER common = gcd( num, den );
IF den < 0 THEN
( -num OVER common, -den OVER common )
ELSE
( num OVER common, den OVER common )
FI
);
OP + = (RATIONAL a, b)RATIONAL: (
INTEGER common = lcm( den OF a, den OF b );
RATIONAL result := ( common OVER den OF a * num OF a + common OVER den OF b * num OF b, common );
num OF result//den OF result
);
OP +:= = (REF RATIONAL a, RATIONAL b)REF RATIONAL: ( a := a + b );
# end code from the Arithmetic/Rational task modified to use LONG LONG INT #
OP / = ( FLOAT f, RATIONAL r )FLOAT: ( f * den OF r ) / num OF r;
INTEGER ag factorial n := 1;
INT ag last factorial := 0;
# returns factorial n, using ag factorial n and ag last factorial to reduce #
# the number of calculations #
PROC ag factorial = ( INT n )INTEGER:
BEGIN
IF n < ag last factorial THEN
ag last factorial := 0;
ag factorial n := 1
FI;
WHILE ag last factorial < n DO
ag factorial n *:= ( ag last factorial +:= 1 )
OD;
ag factorial n
END # ag gfgactorial # ;
# Return the integer portion of the nth term of Almkvist-Giullera sequence. #
PROC almkvist giullera = ( INT n )INTEGER:
ag factorial( 6 * n ) * 32 * ( 532 * n * n + 126 * n + 9 ) OVER ( ( ag factorial( n ) ^ 6 ) * 3 );
BEGIN
print( ( "n | Integer portion of nth term", newline ) );
print( ( "--+---------------------------------------------", newline ) );
FOR n FROM 0 TO 9 DO
print( ( whole( n, 0 ), " | ", whole( almkvist giullera( n ), -44 ), newline ) )
OD;
FLOAT epsilon = FLOAT(10) ^ -70;
FLOAT prev := 0, pi := 0;
RATIONAL sum := zero // one;
FOR n FROM 0
WHILE
RATIONAL term = almkvist giullera( n ) // ( ten ^ ( 6 * n + 3 ) );
sum +:= term;
pi := long long sqrt( FLOAT(1) / sum );
ABS ( pi - prev ) >= epsilon
DO
prev := pi
OD;
print( ( newline, "Pi to 70 decimal places is:", newline ) );
print( ( fixed( pi, -72, 70 ), newline ) )
END
</syntaxhighlight>
{{out}}
<pre>
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
</pre>
=={{header|ARM Assembly}}==
{{works with|as|Raspberry Pi <br> or android 32 bits with application Termux}}
<syntaxhighlight 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"
</syntaxhighlight>
<pre>
96
5122560
190722470400
7574824857600000
312546150372456000000
13207874703225491420651520
567273919793089083292259942400
24650600248172987140112763715584000
1080657854354639453670407474439566400000
47701779391594966287470570490839978880000000
PI =
3.1415926535897932384626433832795028841971693993751058209749445923078164
</pre>
=={{header|C sharp|C#}}==
{{libheader|System.Numerics}}
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).
<
using BI = System.Numerics.BigInteger;
using static System.Console;
Line 147 ⟶ 874:
static void Main(string[] args) {
WriteLine(dump(70, true)); }
}</
{{out}}
<pre> 0 9600000000000000000000000000000000000000000000000000000000000
Line 162 ⟶ 889:
53 iterations required for 70 digits after the decimal point.
3.1415926535897932384626433832795028841971693993751058209749445923078164
</pre>
=={{header|C++}}==
{{libheader|Boost}}
{{libheader|GMP}}
<syntaxhighlight 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';
}</syntaxhighlight>
{{out}}
<pre>
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
</pre>
Line 168 ⟶ 962:
{{libheader|computable-reals}}
{{trans|Raku}}
<
(use-package :computable-reals)
(setq *print-prec* 70)
Line 212 ⟶ 1,006:
(format t "~%~%Pi after ~a iterations: " *iterations*)
(print-r (almkvist-giullera-pi *iterations*) *print-prec*)
</syntaxhighlight>
{{Out}}
Line 233 ⟶ 1,027:
=={{header|dc}}==
{{trans|Common Lisp}}
<
[ 1 Sp [ d lp * sp 1 - d 1 <f ]Sf d 1 <f Lfsz sz Lp ]sF
Line 275 ⟶ 1,069:
[* print resulting value of pi to 70 places *]sz
[Pi after ]n 52n [ iterations:]p
99k 52 lPx 70k 1 / p</
{{Out}}
Line 298 ⟶ 1,092:
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.
<syntaxhighlight lang="erlang">
-mode(compile).
Line 403 ⟶ 1,197:
io:format("~npi to 70 decimal places:~n"),
io:format("~s~n", [Pi70]).
</syntaxhighlight>
{{Out}}
<pre>
Line 422 ⟶ 1,216:
pi to 70 decimal places:
3.1415926535897932384626433832795028841971693993751058209749445923078164
</pre>
=={{header|F_Sharp|F#}}==
This task uses [[Isqrt_(integer_square_root)_of_X#F.23]]
<syntaxhighlight 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,".")}"""
</syntaxhighlight>
{{out}}
<pre>
96
5122560
190722470400
7574824857600000
312546150372456000000
13207874703225491420651520
567273919793089083292259942400
24650600248172987140112763715584000
1080657854354639453670407474439566400000
47701779391594966287470570490839978880000000
pi to 70 decimal places is 3.14159265358979323846264338327950288419716939937510582097494459230781640
</pre>
=={{header|Factor}}==
{{works with|Factor|0.99 2020-08-14}}
<
math.factorials math.functions sequences ;
Line 468 ⟶ 1,288:
"%d %44d %3d %.33f\n" printf
] each-integer nl
"Pi to 70 decimal places:" print pi</
{{out}}
<pre>
Line 490 ⟶ 1,310:
=={{header|Go}}==
{{trans|Wren}}
<
import (
Line 556 ⟶ 1,376:
fmt.Println("\nPi to 70 decimal places is:")
fmt.Println(pi.Text('f', 70))
}</
{{out}}
Line 580 ⟶ 1,400:
{{libheader|numbers}}
{{trans|Common Lisp}}
<
import Data.Number.CReal
import GHC.Integer
Line 626 ⟶ 1,446:
powInteger b 1 = b
powInteger b e = b `timesInteger` powInteger b (e `minusInteger` 1)
</syntaxhighlight>
{{Out}}
Line 643 ⟶ 1,463:
Pi after 52 iterations:
3.1415926535897932384626433832795028841971693993751058209749445923078164</pre>
=={{header|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.
<syntaxhighlight 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 ''
</syntaxhighlight>
{{Out}}
<pre>
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
</pre>
=={{header|Java}}==
<syntaxhighlight lang="java">
import java.math.BigDecimal;
import java.math.BigInteger;
import java.math.MathContext;
import java.math.RoundingMode;
public final class AlmkvistGiulleraFormula {
public static void main(String[] aArgs) {
System.out.println("n Integer part");
System.out.println("================================================");
for ( int n = 0; n <= 9; n++ ) {
System.out.println(String.format("%d%47s", n, almkvistGiullera(n).toString()));
}
final int decimalPlaces = 70;
final MathContext mathContext = new MathContext(decimalPlaces + 1, RoundingMode.HALF_EVEN);
final BigDecimal epsilon = BigDecimal.ONE.divide(BigDecimal.TEN.pow(decimalPlaces));
BigDecimal previous = BigDecimal.ONE;
BigDecimal sum = BigDecimal.ZERO;
BigDecimal pi = BigDecimal.ZERO;
int n = 0;
while ( pi.subtract(previous).abs().compareTo(epsilon) >= 0 ) {
BigDecimal nextTerm = new BigDecimal(almkvistGiullera(n)).divide(BigDecimal.TEN.pow(6 * n + 3));
sum = sum.add(nextTerm);
previous = pi;
n += 1;
pi = BigDecimal.ONE.divide(sum, mathContext).sqrt(mathContext);
}
System.out.println(System.lineSeparator() + "pi to " + decimalPlaces + " decimal places:");
System.out.println(pi);
}
// The integer part of the n'th term of Almkvist-Giullera series.
private static BigInteger almkvistGiullera(int aN) {
BigInteger term1 = factorial(6 * aN).multiply(BigInteger.valueOf(32));
BigInteger term2 = BigInteger.valueOf(532 * aN * aN + 126 * aN + 9);
BigInteger term3 = factorial(aN).pow(6).multiply(BigInteger.valueOf(3));
return term1.multiply(term2).divide(term3);
}
private static BigInteger factorial(int aNumber) {
BigInteger result = BigInteger.ONE;
for ( int i = 2; i <= aNumber; i++ ) {
result = result.multiply(BigInteger.valueOf(i));
}
return result;
}
}
</syntaxhighlight>
{{ out }}
<pre>
n Integer part
================================================
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
</pre>
=={{header|JavaScript}}==
{{trans|Common Lisp}}
{{works with|Node.js|13+}}
{{libheader|bigfloat-esnext}}
{{libheader|es-main}} to support use of module as main code
<syntaxhighlight 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();
</syntaxhighlight>
{{Out}}
<pre>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</pre>
=={{header|jq}}==
'''Adapted from [[#Wren|Wren]]'''
'''Works with gojq, the Go implementation of jq'''
This entry uses the "rational" module, which can be found at [[Arithmetic/Rational#jq]].
'''Preliminaries'''
<syntaxhighlight 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; </syntaxhighlight>
'''Almkvist-Giullera Formula'''
<syntaxhighlight 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; </syntaxhighlight>
'''The Tasks'''
<syntaxhighlight 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)</syntaxhighlight>
{{out}}
<pre>
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
</pre>
=={{header|Julia}}==
<
setprecision(BigFloat, 300)
Line 679 ⟶ 1,785:
println("Computer π is ", format(π + big"0.0", precision=70))
</
<pre>
N Integer Term Power of 10 Nth Term
Line 696 ⟶ 1,802:
π to 70 digits is 3.1415926535897932384626433832795028841971693993751058209749445923078164
Computer π is 3.1415926535897932384626433832795028841971693993751058209749445923078164
</pre>
=={{header|Kotlin}}==
{{trans|Java}}
<syntaxhighlight lang="Kotlin">
import java.math.BigDecimal
import java.math.BigInteger
import java.math.MathContext
import java.math.RoundingMode
object CodeKt{
@JvmStatic
fun main(args: Array<String>) {
println("n Integer part")
println("================================================")
for (n in 0..9) {
println(String.format("%d%47s", n, almkvistGiullera(n).toString()))
}
val decimalPlaces = 70
val mathContext = MathContext(decimalPlaces + 1, RoundingMode.HALF_EVEN)
val epsilon = BigDecimal.ONE.divide(BigDecimal.TEN.pow(decimalPlaces))
var previous = BigDecimal.ONE
var sum = BigDecimal.ZERO
var pi = BigDecimal.ZERO
var n = 0
while (pi.subtract(previous).abs().compareTo(epsilon) >= 0) {
val nextTerm = BigDecimal(almkvistGiullera(n)).divide(BigDecimal.TEN.pow(6 * n + 3), mathContext)
sum = sum.add(nextTerm)
previous = pi
n += 1
pi = BigDecimal.ONE.divide(sum, mathContext).sqrt(mathContext)
}
println("\npi to $decimalPlaces decimal places:")
println(pi)
}
private fun almkvistGiullera(aN: Int): BigInteger {
val term1 = factorial(6 * aN) * BigInteger.valueOf(32)
val term2 = BigInteger.valueOf(532L * aN * aN + 126 * aN + 9)
val term3 = factorial(aN).pow(6) * BigInteger.valueOf(3)
return term1 * term2 / term3
}
private fun factorial(aNumber: Int): BigInteger {
var result = BigInteger.ONE
for (i in 2..aNumber) {
result *= BigInteger.valueOf(i.toLong())
}
return result
}
private fun BigDecimal.sqrt(context: MathContext): BigDecimal {
var x = BigDecimal(Math.sqrt(this.toDouble()), context)
if (this == BigDecimal.ZERO) return BigDecimal.ZERO
val two = BigDecimal.valueOf(2)
while (true) {
val y = this.divide(x, context)
x = x.add(y).divide(two, context)
val nextY = this.divide(x, context)
if (y == nextY || y == nextY.add(BigDecimal.ONE.divide(BigDecimal.TEN.pow(context.precision), context))) {
break
}
}
return x
}
}
</syntaxhighlight>
{{out}}
<pre>
n Integer part
================================================
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
</pre>
=={{header|Mathematica}}/{{header|Wolfram Language}}==
<
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]</
{{out}}
<pre>{96,5122560,190722470400,7574824857600000,312546150372456000000,13207874703225491420651520,567273919793089083292259942400,24650600248172987140112763715584000,1080657854354639453670407474439566400000,47701779391594966287470570490839978880000000}
Line 712 ⟶ 1,907:
{{libheader|nim-decimal}}
Derived from Wren version with some simplifications.
<
import decimal
Line 749 ⟶ 1,944:
inc n
let pi = 1 / sqrt(sum)
echo ($pi)[0..71]</
{{out}}
Line 770 ⟶ 1,965:
=={{header|Perl}}==
{{trans|Raku}}
<
use warnings;
use feature qw(say);
Line 807 ⟶ 2,002:
printf("π to %s decimal places is:\n%s\n",
$precision, almkvist_giullera_pi($precision));</
{{out}}
<pre>First 10 integer portions:
Line 824 ⟶ 2,019:
=={{header|Phix}}==
<!--<
<span style="color: #008080;">with</span> <span style="color: #008080;">javascript_semantics</span>
<span style="color: #7060A8;">requires</span><span style="color: #0000FF;">(</span><span style="color: #008000;">"1.0.0"</span><span style="color: #0000FF;">)</span>
Line 876 ⟶ 2,071:
<span style="color: #7060A8;">mpfr_const_pi</span><span style="color: #0000FF;">(</span><span style="color: #000000;">pi</span><span style="color: #0000FF;">)</span>
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"Pi (builtin) : %s\n"</span><span style="color: #0000FF;">,</span><span style="color: #7060A8;">mpfr_get_fixed</span><span style="color: #0000FF;">(</span><span style="color: #000000;">pi</span><span style="color: #0000FF;">,</span><span style="color: #000000;">70</span><span style="color: #0000FF;">))</span>
<!--</
{{out}}
<pre>
Line 896 ⟶ 2,091:
Pi to 70 d.p.: 3.1415926535897932384626433832795028841971693993751058209749445923078164
Pi (builtin) : 3.1415926535897932384626433832795028841971693993751058209749445923078164
</pre>
=={{header|PicoLisp}}==
<syntaxhighlight lang="picolisp">(scl 70)
(de fact (N)
(if (=0 N)
1
(* N (fact (dec N))) ) )
(de almkvist (N)
(let
(A (* 32 (fact (* 6 N)))
B (+ (* 532 N N) (* 126 N) 9)
C (* (** (fact N) 6) 3) )
(/ (* A B) C) ) )
(de integral (N)
(*/
1.0
(almkvist N)
(** 10 (+ 3 (* 6 N))) ) )
(let (S 0 N -1)
(do 10
(println (inc 'N) (almkvist N)) )
(prinl)
(setq N -1)
(while (gt0 (integral (inc 'N)))
(inc 'S @) )
(setq S (sqrt (*/ 1.0 1.0 S) 1.0))
(prinl "Pi to 70 decimal places is:")
(prinl (format S *Scl)) )</syntaxhighlight>
{{out}}
<pre>
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.1415926535897932384626433832795028841971693993751058209749445923078152
</pre>
=={{header|Python}}==
<
with mp.workdps(72):
Line 934 ⟶ 2,173:
print('mpmath π is ', end='')
mp.nprint(mp.pi, 71)
</
<pre>
Term 0 96
Line 953 ⟶ 2,192:
=={{header|Quackery}}==
<
[ 1 swap times [ i^ 1+ * ] ] is ! ( n --> n )
Line 971 ⟶ 2,210:
53 times [ i^ vterm v+ ]
1/v 70 vsqrt drop
70 point$ echo$ cr</
{{out}}
Line 990 ⟶ 2,229:
=={{header|Raku}}==
<syntaxhighlight lang="raku"
use BigRoot;
Line 1,021 ⟶ 2,260:
loop { $target eq ( my $next = Pi ++$Nth ) ?? ( last ) !! $target = $next }
say "π to $Precision decimal places is :\n$target"</
{{out}}
<pre>First 10 integer portions :
Line 1,038 ⟶ 2,277:
=={{header|REXX}}==
<
numeric digits length( pi() ) + length(.); w= 102
say $( , 3) $( , w%2) $('power', 5) $( , w)
Line 1,075 ⟶ 2,314:
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</
{{out|output|text= when using the internal default input:}}
(Shown at two─thirds size.)
Line 1,105 ⟶ 2,344:
3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821480865132823066470938446095505822317253594081284811174503
────────────────────────────────────────────── ↑↑↑ the true pi to 160 fractional decimal digits (above) is ↑↑↑ ───────────────────────────────────────────────
</pre>
=={{header|RPL}}==
{{works with|HP|49}}
The starred formula can be implemented like this:
« → n
« 32 6 n * FACT *
532 n SQ * 126 n * + 9 + *
3 / n FACT 6 ^ /
» » '<span style="color:blue">ALMKV</span>' STO
or like that:
« → n '2^5*(6*n)!*(532*SQ(n)+126*n+9)/(3*n!^6)'
» '<span style="color:blue">ALMKV</span>' STO
which is more readable, although a little bit (0.6%) slower than the pure reverse Polish approach.
<code>LDIVN</code> is defined at [[Metallic ratios#RPL|Metallic ratios]]
« 0 → x t
« 0 1
'''WHILE''' DUP x ≤ '''REPEAT''' 4 * '''END'''
'''WHILE''' DUP 1 >
'''REPEAT''' 4 IQUOT
DUP2 + x SWAP - 't' STO
SWAP 2 IQUOT SWAP
'''IF''' t 0 ≥ '''THEN''' t 'x' STO SWAP OVER + SWAP '''END'''
'''END''' DROP
» » '<span style="color:blue">ISQRT</span>' STO
« -105 CF <span style="color:grey">@ set exact mode</span>
« n <span style="color:blue">ALMKV</span> » 'n' 0 9 1 SEQ
-1 → j
« 0 ""
'''DO''' SWAP 'j' INCR <span style="color:blue">ALMKV</span> 10 6 j * 3 + ^ / + EVAL
'''UNTIL''' DUP FXND <span style="color:blue">ISQRT</span> SWAP <span style="color:blue">ISQRT</span> 70 <span style="color:blue">LDIVN</span> ROT OVER ==
'''END'''
"π" →TAG NIP j "iterations" →TAG
» » '<span style="color:blue">TASK</span>' STO
{{out}}
<pre>
3: { 96 5122560 190722470400 7574824857600000 312546150372456000000 13207874703225491420651520 567273919793089083292259942400 24650600248172987140112763715584000 1080657854354639453670407474439566400000 47701779391594966287470570490839978880000000 }
2: π:"3.1415926535897932384626433832795028841971693993751058209749445923078164"
1: iterations:53
</pre>
=={{header|Scala}}==
{{trans|Java}}
<syntaxhighlight lang="Scala">
import java.math.{BigDecimal, BigInteger, MathContext, RoundingMode}
object AlmkvistGiulleraFormula extends App {
println("n Integer part")
println("================================================")
for (n <- 0 to 9) {
val term = almkvistGiullera(n).toString
println(f"$n%1d" + " " * (47 - term.length) + term)
}
val decimalPlaces = 70
val mathContext = new MathContext(decimalPlaces + 1, RoundingMode.HALF_EVEN)
val epsilon = BigDecimal.ONE.divide(BigDecimal.TEN.pow(decimalPlaces))
var previous = BigDecimal.ONE
var sum = BigDecimal.ZERO
var pi = BigDecimal.ZERO
var n = 0
while (pi.subtract(previous).abs.compareTo(epsilon) >= 0) {
val nextTerm = new BigDecimal(almkvistGiullera(n)).divide(BigDecimal.TEN.pow(6 * n + 3))
sum = sum.add(nextTerm)
previous = pi
n += 1
pi = BigDecimal.ONE.divide(sum, mathContext).sqrt(mathContext)
}
println("\npi to " + decimalPlaces + " decimal places:")
println(pi)
def almkvistGiullera(aN: Int): BigInteger = {
val term1 = factorial(6 * aN).multiply(BigInteger.valueOf(32))
val term2 = BigInteger.valueOf(532 * aN * aN + 126 * aN + 9)
val term3 = factorial(aN).pow(6).multiply(BigInteger.valueOf(3))
term1.multiply(term2).divide(term3)
}
def factorial(aNumber: Int): BigInteger = {
var result = BigInteger.ONE
for (i <- 2 to aNumber) {
result = result.multiply(BigInteger.valueOf(i))
}
result
}
}
</syntaxhighlight>
{{out}}
<pre>
n Integer part
================================================
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
</pre>
=={{header|Sidef}}==
<
(32 * (14*n * (38*n + 9) + 9) * (6*n)!) / (3 * n!**6)
}
Line 1,136 ⟶ 2,484:
say "π to #{n} decimal places is:"
say almkvist_giullera_pi(n)
}</
{{out}}
<pre>
Line 1,157 ⟶ 2,505:
{{libheader|System.Numerics}}
{{trans|C#}}
<
Module Module1
Line 1,197 ⟶ 2,545:
End Sub
End Module</
{{out}}
<pre> 0 9600000000000000000000000000000000000000000000000000000000000
Line 1,217 ⟶ 2,565:
{{libheader|Wren-big}}
{{libheader|Wren-fmt}}
<
import "./fmt" for Fmt
var factorial = Fn.new { |n|
Line 1,260 ⟶ 2,608:
var pi = BigRat.one/sum.sqrt(70)
System.print("\nPi to 70 decimal places is:")
System.print(pi.toDecimal(70))</
{{out}}
| |||