Conference turris::fortran

    
            Note that there is no "absolute accuracy" in the whole range
            of floating point. Real*8 in G- or T-floating has a 53-bit
            fraction, that's it.
    
            Look at your example:
    
            x=10.0
    
            This is exactely represented, the fraction starts with 101,
            the rest is zero. (Forget about the hidden most significant
            bit for the moment.)
    
            y=y*x
    
            This will apparently calulate 100,1000,10000 etc.
            But each time you multiply with 10 the fraction will get around
            2.3 more non-zero bits. 53/2.3 = 23, so somewhere around
            10**23 you get underflow, the stored datum is not exactely 10**k.
            After that you get more underflow every round in the loop.
    
            I can't say precisely how z**i does it, some accumulated
            underflow is probably avoided, but still there is the
            limitation  of the 53-bit fraction.
    
            Finally, your program does not write out the eaxct values
            of y and z. Write(6,*) writes the values rounded to 16 digits.
            When investigating behaviours like this it's often more
            practical to type out the hexadecimal representation.
    
    
    /�ke
    

    "an accuracy problem" - um, not really.  Just an accuracy limitation.  
    Floating point numbers are stored in computers with a limited number of
    binary digits.
    
    "absolute accuracy" - then don't count anything that can't be counted
    as a reasonable size integer.  (In meterology??)
    
    "overhead of function calls" - inline code is sometimes faster than a
    function call, but not always.  If the function is decently written, it
    can more than pay back the overhead of its call.  A well-written
    function call will also often be MORE accurate than spelling out the
    steps in detail in your own program, because the folks writing the
    functions have carefully picked algorithms that tend to NOT accumulate
    error as they go.
    
    There's a whole field here - Numerical Analysis, Scientific Computing, 
    Algorithms, How To Program Good In Fortran - your customer may want to
    find a book with a title something like that, or a good course at a
    community college.

By the way, the standard way to implement x**i for i not a constant, is to use
steps of squaring and multiplying by x, rather than simply multiplying by x.
I'm sure this is the method the RTL uses.  For example, to compute x**35, it
would follow a procedure such as
	y = x
	y = y*y = x**2
	y = y*y = x**4
	y = y*y = x**8
	y = y*y = x**16
	y = y*x = x**17
	y = y*y = x**34
	y = y*x = x**35
Only 7 multiplies, compared to 34 or 35 for the naive method in your program.
Fewer multiplies means it's faster, and also that there are fewer opportunities
for round-off error to accumulate.

For some values of i, there are faster ways -- for example, computing x**15
this way takes 6 multiplies, but you can do it with only 5:
	z = x*x = x**2
	z = z*z = x**4
	z = z*x = x**5
	y = z*z = x**10
	y = y*z = x**15
but I don't think the RTL does this.

    Re : .1-.3
    
    	Forgot one point.  The customer said, he had run the same program
    on a Sun UltraSparc system and he got exact results (meaning 10.0**23
    or 10.0**24, etc.)
    
    	So, where does Alpha differ ?
    
    	Regards.
    
    	R. Devarajan
    
    

I notice that Digital chose to use 16 digits of precision in our list
directed output.  Is it possible that Sun chose to use only 15 digits
of precision?  You would not see this behavior on our system if we had
chosen to print the results with only 15 digits.

The Sun Fortran docs I have do not specify what format is used for list-directed
output, only to say that it chooses a "simple" representation over an accurate
one.

As .-1 noted, if you print this stuff out E20.15, they all are ".1000...000Enn"
so they look accurate... they're really just "simple".

/Stan

    Thanks for all the inputs.  I have sent a mail to the customer that
    they use a format statement with E20.15 for a REAL *8 variable and
    E21.16 for a REAL *16 variable.  I have checked with the program that
    if the precision goes beyond 15 digit for a R*8 variable, "accuracy" is
    lost (basically some calculations give 9s).  In these cases, R*16 is
    helpful.
    
    Thanks and regards,
    
    R. Devarajan