Arithmetic-geometric mean/Calculate Pi: Difference between revisions

Content added Content deleted
m (syntax highlighting fixup automation)
m (Automated syntax highlighting fixup (second round - minor fixes))
Line 1: Line 1:
[[Category:Geometry]]
{{task}}
{{task}}
[http://www.maa.org/sites/default/files/pdf/upload_library/22/Ford/Almkvist-Berndt585-608.pdf Almkvist Berndt 1988] begins with an investigation of why the agm is such an efficient algorithm, and proves that it converges quadratically. This is an efficient method to calculate <math>\pi</math>.
[http://www.maa.org/sites/default/files/pdf/upload_library/22/Ford/Almkvist-Berndt585-608.pdf Almkvist Berndt 1988] begins with an investigation of why the agm is such an efficient algorithm, and proves that it converges quadratically. This is an efficient method to calculate <math>\pi</math>.
Line 22: Line 23:
<!-- Example by Nigel Galloway, Feb 8, 2012 -->
<!-- Example by Nigel Galloway, Feb 8, 2012 -->
{{libheader|GMP}}
{{libheader|GMP}}
<syntaxhighlight lang=c>#include "gmp.h"
<syntaxhighlight lang="c">#include "gmp.h"


void agm (const mpf_t in1, const mpf_t in2, mpf_t out1, mpf_t out2) {
void agm (const mpf_t in1, const mpf_t in2, mpf_t out1, mpf_t out2) {
Line 74: Line 75:


=={{header|C sharp|C#}}==
=={{header|C sharp|C#}}==
{{libheader|System.Numerics}}Can specify the number of desired digits on the command line, default is 25000, which takes a few seconds (depending on your system's performance).<syntaxhighlight lang=csharp>using System;
{{libheader|System.Numerics}}Can specify the number of desired digits on the command line, default is 25000, which takes a few seconds (depending on your system's performance).<syntaxhighlight lang="csharp">using System;
using System.Numerics;
using System.Numerics;


Line 140: Line 141:
{{trans|C}}
{{trans|C}}
{{libheader|GMP}}
{{libheader|GMP}}
<syntaxhighlight lang=cpp>#include <gmpxx.h>
<syntaxhighlight lang="cpp">#include <gmpxx.h>
#include <chrono>
#include <chrono>


Line 190: Line 191:
=={{header|Clojure}}==
=={{header|Clojure}}==
{{trans|Ruby}}
{{trans|Ruby}}
<syntaxhighlight lang=lisp>(ns async-example.core
<syntaxhighlight lang="lisp">(ns async-example.core
(:use [criterium.core])
(:use [criterium.core])
(:gen-class))
(:gen-class))
Line 266: Line 267:
{{libheader|MMA}}
{{libheader|MMA}}
<p>This is an example that uses the Common Lisp Bigfloat Package (http://www.cs.berkeley.edu/~fateman/lisp/mma4max/more/bf.lisp)</p>
<p>This is an example that uses the Common Lisp Bigfloat Package (http://www.cs.berkeley.edu/~fateman/lisp/mma4max/more/bf.lisp)</p>
<syntaxhighlight lang=lisp>(load "bf.fasl")
<syntaxhighlight lang="lisp">(load "bf.fasl")


;;(setf mma::bigfloat-bin-prec 1000)
;;(setf mma::bigfloat-bin-prec 1000)
Line 289: Line 290:
=={{header|D}}==
=={{header|D}}==
{{trans|C#}}
{{trans|C#}}
<syntaxhighlight lang=d>import std.bigint;
<syntaxhighlight lang="d">import std.bigint;
import std.conv;
import std.conv;
import std.math;
import std.math;
Line 372: Line 373:
Thanks for Velthuis BigIntegers library[https://github.com/rvelthuis/DelphiBigNumbers].
Thanks for Velthuis BigIntegers library[https://github.com/rvelthuis/DelphiBigNumbers].
{{Trans|C#}}
{{Trans|C#}}
<syntaxhighlight lang=Delphi>
<syntaxhighlight lang="delphi">
program Calculate_Pi;
program Calculate_Pi;


Line 490: Line 491:
=={{header|Erlang}}==
=={{header|Erlang}}==
{{trans|python}}
{{trans|python}}
<syntaxhighlight lang=erlang>
<syntaxhighlight lang="erlang">
-module(pi).
-module(pi).
-export([agmPi/1, agmPiBody/5]).
-export([agmPi/1, agmPiBody/5]).
Line 522: Line 523:
{{libheader|iso_fortran_env}}
{{libheader|iso_fortran_env}}
{{trans|Julia}}
{{trans|Julia}}
<syntaxhighlight lang=fortran>program CalcPi
<syntaxhighlight lang="fortran">program CalcPi
! Use real128 numbers: (append '_rf')
! Use real128 numbers: (append '_rf')
use iso_fortran_env, only: rf => real128
use iso_fortran_env, only: rf => real128
Line 571: Line 572:


=={{header|Go}}==
=={{header|Go}}==
<syntaxhighlight lang=go>package main
<syntaxhighlight lang="go">package main


import (
import (
Line 618: Line 619:
=={{header|Groovy}}==
=={{header|Groovy}}==
{{trans|Java}}
{{trans|Java}}
<syntaxhighlight lang=groovy>import java.math.MathContext
<syntaxhighlight lang="groovy">import java.math.MathContext


class CalculatePi {
class CalculatePi {
Line 658: Line 659:
{{libheader|MPFR}}
{{libheader|MPFR}}
{{libheader|hmpfr}}
{{libheader|hmpfr}}
<syntaxhighlight lang=Haskell>import Prelude hiding (pi)
<syntaxhighlight lang="haskell">import Prelude hiding (pi)
import Data.Number.MPFR hiding (sqrt, pi, div)
import Data.Number.MPFR hiding (sqrt, pi, div)
import Data.Number.MPFR.Instances.Near ()
import Data.Number.MPFR.Instances.Near ()
Line 698: Line 699:
[https://code.jsoftware.com/wiki/Essays/Extended%20Precision%20Functions Extended precision functions] and [https://code.jsoftware.com/wiki/Essays/Chudnovsky_Algorithm Pi].
[https://code.jsoftware.com/wiki/Essays/Extended%20Precision%20Functions Extended precision functions] and [https://code.jsoftware.com/wiki/Essays/Chudnovsky_Algorithm Pi].
Translated from python:
Translated from python:
<syntaxhighlight lang=J>DP=: 100
<syntaxhighlight lang="j">DP=: 100


round=: DP&$: : (4 : 0)
round=: DP&$: : (4 : 0)
Line 745: Line 746:
1583455951826865080542496790424338362837978447536228171662934224565463064033895909488933268392567279887495006936541219489670405121434573776487989539520749180843985094860051126840117004097133550161882511486508109869673199973040182062140382647367514024790194...</pre>
1583455951826865080542496790424338362837978447536228171662934224565463064033895909488933268392567279887495006936541219489670405121434573776487989539520749180843985094860051126840117004097133550161882511486508109869673199973040182062140382647367514024790194...</pre>


That said, note that J offers a more direct approach here:<syntaxhighlight lang=J> 102j100":<.@o.&.(*&(10^100x))1
That said, note that J offers a more direct approach here:<syntaxhighlight lang="j"> 102j100":<.@o.&.(*&(10^100x))1
3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679</syntaxhighlight>
3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679</syntaxhighlight>


However, a limitation of this approach is that the floor function -- which we use here to signal the desire for an exact approximation of pi -- does not round. But we can ask for extra digits, for example:<syntaxhighlight lang=J> 113j111":<.@o.&.(*&(10^111x))1
However, a limitation of this approach is that the floor function -- which we use here to signal the desire for an exact approximation of pi -- does not round. But we can ask for extra digits, for example:<syntaxhighlight lang="j"> 113j111":<.@o.&.(*&(10^111x))1
3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117067982148086513</syntaxhighlight>
3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117067982148086513</syntaxhighlight>


Line 756: Line 757:
{{trans|Kotlin}}
{{trans|Kotlin}}
Used features of Java 8
Used features of Java 8
<syntaxhighlight lang=Java>import java.math.BigDecimal;
<syntaxhighlight lang="java">import java.math.BigDecimal;
import java.math.MathContext;
import java.math.MathContext;
import java.util.Objects;
import java.util.Objects;
Line 803: Line 804:
`rsqrt` function
`rsqrt` function
for taking square roots of rationals; such a module can be found at [[Arithmetic/Rational]].
for taking square roots of rationals; such a module can be found at [[Arithmetic/Rational]].
<syntaxhighlight lang=jq># include "rational"; # a reminder
<syntaxhighlight lang="jq"># include "rational"; # a reminder


def pi(precision):
def pi(precision):
Line 830: Line 831:


=={{header|Julia}}==
=={{header|Julia}}==
<syntaxhighlight lang=julia>using Printf
<syntaxhighlight lang="julia">using Printf


agm1step(x, y) = (x + y) / 2, sqrt(x * y)
agm1step(x, y) = (x + y) / 2, sqrt(x * y)
Line 879: Line 880:


=={{header|Kotlin}}==
=={{header|Kotlin}}==
<syntaxhighlight lang=scala>import java.math.BigDecimal
<syntaxhighlight lang="scala">import java.math.BigDecimal
import java.math.MathContext
import java.math.MathContext


Line 919: Line 920:
With 16 steps we have 89415 correct digits:
With 16 steps we have 89415 correct digits:


<syntaxhighlight lang=maple>agm:=proc(n)
<syntaxhighlight lang="maple">agm:=proc(n)
local a:=1,g:=evalf(sqrt(1/2)),s:=0,p:=4,i;
local a:=1,g:=evalf(sqrt(1/2)),s:=0,p:=4,i;
for i to n do
for i to n do
Line 935: Line 936:


=={{header|Mathematica}}/{{header|Wolfram Language}}==
=={{header|Mathematica}}/{{header|Wolfram Language}}==
<syntaxhighlight lang=mathematica>pi[n_, prec_] :=
<syntaxhighlight lang="mathematica">pi[n_, prec_] :=
Module[{a = 1, g = N[1/Sqrt[2], prec], k, s = 0, p = 4},
Module[{a = 1, g = N[1/Sqrt[2], prec], k, s = 0, p = 4},
For[k = 1, k < n, k++,
For[k = 1, k < n, k++,
Line 950: Line 951:
=={{header|МК-61/52}}==
=={{header|МК-61/52}}==
{{Output?}}
{{Output?}}
<lang>3 П0 1 П1 П4 2 КвКор 1/x П2 1
<syntaxhighlight lang="text">3 П0 1 П1 П4 2 КвКор 1/x П2 1
^ 4 / П3 ИП3 ИП1 ИП2 + 2 /
^ 4 / П3 ИП3 ИП1 ИП2 + 2 /
П5 ИП1 - x^2 ИП4 * - П3 ИП1 ИП2
П5 ИП1 - x^2 ИП4 * - П3 ИП1 ИП2
Line 959: Line 960:
{{libheader|bignum}}
{{libheader|bignum}}
{{Trans|Delphy}}
{{Trans|Delphy}}
<syntaxhighlight lang=Nim>from math import sqrt
<syntaxhighlight lang="nim">from math import sqrt
import times
import times


Line 1,057: Line 1,058:
=={{header|OCaml}}==
=={{header|OCaml}}==
program for calculating digits of pi
program for calculating digits of pi
<syntaxhighlight lang=OCaml>let limit = 10000 and n = 2800
<syntaxhighlight lang="ocaml">let limit = 10000 and n = 2800
let x = Array.make (n+1) 2000
let x = Array.make (n+1) 2000


Line 1,080: Line 1,081:


=={{header|PARI/GP}}==
=={{header|PARI/GP}}==
<syntaxhighlight lang=parigp>pi(n)=my(a=1,g=2^-.5);(1-2*sum(k=1,n,[a,g]=[(a+g)/2,sqrt(a*g)];(a^2-g^2)<<k))^-1*4*a^2
<syntaxhighlight lang="parigp">pi(n)=my(a=1,g=2^-.5);(1-2*sum(k=1,n,[a,g]=[(a+g)/2,sqrt(a*g)];(a^2-g^2)<<k))^-1*4*a^2
pi(6)</syntaxhighlight>
pi(6)</syntaxhighlight>
{{out}}
{{out}}
Line 1,090: Line 1,091:
The number of steps used is based on the desired accuracy rather than being hard coded, as this is intended to work for 1M digits as well as for 100.
The number of steps used is based on the desired accuracy rather than being hard coded, as this is intended to work for 1M digits as well as for 100.


<syntaxhighlight lang=perl>use Math::BigFloat try => "GMP,Pari";
<syntaxhighlight lang="perl">use Math::BigFloat try => "GMP,Pari";


my $digits = shift || 100; # Get number of digits from command line
my $digits = shift || 100; # Get number of digits from command line
Line 1,117: Line 1,118:
The following is a translation, almost line-for-line, of the Ruby code. It is slower than the above and the last digit or two may not be correct.
The following is a translation, almost line-for-line, of the Ruby code. It is slower than the above and the last digit or two may not be correct.
{{trans|Ruby}}
{{trans|Ruby}}
<syntaxhighlight lang=perl>use strict;
<syntaxhighlight lang="perl">use strict;
use warnings;
use warnings;
use Math::BigFloat;
use Math::BigFloat;
Line 1,140: Line 1,141:
{{libheader|Phix/mpfr}}
{{libheader|Phix/mpfr}}
{{trans|Python}}
{{trans|Python}}
<!--<syntaxhighlight lang=Phix>(phixonline)-->
<!--<syntaxhighlight lang="phix">(phixonline)-->
<span style="color: #008080;">with</span> <span style="color: #008080;">javascript_semantics</span>
<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> <span style="color: #000080;font-style:italic;">-- (mpfr_set_default_prec[ision] has been renamed)</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> <span style="color: #000080;font-style:italic;">-- (mpfr_set_default_prec[ision] has been renamed)</span>
Line 1,204: Line 1,205:
=={{header|PicoLisp}}==
=={{header|PicoLisp}}==
{{trans|Python}}
{{trans|Python}}
<syntaxhighlight lang=PicoLisp>(scl 40)
<syntaxhighlight lang="picolisp">(scl 40)


(de pi ()
(de pi ()
Line 1,230: Line 1,231:
=={{header|Python}}==
=={{header|Python}}==
{{trans|Ruby}}
{{trans|Ruby}}
<syntaxhighlight lang=python>from decimal import *
<syntaxhighlight lang="python">from decimal import *


D = Decimal
D = Decimal
Line 1,251: Line 1,252:
{{libheader|Rmpfr}}
{{libheader|Rmpfr}}


<syntaxhighlight lang=rsplus>library(Rmpfr)
<syntaxhighlight lang="rsplus">library(Rmpfr)


agm <- function(n, prec) {
agm <- function(n, prec) {
Line 1,287: Line 1,288:
=={{header|Racket}}==
=={{header|Racket}}==
{{trans|Ruby}}
{{trans|Ruby}}
<syntaxhighlight lang=Racket>#lang racket
<syntaxhighlight lang="racket">#lang racket
(require math/bigfloat)
(require math/bigfloat)


Line 1,334: Line 1,335:


Notice that we don't get the exact number of decimals required : the last two decimals or so can be wrong. This is because we don't need <math>a_n</math>, but rather <math>a_n^2</math>. Elevating to the square makes us lose a bit of precision. It could be compensated by choosing a slightly higher value of N (in a way that could be precisely calculated), but that would probably be overkill.
Notice that we don't get the exact number of decimals required : the last two decimals or so can be wrong. This is because we don't need <math>a_n</math>, but rather <math>a_n^2</math>. Elevating to the square makes us lose a bit of precision. It could be compensated by choosing a slightly higher value of N (in a way that could be precisely calculated), but that would probably be overkill.
<syntaxhighlight lang=perl6>constant number-of-decimals = 100;
<syntaxhighlight lang="raku" line>constant number-of-decimals = 100;


multi sqrt(Int $n) {
multi sqrt(Int $n) {
Line 1,381: Line 1,382:
Whatever number of digits used, the actual number of digits is five larger than specified, and then the result is rounded to the requested number of digits.
Whatever number of digits used, the actual number of digits is five larger than specified, and then the result is rounded to the requested number of digits.
===version 1===
===version 1===
<syntaxhighlight lang=rexx>/*REXX program calculates the value of pi using the AGM algorithm. */
<syntaxhighlight lang="rexx">/*REXX program calculates the value of pi using the AGM algorithm. */
parse arg d .; if d=='' | d=="," then d= 500 /*D not specified? Then use default. */
parse arg d .; if d=='' | d=="," then d= 500 /*D not specified? Then use default. */
numeric digits d+5 /*set the numeric decimal digits to D+5*/
numeric digits d+5 /*set the numeric decimal digits to D+5*/
Line 1,413: Line 1,414:


For 1,005 decimal digits, &nbsp; it is over &nbsp; '''68''' &nbsp; times faster than version 3.
For 1,005 decimal digits, &nbsp; it is over &nbsp; '''68''' &nbsp; times faster than version 3.
<syntaxhighlight lang=rexx>/*REXX program calculates the AGM (arithmetic─geometric mean) of two (real) numbers. */
<syntaxhighlight lang="rexx">/*REXX program calculates the AGM (arithmetic─geometric mean) of two (real) numbers. */
parse arg a b digs . /*obtain optional numbers from the C.L.*/
parse arg a b digs . /*obtain optional numbers from the C.L.*/
if digs=='' | digs=="," then digs= 100 /*No DIGS specified? Then use default.*/
if digs=='' | digs=="," then digs= 100 /*No DIGS specified? Then use default.*/
Line 1,480: Line 1,481:


===version 3===
===version 3===
<syntaxhighlight lang=rexx>
<syntaxhighlight lang="rexx">
/*REXX*/
/*REXX*/
Line 1,546: Line 1,547:
Using agm.
Using agm.
See [[Talk:Arithmetic-geometric mean]]
See [[Talk:Arithmetic-geometric mean]]
<syntaxhighlight lang=ruby># Calculate Pi using the Arithmetic Geometric Mean of 1 and 1/sqrt(2)
<syntaxhighlight lang="ruby"># Calculate Pi using the Arithmetic Geometric Mean of 1 and 1/sqrt(2)
#
#
#
#
Line 1,602: Line 1,603:


=={{header|Rust}}==
=={{header|Rust}}==
<syntaxhighlight lang=rust>/// calculate pi with algebraic/geometric mean
<syntaxhighlight lang="rust">/// calculate pi with algebraic/geometric mean
pub fn pi(n: usize) -> f64 {
pub fn pi(n: usize) -> f64 {
let mut a : f64 = 1.0;
let mut a : f64 = 1.0;
Line 1,625: Line 1,626:
</syntaxhighlight>
</syntaxhighlight>
Can be invoked like:
Can be invoked like:
<syntaxhighlight lang=rust>
<syntaxhighlight lang="rust">
fn main() {
fn main() {
println!("pi(7): {}", pi(7));
println!("pi(7): {}", pi(7));
Line 1,637: Line 1,638:
=={{header|Scala}}==
=={{header|Scala}}==
===Completely (tail) recursive===
===Completely (tail) recursive===
<syntaxhighlight lang=Scala>import java.math.MathContext
<syntaxhighlight lang="scala">import java.math.MathContext


import scala.annotation.tailrec
import scala.annotation.tailrec
Line 1,676: Line 1,677:


=={{header|Sidef}}==
=={{header|Sidef}}==
<syntaxhighlight lang=ruby>func agm_pi(digits) {
<syntaxhighlight lang="ruby">func agm_pi(digits) {
var acc = (digits + 8);
var acc = (digits + 8);


Line 1,705: Line 1,706:
{{trans|Ruby}}
{{trans|Ruby}}
{{tcllib|math::bigfloat}}
{{tcllib|math::bigfloat}}
<syntaxhighlight lang=tcl>package require math::bigfloat
<syntaxhighlight lang="tcl">package require math::bigfloat
namespace import math::bigfloat::*
namespace import math::bigfloat::*


Line 1,735: Line 1,736:
{{trans|C#}}
{{trans|C#}}
{{Libheader|System.Numerics}}
{{Libheader|System.Numerics}}
<syntaxhighlight lang=vbnet>Imports System, System.Numerics
<syntaxhighlight lang="vbnet">Imports System, System.Numerics


Module Program
Module Program
Line 1,797: Line 1,798:
{{trans|Sidef}}
{{trans|Sidef}}
{{libheader|Wren-big}}
{{libheader|Wren-big}}
<syntaxhighlight lang=ecmascript>import "/big" for BigRat
<syntaxhighlight lang="ecmascript">import "/big" for BigRat


var digits = 500
var digits = 500
Line 1,819: Line 1,820:
3.14159265358979323846264338327950288419716939937510582097494459230781640628620899862803482534211706798214808651328230664709384460955058223172535940812848111745028410270193852110555964462294895493038196442881097566593344612847564823378678316527120190914564856692346034861045432664821339360726024914127372458700660631558817488152092096282925409171536436789259036001133053054882046652138414695194151160943305727036575959195309218611738193261179310511854807446237996274956735188575272489122793818301194912
3.14159265358979323846264338327950288419716939937510582097494459230781640628620899862803482534211706798214808651328230664709384460955058223172535940812848111745028410270193852110555964462294895493038196442881097566593344612847564823378678316527120190914564856692346034861045432664821339360726024914127372458700660631558817488152092096282925409171536436789259036001133053054882046652138414695194151160943305727036575959195309218611738193261179310511854807446237996274956735188575272489122793818301194912
</pre>
</pre>

[[Category:Geometry]]