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

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

atan2.c
#include "errno.h"

#define MAXEXP	1024
#define MINEXP	-1023

#define PI		3.14159265358979323846
#define PIov2	1.57079632679489661923

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

atof.c
/* Copyright (C) 1983 by Manx Software Systems */
#include	<ctype.h>

double
atof(cp)
register char *cp;
{
	double acc, zero = 0.0, ten = 10.0;
	int msign, esign, dpflg;
	int i, dexp;

	while (*cp == ' ' || *cp == '\t')
		++cp;
	if (*cp == '-') {
		++cp;
		msign = 1;
	} else {
		msign = 0;
		if (*cp == '+')
			++cp;
	}
	dpflg = dexp = 0;
	for (acc = zero ; ; ++cp) {
		if (isdigit(*cp)) {
			acc *= ten;
			acc += *cp - '0';
			if (dpflg)
				--dexp;
		} else if (*cp == '.') {
			if (dpflg)
				break;
			dpflg = 1;
		} else
			break;
	}
	if (*cp == 'e' || *cp == 'E') {
		++cp;
		if (*cp == '-') {
			++cp;
			esign = 1;
		} else {
			esign = 0;
			if (*cp == '+')
				++cp;
		}
		for ( i = 0 ; isdigit(*cp) ; i = i*10 + *cp++ - '0' )
			;
		if (esign)
			i = -i;
		dexp += i;
	}
	if (dexp < 0) {
		while (dexp++)
			acc /= ten;
	} else if (dexp > 0) {
		while (dexp--)
			acc *= ten;
	}
	if (msign)
		acc = -acc;
	return acc;
}

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(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;
}
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 };

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

	if ((*(short *)&number&0x7ff0) == 0x7ff0) {
		strcpy(buffer, "NAN");
		return;
	}
	ndig = maxwidth+1;
	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;
		}
	}

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

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

	if (flag) {
		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 (!flag) {
		*buffer++ = 'e';
		if (exp < 0) {
			exp = -exp;
			*buffer++ = '-';
		} else
			*buffer++ = '+';
		if (exp >= 100) {
			*buffer++ = exp/100 + '0';
			exp %= 100;
		}
		*buffer++ = exp/10 + '0';
		*buffer++ = exp%10 + '0';
	}
	*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
.c.o:
	cc +bf8 +x3,5 -o $@ $*.c
.c.o32:
	cc +blf8 +x3,5 -o $@ $*.c
.c.l:
	cc +bcdf8 +x3,5 -o $@ $*.c
.c.l32:
	cc +bpf8 +x3,5 -o $@ $*.c

.a68.o:
	as -o $@ $*.a68
.a68.o32:
	as -eINT32 -o $@ $*.a68
.a68.l:
	as -cdo $@ $*.a68
.a68.l32:
	as -eINT32 -cdo $@ $*.a68

C=atof.o ftoa.o floor.o random.o
ASM=abs.o atan2.o fabs.o frexp.o ldexp.o modf.o pow.o tran.o
C32=ftoa.o32 floor.o32
ASM32=frexp.o32 ldexp.o32 pow.o32
LC=atof.l ftoa.l floor.l random.l
LASM=abs.l atan2.l fabs.l frexp.l ldexp.l modf.l pow.l tran.l
LC32=ftoa.l32 floor.l32
LASM32=frexp.l32 ldexp.l32 pow.l32

all:	16 32 l16 l32

16:		$(ASM) $(C)

32:		$(ASM32) $(C32)

l16:	$(LASM) $(LC)

l32:	$(LASM32) $(LC32)

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(a,b)
double a,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 ran()
{
	static long int iy = 100001;
	
	iy *= 125;
	iy -= (iy/2796203) * 2796203;
	return (double) iy/ 2796203.0;
}

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

	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