*
*   $VER: cosh.s 33.1 (23.1.97)
*
*   hyperbolic cosine
*
*   Version history:
*
*   33.1    23.1.97 (c) Motorola
*
*           - snipped from M68060SP sources
*

    machine 68040
    fpu     1

    XDEF    _cosh
    XDEF    @cosh

    XREF    @exp

*************************************************************************
* cosh():  computes the hyperbolic cosine of a normalized input         *
*                                                                       *
* INPUT *************************************************************** *
*       fp0 = extended precision input                                  *
*                                                                       *
* OUTPUT ************************************************************** *
*       fp0 = cosh(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 *********************************************************** *
*                                                                       *
*       COSH                                                            *
*       1. If |X| > 16380 log2, go to 3.                                *
*                                                                       *
*       2. (|X| <= 16380 log2) Cosh(X) is obtained by the formulae      *
*               y = |X|, z = exp(Y), and                                *
*               cosh(X) = (1/2)*( z + 1/z ).                            *
*               Exit.                                                   *
*                                                                       *
*       3. (|X| > 16380 log2). If |X| > 16480 log2, go to 5.            *
*                                                                       *
*       4. (16380 log2 < |X| <= 16480 log2)                             *
*               cosh(X) = sign(X) * exp(|X|)/2.                         *
*               However, invoking exp(|X|) may cause premature          *
*               overflow. Thus, we calculate sinh(X) as follows:        *
*               Y       := |X|                                          *
*               Fact    :=      2**(16380)                              *
*               Y'      := Y - 16381 log2                               *
*               cosh(X) := Fact * exp(Y').                              *
*               Exit.                                                   *
*                                                                       *
*       5. (|X| > 16480 log2) sinh(X) must overflow. Return             *
*               Huge*Huge to generate overflow and an infinity with     *
*               the appropriate sign. Huge is the largest finite number *
*               in extended format. Exit.                               *
*                                                                       *
*************************************************************************

TWO16380
        dc.l            $7FFB0000,$80000000,$00000000,$00000000
T1      dc.l            $40C62D38,$D3D64634     ; 16381 LOG2 LEAD
T2      dc.l            $3D6F90AE,$B1E75CC7     ; 16381 LOG2 TRAIL

_cosh
        fmove.d         (4,sp),fp0
@cosh
        fmove.x         fp0,-(sp)

        move.l          (sp),d1
        move.w          (4,sp),d1
        and.l           #$7FFFFFFF,d1
        lea             (12,sp),sp
        cmp.l           #$400CB167,d1
        bgt.b           .COSHBIG

;--THIS IS THE USUAL CASE, |X| < 16380 LOG2
;--COSH(X) = (1/2) * ( EXP(X) + 1/EXP(X) )

        fabs.x          fp0                     ; |X|

        jsr             @exp                    ; FP0 IS EXP(|X|)
        fmul.s          #$3F000000,fp0          ; (1/2)EXP(|X|)

        fmove.s         #$3E800000,fp1          ; (1/4)
        fdiv.x          fp0,fp1                 ; 1/(2 EXP(|X|))

        fadd.x          fp1,fp0
        rts

.COSHBIG
        cmp.l           #$400CB2B3,d1
        bgt.b           .COSHHUGE

        fabs.x          fp0
        fsub.d          (T1,pc),fp0             ; (|X|-16381LOG2_LEAD)
        fsub.d          (T2,pc),fp0             ; |X| - 16381 LOG2, ACCURATE

        jsr             @exp
        fmul.x          (TWO16380,pc),fp0
.COSHHUGE
        rts