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 |
{ |
{ |
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 |
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", |