Third party libraries

The final week is spent studying the use of third party libraries; development packages of headers and libraries external to that provided by the compiler and C library cores.

The g2 library

The g2 library is a small, fast and lightweight graphics library that offers a pretty large subset of the functionality of libplot, in a very small package. Its main advantage is that it has C, FORTRAN77 and perl front-ends.
The program below solves the Newtonian equations of motion for a damped oscillator;

v(t) = x(t-+--dt)----x(t)-, x(t + dt) = x(t) + dt ⋅ v(t) dt
and
v(t-+--dt)----v(t) = - x(t) - s⋅ v(t), v(t + dt) = v(t) - x(t) ⋅dt - s ⋅v(t) ⋅dt dt
using the simplistic Euler extrapolation algorithm, and animates the output in an X window. g2 lacks double-buffering, and so there is screen flicker, but a smarter time loop containing fewer calculations can minimize this.
In order to use g2, its header must be included, the library must be linked at compile time, and you must open both a virtual device and attach a display device to it.
#include <stdio.h>  
#include <g2.h>  
#include <g2_X11.h>  
#include <unistd.h> /* so I can use pause() */  
 
int main()  
{  
    int k, dev, dev1;  
    float dt, x, x0, v, v0, s;  
    char str[256];  
    dt=0.003;  
    x0 = 1.0;  
    v0 = 0.0;  
    s=0.15;  
    dev=g2_open_vd(); /* open virtual device */  
    dev1=g2_open_X11(300,300);  
    g2_attach(dev, dev1);  
    /*  x and y-axes */  
    sprintf(str, "%s", "v");  
    g2_line(dev1, 150.0, 0.0, 150.0, 300.0);  
    g2_string(dev1, 155.0, 290.0, str);  
    sprintf(str, "%s", "x");  
    g2_line(dev1, 0.0, 150.0, 300.0, 150.0);  
    g2_string(dev1,290.0, 140.0, str);  
    for(k=1;k<6000;k++){  
x=x0+v0*dt;  
v=v0-x0*dt-s*v0*dt;  
g2_line(dev1, 150.0+100.0*x0, 150.0+100.0*v0, 150.0+100.0*x,\  
                                                      150.0+100.0*v);  
x0=x;  
v0=v;  
    }  
    g2_flush(dev1);  
    pause(); /* send interrupt Crtl+c to kill */  
    g2_close(dev1);  
    return(0);  
}  

Once again we want to use a Makefile to build the project. The corresponding Makefile for a C project is given below.

# Simple C Makefile  
 
CC = gcc  
CPP = /lib/cpp  
 
CFLAGS = -O2  
LDFLAGS = -I/usr/local/include -L/usr/X11R6/lib -L/usr/local/lib -lm -lX11 \  
                                                -lgd -lg2  
all:    phase  
 
phase:  phase.c  
        $(CC) $(C_FLAGS) phase.c $(LDFLAGS) -o phase  
clean:  
        -rm phase

Translating g2 and libplot function calls

The table below may help you translate some function calls.

g2 libplot
dev=g2_open_vd() pl_openpl ()
dev1=g2_open_X11(Width,Height) pl_handle = pl_newpl (”X”, stdin, stdout, stderr)
pl_fspace (LLx, LLy, URx, URy)
g2_attach(dev, dev1) pl_selectpl (pl_handle)
g2_line(dev1, 150.0, 0.0, 150.0, 300.0) pl_fline(150.0, 0.0, 150.0, 300.0)
g2_string(dev1,290.0, 140.0, str) pl_fmove(290.0, 140.0)
pl_label(str)
g2_flush(dev1) pl_flushpl()
g2_close(dev1) pl_deletepl()
g2_rectangle(dev,20,20,150,150) pl_fbox(20.0,20.0,150.0,150.0)
g2_clear(dev1) pl_erase()
g2_filled_rectangle(dev, x1,y1,x2,y2) pl_filltype(1)
pl_fbox(x1,y1,x2,y2)
g2_circle(dev1, x, y, r) pl_fcircle(x,y,r)
c=g2_ink(dev, R,G,B), g2_pen(dev, c)pl_pencolor(R,G,B)
g2_set_font_size(dev, size) pl_ffontsize(size)

Both libraries have many more functions than this, it really is necessary to read the documentation. G2 can also draw triangles and polygons.
My favorite feature of g2 is its ability to grab and read the location of the mouse and the state of its buttons (0 if not pressed, 256 if left-clicked, 1024 if right-clicked) with g2_query_pointer(dev, *x, *y, *button). This can be used to stop, restart or re-initialize animations. Read the g2 documentation to learn more.
One note on colors, g2 has double color bits 0.0 ≤ R,G,B ≤ 1.0, while libplot uses int, 0 ≤ R,G,B ≤ 65535.

An important difference between libplot and g2 is the way in which the graphics window is set up. Libplot uses affine coordinates, and so the user does not need to translate object coordinates into pixels in order to place them in the indow. The g2 function g2_open_X11(int width, int height) opens a window width pixels wide and height pixels high, with no coordinate system. The user must translate particle coordinated into pixel positions in g2 programming.

Another note, g2 has multiple header files g2.h, g2_PS.h g2_GIF.h and the appropriate ones need to be included in your program preamble.

FORTRAN77 example

The simple program below illustrates the basic usage of g2 in FORTRAN77 to create an on-screen figure.
It does so with a do loop. Study the loop syntax carefully; it uses a line label, here the label is the number 100, and instructions to carry out the calculations between the do 100 line and the 100 continue line, while k ≤ 6000. The k variable indexes the loop. Notice that
1. any line beginning with a c or C is a comment,
2. that all instruction lines must begin with a tab,
3. lines can be no longer than 80 characters,
4. there are restrictions on variable names, both in length and in which letter the name can start with,
5. external routines such as functions or subroutines or g2 commands are invoked with call.
Like a C program, we first define our variables. We next initialize some to particular values, and we open up a g2 graphical display device and name it so it can be referred to. The g2 commands are very simple and intuitive; they need to know what to draw and what device to draw it on.

C     Phase space program for oscillator  
C     jrs  
      real phase  
      real x, v, t, x0, v0,a,b,c,d  
      real dt, dev, dev1, q, s  
      character str*256  
      dt = 0.003  
      x0 = 1.0  
      v0 = 0.0  
      s=0.15  
c set up a graphics device  
      dev=g2_open_vd()  
      write (6,*) d  
      dev1=g2_open_X11(300.0, 300.0)  
      write (6,*) d1  
      call g2_attach(dev, dev1)  
c     x and y-axes  
      str="v"//char(0)  
      call g2_line(dev1, 150.0, 0.0, 150.0, 300.0)  
      call g2_string(dev1, 155.0, 290.0, str)  
      str="x"//char(0)  
      call g2_line(dev1, 0.0, 150.0, 300.0, 150.0)  
      call g2_string(dev1,290.0, 140.0, str)  
      do  100 k=1, 6000  
c     Simple Euler integration  
 
         x=x0+v0*dt  
         v=v0-x0*dt-s*v0*dt  
c         call g2_circle(d1, 150.0+100.0*x, 150.0+100.0*v,1.0)  
         a=150.0+100.0*x0  
         b=150.0+100.0*v0  
         c=150.0+100.0*x  
         d=150.0+100.0*v  
         call g2_line(dev1, a,b,c,d)  
c  
         x0=x  
         v0=v  
 
 100     continue  
      call g2_flush(dev1)  
 
c     halt until someone hits a key  
      read (*,*) q  
      stop  
      end

Copy this program and the Makefile into your home directory and build the project. Run the program. Try some modifications to the code and see what you get. Some things to note; the lines

      dev=g2_open_vd()  
      write (6,*) d  
      dev1=g2_open_X11(300.0, 300.0)  
      write (6,*) d1  
      call g2_attach(dev, dev1)

set up a graphics device and attach a 300 × 300 pixel XWindow to it, with the origin in the lower corner. The subsequent graphics calls draw lines in this window, we must relocate the origin to its center and rescale our figure accordingly.

Another FORTRAN77 example

The code below performs a graphical simulation of Conway’s game of life.

c Simple Conway’s game of life  
c FORTRAN77 version  
c use the makefile to build  
c jrs 10.21.2002  
      parameter (N=80)  
      dimension cell(80,80), cell0(80,80)  
      real x,y,xb,yb,q  
 
c initialize using built-in random number function  
 
      do 200 I=1, N  
         do 100 J=1, N  
            if(rand().gt.0.5) then  
               cell0(I,J)=1  
               else  
                  cell0(I,J)=0  
                  end if  
 100              continue  
 200              continue  
 
c set up a graphics device  
 
      dev=g2_open_vd()  
      write (6,*) dev  
      dev1=g2_open_X11(400.0, 400.0)  
      write (6,*) dev1  
      call g2_attach(dev, dev1)  
 
c evolution step first loop is time loop  
 
      do 900 IT=0, 100  
 
c set up correct boundary conditions (that means YOU)  
c I will update the interior cells...  
 
      do 400 I=2, N-1  
         do 300 J=2, N-1  
            IS=cell0(I-1,J+1)+cell0(I, J+1)+cell0(I+1, J+1)  
            IS=IS+cell0(I-1,J)+cell(I+1,J)  
            IS=IS+cell0(I-1,J-1)+cell0(I, J-1)+cell0(I+1, J-1)  
            if((cell0(I,J).eq.0).and.(IS.eq.3)) then  
               cell(I,J)=1  
            else if ((cell0(I,J).eq.1).and.((IS.eq.2).or.(IS.eq.3)))then  
               cell(I,J)=1  
            else  
               cell(I,J)=0  
            end if  
 300        continue  
 400        continue  
 
c then for example... top  
            do 1000 J=2, N-1  
               IS=cell0(J-1,N-1)+cell0(J,N-1)+cell0(J+1,N-1)  
               IS=IS+cell0(J-1,1)+cell0(J,1)+cell0(J+1,1)  
               IS=IS+cell0(J-1,N)+cell0(J+1,N)  
               if((cell0(J,N).eq.0).and.(IS.eq.3)) then  
               cell(J,N)=1  
            else if ((cell0(J,N).eq.1).and.((IS.eq.2).or.(IS.eq.3)))then  
               cell(J,N)=1  
            else  
               cell(J,N)=0  
            end if  
 1000       continue  
 
 
c ... the left wall  
            do 2000 J=2, N-1  
               IS=cell0(2,J-1)+cell0(2, J)+cell0(2, J+1)  
               IS=cell0(N,J-1)+cell0(N, J)+cell0(N, J+1)  
               IS=IS+cell0(1,J-1)+cell0(1, J+1)  
               if((cell0(1,J).eq.0).and.(IS.eq.3)) then  
               cell(1,J)=1  
            else if ((cell0(1,J).eq.1).and.((IS.eq.2).or.(IS.eq.3)))then  
               cell(1,J)=1  
            else  
               cell(1,J)=0  
            end if  
 2000       continue  
 
c right wall  
            do 3000 J=2, N-1  
               IS=cell0(N-1,J-1)+cell0(N-1, J)+cell0(N-1, J+1)  
               IS=cell0(1,J-1)+cell0(1, J)+cell0(1, J+1)  
               IS=IS+cell0(N,J-1)+cell0(N, J+1)  
               if((cell0(N,J).eq.0).and.(IS.eq.3)) then  
               cell(N,J)=1  
            else if ((cell0(N,J).eq.1).and.((IS.eq.2).or.(IS.eq.3)))then  
               cell(N,J)=1  
            else  
               cell(N,J)=0  
            end if  
 3000       continue  
 
c then for example... bottom  
            do 4000 J=2, N-1  
               IS=cell0(J-1,2)+cell0(J,2)+cell0(J+1,2)  
               IS=IS+cell0(J-1,N)+cell0(J,N)+cell0(J+1,N)  
               IS=IS+cell0(J-1,1)+cell0(J+1,1)  
               if((cell0(J,1).eq.0).and.(IS.eq.3)) then  
               cell(J,1)=1  
            else if ((cell0(J,1).eq.1).and.((IS.eq.2).or.(IS.eq.3)))then  
               cell(J,1)=1  
            else  
               cell(J,1)=0  
            end if  
 4000       continue  
c we SHOULD update the four corners, bugger the corners!  
c drawing portion  
 
      do 600 I=1, N  
         do 500 J=1, N  
 
                  x=5.0*float(I)-5.0  
                  y=5.0*float(J)-5.0  
                  xb=x+5.0  
                  yb=y+5.0  
                  call g2_pen(dev1, cell(I,J))  
                  call g2_filled_rectangle(dev1, x,y,xb,yb)  
 
 500              continue  
 600              continue  
 
c resetting portion  
 
      do 800 I=1, N  
         do 700 J=1, N  
            cell0(I,J)=cell(I,J)  
 700        continue  
 800        continue  
 
 900        continue  
 
c time loop ends, now wait for key input  
            read (*,*) q  
            call g2_close(dev)  
            stop  
      end  
 
 
 
 
 

Problems

1. Read through the documentation for g2 and write a program in your choice of language that draws the solar system, with a yellow sun, one elliptical planetary orbit, a blue planet, with axis labels, and name labels for the sun and planet.

2. Animation. The g2 library is not designed for animation, since it lacks double-buffering, but if one does not mind the flickering, it is possible to do animation with g2. To pull it off, you must use g2_clear( device) inside of your loop, draw everything, and flush the graphics buffer with g2_flush(device) within your time loop. Try it out by animating your oscillator in either FORTRAN77 or C. If the animation is too fast, you may need to use some sort of pause statement to slow the loop down. In C we would use

#include <sys/types.h>  
#include <time.h>  
#include <unistd.h>  
...  
/* loop body begins */  
usleep(0);  
...  
/* loop body ends */

Rather than clearing the screen with each cycle, try whiting out the objects that need to be redrawn. This will speed things up considerably, erasing the screen altogether is costly. Use g2_pen(dev, color) to set the color to white, then use g2_filled_rectangle(dev, x1, y1, x2, y2) to overdraw a white box on top of the object you want redrawn or moved.

Very high precision computing

On 32 bit targets (Intel) or even 64 bit machines (Sparc, Alpha, IA64) the representation of double precision numbers may contain too few digits for reliable computations involving summations of large numbers of terms, since this will result in accumulation errors far to the left of what is supposed to be the least reliable digit. In these cases, or in cases in which for example we need 40-500000 digits in a computed number, we can use the GNU multiple precision library gmp-4.1.1. In fact the only limit on th number of decimal places the library is capable of providing are imposed by your hardware.
Gmp is a collection of high-level functions for allocating, initialising, and processing very high precision numbers as new string types, mpz_t, mpf_t, mpq_t and so forth. The library is fairly easy to use, and we will now run over the important aspects of how it is used in a program.
First the library is compiled and installed somewhere in your directory tree. Compilation creates a static library, libgmp.a and a header gmp.h, which should be copied to /usr/local/lib and /usr/local/include respectively. To compile a program using the library functions, its preamble must include the header gmp.h and the compile line should be for example

gcc -o program program.c -I/usr/local/include -lgmp -lm

Declaring GMP types

Variable types provided by the library are integer (mpz_t), rational (mpq_t), and floating point (mpf_t). These can be integrated into ordinary variables, pointers or structures

mpz_t n,m;   /* integers */  
mpz_t *NUM;  /* pointer to integer */  
mpf_t V[N];  /* array of floats */  
struct FT { mpf_t val;  
            mpz_t freq;  
          };

Initialisation of GMP types

All of these variable types must be initialised before they can be used. There are many initialisation functions, as well as combination declaration/initialisations. For example each of the following functions initialises its target to the appropriate mp-type representation of zero

mpz_t I;    /* declare */  
mpz_init(I);/* initialise */  
mpf_t F;  
mpf_init(F);  
mpf_t F2;  
mpf_init2(F2, 256); /* initialise F2 to at least 256 bits */  
mpz_clear(I); /* deallocate I */  
mpf_clear(F2); /* deallocate F2 */

An example of the use of a combination initialise-assign operation is the following

mpz_t pie;  
mpz_init_set_str (pie, "3141592653589793238462643383279502884", 10);

The syntax used in the description of all mp functions outlined in the gmp-info-X files is the following;

void mpz_XXX(mpz_t ROP, mpz_t OP);

in which ROP is a “return” type for the operation. It is not possible to declare functions that actually return an mpf_t type, or mpz_t type. Instead one creates a void function that passes the result back into the calling program by placing its value in a pointer. For example

void square(mpf_t v2,mpf_t v); /* illustrates void mpf_t function */  
void square(mpf_t v2,mpf_t v){  
 mpf_t tmp;  
 mpf_init2(tmp,128);  
 mpf_init_set(tmp,v);  
 mpf_mul(v2, tmp, tmp);  
 }

Assignment of GMP types

There are several useful assignment operators, depending on the desired precision of a data variable. I have found that if you want to retain a very large number of digits in a float, assign it a value with mpf_set_str. Other integer assignments are

 void mpz_init_set (mpz_t ROP, mpz_t OP);  
 void mpz_init_set_ui (mpz_t ROP, unsigned long int OP);  
 void mpz_init_set_si (mpz_t ROP, signed long int OP);  
 void mpz_init_set_d (mpz_t ROP, double OP);  
 int mpz_init_set_str (mpz_t ROP, char *STR, int BASE);

and for floating point types

void mpf_set (mpf_t ROP, mpf_t OP);  
void mpf_set_ui (mpf_t ROP, unsigned long int OP);  
void mpf_set_si (mpf_t ROP, signed long int OP);  
void mpf_set_d (mpf_t ROP, double OP);  
void mpf_set_z (mpf_t ROP, mpz_t OP);  
void mpf_set_q (mpf_t ROP, mpq_t OP);

Examples of these are the following

mpz_t x,y,z,t;  
mpz_init2(x,256);  
mpz_init2(y,256);  
mpz_init2(z,256);  
mpf_set_str(x,’’0.577215664901532860606512’’, 10);  
mpf_set(y,x);  
mpf_set_d(z, log(2.5)); /* this will limit accuracy on 32 bit targets */

GMP arithmetic functions

We begin with a simple double to gmp conversion function. This is useful for feeding computations back into parts of the program that are not mp-aware. For example

 double mpf_get_d (mpf_t OP);

transforms an mpf_t OP back into a double precision number.

Floating point arithmetic
The most useful arithmetic functions are those listed below

void mpf_add (mpf_t ROP, mpf_t OP1, mpf_t OP2); /* Set ROP to OP1 + OP2 */  
void mpf_sub (mpf_t ROP, mpf_t OP1, mpf_t OP2); /* Set ROP to OP1 - OP2 */  
void mpf_mul (mpf_t ROP, mpf_t OP1, mpf_t OP2); /* Set ROP to OP1 * OP2 */  
void mpf_div (mpf_t ROP, mpf_t OP1, mpf_t OP2); /* Set ROP to OP1 / OP2 */  
void mpf_sqrt (mpf_t ROP, mpf_t OP);            /* Set ROP to sqrt (OP) */  
void mpf_neg (mpf_t ROP, mpf_t OP);             /* Set ROP to -OP       */  
void mpf_abs (mpf_t ROP, mpf_t OP);             /* Set ROP to | OP |    */

Comparison and IO functions

We will only list the floating point versions of these. Read through the info pages accompanying the software for information pertaining to the use of the integer and rational functions.

int mpf_cmp (mpf_t OP1, mpf_t OP2); /* Return a positive value if OP1 > OP2,  
 zero if OP1 = OP2, and a negative value if OP1 < OP2 */  
 
int mpf_eq (mpf_t OP1, mpf_t OP2, unsigned long int op3); /*  Return non-zero  
 if the first OP3 bits of OP1 and OP2 are equal, zero otherwise */  
 
size_t mpf_out_str (FILE *STREAM, int BASE, size_t N_DIGITS, mpf_t OP);  
 /* outputs OP as a string with N_DIGITS digits in base BASE */  
 
size_t mpf_inp_str (mpf_t ROP, FILE *STREAM, int BASE); /* inputs a string  
 from a stream pointer and stores it in OP in base BASE. returns number of  
 bytes read */

An example; computing e

Computation of

e ≈ 1 + 1--+ 1--+ 1--+ 1-+ ⋅⋅ ⋅ + -1-- 1! 2! 3! 4! N !
using gmp is a relatively short application.
#include<stdio.h>  
#include<stdlib.h>  
#include<gmp.h>  
#include<math.h>  
int n,N;  
mpf_t one, tmp, term, sum;  
 
main(int argc, char *argv[]){  
    mpf_init2(one,1000);  
    mpf_init2(tmp,1000);  
    mpf_init2(term,1000);  
    mpf_init2(sum,1000);  
    mpf_set_ui(one,1);  
    mpf_set_ui(sum,1);  
    mpf_set_ui(term,1);  
    N=atoi(argv[1]);  
    for(n=2;n<N;n++){  
mpf_set_ui(tmp,n);     /* tmp=n */  
mpf_div(term, term, tmp);    /* term =1/n! */  
mpf_add(sum,sum,term); /* sum= e */  
    }  
    mpf_add(sum,sum,one);  
    mpf_out_str(stdout,10,300, sum);  
    printf("\n");  
 
    mpf_clear(one);  
    mpf_clear(tmp);  
    mpf_clear(term);  
    mpf_clear(sum);  
}  

Note that we declare all types, then initialize them to unsigned integers, work with them, and then de-allocate them. The program accepts N from the command line.

The gmp library is fast; the run time for this program with N = 4000 is under 0.5 sec on a 200 MHz computer.

An example; partitions

Partitions of integers n are defined as sums of positive integers less than n that add up to n, for example

4 = 4, 1 + 3, 1 + 1 + 2, 1 + 1 + 1 + 1, 2 + 2
5 = 5, 1 + 4, 1 + 1 + 3, 1 + 1 + 1 + 2, 1 + 1 + 1 + 1 + 1, 1 + 2 + 2, 2 + 3
and we define p(n) to be the number of partitions of n, so
p(4) = 5, p(5) = 7
Let us define r(n,m) to be the number of partitions of n with one part equal to m, and all other parts equal or greater than m. Then
p(n) = r(n, 1) + r(n, 2) + ⋅⋅ ⋅r(n, n)
For example
r(5, 1) = 5, {1 + 4, 1 + 1 + 3, 1 + 1 + 1 + 2, 1 + 1 + 1 + 1 + 1, 1 + 2 + 2}
r(5, 2) = 1, {2 + 3}
A theorem in number theory is that
r(n, m) = r(n - m, m)+r(n - m, m+1)+ ⋅⋅⋅+r(n - m, n- m), for n - m > m
and it is easy to show that
r(m + k,m) = 0, 0 < k < m
and
r(n, m) = 0 if n < m, r(n, n) = 1
Another very powerful number theoretic result is
∑ k 3k2----k- 3k2-+--k- p(n) = (- 1) (p(n - ) + p(n - )) k=1 2 2
in which we make the further definition p(0) = 1, p(n < 0) = 0. This is not difficult to prove from the Jacobi triple product formula (from the theory of elliptic functions).
One of the greatest triumphs of twentieth century mathematics was the Ramanujan-Hardy-Rademacher formula for the partition of an integer. A simple approximation for the number of partitions derived from this formula is
1 π√ 2n- p(n) → ---√--e 3 4n 3
The GMP program below uses the Jacobi recursion to compute the partition. Try to compare the results to the values gotten from the RHR formula.
#include<stdio.h>  
#include<gmp.h>  
 
long i,j,n,k;  
double d;  
mpz_t S;  
 
void p(long n, mpz_t R);  
 
main(){  
mpz_init(S);  
for(j=2;j<1001;j++){  
p(j,S);  
/*d=mpz_get_d(S);*/  
printf("%d\t",j);  
mpz_out_str(stdout,10,S);  
printf("\n");  
}  
}  
 
void p(long n, mpz_t R)  
{  
 unsigned long l, k,k1,k2;  
mpz_t sgn,sum,tmp,P[n+1];  
for(i=0;i<=n;i++)  
mpz_init(P[i]);  
mpz_init(tmp);  
mpz_init(sum);  
mpz_init(sgn);  
mpz_init_set_ui(P[0],1);  
mpz_init_set_ui(P[1],1);  
mpz_init_set_ui(P[2],2);  
 
 
 for(l=3;l<=n;l++){  
mpz_set_si(sgn,1);  
mpz_set_ui(sum,0);  
   for(k=1;k<=l;k++){  
k1=(3*k*k-k)/2;  
k2=(3*k*k+k)/2;  
if(k1>l)break;  
 else{  
mpz_mul(tmp,P[l-k1],sgn);  
mpz_add(sum,sum,tmp);}  
 if(k2<=l){  
 
mpz_mul(tmp,P[l-k2],sgn);  
mpz_add(sum,sum,tmp);}  
mpz_neg(sgn,sgn);      /* negate the sign */  
}  
mpz_set(P[l],sum);}  
mpz_set(R,P[n]);  
for(i=0;i<=n;i++)  
mpz_clear(P[i]);  
}  
 

for example

n p(n)
1999023036905895619102711927099838059623708152962831993318457458826854491
527110736267250824578187297578414046706043303911639959133346510497979
5374108729994206
1999123245666796169459664402354595717958747780837377246706273194844243524
38991086685230920653300533234255126724400470940647848098209024483153436
37071762908642

The best way to learn how to use these functions is to actually try using them in simple programs. In the following section we will approximate functions on finite domains with Chebyshev series, and we will write gmp-programs to construct the coefficients for the series to very high precision.
One final note; in any program that repeatedly calls GMP-based routines, it is important that those routines deallocate or free the memory used to store any GMP types created and used by the routine with mpf_clear(mpf_t OP), mpz_clear(mpz_t OP) and mpq_clear(mpq_t OP). If you do not do this, you will wish that you had. Trust me on this one.

The gmp library uses CPU-specific instructions. It really should be recompiled natively rather than getting it precompiled from someone else, so that your CPU-type can be detected and the correct instructions be compiled in. If you don’t do this, your operations could result in ILLEGAL INSTRUCTION.

Problems

3. Write a program to compute factorials up to 1000! and write them to the screen as they are created. This is as simple as you can get as far as use of the GMP library is concerned. The GMP library has a function call that computes factorials; do not use it, compute them from

n! = n ⋅ (n - 1)!

4. Modify the sample program for e and use it to compute ex to 200 decimal places using

∑∞ 1 ex = ---xn n! n=0
You should take x from the command line and initialize a gmp-type to the double-precision value of x.

5. S. Ramanujan has produced some astounding expressions for π, that converge very rapidly compared to any ordinary series or iterative expansion. One of the simplest of his results is

√ -- 1 8 ∑∞ (4n)! (1103 + 26390 n) --= ------ -------------------------- π 9801 (n! (396)n)4 n=0
Write a program using the GMP library to compute π to 4096 bit accuracy, or 200 decimal places, whichever comes first. You will need to consult a high precision math table or Mathematica to verify your computation.

6.Use the Borwein formula

∑∞ π = --1---(---4----- ---2---- - ---1----- ---1----) (16)k 8k + 1 8k + 4 8k + 5 8k + 6 k=0
to compute π to 10, 000 decimal places. Is this a more efficient algorithm than the Ramanujan formula? This is well designed for computing digits of π in hexadecimal, it has been used to get the 16 billionth digit.

A short GMP code library

Miller algorithm for Modified Bessel function

These functions are used as coefficients in many Chebyshev accelerated series, and appear in boundary value problems in radiation theory and classical electromagnetism.

#include<stdio.h>  
#include<math.h>  
#include<stdio.h>  
#include<gmp.h>  
 
 
int n,m;  
int MAX=36;  
int TOP, PREC;  
int N;  
mpf_t x;          /* argument of the function */  
mpf_t IN2, IN1,IN0,arg,tmp,y,two;  
mpf_t sum;  
int parity;  
 
main(){  
 
mpf_init2(IN0,256);  
mpf_init2(IN1,256);  
mpf_init2(IN2,256);  
mpf_init2(tmp,256);  
mpf_init2(sum,256);  
mpf_init2(arg,256);  
mpf_init2(two,256);  
mpf_init2(y,256);  
mpf_init2(x,256);  
mpf_set_ui(two,2);  
mpf_set_ui(x,1);  
   /* we will get J_8(x) */  
mpf_set_ui(sum,0);  
mpf_set_ui(IN0,0);  
mpf_set_ui(IN1,1);  
mpf_div(arg,two,x);  
 for(N=1;N<=30;N++){  
parity=0;  
 
/* begin at the top of the recursion at an arbitrary value */  
PREC=floor((sqrt(MAX*(N+2))));  
TOP=2*(N+PREC);  
 for(n=TOP;n>=1;n--){  
mpf_set_ui(tmp,n);  
 
mpf_mul(tmp,tmp,arg);  
 
mpf_mul(IN2,tmp,IN1);  
mpf_add(IN2,IN2,IN0);  
 
mpf_set(IN0,IN1);  
mpf_set(IN1,IN2);  
if(((n-1)/4)*4==n-1) mpf_add(sum,sum,IN1);  
if(((n+1)/4)*4==n+1) mpf_sub(sum,sum,IN1);  
/* only add up the even bessels with the appropriate sign*/  
parity=1-parity;  
if(n==N)mpf_set(y,IN0);}  
mpf_add(sum,sum,sum);  
 mpf_neg(sum,sum);      /* negate sum */  
mpf_add(sum,sum,IN1);  
 mpf_neg(y,y);           /* we have computed -I_n */  
mpf_div(y,y,sum);  
printf("I_%d(",N);  
mpf_out_str(stdout,10,10,x);  
printf(")=\t");  
mpf_out_str(stdout,10,40,y);  
 printf("\n");}  
}  

Legendre polynomials

An example of a callable subroutine. You cannot create functions that return GMP types, but you can create void punctions that write GMP types into locations pointed to by pointers to GMP types.

 
void P(int n,  mpf_t in, mpf_t out)  
{/* high precision Legendre polynomials out = P_n(in)  */  
  /* stores P_i(in) for 0<=i<=n for speed              */  
int i;  
mpf_t tmp1,tmp2,coeff1,coeff2,one,two,index;  
mpf_t legendre[M+1];  
mpf_init2(index,256);  
mpf_init2(one,256);  
mpf_init2(tmp1,256);  
mpf_init2(two,256);  
mpf_init2(tmp2,256);  
mpf_init2(coeff1,256);  
mpf_init2(coeff2,256);  
for(i=0;i<=M;i++)  
mpf_init2(legendre[i],256);  
mpf_set_ui(one,1);  
mpf_set_ui(two,2);  
mpf_set(legendre[0],one);  
mpf_set(legendre[1],in);  
 
  for(i=2;i<=n;i++){  
mpf_set_ui(index,i);  
mpf_mul(coeff1,two,index);  
mpf_sub(coeff1,coeff1,one);  
mpf_div(coeff1, coeff1, index);  
mpf_sub(coeff2,index,one);  
mpf_div(coeff2,coeff2,index);  
mpf_set(tmp1,legendre[i-1]);  
mpf_set(tmp2,legendre[i-2]);  
mpf_mul(tmp1,tmp1,in);  
mpf_mul(tmp1, tmp1,coeff1);  
mpf_mul(tmp2,tmp2,coeff2);  
 mpf_sub(legendre[i],tmp1,tmp2);}  
  mpf_set(out,legendre[n]);  
     /* clean up */  
 mpf_clear(tmp1);        /* serious problems are caused */  
 mpf_clear(tmp2);         /* if you don’t dealloc these */  
 mpf_clear(coeff1);        /* types */  
mpf_clear(coeff2);  
mpf_clear(one);  
mpf_clear(two);  
mpf_clear(index);  
for(i=0;i<M;i++)  
mpf_clear(legendre[i]);  
 
}  

High Precision Computation of π

In subsequent programs it may be very useful to have π to high precision. The following program is inefficient, but computes π with unlimited precision. To 100 decimal places

π = 3.14159265358979323846264338327950288419716939937510582097494
4592307816406286208998628034825342117068
#include<stdio.h>  
#include<stdlib.h>  
#include<math.h>  
#include<gmp.h>  
mpf_t index,eight,c1,c2,root8,denom,threeninesix,tmp,sum,factorial;  
mpf_t factorial4,denom2,ONE;  
int n,m;  
main(){  
mpf_init2(index,2048);  
mpf_init2(eight,2048);  
mpf_init2(c1,2048);  
mpf_init2(c2,2048);  
mpf_init2(root8,2048);  
mpf_init2(denom,2048);  
mpf_init2(denom2,2048);  
mpf_init2(threeninesix,2048);  
mpf_init2(tmp,2048);  
mpf_init2(sum,2048);  
mpf_init2(factorial,2048);  
mpf_init2(factorial4,2048);  
mpf_init2(ONE,2048);  
 mpf_set_ui(eight,8);       /* set constant factors */  
mpf_set_ui(c1,1103);  
mpf_set_ui(c2,26390);  
mpf_set_ui(threeninesix,396);  
mpf_set_ui(denom,9801);  
mpf_mul(threeninesix,threeninesix,threeninesix); /* square */  
mpf_mul(threeninesix,threeninesix,threeninesix); /* (396)ˆ4*/  
 mpf_sqrt(eight,eight); /* compute sqrt(8) */  
mpf_set(denom2,threeninesix);  
mpf_set_ui(ONE,1);  
/* start adding up terms */  
mpf_set(sum,c1);  
 
/* loop, build the factorials */  
 for(n=1;n<19;n++){  
mpf_set_ui(factorial, 1);  
mpf_set_ui(factorial4, 1); /* inefficient, but so what? */  
mpf_set_ui(index,n);  
 for(m=1;m<=n;m++){  
mpf_set_ui(tmp,m);  
 mpf_mul(factorial, factorial,tmp);}  
mpf_out_str(stdout,10,20,factorial);  
printf("\t");  
 mpf_mul(factorial,factorial,factorial); /*square */  
 mpf_mul(factorial,factorial,factorial); /* (n!)ˆ4 */  
 for(m=1;m<=4*n;m++){  
mpf_set_ui(tmp,m);  
 mpf_mul(factorial4, factorial4,tmp);}  
mpf_out_str(stdout,10,20,factorial4);  
printf("\n");  
 mpf_mul(tmp,c2,index);  /* 26390 n */  
mpf_add(tmp,tmp,c1);      /* 1103 + 26390 n */  
mpf_mul(tmp,tmp,factorial4);  
mpf_div(tmp,tmp,factorial);  
mpf_div(tmp,tmp,denom2);  
mpf_add(sum,sum,tmp);  
 mpf_mul(denom2,denom2,threeninesix); /* make (396)ˆ4n */  
 
 } /* finish the sum */  
mpf_mul(sum,sum,eight);  
mpf_div(sum,sum,denom);  
mpf_div(sum,ONE,sum);  
mpf_out_str(stdout,10,100,sum);  
printf("\n");  
}