Elemetary problems using GMP

[Thomas Chassaing]

Hi,

I try to write a program to calculate the energy of the Helium atom 
variationally. For this I need arbitrary precision numbers (quadruple 
precision at least). I wanted to give gmp a try but couldn't even get 
the simplest test program to work. I hope somebody in this group can 
tell me what I'm doing wrong (I read the manual, not very carefully though).

Here's the code:

<<<<<<<<<<<<<<<<<<<<<<<<<<<<<

#include <gmp.h>

#define np 2
#define nbf np*np+1

const double L   = 5.0;

int main()
{
   mpf_t a[nbf], b[nbf], c[nbf], H[nbf][nbf], N[nbf][nbf];
   mpf_t tst;
   int i,j,k,l;
   double d;

   //mpf_set_default_prec(256);

   mpf_init_set_d(tst,(double)(5.0));
   gmp_printf("%20.15F\n",tst);

   for (i=0;i<nbf;i++) {
       mpf_init(a[i]);
       mpf_init(b[i]);
           mpf_init(c[i]);
   }

   mpf_set_d(a[0],1.858058);
   mpf_set_d(b[0],0.0);
   mpf_set_d(c[0],1.0);
   gmp_printf("%4d: %.10f  %F  %F\n", l, a[0], b[0], c[0]);
   l = 1;
   for (i=0; i<np; i++) {
         for (j=0; j<np; j++) {
           d = L*i/(np+1);
           mpf_set_d(a[l],d);
           mpf_set_d(b[l],d);
       mpf_set_d(c[l],0.0);
           gmp_printf("%4d: %F  %F  %F\n", l, a[l], b[l], c[l]);
           l++;
         }
   }
   return 0;
}

 >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>

And here's what I get:

    5.000002861018331
    0: -1.9958961010  0.000000  0.000000
    1: -1.995911  0.000000  0.000000
    2: -1.995927  0.000000  0.000000
    3: -1.995942  0.000000  0.000000
    4: -1.995957  0.000000  0.000000

Already it seems strange that the 5.0 is not exactly represented (since 
this number can be represented exactly with double).

Thanks in advance
T. Chassaing

-- 
Thomas Chassaing
Physical Chemistry Institute
University of Zurich
Winterthurerstrasse 190
CH-8057 Zurich, Switzerland

Phone : ++41 1 635 4490		
E-mail: tchassai at pci.unizh.ch