Almkvist-Giullera formula for pi: Difference between revisions
Add Kotlin implementation
Thundergnat (talk | contribs) m (argh.) |
(Add Kotlin implementation) |
||
| (10 intermediate revisions by 6 users not shown) | |||
Line 34:
{{trans|C#}}
<syntaxhighlight lang="11l">F isqrt(BigInt =x)
BigInt q = 1
BigInt r = 0
Line 113:
=={{header|AArch64 Assembly}}==
{{works with|as|Raspberry Pi 3B version Buster 64 bits <br> or android 64 bits with application Termux }}
<syntaxhighlight lang=
/* ARM assembly AARCH64 Raspberry PI 3B */
/* program calculPi64.s */
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}}
<syntaxhighlight lang=
/* ARM assembly Raspberry PI */
/* program calculPi.s */
Line 724 ⟶ 838:
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).
<syntaxhighlight lang="csharp">using System;
using BI = System.Numerics.BigInteger;
using static System.Console;
Line 780 ⟶ 895:
{{libheader|Boost}}
{{libheader|GMP}}
<syntaxhighlight lang="cpp">#include <boost/multiprecision/cpp_dec_float.hpp>
#include <boost/multiprecision/gmp.hpp>
#include <iomanip>
Line 847 ⟶ 962:
{{libheader|computable-reals}}
{{trans|Raku}}
<syntaxhighlight lang="lisp">(ql:quickload :computable-reals :silent t)
(use-package :computable-reals)
(setq *print-prec* 70)
Line 912 ⟶ 1,027:
=={{header|dc}}==
{{trans|Common Lisp}}
<syntaxhighlight lang="dc">[* factorial *]sz
[ 1 Sp [ d lp * sp 1 - d 1 <f ]Sf d 1 <f Lfsz sz Lp ]sF
Line 977 ⟶ 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=
-mode(compile).
Line 1,105 ⟶ 1,220:
=={{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
Line 1,131 ⟶ 1,246:
=={{header|Factor}}==
{{works with|Factor|0.99 2020-08-14}}
<syntaxhighlight lang="factor">USING: continuations formatting io kernel locals math
math.factorials math.functions sequences ;
Line 1,195 ⟶ 1,310:
=={{header|Go}}==
{{trans|Wren}}
<syntaxhighlight lang="go">package main
import (
Line 1,285 ⟶ 1,400:
{{libheader|numbers}}
{{trans|Common Lisp}}
<syntaxhighlight lang="haskell">import Control.Monad
import Data.Number.CReal
import GHC.Integer
Line 1,353 ⟶ 1,468:
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=
numerator =: monad define "0
(3 * (! x: y)^6) %~ 32 * (!6x*y) * (y*(126 + 532*y)) + 9x
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,411 ⟶ 1,600:
{{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';
Line 1,491 ⟶ 1,680:
'''Preliminaries'''
<syntaxhighlight lang="jq"># A reminder to include the "rational" module:
# include "rational";
Line 1,504 ⟶ 1,693:
end; </syntaxhighlight>
'''Almkvist-Giullera Formula'''
<syntaxhighlight lang="jq">
def almkvistGiullera(print):
. as $n
Line 1,518 ⟶ 1,707:
end; </syntaxhighlight>
'''The Tasks'''
<syntaxhighlight lang="jq">
def task1:
"N Integer Portion Pow Nth Term",
Line 1,562 ⟶ 1,751:
=={{header|Julia}}==
<syntaxhighlight lang="julia">using Formatting
setprecision(BigFloat, 300)
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>
=={{header|Mathematica}}/{{header|Wolfram Language}}==
<syntaxhighlight lang=
numerator[n_] := (2^5) ((6 n)!) (532 n^2 + 126 n + 9)/(3 (n!)^6)
denominator[n_] := 10^(6 n + 3)
Line 1,629 ⟶ 1,907:
{{libheader|nim-decimal}}
Derived from Wren version with some simplifications.
<syntaxhighlight lang=
import decimal
Line 1,687 ⟶ 1,965:
=={{header|Perl}}==
{{trans|Raku}}
<syntaxhighlight lang="perl">use strict;
use warnings;
use feature qw(say);
Line 1,741 ⟶ 2,019:
=={{header|Phix}}==
<!--<syntaxhighlight lang=
<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 1,816 ⟶ 2,094:
=={{header|PicoLisp}}==
<syntaxhighlight lang=
(de fact (N)
(if (=0 N)
Line 1,860 ⟶ 2,138:
=={{header|Python}}==
<syntaxhighlight lang="python">import mpmath as mp
with mp.workdps(72):
Line 1,914 ⟶ 2,192:
=={{header|Quackery}}==
<syntaxhighlight lang=
[ 1 swap times [ i^ 1+ * ] ] is ! ( n --> n )
Line 1,951 ⟶ 2,229:
=={{header|Raku}}==
<syntaxhighlight lang=
use BigRoot;
Line 1,999 ⟶ 2,277:
=={{header|REXX}}==
<syntaxhighlight 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)
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>
=={{header|Sidef}}==
<syntaxhighlight lang="ruby">func almkvist_giullera(n) {
(32 * (14*n * (38*n + 9) + 9) * (6*n)!) / (3 * n!**6)
}
Line 2,118 ⟶ 2,505:
{{libheader|System.Numerics}}
{{trans|C#}}
<syntaxhighlight lang="vbnet">Imports System, BI = System.Numerics.BigInteger, System.Console
Module Module1
Line 2,178 ⟶ 2,565:
{{libheader|Wren-big}}
{{libheader|Wren-fmt}}
<syntaxhighlight lang=
import "./fmt" for Fmt
var factorial = Fn.new { |n|
| |||