projects / unix-history / blame
| Commit | Line | Data |
|---|---|---|
| 89b330d0 | 1 | # Copyright (c) 1985 Regents of the University of California. |
| fe5be67b KB |
2 | # All rights reserved. |
| 3 | # | |
| 2bc65d90 KB |
4 | # %sccs.include.redist.sh% |
| 5 | # | |
| 6 | # @(#)argred.s 5.4 (Berkeley) %G% | |
| 89b330d0 | 7 | # |
| 7c0a3811 GK |
8 | .data |
| 9 | .align 2 | |
| 10 | _sccsid: | |
| 2bc65d90 | 11 | .asciz "@(#)argred.s 1.1 (Berkeley) 8/21/85; 5.4 (ucb.elefunt) %G%" |
| 89b330d0 ZAL |
12 | |
| 13 | # libm$argred implements Bob Corbett's argument reduction and | |
| 14 | # libm$sincos implements Peter Tang's double precision sin/cos. | |
| 15 | # | |
| 16 | # Note: The two entry points libm$argred and libm$sincos are meant | |
| 17 | # to be used only by _sin, _cos and _tan. | |
| 18 | # | |
| 19 | # method: true range reduction to [-pi/4,pi/4], P. Tang & B. Corbett | |
| 20 | # S. McDonald, April 4, 1985 | |
| 21 | # | |
| 22 | .globl libm$argred | |
| 23 | .globl libm$sincos | |
| 24 | .text | |
| 25 | .align 1 | |
| 26 | ||
| 27 | libm$argred: | |
| 28 | # | |
| 29 | # Compare the argument with the largest possible that can | |
| 30 | # be reduced by table lookup. r3 := |x| will be used in table_lookup . | |
| 31 | # | |
| 32 | movd r0,r3 | |
| 33 | bgeq abs1 | |
| 34 | mnegd r3,r3 | |
| 35 | abs1: | |
| 36 | cmpd r3,$0d+4.55530934770520019583e+01 | |
| 37 | blss small_arg | |
| 38 | jsb trigred | |
| 39 | rsb | |
| 40 | small_arg: | |
| 41 | jsb table_lookup | |
| 42 | rsb | |
| 43 | # | |
| 44 | # At this point, | |
| 45 | # r0 contains the quadrant number, 0, 1, 2, or 3; | |
| 46 | # r2/r1 contains the reduced argument as a D-format number; | |
| 47 | # r3 contains a F-format extension to the reduced argument; | |
| 48 | # r4 contains a 0 or 1 corresponding to a sin or cos entry. | |
| 49 | # | |
| 50 | libm$sincos: | |
| 51 | # | |
| 52 | # Compensate for a cosine entry by adding one to the quadrant number. | |
| 53 | # | |
| 54 | addl2 r4,r0 | |
| 55 | # | |
| 56 | # Polyd clobbers r5-r0 ; save X in r7/r6 . | |
| 57 | # This can be avoided by rewriting trigred . | |
| 58 | # | |
| 59 | movd r1,r6 | |
| 60 | # | |
| 61 | # Likewise, save alpha in r8 . | |
| 62 | # This can be avoided by rewriting trigred . | |
| 63 | # | |
| 64 | movf r3,r8 | |
| 65 | # | |
| 66 | # Odd or even quadrant? cosine if odd, sine otherwise. | |
| 67 | # Save floor(quadrant/2) in r9 ; it determines the final sign. | |
| 68 | # | |
| 69 | rotl $-1,r0,r9 | |
| 70 | blss cosine | |
| 71 | sine: | |
| 72 | muld2 r1,r1 # Xsq = X * X | |
| 317a83dc ZAL |
73 | cmpw $0x2480,r1 # [zl] Xsq > 2^-56? |
| 74 | blss 1f # [zl] yes, go ahead and do polyd | |
| 75 | clrq r1 # [zl] work around 11/780 FPA polyd bug | |
| 76 | 1: | |
| 89b330d0 ZAL |
77 | polyd r1,$7,sin_coef # Q = P(Xsq) , of deg 7 |
| 78 | mulf3 $0f3.0,r8,r4 # beta = 3 * alpha | |
| 79 | mulf2 r0,r4 # beta = Q * beta | |
| 80 | addf2 r8,r4 # beta = alpha + beta | |
| 81 | muld2 r6,r0 # S(X) = X * Q | |
| 82 | # cvtfd r4,r4 ... r5 = 0 after a polyd. | |
| 83 | addd2 r4,r0 # S(X) = beta + S(X) | |
| 84 | addd2 r6,r0 # S(X) = X + S(X) | |
| 85 | brb done | |
| 86 | cosine: | |
| 87 | muld2 r6,r6 # Xsq = X * X | |
| 88 | beql zero_arg | |
| 89 | mulf2 r1,r8 # beta = X * alpha | |
| 90 | polyd r6,$7,cos_coef # Q = P'(Xsq) , of deg 7 | |
| 91 | subd3 r0,r8,r0 # beta = beta - Q | |
| 92 | subw2 $0x80,r6 # Xsq = Xsq / 2 | |
| 93 | addd2 r0,r6 # Xsq = Xsq + beta | |
| 94 | zero_arg: | |
| 95 | subd3 r6,$0d1.0,r0 # C(X) = 1 - Xsq | |
| 96 | done: | |
| 97 | blbc r9,even | |
| 98 | mnegd r0,r0 | |
| 99 | even: | |
| 100 | rsb | |
| 101 | ||
| 102 | .data | |
| 103 | .align 2 | |
| 104 | ||
| 105 | sin_coef: | |
| 106 | .double 0d-7.53080332264191085773e-13 # s7 = 2^-29 -1.a7f2504ffc49f8.. | |
| 107 | .double 0d+1.60573519267703489121e-10 # s6 = 2^-21 1.611adaede473c8.. | |
| 108 | .double 0d-2.50520965150706067211e-08 # s5 = 2^-1a -1.ae644921ed8382.. | |
| 109 | .double 0d+2.75573191800593885716e-06 # s4 = 2^-13 1.71de3a4b884278.. | |
| 110 | .double 0d-1.98412698411850507950e-04 # s3 = 2^-0d -1.a01a01a0125e7d.. | |
| 111 | .double 0d+8.33333333333325688985e-03 # s2 = 2^-07 1.11111111110e50 | |
| 112 | .double 0d-1.66666666666666664354e-01 # s1 = 2^-03 -1.55555555555554 | |
| 113 | .double 0d+0.00000000000000000000e+00 # s0 = 0 | |
| 114 | ||
| 115 | cos_coef: | |
| 116 | .double 0d-1.13006966202629430300e-11 # s7 = 2^-25 -1.8D9BA04D1374BE.. | |
| 117 | .double 0d+2.08746646574796004700e-09 # s6 = 2^-1D 1.1EE632650350BA.. | |
| 118 | .double 0d-2.75573073031284417300e-07 # s5 = 2^-16 -1.27E4F31411719E.. | |
| 119 | .double 0d+2.48015872682668025200e-05 # s4 = 2^-10 1.A01A0196B902E8.. | |
| 120 | .double 0d-1.38888888888464709200e-03 # s3 = 2^-0A -1.6C16C16C11FACE.. | |
| 121 | .double 0d+4.16666666666664761400e-02 # s2 = 2^-05 1.5555555555539E | |
| 122 | .double 0d+0.00000000000000000000e+00 # s1 = 0 | |
| 123 | .double 0d+0.00000000000000000000e+00 # s0 = 0 | |
| 124 | ||
| 125 | #\f | |
| 126 | # Multiples of pi/2 expressed as the sum of three doubles, | |
| 127 | # | |
| 128 | # trailing: n * pi/2 , n = 0, 1, 2, ..., 29 | |
| 129 | # trailing[n] , | |
| 130 | # | |
| 131 | # middle: n * pi/2 , n = 0, 1, 2, ..., 29 | |
| 132 | # middle[n] , | |
| 133 | # | |
| 134 | # leading: n * pi/2 , n = 0, 1, 2, ..., 29 | |
| 135 | # leading[n] , | |
| 136 | # | |
| 137 | # where | |
| 138 | # leading[n] := (n * pi/2) rounded, | |
| 139 | # middle[n] := (n * pi/2 - leading[n]) rounded, | |
| 140 | # trailing[n] := (( n * pi/2 - leading[n]) - middle[n]) rounded . | |
| 141 | ||
| 142 | trailing: | |
| 143 | .double 0d+0.00000000000000000000e+00 # 0 * pi/2 trailing | |
| 144 | .double 0d+4.33590506506189049611e-35 # 1 * pi/2 trailing | |
| 145 | .double 0d+8.67181013012378099223e-35 # 2 * pi/2 trailing | |
| 146 | .double 0d+1.30077151951856714215e-34 # 3 * pi/2 trailing | |
| 147 | .double 0d+1.73436202602475619845e-34 # 4 * pi/2 trailing | |
| 148 | .double 0d-1.68390735624352669192e-34 # 5 * pi/2 trailing | |
| 149 | .double 0d+2.60154303903713428430e-34 # 6 * pi/2 trailing | |
| 150 | .double 0d-8.16726343231148352150e-35 # 7 * pi/2 trailing | |
| 151 | .double 0d+3.46872405204951239689e-34 # 8 * pi/2 trailing | |
| 152 | .double 0d+3.90231455855570147991e-34 # 9 * pi/2 trailing | |
| 153 | .double 0d-3.36781471248705338384e-34 # 10 * pi/2 trailing | |
| 154 | .double 0d-1.06379439835298071785e-33 # 11 * pi/2 trailing | |
| 155 | .double 0d+5.20308607807426856861e-34 # 12 * pi/2 trailing | |
| 156 | .double 0d+5.63667658458045770509e-34 # 13 * pi/2 trailing | |
| 157 | .double 0d-1.63345268646229670430e-34 # 14 * pi/2 trailing | |
| 158 | .double 0d-1.19986217995610764801e-34 # 15 * pi/2 trailing | |
| 159 | .double 0d+6.93744810409902479378e-34 # 16 * pi/2 trailing | |
| 160 | .double 0d-8.03640094449267300110e-34 # 17 * pi/2 trailing | |
| 161 | .double 0d+7.80462911711140295982e-34 # 18 * pi/2 trailing | |
| 162 | .double 0d-7.16921993148029483506e-34 # 19 * pi/2 trailing | |
| 163 | .double 0d-6.73562942497410676769e-34 # 20 * pi/2 trailing | |
| 164 | .double 0d-6.30203891846791677593e-34 # 21 * pi/2 trailing | |
| 165 | .double 0d-2.12758879670596143570e-33 # 22 * pi/2 trailing | |
| 166 | .double 0d+2.53800212047402350390e-33 # 23 * pi/2 trailing | |
| 167 | .double 0d+1.04061721561485371372e-33 # 24 * pi/2 trailing | |
| 168 | .double 0d+6.11729905311472319056e-32 # 25 * pi/2 trailing | |
| 169 | .double 0d+1.12733531691609154102e-33 # 26 * pi/2 trailing | |
| 170 | .double 0d-3.70049587943078297272e-34 # 27 * pi/2 trailing | |
| 171 | .double 0d-3.26690537292459340860e-34 # 28 * pi/2 trailing | |
| 172 | .double 0d-1.14812616507957271361e-34 # 29 * pi/2 trailing | |
| 173 | ||
| 174 | middle: | |
| 175 | .double 0d+0.00000000000000000000e+00 # 0 * pi/2 middle | |
| 176 | .double 0d+5.72118872610983179676e-18 # 1 * pi/2 middle | |
| 177 | .double 0d+1.14423774522196635935e-17 # 2 * pi/2 middle | |
| 178 | .double 0d-3.83475850529283316309e-17 # 3 * pi/2 middle | |
| 179 | .double 0d+2.28847549044393271871e-17 # 4 * pi/2 middle | |
| 180 | .double 0d-2.69052076007086676522e-17 # 5 * pi/2 middle | |
| 181 | .double 0d-7.66951701058566632618e-17 # 6 * pi/2 middle | |
| 182 | .double 0d-1.54628301484890040587e-17 # 7 * pi/2 middle | |
| 183 | .double 0d+4.57695098088786543741e-17 # 8 * pi/2 middle | |
| 184 | .double 0d+1.07001849766246313192e-16 # 9 * pi/2 middle | |
| 185 | .double 0d-5.38104152014173353044e-17 # 10 * pi/2 middle | |
| 186 | .double 0d-2.14622680169080983801e-16 # 11 * pi/2 middle | |
| 187 | .double 0d-1.53390340211713326524e-16 # 12 * pi/2 middle | |
| 188 | .double 0d-9.21580002543456677056e-17 # 13 * pi/2 middle | |
| 189 | .double 0d-3.09256602969780081173e-17 # 14 * pi/2 middle | |
| 190 | .double 0d+3.03066796603896507006e-17 # 15 * pi/2 middle | |
| 191 | .double 0d+9.15390196177573087482e-17 # 16 * pi/2 middle | |
| 192 | .double 0d+1.52771359575124969107e-16 # 17 * pi/2 middle | |
| 193 | .double 0d+2.14003699532492626384e-16 # 18 * pi/2 middle | |
| 194 | .double 0d-1.68853170360202329427e-16 # 19 * pi/2 middle | |
| 195 | .double 0d-1.07620830402834670609e-16 # 20 * pi/2 middle | |
| 196 | .double 0d+3.97700719404595604379e-16 # 21 * pi/2 middle | |
| 197 | .double 0d-4.29245360338161967602e-16 # 22 * pi/2 middle | |
| 198 | .double 0d-3.68013020380794313406e-16 # 23 * pi/2 middle | |
| 199 | .double 0d-3.06780680423426653047e-16 # 24 * pi/2 middle | |
| 200 | .double 0d-2.45548340466059054318e-16 # 25 * pi/2 middle | |
| 201 | .double 0d-1.84316000508691335411e-16 # 26 * pi/2 middle | |
| 202 | .double 0d-1.23083660551323675053e-16 # 27 * pi/2 middle | |
| 203 | .double 0d-6.18513205939560162346e-17 # 28 * pi/2 middle | |
| 204 | .double 0d-6.18980636588357585202e-19 # 29 * pi/2 middle | |
| 205 | ||
| 206 | leading: | |
| 207 | .double 0d+0.00000000000000000000e+00 # 0 * pi/2 leading | |
| 208 | .double 0d+1.57079632679489661351e+00 # 1 * pi/2 leading | |
| 209 | .double 0d+3.14159265358979322702e+00 # 2 * pi/2 leading | |
| 210 | .double 0d+4.71238898038468989604e+00 # 3 * pi/2 leading | |
| 211 | .double 0d+6.28318530717958645404e+00 # 4 * pi/2 leading | |
| 212 | .double 0d+7.85398163397448312306e+00 # 5 * pi/2 leading | |
| 213 | .double 0d+9.42477796076937979208e+00 # 6 * pi/2 leading | |
| 214 | .double 0d+1.09955742875642763501e+01 # 7 * pi/2 leading | |
| 215 | .double 0d+1.25663706143591729081e+01 # 8 * pi/2 leading | |
| 216 | .double 0d+1.41371669411540694661e+01 # 9 * pi/2 leading | |
| 217 | .double 0d+1.57079632679489662461e+01 # 10 * pi/2 leading | |
| 218 | .double 0d+1.72787595947438630262e+01 # 11 * pi/2 leading | |
| 219 | .double 0d+1.88495559215387595842e+01 # 12 * pi/2 leading | |
| 220 | .double 0d+2.04203522483336561422e+01 # 13 * pi/2 leading | |
| 221 | .double 0d+2.19911485751285527002e+01 # 14 * pi/2 leading | |
| 222 | .double 0d+2.35619449019234492582e+01 # 15 * pi/2 leading | |
| 223 | .double 0d+2.51327412287183458162e+01 # 16 * pi/2 leading | |
| 224 | .double 0d+2.67035375555132423742e+01 # 17 * pi/2 leading | |
| 225 | .double 0d+2.82743338823081389322e+01 # 18 * pi/2 leading | |
| 226 | .double 0d+2.98451302091030359342e+01 # 19 * pi/2 leading | |
| 227 | .double 0d+3.14159265358979324922e+01 # 20 * pi/2 leading | |
| 228 | .double 0d+3.29867228626928286062e+01 # 21 * pi/2 leading | |
| 229 | .double 0d+3.45575191894877260523e+01 # 22 * pi/2 leading | |
| 230 | .double 0d+3.61283155162826226103e+01 # 23 * pi/2 leading | |
| 231 | .double 0d+3.76991118430775191683e+01 # 24 * pi/2 leading | |
| 232 | .double 0d+3.92699081698724157263e+01 # 25 * pi/2 leading | |
| 233 | .double 0d+4.08407044966673122843e+01 # 26 * pi/2 leading | |
| 234 | .double 0d+4.24115008234622088423e+01 # 27 * pi/2 leading | |
| 235 | .double 0d+4.39822971502571054003e+01 # 28 * pi/2 leading | |
| 236 | .double 0d+4.55530934770520019583e+01 # 29 * pi/2 leading | |
| 237 | ||
| 238 | twoOverPi: | |
| 239 | .double 0d+6.36619772367581343076e-01 | |
| 240 | .text | |
| 241 | .align 1 | |
| 242 | ||
| 243 | table_lookup: | |
| 244 | muld3 r3,twoOverPi,r0 | |
| 245 | cvtrdl r0,r0 # n = nearest int to ((2/pi)*|x|) rnded | |
| 246 | mull3 $8,r0,r5 | |
| 247 | subd2 leading(r5),r3 # p = (|x| - leading n*pi/2) exactly | |
| 248 | subd3 middle(r5),r3,r1 # q = (p - middle n*pi/2) rounded | |
| 249 | subd2 r1,r3 # r = (p - q) | |
| 250 | subd2 middle(r5),r3 # r = r - middle n*pi/2 | |
| 251 | subd2 trailing(r5),r3 # r = r - trailing n*pi/2 rounded | |
| 252 | # | |
| 253 | # If the original argument was negative, | |
| 254 | # negate the reduce argument and | |
| 255 | # adjust the octant/quadrant number. | |
| 256 | # | |
| 257 | tstw 4(ap) | |
| 258 | bgeq abs2 | |
| 259 | mnegf r1,r1 | |
| 260 | mnegf r3,r3 | |
| 261 | # subb3 r0,$8,r0 ...used for pi/4 reduction -S.McD | |
| 262 | subb3 r0,$4,r0 | |
| 263 | abs2: | |
| 264 | # | |
| 265 | # Clear all unneeded octant/quadrant bits. | |
| 266 | # | |
| 267 | # bicb2 $0xf8,r0 ...used for pi/4 reduction -S.McD | |
| 268 | bicb2 $0xfc,r0 | |
| 269 | rsb | |
| 270 | #\f | |
| 271 | # p.0 | |
| 272 | .text | |
| 273 | .align 2 | |
| 274 | # | |
| 275 | # Only 256 (actually 225) bits of 2/pi are needed for VAX double | |
| 276 | # precision; this was determined by enumerating all the nearest | |
| 277 | # machine integer multiples of pi/2 using continued fractions. | |
| 278 | # (8a8d3673775b7ff7 required the most bits.) -S.McD | |
| 279 | # | |
| 280 | .long 0 | |
| 281 | .long 0 | |
| 282 | .long 0xaef1586d | |
| 283 | .long 0x9458eaf7 | |
| 284 | .long 0x10e4107f | |
| 285 | .long 0xd8a5664f | |
| 286 | .long 0x4d377036 | |
| 287 | .long 0x09d5f47d | |
| 288 | .long 0x91054a7f | |
| 289 | .long 0xbe60db93 | |
| 290 | bits2opi: | |
| 291 | .long 0x00000028 | |
| 292 | .long 0 | |
| 293 | # | |
| 294 | # Note: wherever you see the word `octant', read `quadrant'. | |
| 295 | # Currently this code is set up for pi/2 argument reduction. | |
| 296 | # By uncommenting/commenting the appropriate lines, it will | |
| 297 | # also serve as a pi/4 argument reduction code. | |
| 298 | # | |
| 299 | ||
| 300 | # p.1 | |
| 301 | # Trigred preforms argument reduction | |
| 302 | # for the trigonometric functions. It | |
| 303 | # takes one input argument, a D-format | |
| 304 | # number in r1/r0 . The magnitude of | |
| 305 | # the input argument must be greater | |
| 306 | # than or equal to 1/2 . Trigred produces | |
| 307 | # three results: the number of the octant | |
| 308 | # occupied by the argument, the reduced | |
| 309 | # argument, and an extension of the | |
| 310 | # reduced argument. The octant number is | |
| 311 | # returned in r0 . The reduced argument | |
| 312 | # is returned as a D-format number in | |
| 313 | # r2/r1 . An 8 bit extension of the | |
| 314 | # reduced argument is returned as an | |
| 315 | # F-format number in r3. | |
| 316 | # p.2 | |
| 317 | trigred: | |
| 318 | # | |
| 319 | # Save the sign of the input argument. | |
| 320 | # | |
| 321 | movw r0,-(sp) | |
| 322 | # | |
| 323 | # Extract the exponent field. | |
| 324 | # | |
| 325 | extzv $7,$7,r0,r2 | |
| 326 | # | |
| 327 | # Convert the fraction part of the input | |
| 328 | # argument into a quadword integer. | |
| 329 | # | |
| 330 | bicw2 $0xff80,r0 | |
| 331 | bisb2 $0x80,r0 # -S.McD | |
| 332 | rotl $16,r0,r0 | |
| 333 | rotl $16,r1,r1 | |
| 334 | # | |
| 335 | # If r1 is negative, add 1 to r0 . This | |
| 336 | # adjustment is made so that the two's | |
| 337 | # complement multiplications done later | |
| 338 | # will produce unsigned results. | |
| 339 | # | |
| 340 | bgeq posmid | |
| 341 | incl r0 | |
| 342 | posmid: | |
| 343 | # p.3 | |
| 344 | # | |
| 345 | # Set r3 to the address of the first quadword | |
| 346 | # used to obtain the needed portion of 2/pi . | |
| 347 | # The address is longword aligned to ensure | |
| 348 | # efficient access. | |
| 349 | # | |
| 350 | ashl $-3,r2,r3 | |
| 351 | bicb2 $3,r3 | |
| 352 | subl3 r3,$bits2opi,r3 | |
| 353 | # | |
| 354 | # Set r2 to the size of the shift needed to | |
| 355 | # obtain the correct portion of 2/pi . | |
| 356 | # | |
| 357 | bicb2 $0xe0,r2 | |
| 358 | # p.4 | |
| 359 | # | |
| 360 | # Move the needed 128 bits of 2/pi into | |
| 361 | # r11 - r8 . Adjust the numbers to allow | |
| 362 | # for unsigned multiplication. | |
| 363 | # | |
| 364 | ashq r2,(r3),r10 | |
| 365 | ||
| 366 | subl2 $4,r3 | |
| 367 | ashq r2,(r3),r9 | |
| 368 | bgeq signoff1 | |
| 369 | incl r11 | |
| 370 | signoff1: | |
| 371 | subl2 $4,r3 | |
| 372 | ashq r2,(r3),r8 | |
| 373 | bgeq signoff2 | |
| 374 | incl r10 | |
| 375 | signoff2: | |
| 376 | subl2 $4,r3 | |
| 377 | ashq r2,(r3),r7 | |
| 378 | bgeq signoff3 | |
| 379 | incl r9 | |
| 380 | signoff3: | |
| 381 | # p.5 | |
| 382 | # | |
| 383 | # Multiply the contents of r0/r1 by the | |
| 384 | # slice of 2/pi in r11 - r8 . | |
| 385 | # | |
| 386 | emul r0,r8,$0,r4 | |
| 387 | emul r0,r9,r5,r5 | |
| 388 | emul r0,r10,r6,r6 | |
| 389 | ||
| 390 | emul r1,r8,$0,r7 | |
| 391 | emul r1,r9,r8,r8 | |
| 392 | emul r1,r10,r9,r9 | |
| 393 | emul r1,r11,r10,r10 | |
| 394 | ||
| 395 | addl2 r4,r8 | |
| 396 | adwc r5,r9 | |
| 397 | adwc r6,r10 | |
| 398 | # p.6 | |
| 399 | # | |
| 400 | # If there are more than five leading zeros | |
| 401 | # after the first two quotient bits or if there | |
| 402 | # are more than five leading ones after the first | |
| 403 | # two quotient bits, generate more fraction bits. | |
| 404 | # Otherwise, branch to code to produce the result. | |
| 405 | # | |
| 406 | bicl3 $0xc1ffffff,r10,r4 | |
| 407 | beql more1 | |
| 408 | cmpl $0x3e000000,r4 | |
| 409 | bneq result | |
| 410 | more1: | |
| 411 | # p.7 | |
| 412 | # | |
| 413 | # generate another 32 result bits. | |
| 414 | # | |
| 415 | subl2 $4,r3 | |
| 416 | ashq r2,(r3),r5 | |
| 417 | bgeq signoff4 | |
| 418 | ||
| 419 | emul r1,r6,$0,r4 | |
| 420 | addl2 r1,r5 | |
| 421 | emul r0,r6,r5,r5 | |
| 422 | addl2 r0,r6 | |
| 423 | brb addbits1 | |
| 424 | ||
| 425 | signoff4: | |
| 426 | emul r1,r6,$0,r4 | |
| 427 | emul r0,r6,r5,r5 | |
| 428 | ||
| 429 | addbits1: | |
| 430 | addl2 r5,r7 | |
| 431 | adwc r6,r8 | |
| 432 | adwc $0,r9 | |
| 433 | adwc $0,r10 | |
| 434 | # p.8 | |
| 435 | # | |
| 436 | # Check for massive cancellation. | |
| 437 | # | |
| 438 | bicl3 $0xc0000000,r10,r6 | |
| 439 | # bneq more2 -S.McD Test was backwards | |
| 440 | beql more2 | |
| 441 | cmpl $0x3fffffff,r6 | |
| 442 | bneq result | |
| 443 | more2: | |
| 444 | # p.9 | |
| 445 | # | |
| 446 | # If massive cancellation has occurred, | |
| 447 | # generate another 24 result bits. | |
| 448 | # Testing has shown there will always be | |
| 449 | # enough bits after this point. | |
| 450 | # | |
| 451 | subl2 $4,r3 | |
| 452 | ashq r2,(r3),r5 | |
| 453 | bgeq signoff5 | |
| 454 | ||
| 455 | emul r0,r6,r4,r5 | |
| 456 | addl2 r0,r6 | |
| 457 | brb addbits2 | |
| 458 | ||
| 459 | signoff5: | |
| 460 | emul r0,r6,r4,r5 | |
| 461 | ||
| 462 | addbits2: | |
| 463 | addl2 r6,r7 | |
| 464 | adwc $0,r8 | |
| 465 | adwc $0,r9 | |
| 466 | adwc $0,r10 | |
| 467 | # p.10 | |
| 468 | # | |
| 469 | # The following code produces the reduced | |
| 470 | # argument from the product bits contained | |
| 471 | # in r10 - r7 . | |
| 472 | # | |
| 473 | result: | |
| 474 | # | |
| 475 | # Extract the octant number from r10 . | |
| 476 | # | |
| 477 | # extzv $29,$3,r10,r0 ...used for pi/4 reduction -S.McD | |
| 478 | extzv $30,$2,r10,r0 | |
| 479 | # | |
| 480 | # Clear the octant bits in r10 . | |
| 481 | # | |
| 482 | # bicl2 $0xe0000000,r10 ...used for pi/4 reduction -S.McD | |
| 483 | bicl2 $0xc0000000,r10 | |
| 484 | # | |
| 485 | # Zero the sign flag. | |
| 486 | # | |
| 487 | clrl r5 | |
| 488 | # p.11 | |
| 489 | # | |
| 490 | # Check to see if the fraction is greater than | |
| 491 | # or equal to one-half. If it is, add one | |
| 492 | # to the octant number, set the sign flag | |
| 493 | # on, and replace the fraction with 1 minus | |
| 494 | # the fraction. | |
| 495 | # | |
| 496 | # bitl $0x10000000,r10 ...used for pi/4 reduction -S.McD | |
| 497 | bitl $0x20000000,r10 | |
| 498 | beql small | |
| 499 | incl r0 | |
| 500 | incl r5 | |
| 501 | # subl3 r10,$0x1fffffff,r10 ...used for pi/4 reduction -S.McD | |
| 502 | subl3 r10,$0x3fffffff,r10 | |
| 503 | mcoml r9,r9 | |
| 504 | mcoml r8,r8 | |
| 505 | mcoml r7,r7 | |
| 506 | small: | |
| 507 | # p.12 | |
| 508 | # | |
| 509 | ## Test whether the first 29 bits of the ...used for pi/4 reduction -S.McD | |
| 510 | # Test whether the first 30 bits of the | |
| 511 | # fraction are zero. | |
| 512 | # | |
| 513 | tstl r10 | |
| 514 | beql tiny | |
| 515 | # | |
| 516 | # Find the position of the first one bit in r10 . | |
| 517 | # | |
| 518 | cvtld r10,r1 | |
| 519 | extzv $7,$7,r1,r1 | |
| 520 | # | |
| 521 | # Compute the size of the shift needed. | |
| 522 | # | |
| 523 | subl3 r1,$32,r6 | |
| 524 | # | |
| 525 | # Shift up the high order 64 bits of the | |
| 526 | # product. | |
| 527 | # | |
| 528 | ashq r6,r9,r10 | |
| 529 | ashq r6,r8,r9 | |
| 530 | brb mult | |
| 531 | # p.13 | |
| 532 | # | |
| 533 | # Test to see if the sign bit of r9 is on. | |
| 534 | # | |
| 535 | tiny: | |
| 536 | tstl r9 | |
| 537 | bgeq tinier | |
| 538 | # | |
| 539 | # If it is, shift the product bits up 32 bits. | |
| 540 | # | |
| 541 | movl $32,r6 | |
| 542 | movq r8,r10 | |
| 543 | tstl r10 | |
| 544 | brb mult | |
| 545 | # p.14 | |
| 546 | # | |
| 547 | # Test whether r9 is zero. It is probably | |
| 548 | # impossible for both r10 and r9 to be | |
| 549 | # zero, but until proven to be so, the test | |
| 550 | # must be made. | |
| 551 | # | |
| 552 | tinier: | |
| 553 | beql zero | |
| 554 | # | |
| 555 | # Find the position of the first one bit in r9 . | |
| 556 | # | |
| 557 | cvtld r9,r1 | |
| 558 | extzv $7,$7,r1,r1 | |
| 559 | # | |
| 560 | # Compute the size of the shift needed. | |
| 561 | # | |
| 562 | subl3 r1,$32,r1 | |
| 563 | addl3 $32,r1,r6 | |
| 564 | # | |
| 565 | # Shift up the high order 64 bits of the | |
| 566 | # product. | |
| 567 | # | |
| 568 | ashq r1,r8,r10 | |
| 569 | ashq r1,r7,r9 | |
| 570 | brb mult | |
| 571 | # p.15 | |
| 572 | # | |
| 573 | # The following code sets the reduced | |
| 574 | # argument to zero. | |
| 575 | # | |
| 576 | zero: | |
| 577 | clrl r1 | |
| 578 | clrl r2 | |
| 579 | clrl r3 | |
| 580 | brw return | |
| 581 | # p.16 | |
| 582 | # | |
| 583 | # At this point, r0 contains the octant number, | |
| 584 | # r6 indicates the number of bits the fraction | |
| 585 | # has been shifted, r5 indicates the sign of | |
| 586 | # the fraction, r11/r10 contain the high order | |
| 587 | # 64 bits of the fraction, and the condition | |
| 588 | # codes indicate where the sign bit of r10 | |
| 589 | # is on. The following code multiplies the | |
| 590 | # fraction by pi/2 . | |
| 591 | # | |
| 592 | mult: | |
| 593 | # | |
| 594 | # Save r11/r10 in r4/r1 . -S.McD | |
| 595 | movl r11,r4 | |
| 596 | movl r10,r1 | |
| 597 | # | |
| 598 | # If the sign bit of r10 is on, add 1 to r11 . | |
| 599 | # | |
| 600 | bgeq signoff6 | |
| 601 | incl r11 | |
| 602 | signoff6: | |
| 603 | # p.17 | |
| 604 | # | |
| 605 | # Move pi/2 into r3/r2 . | |
| 606 | # | |
| 607 | movq $0xc90fdaa22168c235,r2 | |
| 608 | # | |
| 609 | # Multiply the fraction by the portion of pi/2 | |
| 610 | # in r2 . | |
| 611 | # | |
| 612 | emul r2,r10,$0,r7 | |
| 613 | emul r2,r11,r8,r7 | |
| 614 | # | |
| 615 | # Multiply the fraction by the portion of pi/2 | |
| 616 | # in r3 . | |
| 617 | emul r3,r10,$0,r9 | |
| 618 | emul r3,r11,r10,r10 | |
| 619 | # | |
| 620 | # Add the product bits together. | |
| 621 | # | |
| 622 | addl2 r7,r9 | |
| 623 | adwc r8,r10 | |
| 624 | adwc $0,r11 | |
| 625 | # | |
| 626 | # Compensate for not sign extending r8 above.-S.McD | |
| 627 | # | |
| 628 | tstl r8 | |
| 629 | bgeq signoff6a | |
| 630 | decl r11 | |
| 631 | signoff6a: | |
| 632 | # | |
| 633 | # Compensate for r11/r10 being unsigned. -S.McD | |
| 634 | # | |
| 635 | addl2 r2,r10 | |
| 636 | adwc r3,r11 | |
| 637 | # | |
| 638 | # Compensate for r3/r2 being unsigned. -S.McD | |
| 639 | # | |
| 640 | addl2 r1,r10 | |
| 641 | adwc r4,r11 | |
| 642 | # p.18 | |
| 643 | # | |
| 644 | # If the sign bit of r11 is zero, shift the | |
| 645 | # product bits up one bit and increment r6 . | |
| 646 | # | |
| 647 | blss signon | |
| 648 | incl r6 | |
| 649 | ashq $1,r10,r10 | |
| 650 | tstl r9 | |
| 651 | bgeq signoff7 | |
| 652 | incl r10 | |
| 653 | signoff7: | |
| 654 | signon: | |
| 655 | # p.19 | |
| 656 | # | |
| 657 | # Shift the 56 most significant product | |
| 658 | # bits into r9/r8 . The sign extension | |
| 659 | # will be handled later. | |
| 660 | # | |
| 661 | ashq $-8,r10,r8 | |
| 662 | # | |
| 663 | # Convert the low order 8 bits of r10 | |
| 664 | # into an F-format number. | |
| 665 | # | |
| 666 | cvtbf r10,r3 | |
| 667 | # | |
| 668 | # If the result of the conversion was | |
| 669 | # negative, add 1 to r9/r8 . | |
| 670 | # | |
| 671 | bgeq chop | |
| 672 | incl r8 | |
| 673 | adwc $0,r9 | |
| 674 | # | |
| 675 | # If r9 is now zero, branch to special | |
| 676 | # code to handle that possibility. | |
| 677 | # | |
| 678 | beql carryout | |
| 679 | chop: | |
| 680 | # p.20 | |
| 681 | # | |
| 682 | # Convert the number in r9/r8 into | |
| 683 | # D-format number in r2/r1 . | |
| 684 | # | |
| 685 | rotl $16,r8,r2 | |
| 686 | rotl $16,r9,r1 | |
| 687 | # | |
| 688 | # Set the exponent field to the appropriate | |
| 689 | # value. Note that the extra bits created by | |
| 690 | # sign extension are now eliminated. | |
| 691 | # | |
| 692 | subw3 r6,$131,r6 | |
| 693 | insv r6,$7,$9,r1 | |
| 694 | # | |
| 695 | # Set the exponent field of the F-format | |
| 696 | # number in r3 to the appropriate value. | |
| 697 | # | |
| 698 | tstf r3 | |
| 699 | beql return | |
| 700 | # extzv $7,$8,r3,r4 -S.McD | |
| 701 | extzv $7,$7,r3,r4 | |
| 702 | addw2 r4,r6 | |
| 703 | # subw2 $217,r6 -S.McD | |
| 704 | subw2 $64,r6 | |
| 705 | insv r6,$7,$8,r3 | |
| 706 | brb return | |
| 707 | # p.21 | |
| 708 | # | |
| 709 | # The following code generates the appropriate | |
| 710 | # result for the unlikely possibility that | |
| 711 | # rounding the number in r9/r8 resulted in | |
| 712 | # a carry out. | |
| 713 | # | |
| 714 | carryout: | |
| 715 | clrl r1 | |
| 716 | clrl r2 | |
| 717 | subw3 r6,$132,r6 | |
| 718 | insv r6,$7,$9,r1 | |
| 719 | tstf r3 | |
| 720 | beql return | |
| 721 | extzv $7,$8,r3,r4 | |
| 722 | addw2 r4,r6 | |
| 723 | subw2 $218,r6 | |
| 724 | insv r6,$7,$8,r3 | |
| 725 | # p.22 | |
| 726 | # | |
| 727 | # The following code makes an needed | |
| 728 | # adjustments to the signs of the | |
| 729 | # results or to the octant number, and | |
| 730 | # then returns. | |
| 731 | # | |
| 732 | return: | |
| 733 | # | |
| 734 | # Test if the fraction was greater than or | |
| 735 | # equal to 1/2 . If so, negate the reduced | |
| 736 | # argument. | |
| 737 | # | |
| 738 | blbc r5,signoff8 | |
| 739 | mnegf r1,r1 | |
| 740 | mnegf r3,r3 | |
| 741 | signoff8: | |
| 742 | # p.23 | |
| 743 | # | |
| 744 | # If the original argument was negative, | |
| 745 | # negate the reduce argument and | |
| 746 | # adjust the octant number. | |
| 747 | # | |
| 748 | tstw (sp)+ | |
| 749 | bgeq signoff9 | |
| 750 | mnegf r1,r1 | |
| 751 | mnegf r3,r3 | |
| 752 | # subb3 r0,$8,r0 ...used for pi/4 reduction -S.McD | |
| 753 | subb3 r0,$4,r0 | |
| 754 | signoff9: | |
| 755 | # | |
| 756 | # Clear all unneeded octant bits. | |
| 757 | # | |
| 758 | # bicb2 $0xf8,r0 ...used for pi/4 reduction -S.McD | |
| 759 | bicb2 $0xfc,r0 | |
| 760 | # | |
| 761 | # Return. | |
| 762 | # | |
| 763 | rsb |