*
*   $VER: sinh.s 33.1 (22.1.97)
*
*   hyperbolic sine
*
*   Version history:
*
*   33.1    22.1.97 (c) Motorola
*
*           - snipped from M68060SP sources
*

    machine 68040
    fpu     1

    XDEF    _sinh
    XDEF    @sinh

    XREF    @exp
    XREF    @expm1

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

T1      dc.l            $40C62D38,$D3D64634     ; 16381 LOG2 LEAD
T2      dc.l            $3D6F90AE,$B1E75CC7     ; 16381 LOG2 TRAIL

_sinh
        fmove.d         (4,sp),fp0
@sinh
        fmove.x         fp0,-(sp)
        move.l          (sp),d1
        move.w          (4,sp),d1
        move.l          d1,a1                   ; save (compacted) operand
        lea             (12,sp),sp
        and.l           #$7FFFFFFF,d1
        cmp.l           #$400CB167,d1
        bgt.b           .SINHBIG

;--THIS IS THE USUAL CASE, |X| < 16380 LOG2
;--Y = |X|, Z = EXPM1(Y), SINH(X) = SIGN(X)*(1/2)*( Z + Z/(1+Z) )

        fabs.x          fp0

        move.l          a1,-(sp)                ; {a1}
        jsr             @expm1                  ; FP0 IS Z = EXPM1(Y)
        move.l          (sp)+,a1                ; {a1}

        fmove.x         fp0,fp1
        fadd.s          #$3F800000,fp1          ; 1+Z
        fmove.x         fp0,-(sp)
        fdiv.x          fp1,fp0                 ; Z/(1+Z)
        move.l          a1,d1
        and.l           #$80000000,d1
        or.l            #$3F000000,d1
        fadd.x          (sp)+,fp0
        move.l          d1,-(sp)
        fmul.s          (sp)+,fp0               ; last fp inst - possible exceptions set
        rts

.SINHBIG
        cmp.l           #$400CB2B3,d1
        bgt             .t_ovfl

        fabs.x          fp0
        fsub.d          (T1,pc),fp0             ; (|X|-16381LOG2_LEAD)
        clr.l           -(sp)
        move.l          #$80000000,-(sp)
        move.l          a1,d1
        and.l           #$80000000,d1
        or.l            #$7FFB0000,d1
        move.l          d1,-(sp)                ; EXTENDED FMT
        fsub.d          (T2,pc),fp0             ; |X| - 16381 LOG2, ACCURATE
        jsr             @exp
        fmul.x          (sp)+,fp0               ; possible exception
.t_ovfl
        rts