Fortran or Qbasic calculator accuracy

Please excuse me if I am asking a question which has already been answered. I have found that the accuracy of calculation in Fortran and Qbasic is not good enough for my needs. I am using Windows XP, and the inbuilt calculator in that is good, but I cannot get the same accuracy in my programme.
The problem is adding 5.92719e24 to 1,000 (or to 1,000,000) . The calculator does this fine, but both Fortran and Qbasic get it wrong. The correct answer is 5,927,190,000,000,000,000,001,000
(or 5,927,190,000,000,000,001,000,000) .
Fortran gives 5.97219d+24 for both answers, while
Qbasic gives 5.972189887017193d+24 for both answers. Note that Qbasic has given a smaller result than both original numbers. What is going on, why are the calculations not being done correctly?
I am fully (I think) conversant with Qbasic, but not very familiar with Fortran.
If it is not against forum rules, I would appreciate an answer by email as well as on this forum. The reason is I am not familiar with forum usage, and may not be able to find this forum page again.
email carmam "at" tiscali "dot" co "dot" uk
Thank you.

 

Not a reply, an addendum.
In Qbasic, I changed exponent references from eg 1.234e24 to 1.234d24 . This gave me a bit more accuracy, the result is now correct if I add 1d10 or greater to 5.97219d24 . Can I achieve greater accuracy than this?
I can also find this forum again, so no problem with replies.

Thanks.

 

Hi! Welcome to MajorGeeks!

As you may know, the data type for the numbers you are working with in Fortran and Qbasic is called double precision floating point, (REAL*4 in Fortran,) which is accurate to at least 15 decimal digits. Higher precision types are quadrupled, accurate to at least 33 digits, and extended types, which are even more accurate. Then, there's arbitrary precision, sometimes called infinite precision, which expands in size as necessary to preserve accuracy.

According to http://stackoverflow.com/questions/...nough-precision-for-implementing-a-calculator and its links, the Windows calculator uses arbitrary precision for simple operations, (displaying 32 digits accurately, but keeping more in memory,) and less precision for advanced operations. A quad or extended data type can display as or more accurately, but are not as accurate as arbitrary precision in actual computations.

Fortran calls quads REAL*8 and has an extended type called REAL*10, but doesn't support arbitrary precision. You may be able to make arbitrary precision out of strings or arrays of integers, but it could get complicated, so supported types are generally better unless more precision is necessary.

QBasic only has single and double, not quads, extended types, or arbitrary precision. I think I made custom data types in QBasic for these, but I don't remember where they are, how finished they are, or if they are compatible with Fortran. If you're interested, and I can find them, I can probably post them here. Fortran's supported types are easier to use, as QBasic won't use the math operators on custom types, so I had to make custom functions to act as operators.

 

In Fortran, I am using real*8, but may have the declarations wrong. It is many years since I used it. Here is the applicable line :-
REAL*8 Period, Me, Ms, Grav, Pie2, Rad, days
This should set all these variables to quad precision. Me is a fixed value of 5.97219d24, and Ms is input as (in this instance) 1000, so these two lines -
print*,"Ms ",Ms
print*,"Me + Ms " , Me + Ms
should return 1000, and 5.972190000000000000001000 .
Ms is correct at 1000 , but Me + Ms returns 5.97219E+24
If I set Me to 5.97219e24 (or e+24), the result is plainly wrong. it is
5.97218989E+24 , which is smaller than 5.97219d24

I think I found the problem. I set Me to 5.97219q24 and got an error when compiling "quadruple precision floating point unsupported". I am using G77 compiler.

 

You're right. g77 doesn't support quad-precision because it was used before GCC 4.0. The quad-precision library was included in GCC 4.6 and higher versions, which use f77 (FORTRAN 95) instead of g77 (FORTRAN 77). If you only need integers, no fractional part, you might be able to try INTEGER*8, instead: http://www.star.bris.ac.uk/~mbt/docs/ssn73.htx/node6.html I don't know for sure if this nonstandard method will work in g77, but if so, it probably requires a 64-bit system.

 

According to http://www-classes.usc.edu/engr/ce/108/text/fbk01.htm I was wrong about FORTRAN's data types. REAL*8 is a double, so REAL*10 is not really a quad, but probably the 80-bit type used by the math coprocessor. REAL*16 would be a quad. g77 may support REAL*10, but it will likely round to a double in some cases.

 

I found http://crd-legacy.lbl.gov/~dhbailey/mpdist/ has some FORTRAN 77 and 90 packages for high precision arithmetic, including arbitrary precision.

 

I have swapped onto FortranG95, and now I have enough decimal places for precision, but unfortunately, the calculations are still not precise. Here is what is happening. Me is set as 5.97219e24 (or d24 or q24), and Ms is input as 1000 in this instance.

Some pertinant programme lines (in order)
real*16 Me
real*16 Ms

print*,"Me ", Me >>> prints 5.97219000000000039059456E+24
print*,"Ms ",Ms >>> prints 1000
print*,"Me + Ms " , Me + Ms >>> prints 5.97219000000000039059556E+24

5.97219e24 is an exact precise number. Why does Fortran insist on making it larger? The addition is correct, but Me is not.

 

hmm, let's see...
5.97219d24 is stored as 5.97219000000000039059456E+24 in a quad, due to standard floating point rounding. (5.97219e24 may be rounded to 5.972189887017193E+24 in some cases.) 5.97219q24 should be stored accurately enough that the rounded part is below the true decimal point. So, it looks like the quad is converted to or from a double at some point. This may happen when printing the variable. You could try to test the value to see if it equals 5.97219q24 before and after each print, like:
if (Me == 5.97219q24) print *, 'pass'
I'm not positive that that code won't convert both numbers, anyway, but it might show something.

http://ondrejcertik.blogspot.com/2012/01/when-double-precision-is-not-enough.html talks of compiling Fortran to better handle quad precision, in case it helps.

 

<<You could try to test the value to see if it equals 5.97219q24 before and after each print, like:
if (Me == 5.97219q24) print *, 'pass'
I'm not positive that that code won't convert both numbers, anyway, but it might show something.>>

I have got round the problem by subtracting 390594560 from
5.97219000000000039059456e24 . This gives the correct figure for Me, which you may have realised is the mass of the earth. This is fine untill I want to use a different planet, in which case I would have to correct for every one I use. I will try your suggestions.

Thanks for all your help.

 

Good thinking. If you correct Me immediately after setting the value, and it holds the corrected value after printing and everything, then it's probably something that happens when setting the variable that converts it to a double and back.