Iniatialization and getting unexpected values

[Jarvis Loh]

Hi,

There are 2 problems I faced.

(1) Initialization:
const double SInter=0.0000000000000001; (in other words 1E-16)

             (This is actually part of a code on molecular dynamics simulation.)

When I output "SInter" when debugging the code, it gives a value of 0.000000000000000099999999999... (9.9999E-17).
Why is this so? I have also tried SInter=1E-16 before the decimal form, but the same problem occurs.

Is there a way I can input an exact "1E-16"?

(2) Accumulation:

void CollisionTimes() {

 double Prob2[31];                                                          
 Prob2[0]=V;
 double CollisionT;

 for (int i=0;i<N;i++) {
     for (int p=1;p<31;p++) {
         double ct=p*((5E-12)/step);
         Prob2[p]=Prob2[p-1]+(V*exp(-V*ct));
     }

     for (int c=0;c<5000;c++) {

         Randomize=(rand()%11);
         double Randomize_2=Randomize/10.0;
         Graphpos=Randomize_2*Prob2[30];

         for (int j=0;j<30;j++) {
             if ((Graphpos>=Prob2[j]) && (Graphpos<Prob2[j+1])) {
                CollisionT2[i][c]=(j*((5E-12)/step))+((5E-12)/(step*2));           LINE *               
             }
         }
     }
 }

}

[Under another subroutine, i have the following.]
initialize();

 CollisionTimes();                                                                   LINE @
 Actualstep=0;
 counter=0;
 Num=0;
 dt=SInter;
 double CollisionT3=0;
 for (int stepno=0;stepno<timestep;stepno++) { 
     double tottime=(Actualstep)*SInter;// has to be always 10E-16, that is why dt is not used

     AtomNo=0;

         for (int next=1;next<N;next++) {
             int i=AtomNo;
             if (CollisionT2[i][Num]<CollisionT2[next][Num]) {
                AtomNo=i;
             }
             else {
                AtomNo=next;
             }
         }
     printf("AtomNo %d\n", AtomNo); // testing
     printf("Num %d\n", Num); // testing
     printf("CollisionT2[10][0] %.40f\n", CollisionT2[10][0]); // testing
     printf("CollisionT2[AtomNo][Num] %.40f\n", CollisionT2[AtomNo][Num]); // testing*/      LINE #

[this is only an extract, for illustration.]
(I apologise for posting the subroutine here.)

LINE @ will call the subroutine CollisionTimes() after which it will calculate the times and store them in the array. In the CollisionTimes(), when i output the values of CollisionT2[i][c], they are correct, but back at LINE #, CollisionT2[AtomNo][Num] gives a value of zero. What is the reason?

Many thanks!!

Jarvis

[Jarvis Loh]

Thanks Clifford..

I understand what you mean by binary floating point. But here is another similar situation.

It is a physics code and i am supposed to add 0.5at to a previous value of v(velocity). My a(acceleration) will change as time passes, and it will go towards a decreasing trend. At a certain stage, my v will not accumulate anymore. I tried increasing the number of decimal places, but all the rest of them are zero. Example(I output them as the code runs):

v 1.560000000000482100000000.. (STEP 1)
a -0.000000000000001952472081..

v 1.560000000000482100000000.. (STEP 2)
a -0.000000000000001410693948..

v does not increase. That is the problem.

Then I went on to do a little check:

        printf("velo   4.0/3.0  %.120lf\n", 4.0/3.0);
        scanf("%d", nol);

The output is "1.33333333333333330000000000000000000000000.."
Why is that so? It should be 3 all the way, but it seems that the 17th decimal place onwards cannot be used.

Anyone can enlighten me? Thanks alot.

Jarvis

[cpns]

You need to understand how floating point works. A real value as an infinite number of possible states, but a 64bit double has a finite number of discrete values. Large numbers are stored at a lower precision that small values. A double can store a value to approximately 16 significant figures. The formatted output functions do not show any more digits than that because the floating point representation cannot represent them.

This causes a problem when performing operations with very large and very small operands. In your case the value "a" is smaller than the precision of "v", so is essentially zero for the purposes of addition - especially for very small time intervals.

Perhaps a solution for very low accelerations is to use much larger time intervals. Either way I would recommend writing the expression as 0.5 * (a * t), this ensures you don't make 'a' even smaller and loose further precision before multiplying by 't'.

See http://en.wikipedia.org/wiki/Double_precision for details of the double precision representation.

Clifford

[cpns]

... you could also use the "long double" type, but at some point you still run into limitations. Also nothing guanrantees the actual size of long double other than it is at least as long as a double. On x86 it is typically 80bits since that is what the FPU supports. On PowerPC it is 128bits (which incidentally may mean that Apple Macintosh users using their machines for scientific processing may be better off not getting an Intel Mac). See: http://en.wikipedia.org/wiki/Long_double

In newer releases of GCC the the _Decimal128 data type is supported which would also solve your first problem. http://gcc.gnu.org/onlinedocs/gcc-4.3.0/gcc/Decimal-Float.html#Decimal-Float

Clifford

[Jarvis Loh]

I went off to search for libraries that enable me to increase my precision on the values calculated. Here is one site that I was looking at:

http://www.cs.berkeley.edu/~yozo/

I tried to use qd (quad-double) as stated in the above website.

I have done the following:
(1) copied the qd folder of header files into my list of header files
(2) typed #include <qd/dd_real.h> in my code

I know I have to add some commands to the linker command line, and also to add the library (is it a ".a" file?) into my C drive. But I can't find any .a files in the downloaded file from the site, and I do not know what are the commands for the command line. Can anyone tell me what i should do?

The creators of the quad-double algorithms indicated that to build the library, i have to edit the files "make.inc" and "make.config" and then type "make". Where should I be typing this, and again I seem to be unable to find the above stated files from my download.

Have I missed out on anything? Please enlighten me. Many thanks!

Jarvis

[Jarvis Loh]

Hi!

Is there anyone who can advise me on the above?

I appreciate the help. Thanks!

Jarvis

[Jarvis Loh]

I elaborate a bit more on the problem.

I have the following in my code:

include <stdio.h>

include <windows.h>

include <cmath>

include <stdlib.h>

include <cstdlib>

include <fstream>

include <iostream>

include <string>

include <ctime>

include <gl/gl.h>

include <gl/glut.h>

include <qd/dd_real.h>

define PI 3.1415926535898

using namespace std;

The <qd/dd_real.h> was just inserted. What should I add under "Add these commands to the linker command line"?

Instructions on how to build the library are included in the download:
1. You need to figure out which C++, C, (and optionally Fortran 95)
compiler you will be using. Run the configure script by typing

   ./configure CXX=<your C++ compiler> CC=<your C compiler> FC=<your Fortran 95 compiler>

 The script will attempt to automatically detect various 
 system-dependent features used during compilation (such as 
 C++ / C / Fortran compiler characteristics and availability of 
 certain system headers).

 If you want to specify a particular C++ / Fortran-90 compiler 
 flags, you can set them using CXXFLAGS and FCFLAGS.  For example:

   CXX=icpc CXXFLAGS='-O2 -mp' FC=ifort FCFLAGS='-O2 -mp' ./configure

 Some variable of interests are

   CXX       C++ compiler to use
   CXXFLAGS  C++ compiler flags to use
   CC        C compiler to use
   FC        Fortran 90 compiler
   FCFLAGS   Fortran 90 compiler flags to use
   FCLIBS    Fortran 90 libraries needed to to link with C++ code
   LDFLAGS   Linker flags

 For more build options, type "./configure --help".  In particular, 
 if you want to install to a custom path, do something like

   ./configure --prefix=/custom/path

1. Type "make". This will build the qd library.

Where should I type the "./configure...."? I don't really understand these instructions. I read that I have to run the make so that the makefile is activated..(Is this right?) But there is no make.exe in the downloaded files (as in http://www.cs.berkeley.edu/~yozo/). I was told to edit the names of make.inc and make.config to make, but there are no such files too.

Can anyone kindly give a step by step instruction on how to install the quad-double algorithm? Once again, my thanks.

Jarvis

[Wayne Keen]

"make" is the tool that takes all of the information about what files to compile, what files to link, and what things to put it the right search path and compiles everything into an executable or a library.

In the Unix world, there is somewhat of a standard "dance" that once can do to compile and install a library. It involves issuing the following commands at the command line

./configure
make
make install

The first step figures out what you have got, and where it is, in terms of getting your code to compile. Assuming you have all the right stuff, whis process will create a makefile for you.

The second step, as I described above, uses the information in the makefile to compile your code

The last step puts things where you need them.

To do this stuff on Windwos, I find it easiest to have a stand alone installation of MinGW and MSYS. MSYS gives you a nice terminal to work with, and provides the tools to do this sort of build with.

This all probably sounds confusing, but once you have done it a time or two, you will see that it is really easy, and I just did a bad job of explaining it.

Wayne

[Jarvis Loh]

I have installed MinGW and MSYS. But the MSYS website states the following:

Using MSYS with MinGW
The autotools that are installed by MSYS DTK do not work well and can't build DLLs. If you need newer versions of autoconf, automake and libtool then follow these instructions.

We will install autoconf, automake and libtool.

You will need to download updated versions.
autoconf
automake
libtool
Unpack autoconf, automake and libtool to a directory of your choice.
Install them with the 3 with the following command: ./configure --prefix=/mingw && make && make install

Do i unpacked them into ANY directory? I did the above, but got: "No such file or directory"

Then i tried help using ./configure --help, again "No such file or directory", so i presume there is something wrong with my installation?

so after ./configure, do i jus type make, then make install (as 2 separate commands)? My gd is in a folder on the desktop, do i need to copy them into c:\mingw first?

Thanks!

Jarvis

[Wayne Keen]

You can get newer versions of those tools through the MinGW website.

I would try your build with what you already have though, before I built new versions of
your tools.

You can also get newer versions of the tools you mention directly from the MinGW site.
They come as a tar file that has been zipped.

Wayne

[cpns]

(1)
Well that is rather the nature of "binary floating point", it cannot represent all decimal real values precisely. There again decimal real numbers cannot precisely represent all real-world values either (pi for example), so one is not necessarily more or less accurate than another, it is just that the result of a binary floating point calculation may not match that which you might calculate by hand on paper, you are likley to round both values, but decimal rounding and binary rounding result in different approximations of teh same value. You will find that the formatted output routines will perform 'sensible' decimal rounding for display and so you will only see this in the debugger.

It is unlikley the the real-world value you are representing is even that precisley measured. Depending on you application you may find that 9.9999...E-17 is as near as makes no difference to 1.0e16 ;

Note also that you can use exponential notation as literal constants, so for radability you could use 1.0e16 in your initialiser.

(2) The format specifier %f is for single precision (float) types, but you passed it a double. Cast the double to float or use %lf

Note that you can achieve more precision if necessart with the long double data type. Althouigt the formatted output routine may not support it, so you may want to cast such values to float ro double for display.

Clifford