mathhardalib.s

vxworks source code, used for develop vxworks system.

	/* is $f12 or $f14 = +-Inf or Nan ? */
        sll     t2, t0, 1
        sll     t3, t1, 1
        srl     t2, t2, 1
        srl     t3, t3, 1
        li      t4, 255
        sll     t4, 23
        blt     t2, t4, 5f 
	bgt	t2, t4, 3f
        blt     t3, t4, 5f 
	bgt	t3, t4, 4f
        /* x and y are +/- Inf */
        li.s    $f0, TN_pio4                    /* f0 = pi/4 */
        bgez    t0, 1f
        neg.s   $f0                             /* f0 = -pi/4 */
1:      bgez    t1, 2f
        li.s    $f2, 3.0
        mul.s   $f0, $f2, $f0
2:      nop
	j	ra

3:      /* atan2f (Nan, ?) = Nan */
	mov.s	$f0, $f12
	j	ra

4:      /* atan2f (?, Nan) = Nan */
	mov.s	$f0, $f14
	j	ra

5:	move	t3, ra
	bc1t	atan2f1
	div.s	$f0, $f2, $f0
	li.s	$f2, FAT_mpio2
	b	atan2f2
atan2f1:
	beq	t1, 0, atan2fz
	div.s	$f0, $f2
	li.s	$f2, 0.0
atan2f2:
	li.s	$f10, FAT_twomr3
	jal	atan1
	bge	t1, 0, atan2f3
	li.s	$f2, FAT_pi
	sub.s	$f0, $f2, $f0
atan2f3:
	bge	t0, 0, atan2f4
	neg.s	$f0
	j	t3	
atan2f4:
	j	t3

atan2fz:
	/* break 0 */
	j	ra
	.end atan2f

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

* float atan (fltParam)
*     float fltParam;

*/
	.globl atanf
	.ent   atanf
atanf:
	.frame sp, 0, t3
	mfc1	t0, $f12
	move	t3, ra
	/* is $f12 = Nan ? */
        sll     t2, t0, 1
        srl     t2, t2, 1
        li      t4, 255
        sll     t4, 23
        ble     t2, t4, 1f
	/* atanf (Nan) = Nan */
	mov.s	$f0, $f12
	j	t3

1:	abs.s	$f0, $f12
	bge	t0, 0, atanf0
	jal	atanf0
	neg.s	$f0
	j	t3

atanf0:	li.s	$f14, one
	li.s	$f10, FAT_twomr3
	c.le.s	$f0, $f14
	li.s	$f2, 0.0
	bc1t	atan1
	div.s	$f0, $f14, $f0
	li.s	$f2, FAT_mpio2
atan1:	c.le.s	$f0, $f10
	li.s	$f14, FAT_sqrt3m1
	bc1t	latan2		/* branch to local atan2 */
	li.s	$f10, FAT_sqrt3
	mul.s	$f6, $f0, $f14
	add.s	$f4, $f0, $f10
	li.s	$f14, one
	sub.s	$f6, $f14
	add.s	$f0, $f6
	li.s	$f14, FAT_pio6
	div.s	$f0, $f4
	add.s	$f2, $f14
latan2:	mul.s	$f6, $f0, $f0
	li.s	$f14, FAT_p1
	li.s	$f10, FAT_p0
	mul.s	$f4, $f6, $f14
	add.s	$f4, $f10
	li.s	$f14, FAT_q0
	mul.s	$f12, $f4, $f6
	add.s	$f4, $f6, $f14
	div.s	$f4, $f12, $f4
	mul.s	$f4, $f0
	add.s	$f0, $f4
	mfc1	t4, $f2
	add.s	$f0, $f2
	bge	t4, 0, atan4
	neg.s	$f0
	j	ra
atan4:	j	ra
	.end atanf

/*******************************************************************************
*
* fabsf - single precision absolute value
*

* float fabs (fltParam)
*     float fltParam;

*/
        .globl fabsf
        .ent   fabsf
fabsf:
        abs.s   $f0, $f12
        j       ra
        .end fabsf

/*******************************************************************************
*
* fcabs - floating point single precision complex number
*
* RETURN THE ABSOLUTE VALUE OF THE SINGLE PRECISION COMPLEX NUMBER

* float fabs (z)
* struct { float r, i;} z;
* {
* 	return hypotf(z.r,z.i);
* }

*/
	.globl fcabs
	.ent   fcabs
fcabs:
	.frame	sp, 0, ra
	mtc1	$4, $f12
	mtc1	$5, $f14
	j	hypotf
	.end fcabs

/*******************************************************************************
*
* coshf - single precision floating point hyperbolic cosine
*

* float cosh (fltParam)
*     float fltParam;

*/
	.globl coshf
	.ent   coshf
coshf:
	.frame	sp, 0, ra
	li.s	$f10, 88.7228317
	abs.s	$f12
	c.ole.s	$f12, $f10
	SW	ra, 0(sp)
	bc1f	coshf_error
	subu	sp, 16
	jal	expf
	addu	sp, 16
	li.s	$f10, 0.5
	div.s	$f2, $f10, $f0
	LW	ra, 0(sp)
	mul.s	$f0, $f10
	add.s	$f0, $f2
	j	ra

coshf1:
#if	FALSE
	li.s	$f8, 0.69316101074218750000e+0
	li.s	$f6, expfmax
	sub.s	$f12, $f8
	c.lt.s	$f6, $f12
	SW	ra, 0(sp)
	bc1t	coshf_error
	subu	sp, 16
	jal	expf
	addu	sp, 16
	li.s	$f8, 1.0
	li.s	$f6, 0.13830277879601902638e-4
	sub.s	$f2, $f0, $f8
	LW	ra, 0(sp)
	mul.s	$f2, $f6
	add.s	$f0, $f2
	j	ra
#endif	/* FALSE */
coshf_error:
	/* raise Overflow and return +Infinity */
	mfc1	t0, $f12
	srl	t0, 23
	beq	t0, 255, 1f
	li.s	$f0, 2e38
	add.s	$f0, $f0
	j	ra
1:	mov.s	$f0, $f12
	j	ra
	.end coshf

/*******************************************************************************
*
* expf - floating point single precision exponetial function
*

* float expf (fltParam)
*     float fltParam;

*/
        .globl expf
        .ent   expf
expf:

        .frame sp, 0, ra
        /* argument in f12 */
        .set noreorder
        li.s    $f10, 88.7228317		/* expfmax */
        cfc1    t4, $31                         /* read fp control/status */
        c.ole.s $f12, $f10                      /* check if special */
        li.s    $f14, -103.2789299019278	/* expfmin */
        bc1f    90f                             /* - if NaN, +Infinity, */
                                                /* or greater than expmax */
        c.lt.s  $f12, $f14                      /* check for exp(x) = 0 */
        li      t1, -4
        bc1t    80f                             /* - if less than expmin */
        and     t1, t4                          /* rounding mode = nearest */
        .set reorder
        ctc1    t1, $31                         /* write fp control/status */
        /* argument reduction */
        li.s    $f10, 46.166241308446828
        mul.s   $f2, $f12, $f10
        li.s    $f14, 2.1660849390173098e-2
        li.s    $f10, -2.325192846878874e-12
        cvt.w.s $f2, $f2
        mfc1    t0, $f2
        and     t1, t0, 31                      /* region  32<t1<=0 */
        sra     t2, t0, 5                       /* scale t2=t0/32 */
        cvt.s.w $f2, $f2
        mul.s   $f4, $f2, $f14
        mul.s   $f6, $f2, $f10
        sub.s   $f4, $f12, $f4

        add.s   $f2, $f4, $f6
        li.s    $f10, 1.3888949086377719e-3
        li.s    $f14, 8.3333679843421958e-3
        mul.s   $f8, $f2, $f10
        add.s   $f0, $f8, $f14
        li.s    $f10, 4.1666666666226079e-2
        li.s    $f14, 1.6666666666526087e-1
        mul.s   $f0, $f2
        add.s   $f0, $f10
        mul.s   $f0, $f2
        add.s   $f0, $f14
        li.s    $f10, 0.5
        mul.s   $f0, $f2
        add.s   $f0, $f10
        mul.s   $f0, $f2
        mul.s   $f0, $f2
        add.s   $f0, $f6
        add.s   $f0, $f4
        sll     t1, 3
        la      t1, EXPF_DAT(t1)
        l.s     $f10, 0(t1)
        l.s     $f14, 4(t1)
        add.s   $f8, $f10, $f14
        mul.s   $f0, $f8
        add.s   $f0, $f14
        add.s   $f0, $f10
        beq     t2, 0, 60f              /* early out for 0 scale */
        mfc1    t0, $f0                 /* get result high word */
        ctc1    t4, $31                 /* restore control/status */
        sll     t1, t0, 1               /* extract exponent */
        srl     t1, 23+1                /* ... */
        addu    t1, t2                  /* add scale to check for denorm */
        sll     t2, 23
        blez    t1, 70f
        addu    t0, t2                  /* add scale */
        mtc1    t0, $f0                 /* put back in result high word */
        j       ra
60:     /* scale = 0, just restore control/status and return */
        ctc1    t4, $31
        j       ra
70:     /* denorm result */
        addu    t2, 64<<23
        addu    t0, t2
        mtc1    t0, $f1
        li.s    $f2, 5.4210108624275222e-20
        mul.s   $f0, $f2
        j       ra

80:     /* argument < expmin */
        li.s    $f0, 0.0		/* should raise underflow */
        j       ra

90:     /* raise Overflow and return +Infinity */
        mfc1    t0, $f12		/* extract argument exponent */
        sll     t0, 1
        srl     t0, 23+1
        beq     t0, 255, 91f		/* - if NaN or Infinity */
        li.s    $f0, 1.701411735e38
        add.s   $f0, $f0		/* raise Overflow */
        j       ra
91:     mov.s   $f0, $f12		/* dont raise any flags */
                                        /* - if argument is exceptional */ 
					/* just return argument */
        j       ra
        .end expf


/*******************************************************************************
*
* truncf - floating point single precision truncation
*

* float truncf (fltParam)
*     float fltParam;

*/
	.globl truncf
	.ent   truncf
truncf:
	.frame	sp, 0, ra
	mfc1	t0, $f12
	srl	t2, t0, 23
	and	t2, 0xFF
	sub	t2, 127
	bge	t2, 0, truncf1
	mtc1	$0, $f0
	j	ra
truncf1:	
	sub	t2, 23
	bge	t2, 0, truncf2
	neg	t2
	srl	t0, t2
	sll	t0, t2
truncf2:
	mtc1	t0, $f0
	j	ra
	.end truncf

/*******************************************************************************
*
* floorf - single precision floating point floor
*

* float floorf (fltParam)
*     float fltParam;

* XXX WARNING: 
*          following routine incorrectly returns 0 when +/- INF is the input
*          this should be corrected if this routine is to be used in the future
*/
	.globl floorf
	.ent   floorf
floorf:
	.frame	sp, 0, t3
	move	t3, ra
	/* is x = +/- Inf ? */
	mfc1    t0, $f12
        sll     t2, t0, 1
        srl     t2, t2, 1
        li      t4, 255
        sll     t4, 23
        bne     t2, t4, 1f
        /* x is +/- Inf */
        li.s    $f0, 0.0
        j       ra

1:	jal	truncf
	sub.s	$f2, $f12, $f0
	mfc1	t0, $f2
	li.s	$f2, 1.0
	sll	t1, t0, 1
	bge	t0, 0, 2f
	beq	t1, 0, 2f
	sub.s	$f0, $f2
2:	j	t3
	.end floorf

/*******************************************************************************
*
* ceilf - single precision floating point ceil
*

* float ceilf (fltParam)
*     float fltParam;

* XXX WARNING: 
*          following routine incorrectly returns 0 when +/- INF is the input
*          this should be corrected if this routine is to be used in the future
*/
	.globl ceilf
	.ent   ceilf
ceilf:
	.frame	sp, 0, t3
	move	t3, ra
	/* is x = +/- Inf ? */
	mfc1    t0, $f12
        sll     t2, t0, 1
        srl     t2, t2, 1
        li      t4, 255
        sll     t4, 23
        bne     t2, t4, 1f
        /* x is +/- Inf */
        li.s    $f0, 0.0
        j       ra

1: 	jal	truncf
	sub.s	$f2, $f12, $f0
	mfc1	t0, $f2
	li.s	$f2, 1.0
	ble	t0, 0, 2f
	add.s	$f0, $f2
2:	j	t3
	.end ceilf


/*******************************************************************************
*
* fmodf - single precision floating point fmod
*

* float fmodf (fltParamX, fltParamY)
*     float fltParamX;
*     float fltParamY;

*/
        .globl fmodf
        .ent   fmodf
fmodf:
        .frame  sp, 0, t3
	move	t3, ra
	abs.s	$f0, $f14
	mfc1    t0, $f0
	beq     t0, zero, retZero		/* |Y| = 0.0 ? */
	lui	t1, INF_HI_SHORT_SINGLE
	bgt	t0, t1, retY			/* |Y| = Nan ? */
	bge     t0, t1, retX			/* |Y| = Inf ? */
	abs.s   $f2, $f12
	mfc1    t2, $f2				/* |X| -> t2 */
	bge	t2, t1, retZero                 /* |X| = Nan or |X| = Inf ? */

	mov.s	$f16, $f12			/* save fltParamX to f16 */
        div.s   $f12, $f14			/* X/Y */
	jal	truncf				/* f0 = truncf(X/Y) */
	nop
	mul.s	$f2, $f0, $f14			
	sub.s	$f0, $f16, $f2	

	/* set result to sign of fltParamX */

	mfc1    t0, $f16
	srl	t0, 31
	sll	t0, 31
	mfc1    t1, $f0
	sll     t1, 1
	srl	t1, 1
	add	t0, t1
	mtc1	t0, $f0
	j	t3
	nop
retZero:
	li.s	$f0, 0.0
	j	t3
retX:	mov.s 	$f0, $f12
	j	t3
retY:	mov.s 	$f0, $f14
	j	t3
	.end fmodf

/*******************************************************************************
*
* hypotf - floating point single precision Euclidean distance
*

* float hypotf (fltX, fltY)
*     float fltX;
*     float fltY;
*
* NOMANUAL
*/
	.globl hypotf
	.ent   hypotf
hypotf:
	.frame	sp, 0, ra
	mul.s	$f12, $f12
	mul.s	$f14, $f14
	add.s	$f12, $f14
	mfc1	t0, $f13
	li	t2, -(1023<<19)+(1023<<20)
	sra	t1, t0, 20
	li	t3, 2047
	beq	t1, 0, 8f
	srl	t0, 1
	beq	t1, t3, 8f
	srl	t1, t0, 15-2
	and	t1, 31<<2
	lw	t1, _sqrttable(t1)
	addu	t0, t2
	subu	t0, t1
	mtc1	t0, $f1
	mtc1	$0, $f0
	/* 8 -> 18 bits */
	li	t2, (1<<20)
	div.s	$f2, $f12, $f0
	/* 18 cycle interlock */
	add.s	$f0, $f2
	/* 1 cycle interlock (2 cycle stall) */
	mfc1	t0, $f1
	add	t1, t2, 6	/* 17 -> 18 bits */
	subu	t0, t1
	mtc1	t0, $f1
	/* nop */
	/* 18 -> 37 bits */
	div.s	$f2, $f12, $f0
	/* 18 cycle interlock */
	add.s	$f0, $f2
	li.s	$f2, half
	mul.s	$f0, $f2
	j	ra

8:	mov.s	$f0, $f12
	abs.s	$f0
	j	ra
	.end hypotf


/*******************************************************************************