*
*   $VER: asin.s 33.1 (22.1.97)
*
*   Calculates the arcsine of the source
*
*   Version history:
*
*   33.1    22.1.97 (c) Motorola
*
*           - snipped from M68060SP sources
*

    machine 68040
    fpu     1

    XDEF    _asin
    XDEF    @asin

    XREF    @atan

*************************************************************************
* asin():  computes the inverse sine of a normalized input              *
*                                                                       *
* INPUT *************************************************************** *
*       fp0 = extended precision input                                  *
*                                                                       *
* OUTPUT ************************************************************** *
*       fp0 = arcsin(X)                                                 *
*                                                                       *
* ACCURACY and MONOTONICITY ******************************************* *
*       The returned result is within 3 ulps in 64 significant bit,     *
*       i.e. within 0.5001 ulp to 53 bits if the result is subsequently *
*       rounded to double precision. The result is provably monotonic   *
*       in double precision.                                            *
*                                                                       *
* ALGORITHM *********************************************************** *
*                                                                       *
*       ASIN                                                            *
*       1. If |X| >= 1, go to 3.                                        *
*                                                                       *
*       2. (|X| < 1) Calculate asin(X) by                               *
*               z := sqrt( [1-X][1+X] )                                 *
*               asin(X) = atan( x / z ).                                *
*               Exit.                                                   *
*                                                                       *
*       3. If |X| > 1, go to 5.                                         *
*                                                                       *
*       4. (|X| = 1) sgn := sign(X), return asin(X) := sgn * Pi/2. Exit.*
*                                                                       *
*       5. (|X| > 1) Generate an invalid operation by 0 * infinity.     *
*               Exit.                                                   *
*                                                                       *
*************************************************************************

PIBY2   dc.l            $3FFF0000,$C90FDAA2,$2168C235,$00000000

_asin
        fmove.d         (4,sp),fp0
@asin
        fmove.x         fp0,-(sp)

        move.l          (sp),d1
        move.w          (4,sp),d1
        and.l           #$7FFFFFFF,d1
        cmp.l           #$3FFF8000,d1
        bge.b           .ASINBIG

; This catch is added here for the '060 QSP. Originally, the call to
; satan() would handle this case by causing the exception which would
; not be caught until gen_except(). Now, with the exceptions being
; detected inside of satan(), the exception would have been handled there
; instead of inside sasin() as expected.
        cmp.l           #$3FD78000,d1
        blt.w           .ASINTINY

;--THIS IS THE USUAL CASE, |X| < 1
;--ASIN(X) = ATAN( X / SQRT( (1-X)(1+X) ) )

.ASINMAIN
        lea             (12,sp),sp
        fmove.s         #$3F800000,fp1
        fsub.x          fp0,fp1                 ; 1-X
        fmovem.x        fp2,-(sp)               ;  {fp2}
        fmove.s         #$3F800000,fp2
        fadd.x          fp0,fp2                 ; 1+X
        fmul.x          fp2,fp1                 ; (1+X)(1-X)
        fmovem.x        (sp)+,fp2               ;  {fp2}
        fsqrt.x         fp1                     ; SQRT([1-X][1+X])
        fdiv.x          fp1,fp0                 ; X/SQRT([1-X][1+X])
        jsr             @atan
        rts

.ASINBIG
        fabs.x          fp0                     ; |X|
        fcmp.s          #$3F800000,fp0
        fbgt            .operr                  ; cause an operr exception

;--|X| = 1, ASIN(X) = +- PI/2.
.ASINONE
        fmove.x         (PIBY2,pc),fp0
        move.l          (sp),d1
        and.l           #$80000000,d1           ; SIGN BIT OF X
        lea             (12,sp),sp
        or.l            #$3F800000,d1           ; +-1 IN SGL FORMAT
        move.l          d1,-(sp)                ; push SIGN(X) IN SGL-FMT
        fmul.s          (sp)+,fp0
        rts

;--|X| < 2^(-40), ATAN(X) = X
.ASINTINY
.operr
        lea             (12,sp),sp
        rts