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

#define MAXEXP	1024
#define MINEXP	-1023

#define PI		3.14159265358979323846
#define PIov2	1.57079632679489661923

double atan2(double v, double u)
{
	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;
}

fabs.a68
; Copyright (C) 1985 by Manx Software Systems, Inc.
; :ts=8

	mc68881
	public	_fabs
_fabs
	fabs.d	4(sp),fp0
	rts

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.a68
; Copyright (C) 1988 by Manx Software Systems, Inc.
; :ts=8

	mc68881
	public	_fmod
_fmod
	fmove.d	4(sp),fp0	;pick up double
	fmod.d	12(sp),fp0
	rts
frexp.a68
;
;	FLOATING POINT ERROR VALUES
;
UNDER_FLOW	equ	1
OVER_FLOW	equ	2
DIV_BY_ZERO	equ	3
;
;	frexp(d, &i)
;
;	returns 1/2 <= |x| < 1
;	such that d = x *2^i
;
		mc68881
		public	_frexp
_frexp:
		move.l	12(sp),a0		;get address for int
		move.l	4(sp),d0		; get double value
		move.l	d0,d1
		swap	d1
		and.w	#$7ff0,d1
		bne		notzero
		move.l	#0,d1
		beq		done
notzero:
		and.l	#$800fffff,d0	;get rid of old exponent
		lsr.w	#4,d1
		sub.w	#1022,d1
		or.l	#$3fe00000,d0	;change exponent to -1
done:
		if INT32
		ext.l	d1
		move.l	d1,(a0)
		else
		move.w	d1,(a0)
		endc
		move.l	d0,4(sp)		;put back d0
		fmove.d	4(sp),fp0		;get into fp0
		rts
ftoa.c
/* Copyright (C) 1984 by Manx Software Systems, Inc. */

static double round[] = { 10, 1,
	5e-1, 5e-2, 5e-3, 5e-4, 5e-5, 5e-6, 5e-7, 5e-8, 5e-9, 5e-10,
	5e-11, 5e-12, 5e-13, 5e-14, 5e-15, 5e-16 };

void
ftoa(double number, register char *buffer, int maxwidth, int mode)
{
	register int i;
	register int exp, digit, decpos, ndig;
	register int flag;

	if ((*(short *)&number&0x7ff0) == 0x7ff0) {
		exp = (*(short *)&number&0x8000) ? '-' : '+';
		while (maxwidth--)
			*buffer++ = exp;
		*buffer = 0;
		return;
	}

	/*
	 *	Convert the number to a number between 1 and 10 or 0. Calculate the
	 *	exponent at the same time.
	 */

	exp = 0;
	if (number < 0.0) {
		number = -number;
		*buffer++ = '-';
	}
	if (number > 0.0) {
		while (number < round[1]) {
			number *= round[0];
			--exp;
		}
		while (number >= round[0]) {
			number /= round[0];
			++exp;
		}
	}


	/*
	 *	The number of digits to display is equal to the number of digits
	 *	following the decimal point (maxwidth) + digits to the left (exp) +
	 *	a guaranteed digit since the number 10 > x >= 1. For E format, the
	 *	exp digits are not counted, while G format works differently.
	 */

	flag = mode & 3;
	if (flag == 2) {						/* 'g' format */
		if (maxwidth == 0)
			maxwidth = 1;
		if (exp < -4 || exp >= maxwidth)
			flag = -1;						/* switch to 'e' format */
		ndig = maxwidth;
	} else if (flag == 1)					/* 'f' format */
		ndig = maxwidth + exp + 1;
	else									/* 'e' format */
		ndig = maxwidth + 1;

	/*
	 *	Round the number up using 5e-x where x is one more than the
	 *	precision of the number.
	 */

	if (ndig > 0) {
		if ((number += round[(ndig>16?16:ndig)+1]) >= round[0]) {
			number = round[1];
			++exp;
			if (flag > 0)
				++ndig;
		}
	}

	if (flag > 0) {
		if (exp < 0) {
			*buffer++ = '0';
			*buffer++ = '.';
			i = -exp - 1;
			if (ndig <= 0)
				i = maxwidth;
			while (i--)
				*buffer++ = '0';
			decpos = 0;
		} else {
			decpos = exp+1;
		}
	} else {
		decpos = 1;
	}

	if (ndig > 0) {
		for (i = 0 ; ; ++i) {
			if (i < 16) {
				digit = (int)number;
				*buffer++ = digit+'0';
				number = (number - digit) * round[0];
			} else
				*buffer++ = '0';
			if (--ndig == 0)
				break;
			if (decpos && --decpos == 0)
				*buffer++ = '.';
		}
	}

	if (decpos && (mode&0x20))
		*buffer++ = '.';
	if (flag <= 0) {
		*buffer++ = (mode&0x10)?'E':'e';
		if (exp < 0) {
			exp = -exp;
			*buffer++ = '-';
		} else
			*buffer++ = '+';
		*buffer++ = exp/100 + '0';
		exp %= 100;
		*buffer++ = exp/10 + '0';
		*buffer++ = exp%10 + '0';
	}
	if (decpos == 0 && (flag == 2 || flag == -1) && (mode&0x20) == 0) {
		while (*--buffer == '0')
			;
		if (*buffer != '.')
			++buffer;
	}
	*buffer = 0;
}

ldexp.a68
;:ts=8
;
;		ldexp(d, i)
;
;		return x = d * 2^i
;
		mc68881
		public	_ldexp
_ldexp:
	fmove.d	4(sp),fp0		;get the number
	if INT32
	fscale.l  12(sp),fp0		;compute the result
	else
	fscale.w  12(sp),fp0		;compute the result
	endc
	rts

makefile
CC=qcc
AS=as
CFLAGS=
AFLAGS=
O=oss

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

CLIB=\
	atan2.$O floor.$O pow.$O random.$O
ALIB=\
	fabs.$O fmod.$O frexp.$O ldexp.$O modf.$O \
	tran.$O

mlib:	$(CLIB) $(ALIB)

modf.a68
;:ts=8
;
;	modf(d, dptr)
;
;	return fractional part of d and
;	stores integral part in	*dptr
;
	mc68881
	public	_modf
_modf:
	move.l	12(sp),a0	;pick up pointer for storing
	fmove.d	4(sp),fp0	;pick up double
	fmove	fp0,fp1		;copy it
	fintrz	fp1		;throw away fractional part
	fmove.d	fp1,(a0)	;store the integral part
	fsub	fp1,fp0		;throw away the integral part
	rts			;and return the fractional part

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

double pow(double a, double b)
{
	double answer, tmp;
	extern int errno;
	char sign;
	
	if (b == 0)
		return 1.0;		/* anything raised to 0 is 1 */
	if (a == 0) {
		if (b <= 0)
domain:		errno = EDOM;
		return 0.0;
	}
	sign = 0;
	if (modf(b,&answer) == 0) {
		if (a < 0)
			sign = 1, a = -a;
		if (answer < 0)
			answer = -answer;
		modf(answer/2, &tmp);
		if (tmp*2 == answer)
			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());
}

tran.a68
; Copyright (C) 1985 by Manx Software Systems, Inc.
; :ts=8

	mc68881
	public	_acos
_acos
	facos.d	4(sp),fp0
	rts

	public	_asin
_asin
	fasin.d	4(sp),fp0
	rts

	public	_atan
_atan
	fatan.d	4(sp),fp0
	rts

	public	_atanh
_atanh
	fatanh.d	4(sp),fp0
	rts

	public	_cos
_cos
	fcos.d	4(sp),fp0
	rts

	public	_cosh
_cosh
	fcosh.d	4(sp),fp0
	rts

	public	_cotan
_cotan
	fmove.w	#1,fp0
	ftan.d	4(sp),fp1
	fdiv.x	fp1,fp0
	rts

	public	_exp
_exp
	fetox.d	4(sp),fp0
	rts

	public	_log10
_log10
	flog10.d	4(sp),fp0
	rts

	public	_log2
_log2
	flog2.d	4(sp),fp0
	rts

	public	_log
_log
	flogn.d	4(sp),fp0
	rts

	public	_sin
_sin
	fsin.d	4(sp),fp0
	rts

	public	_sinh
_sinh
	fsinh.d	4(sp),fp0
	rts

	public	_sqrt
_sqrt
	fsqrt.d	4(sp),fp0
	rts

	public	_tan
_tan
	ftan.d	4(sp),fp0
	rts

	public	_tanh
_tanh
	ftanh.d	4(sp),fp0
	rts