3 // Copyright (C) 2000, 2001, Intel Corporation
4 // All rights reserved.
6 // Contributed 2/2/2000 by John Harrison, Ted Kubaska, Bob Norin, Shane Story,
7 // and Ping Tak Peter Tang of the Computational Software Lab, Intel Corporation.
9 // Redistribution and use in source and binary forms, with or without
10 // modification, are permitted provided that the following conditions are
13 // * Redistributions of source code must retain the above copyright
14 // notice, this list of conditions and the following disclaimer.
16 // * Redistributions in binary form must reproduce the above copyright
17 // notice, this list of conditions and the following disclaimer in the
18 // documentation and/or other materials provided with the distribution.
20 // * The name of Intel Corporation may not be used to endorse or promote
21 // products derived from this software without specific prior written
24 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
25 // "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
26 // LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
27 // A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL INTEL OR ITS
28 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
29 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
30 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
31 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY
32 // OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY OR TORT (INCLUDING
33 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
34 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 // Intel Corporation is the author of this code, and requests that all
37 // problem reports or change requests be submitted to it directly at
38 // http://developer.intel.com/opensource.
41 //==============================================================
42 // 2/02/00 Initial revision
43 // 4/02/00 Unwind support added.
44 // 6/16/00 Updated tables to enforce symmetry
45 // 8/31/00 Saved 2 cycles in main path, and 9 in other paths.
46 // 9/20/00 The updated tables regressed to an old version, so reinstated them
47 // 10/18/00 Changed one table entry to ensure symmetry
48 // 1/03/01 Improved speed, fixed flag settings for small arguments.
51 //==============================================================
52 // double sin( double x);
53 // double cos( double x);
55 // Overview of operation
56 //==============================================================
60 // Reduce x to region -1/2*pi/2^k ===== 0 ===== +1/2*pi/2^k where k=4
61 // divide x by pi/2^k.
62 // Multiply by 2^k/pi.
63 // nfloat = Round result to integer (round-to-nearest)
65 // r = x - nfloat * pi/2^k
66 // Do this as (x - nfloat * HIGH(pi/2^k)) - nfloat * LOW(pi/2^k) for increased accuracy.
67 // pi/2^k is stored as two numbers that when added make pi/2^k.
68 // pi/2^k = HIGH(pi/2^k) + LOW(pi/2^k)
70 // x = (nfloat * pi/2^k) + r
71 // r is small enough that we can use a polynomial approximation
72 // and is referred to as the reduced argument.
76 // Take the unreduced part and remove the multiples of 2pi.
77 // So nfloat = nfloat (with lower k+1 bits cleared) + lower k+1 bits
79 // nfloat (with lower k+1 bits cleared) is a multiple of 2^(k+1)
81 // nfloat * pi/2^k = N * 2^(k+1) * pi/2^k + (lower k+1 bits) * pi/2^k
82 // nfloat * pi/2^k = N * 2 * pi + (lower k+1 bits) * pi/2^k
83 // nfloat * pi/2^k = N2pi + M * pi/2^k
86 // Sin(x) = Sin((nfloat * pi/2^k) + r)
87 // = Sin(nfloat * pi/2^k) * Cos(r) + Cos(nfloat * pi/2^k) * Sin(r)
89 // Sin(nfloat * pi/2^k) = Sin(N2pi + Mpi/2^k)
90 // = Sin(N2pi)Cos(Mpi/2^k) + Cos(N2pi)Sin(Mpi/2^k)
93 // Cos(nfloat * pi/2^k) = Cos(N2pi + Mpi/2^k)
94 // = Cos(N2pi)Cos(Mpi/2^k) + Sin(N2pi)Sin(Mpi/2^k)
97 // Sin(x) = Sin(Mpi/2^k) Cos(r) + Cos(Mpi/2^k) Sin(r)
103 // There are 2^(k+1) Sin entries in a table.
104 // There are 2^(k+1) Cos entries in a table.
106 // Get Sin(Mpi/2^k) and Cos(Mpi/2^k) by table lookup.
111 // Calculate Cos(r) and Sin(r) by polynomial approximation.
113 // Cos(r) = 1 + r^2 q1 + r^4 q2 + r^6 q3 + ... = Series for Cos
114 // Sin(r) = r + r^3 p1 + r^5 p2 + r^7 p3 + ... = Series for Sin
116 // and the coefficients q1, q2, ... and p1, p2, ... are stored in a table
120 // Sin(x) = Sin(Mpi/2^k) Cos(r) + Cos(Mpi/2^k) Sin(r)
124 // Sm = Sin(Mpi/2^k) and Cm = Cos(Mpi/2^k)
128 // P = p1 + r^2p2 + r^4p3 + r^6p4
129 // Q = q1 + r^2q2 + r^4q3 + r^6q4
132 // Sin(r) = r + rcub * P
133 // = r + r^3p1 + r^5p2 + r^7p3 + r^9p4 + ... = Sin(r)
135 // The coefficients are not exactly these values, but almost.
139 // p3 = -1/5040 = -1/7!
140 // p4 = 1/362889 = 1/9!
144 // Answer = Sm Cos(r) + Cm P
146 // Cos(r) = 1 + rsq Q
147 // Cos(r) = 1 + r^2 Q
148 // Cos(r) = 1 + r^2 (q1 + r^2q2 + r^4q3 + r^6q4)
149 // Cos(r) = 1 + r^2q1 + r^4q2 + r^6q3 + r^8q4 + ...
151 // Sm Cos(r) = Sm(1 + rsq Q)
152 // Sm Cos(r) = Sm + Sm rsq Q
153 // Sm Cos(r) = Sm + s_rsq Q
160 #include "libm_support.h"
163 //==============================================================
164 // general input registers:
168 // predicate registers used:
171 // floating-point registers used
176 //==============================================================
179 sind_int_Nfloat = f11
186 sind_Inv_Pi_by_16 = f32
187 sind_Pi_by_16_hi = f33
188 sind_Pi_by_16_lo = f34
190 sind_Inv_Pi_by_64 = f35
191 sind_Pi_by_64_hi = f36
192 sind_Pi_by_64_lo = f37
217 sind_SIG_INV_PI_BY_16_2TO61 = f55
218 sind_RSHF_2TO61 = f56
222 sind_W_2TO61_RSH = f60
226 /////////////////////////////////////////////////////////////
232 sind_AD_beta_table = r37
238 sind_GR_sig_inv_pi_by_16 = r14
239 sind_GR_rshf_2to61 = r15
241 sind_GR_exp_2tom61 = r17
260 ASM_TYPE_DIRECTIVE(double_sind_pi,@object)
261 // data8 0xA2F9836E4E44152A, 0x00004001 // 16/pi (significand loaded w/ setf)
263 data8 0xC90FDAA22168C234, 0x00003FFC // pi/16 hi
264 // c4c6628b80dc1cd1 29024e088a
265 data8 0xC4C6628B80DC1CD1, 0x00003FBC // pi/16 lo
266 ASM_SIZE_DIRECTIVE(double_sind_pi)
269 ASM_TYPE_DIRECTIVE(double_sind_pq_k4,@object)
270 data8 0x3EC71C963717C63A // P4
271 data8 0x3EF9FFBA8F191AE6 // Q4
272 data8 0xBF2A01A00F4E11A8 // P3
273 data8 0xBF56C16C05AC77BF // Q3
274 data8 0x3F8111111110F167 // P2
275 data8 0x3FA555555554DD45 // Q2
276 data8 0xBFC5555555555555 // P1
277 data8 0xBFDFFFFFFFFFFFFC // Q1
278 ASM_SIZE_DIRECTIVE(double_sind_pq_k4)
281 double_sin_cos_beta_k4:
282 ASM_TYPE_DIRECTIVE(double_sin_cos_beta_k4,@object)
283 data8 0x0000000000000000 , 0x00000000 // sin( 0 pi/16) S0
284 data8 0x8000000000000000 , 0x00003fff // cos( 0 pi/16) C0
286 data8 0xc7c5c1e34d3055b3 , 0x00003ffc // sin( 1 pi/16) S1
287 data8 0xfb14be7fbae58157 , 0x00003ffe // cos( 1 pi/16) C1
289 data8 0xc3ef1535754b168e , 0x00003ffd // sin( 2 pi/16) S2
290 data8 0xec835e79946a3146 , 0x00003ffe // cos( 2 pi/16) C2
292 data8 0x8e39d9cd73464364 , 0x00003ffe // sin( 3 pi/16) S3
293 data8 0xd4db3148750d181a , 0x00003ffe // cos( 3 pi/16) C3
295 data8 0xb504f333f9de6484 , 0x00003ffe // sin( 4 pi/16) S4
296 data8 0xb504f333f9de6484 , 0x00003ffe // cos( 4 pi/16) C4
299 data8 0xd4db3148750d181a , 0x00003ffe // sin( 5 pi/16) C3
300 data8 0x8e39d9cd73464364 , 0x00003ffe // cos( 5 pi/16) S3
302 data8 0xec835e79946a3146 , 0x00003ffe // sin( 6 pi/16) C2
303 data8 0xc3ef1535754b168e , 0x00003ffd // cos( 6 pi/16) S2
305 data8 0xfb14be7fbae58157 , 0x00003ffe // sin( 7 pi/16) C1
306 data8 0xc7c5c1e34d3055b3 , 0x00003ffc // cos( 7 pi/16) S1
308 data8 0x8000000000000000 , 0x00003fff // sin( 8 pi/16) C0
309 data8 0x0000000000000000 , 0x00000000 // cos( 8 pi/16) S0
312 data8 0xfb14be7fbae58157 , 0x00003ffe // sin( 9 pi/16) C1
313 data8 0xc7c5c1e34d3055b3 , 0x0000bffc // cos( 9 pi/16) -S1
315 data8 0xec835e79946a3146 , 0x00003ffe // sin(10 pi/16) C2
316 data8 0xc3ef1535754b168e , 0x0000bffd // cos(10 pi/16) -S2
318 data8 0xd4db3148750d181a , 0x00003ffe // sin(11 pi/16) C3
319 data8 0x8e39d9cd73464364 , 0x0000bffe // cos(11 pi/16) -S3
321 data8 0xb504f333f9de6484 , 0x00003ffe // sin(12 pi/16) S4
322 data8 0xb504f333f9de6484 , 0x0000bffe // cos(12 pi/16) -S4
325 data8 0x8e39d9cd73464364 , 0x00003ffe // sin(13 pi/16) S3
326 data8 0xd4db3148750d181a , 0x0000bffe // cos(13 pi/16) -C3
328 data8 0xc3ef1535754b168e , 0x00003ffd // sin(14 pi/16) S2
329 data8 0xec835e79946a3146 , 0x0000bffe // cos(14 pi/16) -C2
331 data8 0xc7c5c1e34d3055b3 , 0x00003ffc // sin(15 pi/16) S1
332 data8 0xfb14be7fbae58157 , 0x0000bffe // cos(15 pi/16) -C1
334 data8 0x0000000000000000 , 0x00000000 // sin(16 pi/16) S0
335 data8 0x8000000000000000 , 0x0000bfff // cos(16 pi/16) -C0
338 data8 0xc7c5c1e34d3055b3 , 0x0000bffc // sin(17 pi/16) -S1
339 data8 0xfb14be7fbae58157 , 0x0000bffe // cos(17 pi/16) -C1
341 data8 0xc3ef1535754b168e , 0x0000bffd // sin(18 pi/16) -S2
342 data8 0xec835e79946a3146 , 0x0000bffe // cos(18 pi/16) -C2
344 data8 0x8e39d9cd73464364 , 0x0000bffe // sin(19 pi/16) -S3
345 data8 0xd4db3148750d181a , 0x0000bffe // cos(19 pi/16) -C3
347 data8 0xb504f333f9de6484 , 0x0000bffe // sin(20 pi/16) -S4
348 data8 0xb504f333f9de6484 , 0x0000bffe // cos(20 pi/16) -S4
351 data8 0xd4db3148750d181a , 0x0000bffe // sin(21 pi/16) -C3
352 data8 0x8e39d9cd73464364 , 0x0000bffe // cos(21 pi/16) -S3
354 data8 0xec835e79946a3146 , 0x0000bffe // sin(22 pi/16) -C2
355 data8 0xc3ef1535754b168e , 0x0000bffd // cos(22 pi/16) -S2
357 data8 0xfb14be7fbae58157 , 0x0000bffe // sin(23 pi/16) -C1
358 data8 0xc7c5c1e34d3055b3 , 0x0000bffc // cos(23 pi/16) -S1
360 data8 0x8000000000000000 , 0x0000bfff // sin(24 pi/16) -C0
361 data8 0x0000000000000000 , 0x00000000 // cos(24 pi/16) S0
364 data8 0xfb14be7fbae58157 , 0x0000bffe // sin(25 pi/16) -C1
365 data8 0xc7c5c1e34d3055b3 , 0x00003ffc // cos(25 pi/16) S1
367 data8 0xec835e79946a3146 , 0x0000bffe // sin(26 pi/16) -C2
368 data8 0xc3ef1535754b168e , 0x00003ffd // cos(26 pi/16) S2
370 data8 0xd4db3148750d181a , 0x0000bffe // sin(27 pi/16) -C3
371 data8 0x8e39d9cd73464364 , 0x00003ffe // cos(27 pi/16) S3
373 data8 0xb504f333f9de6484 , 0x0000bffe // sin(28 pi/16) -S4
374 data8 0xb504f333f9de6484 , 0x00003ffe // cos(28 pi/16) S4
377 data8 0x8e39d9cd73464364 , 0x0000bffe // sin(29 pi/16) -S3
378 data8 0xd4db3148750d181a , 0x00003ffe // cos(29 pi/16) C3
380 data8 0xc3ef1535754b168e , 0x0000bffd // sin(30 pi/16) -S2
381 data8 0xec835e79946a3146 , 0x00003ffe // cos(30 pi/16) C2
383 data8 0xc7c5c1e34d3055b3 , 0x0000bffc // sin(31 pi/16) -S1
384 data8 0xfb14be7fbae58157 , 0x00003ffe // cos(31 pi/16) C1
386 data8 0x0000000000000000 , 0x00000000 // sin(32 pi/16) S0
387 data8 0x8000000000000000 , 0x00003fff // cos(32 pi/16) C0
388 ASM_SIZE_DIRECTIVE(double_sin_cos_beta_k4)
398 ////////////////////////////////////////////////////////
399 // There are two entry points: sin and cos
402 // If from sin, p8 is true
403 // If from cos, p9 is true
418 alloc r32=ar.pfs,1,13,0,0
419 movl sind_GR_sig_inv_pi_by_16 = 0xA2F9836E4E44152A // significand of 16/pi
422 addl sind_AD_1 = @ltoff(double_sind_pi), gp
423 movl sind_GR_rshf_2to61 = 0x47b8000000000000 // 1.1000 2^(63+63-2)
428 ld8 sind_AD_1 = [sind_AD_1]
429 fnorm sind_NORM_f8 = f8
430 cmp.eq p8,p9 = r0, r0
433 mov sind_GR_exp_2tom61 = 0xffff-61 // exponent of scaling factor 2^-61
434 mov sind_r_sincos = 0x0
435 br.cond.sptk L(SIND_SINCOS)
440 ASM_SIZE_DIRECTIVE(sin)
455 alloc r32=ar.pfs,1,13,0,0
456 movl sind_GR_sig_inv_pi_by_16 = 0xA2F9836E4E44152A // significand of 16/pi
459 addl sind_AD_1 = @ltoff(double_sind_pi), gp
460 movl sind_GR_rshf_2to61 = 0x47b8000000000000 // 1.1000 2^(63+63-2)
465 ld8 sind_AD_1 = [sind_AD_1]
466 fnorm.s1 sind_NORM_f8 = f8
467 cmp.eq p9,p8 = r0, r0
470 mov sind_GR_exp_2tom61 = 0xffff-61 // exponent of scaling factor 2^-61
471 mov sind_r_sincos = 0x8
472 br.cond.sptk L(SIND_SINCOS)
477 ////////////////////////////////////////////////////////
478 // All entry points end up here.
479 // If from sin, sind_r_sincos is 0 and p8 is true
480 // If from cos, sind_r_sincos is 8 = 2^(k-1) and p9 is true
481 // We add sind_r_sincos to N
486 // Form two constants we need
487 // 16/pi * 2^-2 * 2^63, scaled by 2^61 since we just loaded the significand
488 // 1.1000...000 * 2^(63+63-2) to right shift int(W) into the low significand
489 // fcmp used to set denormal, and invalid on snans
491 setf.sig sind_SIG_INV_PI_BY_16_2TO61 = sind_GR_sig_inv_pi_by_16
492 fcmp.eq.s0 p12,p0=f8,f0
493 mov sind_r_17_ones = 0x1ffff
496 setf.d sind_RSHF_2TO61 = sind_GR_rshf_2to61
497 movl sind_GR_rshf = 0x43e8000000000000 // 1.1000 2^63 for right shift
501 // Form another constant
502 // 2^-61 for scaling Nfloat
503 // 0x10009 is register_bias + 10.
504 // So if f8 > 2^10 = Gamma, go to DBX
506 setf.exp sind_2TOM61 = sind_GR_exp_2tom61
507 fclass.m p13,p0 = f8, 0x23 // Test for x inf
508 mov sind_exp_limit = 0x10009
512 // Load the two pieces of pi/16
513 // Form another constant
514 // 1.1000...000 * 2^63, the right shift constant
516 ldfe sind_Pi_by_16_hi = [sind_AD_1],16
517 setf.d sind_RSHF = sind_GR_rshf
518 fclass.m p14,p0 = f8, 0xc3 // Test for x nan
523 ldfe sind_Pi_by_16_lo = [sind_AD_1],16
524 (p13) frcpa.s0 f8,p12=f0,f0 // force qnan indef for x=inf
528 addl sind_AD_beta_table = @ltoff(double_sin_cos_beta_k4), gp
530 (p13) br.ret.spnt b0 ;; // Exit for x=inf
533 // Start loading P, Q coefficients
536 ldfpd sind_P4,sind_Q4 = [sind_AD_1],16
537 (p8) fclass.m.unc p6,p0 = f8, 0x07 // Test for sin(0)
541 addl sind_AD_beta_table = @ltoff(double_sin_cos_beta_k4), gp
542 (p14) fma.d f8=f8,f1,f0 // qnan for x=nan
543 (p14) br.ret.spnt b0 ;; // Exit for x=nan
549 getf.exp sind_r_signexp = f8
550 (p9) fclass.m.unc p7,p0 = f8, 0x07 // Test for sin(0)
554 ld8 sind_AD_beta_table = [sind_AD_beta_table]
560 ldfpd sind_P3,sind_Q3 = [sind_AD_1],16
561 setf.sig fp_tmp = gr_tmp // Create constant such that fmpy sets inexact
562 (p6) br.ret.spnt b0 ;;
566 and sind_r_exp = sind_r_17_ones, sind_r_signexp
567 (p7) fmerge.s f8 = f1,f1
568 (p7) br.ret.spnt b0 ;;
571 // p10 is true if we must call routines to handle larger arguments
572 // p10 is true if f8 exp is > 0x10009
575 ldfpd sind_P2,sind_Q2 = [sind_AD_1],16
577 cmp.ge p10,p0 = sind_r_exp,sind_exp_limit
581 // sind_W = x * sind_Inv_Pi_by_16
582 // Multiply x by scaled 16/pi and add large const to shift integer part of W to
583 // rightmost bits of significand
585 ldfpd sind_P1,sind_Q1 = [sind_AD_1]
586 fma.s1 sind_W_2TO61_RSH = sind_NORM_f8,sind_SIG_INV_PI_BY_16_2TO61,sind_RSHF_2TO61
590 (p10) cmp.ne.unc p11,p12=sind_r_sincos,r0 // p11 call __libm_cos_double_dbx
591 // p12 call __libm_sin_double_dbx
592 (p11) br.cond.spnt L(COSD_DBX)
593 (p12) br.cond.spnt L(SIND_DBX)
598 // sind_NFLOAT = Round_Int_Nearest(sind_W)
599 // This is done by scaling back by 2^-61 and subtracting the shift constant
602 fms.s1 sind_NFLOAT = sind_W_2TO61_RSH,sind_2TOM61,sind_RSHF
607 // get N = (int)sind_int_Nfloat
609 getf.sig sind_GR_n = sind_W_2TO61_RSH
614 // Add 2^(k-1) (which is in sind_r_sincos) to N
615 // sind_r = -sind_Nfloat * sind_Pi_by_16_hi + x
616 // sind_r = sind_r -sind_Nfloat * sind_Pi_by_16_lo
618 add sind_GR_n = sind_GR_n, sind_r_sincos
619 fnma.s1 sind_r = sind_NFLOAT, sind_Pi_by_16_hi, sind_NORM_f8
624 // Get M (least k+1 bits of N)
626 and sind_GR_m = 0x1f,sind_GR_n ;;
628 shl sind_GR_32m = sind_GR_m,5 ;;
631 // Add 32*M to address of sin_cos_beta table
633 add sind_AD_2 = sind_GR_32m, sind_AD_beta_table
639 ldfe sind_Sm = [sind_AD_2],16
640 (p8) fclass.m.unc p10,p0=f8,0x0b // If sin, note denormal input to set uflow
645 ldfe sind_Cm = [sind_AD_2]
646 fnma.s1 sind_r = sind_NFLOAT, sind_Pi_by_16_lo, sind_r
653 fma.s1 sind_rsq = sind_r, sind_r, f0
658 fmpy.s0 fp_tmp = fp_tmp,fp_tmp // fmpy forces inexact flag
662 // form P and Q series
665 fma.s1 sind_P_temp1 = sind_rsq, sind_P4, sind_P3
671 fma.s1 sind_Q_temp1 = sind_rsq, sind_Q4, sind_Q3
675 // get rcube and sm*rsq
678 fmpy.s1 sind_srsq = sind_Sm,sind_rsq
684 fmpy.s1 sind_rcub = sind_r, sind_rsq
690 fma.s1 sind_Q_temp2 = sind_rsq, sind_Q_temp1, sind_Q2
696 fma.s1 sind_P_temp2 = sind_rsq, sind_P_temp1, sind_P2
702 fma.s1 sind_Q = sind_rsq, sind_Q_temp2, sind_Q1
708 fma.s1 sind_P = sind_rsq, sind_P_temp2, sind_P1
715 fma.s1 sind_Q = sind_srsq,sind_Q, sind_Sm
721 fma.s1 sind_P = sind_rcub,sind_P, sind_r
725 // If sin(denormal), force inexact to be set
728 (p10) fmpy.d.s0 fp_tmp = f8,f8
735 fma.d f8 = sind_Cm, sind_P, sind_Q
739 ASM_SIZE_DIRECTIVE(cos#)
743 .proc __libm_callout_1s
750 .save ar.pfs,GR_SAVE_PFS
751 mov GR_SAVE_PFS=ar.pfs
766 br.call.sptk.many b0=__libm_sin_double_dbx# ;;
780 mov ar.pfs = GR_SAVE_PFS
783 .endp __libm_callout_1s
784 ASM_SIZE_DIRECTIVE(__libm_callout_1s)
787 .proc __libm_callout_1c
794 .save ar.pfs,GR_SAVE_PFS
795 mov GR_SAVE_PFS=ar.pfs
810 br.call.sptk.many b0=__libm_cos_double_dbx# ;;
824 mov ar.pfs = GR_SAVE_PFS
827 .endp __libm_callout_1c
828 ASM_SIZE_DIRECTIVE(__libm_callout_1c)
831 // ====================================================================
832 // ====================================================================
834 // These functions calculate the sin and cos for inputs
836 // __libm_sin_double_dbx# and __libm_cos_double_dbx#
838 // *********************************************************************
839 // *********************************************************************
841 // Function: Combined sin(x) and cos(x), where
843 // sin(x) = sine(x), for double precision x values
844 // cos(x) = cosine(x), for double precision x values
846 // *********************************************************************
848 // Accuracy: Within .7 ulps for 80-bit floating point values
849 // Very accurate for double precision values
851 // *********************************************************************
855 // Floating-Point Registers: f8 (Input and Return Value)
858 // General Purpose Registers:
860 // r44-r45 (Used to pass arguments to pi_by_2 reduce routine)
862 // Predicate Registers: p6-p13
864 // *********************************************************************
866 // IEEE Special Conditions:
868 // Denormal fault raised on denormal inputs
869 // Overflow exceptions do not occur
870 // Underflow exceptions raised when appropriate for sin
871 // (No specialized error handling for this routine)
872 // Inexact raised when appropriate by algorithm
883 // *********************************************************************
885 // Mathematical Description
886 // ========================
888 // The computation of FSIN and FCOS is best handled in one piece of
889 // code. The main reason is that given any argument Arg, computation
890 // of trigonometric functions first calculate N and an approximation
893 // Arg = N pi/2 + alpha, |alpha| <= pi/4.
897 // cos( Arg ) = sin( (N+1) pi/2 + alpha ),
899 // therefore, the code for computing sine will produce cosine as long
900 // as 1 is added to N immediately after the argument reduction
908 // Arg = M pi/2 + alpha, |alpha| <= pi/4,
910 // let I = M mod 4, or I be the two lsb of M when M is represented
911 // as 2's complement. I = [i_0 i_1]. Then
913 // sin( Arg ) = (-1)^i_0 sin( alpha ) if i_1 = 0,
914 // = (-1)^i_0 cos( alpha ) if i_1 = 1.
918 // sin ((-pi/2 + alpha) = (-1) cos (alpha)
920 // sin (alpha) = sin (alpha)
922 // sin (pi/2 + alpha) = cos (alpha)
924 // sin (pi + alpha) = (-1) sin (alpha)
926 // sin ((3/2)pi + alpha) = (-1) cos (alpha)
928 // The value of alpha is obtained by argument reduction and
929 // represented by two working precision numbers r and c where
931 // alpha = r + c accurately.
933 // The reduction method is described in a previous write up.
934 // The argument reduction scheme identifies 4 cases. For Cases 2
935 // and 4, because |alpha| is small, sin(r+c) and cos(r+c) can be
936 // computed very easily by 2 or 3 terms of the Taylor series
937 // expansion as follows:
942 // sin(r + c) = r + c - r^3/6 accurately
943 // cos(r + c) = 1 - 2^(-67) accurately
948 // sin(r + c) = r + c - r^3/6 + r^5/120 accurately
949 // cos(r + c) = 1 - r^2/2 + r^4/24 accurately
951 // The only cases left are Cases 1 and 3 of the argument reduction
952 // procedure. These two cases will be merged since after the
953 // argument is reduced in either cases, we have the reduced argument
954 // represented as r + c and that the magnitude |r + c| is not small
955 // enough to allow the usage of a very short approximation.
957 // The required calculation is either
959 // sin(r + c) = sin(r) + correction, or
960 // cos(r + c) = cos(r) + correction.
964 // sin(r + c) = sin(r) + c sin'(r) + O(c^2)
965 // = sin(r) + c cos (r) + O(c^2)
966 // = sin(r) + c(1 - r^2/2) accurately.
969 // cos(r + c) = cos(r) - c sin(r) + O(c^2)
970 // = cos(r) - c(r - r^3/6) accurately.
972 // We therefore concentrate on accurately calculating sin(r) and
973 // cos(r) for a working-precision number r, |r| <= pi/4 to within
976 // The greatest challenge of this task is that the second terms of
979 // r - r^3/3! + r^r/5! - ...
983 // 1 - r^2/2! + r^4/4! - ...
985 // are not very small when |r| is close to pi/4 and the rounding
986 // errors will be a concern if simple polynomial accumulation is
987 // used. When |r| < 2^-3, however, the second terms will be small
988 // enough (6 bits or so of right shift) that a normal Horner
989 // recurrence suffices. Hence there are two cases that we consider
990 // in the accurate computation of sin(r) and cos(r), |r| <= pi/4.
992 // Case small_r: |r| < 2^(-3)
993 // --------------------------
995 // Since Arg = M pi/4 + r + c accurately, and M mod 4 is [i_0 i_1],
998 // sin(Arg) = (-1)^i_0 * sin(r + c) if i_1 = 0
999 // = (-1)^i_0 * cos(r + c) if i_1 = 1
1001 // can be accurately approximated by
1003 // sin(Arg) = (-1)^i_0 * [sin(r) + c] if i_1 = 0
1004 // = (-1)^i_0 * [cos(r) - c*r] if i_1 = 1
1006 // because |r| is small and thus the second terms in the correction
1007 // are unneccessary.
1009 // Finally, sin(r) and cos(r) are approximated by polynomials of
1010 // moderate lengths.
1012 // sin(r) = r + S_1 r^3 + S_2 r^5 + ... + S_5 r^11
1013 // cos(r) = 1 + C_1 r^2 + C_2 r^4 + ... + C_5 r^10
1015 // We can make use of predicates to selectively calculate
1016 // sin(r) or cos(r) based on i_1.
1018 // Case normal_r: 2^(-3) <= |r| <= pi/4
1019 // ------------------------------------
1021 // This case is more likely than the previous one if one considers
1022 // r to be uniformly distributed in [-pi/4 pi/4]. Again,
1024 // sin(Arg) = (-1)^i_0 * sin(r + c) if i_1 = 0
1025 // = (-1)^i_0 * cos(r + c) if i_1 = 1.
1027 // Because |r| is now larger, we need one extra term in the
1028 // correction. sin(Arg) can be accurately approximated by
1030 // sin(Arg) = (-1)^i_0 * [sin(r) + c(1-r^2/2)] if i_1 = 0
1031 // = (-1)^i_0 * [cos(r) - c*r*(1 - r^2/6)] i_1 = 1.
1033 // Finally, sin(r) and cos(r) are approximated by polynomials of
1034 // moderate lengths.
1036 // sin(r) = r + PP_1_hi r^3 + PP_1_lo r^3 +
1037 // PP_2 r^5 + ... + PP_8 r^17
1039 // cos(r) = 1 + QQ_1 r^2 + QQ_2 r^4 + ... + QQ_8 r^16
1041 // where PP_1_hi is only about 16 bits long and QQ_1 is -1/2.
1042 // The crux in accurate computation is to calculate
1044 // r + PP_1_hi r^3 or 1 + QQ_1 r^2
1046 // accurately as two pieces: U_hi and U_lo. The way to achieve this
1047 // is to obtain r_hi as a 10 sig. bit number that approximates r to
1048 // roughly 8 bits or so of accuracy. (One convenient way is
1050 // r_hi := frcpa( frcpa( r ) ).)
1054 // r + PP_1_hi r^3 = r + PP_1_hi r_hi^3 +
1055 // PP_1_hi (r^3 - r_hi^3)
1056 // = [r + PP_1_hi r_hi^3] +
1057 // [PP_1_hi (r - r_hi)
1058 // (r^2 + r_hi r + r_hi^2) ]
1061 // Since r_hi is only 10 bit long and PP_1_hi is only 16 bit long,
1062 // PP_1_hi * r_hi^3 is only at most 46 bit long and thus computed
1063 // exactly. Furthermore, r and PP_1_hi r_hi^3 are of opposite sign
1064 // and that there is no more than 8 bit shift off between r and
1065 // PP_1_hi * r_hi^3. Hence the sum, U_hi, is representable and thus
1066 // calculated without any error. Finally, the fact that
1068 // |U_lo| <= 2^(-8) |U_hi|
1070 // says that U_hi + U_lo is approximating r + PP_1_hi r^3 to roughly
1071 // 8 extra bits of accuracy.
1075 // 1 + QQ_1 r^2 = [1 + QQ_1 r_hi^2] +
1076 // [QQ_1 (r - r_hi)(r + r_hi)]
1079 // Summarizing, we calculate r_hi = frcpa( frcpa( r ) ).
1083 // U_hi := r + PP_1_hi * r_hi^3
1084 // U_lo := PP_1_hi * (r - r_hi) * (r^2 + r*r_hi + r_hi^2)
1085 // poly := PP_1_lo r^3 + PP_2 r^5 + ... + PP_8 r^17
1086 // correction := c * ( 1 + C_1 r^2 )
1090 // U_hi := 1 + QQ_1 * r_hi * r_hi
1091 // U_lo := QQ_1 * (r - r_hi) * (r + r_hi)
1092 // poly := QQ_2 * r^4 + QQ_3 * r^6 + ... + QQ_8 r^16
1093 // correction := -c * r * (1 + S_1 * r^2)
1099 // V := poly + ( U_lo + correction )
1101 // / U_hi + V if i_0 = 0
1103 // \ (-U_hi) - V if i_0 = 1
1105 // It is important that in the last step, negation of U_hi is
1106 // performed prior to the subtraction which is to be performed in
1107 // the user-set rounding mode.
1110 // Algorithmic Description
1111 // =======================
1113 // The argument reduction algorithm is tightly integrated into FSIN
1114 // and FCOS which share the same code. The following is complete and
1115 // self-contained. The argument reduction description given
1116 // previously is repeated below.
1119 // Step 0. Initialization.
1121 // If FSIN is invoked, set N_inc := 0; else if FCOS is invoked,
1124 // Step 1. Check for exceptional and special cases.
1126 // * If Arg is +-0, +-inf, NaN, NaT, go to Step 10 for special
1128 // * If |Arg| < 2^24, go to Step 2 for reduction of moderate
1129 // arguments. This is the most likely case.
1130 // * If |Arg| < 2^63, go to Step 8 for pre-reduction of large
1132 // * If |Arg| >= 2^63, go to Step 10 for special handling.
1134 // Step 2. Reduction of moderate arguments.
1136 // If |Arg| < pi/4 ...quick branch
1137 // N_fix := N_inc (integer)
1140 // Branch to Step 4, Case_1_complete
1141 // Else ...cf. argument reduction
1142 // N := Arg * two_by_PI (fp)
1143 // N_fix := fcvt.fx( N ) (int)
1144 // N := fcvt.xf( N_fix )
1145 // N_fix := N_fix + N_inc
1146 // s := Arg - N * P_1 (first piece of pi/2)
1147 // w := -N * P_2 (second piece of pi/2)
1149 // If |s| >= 2^(-33)
1150 // go to Step 3, Case_1_reduce
1152 // go to Step 7, Case_2_reduce
1156 // Step 3. Case_1_reduce.
1159 // c := (s - r) + w ...observe order
1161 // Step 4. Case_1_complete
1163 // ...At this point, the reduced argument alpha is
1164 // ...accurately represented as r + c.
1165 // If |r| < 2^(-3), go to Step 6, small_r.
1167 // Step 5. Normal_r.
1169 // Let [i_0 i_1] by the 2 lsb of N_fix.
1171 // r_hi := frcpa( frcpa( r ) )
1175 // poly := r*FR_rsq*(PP_1_lo + FR_rsq*(PP_2 + ... FR_rsq*PP_8))
1176 // U_hi := r + PP_1_hi*r_hi*r_hi*r_hi ...any order
1177 // U_lo := PP_1_hi*r_lo*(r*r + r*r_hi + r_hi*r_hi)
1178 // correction := c + c*C_1*FR_rsq ...any order
1180 // poly := FR_rsq*FR_rsq*(QQ_2 + FR_rsq*(QQ_3 + ... + FR_rsq*QQ_8))
1181 // U_hi := 1 + QQ_1 * r_hi * r_hi ...any order
1182 // U_lo := QQ_1 * r_lo * (r + r_hi)
1183 // correction := -c*(r + S_1*FR_rsq*r) ...any order
1186 // V := poly + (U_lo + correction) ...observe order
1188 // result := (i_0 == 0? 1.0 : -1.0)
1190 // Last instruction in user-set rounding mode
1192 // result := (i_0 == 0? result*U_hi + V :
1199 // ...Use flush to zero mode without causing exception
1200 // Let [i_0 i_1] be the two lsb of N_fix.
1205 // z := FR_rsq*FR_rsq; z := FR_rsq*z *r
1206 // poly_lo := S_3 + FR_rsq*(S_4 + FR_rsq*S_5)
1207 // poly_hi := r*FR_rsq*(S_1 + FR_rsq*S_2)
1211 // z := FR_rsq*FR_rsq; z := FR_rsq*z
1212 // poly_lo := C_3 + FR_rsq*(C_4 + FR_rsq*C_5)
1213 // poly_hi := FR_rsq*(C_1 + FR_rsq*C_2)
1214 // correction := -c*r
1218 // poly := poly_hi + (z * poly_lo + correction)
1220 // If i_0 = 1, result := -result
1222 // Last operation. Perform in user-set rounding mode
1224 // result := (i_0 == 0? result + poly :
1228 // Step 7. Case_2_reduce.
1230 // ...Refer to the write up for argument reduction for
1231 // ...rationale. The reduction algorithm below is taken from
1232 // ...argument reduction description and integrated this.
1235 // U_1 := N*P_2 + w ...FMA
1236 // U_2 := (N*P_2 - U_1) + w ...2 FMA
1237 // ...U_1 + U_2 is N*(P_2+P_3) accurately
1240 // c := ( (s - r) - U_1 ) - U_2
1242 // ...The mathematical sum r + c approximates the reduced
1243 // ...argument accurately. Note that although compared to
1244 // ...Case 1, this case requires much more work to reduce
1245 // ...the argument, the subsequent calculation needed for
1246 // ...any of the trigonometric function is very little because
1247 // ...|alpha| < 1.01*2^(-33) and thus two terms of the
1248 // ...Taylor series expansion suffices.
1251 // poly := c + S_1 * r * r * r ...any order
1258 // If i_0 = 1, result := -result
1260 // Last operation. Perform in user-set rounding mode
1262 // result := (i_0 == 0? result + poly :
1268 // Step 8. Pre-reduction of large arguments.
1270 // ...Again, the following reduction procedure was described
1271 // ...in the separate write up for argument reduction, which
1272 // ...is tightly integrated here.
1274 // N_0 := Arg * Inv_P_0
1275 // N_0_fix := fcvt.fx( N_0 )
1276 // N_0 := fcvt.xf( N_0_fix)
1278 // Arg' := Arg - N_0 * P_0
1280 // N := Arg' * two_by_PI
1281 // N_fix := fcvt.fx( N )
1282 // N := fcvt.xf( N_fix )
1283 // N_fix := N_fix + N_inc
1285 // s := Arg' - N * P_1
1288 // If |s| >= 2^(-14)
1294 // Step 9. Case_4_reduce.
1296 // ...first obtain N_0*d_1 and -N*P_2 accurately
1297 // U_hi := N_0 * d_1 V_hi := -N*P_2
1298 // U_lo := N_0 * d_1 - U_hi V_lo := -N*P_2 - U_hi ...FMAs
1300 // ...compute the contribution from N_0*d_1 and -N*P_3
1303 // t := U_lo + V_lo + w ...any order
1305 // ...at this point, the mathematical value
1306 // ...s + U_hi + V_hi + t approximates the true reduced argument
1307 // ...accurately. Just need to compute this accurately.
1309 // ...Calculate U_hi + V_hi accurately:
1311 // if |U_hi| >= |V_hi| then
1312 // a := (U_hi - A) + V_hi
1314 // a := (V_hi - A) + U_hi
1316 // ...order in computing "a" must be observed. This branch is
1317 // ...best implemented by predicates.
1318 // ...A + a is U_hi + V_hi accurately. Moreover, "a" is
1319 // ...much smaller than A: |a| <= (1/2)ulp(A).
1321 // ...Just need to calculate s + A + a + t
1322 // C_hi := s + A t := t + a
1323 // C_lo := (s - C_hi) + A
1326 // ...Final steps for reduction
1328 // c := (C_hi - r) + C_lo
1330 // ...At this point, we have r and c
1331 // ...And all we need is a couple of terms of the corresponding
1332 // ...Taylor series.
1335 // poly := c + r*FR_rsq*(S_1 + FR_rsq*S_2)
1338 // poly := FR_rsq*(C_1 + FR_rsq*C_2)
1342 // If i_0 = 1, result := -result
1344 // Last operation. Perform in user-set rounding mode
1346 // result := (i_0 == 0? result + poly :
1350 // Large Arguments: For arguments above 2**63, a Payne-Hanek
1351 // style argument reduction is used and pi_by_2 reduce is called.
1363 ASM_TYPE_DIRECTIVE(FSINCOS_CONSTANTS,@object)
1364 data4 0x4B800000, 0xCB800000, 0x00000000,0x00000000 // two**24, -two**24
1365 data4 0x4E44152A, 0xA2F9836E, 0x00003FFE,0x00000000 // Inv_pi_by_2
1366 data4 0xCE81B9F1, 0xC84D32B0, 0x00004016,0x00000000 // P_0
1367 data4 0x2168C235, 0xC90FDAA2, 0x00003FFF,0x00000000 // P_1
1368 data4 0xFC8F8CBB, 0xECE675D1, 0x0000BFBD,0x00000000 // P_2
1369 data4 0xACC19C60, 0xB7ED8FBB, 0x0000BF7C,0x00000000 // P_3
1370 data4 0x5F000000, 0xDF000000, 0x00000000,0x00000000 // two_to_63, -two_to_63
1371 data4 0x6EC6B45A, 0xA397E504, 0x00003FE7,0x00000000 // Inv_P_0
1372 data4 0xDBD171A1, 0x8D848E89, 0x0000BFBF,0x00000000 // d_1
1373 data4 0x18A66F8E, 0xD5394C36, 0x0000BF7C,0x00000000 // d_2
1374 data4 0x2168C234, 0xC90FDAA2, 0x00003FFE,0x00000000 // pi_by_4
1375 data4 0x2168C234, 0xC90FDAA2, 0x0000BFFE,0x00000000 // neg_pi_by_4
1376 data4 0x3E000000, 0xBE000000, 0x00000000,0x00000000 // two**-3, -two**-3
1377 data4 0x2F000000, 0xAF000000, 0x9E000000,0x00000000 // two**-33, -two**-33, -two**-67
1378 data4 0xA21C0BC9, 0xCC8ABEBC, 0x00003FCE,0x00000000 // PP_8
1379 data4 0x720221DA, 0xD7468A05, 0x0000BFD6,0x00000000 // PP_7
1380 data4 0x640AD517, 0xB092382F, 0x00003FDE,0x00000000 // PP_6
1381 data4 0xD1EB75A4, 0xD7322B47, 0x0000BFE5,0x00000000 // PP_5
1382 data4 0xFFFFFFFE, 0xFFFFFFFF, 0x0000BFFD,0x00000000 // C_1
1383 data4 0x00000000, 0xAAAA0000, 0x0000BFFC,0x00000000 // PP_1_hi
1384 data4 0xBAF69EEA, 0xB8EF1D2A, 0x00003FEC,0x00000000 // PP_4
1385 data4 0x0D03BB69, 0xD00D00D0, 0x0000BFF2,0x00000000 // PP_3
1386 data4 0x88888962, 0x88888888, 0x00003FF8,0x00000000 // PP_2
1387 data4 0xAAAB0000, 0xAAAAAAAA, 0x0000BFEC,0x00000000 // PP_1_lo
1388 data4 0xC2B0FE52, 0xD56232EF, 0x00003FD2,0x00000000 // QQ_8
1389 data4 0x2B48DCA6, 0xC9C99ABA, 0x0000BFDA,0x00000000 // QQ_7
1390 data4 0x9C716658, 0x8F76C650, 0x00003FE2,0x00000000 // QQ_6
1391 data4 0xFDA8D0FC, 0x93F27DBA, 0x0000BFE9,0x00000000 // QQ_5
1392 data4 0xAAAAAAAA, 0xAAAAAAAA, 0x0000BFFC,0x00000000 // S_1
1393 data4 0x00000000, 0x80000000, 0x0000BFFE,0x00000000 // QQ_1
1394 data4 0x0C6E5041, 0xD00D00D0, 0x00003FEF,0x00000000 // QQ_4
1395 data4 0x0B607F60, 0xB60B60B6, 0x0000BFF5,0x00000000 // QQ_3
1396 data4 0xAAAAAA9B, 0xAAAAAAAA, 0x00003FFA,0x00000000 // QQ_2
1397 data4 0xFFFFFFFE, 0xFFFFFFFF, 0x0000BFFD,0x00000000 // C_1
1398 data4 0xAAAA719F, 0xAAAAAAAA, 0x00003FFA,0x00000000 // C_2
1399 data4 0x0356F994, 0xB60B60B6, 0x0000BFF5,0x00000000 // C_3
1400 data4 0xB2385EA9, 0xD00CFFD5, 0x00003FEF,0x00000000 // C_4
1401 data4 0x292A14CD, 0x93E4BD18, 0x0000BFE9,0x00000000 // C_5
1402 data4 0xAAAAAAAA, 0xAAAAAAAA, 0x0000BFFC,0x00000000 // S_1
1403 data4 0x888868DB, 0x88888888, 0x00003FF8,0x00000000 // S_2
1404 data4 0x055EFD4B, 0xD00D00D0, 0x0000BFF2,0x00000000 // S_3
1405 data4 0x839730B9, 0xB8EF1C5D, 0x00003FEC,0x00000000 // S_4
1406 data4 0xE5B3F492, 0xD71EA3A4, 0x0000BFE5,0x00000000 // S_5
1407 data4 0x38800000, 0xB8800000, 0x00000000 // two**-14, -two**-14
1408 ASM_SIZE_DIRECTIVE(FSINCOS_CONSTANTS)
1411 FR_Neg_Two_to_M3 = f32
1417 FR_Neg_Two_to_24 = f36
1418 FR_Neg_Pi_by_4 = f36
1419 FR_Neg_Two_to_M14 = f37
1420 FR_Neg_Two_to_M33 = f38
1421 FR_Neg_Two_to_M67 = f39
1422 FR_Inv_pi_by_2 = f40
1485 FR_Neg_Two_to_63 = f94
1493 GR_Table_Base1 = r33
1504 .proc __libm_sin_double_dbx#
1506 __libm_sin_double_dbx:
1509 alloc GR_Table_Base = ar.pfs,0,12,2,0
1510 movl GR_Sin_or_Cos = 0x0 ;;
1515 addl GR_Table_Base = @ltoff(FSINCOS_CONSTANTS#), gp
1521 ld8 GR_Table_Base = [GR_Table_Base]
1531 br.cond.sptk L(SINCOS_CONTINUE) ;;
1534 .endp __libm_sin_double_dbx#
1535 ASM_SIZE_DIRECTIVE(__libm_sin_double_dbx)
1538 .proc __libm_cos_double_dbx#
1539 __libm_cos_double_dbx:
1542 alloc GR_Table_Base= ar.pfs,0,12,2,0
1543 movl GR_Sin_or_Cos = 0x1 ;;
1548 addl GR_Table_Base = @ltoff(FSINCOS_CONSTANTS#), gp
1554 ld8 GR_Table_Base = [GR_Table_Base]
1561 // Load Table Address
1566 add GR_Table_Base1 = 96, GR_Table_Base
1567 ldfs FR_Two_to_24 = [GR_Table_Base], 4
1575 // Load 2**24, load 2**63.
1577 ldfs FR_Neg_Two_to_24 = [GR_Table_Base], 12
1582 ldfs FR_Two_to_63 = [GR_Table_Base1], 4
1584 // Check for unnormals - unsupported operands. We do not want
1585 // to generate denormal exception
1586 // Check for NatVals, QNaNs, SNaNs, +/-Infs
1587 // Check for EM unsupporteds
1590 fclass.m.unc p6, p8 = FR_Input_X, 0x1E3
1596 fclass.nm.unc p8, p0 = FR_Input_X, 0x1FF
1597 // GR_Sin_or_Cos denotes
1602 ldfs FR_Neg_Two_to_63 = [GR_Table_Base1], 12
1603 fclass.m.unc p10, p0 = FR_Input_X, 0x007
1604 (p6) br.cond.spnt L(SINCOS_SPECIAL) ;;
1610 (p8) br.cond.spnt L(SINCOS_SPECIAL) ;;
1617 // Branch if +/- NaN, Inf.
1618 // Load -2**24, load -2**63.
1620 (p10) br.cond.spnt L(SINCOS_ZERO) ;;
1624 ldfe FR_Inv_pi_by_2 = [GR_Table_Base], 16
1625 ldfe FR_Inv_P_0 = [GR_Table_Base1], 16
1631 ldfe FR_d_1 = [GR_Table_Base1], 16
1635 // Raise possible denormal operand flag with useful fcmp
1637 // Load Inv_P_0 for pre-reduction
1642 ldfe FR_P_0 = [GR_Table_Base], 16
1643 ldfe FR_d_2 = [GR_Table_Base1], 16
1654 ldfe FR_P_1 = [GR_Table_Base], 16 ;;
1660 ldfe FR_P_2 = [GR_Table_Base], 16
1666 ldfe FR_P_3 = [GR_Table_Base], 16
1667 fcmp.le.unc.s1 p7, p8 = FR_Input_X, FR_Neg_Two_to_24
1673 // Branch if +/- zero.
1674 // Decide about the paths to take:
1675 // If -2**24 < FR_Input_X < 2**24 - CASE 1 OR 2
1676 // OTHERWISE - CASE 3 OR 4
1678 fcmp.le.unc.s0 p10, p11 = FR_Input_X, FR_Neg_Two_to_63
1684 (p8) fcmp.ge.s1 p7, p0 = FR_Input_X, FR_Two_to_24
1689 ldfe FR_Pi_by_4 = [GR_Table_Base1], 16
1690 (p11) fcmp.ge.s1 p10, p0 = FR_Input_X, FR_Two_to_63
1695 ldfe FR_Neg_Pi_by_4 = [GR_Table_Base1], 16 ;;
1696 ldfs FR_Two_to_M3 = [GR_Table_Base1], 4
1701 ldfs FR_Neg_Two_to_M3 = [GR_Table_Base1], 12
1711 (p10) br.cond.spnt L(SINCOS_ARG_TOO_LARGE) ;;
1718 // Branch out if x >= 2**63. Use Payne-Hanek Reduction
1720 (p7) br.cond.spnt L(SINCOS_LARGER_ARG) ;;
1726 // Branch if Arg <= -2**24 or Arg >= 2**24 and use pre-reduction.
1728 fma.s1 FR_N_float = FR_Input_X, FR_Inv_pi_by_2, f0
1734 fcmp.lt.unc.s1 p6, p7 = FR_Input_X, FR_Pi_by_4
1741 // Select the case when |Arg| < pi/4
1742 // Else Select the case when |Arg| >= pi/4
1744 fcvt.fx.s1 FR_N_fix = FR_N_float
1752 // Check if Arg < pi/4
1754 (p6) fcmp.gt.s1 p6, p7 = FR_Input_X, FR_Neg_Pi_by_4
1758 // Case 2: Convert integer N_fix back to normalized floating-point value.
1759 // Case 1: p8 is only affected when p6 is set
1763 (p7) ldfs FR_Two_to_M33 = [GR_Table_Base1], 4
1765 // Grab the integer part of N and call it N_fix
1767 (p6) fmerge.se FR_r = FR_Input_X, FR_Input_X
1768 // If |x| < pi/4, r = x and c = 0
1769 // lf |x| < pi/4, is x < 2**(-3).
1772 (p6) mov GR_N_Inc = GR_Sin_or_Cos ;;
1777 (p7) ldfs FR_Neg_Two_to_M33 = [GR_Table_Base1], 4
1778 (p6) fmerge.se FR_c = f0, f0
1783 (p6) fcmp.lt.unc.s1 p8, p9 = FR_Input_X, FR_Two_to_M3
1790 // lf |x| < pi/4, is -2**(-3)< x < 2**(-3) - set p8.
1792 // Create the right N for |x| < pi/4 and otherwise
1793 // Case 2: Place integer part of N in GP register
1795 (p7) fcvt.xf FR_N_float = FR_N_fix
1801 (p7) getf.sig GR_N_Inc = FR_N_fix
1802 (p8) fcmp.gt.s1 p8, p0 = FR_Input_X, FR_Neg_Two_to_M3 ;;
1809 // Load 2**(-33), -2**(-33)
1811 (p8) br.cond.spnt L(SINCOS_SMALL_R) ;;
1817 (p6) br.cond.sptk L(SINCOS_NORMAL_R) ;;
1820 // if |x| < pi/4, branch based on |x| < 2**(-3) or otherwise.
1823 // In this branch, |x| >= pi/4.
1827 ldfs FR_Neg_Two_to_M67 = [GR_Table_Base1], 8
1831 fnma.s1 FR_s = FR_N_float, FR_P_1, FR_Input_X
1834 // s = -N * P_1 + Arg
1836 add GR_N_Inc = GR_N_Inc, GR_Sin_or_Cos
1841 fma.s1 FR_w = FR_N_float, FR_P_2, f0
1848 // Adjust N_fix by N_inc to determine whether sine or
1849 // cosine is being calculated
1851 fcmp.lt.unc.s1 p7, p6 = FR_s, FR_Two_to_M33
1857 (p7) fcmp.gt.s1 p7, p6 = FR_s, FR_Neg_Two_to_M33
1863 // Remember x >= pi/4.
1864 // Is s <= -2**(-33) or s >= 2**(-33) (p6)
1865 // or -2**(-33) < s < 2**(-33) (p7)
1866 (p6) fms.s1 FR_r = FR_s, f1, FR_w
1872 (p7) fma.s1 FR_w = FR_N_float, FR_P_3, f0
1878 (p7) fma.s1 FR_U_1 = FR_N_float, FR_P_2, FR_w
1884 (p6) fms.s1 FR_c = FR_s, f1, FR_r
1891 // For big s: r = s - w: No futher reduction is necessary
1892 // For small s: w = N * P_3 (change sign) More reduction
1894 (p6) fcmp.lt.unc.s1 p8, p9 = FR_r, FR_Two_to_M3
1900 (p8) fcmp.gt.s1 p8, p9 = FR_r, FR_Neg_Two_to_M3
1906 (p7) fms.s1 FR_r = FR_s, f1, FR_U_1
1913 // For big s: Is |r| < 2**(-3)?
1914 // For big s: c = S - r
1915 // For small s: U_1 = N * P_2 + w
1917 // If p8 is set, prepare to branch to Small_R.
1918 // If p9 is set, prepare to branch to Normal_R.
1919 // For big s, r is complete here.
1921 (p6) fms.s1 FR_c = FR_c, f1, FR_w
1923 // For big s: c = c + w (w has not been negated.)
1924 // For small s: r = S - U_1
1926 (p8) br.cond.spnt L(SINCOS_SMALL_R) ;;
1932 (p9) br.cond.sptk L(SINCOS_NORMAL_R) ;;
1936 (p7) add GR_Table_Base1 = 224, GR_Table_Base1
1938 // Branch to SINCOS_SMALL_R or SINCOS_NORMAL_R
1940 (p7) fms.s1 FR_U_2 = FR_N_float, FR_P_2, FR_U_1
1946 (p7) extr.u GR_i_1 = GR_N_Inc, 0, 1
1952 // Get [i_0,i_1] - two lsb of N_fix_gr.
1953 // Do dummy fmpy so inexact is always set.
1955 (p7) cmp.eq.unc p9, p10 = 0x0, GR_i_1
1956 (p7) extr.u GR_i_0 = GR_N_Inc, 1, 1 ;;
1959 // For small s: U_2 = N * P_2 - U_1
1960 // S_1 stored constant - grab the one stored with the
1965 (p7) ldfe FR_S_1 = [GR_Table_Base1], 16
1967 // Check if i_1 and i_0 != 0
1969 (p10) fma.s1 FR_poly = f0, f1, FR_Neg_Two_to_M67
1970 (p7) cmp.eq.unc p11, p12 = 0x0, GR_i_0 ;;
1975 (p7) fms.s1 FR_s = FR_s, f1, FR_r
1986 (p7) fma.s1 FR_rsq = FR_r, FR_r, f0
1992 (p7) fma.s1 FR_U_2 = FR_U_2, f1, FR_w
1998 (p7) fmerge.se FR_Input_X = FR_r, FR_r
2004 (p10) fma.s1 FR_Input_X = f0, f1, f1
2012 // Save r as the result.
2014 (p7) fms.s1 FR_c = FR_s, f1, FR_U_1
2021 // if ( i_1 ==0) poly = c + S_1*r*r*r
2024 (p12) fnma.s1 FR_Input_X = FR_Input_X, f1, f0
2030 (p7) fma.s1 FR_r = FR_S_1, FR_r, f0
2036 (p7) fma.d.s0 FR_S_1 = FR_S_1, FR_S_1, f0
2043 // If i_1 != 0, poly = 2**(-67)
2045 (p7) fms.s1 FR_c = FR_c, f1, FR_U_2
2054 (p9) fma.s1 FR_poly = FR_r, FR_rsq, FR_c
2061 // i_0 != 0, so Result = -Result
2063 (p11) fma.d.s0 FR_Input_X = FR_Input_X, f1, FR_poly
2069 (p12) fms.d.s0 FR_Input_X = FR_Input_X, f1, FR_poly
2071 // if (i_0 == 0), Result = Result + poly
2072 // else Result = Result - poly
2076 L(SINCOS_LARGER_ARG):
2080 fma.s1 FR_N_0 = FR_Input_X, FR_Inv_P_0, f0
2085 // This path for argument > 2*24
2086 // Adjust table_ptr1 to beginning of table.
2091 addl GR_Table_Base = @ltoff(FSINCOS_CONSTANTS#), gp
2097 ld8 GR_Table_Base = [GR_Table_Base]
2106 // N_0 = Arg * Inv_P_0
2110 add GR_Table_Base = 688, GR_Table_Base ;;
2111 ldfs FR_Two_to_M14 = [GR_Table_Base], 4
2116 ldfs FR_Neg_Two_to_M14 = [GR_Table_Base], 0
2124 // Load values 2**(-14) and -2**(-14)
2126 fcvt.fx.s1 FR_N_0_fix = FR_N_0
2133 // N_0_fix = integer part of N_0
2135 fcvt.xf FR_N_0 = FR_N_0_fix
2142 // Make N_0 the integer part
2144 fnma.s1 FR_ArgPrime = FR_N_0, FR_P_0, FR_Input_X
2150 fma.s1 FR_w = FR_N_0, FR_d_1, f0
2157 // Arg' = -N_0 * P_0 + Arg
2160 fma.s1 FR_N_float = FR_ArgPrime, FR_Inv_pi_by_2, f0
2169 fcvt.fx.s1 FR_N_fix = FR_N_float
2176 // N_fix is the integer part
2178 fcvt.xf FR_N_float = FR_N_fix
2183 getf.sig GR_N_Inc = FR_N_fix
2191 add GR_N_Inc = GR_N_Inc, GR_Sin_or_Cos ;;
2197 // N is the integer part of the reduced-reduced argument.
2198 // Put the integer in a GP register
2200 fnma.s1 FR_s = FR_N_float, FR_P_1, FR_ArgPrime
2206 fnma.s1 FR_w = FR_N_float, FR_P_2, FR_w
2213 // s = -N*P_1 + Arg'
2215 // N_fix_gr = N_fix_gr + N_inc
2217 fcmp.lt.unc.s1 p9, p8 = FR_s, FR_Two_to_M14
2223 (p9) fcmp.gt.s1 p9, p8 = FR_s, FR_Neg_Two_to_M14
2230 // For |s| > 2**(-14) r = S + w (r complete)
2231 // Else U_hi = N_0 * d_1
2233 (p9) fma.s1 FR_V_hi = FR_N_float, FR_P_2, f0
2239 (p9) fma.s1 FR_U_hi = FR_N_0, FR_d_1, f0
2246 // Either S <= -2**(-14) or S >= 2**(-14)
2247 // or -2**(-14) < s < 2**(-14)
2249 (p8) fma.s1 FR_r = FR_s, f1, FR_w
2255 (p9) fma.s1 FR_w = FR_N_float, FR_P_3, f0
2262 // We need abs of both U_hi and V_hi - don't
2263 // worry about switched sign of V_hi.
2265 (p9) fms.s1 FR_A = FR_U_hi, f1, FR_V_hi
2272 // Big s: finish up c = (S - r) + w (c complete)
2273 // Case 4: A = U_hi + V_hi
2274 // Note: Worry about switched sign of V_hi, so subtract instead of add.
2276 (p9) fnma.s1 FR_V_lo = FR_N_float, FR_P_2, FR_V_hi
2282 (p9) fms.s1 FR_U_lo = FR_N_0, FR_d_1, FR_U_hi
2288 (p9) fmerge.s FR_V_hiabs = f0, FR_V_hi
2294 // For big s: c = S - r
2295 // For small s do more work: U_lo = N_0 * d_1 - U_hi
2297 (p9) fmerge.s FR_U_hiabs = f0, FR_U_hi
2304 // For big s: Is |r| < 2**(-3)
2305 // For big s: if p12 set, prepare to branch to Small_R.
2306 // For big s: If p13 set, prepare to branch to Normal_R.
2308 (p8) fms.s1 FR_c = FR_s, f1, FR_r
2315 // For small S: V_hi = N * P_2
2317 // Note the product does not include the (-) as in the writeup
2318 // so (-) missing for V_hi and w.
2320 (p8) fcmp.lt.unc.s1 p12, p13 = FR_r, FR_Two_to_M3
2326 (p12) fcmp.gt.s1 p12, p13 = FR_r, FR_Neg_Two_to_M3
2332 (p8) fma.s1 FR_c = FR_c, f1, FR_w
2338 (p9) fms.s1 FR_w = FR_N_0, FR_d_2, FR_w
2339 (p12) br.cond.spnt L(SINCOS_SMALL_R) ;;
2345 (p13) br.cond.sptk L(SINCOS_NORMAL_R) ;;
2351 // Big s: Vector off when |r| < 2**(-3). Recall that p8 will be true.
2352 // The remaining stuff is for Case 4.
2353 // Small s: V_lo = N * P_2 + U_hi (U_hi is in place of V_hi in writeup)
2354 // Note: the (-) is still missing for V_lo.
2355 // Small s: w = w + N_0 * d_2
2356 // Note: the (-) is now incorporated in w.
2358 (p9) fcmp.ge.unc.s1 p10, p11 = FR_U_hiabs, FR_V_hiabs
2359 extr.u GR_i_1 = GR_N_Inc, 0, 1 ;;
2367 (p9) fma.s1 FR_t = FR_U_lo, f1, FR_V_lo
2368 extr.u GR_i_0 = GR_N_Inc, 1, 1 ;;
2377 (p10) fms.s1 FR_a = FR_U_hi, f1, FR_A
2383 (p11) fma.s1 FR_a = FR_V_hi, f1, FR_A
2390 addl GR_Table_Base = @ltoff(FSINCOS_CONSTANTS#), gp
2396 ld8 GR_Table_Base = [GR_Table_Base]
2404 add GR_Table_Base = 528, GR_Table_Base
2406 // Is U_hiabs >= V_hiabs?
2408 (p9) fma.s1 FR_C_hi = FR_s, f1, FR_A
2413 ldfe FR_C_1 = [GR_Table_Base], 16 ;;
2414 ldfe FR_C_2 = [GR_Table_Base], 64
2421 // c = c + C_lo finished.
2424 ldfe FR_S_1 = [GR_Table_Base], 16
2428 fma.s1 FR_t = FR_t, f1, FR_w ;;
2431 // r and c have been computed.
2432 // Make sure ftz mode is set - should be automatic when using wre
2434 // Get [i_0,i_1] - two lsb of N_fix.
2439 ldfe FR_S_2 = [GR_Table_Base], 64
2443 (p10) fms.s1 FR_a = FR_a, f1, FR_V_hi
2444 cmp.eq.unc p9, p10 = 0x0, GR_i_0
2450 // For larger u than v: a = U_hi - A
2451 // Else a = V_hi - A (do an add to account for missing (-) on V_hi
2453 fms.s1 FR_C_lo = FR_s, f1, FR_C_hi
2459 (p11) fms.s1 FR_a = FR_U_hi, f1, FR_a
2460 cmp.eq.unc p11, p12 = 0x0, GR_i_1
2466 // If u > v: a = (U_hi - A) + V_hi
2467 // Else a = (V_hi - A) + U_hi
2468 // In each case account for negative missing from V_hi.
2470 fma.s1 FR_C_lo = FR_C_lo, f1, FR_A
2477 // C_lo = (S - C_hi) + A
2479 fma.s1 FR_t = FR_t, f1, FR_a
2488 fma.s1 FR_C_lo = FR_C_lo, f1, FR_t
2496 // Adjust Table_Base to beginning of table
2498 fma.s1 FR_r = FR_C_hi, f1, FR_C_lo
2507 fma.s1 FR_rsq = FR_r, FR_r, f0
2514 // Table_Base points to C_1
2517 fms.s1 FR_c = FR_C_hi, f1, FR_r
2524 // if i_1 ==0: poly = S_2 * FR_rsq + S_1
2525 // else poly = C_2 * FR_rsq + C_1
2527 (p11) fma.s1 FR_Input_X = f0, f1, FR_r
2533 (p12) fma.s1 FR_Input_X = f0, f1, f1
2540 // Compute r_cube = FR_rsq * r
2542 (p11) fma.s1 FR_poly = FR_rsq, FR_S_2, FR_S_1
2548 (p12) fma.s1 FR_poly = FR_rsq, FR_C_2, FR_C_1
2555 // Compute FR_rsq = r * r
2558 fma.s1 FR_r_cubed = FR_rsq, FR_r, f0
2568 fma.s1 FR_c = FR_c, f1, FR_C_lo
2575 // if i_1 ==0: poly = r_cube * poly + c
2576 // else poly = FR_rsq * poly
2578 (p10) fms.s1 FR_Input_X = f0, f1, FR_Input_X
2585 // if i_1 ==0: Result = r
2586 // else Result = 1.0
2588 (p11) fma.s1 FR_poly = FR_r_cubed, FR_poly, FR_c
2594 (p12) fma.s1 FR_poly = FR_rsq, FR_poly, f0
2601 // if i_0 !=0: Result = -Result
2603 (p9) fma.d.s0 FR_Input_X = FR_Input_X, f1, FR_poly
2609 (p10) fms.d.s0 FR_Input_X = FR_Input_X, f1, FR_poly
2611 // if i_0 == 0: Result = Result + poly
2612 // else Result = Result - poly
2620 extr.u GR_i_1 = GR_N_Inc, 0, 1 ;;
2623 // Compare both i_1 and i_0 with 0.
2624 // if i_1 == 0, set p9.
2625 // if i_0 == 0, set p11.
2627 cmp.eq.unc p9, p10 = 0x0, GR_i_1 ;;
2632 fma.s1 FR_rsq = FR_r, FR_r, f0
2633 extr.u GR_i_0 = GR_N_Inc, 1, 1 ;;
2641 (p10) fnma.s1 FR_c = FR_c, FR_r, f0
2642 cmp.eq.unc p11, p12 = 0x0, GR_i_0
2646 // ******************************************************************
2647 // ******************************************************************
2648 // ******************************************************************
2649 // r and c have been computed.
2650 // We know whether this is the sine or cosine routine.
2651 // Make sure ftz mode is set - should be automatic when using wre
2654 // Set table_ptr1 to beginning of constant table.
2655 // Get [i_0,i_1] - two lsb of N_fix_gr.
2660 addl GR_Table_Base = @ltoff(FSINCOS_CONSTANTS#), gp
2666 ld8 GR_Table_Base = [GR_Table_Base]
2674 // Set table_ptr1 to point to S_5.
2675 // Set table_ptr1 to point to C_5.
2676 // Compute FR_rsq = r * r
2680 (p9) add GR_Table_Base = 672, GR_Table_Base
2681 (p10) fmerge.s FR_r = f1, f1
2682 (p10) add GR_Table_Base = 592, GR_Table_Base ;;
2685 // Set table_ptr1 to point to S_5.
2686 // Set table_ptr1 to point to C_5.
2690 (p9) ldfe FR_S_5 = [GR_Table_Base], -16 ;;
2692 // if (i_1 == 0) load S_5
2693 // if (i_1 != 0) load C_5
2695 (p9) ldfe FR_S_4 = [GR_Table_Base], -16
2700 (p10) ldfe FR_C_5 = [GR_Table_Base], -16
2702 // Z = FR_rsq * FR_rsq
2704 (p9) ldfe FR_S_3 = [GR_Table_Base], -16
2706 // Compute FR_rsq = r * r
2707 // if (i_1 == 0) load S_4
2708 // if (i_1 != 0) load C_4
2710 fma.s1 FR_Z = FR_rsq, FR_rsq, f0 ;;
2713 // if (i_1 == 0) load S_3
2714 // if (i_1 != 0) load C_3
2718 (p9) ldfe FR_S_2 = [GR_Table_Base], -16 ;;
2720 // if (i_1 == 0) load S_2
2721 // if (i_1 != 0) load C_2
2723 (p9) ldfe FR_S_1 = [GR_Table_Base], -16
2728 (p10) ldfe FR_C_4 = [GR_Table_Base], -16 ;;
2729 (p10) ldfe FR_C_3 = [GR_Table_Base], -16
2734 (p10) ldfe FR_C_2 = [GR_Table_Base], -16 ;;
2735 (p10) ldfe FR_C_1 = [GR_Table_Base], -16
2743 // poly_lo = FR_rsq * C_5 + C_4
2744 // poly_hi = FR_rsq * C_2 + C_1
2746 (p9) fma.s1 FR_Z = FR_Z, FR_r, f0
2753 // if (i_1 == 0) load S_1
2754 // if (i_1 != 0) load C_1
2756 (p9) fma.s1 FR_poly_lo = FR_rsq, FR_S_5, FR_S_4
2764 // dummy fmpy's to flag inexact.
2766 (p9) fma.d.s0 FR_S_4 = FR_S_4, FR_S_4, f0
2773 // poly_lo = FR_rsq * poly_lo + C_3
2774 // poly_hi = FR_rsq * poly_hi
2776 fma.s1 FR_Z = FR_Z, FR_rsq, f0
2782 (p9) fma.s1 FR_poly_hi = FR_rsq, FR_S_2, FR_S_1
2790 // poly_lo = FR_rsq * S_5 + S_4
2791 // poly_hi = FR_rsq * S_2 + S_1
2793 (p10) fma.s1 FR_poly_lo = FR_rsq, FR_C_5, FR_C_4
2801 // Z = Z * r for only one of the small r cases - not there
2802 // in original implementation notes.
2804 (p9) fma.s1 FR_poly_lo = FR_rsq, FR_poly_lo, FR_S_3
2810 (p10) fma.s1 FR_poly_hi = FR_rsq, FR_C_2, FR_C_1
2816 (p10) fma.d.s0 FR_C_1 = FR_C_1, FR_C_1, f0
2822 (p9) fma.s1 FR_poly_hi = FR_poly_hi, FR_rsq, f0
2829 // poly_lo = FR_rsq * poly_lo + S_3
2830 // poly_hi = FR_rsq * poly_hi
2832 (p10) fma.s1 FR_poly_lo = FR_rsq, FR_poly_lo, FR_C_3
2838 (p10) fma.s1 FR_poly_hi = FR_poly_hi, FR_rsq, f0
2845 // if (i_1 == 0): dummy fmpy's to flag inexact
2848 (p9) fma.s1 FR_poly_hi = FR_r, FR_poly_hi, f0
2855 // poly_hi = r * poly_hi
2857 fma.s1 FR_poly = FR_Z, FR_poly_lo, FR_c
2863 (p12) fms.s1 FR_r = f0, f1, FR_r
2870 // poly_hi = Z * poly_lo + c
2871 // if i_0 == 1: r = -r
2873 fma.s1 FR_poly = FR_poly, f1, FR_poly_hi
2879 (p12) fms.d.s0 FR_Input_X = FR_r, f1, FR_poly
2886 // poly = poly + poly_hi
2888 (p11) fma.d.s0 FR_Input_X = FR_r, f1, FR_poly
2890 // if (i_0 == 0) Result = r + poly
2891 // if (i_0 != 0) Result = r - poly
2899 extr.u GR_i_1 = GR_N_Inc, 0, 1 ;;
2901 // Set table_ptr1 and table_ptr2 to base address of
2903 cmp.eq.unc p9, p10 = 0x0, GR_i_1 ;;
2908 fma.s1 FR_rsq = FR_r, FR_r, f0
2909 extr.u GR_i_0 = GR_N_Inc, 1, 1 ;;
2914 frcpa.s1 FR_r_hi, p6 = f1, FR_r
2915 cmp.eq.unc p11, p12 = 0x0, GR_i_0
2919 // ******************************************************************
2920 // ******************************************************************
2921 // ******************************************************************
2923 // r and c have been computed.
2924 // We known whether this is the sine or cosine routine.
2925 // Make sure ftz mode is set - should be automatic when using wre
2926 // Get [i_0,i_1] - two lsb of N_fix_gr alone.
2931 addl GR_Table_Base = @ltoff(FSINCOS_CONSTANTS#), gp
2937 ld8 GR_Table_Base = [GR_Table_Base]
2945 (p10) add GR_Table_Base = 384, GR_Table_Base
2946 (p12) fms.s1 FR_Input_X = f0, f1, f1
2947 (p9) add GR_Table_Base = 224, GR_Table_Base ;;
2952 (p10) ldfe FR_QQ_8 = [GR_Table_Base], 16
2954 // if (i_1==0) poly = poly * FR_rsq + PP_1_lo
2955 // else poly = FR_rsq * poly
2957 (p11) fma.s1 FR_Input_X = f0, f1, f1 ;;
2961 (p10) ldfe FR_QQ_7 = [GR_Table_Base], 16
2963 // Adjust table pointers based on i_0
2964 // Compute rsq = r * r
2966 (p9) ldfe FR_PP_8 = [GR_Table_Base], 16
2967 fma.s1 FR_r_cubed = FR_r, FR_rsq, f0 ;;
2971 (p9) ldfe FR_PP_7 = [GR_Table_Base], 16
2972 (p10) ldfe FR_QQ_6 = [GR_Table_Base], 16
2974 // Load PP_8 and QQ_8; PP_7 and QQ_7
2976 frcpa.s1 FR_r_hi, p6 = f1, FR_r_hi ;;
2979 // if (i_1==0) poly = PP_7 + FR_rsq * PP_8.
2980 // else poly = QQ_7 + FR_rsq * QQ_8.
2984 (p9) ldfe FR_PP_6 = [GR_Table_Base], 16
2985 (p10) ldfe FR_QQ_5 = [GR_Table_Base], 16
2990 (p9) ldfe FR_PP_5 = [GR_Table_Base], 16
2991 (p10) ldfe FR_S_1 = [GR_Table_Base], 16
2996 (p10) ldfe FR_QQ_1 = [GR_Table_Base], 16
2997 (p9) ldfe FR_C_1 = [GR_Table_Base], 16
3002 (p10) ldfe FR_QQ_4 = [GR_Table_Base], 16 ;;
3003 (p9) ldfe FR_PP_1 = [GR_Table_Base], 16
3008 (p10) ldfe FR_QQ_3 = [GR_Table_Base], 16
3010 // if (i_1=0) corr = corr + c*c
3011 // else corr = corr * c
3013 (p9) ldfe FR_PP_4 = [GR_Table_Base], 16
3014 (p10) fma.s1 FR_poly = FR_rsq, FR_QQ_8, FR_QQ_7 ;;
3017 // if (i_1=0) poly = rsq * poly + PP_5
3018 // else poly = rsq * poly + QQ_5
3019 // Load PP_4 or QQ_4
3023 (p9) ldfe FR_PP_3 = [GR_Table_Base], 16
3024 (p10) ldfe FR_QQ_2 = [GR_Table_Base], 16
3026 // r_hi = frcpa(frcpa(r)).
3027 // r_cube = r * FR_rsq.
3029 (p9) fma.s1 FR_poly = FR_rsq, FR_PP_8, FR_PP_7 ;;
3032 // Do dummy multiplies so inexact is always set.
3036 (p9) ldfe FR_PP_2 = [GR_Table_Base], 16
3040 (p9) fma.s1 FR_U_lo = FR_r_hi, FR_r_hi, f0
3046 (p9) ldfe FR_PP_1_lo = [GR_Table_Base], 16
3047 (p10) fma.s1 FR_corr = FR_S_1, FR_r_cubed, FR_r
3052 (p10) fma.s1 FR_poly = FR_rsq, FR_poly, FR_QQ_6
3059 // if (i_1=0) U_lo = r_hi * r_hi
3060 // else U_lo = r_hi + r
3062 (p9) fma.s1 FR_corr = FR_C_1, FR_rsq, f0
3069 // if (i_1=0) corr = C_1 * rsq
3070 // else corr = S_1 * r_cubed + r
3072 (p9) fma.s1 FR_poly = FR_rsq, FR_poly, FR_PP_6
3078 (p10) fma.s1 FR_U_lo = FR_r_hi, f1, FR_r
3085 // if (i_1=0) U_hi = r_hi + U_hi
3086 // else U_hi = QQ_1 * U_hi + 1
3088 (p9) fma.s1 FR_U_lo = FR_r, FR_r_hi, FR_U_lo
3095 // U_hi = r_hi * r_hi
3097 fms.s1 FR_r_lo = FR_r, f1, FR_r_hi
3104 // Load PP_1, PP_6, PP_5, and C_1
3105 // Load QQ_1, QQ_6, QQ_5, and S_1
3107 fma.s1 FR_U_hi = FR_r_hi, FR_r_hi, f0
3113 (p10) fma.s1 FR_poly = FR_rsq, FR_poly, FR_QQ_5
3119 (p10) fnma.s1 FR_corr = FR_corr, FR_c, f0
3126 // if (i_1=0) U_lo = r * r_hi + U_lo
3127 // else U_lo = r_lo * U_lo
3129 (p9) fma.s1 FR_corr = FR_corr, FR_c, FR_c
3135 (p9) fma.s1 FR_poly = FR_rsq, FR_poly, FR_PP_5
3142 // if (i_1 =0) U_hi = r + U_hi
3143 // if (i_1 =0) U_lo = r_lo * U_lo
3146 (p9) fma.d.s0 FR_PP_5 = FR_PP_5, FR_PP_4, f0
3152 (p9) fma.s1 FR_U_lo = FR_r, FR_r, FR_U_lo
3158 (p10) fma.s1 FR_U_lo = FR_r_lo, FR_U_lo, f0
3165 // if (i_1=0) poly = poly * rsq + PP_6
3166 // else poly = poly * rsq + QQ_6
3168 (p9) fma.s1 FR_U_hi = FR_r_hi, FR_U_hi, f0
3174 (p10) fma.s1 FR_poly = FR_rsq, FR_poly, FR_QQ_4
3180 (p10) fma.s1 FR_U_hi = FR_QQ_1, FR_U_hi, f1
3186 (p10) fma.d.s0 FR_QQ_5 = FR_QQ_5, FR_QQ_5, f0
3193 // if (i_1!=0) U_hi = PP_1 * U_hi
3194 // if (i_1!=0) U_lo = r * r + U_lo
3195 // Load PP_3 or QQ_3
3197 (p9) fma.s1 FR_poly = FR_rsq, FR_poly, FR_PP_4
3203 (p9) fma.s1 FR_U_lo = FR_r_lo, FR_U_lo, f0
3209 (p10) fma.s1 FR_U_lo = FR_QQ_1,FR_U_lo, f0
3215 (p9) fma.s1 FR_U_hi = FR_PP_1, FR_U_hi, f0
3221 (p10) fma.s1 FR_poly = FR_rsq, FR_poly, FR_QQ_3
3230 (p9) fma.s1 FR_poly = FR_rsq, FR_poly, FR_PP_3
3237 // if (i_1==0) poly = FR_rsq * poly + PP_3
3238 // else poly = FR_rsq * poly + QQ_3
3241 (p9) fma.s1 FR_U_lo = FR_PP_1, FR_U_lo, f0
3248 // if (i_1 =0) poly = poly * rsq + pp_r4
3249 // else poly = poly * rsq + qq_r4
3251 (p9) fma.s1 FR_U_hi = FR_r, f1, FR_U_hi
3257 (p10) fma.s1 FR_poly = FR_rsq, FR_poly, FR_QQ_2
3264 // if (i_1==0) U_lo = PP_1_hi * U_lo
3265 // else U_lo = QQ_1 * U_lo
3267 (p9) fma.s1 FR_poly = FR_rsq, FR_poly, FR_PP_2
3274 // if (i_0==0) Result = 1
3277 fma.s1 FR_V = FR_U_lo, f1, FR_corr
3283 (p10) fma.s1 FR_poly = FR_rsq, FR_poly, f0
3290 // if (i_1==0) poly = FR_rsq * poly + PP_2
3291 // else poly = FR_rsq * poly + QQ_2
3293 (p9) fma.s1 FR_poly = FR_rsq, FR_poly, FR_PP_1_lo
3299 (p10) fma.s1 FR_poly = FR_rsq, FR_poly, f0
3308 (p9) fma.s1 FR_poly = FR_r_cubed, FR_poly, f0
3315 // if (i_1==0) poly = r_cube * poly
3316 // else poly = FR_rsq * poly
3318 fma.s1 FR_V = FR_poly, f1, FR_V
3324 (p12) fms.d.s0 FR_Input_X = FR_Input_X, FR_U_hi, FR_V
3333 (p11) fma.d.s0 FR_Input_X = FR_Input_X, FR_U_hi, FR_V
3335 // if (i_0==0) Result = Result * U_hi + V
3336 // else Result = Result * U_hi - V
3342 // If cosine, FR_Input_X = 1
3343 // If sine, FR_Input_X = +/-Zero (Input FR_Input_X)
3344 // Results are exact, no exceptions
3349 cmp.eq.unc p6, p7 = 0x1, GR_Sin_or_Cos
3356 (p7) fmerge.s FR_Input_X = FR_Input_X, FR_Input_X
3362 (p6) fmerge.s FR_Input_X = f1, f1
3369 // Path for Arg = +/- QNaN, SNaN, Inf
3370 // Invalid can be raised. SNaNs
3376 fmpy.d.s0 FR_Input_X = FR_Input_X, f0
3379 .endp __libm_cos_double_dbx#
3380 ASM_SIZE_DIRECTIVE(__libm_cos_double_dbx#)
3385 // Call int pi_by_2_reduce(double* x, double *y)
3386 // for |arguments| >= 2**63
3387 // Address to save r and c as double
3392 // r45 sp+32 -> f0 r
3393 // r44 -> sp+16 -> InputX
3394 // sp sp -> scratch provided to callee
3398 .proc __libm_callout_2
3400 L(SINCOS_ARG_TOO_LARGE):
3404 add r45=-32,sp // Parameter: r address
3406 .save ar.pfs,GR_SAVE_PFS
3407 mov GR_SAVE_PFS=ar.pfs // Save ar.pfs
3411 add sp=-64,sp // Create new stack
3413 mov GR_SAVE_GP=gp // Save gp
3416 stfe [r45] = f0,16 // Clear Parameter r on stack
3417 add r44 = 16,sp // Parameter x address
3418 .save b0, GR_SAVE_B0
3419 mov GR_SAVE_B0=b0 // Save b0
3423 stfe [r45] = f0,-16 // Clear Parameter c on stack
3428 stfe [r44] = FR_Input_X // Store Parameter x on stack
3430 br.call.sptk b0=__libm_pi_by_2_reduce# ;;
3435 ldfe FR_Input_X =[r44],16
3437 // Get r and c off stack
3439 adds GR_Table_Base1 = -16, GR_Table_Base1
3441 // Get r and c off stack
3443 add GR_N_Inc = GR_Sin_or_Cos,r8 ;;
3448 // Get X off the stack
3449 // Readjust Table ptr
3451 ldfs FR_Two_to_M3 = [GR_Table_Base1],4
3455 ldfs FR_Neg_Two_to_M3 = [GR_Table_Base1],0
3462 add sp = 64,sp // Restore stack pointer
3463 fcmp.lt.unc.s1 p6, p0 = FR_r, FR_Two_to_M3
3464 mov b0 = GR_SAVE_B0 // Restore return address
3467 mov gp = GR_SAVE_GP // Restore gp
3468 mov ar.pfs = GR_SAVE_PFS // Restore ar.pfs
3475 (p6) fcmp.gt.unc.s1 p6, p0 = FR_r, FR_Neg_Two_to_M3
3482 (p6) br.cond.spnt L(SINCOS_SMALL_R) ;;
3488 br.cond.sptk L(SINCOS_NORMAL_R) ;;
3491 .endp __libm_callout_2
3492 ASM_SIZE_DIRECTIVE(__libm_callout_2)
3494 .type __libm_pi_by_2_reduce#,@function
3495 .global __libm_pi_by_2_reduce#
3498 .type __libm_sin_double_dbx#,@function
3499 .global __libm_sin_double_dbx#
3500 .type __libm_cos_double_dbx#,@function
3501 .global __libm_cos_double_dbx#