home *** CD-ROM | disk | FTP | other *** search
/ Fresh Fish 4 / FreshFish_May-June1994.bin / bsd / src / libm / libm-amiga / vax / argred.s next >
Text File  |  1993-09-23  |  21KB  |  790 lines
  1. # Copyright (c) 1985 Regents of the University of California.
  2. # All rights reserved.
  3. #
  4. # Redistribution and use in source and binary forms, with or without
  5. # modification, are permitted provided that the following conditions
  6. # are met:
  7. # 1. Redistributions of source code must retain the above copyright
  8. #    notice, this list of conditions and the following disclaimer.
  9. # 2. Redistributions in binary form must reproduce the above copyright
  10. #    notice, this list of conditions and the following disclaimer in the
  11. #    documentation and/or other materials provided with the distribution.
  12. # 3. All advertising materials mentioning features or use of this software
  13. #    must display the following acknowledgement:
  14. #    This product includes software developed by the University of
  15. #    California, Berkeley and its contributors.
  16. # 4. Neither the name of the University nor the names of its contributors
  17. #    may be used to endorse or promote products derived from this software
  18. #    without specific prior written permission.
  19. #
  20. # THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
  21. # ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
  22. # IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
  23. # ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
  24. # FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
  25. # DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
  26. # OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
  27. # HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
  28. # LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
  29. # OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
  30. # SUCH DAMAGE.
  31. #
  32. #    @(#)argred.s    5.4 (Berkeley) 10/9/90
  33. #
  34.     .data
  35.     .align    2
  36. _sccsid:
  37. .asciz    "@(#)argred.s    1.1 (Berkeley) 8/21/85; 5.4 (ucb.elefunt) 10/9/90"
  38.  
  39. #  libm$argred implements Bob Corbett's argument reduction and
  40. #  libm$sincos implements Peter Tang's double precision sin/cos.
  41. #  
  42. #  Note: The two entry points libm$argred and libm$sincos are meant
  43. #        to be used only by _sin, _cos and _tan.
  44. #
  45. # method: true range reduction to [-pi/4,pi/4], P. Tang  &  B. Corbett
  46. # S. McDonald, April 4,  1985
  47. #
  48.     .globl    libm$argred
  49.     .globl    libm$sincos
  50.     .text
  51.     .align    1
  52.  
  53. libm$argred:
  54. #
  55. #  Compare the argument with the largest possible that can
  56. #  be reduced by table lookup.  r3 := |x|  will be used in  table_lookup .
  57. #
  58.     movd    r0,r3
  59.     bgeq    abs1
  60.     mnegd    r3,r3
  61. abs1:
  62.     cmpd    r3,$0d+4.55530934770520019583e+01
  63.     blss    small_arg
  64.     jsb    trigred
  65.     rsb
  66. small_arg:
  67.     jsb    table_lookup
  68.     rsb
  69. #
  70. #  At this point,
  71. #       r0  contains the quadrant number, 0, 1, 2, or 3;
  72. #    r2/r1  contains the reduced argument as a D-format number;
  73. #         r3  contains a F-format extension to the reduced argument;
  74. #          r4  contains a  0 or 1  corresponding to a  sin or cos  entry.
  75. #
  76. libm$sincos:
  77. #
  78. #  Compensate for a cosine entry by adding one to the quadrant number.
  79. #
  80.     addl2    r4,r0
  81. #
  82. #  Polyd clobbers  r5-r0 ;  save  X  in  r7/r6 .
  83. #  This can be avoided by rewriting  trigred .
  84. #
  85.     movd    r1,r6
  86. #
  87. #  Likewise, save  alpha  in  r8 .
  88. #  This can be avoided by rewriting  trigred .
  89. #
  90.     movf    r3,r8
  91. #
  92. #  Odd or even quadrant?  cosine if odd, sine otherwise.
  93. #  Save  floor(quadrant/2) in  r9  ; it determines the final sign.
  94. #
  95.     rotl    $-1,r0,r9
  96.     blss    cosine
  97. sine:
  98.     muld2    r1,r1        # Xsq = X * X
  99.     cmpw    $0x2480,r1    # [zl] Xsq > 2^-56?
  100.     blss    1f        # [zl] yes, go ahead and do polyd
  101.     clrq    r1        # [zl] work around 11/780 FPA polyd bug
  102. 1:
  103.     polyd    r1,$7,sin_coef    # Q = P(Xsq) , of deg 7
  104.     mulf3    $0f3.0,r8,r4    # beta = 3 * alpha
  105.     mulf2    r0,r4        # beta = Q * beta
  106.     addf2    r8,r4        # beta = alpha + beta
  107.     muld2    r6,r0        # S(X) = X * Q
  108. #    cvtfd    r4,r4        ... r5 = 0 after a polyd.
  109.     addd2    r4,r0        # S(X) = beta + S(X)
  110.     addd2    r6,r0        # S(X) = X + S(X)
  111.     brb    done
  112. cosine:
  113.     muld2    r6,r6        # Xsq = X * X
  114.     beql    zero_arg
  115.     mulf2    r1,r8        # beta = X * alpha
  116.     polyd    r6,$7,cos_coef    # Q = P'(Xsq) , of deg 7
  117.     subd3    r0,r8,r0    # beta = beta - Q
  118.     subw2    $0x80,r6    # Xsq = Xsq / 2
  119.     addd2    r0,r6        # Xsq = Xsq + beta
  120. zero_arg:
  121.     subd3    r6,$0d1.0,r0    # C(X) = 1 - Xsq
  122. done:
  123.     blbc    r9,even
  124.     mnegd    r0,r0
  125. even:
  126.     rsb
  127.  
  128. .data
  129. .align    2
  130.  
  131. sin_coef:
  132.     .double    0d-7.53080332264191085773e-13    # s7 = 2^-29 -1.a7f2504ffc49f8..
  133.     .double    0d+1.60573519267703489121e-10    # s6 = 2^-21  1.611adaede473c8..
  134.     .double    0d-2.50520965150706067211e-08    # s5 = 2^-1a -1.ae644921ed8382..
  135.     .double    0d+2.75573191800593885716e-06    # s4 = 2^-13  1.71de3a4b884278..
  136.     .double    0d-1.98412698411850507950e-04    # s3 = 2^-0d -1.a01a01a0125e7d..
  137.     .double    0d+8.33333333333325688985e-03    # s2 = 2^-07  1.11111111110e50
  138.     .double    0d-1.66666666666666664354e-01    # s1 = 2^-03 -1.55555555555554
  139.     .double    0d+0.00000000000000000000e+00    # s0 = 0
  140.  
  141. cos_coef:
  142.     .double    0d-1.13006966202629430300e-11    # s7 = 2^-25 -1.8D9BA04D1374BE..
  143.     .double    0d+2.08746646574796004700e-09    # s6 = 2^-1D  1.1EE632650350BA..
  144.     .double    0d-2.75573073031284417300e-07    # s5 = 2^-16 -1.27E4F31411719E..
  145.     .double    0d+2.48015872682668025200e-05    # s4 = 2^-10  1.A01A0196B902E8..
  146.     .double    0d-1.38888888888464709200e-03    # s3 = 2^-0A -1.6C16C16C11FACE..
  147.     .double    0d+4.16666666666664761400e-02    # s2 = 2^-05  1.5555555555539E
  148.     .double    0d+0.00000000000000000000e+00    # s1 = 0
  149.     .double    0d+0.00000000000000000000e+00    # s0 = 0
  150.  
  151. #
  152. #  Multiples of  pi/2  expressed as the sum of three doubles,
  153. #
  154. #  trailing:    n * pi/2 ,  n = 0, 1, 2, ..., 29
  155. #            trailing[n] ,
  156. #
  157. #  middle:    n * pi/2 ,  n = 0, 1, 2, ..., 29
  158. #            middle[n]   ,
  159. #
  160. #  leading:    n * pi/2 ,  n = 0, 1, 2, ..., 29
  161. #            leading[n]  ,
  162. #
  163. #    where
  164. #        leading[n]  := (n * pi/2)  rounded,
  165. #        middle[n]   := (n * pi/2  -  leading[n])  rounded,
  166. #        trailing[n] := (( n * pi/2 - leading[n]) - middle[n])  rounded .
  167.  
  168. trailing:
  169.     .double    0d+0.00000000000000000000e+00    #  0 * pi/2  trailing
  170.     .double    0d+4.33590506506189049611e-35    #  1 * pi/2  trailing
  171.     .double    0d+8.67181013012378099223e-35    #  2 * pi/2  trailing
  172.     .double    0d+1.30077151951856714215e-34    #  3 * pi/2  trailing
  173.     .double    0d+1.73436202602475619845e-34    #  4 * pi/2  trailing
  174.     .double    0d-1.68390735624352669192e-34    #  5 * pi/2  trailing
  175.     .double    0d+2.60154303903713428430e-34    #  6 * pi/2  trailing
  176.     .double    0d-8.16726343231148352150e-35    #  7 * pi/2  trailing
  177.     .double    0d+3.46872405204951239689e-34    #  8 * pi/2  trailing
  178.     .double    0d+3.90231455855570147991e-34    #  9 * pi/2  trailing
  179.     .double    0d-3.36781471248705338384e-34    # 10 * pi/2  trailing
  180.     .double    0d-1.06379439835298071785e-33    # 11 * pi/2  trailing
  181.     .double    0d+5.20308607807426856861e-34    # 12 * pi/2  trailing
  182.     .double    0d+5.63667658458045770509e-34    # 13 * pi/2  trailing
  183.     .double    0d-1.63345268646229670430e-34    # 14 * pi/2  trailing
  184.     .double    0d-1.19986217995610764801e-34    # 15 * pi/2  trailing
  185.     .double    0d+6.93744810409902479378e-34    # 16 * pi/2  trailing
  186.     .double    0d-8.03640094449267300110e-34    # 17 * pi/2  trailing
  187.     .double    0d+7.80462911711140295982e-34    # 18 * pi/2  trailing
  188.     .double    0d-7.16921993148029483506e-34    # 19 * pi/2  trailing
  189.     .double    0d-6.73562942497410676769e-34    # 20 * pi/2  trailing
  190.     .double    0d-6.30203891846791677593e-34    # 21 * pi/2  trailing
  191.     .double    0d-2.12758879670596143570e-33    # 22 * pi/2  trailing
  192.     .double    0d+2.53800212047402350390e-33    # 23 * pi/2  trailing
  193.     .double    0d+1.04061721561485371372e-33    # 24 * pi/2  trailing
  194.     .double    0d+6.11729905311472319056e-32    # 25 * pi/2  trailing
  195.     .double    0d+1.12733531691609154102e-33    # 26 * pi/2  trailing
  196.     .double    0d-3.70049587943078297272e-34    # 27 * pi/2  trailing
  197.     .double    0d-3.26690537292459340860e-34    # 28 * pi/2  trailing
  198.     .double    0d-1.14812616507957271361e-34    # 29 * pi/2  trailing
  199.  
  200. middle:
  201.     .double    0d+0.00000000000000000000e+00    #  0 * pi/2  middle
  202.     .double    0d+5.72118872610983179676e-18    #  1 * pi/2  middle
  203.     .double    0d+1.14423774522196635935e-17    #  2 * pi/2  middle
  204.     .double    0d-3.83475850529283316309e-17    #  3 * pi/2  middle
  205.     .double    0d+2.28847549044393271871e-17    #  4 * pi/2  middle
  206.     .double    0d-2.69052076007086676522e-17    #  5 * pi/2  middle
  207.     .double    0d-7.66951701058566632618e-17    #  6 * pi/2  middle
  208.     .double    0d-1.54628301484890040587e-17    #  7 * pi/2  middle
  209.     .double    0d+4.57695098088786543741e-17    #  8 * pi/2  middle
  210.     .double    0d+1.07001849766246313192e-16    #  9 * pi/2  middle
  211.     .double    0d-5.38104152014173353044e-1