*
* $VER: 64bit.a 2.2 (15.2.95)
*
* This file contains several routines as replacements for `inner loop'
* algorithms of MetaFont 2.71 and MetaPost 0.62 written in assembly
* language according to suggestions made by Donald Knuth [Computers &
* Typesetting D, pp.39-45].  These routines are appropriate for the
* Motorola family of 680X0 processors, where `X' stands for `2, 3, 4,
* and higher'.
*
* Version 1 of these routines was written by Stefan Scherer.
* Version 2 was modified for use of registers instead of stack parameters.
*
* Copyright (c) 1993,1994 Stefan Scherer.
* Copyright (c) 1994,1995 Andreas Scherer.
*
* Changes:
*    30.11.94: Function arguments are passed in registers instead of via
*              the stack.  This change was suggested by Giuseppe Ghibò
*              one day after the archive was shipped to the Internet.
*              Sorry.  But I promise you don't loose much time anyway.
*              As the complete interface has changed, new version (V2.0)
*    26.01.95: And I really should tell you in this file, what the new
*              interface looks like, so (V2.1)
*    15.02.95: We use these routines in John Hobby's MetaPost as well.
*
		SECTION	text,CODE

		XREF	_aritherror

EL_GORDO	EQU	$7fffffff

******* 64bit/makef *********************************************************
*
*   NAME
*	makef -- machine dependent replacement for the tight loop in
*	the make_fraction() algorithm [module 107].
*
*   SYNOPSIS
*	f = makef( p, q )
*
*   C PROTOTYPE
*	fraction __asm makef(
*          register __d0 integer,
*          register __d1 integer );
*
*   FUNCTION
*	Produces the fraction equivalent of p/q, given integers p and q;
*	it computes the integer f according to the formula
*
*		f = floor ( 2^28 · p / q + 0.5),
*
*	when p and q both are positive.
*
*   INPUTS
*	p - unsigned long integer (32 bit), expected in register D0
*	q - unsigned long integer (32 bit), expected in register D1
*
*   RESULT
*	f - fraction equivalent of p/q (32 bit)
*
*   NOTE
*	This code uses 68020 instructions!
*
*   SEE ALSO
*	takef() -- ``backward'' algorithm
*
*****************************************************************************

		XDEF	_makef

_makef		move.l	d2,-(sp)		; save register

		move.l	d0,d2			; 2^28 · p (64 bit)
		lsr.l	#4,d2			; calculate 32 high bits
		and.l	#%1111,d0
		ror.l	#4,d0			; calculate 32 low bits

		divu.l	d1,d2:d0		; (2^28·p) / q (32 bit)

		lsl.l	#1,d2			; double division rest
		cmp.l	d1,d2			; division rest >= q/2 ?
		bcs.s	make_f_done		; no, we're finished
		addq.l	#1,d0			; yes, increase result

make_f_done	move.l	(sp)+,d2		; restore register
		rts				; return to sender

******* 64bit/takefraction **************************************************
*
*   NAME
*	takefraction -- machine dependent replacement for the
*	take_fraction() algorithm [module 109].
*
*   SYNOPSIS
*	p = takefraction( q, f )
*
*   C PROTOTYPE
*	integer __asm takefraction(
*          register __d1 integer,
*          register __d0 fraction );
*
*       Yes, __d0 and __d1 are twisted!
*
*   FUNCTION
*	Multiplies a given integer q by a fraction f.  When the operands
*	are positive, it computes the integer p according to the formula
*
*		p = floor ( q · f / 2^28 + 0.5),
*
*	a symmetric function of q and f.
*
*   INPUTS
*	q - unsigned long integer (32 bit), expected in register D1
*	f - unsigned fraction value (32 bit), expected in register D0
*
*   RESULT
*	p - unsigned long integer equivalent of qf (32 bit)
*
*   NOTES
*	This code uses 68020 instructions!
*
*	Used registers:
*	d0 - f and result p
*	d1 - q
*	d2 - n
*	d3 - negative
*	d4 - scratch register
*
*   SEE ALSO
*	makef() -- ``forward'' algorithm
*
*****************************************************************************

		XDEF	_takefraction

_takefraction	movem.l	d2-d4,-(sp)		; save register

		moveq	#0,d3			; negative = false
		tst.l	d0			; f >= 0 ?
		bge.s	take_f_pos		; yes
		neg.l	d0			; else, negate f
		moveq	#1,d3			; negative = true
		
take_f_pos	tst.l	d1			; q >= 0 ?
		bge.s	take_q_pos		; yes
		neg.l	d1			; else, negate q
		bchg	#0,d3			; negative = ! negative

take_q_pos	cmp.l	#$10000000,d0		; f < $10000000 ?
		bcs.s	take_f_less		; yes
		move.l	d0,d2			; n = f
		rol.l	#4,d2
		and.l	#%1111,d2		; n = f / $10000000
		and.l	#$0fffffff,d0		; f = f % $10000000
		move.l	#EL_GORDO,d4
		divu.l	d2,d4			; EL_GORDO / n
		cmp.l	d4,d1			; q <= (EL_GORDO/n) ?
		bls.s	take_q_less
		move.b	#1,_aritherror(a4)	; aritherror = true !
		move.l	#EL_GORDO,d2		; n = EL_GORDO
		bra.s	take_if_f_done

take_q_less	mulu.l	d1,d2			; n = n · q
		bra.s	take_if_f_done

take_f_less	moveq	#0,d2			; n = 0

take_if_f_done	mulu.l	d1,d4:d0		; multiply q by f
		divu.l	#$10000000,d4:d0	; (q·f) / 2^28
		cmp.l	#$08000000,d4		; division rest >= (2^28)/2 ?
		bcs.s	takef_f_done		; no, we're finished
		addq.l	#1,d0			; yes, increase result by 1

takef_f_done	move.l	d2,d4
		sub.l	#EL_GORDO,d4		; becareful = n - EL_GORDO
		add.l	d0,d4			; becareful + p > 0 ?
		tst.l	d4
		ble.s	take_careful		; no
		move.b	#1,_aritherror(a4)	; aritherror = true
		move.l	#EL_GORDO,d0		; return EL_GORDO
		bra.s	take_pos

take_careful	add.l	d2,d0			; result = n + p
		tst.l	d3			; negative ?
		beq.s	take_pos		; no, positive
		neg.l	d0			; yes, invert result

take_pos	movem.l	(sp)+,d2-d4		; restore register
		rts				; return to sender

******* 64bit/makes *********************************************************
*
*   NAME
*	makes -- machine dependent replacement for the tight loop in
*	the make_scaled() algorithm [module 114].
*
*   SYNOPSIS
*	s = makes( p, q )
*
*   C PROTOTYPE
*	scaled __asm makes(
*          register __d0 integer,
*          register __d1 integer );
*
*   FUNCTION
*	Produces the scaled number equivalent of p/q, given integers p and q;
*	it computes the integer s according to the formula
*
*		s = floor ( 2^16 · p / q + 0.5),
*
*	when p and q both are positive.
*
*   INPUTS
*	p - unsigned long integer (32 bit), expected in register D0
*	q - unsigned long integer (32 bit), expected in register D1
*
*   RESULT
*	s - scaled number equivalent of p/q (32 bit)
*
*   NOTE
*	This code uses 68020 instructions!
*
*   SEE ALSO
*	takes() -- ``backward'' algorithm
*
*****************************************************************************

		XDEF	_makes

_makes		move.l	d2,-(sp)		; save register

		swap	d0			; 2^16 · p (64 bit)
		moveq	#0,d2
		move.w	d0,d2
		clr.w	d0

		divu.l	d1,d2:d0		; (2^16·p) / q (32 bit)

		lsl.l	#1,d2			; double division rest
		cmp.l	d1,d2			; division rest >= q/2 ?
		bcs.s	make_s_done		; no, we're finished
		addq.l	#1,d0			; yes, increase result

make_s_done	move.l	(sp)+,d2		; restore register
		rts				; return to sender

******* 64bit/takes **********************************************************
*
*   NAME
*	takes -- machine dependent replacement for the tight loop in
*	the take_scaled() algorithm [module 112].
*
*   SYNOPSIS
*	p = takes( q, s )
*
*   C PROTOTYPE
*	integer __asm takes(
*          register __d0 integer,
*          register __d1 scaled );
*
*   FUNCTION
*	Multiplies a given integer q by a scaled quantity s.  When the
*	operands are positive, it computes the integer p according to
*	the formula
*
*		p = floor ( q · s / 2^16 + 0.5),
*
*	a symmetric function of q and s.
*
*   INPUTS
*	q - unsigned long integer (32 bit), expected in register D0
*	s - unsigned scaled quantity (32 bit), expected in register D1
*
*   RESULT
*	p - unsigned long integer equivalent of qs (32 bit)
*
*   NOTE
*	This code uses 68020 instructions!
*
*   SEE ALSO
*	makes() -- ``forward'' algorithm
*
*****************************************************************************

		XDEF	_takes

_takes		move.l	d2,-(sp)		; save register

		mulu.l	d1,d2:d0		; multiply q by f
		divu.l	#$00010000,d2:d0	; (q·f) / 2^16

		cmp.l	#$00008000,d2		; division rest >= (2^16)/2 ?
		bcs.s	take_s_done		; no, we're finished
		addq.l	#1,d0			; yes, increase result

take_s_done	move.l	(sp)+,d2		; restore register
		rts				; return to sender

		END