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