114 lines
3.3 KiB
C
114 lines
3.3 KiB
C
|
|
#ifndef lint
|
|
static char sccsid[] = "@(#)fmod.c 1.1 92/07/30 SMI";
|
|
#endif
|
|
|
|
/*
|
|
* Copyright (c) 1988 by Sun Microsystems, Inc.
|
|
*/
|
|
|
|
#include <math.h>
|
|
|
|
static long
|
|
is = 0x80000000,
|
|
im = 0x000fffff,
|
|
iu = 0x00100000;
|
|
|
|
static double
|
|
zero = 0.0,
|
|
two110 = 1.29807421463370690713e+33, /* 46D00000 00000000 */
|
|
twom110 = 7.70371977754894341222e-34, /* 39100000 00000000 */
|
|
twom969 = 2.00416836000897277800e-292; /* 03600000 00000000 */
|
|
|
|
double fmod(x,y)
|
|
double x,y ;
|
|
{
|
|
double a,b;
|
|
long *px = (long*)&x;
|
|
long *pa = (long*)&a;
|
|
long *pb = (long*)&b;
|
|
int n0,n1,n,k,ia,ib;
|
|
int x0,y0,z0;
|
|
unsigned x1,y1,z1;
|
|
|
|
/* purge off exception values */
|
|
if(!finite(x)||y!=y||y==zero) return (x*y)/(x*y);
|
|
a = fabs(x); b = fabs(y);
|
|
if(a<=b) { if(a<b) return x; else return zero*x;}
|
|
if(b<twom969) {
|
|
b *= two110;
|
|
return twom110*fmod(two110*fmod(x,b),b);
|
|
}
|
|
|
|
/* check out ordering of floating point number */
|
|
if ((*(int*)&two110)==0x46d00000) {n0=0;n1=1;} else {n0=1;n1=0;}
|
|
|
|
/* a,b is in reasonable range. Now prepare for fix point arithmetic */
|
|
ia = pa[n0]&0x7ff00000; ib = pb[n0]&0x7ff00000;
|
|
n = ((ia-ib)>>20);
|
|
x0 = iu|(im&pa[n0]); x1 = pa[n1];
|
|
y0 = iu|(im&pb[n0]); y1 = pb[n1];
|
|
|
|
/* fix point fmod */
|
|
/*
|
|
while(n--) {
|
|
z0=x0-y0;z1=x1-y1; if(x1<y1) z0 -= 1;
|
|
if(z0<0){x0 = x0+x0+((x1&is)!=0); x1 = x1+x1;}
|
|
else {
|
|
if((z0|z1)==0){pa[n0]=px[n0]&is;pa[n1]=0;return a;}
|
|
x0 = z0+z0+((z1&is)!=0); x1 = z1+z1;
|
|
}
|
|
}
|
|
*/
|
|
/* unroll the above loop 4 times to gain performance */
|
|
k = n&3; n >>= 2;
|
|
while(n--) {
|
|
z0=x0-y0;z1=x1-y1; if(x1<y1) z0 -= 1;
|
|
if(z0<0){x0 = x0+x0+((x1&is)!=0); x1 = x1+x1;}
|
|
else { if((z0|z1)==0){pa[n0]=px[n0]&is;pa[n1]=0;return a;}
|
|
x0 = z0+z0+((z1&is)!=0); x1 = z1+z1; }
|
|
z0=x0-y0;z1=x1-y1; if(x1<y1) z0 -= 1;
|
|
if(z0<0){x0 = x0+x0+((x1&is)!=0); x1 = x1+x1;}
|
|
else { if((z0|z1)==0){pa[n0]=px[n0]&is;pa[n1]=0;return a;}
|
|
x0 = z0+z0+((z1&is)!=0); x1 = z1+z1; }
|
|
z0=x0-y0;z1=x1-y1; if(x1<y1) z0 -= 1;
|
|
if(z0<0){x0 = x0+x0+((x1&is)!=0); x1 = x1+x1;}
|
|
else { if((z0|z1)==0){pa[n0]=px[n0]&is;pa[n1]=0;return a;}
|
|
x0 = z0+z0+((z1&is)!=0); x1 = z1+z1; }
|
|
z0=x0-y0;z1=x1-y1; if(x1<y1) z0 -= 1;
|
|
if(z0<0){x0 = x0+x0+((x1&is)!=0); x1 = x1+x1;}
|
|
else { if((z0|z1)==0){pa[n0]=px[n0]&is;pa[n1]=0;return a;}
|
|
x0 = z0+z0+((z1&is)!=0); x1 = z1+z1; }
|
|
}
|
|
switch(k) {
|
|
case 3:
|
|
z0=x0-y0;z1=x1-y1; if(x1<y1) z0 -= 1;
|
|
if(z0<0){x0 = x0+x0+((x1&is)!=0); x1 = x1+x1;}
|
|
else { if((z0|z1)==0){pa[n0]=px[n0]&is;pa[n1]=0;return a;}
|
|
x0 = z0+z0+((z1&is)!=0); x1 = z1+z1; }
|
|
case 2:
|
|
z0=x0-y0;z1=x1-y1; if(x1<y1) z0 -= 1;
|
|
if(z0<0){x0 = x0+x0+((x1&is)!=0); x1 = x1+x1;}
|
|
else { if((z0|z1)==0){pa[n0]=px[n0]&is;pa[n1]=0;return a;}
|
|
x0 = z0+z0+((z1&is)!=0); x1 = z1+z1; }
|
|
case 1:
|
|
z0=x0-y0;z1=x1-y1; if(x1<y1) z0 -= 1;
|
|
if(z0<0){x0 = x0+x0+((x1&is)!=0); x1 = x1+x1;}
|
|
else { if((z0|z1)==0){pa[n0]=px[n0]&is;pa[n1]=0;return a;}
|
|
x0 = z0+z0+((z1&is)!=0); x1 = z1+z1; }
|
|
}
|
|
z0=x0-y0;z1=x1-y1; if(x1<y1) z0 -= 1;
|
|
if(z0>=0) {x0=z0;x1=z1;}
|
|
|
|
/* convert back to floating value and restore the sign */
|
|
if((x0|x1)==0) {pa[n0]=px[n0]&is;pa[n1]=0;return a;}
|
|
while((x0&iu)==0) {
|
|
x0 = x0+x0+((x1&is)!=0); x1 = x1+x1;
|
|
ib -= iu;
|
|
}
|
|
ib |= px[n0]&is;
|
|
pa[n0] = ib|(x0&im);
|
|
pa[n1] = x1;
|
|
return a;
|
|
}
|