sys/arch/m68k/fpsp/satanh.sa

   1 *       $NetBSD$
   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 *       satanh.sa 3.3 12/19/90
  35 *
  36 *       The entry point satanh computes the inverse
  37 *       hyperbolic tangent of
  38 *       an input argument; satanhd does the same except for denormalized
  39 *       input.
  40 *
  41 *       Input: Double-extended number X in location pointed to
  42 *               by address register a0.
  43 *
  44 *       Output: The value arctanh(X) returned in floating-point register Fp0.
  45 *
  46 *       Accuracy and Monotonicity: The returned result is within 3 ulps in
  47 *               64 significant bit, i.e. within 0.5001 ulp to 53 bits if the
  48 *               result is subsequently rounded to double precision. The
  49 *               result is provably monotonic in double precision.
  50 *
  51 *       Speed: The program satanh takes approximately 270 cycles.
  52 *
  53 *       Algorithm:
  54 *
  55 *       ATANH
  56 *       1. If |X| >= 1, go to 3.
  57 *
  58 *       2. (|X| < 1) Calculate atanh(X) by
  59 *               sgn := sign(X)
  60 *               y := |X|
  61 *               z := 2y/(1-y)
  62 *               atanh(X) := sgn * (1/2) * logp1(z)
  63 *               Exit.
  64 *
  65 *       3. If |X| > 1, go to 5.
  66 *
  67 *       4. (|X| = 1) Generate infinity with an appropriate sign and
  68 *               divide-by-zero by
  69 *               sgn := sign(X)
  70 *               atan(X) := sgn / (+0).
  71 *               Exit.
  72 *
  73 *       5. (|X| > 1) Generate an invalid operation by 0 * infinity.
  74 *               Exit.
  75 *
  76
  77 satanh  IDNT    2,1 Motorola 040 Floating Point Software Package
  78
  79         section 8
  80
  81         xref    t_dz
  82         xref    t_operr
  83         xref    t_frcinx
  84         xref    t_extdnrm
  85         xref    slognp1
  86
  87         xdef    satanhd
  88 satanhd:
  89 *--ATANH(X) = X FOR DENORMALIZED X
  90
  91         bra             t_extdnrm
  92
  93         xdef    satanh
  94 satanh:
  95         move.l          (a0),d0
  96         move.w          4(a0),d0
  97         ANDI.L          #$7FFFFFFF,D0
  98         CMPI.L          #$3FFF8000,D0
  99         BGE.B           ATANHBIG
 100
 101 *--THIS IS THE USUAL CASE, |X| < 1
 102 *--Y = |X|, Z = 2Y/(1-Y), ATANH(X) = SIGN(X) * (1/2) * LOG1P(Z).
 103
 104         FABS.X          (a0),FP0        ...Y = |X|
 105         FMOVE.X         FP0,FP1
 106         FNEG.X          FP1             ...-Y
 107         FADD.X          FP0,FP0         ...2Y
 108         FADD.S          #:3F800000,FP1  ...1-Y
 109         FDIV.X          FP1,FP0         ...2Y/(1-Y)
 110         move.l          (a0),d0
 111         ANDI.L          #$80000000,D0
 112         ORI.L           #$3F000000,D0   ...SIGN(X)*HALF
 113         move.l          d0,-(sp)
 114
 115         fmovem.x        fp0,(a0)        ...overwrite input
 116         move.l          d1,-(sp)
 117         clr.l           d1
 118         bsr             slognp1         ...LOG1P(Z)
 119         fmove.l         (sp)+,fpcr
 120         FMUL.S          (sp)+,FP0
 121         bra             t_frcinx
 122
 123 ATANHBIG:
 124         FABS.X          (a0),FP0        ...|X|
 125         FCMP.S          #:3F800000,FP0
 126         fbgt            t_operr
 127         bra             t_dz
 128
 129         end