<- previous    index    next ->

Lecture 15, Multiple precision 


When 64-bit floating point is not accurate enough 
When 64-bit integers are way too small 

You will need this lecture for Homework 4 and
cs455 Project

"C" has gmp, gnu multiple precision. 

Java has a bignum package.
Java has apfloat.jar

Ada has arbitrary decimal digit precision floating point.

Fortran has a multiple precision library.

Python has arbitrary precision integer arithmetic
Python has mpmath and other packages

Hmmm? Multiple precision must be important and have a use.

Computing Pi to a million places is a demonstration, but there
are more reasonable uses.

Types of multiple precision:
  Unlimited length integers
  Unlimited length fractions  1/integer
  Unlimited length rationals  integer/integer
  Unlimited length floating point
  Arbitrary yet fixed number of digits floating point


for "C" get gmp, GNU Multiple Precision!
download gmp from https://gmplib.org
gmp.guide

Here are a few simple gmp samples

test_mpf.c
test_mpf.out
test_mpq.c
test_mpq.out
test_mpz.c
test_mpz.out
gmp fact.c
fact_gmp.out


Java BigDecimal that is BigInteger with a scale factor
Big_pi.java test program
Big_pi_java.out test results
Big_pi_check.java test program
Big_pi_check_java.out test results
test_apfloat.java test program
test_apfloat.out test results
mytest_apfloat.java test program
mytest_apfloat.out test results
Big_math.java various functions exp, sin
Big_simeq.java simultaneous equations
Big_inverse.java invert matrix

Fortran 95 module that implements big integers
big_integers_module.f90
test_big.f90  test program
test_big_f90.out  test results

Ada 95 precise, rational and digits arithmetic
directory of Ada 95 files

Python simple factorial(52) is
factorial.py  program
factorial_py.out  output
factbig.py3  really big 1,000,000!


Python simple 2^200 integer, yet floating point is IEEE 64 bit
power2.py  program
power2_py.out  output
test_mpmath.py  program
test_mpmath_py.out  output
test_mpmath.py3  program
test_mpmath_py3.out  output
mpmath_example.py  program
mpmath_example_py.out  output
mpmath_example.py3  program
mpmath_example_py3.out  output


A quick conversion of simeq.c to mpf_simeq.c solves simultaneous
equations with 50 digits, could be many more digits using gmp mpf_t.

Using the very difficult to get accurate answers matrix:

test_mpf_simeq.c
mpf_simeq.c
mpf_simeq.h
test_mpf_simeq.out
test_mpf_simeq_300.out


Can irrational numbers be combined to produce an integer?

It appears that  e^(Pi*sqrt(163)) is the integer 262,537,412,640,768,744
and that is rather amazing to me.

How might this be validated? It has been tried using higher and
higher precision and the value computes to about
262,537,412,640,768,743.999,999,999,999,999,999,999

We know 1.5 base 10 can be represented as  1.1 base 2 exactly.
We know 1/3 can not be represented exactly in a finite number of
decimal digits or a finite number of binary digits
   1/3 = 0.3333333333333333333333333333333  decimal approximately
   1/3 = 0.0101010101010101010101010101010  binary approximately
And, for what it is worth
We know 1/3 can be represented exactly as  0.1 base 3.

A quick "C" program gives the first 14 digits, using IEEE 64-bit floating point

 /* irrational.c  e^(Pi*sqrt(163)) is an integer?  */
 #include <math.h>
 #include <stdio.h>

 #define e  2.7182818284590452353602874713526624977572
 #define Pi 3.1415926535897932384626433832795028841971

 int main(int argc, char * argv[])
 {
   printf("e^(Pi*sqrt(163)) is an integer? \n");
   printf("  262537412640768744.0000000  check \n");
   printf("%28.7f  using 15 digit arithmetic \n",pow(e,Pi*sqrt(163.0)));
   return 0;
 }
With output:
 e^(Pi*sqrt(163)) is an integer? 
   262537412640768744.0000000  check 
   262537412640767712.0000000  using 15 digit arithmetic 
   262537412640768743.99999999999925007259719818568887935385633733699 using 300 digit

What is needed is to show convergence from above and below,
and it would be nice if the convergence was uniform. This
could use 50, 100, 150, 200 digit precision for the computation.
Pick a number of digits. Increment the bottom digit of e, Pi and 163
then do the computation. Decrement the bottom digit of e, Pi and 163
then do the computation. Check if the upper and lower values of
the fraction are converging toward zero. Check if the convergence
is uniform, balanced, for the upper and lower values.

This is left as an exercise to the student.


Computing Pi to arbitrary precision


One method of computing Pi is to compute 4*atan(1) or
24*atan(b2)  b2=(2*sqrt(b1)-1)/b1  b1=2-sqrt(3)
get_pi.c
get_pi.out
get_pi.c using double and math.h atan 
pi=24*atan(b2) b2=0.131652, b1=0.267949 
    Pi= 3.1415926535897932384626433832795028841971 
get_pi= 3.14159265358980333 
 pi-Pi= 0.00000000000001021 

get_mpf_pi.c
get_mpf_pi.out
get_mpf_pi.c using digits=50 
b1=    0.26794919243112270647255365849412763305719474618962 
b2=    0.13165249758739585347152645740971710359281410222324 
mpf_pi= 3.1415926535897932384626433832795028841971693999457 
    Pi= 3.1415926535897932384626433832795028841971 


Now you can do Homework 4




Using gmp to solve higher power methods


For reference in later lectures on PDE, when "double" is
just not accurate enough to allow higher power methods:

Solving a relatively simple partial differential equation

    du/dx + du/dy = 0

using a uniform grid on rectangle 0 <= X <= 2Pi, 0 <= Y <= 2Pi
with boundary values (and analytic solution)

    u(x,y) = sin(x-y)

Due to the symmetry of the problem, the number of boundary
points in X and Y can not be the same. If nx = ny then
the solution of simultaneous equations becomes unsolvable
because of a singular matrix.

Solutions above nx or ny equal 13 cause unstable values
for the derivative coefficients with standard nderiv
computation. Actually, integer overflow occurs first.

By using gmp with high precision integer, mpz, and
high precision floating point, mpf, then high
accuracy solutions can be computed.

Setting "digits" at 200, "bits" becomes 200*3.32 = 664.
This is used for computing boundary values and the
numeric checking against the analytic solution.

The results are:
  nx  ny  maxerror
  13  12  1.40e-4     same as nx=12, ny=13, they can not be equal
  15  14  4.38e-6
  17  16  1.02e-7
  19  18  1.84e-9
  21  20  2.70e-11
  23  22  3.23e-13
  25  24  3.23e-15
  31  30  1.25e-21
  33  32  6.86e-24
  35  34  3.33e-26
  37  36  1.44e-28

  37  36  8.59e-18  0  to  4Pi
  37  36  3.40e-7   0  to  8Pi

The basic "C" code with lots of printout is:
source code pde2sin_eq.c
output pde2sin_eq_c.out
source code pde2sin_sparse.c
output pde2sin_sparse_c.out

The same programs, converting "double" to gmp "mpf_t"
and calling gmp versions of functions is:
source code pde2sin_eq_gmp.c
output pde2sin_eq_gmp.out
source code pde2sin_sparse_gmp.c
output pde2sin_sparse_gmp.out
gmp function mpf_deriv.h
gmp function mpf_deriv.c
gmp function mpf_simeq.h
gmp function mpf_simeq.c
gmp function mpf_sparse.h
gmp function mpf_sparse.c

Then, using non uniform derivative (on same uniform grid):
source code pde2sin_nueq_gmp.c
output pde2sin_nueq_gmp.out
source code pde2sin8_nueq_gmp.c
output pde2sin8_nueq_gmp.out
source code pde4sin_nueqsp.c
output pde4sin_nueqsp.out
gmp function mpf_nuderiv.h
gmp function mpf_nuderiv.c
gmp function mpf_inverse.h
gmp function mpf_inverse.c
gmp function mpf_simeq.h
gmp function mpf_simeq.c
gmp function mpf_sparse.h
gmp function mpf_sparse.c

binary.txt convert binary to decimal

multiple precision HW4 is assigned

    <- previous    index    next ->

Other links

Go to top