#if !_UWIN || (_lib__copysign||_lib_copysign) && _lib_logb && (_lib__finite||_lib_finite) && (_lib_drem||_lib_remainder) && _lib_sqrt && _lib_ilogb && (_lib__scalb||_lib_scalb)
void _STUB_support(){}
#else
/*
* Copyright (c) 1985, 1993
* The Regents of the University of California. All rights reserved.
*
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions
* are met:
* 1. Redistributions of source code must retain the above copyright
* notice, this list of conditions and the following disclaimer.
* 2. Redistributions in binary form must reproduce the above copyright
* notice, this list of conditions and the following disclaimer in the
* 3. Neither the name of the University nor the names of its contributors
* may be used to endorse or promote products derived from this software
* without specific prior written permission.
*
* THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
* ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
* IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
* ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
* FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
* DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
* OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
* HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
* LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
* OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
* SUCH DAMAGE.
*/
#ifndef lint
#endif /* not lint */
/*
* 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
* its own assembler.
******************************************************************************
*
* IEEE 754 required operations:
* 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
* the rounding mod.
*
* IEEE 754 recommended functions:
* (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)
* is +INF, and logb(NAN) is that NAN.
* (d) finite(x)
* returns the value TRUE if -INF < x < +INF and returns
* FALSE otherwise.
*
*
* CODED IN C BY K.C. NG, 11/25/84;
* REVISED BY K.C. NG on 1/22/85, 2/13/85, 3/24/85.
*/
#include "mathimpl.h"
#include <errno.h>
#else /* defined(vax)||defined(tahoe) */
#endif /* defined(vax)||defined(tahoe) */
#if !_lib__scalb || !_lib_scalb
extern double _scalb(x,N)
double x; double N;
{
int k;
#ifdef national
#else /* national */
unsigned short *px=(unsigned short *) &x;
#endif /* national */
if( x == zero ) return(x);
if (N < -260)
else if (N > 260) {
}
#else /* defined(vax)||defined(tahoe) */
if( k == 0 ) {
#endif /* defined(vax)||defined(tahoe) */
if((k = (k>>gap)+ N) > 0 )
else
if( k > -prep1 )
/* gradual underflow */
else
}
return(x);
}
#endif
#if !_lib_scalb
extern double scalb(x,N)
double x; double N;
{
return _scalb(x, N);
}
#endif
#if !_lib__copysign
extern double _copysign(x,y)
double x,y;
{
#ifdef national
#else /* national */
unsigned short *px=(unsigned short *) &x,
*py=(unsigned short *) &y;
#endif /* national */
#endif /* defined(vax)||defined(tahoe) */
return(x);
}
#endif
#if !_lib_copysign
extern double copysign(x,y)
double x,y;
{
return _copysign(x,y);
}
#endif
#if !_lib_logb
extern double logb(x)
double x;
{
#ifdef national
#else /* national */
short *px=(short *) &x, k;
#endif /* national */
#else /* defined(vax)||defined(tahoe) */
if ( k != 0 )
else if( x != zero)
return ( -1022.0 );
else
return(-(1.0/zero));
else if(x != x)
return(x);
else
#endif /* defined(vax)||defined(tahoe) */
}
#endif
#if !_lib__finite
extern int _finite(x)
double x;
{
return(1);
#else /* defined(vax)||defined(tahoe) */
#ifdef national
#else /* national */
#endif /* national */
#endif /* defined(vax)||defined(tahoe) */
}
#endif
#if !_lib_finite
extern int finite(x)
double x;
{
return _finite(x);
}
#endif
#if !_lib_drem
extern double drem(x,p)
double x,p;
{
#if _lib_remainder
return remainder(x,p);
#else
short sign;
unsigned short k;
#ifdef national
unsigned short
*px=(unsigned short *) &x +3,
*pp=(unsigned short *) &p +3,
#else /* national */
unsigned short
*px=(unsigned short *) &x ,
*pp=(unsigned short *) &p ,
#endif /* national */
#else /* defined(vax)||defined(tahoe) */
#endif /* defined(vax)||defined(tahoe) */
return (x-p)-(x-p); /* create nan if x is inf */
if (p == zero) {
#else /* defined(vax)||defined(tahoe) */
#endif /* defined(vax)||defined(tahoe) */
}
#else /* defined(vax)||defined(tahoe) */
#endif /* defined(vax)||defined(tahoe) */
{ if (p != p) return p; else return x;}
/* subnormal p, or almost subnormal p */
else if ( p >= novf/2)
else
{
while ( x > dp )
{
*pt += k ;
#else /* defined(vax)||defined(tahoe) */
#endif /* defined(vax)||defined(tahoe) */
x -= tmp ;
}
if ( x > hp )
{ x -= p ; if ( x >= hp ) x -= p ; }
if (x)
#endif /* defined(vax)||defined(tahoe) */
return( x);
}
#endif
}
#endif
#if !_lib_remainder
extern double remainder(x,p)
double x,p;
{
return drem(x,p);
}
#endif
#if !_lib_sqrt
extern double sqrt(x)
double x;
{
double q,s,b,r;
double t;
int m,n,i;
int k=54;
#else /* defined(vax)||defined(tahoe) */
int k=51;
#endif /* defined(vax)||defined(tahoe) */
/* sqrt(NaN) is NaN, sqrt(+-0) = +-0 */
if(x!=x||x==zero) return(x);
/* sqrt(negative) is invalid */
if(x<zero) {
#else /* defined(vax)||defined(tahoe) */
#endif /* defined(vax)||defined(tahoe) */
}
/* sqrt(INF) is INF */
if(!finite(x)) return(x);
/* scale x to [1,4) */
n=logb(x);
x=scalb(x,-n);
m += n;
n = m/2;
if((n+n)!=m) {x *= 2; m -=1; n=m/2;}
/* generate sqrt(x) bit by bit (accumulating in q) */
q=1.0; s=4.0; x -= 1.0; r=1;
for(i=1;i<=k;i++) {
t=s+1; x *= 4; r /= 2;
if(t<=x) {
s=t+t+2, x -= t; q += r;}
else
s *= 2;
}
/* generate the last bit and determine the final rounding */
r/=2; x *= 4;
if(s<x) {
q+=r; x -=s; s += 2; s *= 2; x *= 4;
t = (x-s)-5;
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;
b=1.0+r/4; if(b>1.0) t=1;
if(t>=0) q+=r; }
}
#endif
#if 0
/* DREM(X,Y)
* RETURN X REM Y =X-N*Y, N=[X/Y] ROUNDED (ROUNDED TO EVEN IN THE HALF WAY CASE)
* DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS)
* INTENDED FOR ASSEMBLY LANGUAGE
* CODED IN C BY K.C. NG, 3/23/85, 4/8/85.
*
* 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"
*/
extern double drem(x,y)
double x,y;
{
#ifdef national /* order of words in floating point number */
#else /* VAX, SUN, ZILOG, TAHOE */
#endif
static const double zero=0.0;
short k;
long n;
int i,e;
/* return NaN if x is NaN, or y is NaN, or x is INF, or y is zero */
if(x!=x) return(x); if(y!=y) return(y); /* x or y is NaN */
if(y==zero) return(y/y);
/* save the inexact flag and inexact enable in i and e respectively
* and reset them to zero
*/
/* subnormal number */
nx=0;
/* if y is tiny (biased exponent <= 57), scale up y to y*2**57 */
/* mask off the least significant 27 bits of y */
/* LOOP: argument reduction on x whenever x > y */
loop:
while ( x > y )
{
t=y;
if(k>0) /* if x/y >= 2**26, scale up y so that x/y < 2**26 */
}
/* end while (x > y) */
/* final adjustment */
hy=y/2.0;
/* restore inexact flag and inexact enable */
return(x);
}
#endif
#if 0
/* SQRT
* RETURN CORRECTLY ROUNDED (ACCORDING TO THE ROUNDING MODE) SQRT
* FOR IEEE DOUBLE PRECISION ONLY, INTENDED FOR ASSEMBLY LANGUAGE
* CODED IN C BY K.C. NG, 3/22/85.
*
* 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"
* swapENI(e) ...return the status of inexact enable and reset it to "e"
* addc(t) ...perform t=t+1 regarding t as a 64 bit unsigned integer
* subc(t) ...perform t=t-1 regarding t as a 64 bit unsigned integer
*/
static const unsigned long table[] = {
0, 1204, 3062, 5746, 9193, 13348, 18162, 23592, 29598, 36145, 43202, 50740,
58733, 67158, 75992, 85215, 83599, 71378, 60428, 50647, 41945, 34246, 27478,
21581, 16499, 12183, 8588, 5674, 3403, 1742, 661, 130, };
extern double newsqrt(x)
double x;
{
long const mexp=0x7ff00000;
unsigned long *py=(unsigned long *) &y ,
*pt=(unsigned long *) &t ,
*px=(unsigned long *) &x ;
#ifdef national /* ordering of word in a floating point number */
#else
#endif
/* Rounding Mode: RN ...round-to-nearest
* RZ ...round-towards 0
* RP ...round-towards +INF
* RM ...round-towards -INF
*/
/* machine dependent: work on a Zilog Z8070
* and a National 32081 & 16081
*/
/* exceptions */
if(x!=x||x==0.0) return(x); /* sqrt(NaN) is NaN, sqrt(+-0) = +-0 */
if(x<0) return((x-x)/(x-x)); /* sqrt(negative) is invalid */
/* save, reset, initialize */
e=swapENI(0); /* ...save and reset the inexact enable */
i=swapINX(0); /* ...save INEXACT flag */
scalx=0;
/* subnormal number, scale up x to x*2**54 */
* if (x > 2**512) x=x/2**512; if (x < 2**-512) x=x*2**512 */
/* magic initial approximation to almost 8 sig. bits */
/* Heron's rule once with correction to improve y to almost 18 sig. bits */
/* triple to almost 56 sig. bits; now y approx. sqrt(x) to within 1 ulp */
/* twiddle last bit to force y correctly rounded */
swapINX(0); /* ...clear INEXACT flag */
swapENI(e); /* ...restore inexact enable status */
t=x/y; /* ...chopped quotient, possibly inexact */
j=swapINX(i); /* ...read and restore inexact flag */
y=y+t; /* ...chopped sum */
swapRM(r); /* ...restore Rounding Mode */
return(y);
}
#endif
#if !_lib_ilogb
extern int ilogb(double x)
{
return((int)logb(x));
}
#endif
#endif