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

   1 *       $NetBSD: sasin.sa,v 1.2 1994/10/26 07:49:29 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 *       sasin.sa 3.3 12/19/90
  35 *
  36 *       Description: The entry point sAsin computes the inverse sine of
  37 *               an input argument; sAsind 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 arcsin(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 sASIN takes approximately 310 cycles.
  51 *
  52 *       Algorithm:
  53 *
  54 *       ASIN
  55 *       1. If |X| >= 1, go to 3.
  56 *
  57 *       2. (|X| < 1) Calculate asin(X) by
  58 *               z := sqrt( [1-X][1+X] )
  59 *               asin(X) = atan( x / z ).
  60 *               Exit.
  61 *
  62 *       3. If |X| > 1, go to 5.
  63 *
  64 *       4. (|X| = 1) sgn := sign(X), return asin(X) := sgn * Pi/2. Exit.
  65 *
  66 *       5. (|X| > 1) Generate an invalid operation by 0 * infinity.
  67 *               Exit.
  68 *
  69
  70 SASIN   IDNT    2,1 Motorola 040 Floating Point Software Package
  71
  72         section 8
  73
  74 PIBY2   DC.L $3FFF0000,$C90FDAA2,$2168C235,$00000000
  75
  76         xref    t_operr
  77         xref    t_frcinx
  78         xref    t_extdnrm
  79         xref    satan
  80
  81         xdef    sasind
  82 sasind:
  83 *--ASIN(X) = X FOR DENORMALIZED X
  84
  85         bra             t_extdnrm
  86
  87         xdef    sasin
  88 sasin:
  89         FMOVE.X         (a0),FP0        ...LOAD INPUT
  90
  91         move.l          (a0),d0
  92         move.w          4(a0),d0
  93         ANDI.L          #$7FFFFFFF,D0
  94         CMPI.L          #$3FFF8000,D0
  95         BGE.B           asinbig
  96
  97 *--THIS IS THE USUAL CASE, |X| < 1
  98 *--ASIN(X) = ATAN( X / SQRT( (1-X)(1+X) ) )
  99
 100         FMOVE.S         #:3F800000,FP1
 101         FSUB.X          FP0,FP1         ...1-X
 102         fmovem.x        fp2,-(a7)
 103         FMOVE.S         #:3F800000,FP2
 104         FADD.X          FP0,FP2         ...1+X
 105         FMUL.X          FP2,FP1         ...(1+X)(1-X)
 106         fmovem.x        (a7)+,fp2
 107         FSQRT.X         FP1             ...SQRT([1-X][1+X])
 108         FDIV.X          FP1,FP0         ...X/SQRT([1-X][1+X])
 109         fmovem.x        fp0,(a0)
 110         bsr             satan
 111         bra             t_frcinx
 112
 113 asinbig:
 114         FABS.X          FP0      ...|X|
 115         FCMP.S          #:3F800000,FP0
 116         fbgt            t_operr         ;cause an operr exception
 117
 118 *--|X| = 1, ASIN(X) = +- PI/2.
 119
 120         FMOVE.X         PIBY2,FP0
 121         move.l          (a0),d0
 122         ANDI.L          #$80000000,D0   ...SIGN BIT OF X
 123         ORI.L           #$3F800000,D0   ...+-1 IN SGL FORMAT
 124         MOVE.L          D0,-(sp)        ...push SIGN(X) IN SGL-FMT
 125         FMOVE.L         d1,FPCR
 126         FMUL.S          (sp)+,FP0
 127         bra             t_frcinx
 128
 129         end