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