Rosetta Code

Almkvist-Giullera formula for pi: Difference between revisions

Page Discussion

Almkvist-Giullera formula for pi (view source)

Revision as of 04:45, 8 February 2024

13,853 bytes added ,  3 months ago
Add Kotlin implementation
m (Automated syntax highlighting fixup (second round - minor fixes))
(Add Kotlin implementation)
 
(9 intermediate revisions by 5 users not shown)
Line 419:
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}}
Line 724 ⟶ 838:
3.1415926535897932384626433832795028841971693993751058209749445923078164
</pre>
 
=={{header|C sharp|C#}}==
{{libheader|System.Numerics}}
Line 1,403 ⟶ 1,518:
 
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>
 
Line 1,613 ⟶ 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>
 
Line 2,066 ⟶ 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>
 
Line 2,178 ⟶ 2,565:
{{libheader|Wren-big}}
{{libheader|Wren-fmt}}
<syntaxhighlight lang="ecmascriptwren">import "./big" for BigInt, BigRat
import "./fmt" for Fmt
 
var factorial = Fn.new { |n|
Aspen138
337

edits

Retrieved from "https://rosettacode.org/wiki/Special:MobileDiff/328939...354322"