Increasing accuracy of “sum()” command in Fortran

user_10412550 0
Viewed 155

I am trying to perform summation of an array.

There are 1024 elements in that array and on applying the command "sum(a(:))", where a(:) is the array, I get the answer 1981.9072.

If I do the same summation of 1024 elements in Excel sheet the answer is 1981.93530 which is the right answer. So, a difference of 0.0281 is observed between the above two values. As I increase the number of elements in a(:) the difference in value obtained from "sum()" and Excel sheet increases.

I think the "sum()" value is different due to rounding off error. How do I get the true value (Excel value) using "sum()" without any rounding off error ?

fortran
1 Answers

REAL(4) is a 32-bit IEEE-754 floating-point number which is imprecise.

You absolutely need to read and understand this important QA on the matter: Why are floating point numbers inaccurate?

Here's a page on the same topic from the GCC Fortran compiler docs: https://gcc.gnu.org/onlinedocs/gcc-3.3.6/g77/Floating_002dpoint-Errors.html

You need to change the type of a to REAL(8) - but note that even if you do, the sum still won't be completely accurate (it will just be less inaccurate).

This can be reproduced consistently in any language that implements IEEE-754 32-bit floating-point numbers, like C# (I'm using C# as I don't have access to a Fortran90 compiler right now):

static void Main()
{
    Single t = 0f;

    for( Int32 i = 0; i < 1024; i++ )
    {
        t += 1.93548369f; // This is a 32-bit float literal
    }

    Console.WriteLine( t ); // "1981.9072"
}

If you sum the numbers using a double-precision (64-bit) floating number type you get the same result as Excel:

static void Main()
{
    Double d = 0d;

    for( Int32 i = 0; i < 1024; i++ )
    {
        d += 1.93548369d; // This is a 64-bit float literal
    }

    Console.WriteLine( d ); // "1981.935302734375" --> "1981.9350"
}

If you want an accurate and precise answer then you need to use a data-type that doesn't use IEEE-754 - I'm not a Fortran user so I don't know what/if Fortran's equivalent of C#'s Decimal is, but when I do the calculation with Decimal I get 1981.93529856 which is different to Excel's answer (this is because Excel uses 64-bit IEEE-754 instead of a real Decimal type):

static void Main()
{
    Decimal dec = 0M;
    
    for( Int32 i = 0; i < 1024; i++ )
    {
        dec += 1.93548369M; // This is a decimal literal, note the `M`.
    }

    Console.WriteLine( dec ); // "1981.93529856" --> "1981.9352"
}
user_1591450
Related