-ARCHIVE- asin.c 560
/*	asin and acos - quick and dirty;
*/

#include "errno.h"

static double halfpi      = 1.570796326794896619;

double asin(x)
double x;
{
  int sign;
  double temp,atan(),sqrt();

  if(x<0.0){
    x=-x;
    sign=-1;
  }else sign=0;
  if(x>1.0){
    errno=EDOM;
    return 0.0;
  }
  temp=sqrt(1.0-x*x);
  if(x>0.7)temp=halfpi-atan(temp/x);
  else temp=atan(x/temp);
  if(sign)temp=-temp;
  return temp;
}

double acos(x)
double x;
{

  if((x>1.0)||(x<-1.0)){
    errno=EDOM;
    return 0.0;
  }
  return halfpi-asin(x);
}
-ARCHIVE- atan.c 1149
/*	atan and atan2
*/

static double a[]={
0.0,
0.523598775598298,
1.570796326794896,
1.047197551196597
};

static double _atan(f)
double f;
{
  int n;
  double g,gpg,qg,r,fabs();

  if(f>1.0){
    f=1.0/f;
    n=2;
  } else n=0;
  if(f>0.267949192431122){
    g=0.732050807568877*f-0.5;
    g-=0.5;
    f=(g+f)/(f+1.732050807568877);
    ++n;
  }
  if(fabs(f)>1.0e-8){
    g=f*f;
    gpg=(((-0.837582993681500*g-0.849462403513206e1)*g-0.205058551958616e2)
*g-0.136887688941919e2)*g;
    qg=(((g+0.150240011600285e2)*g+0.595784361425973e2)
*g+0.861573495971302e2)*g+0.410663066825757e2;
    r=gpg/qg;
    f+=f*r;
  }
  if(n>1)f=-f;
  f+=a[n];
  return f;
}

double atan(x)
double x;
{
  double r;

  r=_atan(fabs(x));
  if(x<0.0)r=-r;
  return r;
}

double atan2(v,u)
double v,u;
{
  int j,*ip;
  double r,fabs();

  r=a[2];
  if(u==0.0){
    if(v==0.0)return 0.0;
    goto C;
  }
  ip=&v;
  j=(*ip>>4)&0x3ff;
  ip=&u;
  j-=(*ip>>4)*0x3ff;
  if(j>0x3fc)goto C;
  if(j<-0x3fc)r=0.0;
  else r=_atan(fabs(v/u));
  if(u<0.0)r=3.141592653589793-r;
C:
  if(v<0.0)r=-r;
  return r;
}
-ARCHIVE- ceil.c 104
/*	return ceil
*/

double ceil(val)
double val;
{
  double floor();

  return -floor(-val);
}
-ARCHIVE- errno.c 53
/*	define the error control field
*/

int errno;
-ARCHIVE- errno.h 209
/*	error return definitions for math functions
*/

extern int errno;		/* contains the error indicator */

#define EDOM 33			/* arg not in domain of function */
#define ERANGE 34		/* result too large */
-ARCHIVE- exp.c 617
/*	exp function
*/

#include "errno.h"

double exp(x)
double x;
{
  int n;
  static double bigval=500.0;
  double temp,xn,g,z,gpz,qz,rg,floor(),ldexp();

  if(x==0.0)return 1.0;
  if(x<-bigval)return 0.0;
  if(x>bigval){
    errno=ERANGE;
    return 1.0e300;
  }
  temp=x*1.44269504088896;
  xn=floor(temp);
  if(temp-xn>0.5)xn+=1.0;
  n=xn;
  g=x-xn*0.693359375;
  g+=xn*2.12194440054690e-4;
  z=g*g;
  gpz=((0.165203300268279e-4*z+0.694360001511792e-2)*z+0.24999999999992)*g;
  qz=  (0.495862884905441e-3*z+0.555538666969001e-1)*z+0.5;
  rg=0.5+gpz/(qz-gpz);
  return ldexp(rg,n+1);
}
-ARCHIVE- fabs.c 113
/*	return the absolute value of a float
*/

double fabs(x)
double x;
{

  if(x<0.0)x=-x;
  return x;
}
-ARCHIVE- floor.c 201
/*	return integer <= val
*/

double floor(val)
double val;
{
  double modf();

  if(val<0.0){
    if(modf(-val,&val)!=0.0)val+=1.0;
    val=-val;
  } else modf(val,&val);
  return val;
}
-ARCHIVE- frexp.c 244
/*	return exponent and fraction
*/

double frexp(val,ip)
double val;
int *ip;
{
  int *tip;

  if(val==0.0)*ip=0;
  else {
    tip=&val;
    *ip=((tip[3]>>4)&0x7ff)-0x3fe;
    tip[3]=(tip[3]&0x800f)+0x3fe0;
  }
  return val;
}
-ARCHIVE- ldexp.c 352
/*	rebuild a floating point number
*/

#include "errno.h"

double ldexp(val,exp)
double val;
int exp;
{
  int *ip;

  ip=&val;
  exp+=(ip[3]>>4)&0x3ff;
  if(exp<1)val=0.0;		/* it underflowed */
  else if(exp>0x7ff) {
    ip[3]|=0x7fff;
    ip[2]=ip[1]=ip[0]=-1;
    errno=ERANGE;
  } ip[3]=(ip[3]&0x800f)|(exp<<4);
  return val;
}
-ARCHIVE- log.c 726
/*	log base e and base 10
*/

#include "errno.h"

double log(val)
double val;
{
  int n;
  double f,znum,zden,z,w,aw,bw,rz,frexp();

  if(val<=0.0){
    errno=EDOM;
    return -1.0e300;
  }
  f=frexp(val,&n);
  znum=f-0.5;
  if(f>0.707106781186547){
    znum-=0.5;
    zden=f*0.5+0.5;
  } else {
    n-=1;
    zden=znum*0.5+0.5;
  }
  z=znum/zden;
  w=z*z;
  bw=((w-0.356679777390346e2)*w+0.312032220919245e3)*w-0.769499321084948e3;
  aw=(  -0.789561128874912e0 *w+0.163839435630215e2)*w-0.641249434237455e2;
  rz=z+z*(w*aw/bw);
  val=355;
  val/=512;
  val*=n;
  val+=rz-n*2.1219444005469e-4;
  return val;
}

double log10(val)
double val;
{

   return log(val)*0.434294481903252;
}
-ARCHIVE- modf.c 464
/*	return the integer and fractional part of a number
*/

double modf(val,ivp)
double val;		/* the input value */
double *ivp;		/* where to put integer part */
{
  int exp,*ip;

  ip=&val;
  exp=((ip[3]>>4)&0x7ff)-0x3ff;
  if(exp<0)*ivp=0.0;	/* very small */
  else {
    *ivp=val;
    if(exp<52){
      exp=52-exp;
      for(ip=ivp;exp>0;exp-=16){
	if(exp>=16)*ip++=0;
        else *ip++&=(-1)<<exp;
      }
    }
  }
  return val-*ivp;
}
-ARCHIVE- pow.c 373
/*	returns x**y
*/

#include "errno.h"

double pow(x,y)
double x, y;
{
  long l;
  double temp,log(),exp();

  if(x<=0.0){
    if(x==0.0){
      if(y<=0.0)goto bad;
      return 0.0;
    }
    l=y;
    if(l!=y)goto bad;
    temp=exp(y*log(-x));
    if(l&1)temp=-temp;
  } else temp=exp(y*log(x));
  return temp;
bad:
  errno=EDOM;
  return 0.0;
}
-ARCHIVE- sin.c 1068
/*	sin/cos routines
*/

static double r[]={
  -0.16666666666666665052e0,
   0.83333333333331650314e-2,
  -0.19841269841201840457e-3,
   0.27557319210152756119e-5,
  -0.25052106798274584544e-7,
   0.16058936490371589114e-9,
  -0.76429178068910467734e-12,
   0.27204790957888846175e-14
};

static double _sincos(x,y,sgn,cosflag)
double x;
double y;
{
  int j,n;
  double modf(),fabs(),xn,f,g,rg;

  if(y>=1.0e8)return 0.0;	/* too large */
  if(modf(y*0.318309886183790,&xn)>=0.5)xn+=1.0;
  n=xn;
  if(n&1)sgn=!sgn;
  if(cosflag)xn-=0.5;
  f=fabs(x)-3.1416015625*xn;
  f+=8.90891020676154e-6*xn;
  if(fabs(f)>1.0e-8){
    g=f*f;
    for(rg=0.0,j=8;j--;)rg=(rg+r[j])*g;
    f=f+f*rg;
  }
  if(sgn)f=-f;
  return f;
}

double cos(arg)
double arg;
{
  double t;

  if(arg<0.0)t=-arg; else t=arg;
  return _sincos(arg,t+1.57079632679489661923,0,1);
}

double sin(arg)
double arg;
{
  int sgn;
  double t;

  if(arg<0.0){
    sgn=1;
    t=-arg;
  } else {
    sgn=0;
    t=arg;
  }
  return _sincos(arg,t,sgn,0);
}
-ARCHIVE- sqrt.c 387
/*	square root, Newton Rapson
*/

#include "errno.h"

double sqrt(arg)
double arg;
{
  int n,j;
  double f,y,frexp(),ldexp();

  if(arg<=0.0){
    if(arg<0.0)errno=EDOM;
    return 0.0;
  }
  f=frexp(arg,&n);
  y=0.41731+0.59016*f;	/* approx of root */
  for(j=0;j<4;++j)y=0.5*(y+f/y);
  if(n&1){
    y=y*0.707106781186547;
    ++n;
  }
  return ldexp(y,n>>1);
}
-ARCHIVE- tan.c 894
/*	tan and cotan
*/

#include "errno.h"

static double _tancot(x,iflag)
double x;
int iflag;
{
  int n;
  double fabs(),modf(),y,xn,f,g,xnum,xden;

  y=fabs(x);
  if(iflag && y<1.0e-300){
    if(x<0.0)return -1.0e300;
    return 1.0e300;
  }
  if(y>1.0e8)return 0.0;
  if(modf(x*0.63661977236758134308,&xn)>=0.5)xn+=1.0;
  n=xn;
  f=x-xn*1.57080078125;
  f+=xn*4.454455103380768678308e-6;
  if(fabs(f)<1.0e-8){
    xnum=f;
    xden=1.0;
  } else {
    g=f*f;
    xnum=((-0.178617073422544e-4*g+0.342488782358906e-2)*g-0.133383500064220)*g*f+f;
    xden=(((0.498194339937865e-6*g-0.311815319070100e-3)*g+0.256638322894401e-1)*g-0.466716833397553)*g+0.1e1;
  }
  if(n&1)xnum=-xnum;
  if(iflag==(n&1))return xnum/xden;
  return xden/xnum;
}

double tan(x)
double x;
{

  return _tancot(x,0);
}

double cotan(x)
double x;
{

  return _tancot(x,1);
}