1# Copyright (c) 1985 Regents of the University of California. 2# All rights reserved. 3# 4# %sccs.include.redist.sh% 5# 6# @(#)argred.s 5.4 (Berkeley) 10/09/90 7# 8 .data 9 .align 2 10_sccsid: 11.asciz "@(#)argred.s 1.1 (Berkeley) 8/21/85; 5.4 (ucb.elefunt) 10/09/90" 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 27libm$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 35abs1: 36 cmpd r3,$0d+4.55530934770520019583e+01 37 blss small_arg 38 jsb trigred 39 rsb 40small_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# 50libm$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 71sine: 72 muld2 r1,r1 # Xsq = X * X 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 761: 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 86cosine: 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 94zero_arg: 95 subd3 r6,$0d1.0,r0 # C(X) = 1 - Xsq 96done: 97 blbc r9,even 98 mnegd r0,r0 99even: 100 rsb 101 102.data 103.align 2 104 105sin_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 115cos_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# 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 142trailing: 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 174middle: 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 206leading: 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 238twoOverPi: 239 .double 0d+6.36619772367581343076e-01 240 .text 241 .align 1 242 243table_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 263abs2: 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# 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 290bits2opi: 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 317trigred: 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 342posmid: 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 370signoff1: 371 subl2 $4,r3 372 ashq r2,(r3),r8 373 bgeq signoff2 374 incl r10 375signoff2: 376 subl2 $4,r3 377 ashq r2,(r3),r7 378 bgeq signoff3 379 incl r9 380signoff3: 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 410more1: 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 425signoff4: 426 emul r1,r6,$0,r4 427 emul r0,r6,r5,r5 428 429addbits1: 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 443more2: 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 459signoff5: 460 emul r0,r6,r4,r5 461 462addbits2: 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# 473result: 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 506small: 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# 535tiny: 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# 552tinier: 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# 576zero: 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# 592mult: 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 602signoff6: 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 631signoff6a: 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 653signoff7: 654signon: 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 679chop: 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# 714carryout: 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# 732return: 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 741signoff8: 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 754signoff9: 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 764