⟦db21c941d⟧

TextFile

#ifndef lint
static	char RCSid[] = "$Header: mf_arith.c,v 1.0 86/01/31 14:59:11 richards Released $";
#endif

/*
 * mf_arith.c -- optimized arithmetic routines for UNIX METAFONT 84
 *
 *	These functions account for a significant percentage of the processing
 *	time used by METAFONT, and have been recoded in assembly language for
 *	some classes of machines.
 *
 *	User BEWARE!  These codings make assumptions as to how the Pascal
 *	compiler passes arguments to external subroutines, which should not
 *	be ignored if you take this code to a different environment than
 *	it was originally designed and tested on!  You have been warned!
 *
 *	function makefraction(p, q: integer): fraction;
 *		{ Calculate the function floor( (p * 2^28) / q + 0.5 )	}
 *		{ if both p and q are positive.  If not, then the value	}
 *		{ of makefraction is calculated as though both *were*	}
 *		{ positive, then the result sign adjusted.		}
 *		{ (e.g. makefraction ALWAYS rounds away from zero)	}
 *		{ In case of an overflow, return the largest possible	}
 *		{ value (2^32-1) with the correct sign, and set global	}
 *		{ variable "aritherror" to 1.  Note that -2^32 is 	}
 *		{ considered to be an illegal product for this type of	}
 *		{ arithmetic!						}
 *
 *	function takefraction(f: fraction; q: integer): integer;
 *		{ Calculate the function floor( (p * q) / 2^28 + 0.5 )	}
 *		{ takefraction() rounds in a similar manner as		}
 *		{   makefraction() above.				}
 *
 */

#define	EXTERN	extern
#include "../mfd.h"

extern	boolean	aritherror;		/* to be set on overflow */

#define	EL_GORDO	0x7fffffff	/* 2^31-1		*/
#define	FRACTION_ONE	0x10000000	/* 2^28			*/
#define	FRACTION_HALF	0x08000000	/* 2^27			*/
#define	FRACTION_FOUR	0x40000000	/* 2^30			*/

#if defined(VAX4_2) || defined(VAX4_3) || defined(PYRAMID)

#if	defined(VAX4_2) || defined(VAX4_3)

/*
 * DEC VAX-11 class systems running 4.2 or 4.3 BSD Unix and Berkeley PC
 *
 * In the VAX assembler code that follows, remember, REMEMBER that
 * a quad result addressed in Rn has most significant longword
 * in Rn+1, NOT Rn!
 */

fraction zmakefraction(p, q)
	integer		p, q;
{
	register	r11;

	asm("	movl	8(ap),r0		"); /* q		      */
	asm("	xorl3	4(ap),r0,r11		"); /* sign of prod	      */
	asm("	jgeq	Lmf1			");
	asm("	mnegl	r0,r0			"); /* abs(q)*sign(p)	      */
	asm("Lmf1:				");
	asm("	divl2	$2,r0			"); /* 	 div 2 (signed)       */
	asm("	emul	4(ap),$0x10000000,r0,r0	"); /* p*2^28+(sign(p)*(q/2)) */
	asm("	ediv	8(ap),r0,r0,r1		"); /*((p*2^28)+	      */
	asm("	jvs	Lmf2			"); /*	sign(p)*(q/2))/q      */
	asm("	mnegl	r0,r1			"); /* check for 2^32	      */
	asm("	jvc	Lmf3			");
	asm("Lmf2:				");
	asm("	movl	$1,_aritherror		");
	asm("	movl	$0x7fffffff,r0		"); /* EL_GORDO		      */
	asm("	tstl	r11			");
	asm("	jgeq	Lmf3			");
	asm("	mnegl	r0,r0			");
	asm("Lmf3:				");
	asm("	ret				");

}

integer ztakefraction(q, f)
	integer		q;
	fraction	f;
{
	register	r11;

	asm("	movl	$0x8000000,r0		"); /* 2^27		    */
	asm("	xorl3	4(ap),8(ap),r11		"); /* sign of prod	    */
	asm("	jgeq	Ltf1			");
	asm("	mnegl	r0,r0			"); /* 2^27*sgn(prod)	    */
	asm("Ltf1:				");
	asm("	emul	4(ap),8(ap),r0,r0	"); /* q*f+sgn(prod)*2^27   */
	asm("	ediv	$0x10000000,r0,r0,r1	"); /*	div 2^28 (signed!)  */
	asm("	jvs	Ltf2			");
	asm("	mnegl	r0,r1			"); /* check for -2^32	    */
	asm("	jvc	Ltf3			");
	asm("Ltf2:				");
	asm("	movl	$1,_aritherror		");
	asm("	movl	$0x7fffffff,r0		"); /* EL_GORDO		    */
	asm("	tstl	r11			");
	asm("	jgeq	Ltf3			");
	asm("	mnegl	r0,r0			");
	asm("Ltf3:				");
	asm("	ret				");
}

#endif	VAX

#if	defined(PYRAMID)

/*
 * Pyramid 90x class of machines, running Pyramid's PASCAL
 */

fraction zmakefraction(p, q)
	integer		p, q;
{

	asm("	movw	$0,lr0		"); /* high word of 64 bit accum    */
	asm("	mabsw	pr0,lr1		"); /* lr1 = abs(p)		    */
	asm("	mabsw	pr1,lr2		"); /* lr2 = abs(q)		    */
	asm("	ashll	$29,lr0		"); /* lr0'lr1 = abs(p)*2^29	    */
	asm("	addw	lr2,lr1		"); /*  abs(p)*2^29 + abs(q)	    */
	asm("	addwc	$0,lr0		");
	asm("	ashrl	$1,lr0		"); /*  (abs(p)*2^28 + abs(q/2))    */
	asm("	ediv	lr2,lr0		"); /*  (abs(p)*2^28 / abs(q)) + .5 */
	asm("	bvc	L1		");
	asm("	movw	$1,_aritherror	");
	asm("	movw	$0x7fffffff,lr0	"); /* EL_GORDO			    */
	asm("L1:			");
	asm("	xorw	pr0,pr1		"); /* determine sign of result	    */
	asm("	blt	L2		");
	asm("	movw	lr0,pr0		");
	asm("	ret			");
	asm("L2:			");
	asm("	mnegw	lr0,pr0		"); /* negative			    */
	asm("	ret			");
}


integer ztakefraction(q, f)
	integer		q;
	fraction	f;
{
	asm("	movw	$0,lr0		"); /* high word of 64 bit accum    */
	asm("	mabsw	pr0,lr1		"); /* lr1 = abs(q)		    */
	asm("	mabsw	pr1,lr2		"); /* lr2 = abs(f)		    */
	asm("	emul	lr2,lr0		"); /* lr0'lr1 = abs(q) * abs(f)    */
	asm("	addw	$0x08000000,lr1	"); /*  (abs(q)*abs(f))+2^27	    */
	asm("	addwc	$0,lr0		");
	asm("	ashll	$4,lr0		"); /* lr0=(abs(q)*abs(f))/2^28 + .5 */
	asm("	bvc	L3		");
	asm("	movw	$1,_aritherror	");
	asm("	movw	$0x7fffffff,lr0	"); /* EL_GORDO			    */
	asm("L3:			");
	asm("	xorw	pr0,pr1		"); /* construct sign of product    */
	asm("	blt	L4		");
	asm("	movw	lr0,pr0		");
	asm("	ret			");
	asm("L4:			");
	asm("	mnegw	lr0,pr0		"); /* negative			    */
	asm("	ret			");
}

#endif	PYRAMID

#else


/*
 * Generic Routines for other machine classes:
 *
 *	essentially a verbatim copy of the WEB
 *	routines, with the hope that the C compiler
 *	may do better than the Pascal on code quality.
 */

fraction zmakefraction(p, q)
    register	integer p, q;
{
    int		negative;
    register	int be_careful;
    register	fraction f;
    int		n;

    if (p < 0) {
	negative = 1;
	p = -p;
    } else negative = 0;

    if (q < 0) {
	negative = !negative;
	q = -q;
    }

    n = p / q;
    p = p % q;

    if (n >= 8) {
	aritherror = true;
	return (negative? -EL_GORDO : EL_GORDO);
    }

    n = (n-1)*FRACTION_ONE;
    f = 1;
    do {
	be_careful = p - q;
	p = be_careful + p;
	if (p >= 0) {
	    f = f + f + 1;
	} else {
	    f <<= 1;
	    p += q;
	}
    } while (f < FRACTION_ONE);

    be_careful = p - q;
    if ((be_careful + p) >= 0)
	f += 1;
    return (negative? -(f+n) : (f+n));

}

integer ztakefraction(q, f)
    register	integer q;
    register	fraction f;
{
    int		n, negative;
    register	int p, be_careful;

    if (f < 0) {
	negative = 1;
	f = -f;
    } else negative = 0;
    if (q < 0) {
	negative = !negative;
	q = -q;
    }

    if (f < FRACTION_ONE)
	n = 0;
    else {
	n = f / FRACTION_ONE;
	f = f % FRACTION_ONE;
    	if (q < (EL_GORDO /  n))
	    n = n * q;
	else {
	    aritherror = true;
	    n = EL_GORDO;
	}
    }

    f += FRACTION_ONE;
    p = FRACTION_HALF;
    if (q < FRACTION_FOUR) 
	do {
	    if (f & 0x1)
		p = (p+q) >> 1;
	    else
		p >>= 1;
	    f >>= 1;
	} while (f > 1);
    else
	do {
	    if (f & 0x1)
		p = p + ((q-p) >> 1);
	    else
		p >>= 1;
	    f >>= 1;
	} while (f > 1);

    be_careful = n - EL_GORDO;
    if ((be_careful + p) > 0) {
	aritherror = true;
	n = EL_GORDO - p;
    }

    return(negative? -(n+p) : (n+p));
}

#endif VAX4_2
DataMuseum.dk

DKUUG/EUUG Conference tapes

⟦db21c941d⟧ TextFile

Derivation

TextFile
