Re: [Gcl-devel] Re: [Maxima] Re: Epsilon calculation

[Camm Maguire]

Greetings!  Just a followup note to state that I think this issue has
now been resolved in a portable way.  Basically, as others have
stated, the FMEM, etc. bits were not succeeding in forcing the
intermediate results out of the register and into memory.  'volatile'
does do this, and the following code is giving the correct answers now
on all 11 Debian architectures:

        {
          
          volatile double rd,dd,td,td1;
          volatile float  rf,df,tf,tf1;
          int i,j;
#define MAX 500
          
          for (rf=1.0f,df=0.5f,i=j=0;i<MAX && j<MAX && 
df!=1.0f;i++,df=1.0f-(0.5f*(1.0f-df)))
            for (tf=rf,tf1=tf+1.0f,j=0;j<MAX && 
tf1!=1.0f;j++,rf=tf,tf*=df,tf1=tf+1.0f);
          if (i==MAX||j==MAX)
            printf("WARNING, cannot calculate float_epsilon: %d %d %f   %f %f 
%f\n",i,j,rf,df,tf,tf1);
          float_epsilon=rf;

          for (rf=1.0f,df=0.5f,i=j=0;i<MAX && j<MAX && 
df!=1.0f;i++,df=1.0f-(0.5f*(1.0f-df)))
            for (tf=rf,tf1=1.0f-tf,j=0;j<MAX && 
tf1!=1.0f;j++,rf=tf,tf*=df,tf1=1.0f-tf);
          if (i==MAX||j==MAX)
            printf("WARNING, cannot calculate float_negative_epsilon: %d %d %f  
 %f %f %f\n",i,j,rf,df,tf,tf1);
          float_negative_epsilon=rf;
          
          for (rd=1.0,dd=0.5,i=j=0;i<MAX && j<MAX && 
dd!=1.0;i++,dd=1.0-(0.5*(1.0-dd)))
            for (td=rd,td1=td+1.0,j=0;j<MAX && 
td1!=1.0;j++,rd=td,td*=dd,td1=td+1.0);
          if (i==MAX||j==MAX)
            printf("WARNING, cannot calculate double_epsilon: %d %d %f   %f %f 
%f\n",i,j,rd,dd,td,td1);
          double_epsilon=rd;

          for (rd=1.0,dd=0.5,i=j=0;i<MAX && j<MAX && 
dd!=1.0;i++,dd=1.0-(0.5*(1.0-dd)))
            for (td=rd,td1=1.0-td,j=0;j<MAX && 
td1!=1.0;j++,rd=td,td*=dd,td1=1.0-td);
          if (i==MAX||j==MAX)
            printf("WARNING, cannot calculate double_negative_epsilon: %d %d %f 
  %f %f %f\n",i,j,rd,dd,td,td1);
          double_negative_epsilon=rd;
          
        }


Take care,



Raymond Toy <address@hidden> writes:

> >>>>> "Stavros" == Stavros Macrakis <address@hidden> writes:
> 
>     Stavros> Raymond--
>     Stavros> I understand your theory, but I don't think it's correct.
> 
>     Stavros> Here's my understanding of IEEE floating point.  Suppose we're 
> working
>     Stavros> with a 5-bit stored mantissa and an 8-bit register mantissa.  
> When we
>     Stavros> add .xxxxx and .00010000, the register result is .xxxxx100, 
> which should
>     Stavros> be rounded according to the 'even' rule.  Now what if we add 
> .xxxxx and
>     Stavros> .00010001?  Well, if there were no sticky bit, we would still get
>     Stavros> .xxxxx100 in the register (since the final 1 has been shifted 
> off),
>     Stavros> which again would round to even when stored.  But if we have a 
> sticky
>     Stavros> bit (which remembers if anything non-zero has been shifted off), 
> then
>     Stavros> the register sum is .xxxxx101.  When we store .xxxxx101, we need 
> to
>     Stavros> round to five bits, and since 101 is greater than 100, we round 
> up.  So
>     Stavros> we shouldn't be losing the round-up quality when storing.  After 
> all, it
>     Stavros> would be a pity if despite the elaborate IEEE machinery, you 
> couldn't
>     Stavros> get a simple sum of two floats to round correctly except in 
> intermediate
>     Stavros> register results....
> 
>     Stavros> Or have I gotten something wrong here?
> 
> I think I was wrong.  Your reasoning looks good and I think this is
> what happens if the intermediate results are stored.  If the compiler
> decided to keep everything in registers, things might be different.
> You didn't say what FMEM did and whether it was in the same file or
> not.
> 
>     Stavros> Anyway, your theory does not explain the Lisp facts.  How is it 
> that
>     Stavros> Emacs Lisp and GCL manage to round their results correctly -- 
> remember,
>     Stavros> they're written in C, too -- but native C does not?
> 
> Because they always store the results away in memory after each basic
> operation?  I'm guessing.  I don't really know.  I suppose you could
> look at the C file that gcl produces to see what it's doing.
> 
> All I can say is that since the results you get are different, they
> must be doing something different.  I just don't know exactly what
> that might be.
> 
> Ray
> 
> 
> _______________________________________________
> Gcl-devel mailing list
> address@hidden
> http://mail.gnu.org/mailman/listinfo/gcl-devel
> 
> 

-- 
Camm Maguire                                            address@hidden
==========================================================================
"The earth is but one country, and mankind its citizens."  --  Baha'u'llah