arch/m68k-all/m680x0/fpsp/scosh.sa

   1 *       $NetBSD: scosh.sa,v 1.2 1994/10/26 07:49:39 cgd Exp $
   2
   3 *       MOTOROLA MICROPROCESSOR & MEMORY TECHNOLOGY GROUP
   4 *       M68000 Hi-Performance Microprocessor Division
   5 *       M68040 Software Package
   6 *
   7 *       M68040 Software Package Copyright (c) 1993, 1994 Motorola Inc.
   8 *       All rights reserved.
   9 *
  10 *       THE SOFTWARE is provided on an "AS IS" basis and without warranty.
  11 *       To the maximum extent permitted by applicable law,
  12 *       MOTOROLA DISCLAIMS ALL WARRANTIES WHETHER EXPRESS OR IMPLIED,
  13 *       INCLUDING IMPLIED WARRANTIES OF MERCHANTABILITY OR FITNESS FOR A
  14 *       PARTICULAR PURPOSE and any warranty against infringement with
  15 *       regard to the SOFTWARE (INCLUDING ANY MODIFIED VERSIONS THEREOF)
  16 *       and any accompanying written materials.
  17 *
  18 *       To the maximum extent permitted by applicable law,
  19 *       IN NO EVENT SHALL MOTOROLA BE LIABLE FOR ANY DAMAGES WHATSOEVER
  20 *       (INCLUDING WITHOUT LIMITATION, DAMAGES FOR LOSS OF BUSINESS
  21 *       PROFITS, BUSINESS INTERRUPTION, LOSS OF BUSINESS INFORMATION, OR
  22 *       OTHER PECUNIARY LOSS) ARISING OF THE USE OR INABILITY TO USE THE
  23 *       SOFTWARE.  Motorola assumes no responsibility for the maintenance
  24 *       and support of the SOFTWARE.
  25 *
  26 *       You are hereby granted a copyright license to use, modify, and
  27 *       distribute the SOFTWARE so long as this entire notice is retained
  28 *       without alteration in any modified and/or redistributed versions,
  29 *       and that such modified versions are clearly identified as such.
  30 *       No licenses are granted by implication, estoppel or otherwise
  31 *       under any patents or trademarks of Motorola, Inc.
  32
  33 *
  34 *       scosh.sa 3.1 12/10/90
  35 *
  36 *       The entry point sCosh computes the hyperbolic cosine of
  37 *       an input argument; sCoshd does the same except for denormalized
  38 *       input.
  39 *
  40 *       Input: Double-extended number X in location pointed to
  41 *               by address register a0.
  42 *
  43 *       Output: The value cosh(X) returned in floating-point register Fp0.
  44 *
  45 *       Accuracy and Monotonicity: The returned result is within 3 ulps in
  46 *               64 significant bit, i.e. within 0.5001 ulp to 53 bits if the
  47 *               result is subsequently rounded to double precision. The
  48 *               result is provably monotonic in double precision.
  49 *
  50 *       Speed: The program sCOSH takes approximately 250 cycles.
  51 *
  52 *       Algorithm:
  53 *
  54 *       COSH
  55 *       1. If |X| > 16380 log2, go to 3.
  56 *
  57 *       2. (|X| <= 16380 log2) Cosh(X) is obtained by the formulae
  58 *               y = |X|, z = exp(Y), and
  59 *               cosh(X) = (1/2)*( z + 1/z ).
  60 *               Exit.
  61 *
  62 *       3. (|X| > 16380 log2). If |X| > 16480 log2, go to 5.
  63 *
  64 *       4. (16380 log2 < |X| <= 16480 log2)
  65 *               cosh(X) = sign(X) * exp(|X|)/2.
  66 *               However, invoking exp(|X|) may cause premature overflow.
  67 *               Thus, we calculate sinh(X) as follows:
  68 *               Y       := |X|
  69 *               Fact    :=      2**(16380)
  70 *               Y'      := Y - 16381 log2
  71 *               cosh(X) := Fact * exp(Y').
  72 *               Exit.
  73 *
  74 *       5. (|X| > 16480 log2) sinh(X) must overflow. Return
  75 *               Huge*Huge to generate overflow and an infinity with
  76 *               the appropriate sign. Huge is the largest finite number in
  77 *               extended format. Exit.
  78 *
  79
  80 SCOSH   IDNT    2,1 Motorola 040 Floating Point Software Package
  81
  82         section 8
  83
  84         xref    t_ovfl
  85         xref    t_frcinx
  86         xref    setox
  87
  88 T1      DC.L $40C62D38,$D3D64634 ... 16381 LOG2 LEAD
  89 T2      DC.L $3D6F90AE,$B1E75CC7 ... 16381 LOG2 TRAIL
  90
  91 TWO16380 DC.L $7FFB0000,$80000000,$00000000,$00000000
  92
  93         xdef    scoshd
  94 scoshd:
  95 *--COSH(X) = 1 FOR DENORMALIZED X
  96
  97         FMOVE.S         #:3F800000,FP0
  98
  99         FMOVE.L         d1,FPCR
 100         FADD.S          #:00800000,FP0
 101         bra             t_frcinx
 102
 103         xdef    scosh
 104 scosh:
 105         FMOVE.X         (a0),FP0        ...LOAD INPUT
 106
 107         move.l          (a0),d0
 108         move.w          4(a0),d0
 109         ANDI.L          #$7FFFFFFF,d0
 110         CMPI.L          #$400CB167,d0
 111         BGT.B           COSHBIG
 112
 113 *--THIS IS THE USUAL CASE, |X| < 16380 LOG2
 114 *--COSH(X) = (1/2) * ( EXP(X) + 1/EXP(X) )
 115
 116         FABS.X          FP0             ...|X|
 117
 118         move.l          d1,-(sp)
 119         clr.l           d1
 120         fmovem.x        fp0,(a0)        ;pass parameter to setox
 121         bsr             setox           ...FP0 IS EXP(|X|)
 122         FMUL.S          #:3F000000,FP0  ...(1/2)EXP(|X|)
 123         move.l          (sp)+,d1
 124
 125         FMOVE.S         #:3E800000,FP1  ...(1/4)
 126         FDIV.X          FP0,FP1         ...1/(2 EXP(|X|))
 127
 128         FMOVE.L         d1,FPCR
 129         FADD.X          fp1,FP0
 130
 131         bra             t_frcinx
 132
 133 COSHBIG:
 134         CMPI.L          #$400CB2B3,d0
 135         BGT.B           COSHHUGE
 136
 137         FABS.X          FP0
 138         FSUB.D          T1(pc),FP0              ...(|X|-16381LOG2_LEAD)
 139         FSUB.D          T2(pc),FP0              ...|X| - 16381 LOG2, ACCURATE
 140
 141         move.l          d1,-(sp)
 142         clr.l           d1
 143         fmovem.x        fp0,(a0)
 144         bsr             setox
 145         fmove.l         (sp)+,fpcr
 146
 147         FMUL.X          TWO16380(pc),FP0
 148         bra             t_frcinx
 149
 150 COSHHUGE:
 151         fmove.l         #0,fpsr         ;clr N bit if set by source
 152         bclr.b          #7,(a0)         ;always return positive value
 153         fmovem.x        (a0),fp0
 154         bra             t_ovfl
 155
 156         end