Arithmetic-geometric mean/Calculate Pi: Difference between revisions
Content added Content deleted
Thundergnat (talk | contribs) m (syntax highlighting fixup automation) |
Thundergnat (talk | contribs) m (Automated syntax highlighting fixup (second round - minor fixes)) |
||
Line 1: | Line 1: | ||
⚫ | |||
{{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= |
<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= |
<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= |
<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= |
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= |
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= |
<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= |
<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= |
<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= |
<!--<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= |
<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= |
<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= |
<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, it is over '''68''' times faster than version 3. |
For 1,005 decimal digits, it is over '''68''' 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= |
<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> |
||
⚫ |