mathhardalib.s

vxworks source code, used for develop vxworks system.

	j	ra

90:	/* raise Overflow and return +Infinity */
	mfc1	t0, $f13			/* extract argument exponent */
	sll	t0, 1
	srl	t0, 20+1
	beq	t0, 2047, 91f			/* - if NaN or Infinity */
	li.d	$f0, 0.898846567431158e308
	add.d	$f0, $f0			/* raise Overflow */
	j	ra
91:	mov.d	$f0, $f12			/* dont raise any flags */
						/* - if argument is exceptional */
						/* just return argument */
	j	ra
	.end expm1

/*******************************************************************************
*
* trunc - floating point truncation
*

* INTERNAL
* An alternate algorithm would to check for numbers < 2**53,
* set the rounding mode, add 2**53, and subtract 2**53. 

* double trunc (dblParam)
*     double dblParam;

*/
	.globl trunc
	.ent   trunc
trunc:
	.frame	sp, 0, ra
	mfc1	t1, $f13
	mfc1	t0, $f12
	srl	t2, t1, 20
	and	t2, 0x7FF
	sub	t2, 1023
	bge	t2, 0, trunc1
	mtc1	$0, $f0
	mtc1	$0, $f1
	j	ra
trunc1:
	sub	t2, 20
	bgt	t2, 0, trunc2
	neg	t2
	srl	t1, t2
	sll	t1, t2
	mtc1	$0, $f0
	mtc1	t1, $f1
	j	ra
trunc2:
	sub	t2, 32
	bge	t2, 0, trunc3
	neg	t2
	srl	t0, t2
	sll	t0, t2
trunc3:
	mtc1	t0, $f0
	mtc1	t1, $f1
	j	ra
	.end trunc


/*******************************************************************************
*
* z_abs_ - 
*
* NOMANUAL
*
* This entrypoint provided for FCOM for ZABS (complex absolute value),
* because FCOM has a hard time calling HYPOT directly.  Also used by
* FCOM when user writes INTRINSIC ZABS.  The latter must use pass by
* reference, of course. 
*/
	.globl z_abs_
	.ent   z_abs_
z_abs_:
	.frame	sp, 0, ra
	l.d	$f12, 0(a0)
	l.d	$f14, 8(a0)
	/* just fall through */
	.end z_abs_

/*******************************************************************************
*
* hypot - 
*
* NOMANUAL
*
*/
	.globl hypot
	.ent   hypot
hypot:
	.frame	sp, 0, ra
	.set noreorder
	abs.d	$f2, $f12
	abs.d	$f4, $f14
	mfc1	t0, $f3
	mfc1	t1, $f5
	li	t7, 2047
	srl	t2, t0, 20
	srl	t3, t1, 20
	beq	t2, t7, 70f
	 c.lt.d	$f2, $f4
	beq	t3, t7, 75f
	 subu	t4, t2, t3
	bc1f	10f
	 slt	t5, t4, 31
	abs.d	$f2, $f14
	abs.d	$f4, $f12
	move	t1, t0
	subu	t4, t3, t2
	slt	t5, t4, 31
10:	beq	t5, 0, 20f
	 mfc1	t0, $f4
11:	bne	t1, 0, 12f
	 sub.d	$f6, $f2, $f4
	beq	t0, 0, 20f
	 nop
12:
	.set reorder
	SETFRAME(hypot,4)
	subu	sp, FRAMESZ(hypot)
	SW	ra, FRAMERA(hypot)(sp)
	s.d	$f2, FRAMEA0(hypot)(sp)
	s.d	$f4, FRAMEA1(hypot)(sp)
	c.lt.d	$f4, $f6
	bc1f	13f

 	div.d	$f6, $f2, $f4
	li.d	$f2, one
	mul.d	$f12, $f6, $f6
	s.d	$f6, FRAMER0(hypot)(sp)
	add.d	$f12, $f2
	jal	sqrt
	l.d	$f6, FRAMER0(hypot)(sp)
	add.d	$f6, $f0
	b	14f
13:
	li.d	$f10, two
	div.d	$f6, $f4
	add.d	$f8, $f6, $f10
	mul.d	$f8, $f6
	s.d	$f6, FRAMER0(hypot)(sp)
	add.d	$f12, $f8, $f10
	s.d	$f8, FRAMER2(hypot)(sp)
	jal	sqrt
	li.d	$f10, HYPT_sqrt2
	l.d	$f6, FRAMER0(hypot)(sp)
	l.d	$f8, FRAMER2(hypot)(sp)
	add.d	$f0, $f10
	div.d	$f8, $f0
	add.d	$f6, $f8
	li.d	$f10, HYPT_r2p1lo
	li.d	$f12, HYPT_r2p1hi
	add.d	$f6, $f10
	add.d	$f6, $f12
14:
	l.d	$f4, FRAMEA2(hypot)(sp)
	LW	ra, FRAMERA(hypot)(sp)
	div.d	$f6, $f4, $f6
	l.d	$f2, FRAMEA0(hypot)(sp)
	addu	sp, FRAMESZ(hypot)
	add.d	$f0, $f6, $f2
	j	ra

20:	mov.d	$f0, $f2
	j	ra

22:	mov.d	$f0, $f4
	j	ra

70:	c.eq.d	$f12, $f12
	bc1t	20b
	beq	t3, t7, 75f
	mov.d	$f0, $f12
	j	ra

75:	c.eq.d	$f14, $f14
	bc1t	22b
	mov.d	$f0, $f14
	j	ra
	.end hypot


/*******************************************************************************
*
* log1p - 
*

* double log1p (dblParam)
*     double dblParam;

*/
	.globl log1p
	.ent   log1p
log1p:
	.frame	sp, 0, ra
	/* argument in f12 */
	.set noreorder
	cfc1	t5, $31
	and	t6, t5, -4
	li.d	$f10,  6.4494458917859432e-02	# ceil( exp( 1/16)-1)
	li.d	$f16, -6.0586937186524220e-02	# floor(exp(-1/16)-1)
	c.ult.d	$f12, $f10
	li.d	$f10, 1.0
	bc1f	1f
	 c.olt.d $f16, $f12
	neg.d	$f16, $f10
	bc1t	5f
	 c.ule.d $f12, $f16
	nop
	bc1t	8f
	 nop
1:	add.d	$f14, $f12, $f10
	li	t4, 2047
	mfc1	t0, $f15
	ctc1	t6, $31
	srl	t1, t0, 20
	beq	t1, t4, 7f
	 subu	t1, 1023
	.set reorder
	sll	t2, t1, 20
	subu	t9, t0, t2
	mtc1	t9, $f15
	li.d	$f16, 3.5184372088832000e+13	# 2^(53-8)
	mtc1	t1, $f8
	add.d	$f0, $f14, $f16
	la	t4, LOGC_TBL
	sub.d	$f18, $f0, $f16
	mfc1	t3, $f0

	slt	v0, t1, -1
	slt	v1, t1, 53
	beq	v0, 0, 2f
	sub.d	$f12, $f14, $f18
	b	4f
2:
	mfc1	t9, $f13
	subu	t9, t2
	mtc1	t9, $f13
	li	t8, 1023<<20
	subu	t8, t2
	mtc1	t8, $f1
	mtc1	$0, $f0
	bne	v1, 0, 3f
	sub.d	$f12, $f18
	add.d	$f12, $f0
	b	4f
3:
	sub.d	$f0, $f18
	add.d	$f12, $f0
4:
	add.d	$f12, $f12
	add.d	$f18, $f14
	div.d	$f12, $f18

	cvt.d.w	$f8, $f8
	l.d	$f10, 128*16+0(t4)	# log2head
	l.d	$f16, 128*16+8(t4)	# log2trail
	mul.d	$f0, $f8, $f10
	sll	t3, 4
	addu	t3, t4
	mul.d	$f2, $f8, $f16
	l.d	$f4, -128*16+0(t3)
	l.d	$f6, -128*16+8(t3)
	add.d	$f0, $f4
	add.d	$f2, $f6

	mul.d	$f18, $f12, $f12
	li.d	$f10, 1.2500053168098584e-02
	li.d	$f16, 8.3333333333039133e-02
	mul.d	$f14, $f18, $f10
	add.d	$f14, $f16
	mul.d	$f14, $f18
	mul.d	$f14, $f12
	add.d	$f14, $f2
	add.d	$f14, $f12
	ctc1	t5, $31
	add.d	$f0, $f14
	j	ra

5:	/* exp(-1/16)-1 < x < exp(1/16)-1 */
	/* use special approximation */
	li.d	$f16, 1.57e-16
	abs.d	$f14, $f12
	c.lt.d	$f14, $f16
	bc1t	7f
	ctc1	t6, $31
	add.d	$f16, $f10, $f10	# 2.0
	add.d	$f14, $f12, $f16
	div.d	$f14, $f10, $f14
	mul.d	$f2, $f12, $f14
	add.d	$f2, $f2
	mul.d	$f4, $f2, $f2
	li.d	$f10, 4.3488777770761457e-04
	li.d	$f16, 2.2321399879194482e-03
	mul.d	$f6, $f4, $f10
	add.d	$f6, $f16
	li.d	$f10, 1.2500000003771751e-02
	mul.d	$f6, $f4
	add.d	$f6, $f10
	li.d	$f16, 8.3333333333331788e-02
	mul.d	$f6, $f4
	add.d	$f6, $f16
	mul.d	$f6, $f4
	mul.d	$f6, $f2
	cvt.s.d	$f4, $f2
	cvt.d.s	$f4, $f4
	cvt.s.d	$f2, $f12
	cvt.d.s	$f2, $f2
	sub.d	$f8, $f12, $f2
	sub.d	$f0, $f12, $f4
	add.d	$f0, $f0
	mul.d	$f2, $f4
	sub.d	$f0, $f2
	mul.d	$f8, $f4
	sub.d	$f0, $f8
	mul.d	$f0, $f14
	add.d	$f0, $f6
	ctc1	t5, $31
	add.d	$f0, $f4
	j	ra

7:	/* log(+Infinity) = +Infinity */
	/* log(NaN) = NaN */
	mov.d	$f0, $f12
	j	ra

8:	/* x <= -1 or x = NaN */
	c.eq.d	$f12, $f16
	li.d	$f10, 0.0
	bc1f	9f
	div.d	$f0, $f16, $f10
	j	ra

9:	c.eq.d	$f12, $f12
	bc1f	7b
	div.d	$f0, $f10, $f10
	j	ra
	.end log1p

/********************************************************************************
* rint - convert double to integer
*
 
* int rint (dblParam)
*     double dblParam;
 
*/
	.globl rint
	.ent   rint
rint:
	.frame	sp, 0, ra
	li.d	$f4, 4503599627370496.0	/* 2^52 */
	abs.d	$f2, $f12		/* |arg| */
	c.olt.d	$f2, $f4		/* if |arg| >= 2^52 or arg is NaN */
	mfc1	t0, $f13
	mov.d	$f0, $f12
	bc1f	4f			/* then done */
	/* < 2^52 */
	sll	t1, t0, 1
	bgez	t0, 2f			/* if input negative, negate result */
	/* negative */
	beq	t1, 0, 3f		/* possible -0 */
1:	sub.d	$f0, $f12, $f4
	add.d	$f0, $f4
	j	ra
2:	/* positive */
	add.d	$f0, $f12, $f4		/* bias by 2^52 to force non-integer
					   bits off end */
	sub.d	$f0, $f4		/* unbias */
	j	ra

3:	/* msw = 80000000 */
	mfc1	t1, $f12		/* if -0, return -0 */
	bne	t1, 0, 1b		/* if negative denorm, process that */
4:	j	ra
	.end rint
		
/******************************************************************************
*
* isNaN - 
*
* This routine takes an input double-precision floating point
* parameter, <dbl>, and returns TRUE if <dbl> contains an IEEE not
* a number value.
*
*
* BOOL isNaN (dbl)
*     double	dbl;
*
* RETURNS: TRUE if <dbl> is NaN, else FALSE
*/
	.ent	isNaN
isNaN:
	mfc1	v0, fp13
	srl	v0, v0, 16
	li	t0, INF_HI_SHORT
	andi	v0, v0, 0x7fff
	swc1	fp12, 4(sp)
	slt	v0, t0, v0
	j	ra
	.end	isNaN

#endif	/* INCLUDE_DOUBLE_PRECISION */


#ifdef  INCLUDE_SINGLE_PRECISION
/*******************************************************************************
*
* acosf - floating point arc-cosine
*

* float acos (fltParam)
*     float fltParam;

*/
	.globl acosf
	.ent   acosf
acosf:
	.frame	sp, 0, ra
	li.s	$f8, half
	abs.s	$f14, $f12
	c.le.s	$f14, $f8
	li.s	$f10, FACS_eps
	move	t9, ra
	bc1f	acosf2
	c.lt.s	$f14, $f10
	mov.s	$f0, $f12
	bc1t	1f
	mul.s	$f2, $f12, $f12
	jal	asincosf2
1:	li.s	$f8, FACS_pio2
	sub.s	$f0, $f8, $f0
	j	t9
acosf2:
	jal	asincosf1
	/* nop */
	bltz	t0, 1f
	/* nop */
	neg.s	$f0
	j	t9
1:	li.s	$f8, FACS_pi
	add.s	$f0, $f8
	j	t9
	.end acosf

/*******************************************************************************
*
* asinf - floating point arc-sine
*

* float asinf (fltParam)
*     float fltParam;

*/
	.globl asinf
	.ent   asinf
asinf:
	.frame	sp, 0, ra
	li.s	$f8, half
	abs.s	$f14, $f12
	c.ole.s	$f14, $f8
	li.s	$f10, FACS_eps
	move	t9, ra
	bc1f	asinf2
	c.lt.s	$f14, $f10
	mul.s	$f2, $f12, $f12
	mov.s	$f0, $f12
	bc1f	asincosf2
	j	ra
	/* nop */

asinf2:
	jal	asincosf1
	/* nop */
	li.s	$f8, FACS_pio2
	add.s	$f0, $f8
	bltz	t0, 1f
	j	t9
	/* nop */
1:	neg.s	$f0
	j	t9


asincosf1:
	li.s	$f10, one
	sw	ra, 0(sp)
	c.ole.s	$f14, $f10
	s.s	$f12, 4(sp)
	sub.s	$f0, $f10, $f14
	bc1f	asinf_error
	mul.s	$f12, $f0, $f8
	s.s	$f12, 8(sp)
	subu	sp, 16
	jal	sqrtf
	addu	sp, 16
	l.s	$f2, 8(sp)
	lw	ra, 0(sp)
	lw	t0, 4(sp)
	add.s	$f0, $f0
	neg.s	$f0
	/* fall through */
asincosf2:
	li.s	$f8, FACS_p2
	li.s	$f10, FACS_q1
	mul.s	$f4, $f2, $f8
	add.s	$f6, $f2, $f10
	li.s	$f8, FACS_p1
	mul.s	$f6, $f2
	add.s	$f4, $f8
	li.s	$f10, FACS_q0
	mul.s	$f4, $f2
	add.s	$f6, $f10

	div.s	$f4, $f6
	mul.s	$f4, $f0
	add.s	$f0, $f4
	j	ra

asinf_error:	/* |x| > 1 */
	c.un.s	$f12, $f12	/* - if x = NaN, return x */
	li.s	$f0, 0.0	/* - else generate a NaN */
	bc1t	1f
	div.s	$f0, $f0
	j	t9
1:	mov.s	$f0, $f12
	j	t9
	.end asinf

/*******************************************************************************
*
* atan2f - single precision floating point arctangent
*

* float atan2f (fltParam)
*     float fltParam;

*/
	.globl atan2f
	.ent   atan2f
atan2f:
	.frame sp, 0, t3
	abs.s	$f0, $f12
	abs.s	$f2, $f14
	c.le.s	$f0, $f2
	mfc1	t0, $f12		/* save signs of both operands */
	mfc1	t1, $f14		/* ... */