<- previous index next ->
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
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
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
mpmath_example.py program
mpmath_example_py.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
<- previous index next ->