[Gcl-devel] Re: Epsilon calculation

[Camm Maguire]

"Stavros Macrakis" <address@hidden> writes:

> Camm,
> 
> The new Epsilon calculations look excellent.
> 
> But I am puzzled by the value 1.1107651257113995E-16 (call that Q765)
> for double-float-epsilon.  When running in Lisp, the correct value is in
> fact 1.1102230246251568d-16 (call that Q223).  That is, 1+Q223-1 is
> non-zero.
> 
> I have implemented your algorithm both in C and in Lisp.  In Lisp, it
> gives the correct result (Q223).  But in C, it gives Q765, UNLESS you
> don't force intermediate results to memory.  But surely in Lisp, all the
> intermediate results are in fact going to memory!  And the binary
> representation of Q765 looks bizarre.
> 
> It is not your code that's the problem -- I tried a different algorithm
> and got the same result.
> 
> Any clue what is going on here?  Or am I just hallucinating?
> 
>       -s
> 
> 
> 
Greetings!  Beyond this problem, another user reported that the
algorithm looped forever on his Athlon XP.  So I've redone the
calculation again just in the IEEE case.  There is still a
computational anomaly in the double precision case, which I've
temporarily overriden to give the right answers for now.  Comments
welcome!

=============================================================================
Index: num_co.c
===================================================================
RCS file: /cvsroot/gcl/gcl/o/num_co.c,v
retrieving revision 1.6
retrieving revision 1.7
diff -u -r1.6 -r1.7
--- num_co.c    17 Oct 2002 01:40:56 -0000      1.6
+++ num_co.c    18 Oct 2002 15:00:56 -0000      1.7
@@ -1203,6 +1203,63 @@
 void  _assure_in_memory (void *p);
 
 #define FMEM(f) _assure_in_memory(&f)
+
+#ifdef IEEEFLOAT
+/* from ieee754.h */
+
+typedef union {
+  float f;
+  
+  /* This is the IEEE 754 single-precision format.  */
+  struct float_bits
+  {
+#ifndef LITTLE_END
+    unsigned int negative:1;
+    unsigned int exponent:8;
+    unsigned int mantissa:23;
+#else                          /* Big endian.  */
+    unsigned int mantissa:23;
+    unsigned int exponent:8;
+    unsigned int negative:1;
+#endif                         /* Little endian.  */
+  } ieee;
+} IEEE_float;
+  
+typedef union  {
+
+  double d;
+  
+  /* This is the IEEE 754 double-precision format.  */
+  struct double_bits
+  {
+#ifndef LITTLE_END
+    unsigned int negative:1;
+    unsigned int exponent:11;
+    /* Together these comprise the mantissa.  */
+    unsigned int mantissa0:20;
+    unsigned int mantissa1:32;
+#else                          /* Big endian.  */
+    /* # if    __FLOAT_WORD_ORDER == BIG_ENDIAN */
+    /*     unsigned int mantissa0:20; */
+    /*     unsigned int exponent:11; */
+    /*     unsigned int negative:1; */
+    /*     unsigned int mantissa1:32; */
+    /* # else */
+    /* Together these comprise the mantissa.  */
+    unsigned int mantissa1:32;
+    unsigned int mantissa0:20;
+    unsigned int exponent:11;
+    unsigned int negative:1;
+    /* # endif */
+#endif                         /* Little endian.  */
+  } ieee;
+} IEEE_double;
+  
+
+#endif
+
+
+
 void
 init_num_co(void)
 {
@@ -1289,33 +1346,6 @@
 #endif
 #endif
 
-#ifdef MV
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-#endif
-
 #if defined(S3000) && ~defined(DBL_MAX_10_EXP)
        l[0] = 0x7fffffff;
        l[1] = 0xffffffff;
@@ -1351,37 +1381,50 @@
        biggest_float = FLT_MAX;
 #endif
        
-
-/* We want the smallest number not satisfying something,
-   and so we go quickly down, and then back up.  We have
-   to use a function call for test, since in line code may keep
-   too much precision, while the usual lisp eql,is not
-   in line.
-   We use SMALL as a multiple to come back up by.
-   We use FMEM(double_negative_epsilon)
-   to force the quantity into memory by taking its address
-   and then passing it to a function.
-*/
-
+       
 #ifdef IEEEFLOAT
  { 
-   double div,td;
-   float ts;
-   for (float_epsilon=1.0,div=0.5;(float)div!=1.0;div=1.0-(0.5*(1.0-div)))
-     for 
(ts=float_epsilon;FMEM(ts),!SF_EQL((float)(1.0+ts),(float)1.0);float_epsilon=ts,ts*=div);
-
-   for 
(float_negative_epsilon=1.0,div=0.5;(float)div!=1.0;div=1.0-(0.5*(1.0-div)))
-     for 
(ts=float_negative_epsilon;FMEM(ts),!SF_EQL((float)(1.0-ts),(float)1.0);float_negative_epsilon=ts,ts*=div);
-
-   for (double_epsilon=1.0,div=0.5;(double)div!=1.0;div=1.0-(0.5*(1.0-div)))
-     for 
(td=double_epsilon;FMEM(td),!LF_EQL((double)(1.0+td),(double)1.0);double_epsilon=td,td*=div);
+   IEEE_float ief;
+   IEEE_double ied;
 
-   for 
(double_negative_epsilon=1.0,div=0.5;(double)div!=1.0;div=1.0-(0.5*(1.0-div)))
-     for 
(td=double_negative_epsilon;FMEM(td),!LF_EQL((double)(1.0-td),(double)1.0);double_negative_epsilon=td,td*=div);
+   if (sizeof(ief)!=sizeof(ief.f))
+     FEerror("Bad ieee float definition\n");
+   if (sizeof(ied)!=sizeof(ied.d))
+     FEerror("Bad ieee float definition\n");
+
+   for (float_epsilon=ief.f=1.0,ief.ieee.mantissa=1;
+       FMEM(ief.f),!SF_EQL((float)(1.0+ief.f),(float)1.0);
+       float_epsilon=ief.f,ief.ieee.exponent--);
+
+   for (float_negative_epsilon=ief.f=1.0,ief.ieee.mantissa=1;
+       FMEM(ief.f),!SF_EQL((float)(1.0-ief.f),(float)1.0);
+       float_negative_epsilon=ief.f,ief.ieee.exponent--);
+
+   for (double_epsilon=ied.d=1.0,ied.ieee.mantissa1=1;
+       FMEM(ied.d),!LF_EQL((1.0+ied.d),1.0);
+       double_epsilon=ied.d,ied.ieee.exponent--);
+   /* FIXME double calculations end one iteration too early */
+   double_epsilon=ied.d;
+
+   for (double_negative_epsilon=ied.d=1.0,ied.ieee.mantissa1=1;
+       FMEM(ied.d),!LF_EQL((1.0-ied.d),1.0);
+       double_negative_epsilon=ied.d,ied.ieee.exponent--);
+   double_negative_epsilon=ied.d;
 
  }
 #else
 
+ /* We want the smallest number not satisfying something,
+    and so we go quickly down, and then back up.  We have
+    to use a function call for test, since in line code may keep
+    too much precision, while the usual lisp eql,is not
+    in line.
+    We use SMALL as a multiple to come back up by.
+    We use FMEM(double_negative_epsilon)
+    to force the quantity into memory by taking its address
+    and then passing it to a function.
+ */
+ 
 #define SMALL 1.05     
        for (float_epsilon = 1.0;
             FMEM(float_epsilon),!SF_EQL((float)(1.0 + 
float_epsilon),(float)1.0);
=============================================================================


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