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