#include "FEATURE/uwin"
#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
#ifndef lint
static char sccsid[] = "@(#)support.c 8.1 (Berkeley) 6/4/93";
#endif
#include "mathimpl.h"
#if defined(vax)||defined(tahoe)
#include <errno.h>
static const unsigned short msign=0x7fff , mexp =0x7f80 ;
static const short prep1=57, gap=7, bias=129 ;
static const double novf=1.7E38, nunf=3.0E-39, zero=0.0 ;
#else
static const unsigned short msign=0x7fff, mexp =0x7ff0 ;
static const short prep1=54, gap=4, bias=1023 ;
static const double novf=1.7E308, nunf=3.0E-308,zero=0.0;
#endif
#if !_lib__scalb || !_lib_scalb
extern double _scalb(x,N)
double x; double N;
{
int k;
#ifdef national
unsigned short *px=(unsigned short *) &x + 3;
#else
unsigned short *px=(unsigned short *) &x;
#endif
if( x == zero ) return(x);
#if defined(vax)||defined(tahoe)
if( (k= *px & mexp ) != ~msign ) {
if (N < -260)
return(nunf*nunf);
else if (N > 260) {
return(copysign(infnan(ERANGE),x));
}
#else
if( (k= *px & mexp ) != mexp ) {
if( N<-2100) return(nunf*nunf); else if(N>2100) return(novf+novf);
if( k == 0 ) {
x *= scalb(1.0,prep1); N -= prep1; return(scalb(x,N));}
#endif
if((k = (k>>gap)+ N) > 0 )
if( k < (mexp>>gap) ) *px = (*px&~mexp) | (k<<gap);
else x=novf+novf;
else
if( k > -prep1 )
{*px=(*px&~mexp)|(short)(1<<gap); x *= scalb(1.0,k-1);}
else
return(nunf*nunf);
}
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
unsigned short *px=(unsigned short *) &x+3,
*py=(unsigned short *) &y+3;
#else
unsigned short *px=(unsigned short *) &x,
*py=(unsigned short *) &y;
#endif
#if defined(vax)||defined(tahoe)
if ( (*px & mexp) == 0 ) return(x);
#endif
*px = ( *px & msign ) | ( *py & ~msign );
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
short *px=(short *) &x+3, k;
#else
short *px=(short *) &x, k;
#endif
#if defined(vax)||defined(tahoe)
return (int)(((*px&mexp)>>gap)-bias);
#else
if( (k= *px & mexp ) != mexp )
if ( k != 0 )
return ( (k>>gap) - bias );
else if( x != zero)
return ( -1022.0 );
else
return(-(1.0/zero));
else if(x != x)
return(x);
else
{*px &= msign; return(x);}
#endif
}
#endif
#if !_lib__finite
extern int _finite(x)
double x;
{
#if defined(vax)||defined(tahoe)
return(1);
#else
#ifdef national
return( (*((short *) &x+3 ) & mexp ) != mexp );
#else
return( (*((short *) &x ) & mexp ) != mexp );
#endif
#endif
}
#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;
double hp,dp,tmp;
unsigned short k;
#ifdef national
unsigned short
*px=(unsigned short *) &x +3,
*pp=(unsigned short *) &p +3,
*pd=(unsigned short *) &dp +3,
*pt=(unsigned short *) &tmp+3;
#else
unsigned short
*px=(unsigned short *) &x ,
*pp=(unsigned short *) &p ,
*pd=(unsigned short *) &dp ,
*pt=(unsigned short *) &tmp;
#endif
*pp &= msign ;
#if defined(vax)||defined(tahoe)
if( ( *px & mexp ) == ~msign )
#else
if( ( *px & mexp ) == mexp )
#endif
return (x-p)-(x-p);
if (p == zero) {
#if defined(vax)||defined(tahoe)
return(infnan(EDOM));
#else
return zero/zero;
#endif
}
#if defined(vax)||defined(tahoe)
if( ( *pp & mexp ) == ~msign )
#else
if( ( *pp & mexp ) == mexp )
#endif
{ if (p != p) return p; else return x;}
else if ( ((*pp & mexp)>>gap) <= 1 )
{ 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
{
dp=p+p; hp=p/2;
sign= *px & ~msign ;
*px &= msign ;
while ( x > dp )
{
k=(*px & mexp) - (*pd & mexp) ;
tmp = dp ;
*pt += k ;
#if defined(vax)||defined(tahoe)
if( x < tmp ) *pt -= 128 ;
#else
if( x < tmp ) *pt -= 16 ;
#endif
x -= tmp ;
}
if ( x > hp )
{ x -= p ; if ( x >= hp ) x -= p ; }
#if defined(vax)||defined(tahoe)
if (x)
#endif
*px ^= sign;
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;
double const zero=0.0;
int m,n,i;
#if defined(vax)||defined(tahoe)
int k=54;
#else
int k=51;
#endif
if(x!=x||x==zero) return(x);
if(x<zero) {
#if defined(vax)||defined(tahoe)
return (infnan(EDOM));
#else
return(zero/zero);
#endif
}
if(!finite(x)) return(x);
n=logb(x);
x=scalb(x,-n);
if((m=logb(x))!=0) x=scalb(x,-m);
m += n;
n = m/2;
if((n+n)!=m) {x *= 2; m -=1; n=m/2;}
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;
}
r/=2; x *= 4;
if(x==zero) goto end; 100+r;
if(s<x) {
q+=r; x -=s; s += 2; s *= 2; x *= 4;
t = (x-s)-5;
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; }
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));
}
#endif
#if 0
extern double drem(x,y)
double x,y;
{
#ifdef national
static const n0=3,n1=2,n2=1,n3=0;
#else
static const n0=0,n1=1,n2=2,n3=3;
#endif
static const unsigned short mexp =0x7ff0, m25 =0x0190, m57 =0x0390;
static const double zero=0.0;
double hy,y1,t,t1;
short k;
long n;
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 ;
xexp = px[n0] & mexp ;
yexp = py[n0] & mexp ;
sign = px[n0] &0x8000;
if(x!=x) return(x); if(y!=y) return(y);
if( xexp == mexp ) return(zero/zero);
if(y==zero) return(y/y);
i=swapINX(0); e=swapENI(0);
nx=0;
if(yexp==0) {t=1.0,pt[n0]+=m57; y*=t; nx=m57;}
if( yexp <= m57 ) {py[n0]+=m57; nx+=m57; yexp+=m57;}
nf=nx;
py[n0] &= 0x7fff;
px[n0] &= 0x7fff;
t=y; pt[n3]=0; pt[n2]&=0xf800; y1=t;
loop:
while ( x > y )
{
t=y;
t1=y1;
xexp=px[n0]&mexp;
k=xexp-yexp-m25;
if(k>0)
{pt[n0]+=k;pt1[n0]+=k;}
n=x/t; x=(x-n*t1)-n*(t-t1);
}
if(nx!=0) {t=1.0; pt[n0]+=nx; x*=t; nx=0; goto loop;}
hy=y/2.0;
if(x>hy||((x==hy)&&n%2==1)) x-=y;
px[n0] ^= sign;
if(nf!=0) { t=1.0; pt[n0]-=nf; x*=t;}
swapINX(i); swapENI(e);
return(x);
}
#endif
#if 0
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;
{
double y,z,t,addc(),subc()
double const b54=134217728.*134217728.;
long mx,scalx;
long const mexp=0x7ff00000;
int i,j,r,e,swapINX(),swapRM(),swapENI();
unsigned long *py=(unsigned long *) &y ,
*pt=(unsigned long *) &t ,
*px=(unsigned long *) &x ;
#ifdef national
const int n0=1, n1=0;
#else
const int n0=0, n1=1;
#endif
const int RN=0,RZ=1,RP=2,RM=3;
if(x!=x||x==0.0) return(x);
if(x<0) return((x-x)/(x-x));
if((mx=px[n0]&mexp)==mexp) return(x);
e=swapENI(0);
i=swapINX(0);
r=swapRM(RN);
scalx=0;
if(mx==0) {x *= b54 ; scalx-=0x01b00000;}
if(mx>0x5ff00000) {px[n0] -= 0x20000000; scalx+= 0x10000000;}
if(mx<0x1ff00000) {px[n0] += 0x20000000; scalx-= 0x10000000;}
py[n0]=(px[n0]>>1)+0x1ff80000;
py[n0]=py[n0]-table[(py[n0]>>15)&31];
t=x/y; y=y+t; py[n0]=py[n0]-0x00100006; py[n1]=0;
t=y*y; z=t; pt[n0]+=0x00100000; t+=z; z=(x-z)*y;
t=z/(t+x) ; pt[n0]+=0x00100000; y+=t;
swapRM(RZ);
swapINX(0);
swapENI(e);
t=x/y;
j=swapINX(i);
if(j==0) { if(t==y) goto end; else t=subc(t); }
b54+0.1;
if(r==RN) t=addc(t);
else if(r==RP) { t=addc(t);y=addc(y);}
y=y+t;
py[n0]=py[n0]-0x00100000;
end: py[n0]=py[n0]+scalx;
swapRM(r);
return(y);
}
#endif
#if !_lib_ilogb
extern int ilogb(double x)
{
return((int)logb(x));
}
#endif
#endif