asin.c
#include "math.h"
#include "errno.h"

double arcsine();

double asin(x)
double x;
{
	return arcsine(x,0);
}

double acos(x)
double x;
{
	return arcsine(x,1);
}

#define P1 -0.27368494524164255994e+2
#define P2 +0.57208227877891731407e+2
#define P3 -0.39688862997504877339e+2
#define P4 +0.10152522233806463645e+2
#define P5 -0.69674573447350646411
#define Q0 -0.16421096714498560795e+3
#define Q1 +0.41714430248260412556e+3
#define Q2 -0.38186303361750149284e+3
#define Q3 +0.15095270841030604719e+3
#define Q4 -0.23823859153670238830e+2

#define P(g) ((((P5*g P4)*g P3)*g P2)*g P1)
#define Q(g) (((((g Q4)*g Q3)*g Q2)*g Q1)*g Q0)

double arcsine(x,flg)
double x;
{
	double y, g, r;
	register int i;
	extern int errno;
	static double a[2] = { 0.0, 0.78539816339744830962 };
	static double b[2] = { 1.57079632679489661923, 0.78539816339744830962 };

	y = fabs(x);
	i = flg;
	if (y < 2.3e-10)
		r = y;
	else {
		if (y > 0.5) {
			i = 1-i;
			if (y > 1.0) {
				errno = EDOM;
				return 0.0;
			}
			g = (0.5-y)+0.5;
			g = ldexp(g,-1);
			y = sqrt(g);
			y = -(y+y);
		} else
			g = y*y;
		r = y + y*
				((P(g)*g)
				/Q(g));
	}
	if (flg) {
		if (x < 0.0)
			r = (b[i] + r) + b[i];
		else
			r = (a[i] - r) + a[i];
	} else {
		r = (a[i] + r) + a[i];
		if (x < 0.0)
			r = -r;
	}
	return r;
}
atan.c
#include "math.h"
#include "errno.h"

#ifdef MPU8086
#define MAXEXP	1024
#define MINEXP	-1023
#else
#define MAXEXP	504
#define MINEXP	-512
#endif

#define PI		3.14159265358979323846
#define PIov2	1.57079632679489661923

double atan2(v,u)
double u,v;
{
	double f, frexp();
	int vexp, uexp;
	extern int flterr;
	extern int errno;

	if (u == 0.0) {
		if (v == 0.0) {
			errno = EDOM;
			return 0.0;
		} else if (v > 0.0 )
			return PIov2;
		return -PIov2;
	}

	frexp(v, &vexp);
	frexp(u, &uexp);
	if (vexp-uexp > MAXEXP-3)	/* overflow */
		f = PIov2;
	else {
		if (vexp-uexp < MINEXP+3)	/* underflow */
			f = 0.0;
		else
			f = atan(fabs(v/u));
		if (u < 0.0)
			f = PI - f;
	}
	if (v < 0.0)
		f = -f;
	return f;
}

#define P0 -0.13688768894191926929e+2
#define P1 -0.20505855195861651981e+2
#define P2 -0.84946240351320683534e+1
#define P3 -0.83758299368150059274e+0
#define Q0 +0.41066306682575781263e+2
#define Q1 +0.86157349597130242515e+2
#define Q2 +0.59578436142597344465e+2
#define Q3 +0.15024001160028576121e+2

#define P(g) (((P3*g P2)*g P1)*g P0)
#define Q(g) ((((g Q3)*g Q2)*g Q1)*g Q0)

double atan(x)
double x;
{
	double f, r, g;
	int n;
	static double Avals[4] = {
		0.0,
		0.52359877559829887308,
		1.57079632679489661923,
		1.04719755119659774615
	};
	
	n = 0;
	f = fabs(x);
	if (f > 1.0) {
		f = 1.0/f;
		n = 2;
	}
	if (f > 0.26794919243112270647) {
		f = (((0.73205080756887729353*f - 0.5) - 0.5) + f) /
				(1.73205080756887729353 + f);
		++n;
	}
	if (fabs(f) < 2.3e-10)
		r = f;
	else {
		g = f*f;
		r = f + f *
			((P(g)*g)
			/Q(g));
	}
	if (n > 1)
		r = -r;
	r += Avals[n];
	if (x < 0.0)
		r = -r;
	return r;
}
exp.c
#include "math.h"
#include "errno.h"

#define P0 0.25000000000000000000e+0
#define P1 0.75753180159422776666e-2
#define P2 0.31555192765684646356e-4
#define Q0 0.50000000000000000000e+0
#define Q1 0.56817302698551221787e-1
#define Q2 0.63121894374398503557e-3
#define Q3 0.75104028399870046114e-6

#define P(z) ((P2*z + P1)*z + P0)
#define Q(z) (((Q3*z + Q2)*z + Q1)*z + Q0)

#define EPS	2.710505e-20

double
exp(x)
double x;
{
	int n;
	double xn, g, r, z;
	extern int errno;
	
	if (x > LOGHUGE) {
		errno = ERANGE;
		return HUGE_VAL;
	}
	if (x < LOGTINY) {
		errno = ERANGE;
		return 0.0;
	}
	if (fabs(x) < EPS)
		return 1.0;
	n = z = x * 1.4426950408889634074;
	if (n < 0)
		--n;
	if (z-n >= 0.5)
		++n;
	xn = n;
	g = ((x - xn*0.693359375)) + xn*2.1219444005469058277e-4;
	z = g*g;
	r = P(z)*g;
	r = 0.5 + r/(Q(z)-r);
	return ldexp(r,n+1);
}
floor.c
#include "math.h"

double floor(d)
double d;
{
	if (d < 0.0)
		return -ceil(-d);
	modf(d, &d);
	return d;
}

double ceil(d)
double d;
{
	if (d < 0.0)
		return -floor(-d);
	if (modf(d, &d) > 0.0)
		++d;
	return d;
}
log.c
#include "math.h"
#include "errno.h"

double log10(x)
double x;
{
	return log(x)*0.43429448190325182765;
}

#define A0 -0.64124943423745581147e+2
#define A1 +0.16383943563021534222e+2
#define A2 -0.78956112887491257267e+0
#define A(w) ((A2*w A1)*w A0)

#define B0 -0.76949932108494879777e+3
#define B1 +0.31203222091924532844e+3
#define B2 -0.35667977739034646171e+2
#define B(w) (((w B2)*w B1)*w B0)

#define C0 0.70710678118654752440
#define C1 0.693359375
#define C2 -2.121944400546905827679e-4

double log(x)
double x;
{
	double Rz, f, z, w, znum, zden, xn;
	int n;
	extern int errno;
	
	if (x <= 0.0) {
		errno = EDOM;
		return -HUGE_VAL;
	}
	f = frexp(x, &n);
	if (f > C0) {
		znum = (znum = f-0.5) - 0.5; /* the assignment prevents const. eval */
		zden = f*0.5 + 0.5;
	} else {
		--n;
		znum = f - 0.5;
		zden = znum*0.5 + 0.5;
	}
	z = znum/zden;
	w = z*z;
/* the lines below are split up to allow expansion of A(w) and B(w) */
	Rz = z + z * (w *
			 A(w)
			/B(w));
	xn = n;
	return (xn*C2 + Rz) + xn*C1;
}
makefile
.c.o:
	cc +bfi +x3,5 -o $@ $*.c
.c.o32:
	cc +blfi +x3,5 -o $@ $*.c
.c.l:
	cc +bcdfi +x3,5 -o $@ $*.c
.c.l32:
	cc +bpfi +x3,5 -o $@ $*.c

C=asin.o atan.o exp.o floor.o log.o\
	pow.o random.o sin.o sinh.o sqrt.o\
	tan.o tanh.o
C32=asin.o32 atan.o32 exp.o32 floor.o32 log.o32\
	pow.o32 random.o32 sin.o32 sinh.o32 sqrt.o32\
	tan.o32 tanh.o32
LC=asin.l atan.l exp.l floor.l log.l\
	pow.l random.l sin.l sinh.l sqrt.l\
	tan.l tanh.l
LC32=asin.l32 atan.l32 exp.l32 floor.l32 log.l32\
	pow.l32 random.l32 sin.l32 sinh.l32 sqrt.l32\
	tan.l32 tanh.l32

all:	16 32 l16 l32

16:		$(C)

32:		$(C32)

l16:	$(LC)

l32:	$(LC32)

npow.c
#include "math.h"
#include "errno.h"

double pow(a,b)
double a,b;
{
	double loga;
	extern int errno;
	double whole, frac;
	register int i;
	int sgn = 0;

	frac = modf(b,&whole);
	if (!frac) {
		if (whole < 0) {
			sgn = 1;
			whole = -whole;
		}
		i = whole;
		while (--i)
			a = a * a;
		if (sgn)
			a = 1.0/a;
		return a;
	}
	if (a<=0.0) {
		if (a<0.0 || a==0.0 && b<=0.0) {
			errno = EDOM;
			return -HUGE_VAL;
		}
		else return 0.0;
	}
	loga = log(a);
	loga *= b;
	if (loga > LOGHUGE) {
		errno = ERANGE;
		return HUGE_VAL;
	}
	if (loga < LOGTINY) {
		errno = ERANGE;
		return 0.0;
	}
	return exp(loga);
}
pow.c
#include "math.h"
#include "errno.h"

double pow(a,b)
double a,b;
{
	double answer;
	extern int errno;
	register long count;
	char sign, inverse;
	
	if (a == 0) {
		if (b <= 0)
domain:		errno = EDOM;
		return 0.0;
	}
	if (b == 0)
		return 1.0;		/* anything raised to 0 is 1 */
	inverse = sign = 0;
	if (modf(b,&answer) == 0) {
		if (answer < 0)
			inverse = 1, answer = -answer;
		if ((count = answer) == answer) {
			for (answer = 1.0 ; count ; count >>= 1, a *= a)
				if ((int)count & 1)
					answer *= a;
			if (inverse)
				answer = 1.0/answer;
			return answer;
		}
		if (a < 0)
			sign = 1, a = -a;
		if ((count&1) == 0)
			sign = 0;		/* number is even so sign is positive */

	} else if (a < 0)
		goto domain;

	answer = exp(log(a)*b);
	return sign ? -answer : answer;
}
random.c
/*
 * Random number generator -
 * adapted from the FORTRAN version 
 * in "Software Manual for the Elementary Functions"
 * by W.J. Cody, Jr and William Waite.
 */

static long int iy = 100001;

sran(seed)
long seed;
{
	iy = seed;
}

double ran()
{
	iy *= 125;
	iy -= (iy/2796203) * 2796203;
	return (double) iy/ 2796203.0;
}

double randl(x)
double x;
{
	double exp();

	return exp(x*ran());
}

sin.c
#include "math.h"
#include "errno.h"

double cos(x)
double x;
{
	double sincos();

	return sincos(x, fabs(x) + 1.57079632679489661923, 0);
}

double sin(x)
double x;
{
	double sincos();
	
	if (x < 0.0)
		return sincos(x,-x,1);
	else
		return sincos(x,x,0);
}

#define R1 -0.16666666666666665052e+00
#define R2 +0.83333333333331650314e-02
#define R3 -0.19841269841201840457e-03
#define R4 +0.27557319210152756119e-05
#define R5 -0.25052106798274584544e-07
#define R6 +0.16058936490371589114e-09
#define R7 -0.76429178068910467734e-12
#define R8 +0.27204790957888846175e-14

#define YMAX 6.7465e09

static double sincos(x,y,sgn)
double x,y;
{
	double f, xn, g;
	extern int errno;

	if (y >= YMAX) {
		errno = ERANGE;
		return 0.0;
	}
	if (modf(y * 0.31830988618379067154, &xn) >= 0.5)
		++xn;
	if ((int)xn & 1)
		sgn = !sgn;
	if (fabs(x) != y)
		xn -= 0.5;
	g = modf(fabs(x), &x);		/* break into fraction and integer parts */
	f = ((x - xn*(3217.0/1024)) + g) - xn*-8.9089102067615373566e-6;
	if (fabs(f) > 2.3283e-10) {
		g = f*f;
		f = (((((((R8*g R7)*g R6)*g R5)*g
				R4)*g R3)*g R2)*g R1)*g*f+f;
	}
	if (sgn)
		f = -f;
	return f;
}
sinh.c
#include "math.h"
#include "errno.h"

extern int errno;

#define P0 -0.35181283430177117881e+6
#define P1 -0.11563521196851768270e+5
#define P2 -0.16375798202630751372e+3
#define P3 -0.78966127417357099479e+0
#define Q0 -0.21108770058106271242e+7
#define Q1 +0.36162723109421836460e+5
#define Q2 -0.27773523119650701667e+3

#define PS(x) (((P3*x P2)*x P1)*x P0)
#define QS(x) (((x Q2)*x Q1)*x Q0)

double sinh(x)
double x;
{
	double y, w, z;
	int sign;
	
	y = x;
	sign = 0;
	if (x < 0.0) {
		y = -x;
		sign = 1;
	}
	if (y > 1.0) {
		w = y - 0.6931610107421875000;
		if (w > LOGHUGE) {
			errno = ERANGE;
			z = HUGE_VAL;
		} else {
			z = exp(w);
			if (w < 19.95)
				z -= 0.24999308500451499336 / z;
			z += 0.13830277879601902638e-4 * z;
		}
		if (sign)
			z = -z;
	} else if (y < 2.3e-10)
		z = x;
	else {
		z = x*x;
		z = x + x *
				(z*(PS(z)
				/QS(z)));
	}
	return z;
}

double cosh(x)
double x;
{
	double y, w, z;
	
	y = fabs(x);
	if (y > 1.0) {
		w = y - 0.6931610107421875000;
		if (w > LOGHUGE) {
			errno = ERANGE;
			return HUGE_VAL;
		}
		z = exp(w);
		if (w < 19.95)
			z += 0.24999308500451499336 / z;
		z += 0.13830277879601902638e-4 * z;
	} else {
		z = exp(y);
		z = z*0.5 + 0.5/z;
	}
	return z;
}
sqrt.c
#include "math.h"
#include "errno.h"

double sqrt(x)
double x;
{
	double f, y;
	int n;
	extern int errno;
	
	if (x == 0.0)
		return x;
	if (x < 0.0) {
		errno = EDOM;
		return 0.0;
	}
	f = frexp(x, &n);
	y = 0.41731 + 0.59016 * f;
	y = (y + f/y);
	y = ldexp(y,-2) + f/y;	/* fast calculation of y2 */
	y = ldexp(y + f/y, -1);
	y = ldexp(y + f/y, -1);
	
	if (n&1) {
		y *= 0.70710678118654752440;
		++n;
	}
	return ldexp(y,n/2);
}
tan.c
#include "math.h"
#include "errno.h"

extern int errno;

static double tansub();

#if MPU8080 || MPUZ80 || MPU6502
#define TOOSMALL	(1.0/HUGE_VAL)
#else
#define TOOSMALL	TINY_VAL
#endif

double cotan(x)
double x;
{
	double y;
	
	y = fabs(x);
	if (y < TOOSMALL) {
		errno = ERANGE;
		if (x < 0.0)
			return -HUGE_VAL;
		else
			return HUGE_VAL;
	}
	return tansub(x,y,2);
}

double tan(x)
double x;
{
	return tansub(x, fabs(x), 0);
}

#define P1 -0.13338350006421960681e+0
#define P2 +0.34248878235890589960e-2
#define P3 -0.17861707342254426711e-4
#define Q0 +1.0
#define Q1 -0.46671683339755294240e+0
#define Q2 +0.25663832289440112864e-1
#define Q3 -0.31181531907010027307e-3
#define Q4 +0.49819433993786512270e-6

#define P(f,g) (((P3*g P2)*g P1)*g*f + f)
#define Q(g) ((((Q4*g Q3)*g Q2)*g Q1)*g Q0)

#define YMAX 6.74652e09

static double tansub(x, y, flag)
double x,y;
{
	double f, g, xn;
	double xnum, xden;
	
	if (y > YMAX) {
		errno = ERANGE;
		return 0.0;
	}
	if (fabs(modf(x*0.63661977236758134308, &xn)) >= 0.5)
		xn += (x < 0.0) ? -1.0 : 1.0;
	f = modf(x, &g);
	f = ((g - xn*(3217.0/2048)) + f) - xn*-4.454455103380768678308e-6;
	if (fabs(f) < 2.33e-10) {
		xnum = f;
		xden = 1.0;
	} else {
		g = f*f;
		xnum = P(f,g);
		xden = Q(g);
	}
	flag |= ((int)xn & 1);
	switch (flag) {
	case 1:		/* A: tan, xn odd */
		xnum = -xnum;
	case 2:		/* B: cotan, xn even */
		return xden/xnum;
		
	case 3:		/* C: cotan, xn odd */
		xnum = -xnum;
	case 0:		/* D: tan, xn even */
		return xnum/xden;
	}
	return 0.0;
}
tanh.c
#include "math.h"

#define P0 -0.16134119023996228053e+4
#define P1 -0.99225929672236083313e+2
#define P2 -0.96437492777225469787e+0
#define Q0 +0.48402357071988688686e+4
#define Q1 +0.22337720718962312926e+4
#define Q2 +0.11274474380534949335e+3

#define gP(g) (((P2*g P1)*g P0)*g)
#define Q(g) (((g Q2)*g Q1)*g Q0)

double tanh(x)
double x;
{
	double f,g,r;
	
	f = fabs(x);
	if (f > 25.3)
		r = 1.0;
	else if (f > 0.54930614433405484570) {
		r = 0.5 - 1.0/(exp(f+f)+1.0);
		r += r;
	} else if (f < 2.3e-10)
		r = f;
	else {
		g = f*f;
		r = f + f*
			(gP(g)
			/Q(g));
	}
	if (x < 0.0)
		r = -r;
	return r;
}