summaryrefslogtreecommitdiff
path: root/lib/libm/ieee/support.c
diff options
context:
space:
mode:
Diffstat (limited to 'lib/libm/ieee/support.c')
-rw-r--r--lib/libm/ieee/support.c106
1 files changed, 53 insertions, 53 deletions
diff --git a/lib/libm/ieee/support.c b/lib/libm/ieee/support.c
index e97683942188..0ae10bb1e4c8 100644
--- a/lib/libm/ieee/support.c
+++ b/lib/libm/ieee/support.c
@@ -35,37 +35,37 @@
static char sccsid[] = "@(#)support.c 8.1 (Berkeley) 6/4/93";
#endif /* not lint */
-/*
- * Some IEEE standard 754 recommended functions and remainder and sqrt for
+/*
+ * Some IEEE standard 754 recommended functions and remainder and sqrt for
* supporting the C elementary functions.
******************************************************************************
* WARNING:
* These codes are developed (in double) to support the C elementary
* functions temporarily. They are not universal, and some of them are very
- * slow (in particular, drem and sqrt is extremely inefficient). Each
- * computer system should have its implementation of these functions using
+ * slow (in particular, drem and sqrt is extremely inefficient). Each
+ * computer system should have its implementation of these functions using
* its own assembler.
******************************************************************************
*
* IEEE 754 required operations:
- * drem(x,p)
+ * drem(x,p)
* returns x REM y = x - [x/y]*y , where [x/y] is the integer
* nearest x/y; in half way case, choose the even one.
- * sqrt(x)
- * returns the square root of x correctly rounded according to
+ * sqrt(x)
+ * returns the square root of x correctly rounded according to
* the rounding mod.
*
* IEEE 754 recommended functions:
- * (a) copysign(x,y)
- * returns x with the sign of y.
- * (b) scalb(x,N)
+ * (a) copysign(x,y)
+ * returns x with the sign of y.
+ * (b) scalb(x,N)
* returns x * (2**N), for integer values N.
- * (c) logb(x)
- * returns the unbiased exponent of x, a signed integer in
- * double precision, except that logb(0) is -INF, logb(INF)
+ * (c) logb(x)
+ * returns the unbiased exponent of x, a signed integer in
+ * double precision, except that logb(0) is -INF, logb(INF)
* is +INF, and logb(NAN) is that NAN.
- * (d) finite(x)
- * returns the value TRUE if -INF < x < +INF and returns
+ * (d) finite(x)
+ * returns the value TRUE if -INF < x < +INF and returns
* FALSE otherwise.
*
*
@@ -78,7 +78,7 @@ static char sccsid[] = "@(#)support.c 8.1 (Berkeley) 6/4/93";
#if defined(vax)||defined(tahoe) /* VAX D format */
#include <errno.h>
static const unsigned short msign=0x7fff , mexp =0x7f80 ;
- static const short prep1=57, gap=7, bias=129 ;
+ static const short prep1=57, gap=7, bias=129 ;
static const double novf=1.7E38, nunf=3.0E-39, zero=0.0 ;
#else /* defined(vax)||defined(tahoe) */
static const unsigned short msign=0x7fff, mexp =0x7ff0 ;
@@ -97,7 +97,7 @@ double x; int N;
unsigned short *px=(unsigned short *) &x;
#endif /* national */
- if( x == zero ) return(x);
+ if( x == zero ) return(x);
#if defined(vax)||defined(tahoe)
if( (k= *px & mexp ) != ~msign ) {
@@ -117,7 +117,7 @@ double x; int N;
if( k < (mexp>>gap) ) *px = (*px&~mexp) | (k<<gap);
else x=novf+novf; /* overflow */
else
- if( k > -prep1 )
+ if( k > -prep1 )
/* gradual underflow */
{*px=(*px&~mexp)|(short)(1<<gap); x *= scalb(1.0,k-1);}
else
@@ -147,7 +147,7 @@ double x,y;
}
double logb(x)
-double x;
+double x;
{
#ifdef national
@@ -164,8 +164,8 @@ double x;
return ( (k>>gap) - bias );
else if( x != zero)
return ( -1022.0 );
- else
- return(-(1.0/zero));
+ else
+ return(-(1.0/zero));
else if(x != x)
return(x);
else
@@ -174,7 +174,7 @@ double x;
}
finite(x)
-double x;
+double x;
{
#if defined(vax)||defined(tahoe)
return(1);
@@ -192,16 +192,16 @@ double x,p;
{
short sign;
double hp,dp,tmp;
- unsigned short k;
+ unsigned short k;
#ifdef national
unsigned short
- *px=(unsigned short *) &x +3,
+ *px=(unsigned short *) &x +3,
*pp=(unsigned short *) &p +3,
*pd=(unsigned short *) &dp +3,
*pt=(unsigned short *) &tmp+3;
#else /* national */
unsigned short
- *px=(unsigned short *) &x ,
+ *px=(unsigned short *) &x ,
*pp=(unsigned short *) &p ,
*pd=(unsigned short *) &dp ,
*pt=(unsigned short *) &tmp;
@@ -230,13 +230,13 @@ double x,p;
#endif /* defined(vax)||defined(tahoe) */
{ if (p != p) return p; else return x;}
- else if ( ((*pp & mexp)>>gap) <= 1 )
+ else if ( ((*pp & mexp)>>gap) <= 1 )
/* subnormal p, or almost subnormal p */
{ double b; b=scalb(1.0,(int)prep1);
p *= b; x = drem(x,p); x *= b; return(drem(x,p)/b);}
else if ( p >= novf/2)
{ p /= 2 ; x /= 2; return(drem(x,p)*2);}
- else
+ else
{
dp=p+p; hp=p/2;
sign= *px & ~msign ;
@@ -294,13 +294,13 @@ double x;
}
/* sqrt(INF) is INF */
- if(!finite(x)) return(x);
+ if(!finite(x)) return(x);
/* scale x to [1,4) */
n=logb(x);
x=scalb(x,-n);
if((m=logb(x))!=0) x=scalb(x,-m); /* subnormal number */
- m += n;
+ m += n;
n = m/2;
if((n+n)!=m) {x *= 2; m -=1; n=m/2;}
@@ -313,23 +313,23 @@ double x;
else
s *= 2;
}
-
+
/* generate the last bit and determine the final rounding */
- r/=2; x *= 4;
+ r/=2; x *= 4;
if(x==zero) goto end; 100+r; /* trigger inexact flag */
if(s<x) {
q+=r; x -=s; s += 2; s *= 2; x *= 4;
- t = (x-s)-5;
+ t = (x-s)-5;
b=1.0+3*r/4; if(b==1.0) goto end; /* b==1 : Round-to-zero */
b=1.0+r/4; if(b>1.0) t=1; /* b>1 : Round-to-(+INF) */
if(t>=0) q+=r; } /* else: Round-to-nearest */
- else {
- s *= 2; x *= 4;
- t = (x-s)-1;
+ else {
+ s *= 2; x *= 4;
+ t = (x-s)-1;
b=1.0+3*r/4; if(b==1.0) goto end;
b=1.0+r/4; if(b>1.0) t=1;
if(t>=0) q+=r; }
-
+
end: return(scalb(q,n));
}
@@ -342,13 +342,13 @@ end: return(scalb(q,n));
*
* Warning: this code should not get compiled in unless ALL of
* the following machine-dependent routines are supplied.
- *
+ *
* Required machine dependent functions (not on a VAX):
* swapINX(i): save inexact flag and reset it to "i"
* swapENI(e): save inexact enable and reset it to "e"
*/
-double drem(x,y)
+double drem(x,y)
double x,y;
{
@@ -363,8 +363,8 @@ double x,y;
double hy,y1,t,t1;
short k;
long n;
- int i,e;
- unsigned short xexp,yexp, *px =(unsigned short *) &x ,
+ int i,e;
+ unsigned short xexp,yexp, *px =(unsigned short *) &x ,
nx,nf, *py =(unsigned short *) &y ,
sign, *pt =(unsigned short *) &t ,
*pt1 =(unsigned short *) &t1 ;
@@ -381,7 +381,7 @@ double x,y;
/* save the inexact flag and inexact enable in i and e respectively
* and reset them to zero
*/
- i=swapINX(0); e=swapENI(0);
+ i=swapINX(0); e=swapENI(0);
/* subnormal number */
nx=0;
@@ -391,7 +391,7 @@ double x,y;
if( yexp <= m57 ) {py[n0]+=m57; nx+=m57; yexp+=m57;}
nf=nx;
- py[n0] &= 0x7fff;
+ py[n0] &= 0x7fff;
px[n0] &= 0x7fff;
/* mask off the least significant 27 bits of y */
@@ -408,7 +408,7 @@ loop:
if(k>0) /* if x/y >= 2**26, scale up y so that x/y < 2**26 */
{pt[n0]+=k;pt1[n0]+=k;}
n=x/t; x=(x-n*t1)-n*(t-t1);
- }
+ }
/* end while (x > y) */
if(nx!=0) {t=1.0; pt[n0]+=nx; x*=t; nx=0; goto loop;}
@@ -416,14 +416,14 @@ loop:
/* final adjustment */
hy=y/2.0;
- if(x>hy||((x==hy)&&n%2==1)) x-=y;
+ if(x>hy||((x==hy)&&n%2==1)) x-=y;
px[n0] ^= sign;
if(nf!=0) { t=1.0; pt[n0]-=nf; x*=t;}
/* restore inexact flag and inexact enable */
- swapINX(i); swapENI(e);
+ swapINX(i); swapENI(e);
- return(x);
+ return(x);
}
#endif
@@ -435,7 +435,7 @@ loop:
*
* Warning: this code should not get compiled in unless ALL of
* the following machine-dependent routines are supplied.
- *
+ *
* Required machine dependent functions:
* swapINX(i) ...return the status of INEXACT flag and reset it to "i"
* swapRM(r) ...return the current Rounding Mode and reset it to "r"
@@ -456,16 +456,16 @@ double x;
double const b54=134217728.*134217728.; /* b54=2**54 */
long mx,scalx;
long const mexp=0x7ff00000;
- int i,j,r,e,swapINX(),swapRM(),swapENI();
+ int i,j,r,e,swapINX(),swapRM(),swapENI();
unsigned long *py=(unsigned long *) &y ,
*pt=(unsigned long *) &t ,
*px=(unsigned long *) &x ;
#ifdef national /* ordering of word in a floating point number */
- const int n0=1, n1=0;
+ const int n0=1, n1=0;
#else
- const int n0=0, n1=1;
+ const int n0=0, n1=1;
#endif
-/* Rounding Mode: RN ...round-to-nearest
+/* Rounding Mode: RN ...round-to-nearest
* RZ ...round-towards 0
* RP ...round-towards +INF
* RM ...round-towards -INF
@@ -502,10 +502,10 @@ double x;
t=x/y; y=y+t; py[n0]=py[n0]-0x00100006; py[n1]=0;
/* triple to almost 56 sig. bits; now y approx. sqrt(x) to within 1 ulp */
- t=y*y; z=t; pt[n0]+=0x00100000; t+=z; z=(x-z)*y;
+ t=y*y; z=t; pt[n0]+=0x00100000; t+=z; z=(x-z)*y;
t=z/(t+x) ; pt[n0]+=0x00100000; y+=t;
-/* twiddle last bit to force y correctly rounded */
+/* twiddle last bit to force y correctly rounded */
swapRM(RZ); /* ...set Rounding Mode to round-toward-zero */
swapINX(0); /* ...clear INEXACT flag */
swapENI(e); /* ...restore inexact enable status */
generated by cgit v1.2.3 (git 2.25.1) at 2025-10-31 14:04:10 +0000