how to get single, and quad Precision out of Mathematica which match Fortran result

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


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.0000132;
sum1=0.032;


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?

• One thing going on here is 0.032 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) Jan 24, 2018 at 18:38
• @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? Jan 24, 2018 at 18:45
• right, for anything other than machine precision mathematica is not doing IEEE arithmetic and so you will never get the same results. Jan 24, 2018 at 18:50

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  *)


• (* 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. Apr 13, 2019 at 22:00
$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\ 999999999999999999997945532.  Note that the zero remained exact. • 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. Jan 24, 2018 at 19:59