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

	public	.Pabs
	public	_abs
_abs
	movem.l	4(sp),d0/d1
	jmp	.Pabs

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

	public	_ceil
_ceil:
	movem.l	4(sp),d0/d1
	move.l	#_LVOIEEEDPCeil#,-(sp)
	jmp	amiga_ieee#

dtof.a68
;
;	FLOATING POINT ERROR VALUES
;
ERANGE		equ	13
EDOM		equ	14
;
;	convert double to float
;
;	argument in d0/d1	result in d0
;
		public	.dtof
.dtof:	movem.l	d2/d3/d4,-(sp)
		move.l	d0,d2			;save a copy for exponent extraction
		beq		easy_exit		;no work if zero
		smi		d4				;get sign of result
		swap	d2
		lsr.w	#4,d2
		and.w	#$7ff,d2		;extract exponent
		sub.w	#1023,d2		;unbias
		and.l	#$fffff,d0		;extract mantissa
		or.l	#$100000,d0		;turn hidden bit back on
		move.l	#2,d3
dagain:
		lsl.l	#1,d1			;shift mantissa into place
		roxl.l	#1,d0
		dbra	d3,dagain
		cmp.l	#$80000000,d1	;check for round
		beq		rstar
		bcs		nornd
		add.l	#1,d0			;round up
		cmp.l	#$00ffffff,d0	;check for carry
		bls		nornd
		add.w	#1,d2			;bump exponent
		lsr.l	#1,d0			;and shift result back into place
		bra		nornd
rstar:
		or.w	#1,d0			;rstar case, force low bit to 1
nornd:
		add.w	#127,d2			;add in float bias
		bmi		underflow		;can't represent double in float
		cmp.w	#255,d2
		bgt		overflow		;ditto
		lsl.w	#7,d2			;back in place
		tst.b	d4				;need to set sign?
		beq		ndsign			;no
		or.w	#$8000,d2		;yes, set it
ndsign:
		swap	d2				;get exponent in high word
		clr.w	d2				;clear out
		and.l	#$7fffff,d0		;get rid of hidden bit
		or.l	d2,d0			; put back in number
easy_exit:
		movem.l	(sp)+,d2/d3/d4
		rts
underflow:
		IF INT32
		move.l	#ERANGE,_errno	;under flow
		ELSE
		move.w	#ERANGE,_errno	;under flow
		ENDC
		move.l	#$00800000,d0	;set to smallest number
		bra		dosign
overflow:
		IF INT32
		move.l	#ERANGE,_errno	;over flow
		ELSE
		move.w	#ERANGE,_errno	;over flow
		ENDC
		move.l	#$7fffffff,d0	;set to largest number
dosign:
		tst.b	d4
		beq		dexit
		or.l	#$80000000,d0	;set sign
dexit:
		movem.l	(sp)+,d2/d3/d4
		rts
;
		IF INT32
		global	_errno,4
		ELSE
		global	_errno,2
		ENDC
fabs.a68
; Copyright (C) 1985 by Manx Software Systems, Inc.
; :ts=8

	public	.Pabs
	public	_fabs
_fabs
	movem.l	4(sp),d0/d1
	jmp	.Pabs

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

	public	.Ffix
.Ffix:
	move.l	d1,-(sp)
	jsr	.ftod#
	jsr	.Pfix#
	movem.l	(sp)+,d1
	rts

	public	.Fflt
.Fflt:
	move.l	d1,-(sp)
	jsr	.Pflt#
	jsr	.dtof#
	movem.l	(sp)+,d1
	rts

	public	.Fneg
.Fneg:
	move.l	d1,-(sp)
	jsr	.ftod#
	jsr	.Pneg#
	jsr	.dtof#
	movem.l	(sp)+,d1
	rts

	public	.Ftst
.Ftst:
	movem.l	d0/d1,-(sp)
	jsr	.ftod#
	jsr	.Ptst#
	movem.l	(sp)+,d0/d1
	rts

	public	.Fadd
.Fadd:
	movem.l	d2/d3,-(sp)
	move.l	d1,d2
	jsr	.ftod#
	exg	d0,d2
	move.l	d1,d3
	jsr	.ftod#
	jsr	.Padd#
	jsr	.dtof#
	movem.l	(sp)+,d2/d3
	rts

	public	.Fcmp
.Fcmp:
	movem.l	d2/d3,-(sp)
	move.l	d0,d2
	move.l	d1,d0
	jsr	.ftod#
	exg	d0,d2
	move.l	d1,d3
	jsr	.ftod#
	jsr	.Pcmp#
	movem.l	(sp)+,d2/d3
	rts

	public	.Fdiv
.Fdiv:
	movem.l	d2/d3,-(sp)
	move.l	d0,d2
	move.l	d1,d0
	jsr	.ftod#
	exg	d0,d2
	move.l	d1,d3
	jsr	.ftod#
	jsr	.Pdiv#
	jsr	.dtof#
	movem.l	(sp)+,d2/d3
	rts

	public	.Fmul
.Fmul:
	movem.l	d2/d3,-(sp)
	move.l	d0,d2
	move.l	d1,d0
	jsr	.ftod#
	exg	d0,d2
	move.l	d1,d3
	jsr	.ftod#
	jsr	.Pmul#
	jsr	.dtof#
	movem.l	(sp)+,d2/d3
	rts

	public	.Fsub
.Fsub:
	movem.l	d2/d3,-(sp)
	move.l	d0,d2
	move.l	d1,d0
	jsr	.ftod#
	exg	d0,d2
	move.l	d1,d3
	jsr	.ftod#
	jsr	.Psub#
	jsr	.dtof#
	movem.l	(sp)+,d2/d3
	rts

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

	public	_floor
_floor:
	movem.l	4(sp),d0/d1
	move.l	#_LVOIEEEDPFloor#,-(sp)
	jmp	amiga_ieee#

fmod.c
#include <math.h>

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

	return(modf(x/y,&d)*y);
}
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
;
		public	_frexp
_frexp:	move.l	d2,-(sp)
		move.l	16(sp),a0		;get address for int
		move.l	8(sp),d0		; get double value
		move.l	12(sp),d1
		move.l	d0,d2
		swap	d2
		and.w	#$7ff0,d2
		bne		notzero
		move.l	#0,d2
		beq		done
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
done:
		if INT32
		ext.l	d2
		move.l	d2,(a0)
		else
		move.w	d2,(a0)
		endc
		move.l	(sp)+,d2
		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;
}

ftod.a68
;
;	convert float to double
;
;	argument in d0 result in d0/d1
;
		public	.ftod
.ftod:	movem.l	d3/d4,-(sp)
		tst.l	d0				;save sign
		smi		d3
		clr.l	d1				;clear out LL of result
		and.l	#$7fffffff,d0	;extract exponent and mantissa
		move.l	#2,d4
fagain:
		lsr.l	#1,d0			;shift exponent and mantissa into position
		roxr.l	#1,d1
		dbra	d4,fagain
		move.l	d0,d4
		or.l	d1,d4			;check for 0
		beq		1$
		add.l	#$38000000,d0	;adjust bias for double (1023-127) << 20
1$:
		tst.b	d3				; have to set sign?
		beq		fexit			; no
		or.l	#$80000000,d0	;yes, set it
fexit:
		movem.l	(sp)+,d3/d4
		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
ldtoa.c
void ftoa(double, char *, int, int);

void 
_ldtoa(long double number, char *buffer, int maxwidth, int flag)
{
	ftoa((double)number, buffer, maxwidth, flag);
}
lvo.a68
; Copyright (C) 1986 by Manx Software Systems, Inc.
; :ts=8

_LVOIEEEDPCeil	equ	-96
_LVOIEEEDPFloor	equ	-90
_LVOIEEEDPDiv	equ	-84
_LVOIEEEDPMul	equ	-78
_LVOIEEEDPSub	equ	-72
_LVOIEEEDPAdd	equ	-66
_LVOIEEEDPNeg	equ	-60
_LVOIEEEDPAbs	equ	-54
_LVOIEEEDPTst	equ	-48
_LVOIEEEDPCmp	equ	-42
_LVOIEEEDPFlt	equ	-36
_LVOIEEEDPFix	equ	-30

	public	_LVOIEEEDPCeil
	public	_LVOIEEEDPFloor
	public	_LVOIEEEDPDiv
	public	_LVOIEEEDPMul
	public	_LVOIEEEDPSub
	public	_LVOIEEEDPAdd
	public	_LVOIEEEDPNeg
	public	_LVOIEEEDPAbs
	public	_LVOIEEEDPTst
	public	_LVOIEEEDPCmp
	public	_LVOIEEEDPFlt
	public	_LVOIEEEDPFix

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

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

CLIB=\
	fmod.$O
ALIB=\
	abs.$O fabs.$O frexp.$O \
	ldexp.$O math.$O modf.$O pabs.$O padd.$O \
	pcmp.$O pdiv.$O pfix.$O pflt.$O pmul.$O \
	pneg.$O psub.$O ptst.$O lvo.$O ceil.$O \
	floor.$O floats.$O

mlib:	$(CLIB) $(ALIB)

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

	public	__abort
	public	_SysBase
	public	_DOSBase
	public	_MathIeeeDoubBasBase

	public	amiga_ieee
amiga_ieee:
	tst.l	_MathIeeeDoubBasBase
	bne	skip
	movem.l	d0-d3/a0/a1/a6,-(sp)
	move.l	_SysBase,a6
	move.l	#0,d0
	lea	libnam,a1
	jsr	-552(a6)	;OpenLibrary
	move.l	d0,_MathIeeeDoubBasBase
	bne	done
	move.l	_DOSBase,a6
	jsr	-60(a6)		;Output
	move.l	d0,d1
	lea	msg,a0
	move.l	a0,d2
	move.l	#16,d3
	jsr	-48(a6)		;Write
	move.l	#1,-(sp)
	jsr	__abort
done
	movem.l	(sp)+,d0-d3/a0/a1/a6
skip
	movem.l	a0/a1,-(sp)
	move.l	8(sp),a0
	move.l	a6,8(sp)
	move.l	_MathIeeeDoubBasBase,a6
	jsr	(a6,a0.l)
	movem.l	(sp)+,a0/a1/a6
	rts

libnam	dc.b	'mathieeedoubbas.library',0

msg	dc.b	'no ieee library',$a
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
pabs.a68
; Copyright (C) 1986 by Manx Software Systems, Inc.
; :ts=8

	public	.Pabs
.Pabs:
	move.l	#_LVOIEEEDPAbs#,-(sp)
	jmp	amiga_ieee#

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

	public	.Padd
.Padd:
	move.l	#_LVOIEEEDPAdd#,-(sp)
	jmp	amiga_ieee#

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

	public	.Pcmp
.Pcmp:
	move.l	#_LVOIEEEDPCmp#,-(sp)
	jmp	amiga_ieee#

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

	public	.Pdiv
.Pdiv:
	move.l	#_LVOIEEEDPDiv#,-(sp)
	jmp	amiga_ieee#

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

	public	.Pfix
.Pfix:
	move.l	#_LVOIEEEDPFix#,-(sp)
	jmp	amiga_ieee#

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

	public	.Pflt
.Pflt:
	move.l	#_LVOIEEEDPFlt#,-(sp)
	jmp	amiga_ieee#

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

	public	.Pmul
.Pmul:
	move.l	#_LVOIEEEDPMul#,-(sp)
	jmp	amiga_ieee#

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

	public	.Pneg
.Pneg:
	move.l	#_LVOIEEEDPNeg#,-(sp)
	jmp	amiga_ieee#

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

	public	.Psub
.Psub:
	move.l	#_LVOIEEEDPSub#,-(sp)
	jmp	amiga_ieee#

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

	public	.Ptst
.Ptst:
	move.l	#_LVOIEEEDPTst#,-(sp)
	jmp	amiga_ieee#