/[gcl]/gcl/o/num_co.c
ViewVC logotype

Diff of /gcl/o/num_co.c

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1.7 by camm, Fri Oct 18 15:00:56 2002 UTC revision 1.8 by camm, Fri Nov 1 02:30:43 2002 UTC
# Line 1190  Limagpart(void) Line 1190  Limagpart(void)
1190          }          }
1191  }  }
1192    
 static float sf1,sf2;  
 static int sf_eql(void)  
 {return(sf1==sf2);}  
   
   
   
 int lf_eqlp(double *p, double *q);  
 #define LF_EQL(a,b) (lf1=a,lf2=b,lf_eqlp(&lf1,&lf2))  
 #define SF_EQL(a,b) (sf1=a,sf2=b,sf_eql())  
   
 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  
   
   
   
1193  void  void
1194  init_num_co(void)  init_num_co(void)
1195  {  {
# Line 1267  init_num_co(void) Line 1197  init_num_co(void)
1197          double smallest_double, smallest_norm_double, biggest_double;          double smallest_double, smallest_norm_double, biggest_double;
1198          float float_epsilon, float_negative_epsilon;          float float_epsilon, float_negative_epsilon;
1199          double double_epsilon, double_negative_epsilon;          double double_epsilon, double_negative_epsilon;
         double lf1,lf2;  
   
1200    
1201    
1202  #ifdef VAX  #ifdef VAX
# Line 1381  init_num_co(void) Line 1309  init_num_co(void)
1309          biggest_float = FLT_MAX;          biggest_float = FLT_MAX;
1310  #endif  #endif
1311                    
1312            {
1313              
1314              volatile double rd,dd,td,td1;
1315              volatile float  rf,df,tf,tf1;
1316              int i,j;
1317    #define MAX 500
1318              
1319              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)))
1320                for (tf=rf,tf1=tf+1.0f,j=0;j<MAX && tf1!=1.0f;j++,rf=tf,tf*=df,tf1=tf+1.0f);
1321              if (i==MAX||j==MAX)
1322                printf("WARNING, cannot calculate float_epsilon: %d %d %f   %f %f %f\n",i,j,rf,df,tf,tf1);
1323              float_epsilon=rf;
1324    
1325              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)))
1326                for (tf=rf,tf1=1.0f-tf,j=0;j<MAX && tf1!=1.0f;j++,rf=tf,tf*=df,tf1=1.0f-tf);
1327              if (i==MAX||j==MAX)
1328                printf("WARNING, cannot calculate float_negative_epsilon: %d %d %f   %f %f %f\n",i,j,rf,df,tf,tf1);
1329              float_negative_epsilon=rf;
1330              
1331              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)))
1332                for (td=rd,td1=td+1.0,j=0;j<MAX && td1!=1.0;j++,rd=td,td*=dd,td1=td+1.0);
1333              if (i==MAX||j==MAX)
1334                printf("WARNING, cannot calculate double_epsilon: %d %d %f   %f %f %f\n",i,j,rd,dd,td,td1);
1335              double_epsilon=rd;
1336    
1337              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)))
1338                for (td=rd,td1=1.0-td,j=0;j<MAX && td1!=1.0;j++,rd=td,td*=dd,td1=1.0-td);
1339              if (i==MAX||j==MAX)
1340                printf("WARNING, cannot calculate double_negative_epsilon: %d %d %f   %f %f %f\n",i,j,rd,dd,td,td1);
1341              double_negative_epsilon=rd;
1342              
1343            }
1344    
1345                    
1346  #ifdef IEEEFLOAT  #ifdef IEEEFLOAT
1347   {          /* Maybe check for "right" answer here */
    IEEE_float ief;  
    IEEE_double ied;  
   
    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);  
              float_epsilon /= 2.0)  
                 ;  
         while(SF_EQL((float)(1.0 + float_epsilon),(float)1.0))  
           { FMEM(float_epsilon); float_epsilon=float_epsilon*SMALL;}  
         for (float_negative_epsilon = 1.0;  
             FMEM(float_negative_epsilon), !SF_EQL((float)(1.0 - float_negative_epsilon) ,(float)1.0);  
              float_negative_epsilon /= 2.0)  
                 ;  
         while(SF_EQL((float)(1.0 - float_negative_epsilon) ,(float)1.0))  
           float_negative_epsilon=float_negative_epsilon*SMALL;  
         for (double_epsilon = 1.0;  
              FMEM(double_epsilon), !(LF_EQL(1.0 + double_epsilon, 1.0));  
              double_epsilon /= 2.0)  
                 ;  
         while((LF_EQL(1.0 + double_epsilon, 1.0)))  
           double_epsilon=double_epsilon*SMALL;  
           ;  
         for (double_negative_epsilon = 1.0;  
            FMEM(double_negative_epsilon),  
                !LF_EQL(1.0 - double_negative_epsilon , 1.0);  
              double_negative_epsilon /= 2.0)  
                 ;  
         while(LF_EQL(1.0 - double_negative_epsilon , 1.0))  
           double_negative_epsilon=double_negative_epsilon*SMALL;  
           ;  
1348  #endif  #endif
1349    
1350          make_constant("MOST-POSITIVE-SHORT-FLOAT",          make_constant("MOST-POSITIVE-SHORT-FLOAT",

Legend:
Removed from v.1.7  
changed lines
  Added in v.1.8

savannah-hackers-public@gnu.org
ViewVC Help
Powered by ViewVC 1.1.26