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

    machine 68040
    fpu     1

    XDEF    _acos
    XDEF    @acos

    XREF    @atan

*************************************************************************
* acos():  computes the inverse cosine of a normalized input            *
*                                                                       *
* INPUT *************************************************************** *
*       fp0 = extended precision input                                  *
*                                                                       *
* OUTPUT ************************************************************** *
*       fp0 = arccos(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 *********************************************************** *
*                                                                       *
*       ACOS                                                            *
*       1. If |X| >= 1, go to 3.                                        *
*                                                                       *
*       2. (|X| < 1) Calculate acos(X) by                               *
*               z := (1-X) / (1+X)                                      *
*               acos(X) = 2 * atan( sqrt(z) ).                          *
*               Exit.                                                   *
*                                                                       *
*       3. If |X| > 1, go to 5.                                         *
*                                                                       *
*       4. (|X| = 1) If X > 0, return 0. Otherwise, return Pi. Exit.    *
*                                                                       *
*       5. (|X| > 1) Generate an invalid operation by 0 * infinity.     *
*               Exit.                                                   *
*                                                                       *
*************************************************************************

PI      dc.l            $40000000,$C90FDAA2,$2168C235,$00000000

_acos
        fmove.d         (4,sp),fp0
@acos
        fmove.x         fp0,-(sp)
        move.l          (sp),d1                 ; pack exp w/ upper 16 fraction
        move.w          (4,sp),d1
        and.l           #$7FFFFFFF,d1
        cmp.l           #$3FFF8000,d1
        bge.b           .ACOSBIG

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

.ACOSMAIN
        lea             (12,sp),sp
        fmove.s         #$3F800000,fp1
        fadd.x          fp0,fp1                 ; 1+X
        fneg.x          fp0                     ; -X
        fadd.s          #$3F800000,fp0          ; 1-X
        fdiv.x          fp1,fp0                 ; (1-X)/(1+X)
        fsqrt.x         fp0                     ; SQRT((1-X)/(1+X))
        jsr             @atan                   ; ATAN(SQRT([1-X]/[1+X]))

        fadd.x          fp0,fp0                 ; 2 * ATAN( STUFF )
        rts

.ACOSBIG
        fabs.x          fp0
        fcmp.s          #$3f800000,fp0
        fbgt            .operr

;--|X| = 1, ACOS(X) = 0 OR PI
        tst.b           (sp)                    ; is X positive or negative?
        bpl.b           .ACOSP1

;--X = -1
;Returns PI and inexact exception
.ACOSM1
        fmove.x         (PI,pc),fp0             ; load PI
.operr  lea             (12,sp),sp
        rts
.ACOSP1 fmove.s         #0,fp0
        lea             (12,sp),sp
        rts