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