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

   1 *       $NetBSD: ssinh.sa,v 1.3 1994/10/26 07:50:05 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 *       ssinh.sa 3.1 12/10/90
  35 *
  36 *       The entry point sSinh computes the hyperbolic sine of
  37 *       an input argument; sSinhd 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 sinh(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 sSINH takes approximately 280 cycles.
  51 *
  52 *       Algorithm:
  53 *
  54 *       SINH
  55 *       1. If |X| > 16380 log2, go to 3.
  56 *
  57 *       2. (|X| <= 16380 log2) Sinh(X) is obtained by the formulae
  58 *               y = |X|, sgn = sign(X), and z = expm1(Y),
  59 *               sinh(X) = sgn*(1/2)*( z + z/(1+z) ).
  60 *          Exit.
  61 *
  62 *       3. If |X| > 16480 log2, go to 5.
  63 *
  64 *       4. (16380 log2 < |X| <= 16480 log2)
  65 *               sinh(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 *             sgn     := sign(X)
  70 *             sgnFact := sgn * 2**(16380)
  71 *             Y'      := Y - 16381 log2
  72 *             sinh(X) := sgnFact * exp(Y').
  73 *          Exit.
  74 *
  75 *       5. (|X| > 16480 log2) sinh(X) must overflow. Return
  76 *          sign(X)*Huge*Huge to generate overflow and an infinity with
  77 *          the appropriate sign. Huge is the largest finite number in
  78 *          extended format. Exit.
  79 *
  80
  81 SSINH   IDNT    2,1 Motorola 040 Floating Point Software Package
  82
  83         section 8
  84
  85 T1      DC.L $40C62D38,$D3D64634 ... 16381 LOG2 LEAD
  86 T2      DC.L $3D6F90AE,$B1E75CC7 ... 16381 LOG2 TRAIL
  87
  88         xref    t_frcinx
  89         xref    t_ovfl
  90         xref    t_extdnrm
  91         xref    setox
  92         xref    setoxm1
  93
  94         xdef    ssinhd
  95 ssinhd:
  96 *--SINH(X) = X FOR DENORMALIZED X
  97
  98         bra     t_extdnrm
  99
 100         xdef    ssinh
 101 ssinh:
 102         FMOVE.x (a0),FP0        ...LOAD INPUT
 103
 104         move.l  (a0),d0
 105         move.w  4(a0),d0
 106         move.l  d0,a1           save a copy of original (compacted) operand
 107         AND.L   #$7FFFFFFF,D0
 108         CMP.L   #$400CB167,D0
 109         BGT.B   SINHBIG
 110
 111 *--THIS IS THE USUAL CASE, |X| < 16380 LOG2
 112 *--Y = |X|, Z = EXPM1(Y), SINH(X) = SIGN(X)*(1/2)*( Z + Z/(1+Z) )
 113
 114         FABS.X  FP0             ...Y = |X|
 115
 116         movem.l a1/d1,-(sp)
 117         fmovem.x fp0,(a0)
 118         clr.l   d1
 119         bsr     setoxm1         ...FP0 IS Z = EXPM1(Y)
 120         fmove.l #0,fpcr
 121         movem.l (sp)+,a1/d1
 122
 123         FMOVE.X FP0,FP1
 124         FADD.S  #:3F800000,FP1  ...1+Z
 125         FMOVE.X FP0,-(sp)
 126         FDIV.X  FP1,FP0         ...Z/(1+Z)
 127         MOVE.L  a1,d0
 128         AND.L   #$80000000,D0
 129         OR.L    #$3F000000,D0
 130         FADD.X  (sp)+,FP0
 131         MOVE.L  D0,-(sp)
 132
 133         fmove.l d1,fpcr
 134         fmul.s  (sp)+,fp0       ;last fp inst - possible exceptions set
 135
 136         bra     t_frcinx
 137
 138 SINHBIG:
 139         cmp.l   #$400CB2B3,D0
 140         bgt     t_ovfl
 141         FABS.X  FP0
 142         FSUB.D  T1(pc),FP0      ...(|X|-16381LOG2_LEAD)
 143         clr.l   -(sp)
 144         move.l  #$80000000,-(sp)
 145         move.l  a1,d0
 146         AND.L   #$80000000,D0
 147         OR.L    #$7FFB0000,D0
 148         MOVE.L  D0,-(sp)        ...EXTENDED FMT
 149         FSUB.D  T2(pc),FP0      ...|X| - 16381 LOG2, ACCURATE
 150
 151         move.l  d1,-(sp)
 152         clr.l   d1
 153         fmovem.x fp0,(a0)
 154         bsr     setox
 155         fmove.l (sp)+,fpcr
 156
 157         fmul.x  (sp)+,fp0       ;possible exception
 158         bra     t_frcinx
 159
 160         end