8
$\begingroup$

by default, if a number has decimal . after it, then Mathematica will do computation using machine Precision, which on my PC (intel hardware) running windows 7 64 bit is double Precision.

I'd like to get the computation also but using single and quad precision, to match the following small Fortran program output.

Is there a recommend way to do this? Here is the fortran program and its output

PROGRAM foo
IMPLICIT NONE
INTEGER:: i

INTEGER, PARAMETER ::                              &
    sp = kind(1.0),                                &
    dp = selected_real_kind(2*precision(1.0_sp)),  &
    qp = selected_real_kind(2*precision(1.0_dp))

REAL(kind=sp) ::  sum1,x1
REAL(kind=dp) ::  sum2,x2
REAL(kind=qp) ::  sum3,x3

sum1=0.0_sp; sum2=0.0_dp; sum3=0.0_qp
x1 = 0.00001_sp
x2 = 0.00001_dp
x3 = 0.00001_qp

DO i=1,10**5
   sum1 = sum1 + x1
   sum2 = sum2 + x2
   sum3 = sum3 + x3
END DO

PRINT *,sum1; PRINT *,sum2; PRINT *,sum3;
END

Compiled using gfortran foo.f90 gives

ex1>./a.out 
   1.00099015      ---> SINGLE
  0.99999999999808376      --> DOUBLE
  0.999999999999999999999999999998395508  ----> QUAD

The first line is single, the second is double (which Mathematica matches exactly) and the third is quad.

Here is the Mathematica basic code which gives output that matches the double Precision:

x1=0.00001;
sum1=0.0;
Do[sum1=sum1+x1,{i,1,10^5}]
InputForm[sum1]

Which gives

 0.9999999999980838

I tried doing 

 x1=0.00001`32;
 sum1=0.0`32;

To see if will something close to the quad result, but it had not effect on final output. I still get 0.9999999999980838

How to change the above Mathematica code to make it give the single and quad Precision as shown by Fortran?

Improve this question
$\endgroup$
3
  • $\begingroup$ One thing going on here is 0.0`32 is for some reason MachinePrecision. You can actually initialize sum1=0; and get your desired precision in the calculation. (It does not match the fortran results though) $\endgroup$
    – george2079
    Commented Jan 24, 2018 at 18:38
  • $\begingroup$ @george2079 I do not think this gives same as Fortran. For example, $x1=0.00001\`7$ then $sum=0$ results in $0.9999999999999972025\`6$ while Fortran gives $1.00099015$ I think setting sum=0 instead of sum=0.0 makes it do arbitray precision somewhere internally? $\endgroup$
    – Nasser
    Commented Jan 24, 2018 at 18:45
  • $\begingroup$ right, for anything other than machine precision mathematica is not doing IEEE arithmetic and so you will never get the same results. $\endgroup$
    – george2079
    Commented Jan 24, 2018 at 18:50

2 Answers 2

Reset to default
4
$\begingroup$

The following emulates the various precisions, by directly rounding the operation as done in floating-point. The quad seems off by 1 bit, if Fortran uses IEEE quad precision. Also, Fortran output seems to have an extra digit or two, compared to Mathematica's at the specified precision.

singlebits = 23;
x1 = 0.00001 // Rationalize;
sum1 = 0;
Do[sum1 = Round[sum1 + x1, 2^Floor@Log2[N[(sum1 + x1), 9]*2^-singlebits]], {i, 1, 10^5}]
SetPrecision[sum1, 9]
(*  1.00099015  *)

(* machine precision *)
x1 = 0.00001;
sum1 = 0.0;
Do[sum1 = sum1 + x1, {i, 1, 10^5}]
SetPrecision[sum1, 17]
(*  0.99999999999808376  *)

quadbits = 112;
x1 = 0.00001 // Rationalize;
sum1 = 0;
Do[sum1 = Round[sum1 + x1, 2^Floor@Log2[N[sum1 + x1, 36]*2^-quadbits]], {i, 1, 10^5}]
SetPrecision[sum1, 36]
(*  0.999999999999999999999999999998395508  *)
Improve this answer
$\endgroup$
2
  • $\begingroup$ (* should be 113? *) — No, why? You use 23 for single precision, then 112 is the corresponding number of bits (of number's tail) in 128-bit IEEE format. What seems to be a bit incorrect is the use of N. Namely, N will do some rounding to represent the number in the precision requested, and then you do additional Round, which will result in double rounding and thus potentially increase the error. $\endgroup$
    – Ruslan
    Commented Apr 13, 2019 at 22:00
  • $\begingroup$ @Ruslan Thanks. I don't remember why I wrote that. Perhaps I was reading a description of the format that was unclear whether 113 bits was the size of the significand or the number of bits stored. $\endgroup$
    – Michael E2
    Commented Apr 14, 2019 at 0:10
2
$\begingroup$ 
$Version
x1 = 0.00001 // Rationalize // N[#, 32] &
sum1 = 0 // N[#, 32] &
Do[sum1 = sum1 + x1, {i, 1, 10^5}]
sum1 // InputForm

"11.2.0 for Mac OS X x86 (64-bit) (September 11, 2017)"

0.000010000000000000000000000000000000

0

0.999999999999999999999999999999999\
9999999999999999999979455`32.

Note that the zero remained exact.

Improve this answer
$\endgroup$
3
  • 1
    $\begingroup$ Thanks. But are you sure this is done using quad precision? Mathematica result is 0.9999999999999999999999999999999999999999999999999999979455 but Fortran is 0.999999999999999999999999999998395508 I know Mathematica puts $\'32$ at the end of the number, but this is more than 32 digits precision. $\endgroup$
    – Nasser
    Commented Jan 24, 2018 at 19:59
  • $\begingroup$ To ensure that the result has 32 digit precision, the internal calculations use longer numbers and track the precision. The longer result is shown but -- as indicated -- only 32 digits are considered precise since that is the highest precision of the lowest precision input, and the precision tracking did not see any lose of precision. $\endgroup$
    – Bob Hanlon
    Commented Jan 24, 2018 at 20:09
  • $\begingroup$ @Nasser - I forgot to add your name to my comment above. $\endgroup$
    – Bob Hanlon
    Commented Jan 24, 2018 at 20:36

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.