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