how to get single, and quad Precision out of Mathematica which match Fortran result
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 *)
$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.