home *** CD-ROM | disk | FTP | other *** search
/ Frozen Fish 1: Amiga / FrozenFish-Apr94.iso / bbs / gnu / libm-src.lha / src / amiga / libm / 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-17    # 10 * pi/2  middle
  212.     .double    0d-2.14622680169080983801e-16    # 11 * pi/2  middle
  213.     .double    0d-1.53390340211713326524e-16    # 12 * pi/2  middle
  214.     .double    0d-9.21580002543456677056e-17    # 13 * pi/2  middle
  215.     .double    0d-3.09256602969780081173e-17    # 14 * pi/2  middle
  216.     .double    0d+3.03066796603896507006e-17    # 15 * pi/2  middle
  217.     .double    0d+9.15390196177573087482e-17    # 16 * pi/2  middle
  218.     .double    0d+1.52771359575124969107e-16    # 17 * pi/2  middle
  219.     .double    0d+2.14003699532492626384e-16    # 18 * pi/2  middle
  220.     .double    0d-1.68853170360202329427e-16    # 19 * pi/2  middle
  221.     .double    0d-1.07620830402834670609e-16    # 20 * pi/2  middle
  222.     .double    0d+3.97700719404595604379e-16    # 21 * pi/2  middle
  223.     .double    0d-4.29245360338161967602e-16    # 22 * pi/2  middle
  224.     .double    0d-3.68013020380794313406e-16    # 23 * pi/2  middle
  225.     .double    0d-3.06780680423426653047e-16    # 24 * pi/2  middle
  226.     .double    0d-2.45548340466059054318e-16    # 25 * pi/2  middle
  227.     .double    0d-1.84316000508691335411e-16    # 26 * pi/2  middle
  228.     .double    0d-1.23083660551323675053e-16    # 27 * pi/2  middle
  229.     .double    0d-6.18513205939560162346e-17    # 28 * pi/2  middle
  230.     .double    0d-6.18980636588357585202e-19    # 29 * pi/2  middle
  231.  
  232. leading:
  233.     .double    0d+0.00000000000000000000e+00    #  0 * pi/2  leading
  234.     .double    0d+1.57079632679489661351e+00    #  1 * pi/2  leading
  235.     .double    0d+3.14159265358979322702e+00    #  2 * pi/2  leading
  236.     .double    0d+4.71238898038468989604e+00    #  3 * pi/2  leading
  237.     .double    0d+6.28318530717958645404e+00    #  4 * pi/2  leading
  238.     .double    0d+7.85398163397448312306e+00    #  5 * pi/2  leading
  239.     .double    0d+9.42477796076937979208e+00    #  6 * pi/2  leading
  240.     .double    0d+1.09955742875642763501e+01    #  7 * pi/2  leading
  241.     .double    0d+1.25663706143591729081e+01    #  8 * pi/2  leading
  242.     .double    0d+1.41371669411540694661e+01    #  9 * pi/2  leading
  243.     .double    0d+1.57079632679489662461e+01    # 10 * pi/2  leading
  244.     .double    0d+1.72787595947438630262e+01    # 11 * pi/2  leading
  245.     .double    0d+1.88495559215387595842e+01    # 12 * pi/2  leading
  246.     .double    0d+2.04203522483336561422e+01    # 13 * pi/2  leading
  247.     .double    0d+2.19911485751285527002e+01    # 14 * pi/2  leading
  248.     .double    0d+2.35619449019234492582e+01    # 15 * pi/2  leading
  249.     .double    0d+2.51327412287183458162e+01    # 16 * pi/2  leading
  250.     .double    0d+2.67035375555132423742e+01    # 17 * pi/2  leading
  251.     .double    0d+2.82743338823081389322e+01    # 18 * pi/2  leading
  252.     .double    0d+2.98451302091030359342e+01    # 19 * pi/2  leading
  253.     .double    0d+3.14159265358979324922e+01    # 20 * pi/2  leading
  254.     .double    0d+3.29867228626928286062e+01    # 21 * pi/2  leading
  255.     .double    0d+3.45575191894877260523e+01    # 22 * pi/2  leading
  256.     .double    0d+3.61283155162826226103e+01    # 23 * pi/2  leading
  257.     .double    0d+3.76991118430775191683e+01    # 24 * pi/2  leading
  258.     .double    0d+3.92699081698724157263e+01    # 25 * pi/2  leading
  259.     .double    0d+4.08407044966673122843e+01    # 26 * pi/2  leading
  260.     .double    0d+4.24115008234622088423e+01    # 27 * pi/2  leading
  261.     .double    0d+4.39822971502571054003e+01    # 28 * pi/2  leading
  262.     .double    0d+4.55530934770520019583e+01    # 29 * pi/2  leading
  263.  
  264. twoOverPi:
  265.     .double    0d+6.36619772367581343076e-01
  266.     .text
  267.     .align    1
  268.  
  269. table_lookup:
  270.     muld3    r3,twoOverPi,r0
  271.     cvtrdl    r0,r0            # n = nearest int to ((2/pi)*|x|) rnded
  272.     mull3    $8,r0,r5
  273.     subd2    leading(r5),r3        # p = (|x| - leading n*pi/2) exactly
  274.     subd3    middle(r5),r3,r1    # q = (p - middle  n*pi/2) rounded
  275.     subd2    r1,r3            # r = (p - q)
  276.     subd2    middle(r5),r3        # r =  r - middle  n*pi/2
  277.     subd2    trailing(r5),r3        # r =  r - trailing n*pi/2  rounded
  278. #
  279. #  If the original argument was negative,
  280. #  negate the reduce argument and
  281. #  adjust the octant/quadrant number.
  282. #
  283.     tstw    4(ap)
  284.     bgeq    abs2
  285.     mnegf    r1,r1
  286.     mnegf    r3,r3
  287. #    subb3    r0,$8,r0    ...used for  pi/4  reduction -S.McD
  288.     subb3    r0,$4,r0
  289. abs2:
  290. #
  291. #  Clear all unneeded octant/quadrant bits.
  292. #
  293. #    bicb2    $0xf8,r0    ...used for  pi/4  reduction -S.McD
  294.     bicb2    $0xfc,r0
  295.     rsb
  296. #
  297. #                        p.0
  298.     .text
  299.     .align    2
  300. #
  301. # Only 256 (actually 225) bits of 2/pi are needed for VAX double
  302. # precision; this was determined by enumerating all the nearest
  303. # machine integer multiples of pi/2 using continued fractions.
  304. # (8a8d3673775b7ff7 required the most bits.)        -S.McD
  305. #
  306.     .long    0
  307.     .long    0
  308.     .long    0xaef1586d
  309.     .long    0x9458eaf7
  310.     .long    0x10e4107f
  311.     .long    0xd8a5664f
  312.     .long    0x4d377036
  313.     .long    0x09d5f47d
  314.     .long    0x91054a7f
  315.     .long    0xbe60db93
  316. bits2opi:
  317.     .long    0x00000028
  318.     .long    0
  319. #
  320. #  Note: wherever you see the word `octant', read `quadrant'.
  321. #  Currently this code is set up for  pi/2  argument reduction.
  322. #  By uncommenting/commenting the appropriate lines, it will
  323. #  also serve as a  pi/4  argument reduction code.
  324. #  
  325.  
  326. #                        p.1
  327. #  Trigred  preforms argument reduction
  328. #  for the trigonometric functions.  It
  329. #  takes one input argument, a D-format
  330. #  number in  r1/r0 .  The magnitude of
  331. #  the input argument must be greater
  332. #  than or equal to  1/2 .  Trigred produces
  333. #  three results:  the number of the octant
  334. #  occupied by the argument, the reduced 
  335. #  argument, and an extension of the
  336. #  reduced argument.  The octant number is 
  337. #  returned in  r0 .  The reduced argument 
  338. #  is returned as a D-format number in 
  339. #  r2/r1 .  An 8 bit extension of the 
  340. #  reduced argument is returned as an 
  341. #  F-format number in r3.
  342. #                        p.2
  343. trigred:
  344. #
  345. #  Save the sign of the input argument.
  346. #
  347.     movw    r0,-(sp)
  348. #
  349. #  Extract the exponent field.
  350. #
  351.     extzv    $7,$7,r0,r2
  352. #
  353. #  Convert the fraction part of the input
  354. #  argument into a quadword integer.
  355. #
  356.     bicw2    $0xff80,r0
  357.     bisb2    $0x80,r0    # -S.McD
  358.     rotl    $16,r0,r0
  359.     rotl    $16,r1,r1
  360. #
  361. #  If  r1  is negative, add  1  to  r0 .  This
  362. #  adjustment is made so that the two's
  363. #  complement multiplications done later
  364. #  will produce unsigned results.
  365. #
  366.     bgeq    posmid
  367.     incl    r0
  368. posmid:
  369. #                        p.3
  370. #
  371. #  Set  r3  to the address of the first quadword
  372. #  used to obtain the needed portion of  2/pi .
  373. #  The address is longword aligned to ensure
  374. #  efficient access.
  375. #
  376.     ashl    $-3,r2,r3
  377.     bicb2    $3,r3
  378.     subl3    r3,$bits2opi,r3
  379. #
  380. #  Set  r2  to the size of the shift needed to 
  381. #  obtain the correct portion of  2/pi .
  382. #
  383.     bicb2    $0xe0,r2
  384. #                        p.4
  385. #
  386. #  Move the needed  128  bits of  2/pi  into
  387. #  r11 - r8 .  Adjust the numbers to allow
  388. #  for unsigned multiplication.
  389. #
  390.     ashq    r2,(r3),r10
  391.  
  392.     subl2    $4,r3
  393.     ashq    r2,(r3),r9
  394.     bgeq    signoff1
  395.     incl    r11
  396. signoff1:
  397.     subl2    $4,r3
  398.     ashq    r2,(r3),r8
  399.     bgeq    signoff2
  400.     incl    r10
  401. signoff2:
  402.     subl2    $4,r3
  403.     ashq    r2,(r3),r7
  404.     bgeq    signoff3
  405.     incl    r9
  406. signoff3:
  407. #                        p.5
  408. #
  409. #  Multiply the contents of  r0/r1  by the 
  410. #  slice of  2/pi  in  r11 - r8 .
  411. #
  412.     emul    r0,r8,$0,r4
  413.     emul    r0,r9,r5,r5
  414.     emul    r0,r10,r6,r6
  415.  
  416.     emul    r1,r8,$0,r7
  417.     emul    r1,r9,r8,r8
  418.     emul    r1,r10,r9,r9
  419.     emul    r1,r11,r10,r10
  420.  
  421.     addl2    r4,r8
  422.     adwc    r5,r9
  423.     adwc    r6,r10
  424. #                        p.6
  425. #
  426. #  If there are more than five leading zeros
  427. #  after the first two quotient bits or if there
  428. #  are more than five leading ones after the first
  429. #  two quotient bits, generate more fraction bits.
  430. #  Otherwise, branch to code to produce the result.
  431. #
  432.     bicl3    $0xc1ffffff,r10,r4
  433.     beql    more1
  434.     cmpl    $0x3e000000,r4
  435.     bneq    result
  436. more1:
  437. #                        p.7
  438. #
  439. #  generate another  32  result bits.
  440. #
  441.     subl2    $4,r3
  442.     ashq    r2,(r3),r5
  443.     bgeq    signoff4
  444.  
  445.     emul    r1,r6,$0,r4
  446.     addl2    r1,r5
  447.     emul    r0,r6,r5,r5
  448.     addl2    r0,r6
  449.     brb    addbits1
  450.  
  451. signoff4:
  452.     emul    r1,r6,$0,r4
  453.     emul    r0,r6,r5,r5
  454.  
  455. addbits1:
  456.     addl2    r5,r7
  457.     adwc    r6,r8
  458.     adwc    $0,r9
  459.     adwc    $0,r10
  460. #                        p.8
  461. #
  462. #  Check for massive cancellation.
  463. #
  464.     bicl3    $0xc0000000,r10,r6
  465. #    bneq    more2            -S.McD  Test was backwards
  466.     beql    more2
  467.     cmpl    $0x3fffffff,r6
  468.     bneq    result
  469. more2:
  470. #                        p.9
  471. #
  472. #  If massive cancellation has occurred,
  473. #  generate another  24  result bits.
  474. #  Testing has shown there will always be 
  475. #  enough bits after this point.
  476. #
  477.     subl2    $4,r3
  478.     ashq    r2,(r3),r5
  479.     bgeq    signoff5
  480.  
  481.     emul    r0,r6,r4,r5
  482.     addl2    r0,r6
  483.     brb    addbits2
  484.  
  485. signoff5:
  486.     emul    r0,r6,r4,r5
  487.  
  488. addbits2:
  489.     addl2    r6,r7
  490.     adwc    $0,r8
  491.     adwc    $0,r9
  492.     adwc    $0,r10
  493. #                        p.10
  494. #
  495. #  The following code produces the reduced
  496. #  argument from the product bits contained
  497. #  in  r10 - r7 .
  498. #
  499. result:
  500. #
  501. #  Extract the octant number from  r10 .
  502. #
  503. #    extzv    $29,$3,r10,r0    ...used for  pi/4  reduction -S.McD
  504.     extzv    $30,$2,r10,r0
  505. #
  506. #  Clear the octant bits in  r10 .
  507. #
  508. #    bicl2    $0xe0000000,r10    ...used for  pi/4  reduction -S.McD
  509.     bicl2    $0xc0000000,r10
  510. #
  511. #  Zero the sign flag.
  512. #
  513.     clrl    r5
  514. #                        p.11
  515. #
  516. #  Check to see if the fraction is greater than
  517. #  or equal to one-half.  If it is, add one 
  518. #  to the octant number, set the sign flag
  519. #  on, and replace the fraction with  1 minus
  520. #  the fraction.
  521. #
  522. #    bitl    $0x10000000,r10        ...used for  pi/4  reduction -S.McD
  523.     bitl    $0x20000000,r10
  524.     beql    small
  525.     incl    r0
  526.     incl    r5
  527. #    subl3    r10,$0x1fffffff,r10    ...used for  pi/4  reduction -S.McD
  528.     subl3    r10,$0x3fffffff,r10
  529.     mcoml    r9,r9
  530.     mcoml    r8,r8
  531.     mcoml    r7,r7
  532. small:
  533. #                        p.12
  534. #
  535. ##  Test whether the first  29  bits of the ...used for  pi/4  reduction -S.McD
  536. #  Test whether the first  30  bits of the 
  537. #  fraction are zero.
  538. #
  539.     tstl    r10
  540.     beql    tiny
  541. #
  542. #  Find the position of the first one bit in  r10 .
  543. #
  544.     cvtld    r10,r1
  545.     extzv    $7,$7,r1,r1
  546. #
  547. #  Compute the size of the shift needed.
  548. #
  549.     subl3    r1,$32,r6
  550. #
  551. #  Shift up the high order  64  bits of the
  552. #  product.
  553. #
  554.     ashq    r6,r9,r10
  555.     ashq    r6,r8,r9
  556.     brb    mult
  557. #                        p.13
  558. #
  559. #  Test to see if the sign bit of  r9  is on.
  560. #
  561. tiny:
  562.     tstl    r9
  563.     bgeq    tinier
  564. #
  565. #  If it is, shift the product bits up  32  bits.
  566. #
  567.     movl    $32,r6
  568.     movq    r8,r10
  569.     tstl    r10
  570.     brb    mult
  571. #                        p.14
  572. #
  573. #  Test whether  r9  is zero.  It is probably
  574. #  impossible for both  r10  and  r9  to be
  575. #  zero, but until proven to be so, the test
  576. #  must be made.
  577. #
  578. tinier:
  579.     beql    zero
  580. #
  581. #  Find the position of the first one bit in  r9 .
  582. #
  583.     cvtld    r9,r1
  584.     extzv    $7,$7,r1,r1
  585. #
  586. #  Compute the size of the shift needed.
  587. #
  588.     subl3    r1,$32,r1
  589.     addl3    $32,r1,r6
  590. #
  591. #  Shift up the high order  64  bits of the
  592. #  product.
  593. #
  594.     ashq    r1,r8,r10
  595.     ashq    r1,r7,r9
  596.     brb    mult
  597. #                        p.15
  598. #
  599. #  The following code sets the reduced
  600. #  argument to zero.
  601. #
  602. zero:
  603.     clrl    r1
  604.     clrl    r2
  605.     clrl    r3
  606.     brw    return
  607. #                        p.16
  608. #
  609. #  At this point,  r0  contains the octant number,
  610. #  r6  indicates the number of bits the fraction
  611. #  has been shifted,  r5  indicates the sign of
  612. #  the fraction,  r11/r10  contain the high order
  613. #  64  bits of the fraction, and the condition
  614. #  codes indicate where the sign bit of  r10
  615. #  is on.  The following code multiplies the
  616. #  fraction by  pi/2 .
  617. #
  618. mult:
  619. #
  620. #  Save  r11/r10  in  r4/r1 .        -S.McD
  621.     movl    r11,r4
  622.     movl    r10,r1
  623. #
  624. #  If the sign bit of  r10  is on, add  1  to  r11 .
  625. #
  626.     bgeq    signoff6
  627.     incl    r11
  628. signoff6:
  629. #                        p.17
  630. #
  631. #  Move  pi/2  into  r3/r2 .
  632. #
  633.     movq    $0xc90fdaa22168c235,r2
  634. #
  635. #  Multiply the fraction by the portion of  pi/2
  636. #  in  r2 .
  637. #
  638.     emul    r2,r10,$0,r7
  639.     emul    r2,r11,r8,r7
  640. #
  641. #  Multiply the fraction by the portion of  pi/2 
  642. #  in  r3 .
  643.     emul    r3,r10,$0,r9
  644.     emul    r3,r11,r10,r10
  645. #
  646. #  Add the product bits together.
  647. #
  648.     addl2    r7,r9
  649.     adwc    r8,r10
  650.     adwc    $0,r11
  651. #
  652. #  Compensate for not sign extending  r8  above.-S.McD
  653. #
  654.     tstl    r8
  655.     bgeq    signoff6a
  656.     decl    r11
  657. signoff6a:
  658. #
  659. #  Compensate for  r11/r10  being unsigned.    -S.McD
  660. #
  661.     addl2    r2,r10
  662.     adwc    r3,r11
  663. #
  664. #  Compensate for  r3/r2  being unsigned.    -S.McD
  665. #
  666.     addl2    r1,r10
  667.     adwc    r4,r11
  668. #                        p.18
  669. #
  670. #  If the sign bit of  r11  is zero, shift the
  671. #  product bits up one bit and increment  r6 .
  672. #
  673.     blss    signon
  674.     incl    r6
  675.     ashq    $1,r10,r10
  676.     tstl    r9
  677.     bgeq    signoff7
  678.     incl    r10
  679. signoff7:
  680. signon:
  681. #                        p.19
  682. #
  683. #  Shift the  56  most significant product
  684. #  bits into  r9/r8 .  The sign extension
  685. #  will be handled later.
  686. #
  687.     ashq    $-8,r10,r8
  688. #
  689. #  Convert the low order  8  bits of  r10
  690. #  into an F-format number.
  691. #
  692.     cvtbf    r10,r3
  693. #
  694. #  If the result of the conversion was
  695. #  negative, add  1  to  r9/r8 .
  696. #
  697.     bgeq    chop
  698.     incl    r8
  699.     adwc    $0,r9
  700. #
  701. #  If  r9  is now zero, branch to special
  702. #  code to handle that possibility.
  703. #
  704.     beql    carryout
  705. chop:
  706. #                        p.20
  707. #
  708. #  Convert the number in  r9/r8  into
  709. #  D-format number in  r2/r1 .
  710. #
  711.     rotl    $16,r8,r2
  712.     rotl    $16,r9,r1
  713. #
  714. #  Set the exponent field to the appropriate
  715. #  value.  Note that the extra bits created by
  716. #  sign extension are now eliminated.
  717. #
  718.     subw3    r6,$131,r6
  719.     insv    r6,$7,$9,r1
  720. #
  721. #  Set the exponent field of the F-format
  722. #  number in  r3  to the appropriate value.
  723. #
  724.     tstf    r3
  725.     beql    return
  726. #    extzv    $7,$8,r3,r4    -S.McD
  727.     extzv    $7,$7,r3,r4
  728.     addw2    r4,r6
  729. #    subw2    $217,r6        -S.McD
  730.     subw2    $64,r6
  731.     insv    r6,$7,$8,r3
  732.     brb    return
  733. #                        p.21
  734. #
  735. #  The following code generates the appropriate 
  736. #  result for the unlikely possibility that
  737. #  rounding the number in  r9/r8  resulted in 
  738. #  a carry out.
  739. #
  740. carryout:
  741.     clrl    r1
  742.     clrl    r2
  743.     subw3    r6,$132,r6
  744.     insv    r6,$7,$9,r1
  745.     tstf    r3
  746.     beql    return
  747.     extzv    $7,$8,r3,r4
  748.     addw2    r4,r6
  749.     subw2    $218,r6
  750.     insv    r6,$7,$8,r3
  751. #                        p.22
  752. #
  753. #  The following code makes an needed
  754. #  adjustments to the signs of the 
  755. #  results or to the octant number, and
  756. #  then returns.
  757. #
  758. return:
  759. #
  760. #  Test if the fraction was greater than or 
  761. #  equal to  1/2 .  If so, negate the reduced
  762. #  argument.
  763. #
  764.     blbc    r5,signoff8
  765.     mnegf    r1,r1
  766.     mnegf    r3,r3
  767. signoff8:
  768. #                        p.23
  769. #
  770. #  If the original argument was negative,
  771. #  negate the reduce argument and
  772. #  adjust the octant number.
  773. #
  774.     tstw    (sp)+
  775.     bgeq    signoff9
  776.     mnegf    r1,r1
  777.     mnegf    r3,r3
  778. #    subb3    r0,$8,r0    ...used for  pi/4  reduction -S.McD
  779.     subb3    r0,$4,r0
  780. signoff9:
  781. #
  782. #  Clear all unneeded octant bits.
  783. #
  784. #    bicb2    $0xf8,r0    ...used for  pi/4  reduction -S.McD
  785.     bicb2    $0xfc,r0
  786. #
  787. #  Return.
  788. #
  789.     rsb
  790.