asin.c
/* Copyright 1987 Manx Software Systems, Inc */

#include <math.h>
#include <errno.h>

static double arcsine(double x, int flg);

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

double acos(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)

static double
arcsine(double x, int flg)
{
	register double y;
	register double g;
	register long i;
	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) {
		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;

		y = y + y*
				((P(g)*g)
				/Q(g));
	}
	if (flg) {
		if (x < 0.0)
			y = (b[i] + y) + b[i];
		else
			y = (a[i] - y) + a[i];
	} else {
		y = (a[i] + y) + a[i];
		if (x < 0.0)
			y = -y;
	}
	return y;
}
atan.c
/* Copyright 1987 Manx Software Systems, Inc */

#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(register double v, double u)
{
	register double f;
	int vexp, uexp;

	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(double x)
{
	register double f, g;
	register 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) {
		g = f*f;
		f = f + f *
			((P(g)*g)
			/Q(g));
	}
	if (n > 1)
		f = -f;
	f += Avals[n];
	if (x < 0.0)
		f = -f;
	return f;
}
exp.c
/* Copyright 1987 Manx Software Systems, Inc */

#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(register double x)
{
#	define g (x)
#	define r (g)
	register double z;
#	define xn (z)
	register int n;
	
	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);
#undef g
#undef r
#undef xn
}
floor.c
#include "math.h"

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

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

double
fmod(double x, double y)
{
	double d;

	return(modf(x/y,&d)*y);
}
frexp.a68
; Copyright 1987 Manx Software Systems, Inc
;
;	FLOATING POINT ERROR VALUES
;
UNDER_FLOW	equ	1
OVER_FLOW	equ	2
DIV_BY_ZERO	equ	3
;
;	double frexp(double d, int *i)
;
;	returns 1/2 <= |x| < 1
;	such that d = x *2^i
;
		public	_frexp
_frexp:
		movem.l	4+0(sp),d0/d1/a0	; d0,d1 := d; a0 = &i;
		move.l	d0,d2
		swap	d2
		and.w	#$7ff0,d2
		beq		zero
notzero:
		and.l	#$800fffff,d0	;get rid of old exponent
		lsr.w	#4,d2
		sub.w	#1022,d2
		or.l	#$3fe00000,d0	;change exponent to -1
zero:
		IF INT32
		ext.l	d2
		move.l	d2,(a0)
		ELSE
		move.w	d2,(a0)
		ENDC
		rts
ldexp.a68
;
;	FLOATING POINT ERROR VALUES
;
UNDER_FLOW	equ	1
OVER_FLOW	equ	2
DIV_BY_ZERO	equ	3
;
ERANGE	equ	13
EDOM	equ	14
;
		xdef	_errno
;
;		ldexp(d, i)
;
;		return x = d * 2^i
;
		public	_ldexp
_ldexp	move.l	d2,-(sp)
		movem.l	8(sp),d0/d1		;get the number
		move.l	d0,d2			;get exponent
		and.l	#$800fffff,d0	;get old exponent out
		swap	d2
		and.l	#$7ff0,d2
		lsr.w	#4,d2
		if INT32
		add.w	18(sp),d2
		else
		add.w	16(sp),d2
		endc
		ble		ldunder
		cmp.l	#2047,d2
		bgt		ldover
		lsl.w	#4,d2
		swap	d2
		or.l	d2,d0			;put new exponent in
		move.l	(sp)+,d2
		rts
ldunder:
		move.w	#UNDER_FLOW,_flterr
		if INT32
		move.l	#ERANGE,_errno
		else
		move.w	#ERANGE,_errno
		endc
;		or.l	#$00100000,d0
		move.l	#$c3a99999,d0	;-HUGE_VAL
		move.l	#$8d2c7175,d1
		move.l	(sp)+,d2
		rts
ldover:
		move.w	#OVER_FLOW,_flterr
		if INT32
		move.l	#ERANGE,_errno
		else
		move.w	#ERANGE,_errno
		endc
;		or.l	#$7ff00000,d0
		move.l	#$43a99999,d0
		move.l	#$8d2c7175,d1
		move.l	(sp)+,d2
		rts
;
		global	_flterr,2
log.c
/* Copyright 1987 Manx Software Systems, Inc */

#include <math.h>
#include <errno.h>

double log10(double x)
{
	x = log(x);
	if (x != -HUGE_VAL)
		x *= 0.43429448190325182765;
	return(x);
}

#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(register double x)
{
	register double znum;
#	define z (znum)
#	define Rz (z)

#	define f (x)
#	define zden (f)
#	define w (f)
#	define xn (w)
	int n;
	
	if (x <= 0.0) {
		errno = EDOM;
		if (x == 0.0)
			errno = ERANGE;
		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;
#undef z
#undef Rz
#undef f
#undef zden
#undef w
#undef xn
}
makefile
CC=qcc
AS=as
CFLAGS=
AFLAGS=
O=oss

.c.$O:
	$(CC) -o $@ $(CFLAGS) $*.c
.a68.$O:
	$(AS) -o $@ $(AFLAGS) $*.a68

CLIB=\
	asin.$O atan.$O exp.$O floor.$O fmod.$O \
	log.$O pow.$O random.$O sin.$O sinh.$O \
	sqrt.$O tan.$O tanh.$O
ALIB=\
	frexp.$O ldexp.$O modf.$O

mlib:	$(CLIB) $(ALIB)

modf.a68
;
;	FLOATING POINT ERROR VALUES
;
UNDER_FLOW	equ	1
OVER_FLOW	equ	2
DIV_BY_ZERO	equ	3
;
;		modf(d, dptr)
;
;		return fractional part of d and
;		stores integral part in	*dptr
;
		public	.Pfix
		public	.Pflt
		public	.Psub
		public	_modf
_modf:
		movem.l	d2/d3,-(sp)
		move.l	12(sp),d0		;pick up double
		move.l	16(sp),d1
		move.l	d0,d2			;copy number
		move.l	d1,d3			;
		jsr		.Pfix			;fix the number
;		bvc		fits			;fits in long
		cmp.l	#$7fffffff,d0	;positive overflow
		beq		over1
		cmp.l	#$80000001,d0	;negative overflow
		bne		fits
over1:
		move.l	d2,d0			;overflow
		move.l	d3,d1			;so return original
		bra		over2
fits:
		jsr		.Pflt			;float it again
over2:
		move.l	20(sp),a0
		move.l	d0,(a0)+		;store integral part
		move.l	d1,(a0)
		exg.l	d0,d2
		exg.l	d1,d3
		jsr		.Psub			;get rid of integral part
		movem.l	(sp)+,d2/d3
		rts
oldexp.a68
; Copyright 1987 Manx Software Systems, Inc
;
;	FLOATING POINT ERROR VALUES
;
ERANGE		equ	13
EDOM		equ	14
;
;		double ldexp(double d, int i)
;
;		return x = d * 2^i
;
		public	_ldexp
_ldexp:
		movem.l	4+0(sp),d0/d1	;get the number
		swap	d0
		move.w	d0,d2			;get exponent
		and.w	#$800f,d0	;get old exponent out
		and.w	#$7ff0,d2
		lsr.w	#4,d2
		IF INT32
		add.w	4+8+2(sp),d2
		ELSE
		add.w	4+8(sp),d2
		ENDC
;		bmi		ldunder
;		beq		ldunder
		ble		ldunder
		cmp.w	#2047,d2
		bgt		ldover
		lsl.w	#4,d2
		or.w	d2,d0			;put new exponent in
		swap	d0
		rts
ldunder:
		IF INT32
		move.l	#ERANGE,_errno	;under flow
		ELSE
		move.w	#ERANGE,_errno	;under flow
		ENDC
		swap	d0
		or.l	#$00100000,d0
		rts
ldover:
		IF INT32
		move.l	#ERANGE,_errno	;over flow
		ELSE
		move.w	#ERANGE,_errno	;over flow
		ENDC
		swap	d0
		or.l	#$7ff00000,d0
		rts
;
		IF INT32
		global	_errno,4
		ELSE
		global	_errno,2
		ENDC
pow.c
#include "math.h"
#include "errno.h"

double pow(double a, double b)
{
	double ans;
	register 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,&ans) == 0) {
		if ((answer = ans) < 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.
 */

double exp(double), ran(void);

static long int iy = 100001;

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

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

double
randl(double x)
{

	return exp(x*ran());
}

sin.c
/* Copyright 1987 Manx Software Systems, Inc */

#include <math.h>
#include <errno.h>

static double sincos(double x, register double y, int sgn);

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

double sin(register double x)
{
	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(double x, register double y, int sgn)
{
	double xn;
	register double f;
#	define g (y)

	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;
#undef g
}
sinh.c
/* Copyright 1987 Manx Software Systems, Inc */

#include <math.h>
#include <errno.h>

#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(register double x)
{
#	define w (x)
	register double y;
#	define z (y)
	char 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;
#undef w
#undef z
}

double cosh(double x)
{
	register double y;
#	define w (y)
	register double 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;
#undef w
}
sqrt.c
/* Copyright 1987 Manx Software Systems, Inc */

#include <math.h>
#include <errno.h>

double sqrt(register double x)
{
	register double y;
	int n;
	
	if (x <= 0.0) {
		if (x == 0.0) return x;
		errno = EDOM;
		return 0.0;
	}
    y = 0.4173075996388649989089 + 0.59016206709064458299663 * frexp(x, &n);
	if (n&1)
		y *= 1.41421356237309504880;	/* sqrt(2) */
	y = ldexp(y, n >> 1);

	y = ldexp(y + x/y, -1);				/* y = (y + x/y) / 2.0 */
	y = ldexp(y + x/y, -1);				/* y = (y + x/y) / 2.0 */
 return ldexp(y + x/y, -1);				/* y = (y + x/y) / 2.0 */
}
tan.c
/* Copyright 1987 Manx Software Systems, Inc */

#include <math.h>
#include <errno.h>

static double tansub(double x, int flag);

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

#define YMAX 6.74652e09

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

double tan(double x)
{
	if (fabs(x) > YMAX) {
		errno = ERANGE;
		return 0.0;
	}

	return tansub(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)


static double tansub(register double x, int flag)
{
#	define f (x)
#	define xnum (f)
	register double g;
#	define xden (g)
	double xn, ag;
	
	if (fabs(modf(x*0.63661977236758134308, &xn)) >= 0.5)
		xn += (x < 0.0) ? -1.0 : 1.0;
	f = modf(x, &ag);
	f = ((ag - 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;
#undef f
#undef xnum
#undef xden
}
tanh.c
/* Copyright 1987 Manx Software Systems, Inc */

#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(double x)
{
	register double r;
#	define f (r)
	register double g;
	
	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;
#undef f
}