/* 
 * Copyright (c) 1992 Helsinki University of Technology
 * All Rights Reserved.
 * 
 * Permission to use, copy, modify and distribute this software and its
 * documentation is hereby granted, provided that both the copyright
 * notice and this permission notice appear in all copies of the
 * software, derivative works or modified versions, and any portions
 * thereof, and that both notices appear in supporting documentation.
 * 
 * HELSINKI UNIVERSITY OF TECHNOLOGY ALLOWS FREE USE OF THIS SOFTWARE IN
 * ITS "AS IS" CONDITION. HELSINKI UNIVERSITY OF TECHNOLOGY DISCLAIMS ANY
 * LIABILITY OF ANY KIND FOR ANY DAMAGES WHATSOEVER RESULTING FROM THE
 * USE OF THIS SOFTWARE.
 */
/*
 *   Copyright (c) 1985 Regents of the University of California.
 *   All rights reserved. 
 *   
 *   Redistribution and use in source and binary forms are permitted 
 *   provided that the above copyright notice and this paragraph are 
 *   duplicated in all such forms and that any documentation, 
 *   advertising materials, and other materials related to such 
 *   distribution and use acknowledge that the software was developed 
 *   by the University of California, Berkeley.  The name of the 
 *   University may not be used to endorse or promote products derived 
 *   from this software without specific prior written permission. 
 *   THIS SOFTWARE IS PROVIDED ``AS IS'' AND WITHOUT ANY EXPRESS OR 
 *   IMPLIED WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED 
 *   WARRANTIES OF MERCHANTIBILITY AND FITNESS FOR A PARTICULAR PURPOSE. 
 *   
 *   All recipients should regard themselves as participants in an ongoing 
 *   research project and hence should feel obligated to report their 
 *   experiences (good or bad) with these elementary function codes, using 
 *   the sendbug(8) program, to the authors. 
 *   
 *  	@(#)support.s	5.3 (Berkeley) 6/30/88 
 */

/*   IEEE recommended functions 
 *    
 *   double copysign(x,y) 
 *   double x,y; 
 *   IEEE 754 recommended function, return x*sign(y) 
 *   Coded by K.C. Ng in National 32k assembler, 11/9/85. 
 */
/* 	.vers	2 */
	.text
	.align	2
	.globl	_copysign
_copysign:
	movl	4(sp),f0
	movd	8(sp),r0
	movd	16(sp),r1
	xord	r0,r1
	andd	0x80000000,r1
	cmpqd	0,r1
	beq	end_copysign
	negl	f0,f0
end_copysign:	ret	0

/*   */
/*  double logb(x) */
/*  double x; */
/*  IEEE p854 recommended function, return the exponent of x (return float(N)  */
/*  such that 1 <= x*2**-N < 2, even for subnormal number. */
/*  Coded by K.C. Ng in National 32k assembler, 11/9/85. */
/*  Note: subnormal number (if implemented) will be taken care of.  */
/*  */
/* 	.vers	2 */
	.text
	.align	2
	.globl	_logb
_logb:
/*  */
/*  extract the exponent of x */
/*  glossaries:	r0 = high part of x */
/* 		r1 = unbias exponent of x */
/* 		r2 = 20 (first exponent bit position) */
/*  */
	movd	8(sp),r0
	movd	20,r2
	extd	r2,r0,r1,11	/*  extract the exponent of x */
	cmpqd	0,r1		/*  if exponent bits = 0, goto L3 */
	beq	L3logb
	cmpd	0x7ff,r1
	beq	L2logb		/*  if exponent bits = 0x7ff, goto L2 */
L1logb:	subd	1023,r1		/*  unbias the exponent */
	movdl	r1,f0		/*  convert the exponent to floating value */
	ret	0
/*  */
/*  x is INF or NaN, simply return x */
/*  */
L2logb:
	movl	4(sp),f0	/*  logb(+inf)=+inf, logb(NaN)=NaN */
	ret	0
/*  */
/*  x is 0 or subnormal */
/*  */
L3logb:
	movl	4(sp),f0
	cmpl	0f0,f0
	beq	L5logb		/*  x is 0 , goto L5 (return -inf) */
/*  */
/*  Now x is subnormal */
/*  */
	mull	L64,f0		/*  scale up f0 with 2**64 */
	movl	f0,tos
	movd	tos,r0
	movd	tos,r0		/*  now r0 = new high part of x */
	extd	r2,r0,r1,11	/*  extract the exponent of x to r1 */
	subd	1087,r1		/*  unbias the exponent with correction  */
	movdl	r1,f0		/*  convert the exponent to floating value */
	ret	0
/*  */
/*  x is 0, return logb(0)= -INF */
/*  */
L5logb:
	movl	0f1.0e300,f0
	mull	0f-1.0e300,f0	/*  multiply two big numbers to get -INF */
	ret	0
/*   */
/*  double rint(x) */
/*  double x# */
/*  ... delivers integer nearest x in direction of prevailing rounding */
/*  ... mode */
/*  Coded by K.C. Ng in National 32k assembler, 11/9/85. */
/*  Note: subnormal number (if implemented) will be taken care of.  */
/*  */
/* 	.vers	2 */
	.text
	.align	2
	.globl	_rint
_rint:
/*  */
	movd	8(sp),r0
	movd	20,r2
	extd	r2,r0,r1,11	/*  extract the exponent of x */
	cmpd	0x433,r1
	ble	itself
	movl	L52rint,f2		/*  f2 = L = 2**52 */
	cmpqd	0,r0
	ble	L1rint
	negl	f2,f2		/*  f2 = s = copysign(L,x) */
L1rint:	addl	f2,f0		/*  f0 = x + s */
	subl	f2,f0		/*  f0 = f0 - s */
	ret	0
itself:	movl	4(sp),f0
	ret	0
L52rint:	.double	0x0,0x43300000	/*  L52=2**52 */
/*   */
/*  int finite(x) */
/*  double x# */
/*  IEEE 754 recommended function, return 0 if x is NaN or INF, else 0 */
/*  Coded by K.C. Ng in National 32k assembler, 11/9/85. */
/*  */
/* 	.vers	2 */
	.text
	.align	2
	.globl	_finite
_finite:
	movd	4(sp),r1
	andd	0x800fffff,r1
	cmpd	0x7ff00000,r1
	sned	r0		/*  r0=0 if exponent(x) = 0x7ff */
	ret	0
/*   */
/*  double scalb(x,N) */
/*  double x; int N; */
/*  IEEE 754 recommended function, return x*2**N by adjusting  */
/*  exponent of x. */
/*  Coded by K.C. Ng in National 32k assembler, 11/9/85.  */
/*  Replaced by native 32381 scalb from libc (ldexp.s) jtp */
/*  Note: subnormal number (if implemented) will be taken care of  */
/*  */
/* 	.vers	2 */
	.text
	.align	2
	.globl	_errno
	.globl	_scalb
/* kivinen's code from libc */
_scalb:
        enter   [],8
        movd    12(fp),r1       /* get value, high 32 bit */
        lshd    -20,r1          /* extract exponent */
        andd    0x7ff,r1        /* 11 lower bits */
        subd    1023,r1         /* unbias exponent */
        movd    16(fp),r0       /* get exp */
        addd    r1,r0           /* add exponents */
        cmpd    r0,1023         /* most positive exponent */
        bgt     ldexp_overflow  /* too large */
        cmpd    r0,-1022        /* most negative exponent */
        blt     ldexp_underflow /* too small */
        addd    1023,r0         /* bias exponent */
        lshd    20,r0           /* convert exponent back to its place */
        movd    12(fp),r1       /* get value high 32 bit */
        bicd    0x7ff00000,r1   /* clear exponent */
        ord     r0,r1           /* insert exponent */
        movd    r1,-4(fp)       /* put high 32 bit to local variable */
        movd    8(fp),-8(fp)    /* value low 32 bit to local variable */
        movl    -8(fp),f0       /* return value * 2**exp */
        exit    []
        ret     0               /* neither */
ldexp_underflow:
        movdl   0d0e0,f0        /* if underflow return 0, set errno */
        movd    34,@_errno	/* ERANGE */
        exit    []
        ret     0
ldexp_overflow:
        movl    @huge,f0        /* if overflow return HUGE */
        movdl   0d0e0,f2
        cmpl    f2,8(fp)        /* check original sign */
        bgt     ldexp_positive
        negl    f0,f0           /* if negative, return -HUGE */
ldexp_positive:
        movd    34,@_errno	/* set errno ERANGE */
        exit    []
        ret     0

huge:   .long   0xffffffff      /* the largest number that can */
        .long   0x7fefffff      /* be represented in a long floating */
                                /* number.  This is given in hex in order */
                                /* to avoid floating conversions */
/*  */
/*  if x=0 return 0 */
/*  */
/* 	movl	4(sp),f0 */
/* 	cmpl	0f0,f0 */
/* 	beq	end_scalb		# scalb(0,N) is x itself */
/*  */
/*  extract the exponent of x */
/*  glossaries:	r0 = high part of x,  */
/* 		r1 = unbias exponent of x, */
/* 		r2 = 20 (first exponent bit position). */
/*  */
/* 	movd	8(sp),r0	# r0 = high part of x */
/* 	movd	20,r2		# r2 = 20 */
/* 	extd	r2,r0,r1,11	# extract the exponent of x in r1 */
/* 	cmpd	0x7ff,r1	 */
/*  */
/*  if exponent of x is 0x7ff, then x is NaN or INF# simply return x   */
/*  */
/* 	beq	end_scalb		 */
/* 	cmpqd	0,r1 */
/*  */
/*  if exponent of x is zero, then x is subnormal# goto L19 */
/*  */
/* 	beq	L19scalb		 */
/* 	addd	12(sp),r1	# r1 = (exponent of x) + N */
/* 	bfs	inof		# if integer overflows, goto inof */
/* 	cmpqd	0,r1		# if new exponent <= 0, goto underflow */
/* 	bge	underflow */
/* 	cmpd	2047,r1		# if new exponent >= 2047 goto overflow */
/* 	ble	overflow */
/* 	insd	r2,r1,r0,11	# insert the new exponent  */
/* 	movd	r0,tos */
/* 	movd	8(sp),tos */
/* 	movl	tos,f0		# return x*2**N */
/* end_scalb:	ret	0 */
/* inof:	bcs	underflow	# negative int overflow if Carry bit is set */
/* overflow: */
/* 	andd	0x80000000,r0	# keep the sign of x */
/* 	ord	0x7fe00000,r0	# set x to a huge number */
/* 	movd	r0,tos */
/* 	movqd	0,tos */
/* 	movl	tos,f0 */
/* 	mull	0f1.0e300,f0	# multiply two huge number to get overflow */
/* 	ret	0 */
/* underflow: */
/* 	addd	64,r1		# add 64 to exonent to see if it is subnormal */
/* 	cmpqd	0,r1 */
/* 	bge	zero		# underflow to zero */
/* 	insd	r2,r1,r0,11	# insert the new exponent  */
/* 	movd	r0,tos */
/* 	movd	8(sp),tos */
/* 	movl	tos,f0 */
/* 	mull	L30,f0		# readjust x by multiply it with 2**-64 */
/* 	ret	0 */
/* zero:	andd	0x80000000,r0	# keep the sign of x */
/* 	ord	0x00100000,r0	# set x to a tiny number */
/* 	movd	r0,tos */
/* 	movqd	0,tos */
/* 	movl	tos,f0 */
/* 	mull	0f1.0e-300,f0	# underflow to 0  by multipling two tiny nos. */
/* 	ret	0 */
/* L19scalb:		# subnormal number */
/* 	mull	L32,f0		# scale up x by 2**64 */
/* 	movl	f0,tos */
/* 	movd	tos,r0 */
/* 	movd	tos,r0		# get the high part of new x */
/* 	extd	r2,r0,r1,11	# extract the exponent of x in r1 */
/* 	addd	12(sp),r1	# exponent of x + N */
/* 	subd	64,r1		# adjust it by subtracting 64 */
/* 	cmpqd	0,r1 */
/* 	bge	underflow */
/* 	cmpd	2047,r1 */
/* 	ble	overflow */
/* 	insd	r2,r1,r0,11	# insert back the incremented exponent  */
/* 	movd	r0,tos */
/* 	movd	8(sp),tos */
/* 	movl	tos,f0 */
/* end:	ret	0 */
/* L30:	.double	0x0,0x3bf00000	# floating point 2**-64 */
/* L32:	.double	0x0,0x43f00000	# floating point 2**64 */
/*   */
/*  double drem(x,y) */
/*  double x,y# */
/*  IEEE double remainder function, return x-n*y, where n=x/y rounded to  */
/*  nearest integer (half way case goes to even). Result exact. */
/*  Coded by K.C. Ng in National 32k assembly, 11/19/85. */
/*  */
/* 	.vers	2 */
	.text
	.align	2
	.globl	_drem
_drem:
/*  */
/*  glossaries:	 */
/* 		r2 = high part of x */
/* 		r3 = exponent of x */
/* 		r4 = high part of y */
/* 		r5 = exponent of y */
/* 		r6 = sign of x */
/* 		r7 = constant 0x7ff00000 */
/*  */
/*   16(fp) : y */
/*    8(fp) : x */
/*  -12(fp) : adjustment on y when y is subnormal */
/*  -16(fp) : fsr */
/*  -20(fp) : nx */
/*  -28(fp) : t */
/*  -36(fp) : t1 */
/*  -40(fp) : nf */
/*   */
/*  */
	enter	[r3,r4,r5,r6,r7],40
	movl	f6,tos
	movl	f4,tos
	movl	0f0,-12(fp)
	movd	0,-20(fp)
	movd	0,-40(fp)
	movd	0x7ff00000,r7	/*  initialize r7=0x7ff00000 */
	movd	12(fp),r2	/*  r2 = high(x) */
	movd	r2,r3
	andd	r7,r3		/*  r3 = xexp */
	cmpd	r7,r3
/*  if x is NaN or INF goto L1 */
	beq	L1drem		
	movd	20(fp),r4
	bicd	0x80000000,r4	/*  [31] r4 = high part of |y| */
	movd	r4,20(fp)	/*  y = |y| */
	movd	r4,r5
	andd	r7,r5		/*  r5 = yexp */
	cmpd	r7,r5
	beq	L2drem		/*  if y is NaN or INF goto L2 */
	cmpd	0x04000000,r5	/*   */
	bgt	L3drem		/*  if y is tiny goto L3 */
/*  */
/*  now y != 0 , x is finite */
/*  */
L10:
	movd	r2,r6
	andd	0x80000000,r6	/*  r6 = sign(x) */
	bicd	0x80000000,r2	/*  [31] x <- |x| */
	sfsr	r1
	movd	r1,-16(fp)	/*  save fsr in -16(fp) */
	bicd	0x20,r1		/* [5] */
	lfsr	r1		/*  disable inexact interupt */
	movd	16(fp),r0	/*  r0 = low part of y */
	movd	r0,r1		/*  r1 = r0 = low part of y */
	andd	0xf8000000,r1	/*  mask off the lsb 27 bits of y */

	movd	r2,12(fp)	/*  update x to |x| */
	movd	r0,-28(fp)	/*   */
	movd	r4,-24(fp)	/*  t  = y */
	movd	r4,-32(fp)	/*   */
	movd	r1,-36(fp)	/*  t1 = y with trialing 27 zeros */
	movd	0x01900000,r1	/*  r1 = 25 in exponent field */
LOOP:
	movl	8(fp),f0	/*  f0 = x */
	movl	16(fp),f2	/*  f2 = y */
	cmpl	f0,f2
	ble	fnad		/*  goto fnad (final adjustment) if x <= y */
	movd	r4,-32(fp)
	movd	r3,r0
	subd	r5,r0		/*  xexp - yexp */
	subd	r1,r0		/*  r0 = xexp - yexp - m25 */
	cmpqd	0,r0		/*  r0 > 0 ? */
	bge	1f
	addd	r4,r0		/*  scale up (high) y */
	movd	r0,-24(fp)	/*  scale up t */
	movl	-28(fp),f2	/*  t */
	movd	r0,-32(fp)	/*  scale up t1 */
1:
	movl	-36(fp),f4	/*  t1 */
	movl	f0,f6
	divl	f2,f6		/*  f6 = x/t */
	floorld	f6,r0		/*  r0 = [x/t] */
	movdl	r0,f6		/*  f6 = n */
	subl	f4,f2		/*  t = t - t1 (tail of t1) */
	mull	f6,f4		/*  f4 = n*t1	...exact */
	subl	f4,f0		/*  x = x - n*t1 */
	mull	f6,f2		/*  n*(t-t1)	...exact */
	subl	f2,f0		/*  x = x - n*(t-t1) */
/*  update xexp */
	movl	f0,8(fp)
	movd	12(fp),r3	
	andd	r7,r3
	jump	LOOP
fnad:
	cmpqd	0,-20(fp)	/*  0 = nx? */
	beq	final
	mull	-12(fp),8(fp)	/*  scale up x the same amount as y */
	movd	0,-20(fp)
	movd	12(fp),r2
	movd	r2,r3
	andd	r7,r3		/*  update exponent of x */
	jump	LOOP

final:
	movl	16(fp),f2	/*  f2 = y (f0=x, r0=n) */
	subd	0x100000,r4	/*  high y /2 */
	movd	r4,-24(fp)
	movl	-28(fp),f4	/*  f4 = y/2 */
	cmpl	f0,f4		/*  x > y/2 ? */
	bgt	1f
	bne	2f
	andd	1,r0		/*  n is odd or even */
	cmpqd	0,r0
	beq	2f
1:
	subl	f2,f0		/*  x = x - y */
2:
	cmpqd	0,-40(fp)
	beq	3f
	divl	-12(fp),f0	/*  scale down the answer */
3:
	movl	f0,tos
	xord	r6,tos
	movl	tos,f0
	movd	-16(fp),r0
	lfsr	r0		/*  restore the fsr */

end_drem:	movl	tos,f4
	movl	tos,f6
	exit	[r3,r4,r5,r6,r7]
	ret	0
/*  */
/*  y is NaN or INF */
/*  */
L2drem:	
	movd	16(fp),r0	/*  r0 = low part of y */
	andd	0xfffff,r4	/*  r4 = high part of y & 0x000fffff */
	ord	r4,r0
	cmpqd	0,r0
	beq	L4drem
	movl	16(fp),f0	/*  y is NaN, return y */
	jump	end_drem
L4drem:	movl	8(fp),f0	/*  y is inf, return x */
	jump	end_drem
/*  */
/*  exponent of y is less than 64, y may be zero or subnormal */
/*  */
L3drem:
	movl	16(fp),f0
	cmpl	0f0,f0
	bne	L5drem
	divl	f0,f0		/*  y is 0, return NaN by doing 0/0 */
	jump	end_drem
/*  */
/*  subnormal y or tiny y */
/*  */
L5drem:	
	movd	0x04000000,-20(fp)	/*  nx = 64 in exponent field */
	movl	L64,f2
	movl	f2,-12(fp)
	mull	f2,f0
	cmpl	f0,LTINY
	bgt	L6drem
	mull	f2,f0
	addd	0x04000000,-20(fp)	/*  nx = nx + 64 in exponent field */
	mull	f2,-12(fp)
L6drem:
	movd	-20(fp),-40(fp)
	movl	f0,16(fp)
	movd	20(fp),r4
	movd	r4,r5
	andd	r7,r5		/*  exponent of new y */
	jump	L10drem
/*  */
/*  x is NaN or INF, return x-x */
/*  */
L1drem:
	movl	8(fp),f0
	subl	f0,f0		/*  if x is INF, then INF-INF is NaN */
	ret	0
L64:	.long 0x0,0x43f00000	/*  L64 = 2**64 */
LTINY:	.long 0x0,0x04000000	/*  LTINY = 2**-959 */
