source: rtems/c/src/lib/libcpu/m68k/m68040/fpsp/stanh.s @ 010e7582

//
//      stanh.sa 3.1 12/10/90
//
//      The entry point sTanh computes the hyperbolic tangent of
//      an input argument; sTanhd does the same except for denormalized
//      input.
//
//      Input: Double-extended number X in location pointed to
//              by address register a0.
//
//      Output: The value tanh(X) returned in floating-point register Fp0.
//
//      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.
//
//      Speed: The program stanh takes approximately 270 cycles.
//
//      Algorithm:
//
//      TANH
//      1. If |X| >= (5/2) log2 or |X| <= 2**(-40), go to 3.
//
//      2. (2**(-40) < |X| < (5/2) log2) Calculate tanh(X) by
//              sgn := sign(X), y := 2|X|, z := expm1(Y), and
//              tanh(X) = sgn*( z/(2+z) ).
//              Exit.
//
//      3. (|X| <= 2**(-40) or |X| >= (5/2) log2). If |X| < 1,
//              go to 7.
//
//      4. (|X| >= (5/2) log2) If |X| >= 50 log2, go to 6.
//
//      5. ((5/2) log2 <= |X| < 50 log2) Calculate tanh(X) by
//              sgn := sign(X), y := 2|X|, z := exp(Y),
//              tanh(X) = sgn - [ sgn*2/(1+z) ].
//              Exit.
//
//      6. (|X| >= 50 log2) Tanh(X) = +-1 (round to nearest). Thus, we
//              calculate Tanh(X) by
//              sgn := sign(X), Tiny := 2**(-126),
//              tanh(X) := sgn - sgn*Tiny.
//              Exit.
//
//      7. (|X| < 2**(-40)). Tanh(X) = X.       Exit.
//

//              Copyright (C) Motorola, Inc. 1990
//                      All Rights Reserved
//
//      THIS IS UNPUBLISHED PROPRIETARY SOURCE CODE OF MOTOROLA 
//      The copyright notice above does not evidence any  
//      actual or intended publication of such source code.

//STANH idnt    2,1 | Motorola 040 Floating Point Software Package

        |section        8
        
        .include "fpsp.defs"

        .set    X,FP_SCR5
        .set    XDCARE,X+2
        .set    XFRAC,X+4

        .set    SGN,L_SCR3

        .set    V,FP_SCR6

BOUNDS1:        .long 0x3FD78000,0x3FFFDDCE // ... 2^(-40), (5/2)LOG2

        |xref   t_frcinx
        |xref   t_extdnrm
        |xref   setox
        |xref   setoxm1

        .global stanhd
stanhd:
//--TANH(X) = X FOR DENORMALIZED X

        bra             t_extdnrm

        .global stanh
stanh:
        fmovex          (%a0),%fp0      // ...LOAD INPUT

        fmovex          %fp0,X(%a6)
        movel           (%a0),%d0
        movew           4(%a0),%d0
        movel           %d0,X(%a6)
        andl            #0x7FFFFFFF,%d0
        cmp2l           BOUNDS1(%pc),%d0        // ...2**(-40) < |X| < (5/2)LOG2 ?
        bcss            TANHBORS

//--THIS IS THE USUAL CASE
//--Y = 2|X|, Z = EXPM1(Y), TANH(X) = SIGN(X) * Z / (Z+2).

        movel           X(%a6),%d0
        movel           %d0,SGN(%a6)
        andl            #0x7FFF0000,%d0
        addl            #0x00010000,%d0 // ...EXPONENT OF 2|X|
        movel           %d0,X(%a6)
        andl            #0x80000000,SGN(%a6)
        fmovex          X(%a6),%fp0             // ...FP0 IS Y = 2|X|

        movel           %d1,-(%a7)
        clrl            %d1
        fmovemx %fp0-%fp0,(%a0)
        bsr             setoxm1         // ...FP0 IS Z = EXPM1(Y)
        movel           (%a7)+,%d1

        fmovex          %fp0,%fp1
        fadds           #0x40000000,%fp1        // ...Z+2
        movel           SGN(%a6),%d0
        fmovex          %fp1,V(%a6)
        eorl            %d0,V(%a6)

        fmovel          %d1,%FPCR               //restore users exceptions
        fdivx           V(%a6),%fp0
        bra             t_frcinx

TANHBORS:
        cmpl            #0x3FFF8000,%d0
        blt             TANHSM

        cmpl            #0x40048AA1,%d0
        bgt             TANHHUGE

//-- (5/2) LOG2 < |X| < 50 LOG2,
//--TANH(X) = 1 - (2/[EXP(2X)+1]). LET Y = 2|X|, SGN = SIGN(X),
//--TANH(X) = SGN -     SGN*2/[EXP(Y)+1].

        movel           X(%a6),%d0
        movel           %d0,SGN(%a6)
        andl            #0x7FFF0000,%d0
        addl            #0x00010000,%d0 // ...EXPO OF 2|X|
        movel           %d0,X(%a6)              // ...Y = 2|X|
        andl            #0x80000000,SGN(%a6)
        movel           SGN(%a6),%d0
        fmovex          X(%a6),%fp0             // ...Y = 2|X|

        movel           %d1,-(%a7)
        clrl            %d1
        fmovemx %fp0-%fp0,(%a0)
        bsr             setox           // ...FP0 IS EXP(Y)
        movel           (%a7)+,%d1
        movel           SGN(%a6),%d0
        fadds           #0x3F800000,%fp0        // ...EXP(Y)+1

        eorl            #0xC0000000,%d0 // ...-SIGN(X)*2
        fmoves          %d0,%fp1                // ...-SIGN(X)*2 IN SGL FMT
        fdivx           %fp0,%fp1               // ...-SIGN(X)2 / [EXP(Y)+1 ]

        movel           SGN(%a6),%d0
        orl             #0x3F800000,%d0 // ...SGN
        fmoves          %d0,%fp0                // ...SGN IN SGL FMT

        fmovel          %d1,%FPCR               //restore users exceptions
        faddx           %fp1,%fp0

        bra             t_frcinx

TANHSM:
        movew           #0x0000,XDCARE(%a6)

        fmovel          %d1,%FPCR               //restore users exceptions
        fmovex          X(%a6),%fp0             //last inst - possible exception set

        bra             t_frcinx

TANHHUGE:
//---RETURN SGN(X) - SGN(X)EPS
        movel           X(%a6),%d0
        andl            #0x80000000,%d0
        orl             #0x3F800000,%d0
        fmoves          %d0,%fp0
        andl            #0x80000000,%d0
        eorl            #0x80800000,%d0 // ...-SIGN(X)*EPS

        fmovel          %d1,%FPCR               //restore users exceptions
        fadds           %d0,%fp0

        bra             t_frcinx

        |end