How to calculate numbers to arbitrarily high precision?

Question

I wrote a simple fortran program to compute Gauss's constant :

program main

implicit none

integer :: i, nit
double precision :: u0, v0, ut, vt

nit=60
u0=1.d0
v0=sqrt(2.d0)
print *,1.d0/u0,1.d0/v0

do i=1,nit
  ut=sqrt(u0*v0)
  vt=(u0+v0)/2.d0
  u0=ut
  v0=vt
  print *,1.d0/u0,1.d0/v0
enddo

end program main

Result is 0.83462684167407308 after 4 iterations. Anyway to have better results using the arithmetico-geometric mean method? How do people compute many digits for numbers such as pi, Euler's constant, and so on ? Does each irrational number has a specific algorithm?

Answer 1

If you goal is to insert a constant value into your program, the easiest solution is to look up the value on the web in or a book. Be sure to add a type specification to the numeric value, other Fortran will treat it as the default of single precision. One could write pi as pi_quad = 3.14159265358979323846264338327950288_real128 -- showing the use of a type specifier on a constant.

If you want to do high precision calculations, you could some high precision type available in your compiler. Many compilers now have quadruple precision. If they have the Fortran 2008 version of the ISO_FORTRAN_ENV module, you can request this via the type real128.

Arbitrary precision (user specified number of digits, to very high number of digits) is outside the language and is available in libraries, e.g., MPFUN90, http://crd-legacy.lbl.gov/~dhbailey/mpdist/

Yes, different constants have various algorithms. This is a very large topic.

Answer 2

Solution for pi:

  pi = 4.0d0 * datan(1.0d0)