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.
11 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
12 // "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
13 // LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
14 // A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL INTEL OR ITS
15 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
16 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
17 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
18 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY
19 // OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY OR TORT (INCLUDING
20 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
21 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
23 // Intel Corporation is the author of this code, and requests that all
24 // problem reports or change requests be submitted to it directly at
25 // http://developer.intel.com/opensource.
27 // History: 02/02/00 Initial Version
29 // *********************************************************************
30 // *********************************************************************
32 // Function: __libm_pi_by_two_reduce(x) return r, c, and N where
33 // x = N * pi/4 + (r+c) , where |r+c| <= pi/4.
34 // This function is not designed to be used by the
37 // *********************************************************************
39 // Accuracy: Returns double-precision values
41 // *********************************************************************
45 // Floating-Point Registers: f32-f70
47 // General Purpose Registers:
48 // r8 = return value N
50 // r33 = Address of where to place r and then c
53 // Predicate Registers: p6-p14
55 // *********************************************************************
57 // IEEE Special Conditions:
59 // No condions should be raised.
61 // *********************************************************************
66 // For the forward trigonometric functions sin, cos, sincos, and
67 // tan, the original algorithms for IA 64 handle arguments up to
68 // 1 ulp less than 2^63 in magnitude. For double-extended arguments x,
69 // |x| >= 2^63, this routine returns CASE, N and r_hi, r_lo where
71 // x is accurately approximated by
72 // 2*K*pi + N * pi/2 + r_hi + r_lo, |r_hi+r_lo| <= pi/4.
74 // CASE is 1 unless |r_hi + r_lo| < 2^(-33).
76 // The exact value of K is not determined, but that information is
77 // not required in trigonometric function computations.
79 // We first assume the argument x in question satisfies x >= 2^(63).
80 // In particular, it is positive. Negative x can be handled by symmetry:
82 // -x is accurately approximated by
83 // -2*K*pi + (-N) * pi/2 - (r_hi + r_lo), |r_hi+r_lo| <= pi/4.
85 // The idea of the reduction is that
87 // x * 2/pi = N_big + N + f, |f| <= 1/2
89 // Moreover, for double extended x, |f| >= 2^(-75). (This is an
90 // non-obvious fact found by enumeration using a special algorithm
91 // involving continued fraction.) The algorithm described below
92 // calculates N and an accurate approximation of f.
94 // Roughly speaking, an appropriate 256-bit (4 X 64) portion of
95 // 2/pi is multiplied with x to give the desired information.
97 // II. Representation of 2/PI
98 // ==========================
100 // The value of 2/pi in binary fixed-point is
102 // .101000101111100110......
104 // We store 2/pi in a table, starting at the position corresponding
105 // to bit position 63
107 // bit position 63 62 ... 0 -1 -2 -3 -4 -5 -6 -7 .... -16576
109 // 0 0 ... 0 . 1 0 1 0 1 0 1 .... X
112 // |__ implied binary pt
117 // This describes the algorithm in the most natural way using
118 // unsigned interger multiplication. The implementation section
119 // describes how the integer arithmetic is simulated.
121 // STEP 0. Initialization
122 // ----------------------
124 // Let the input argument x be
126 // x = 2^m * ( 1. b_1 b_2 b_3 ... b_63 ), 63 <= m <= 16383.
128 // The first crucial step is to fetch four 64-bit portions of 2/pi.
129 // To fulfill this goal, we calculate the bit position L of the
130 // beginning of these 256-bit quantity by
134 // Note that -16321 <= L <= -1 because 63 <= m <= 16383; and that
135 // the storage of 2/pi is adequate.
137 // Fetch P_1, P_2, P_3, P_4 beginning at bit position L thus:
139 // bit position L L-1 L-2 ... L-63
143 // each b can be 0 or 1. Also, let P_0 be the two bits correspoding to
144 // bit positions L+2 and L+1. So, when each of the P_j is interpreted
145 // with appropriate scaling, we have
147 // 2/pi = P_big + P_0 + (P_1 + P_2 + P_3 + P_4) + P_small
149 // Note that P_big and P_small can be ignored. The reasons are as follow.
150 // First, consider P_big. If P_big = 0, we can certainly ignore it.
151 // Otherwise, P_big >= 2^(L+3). Now,
153 // P_big * ulp(x) >= 2^(L+3) * 2^(m-63)
154 // >= 2^(65-m + m-63 )
157 // Thus, P_big * x is an integer of the form 4*K. So
159 // x = 4*K * (pi/2) + x*(P_0 + P_1 + P_2 + P_3 + P_4)*(pi/2)
160 // + x*P_small*(pi/2).
162 // Hence, P_big*x corresponds to information that can be ignored for
163 // trigonometic function evaluation.
165 // Next, we must estimate the effect of ignoring P_small. The absolute
166 // error made by ignoring P_small is bounded by
168 // |P_small * x| <= ulp(P_4) * x
169 // <= 2^(L-255) * 2^(m+1)
170 // <= 2^(62-m-255 + m + 1)
173 // Since for double-extended precision, x * 2/pi = integer + f,
174 // 0.5 >= |f| >= 2^(-75), the relative error introduced by ignoring
175 // P_small is bounded by 2^(-192+75) <= 2^(-117), which is acceptable.
177 // Further note that if x is split into x_hi + x_lo where x_lo is the
178 // two bits corresponding to bit positions 2^(m-62) and 2^(m-63); then
182 // is also an integer of the form 4*K; and thus can also be ignored.
183 // Let M := P_0 * x_lo which is a small integer. The main part of the
184 // calculation is really the multiplication of x with the four pieces
185 // P_1, P_2, P_3, and P_4.
187 // Unless the reduced argument is extremely small in magnitude, it
188 // suffices to carry out the multiplication of x with P_1, P_2, and
189 // P_3. x*P_4 will be carried out and added on as a correction only
190 // when it is found to be needed. Note also that x*P_4 need not be
191 // computed exactly. A straightforward multiplication suffices since
192 // the rounding error thus produced would be bounded by 2^(-3*64),
193 // that is 2^(-192) which is small enough as the reduced argument
194 // is bounded from below by 2^(-75).
196 // Now that we have four 64-bit data representing 2/pi and a
197 // 64-bit x. We first need to calculate a highly accurate product
198 // of x and P_1, P_2, P_3. This is best understood as integer
202 // STEP 1. Multiplication
203 // ----------------------
206 // --------- --------- ---------
207 // | P_1 | | P_2 | | P_3 |
208 // --------- --------- ---------
213 // ----------------------------------------------------
215 // --------- ---------
217 // --------- ---------
220 // --------- ---------
222 // --------- ---------
225 // --------- ---------
227 // --------- ---------
229 // ====================================================
230 // --------- --------- --------- ---------
231 // | S_0 | | S_1 | | S_2 | | S_3 |
232 // --------- --------- --------- ---------
236 // STEP 2. Get N and f
237 // -------------------
239 // Conceptually, after the individual pieces S_0, S_1, ..., are obtained,
240 // we have to sum them and obtain an integer part, N, and a fraction, f.
241 // Here, |f| <= 1/2, and N is an integer. Note also that N need only to
242 // be known to module 2^k, k >= 2. In the case when |f| is small enough,
243 // we would need to add in the value x*P_4.
246 // STEP 3. Get reduced argument
247 // ----------------------------
249 // The value f is not yet the reduced argument that we seek. The
252 // x * 2/pi = 4K + N + f
256 // x = 2*K*pi + N * pi/2 + f * (pi/2).
258 // Thus, the reduced argument is given by
260 // reduced argument = f * pi/2.
262 // This multiplication must be performed to extra precision.
264 // IV. Implementation
265 // ==================
267 // Step 0. Initialization
268 // ----------------------
270 // Set sgn_x := sign(x); x := |x|; x_lo := 2 lsb of x.
272 // In memory, 2/pi is stored contigously as
274 // 0x00000000 0x00000000 0xA2F....
276 // |__ implied binary bit
278 // Given x = 2^m * 1.xxxx...xxx; we calculate L := 62 - m. Thus
279 // -1 <= L <= -16321. We fetch from memory 5 integer pieces of data.
281 // P_0 is the two bits corresponding to bit positions L+2 and L+1
282 // P_1 is the 64-bit starting at bit position L
283 // P_2 is the 64-bit starting at bit position L-64
284 // P_3 is the 64-bit starting at bit position L-128
285 // P_4 is the 64-bit starting at bit position L-192
287 // For example, if m = 63, P_0 would be 0 and P_1 would look like
290 // If m = 65, P_0 would be the two msb of 0xA, thus, P_0 is 10 in binary.
291 // P_1 in binary would be 1 0 0 0 1 0 1 1 1 1 ....
293 // Step 1. Multiplication
294 // ----------------------
296 // At this point, P_1, P_2, P_3, P_4 are integers. They are
297 // supposed to be interpreted as
300 // 2^(L-63-64) * P_2;
301 // 2^(L-63-128) * P_3;
302 // 2^(L-63-192) * P_4;
304 // Since each of them need to be multiplied to x, we would scale
305 // both x and the P_j's by some convenient factors: scale each
306 // of P_j's up by 2^(63-L), and scale x down by 2^(L-63).
308 // p_1 := fcvt.xf ( P_1 )
309 // p_2 := fcvt.xf ( P_2 ) * 2^(-64)
310 // p_3 := fcvt.xf ( P_3 ) * 2^(-128)
311 // p_4 := fcvt.xf ( P_4 ) * 2^(-192)
312 // x := replace exponent of x by -1
313 // because 2^m * 1.xxxx...xxx * 2^(L-63)
314 // is 2^(-1) * 1.xxxx...xxx
316 // We are now faced with the task of computing the following
318 // --------- --------- ---------
319 // | P_1 | | P_2 | | P_3 |
320 // --------- --------- ---------
325 // ----------------------------------------------------
327 // --------- ---------
329 // --------- ---------
331 // --------- ---------
333 // --------- ---------
335 // --------- ---------
337 // --------- ---------
339 // ====================================================
340 // ----------- --------- --------- ---------
341 // | S_0 | | S_1 | | S_2 | | S_3 |
342 // ----------- --------- --------- ---------
344 // | |___ binary point
346 // |___ possibly one more bit
348 // Let FPSR3 be set to round towards zero with widest precision
349 // and exponent range. Unless an explicit FPSR is given,
350 // round-to-nearest with widest precision and exponent range is
353 // Define sigma_C := 2^63; sigma_B := 2^(-1); sigma_C := 2^(-65).
355 // Tmp_C := fmpy.fpsr3( x, p_1 );
356 // If Tmp_C >= sigma_C then
358 // C_lo := x*p_1 - C_hi ...fma, exact
360 // C_hi := fadd.fpsr3(sigma_C, Tmp_C) - sigma_C
361 // ...subtraction is exact, regardless
362 // ...of rounding direction
363 // C_lo := x*p_1 - C_hi ...fma, exact
366 // Tmp_B := fmpy.fpsr3( x, p_2 );
367 // If Tmp_B >= sigma_B then
369 // B_lo := x*p_2 - B_hi ...fma, exact
371 // B_hi := fadd.fpsr3(sigma_B, Tmp_B) - sigma_B
372 // ...subtraction is exact, regardless
373 // ...of rounding direction
374 // B_lo := x*p_2 - B_hi ...fma, exact
377 // Tmp_A := fmpy.fpsr3( x, p_3 );
378 // If Tmp_A >= sigma_A then
380 // A_lo := x*p_3 - A_hi ...fma, exact
382 // A_hi := fadd.fpsr3(sigma_A, Tmp_A) - sigma_A
383 // ...subtraction is exact, regardless
384 // ...of rounding direction
385 // A_lo := x*p_3 - A_hi ...fma, exact
388 // ...Note that C_hi is of integer value. We need only the
389 // ...last few bits. Thus we can ensure C_hi is never a big
390 // ...integer, freeing us from overflow worry.
392 // Tmp_C := fadd.fpsr3( C_hi, 2^(70) ) - 2^(70);
393 // ...Tmp_C is the upper portion of C_hi
394 // C_hi := C_hi - Tmp_C
395 // ...0 <= C_hi < 2^7
397 // Step 2. Get N and f
398 // -------------------
400 // At this point, we have all the components to obtain
401 // S_0, S_1, S_2, S_3 and thus N and f. We start by adding
402 // C_lo and B_hi. This sum together with C_hi gives a good
403 // estimation of N and f.
405 // A := fadd.fpsr3( B_hi, C_lo )
406 // B := max( B_hi, C_lo )
407 // b := min( B_hi, C_lo )
409 // a := (B - A) + b ...exact. Note that a is either 0
412 // N := round_to_nearest_integer_value( A );
413 // f := A - N; ...exact because lsb(A) >= 2^(-64)
414 // ...and |f| <= 1/2.
416 // f := f + a ...exact because a is 0 or 2^(-64);
417 // ...the msb of the sum is <= 1/2
418 // ...lsb >= 2^(-64).
420 // N := convert to integer format( C_hi + N );
424 // If sgn_x == 1 (that is original x was negative)
426 // ...this maintains N to be non-negative, but still
427 // ...equivalent to the (negated N) mod 4.
436 // s_lo := (f - s_hi) + g;
442 // A := fadd.fpsr3( A_hi, B_lo )
443 // B := max( A_hi, B_lo )
444 // b := min( A_hi, B_lo )
446 // a := (B - A) + b ...exact. Note that a is either 0
450 // f_lo := (f - f_hi) + A;
452 // ...f-f_hi is exact because either |f| >= |A|, in which
453 // ...case f-f_hi is clearly exact; or otherwise, 0<|f|<|A|
454 // ...means msb(f) <= msb(A) = 2^(-64) => |f| = 2^(-64).
455 // ...If f = 2^(-64), f-f_hi involves cancellation and is
456 // ...exact. If f = -2^(-64), then A + f is exact. Hence
457 // ...f-f_hi is -A exactly, giving f_lo = 0.
461 // If |f| >= 2^(-50) then
465 // f_lo := (f_lo + A_lo) + x*p_4
466 // s_hi := f_hi + f_lo
467 // s_lo := (f_hi - s_hi) + f_lo
472 // Step 3. Get reduced argument
473 // ----------------------------
475 // If sgn_x == 0 (that is original x is positive)
477 // D_hi := Pi_by_2_hi
478 // D_lo := Pi_by_2_lo
479 // ...load from table
483 // D_hi := neg_Pi_by_2_hi
484 // D_lo := neg_Pi_by_2_lo
485 // ...load from table
489 // r_lo := s_hi*D_hi - r_hi ...fma
490 // r_lo := (s_hi*D_lo + r_lo) + s_lo*D_hi
492 // Return CASE, N, r_hi, r_lo
495 #include "libm_support.h"
538 GR_Address_of_Input = r32
539 GR_Address_of_Outputs = r33
578 Constants_Bits_of_2_by_pi:
579 ASM_TYPE_DIRECTIVE(Constants_Bits_of_2_by_pi,@object)
580 data8 0x0000000000000000,0xA2F9836E4E441529
581 data8 0xFC2757D1F534DDC0,0xDB6295993C439041
582 data8 0xFE5163ABDEBBC561,0xB7246E3A424DD2E0
583 data8 0x06492EEA09D1921C,0xFE1DEB1CB129A73E
584 data8 0xE88235F52EBB4484,0xE99C7026B45F7E41
585 data8 0x3991D639835339F4,0x9C845F8BBDF9283B
586 data8 0x1FF897FFDE05980F,0xEF2F118B5A0A6D1F
587 data8 0x6D367ECF27CB09B7,0x4F463F669E5FEA2D
588 data8 0x7527BAC7EBE5F17B,0x3D0739F78A5292EA
589 data8 0x6BFB5FB11F8D5D08,0x56033046FC7B6BAB
590 data8 0xF0CFBC209AF4361D,0xA9E391615EE61B08
591 data8 0x6599855F14A06840,0x8DFFD8804D732731
592 data8 0x06061556CA73A8C9,0x60E27BC08C6B47C4
593 data8 0x19C367CDDCE8092A,0x8359C4768B961CA6
594 data8 0xDDAF44D15719053E,0xA5FF07053F7E33E8
595 data8 0x32C2DE4F98327DBB,0xC33D26EF6B1E5EF8
596 data8 0x9F3A1F35CAF27F1D,0x87F121907C7C246A
597 data8 0xFA6ED5772D30433B,0x15C614B59D19C3C2
598 data8 0xC4AD414D2C5D000C,0x467D862D71E39AC6
599 data8 0x9B0062337CD2B497,0xA7B4D55537F63ED7
600 data8 0x1810A3FC764D2A9D,0x64ABD770F87C6357
601 data8 0xB07AE715175649C0,0xD9D63B3884A7CB23
602 data8 0x24778AD623545AB9,0x1F001B0AF1DFCE19
603 data8 0xFF319F6A1E666157,0x9947FBACD87F7EB7
604 data8 0x652289E83260BFE6,0xCDC4EF09366CD43F
605 data8 0x5DD7DE16DE3B5892,0x9BDE2822D2E88628
606 data8 0x4D58E232CAC616E3,0x08CB7DE050C017A7
607 data8 0x1DF35BE01834132E,0x6212830148835B8E
608 data8 0xF57FB0ADF2E91E43,0x4A48D36710D8DDAA
609 data8 0x425FAECE616AA428,0x0AB499D3F2A6067F
610 data8 0x775C83C2A3883C61,0x78738A5A8CAFBDD7
611 data8 0x6F63A62DCBBFF4EF,0x818D67C12645CA55
612 data8 0x36D9CAD2A8288D61,0xC277C9121426049B
613 data8 0x4612C459C444C5C8,0x91B24DF31700AD43
614 data8 0xD4E5492910D5FDFC,0xBE00CC941EEECE70
615 data8 0xF53E1380F1ECC3E7,0xB328F8C79405933E
616 data8 0x71C1B3092EF3450B,0x9C12887B20AB9FB5
617 data8 0x2EC292472F327B6D,0x550C90A7721FE76B
618 data8 0x96CB314A1679E279,0x4189DFF49794E884
619 data8 0xE6E29731996BED88,0x365F5F0EFDBBB49A
620 data8 0x486CA46742727132,0x5D8DB8159F09E5BC
621 data8 0x25318D3974F71C05,0x30010C0D68084B58
622 data8 0xEE2C90AA4702E774,0x24D6BDA67DF77248
623 data8 0x6EEF169FA6948EF6,0x91B45153D1F20ACF
624 data8 0x3398207E4BF56863,0xB25F3EDD035D407F
625 data8 0x8985295255C06437,0x10D86D324832754C
626 data8 0x5BD4714E6E5445C1,0x090B69F52AD56614
627 data8 0x9D072750045DDB3B,0xB4C576EA17F9877D
628 data8 0x6B49BA271D296996,0xACCCC65414AD6AE2
629 data8 0x9089D98850722CBE,0xA4049407777030F3
630 data8 0x27FC00A871EA49C2,0x663DE06483DD9797
631 data8 0x3FA3FD94438C860D,0xDE41319D39928C70
632 data8 0xDDE7B7173BDF082B,0x3715A0805C93805A
633 data8 0x921110D8E80FAF80,0x6C4BFFDB0F903876
634 data8 0x185915A562BBCB61,0xB989C7BD401004F2
635 data8 0xD2277549F6B6EBBB,0x22DBAA140A2F2689
636 data8 0x768364333B091A94,0x0EAA3A51C2A31DAE
637 data8 0xEDAF12265C4DC26D,0x9C7A2D9756C0833F
638 data8 0x03F6F0098C402B99,0x316D07B43915200C
639 data8 0x5BC3D8C492F54BAD,0xC6A5CA4ECD37A736
640 data8 0xA9E69492AB6842DD,0xDE6319EF8C76528B
641 data8 0x6837DBFCABA1AE31,0x15DFA1AE00DAFB0C
642 data8 0x664D64B705ED3065,0x29BF56573AFF47B9
643 data8 0xF96AF3BE75DF9328,0x3080ABF68C6615CB
644 data8 0x040622FA1DE4D9A4,0xB33D8F1B5709CD36
645 data8 0xE9424EA4BE13B523,0x331AAAF0A8654FA5
646 data8 0xC1D20F3F0BCD785B,0x76F923048B7B7217
647 data8 0x8953A6C6E26E6F00,0xEBEF584A9BB7DAC4
648 data8 0xBA66AACFCF761D02,0xD12DF1B1C1998C77
649 data8 0xADC3DA4886A05DF7,0xF480C62FF0AC9AEC
650 data8 0xDDBC5C3F6DDED01F,0xC790B6DB2A3A25A3
651 data8 0x9AAF009353AD0457,0xB6B42D297E804BA7
652 data8 0x07DA0EAA76A1597B,0x2A12162DB7DCFDE5
653 data8 0xFAFEDB89FDBE896C,0x76E4FCA90670803E
654 data8 0x156E85FF87FD073E,0x2833676186182AEA
655 data8 0xBD4DAFE7B36E6D8F,0x3967955BBF3148D7
656 data8 0x8416DF30432DC735,0x6125CE70C9B8CB30
657 data8 0xFD6CBFA200A4E46C,0x05A0DD5A476F21D2
658 data8 0x1262845CB9496170,0xE0566B0152993755
659 data8 0x50B7D51EC4F1335F,0x6E13E4305DA92E85
660 data8 0xC3B21D3632A1A4B7,0x08D4B1EA21F716E4
661 data8 0x698F77FF2780030C,0x2D408DA0CD4F99A5
662 data8 0x20D3A2B30A5D2F42,0xF9B4CBDA11D0BE7D
663 data8 0xC1DB9BBD17AB81A2,0xCA5C6A0817552E55
664 data8 0x0027F0147F8607E1,0x640B148D4196DEBE
665 data8 0x872AFDDAB6256B34,0x897BFEF3059EBFB9
666 data8 0x4F6A68A82A4A5AC4,0x4FBCF82D985AD795
667 data8 0xC7F48D4D0DA63A20,0x5F57A4B13F149538
668 data8 0x800120CC86DD71B6,0xDEC9F560BF11654D
669 data8 0x6B0701ACB08CD0C0,0xB24855510EFB1EC3
670 data8 0x72953B06A33540C0,0x7BDC06CC45E0FA29
671 data8 0x4EC8CAD641F3E8DE,0x647CD8649B31BED9
672 data8 0xC397A4D45877C5E3,0x6913DAF03C3ABA46
673 data8 0x18465F7555F5BDD2,0xC6926E5D2EACED44
674 data8 0x0E423E1C87C461E9,0xFD29F3D6E7CA7C22
675 data8 0x35916FC5E0088DD7,0xFFE26A6EC6FDB0C1
676 data8 0x0893745D7CB2AD6B,0x9D6ECD7B723E6A11
677 data8 0xC6A9CFF7DF7329BA,0xC9B55100B70DB2E2
678 data8 0x24BA74607DE58AD8,0x742C150D0C188194
679 data8 0x667E162901767A9F,0xBEFDFDEF4556367E
680 data8 0xD913D9ECB9BA8BFC,0x97C427A831C36EF1
681 data8 0x36C59456A8D8B5A8,0xB40ECCCF2D891234
682 data8 0x576F89562CE3CE99,0xB920D6AA5E6B9C2A
683 data8 0x3ECC5F114A0BFDFB,0xF4E16D3B8E2C86E2
684 data8 0x84D4E9A9B4FCD1EE,0xEFC9352E61392F44
685 data8 0x2138C8D91B0AFC81,0x6A4AFBD81C2F84B4
686 data8 0x538C994ECC2254DC,0x552AD6C6C096190B
687 data8 0xB8701A649569605A,0x26EE523F0F117F11
688 data8 0xB5F4F5CBFC2DBC34,0xEEBC34CC5DE8605E
689 data8 0xDD9B8E67EF3392B8,0x17C99B5861BC57E1
690 data8 0xC68351103ED84871,0xDDDD1C2DA118AF46
691 data8 0x2C21D7F359987AD9,0xC0549EFA864FFC06
692 data8 0x56AE79E536228922,0xAD38DC9367AAE855
693 data8 0x3826829BE7CAA40D,0x51B133990ED7A948
694 data8 0x0569F0B265A7887F,0x974C8836D1F9B392
695 data8 0x214A827B21CF98DC,0x9F405547DC3A74E1
696 data8 0x42EB67DF9DFE5FD4,0x5EA4677B7AACBAA2
697 data8 0xF65523882B55BA41,0x086E59862A218347
698 data8 0x39E6E389D49EE540,0xFB49E956FFCA0F1C
699 data8 0x8A59C52BFA94C5C1,0xD3CFC50FAE5ADB86
700 data8 0xC5476243853B8621,0x94792C8761107B4C
701 data8 0x2A1A2C8012BF4390,0x2688893C78E4C4A8
702 data8 0x7BDBE5C23AC4EAF4,0x268A67F7BF920D2B
703 data8 0xA365B1933D0B7CBD,0xDC51A463DD27DDE1
704 data8 0x6919949A9529A828,0xCE68B4ED09209F44
705 data8 0xCA984E638270237C,0x7E32B90F8EF5A7E7
706 data8 0x561408F1212A9DB5,0x4D7E6F5119A5ABF9
707 data8 0xB5D6DF8261DD9602,0x36169F3AC4A1A283
708 data8 0x6DED727A8D39A9B8,0x825C326B5B2746ED
709 data8 0x34007700D255F4FC,0x4D59018071E0E13F
710 data8 0x89B295F364A8F1AE,0xA74B38FC4CEAB2BB
711 ASM_SIZE_DIRECTIVE(Constants_Bits_of_2_by_pi)
713 Constants_Bits_of_pi_by_2:
714 ASM_TYPE_DIRECTIVE(Constants_Bits_of_pi_by_2,@object)
715 data4 0x2168C234,0xC90FDAA2,0x00003FFF,0x00000000
716 data4 0x80DC1CD1,0xC4C6628B,0x00003FBF,0x00000000
717 ASM_SIZE_DIRECTIVE(Constants_Bits_of_pi_by_2)
720 .proc __libm_pi_by_2_reduce#
721 .global __libm_pi_by_2_reduce#
724 __libm_pi_by_2_reduce:
726 // X is at the address in Address_of_Input
727 // Place the two-piece result at the address in Address_of_Outputs
732 alloc r34 = ar.pfs,2,34,0,0
733 (p0) ldfe FR_X = [GR_Address_of_Input]
734 (p0) fsetc.s3 0x00,0x7F ;;
738 (p0) movl GR_BIASL63 = 0x1003E
744 // 0 0 0 0 0. 1 0 1 0
745 // M 0 1 2 .... 63, 64 65 ... 127, 128
746 // ---------------------------------------------
747 // Segment 0. 1 , 2 , 3
748 // START = M - 63 M = 128 becomes 65
749 // LENGTH1 = START & 0x3F 65 become position 1
750 // SEGMENT = shr(START,6) + 1 0 maps to 1, 64 maps to 2,
751 // LENGTH2 = 64 - LENGTH1
752 // Address_BASE = shladd(SEGMENT,3) + BASE
758 (p0) addl GR_BASE = @ltoff(Constants_Bits_of_2_by_pi#), gp
764 ld8 GR_BASE = [GR_BASE]
773 (p0) movl GR_TEMP5 = 0x000000000000FFFE
777 (p0) setf.exp FR_sigma_B = GR_TEMP5
782 (p0) movl GR_TEMP6 = 0x000000000000FFBE ;;
784 // Define sigma_C := 2^63; sigma_B := 2^(-1); sigma_A := 2^(-65).
786 (p0) setf.exp FR_sigma_A = GR_TEMP6
790 // Special Code for testing DE arguments
791 // (p0) movl GR_BIASL63 = 0x0000000000013FFE
792 // (p0) movl GR_x_lo = 0xFFFFFFFFFFFFFFFF
793 // (p0) setf.exp FR_X = GR_BIASL63
794 // (p0) setf.sig FR_ScaleP3 = GR_x_lo
795 // (p0) fmerge.se FR_X = FR_X,FR_ScaleP3
796 // Set sgn_x := sign(x); x := |x|; x_lo := 2 lsb of x.
797 // 2/pi is stored contigously as
798 // 0x00000000 0x00000000.0xA2F....
799 // M = EXP - BIAS ( M >= 63)
800 // Given x = 2^m * 1.xxxx...xxx; we calculate L := 62 - m.
801 // Thus -1 <= L <= -16321.
803 (p0) getf.exp GR_Exp_x = FR_X
804 (p0) getf.sig GR_x_lo = FR_X
805 (p0) fabs FR_X = FR_X ;;
808 (p0) and GR_x_lo = 0x03,GR_x_lo
809 (p0) extr.u GR_M = GR_Exp_x,0,17 ;;
810 (p0) sub GR_START = GR_M,GR_BIASL63
814 (p0) and GR_LENGTH1 = 0x3F,GR_START
815 (p0) shr.u GR_SEGMENT = GR_START,6
819 (p0) add GR_SEGMENT = 0x1,GR_SEGMENT
820 (p0) sub GR_LENGTH2 = 0x40,GR_LENGTH1
822 // P_0 is the two bits corresponding to bit positions L+2 and L+1
823 // P_1 is the 64-bit starting at bit position L
824 // P_2 is the 64-bit starting at bit position L-64
825 // P_3 is the 64-bit starting at bit position L-128
826 // P_4 is the 64-bit starting at bit position L-192
827 // P_1 is made up of Alo and Bhi
828 // P_1 = deposit Alo, position 0, length2 into P_1,position length1
829 // deposit Bhi, position length2, length1 into P_1, position 0
830 // P_2 is made up of Blo and Chi
831 // P_2 = deposit Blo, position 0, length2 into P_2, position length1
832 // deposit Chi, position length2, length1 into P_2, position 0
833 // P_3 is made up of Clo and Dhi
834 // P_3 = deposit Clo, position 0, length2 into P_3, position length1
835 // deposit Dhi, position length2, length1 into P_3, position 0
836 // P_4 is made up of Clo and Dhi
837 // P_4 = deposit Dlo, position 0, length2 into P_4, position length1
838 // deposit Ehi, position length2, length1 into P_4, position 0
840 (p0) cmp.le.unc p6,p7 = 0x2,GR_LENGTH1 ;;
841 (p0) shladd GR_BASE = GR_SEGMENT,3,GR_BASE
842 (p7) cmp.eq.unc p8,p9 = 0x1,GR_LENGTH1 ;;
846 // ld_64 A at Base and increment Base by 8
847 // ld_64 B at Base and increment Base by 8
848 // ld_64 C at Base and increment Base by 8
849 // ld_64 D at Base and increment Base by 8
850 // ld_64 E at Base and increment Base by 8
852 // ---------------------
853 // A, B, C, D, and E look like | length1 | length2 |
854 // ---------------------
856 (p0) ld8 GR_A = [GR_BASE],8
857 (p0) extr.u GR_sgn_x = GR_Exp_x,17,1 ;;
861 (p0) ld8 GR_B = [GR_BASE],8
862 (p0) fmerge.se FR_X = FR_sigma_B,FR_X ;;
865 (p0) ld8 GR_C = [GR_BASE],8
866 (p8) extr.u GR_Temp = GR_A,63,1 ;;
867 (p0) shl GR_TEMP1 = GR_A,GR_LENGTH1
870 (p0) ld8 GR_D = [GR_BASE],8
872 // P_0 = deposit Ahi, position length2, 2 bit into P_0 at position 0.
873 (p6) shr.u GR_P_0 = GR_A,GR_LENGTH2 ;;
874 (p0) shl GR_TEMP2 = GR_B,GR_LENGTH1
877 (p0) ld8 GR_E = [GR_BASE],-40
878 (p0) shr.u GR_P_1 = GR_B,GR_LENGTH2 ;;
879 (p0) shr.u GR_P_2 = GR_C,GR_LENGTH2
882 // Load 16 bit of ASUB from (Base_Address_of_A - 2)
887 // Deposit element 63 from Ahi and place in element 0 of P_0.
891 (p7) ld2 GR_ASUB = [GR_BASE],8
892 (p0) shl GR_TEMP3 = GR_C,GR_LENGTH1 ;;
893 (p0) shl GR_TEMP4 = GR_D,GR_LENGTH1
897 (p0) shr.u GR_P_3 = GR_D,GR_LENGTH2 ;;
898 (p0) shr.u GR_P_4 = GR_E,GR_LENGTH2
901 (p7) and GR_P_0 = 0x03,GR_ASUB
902 (p6) and GR_P_0 = 0x03,GR_P_0 ;;
903 (p0) or GR_P_1 = GR_P_1,GR_TEMP1
906 (p8) and GR_P_0 = 0x1,GR_P_0 ;;
907 (p0) or GR_P_2 = GR_P_2,GR_TEMP2
908 (p8) shl GR_P_0 = GR_P_0,0x1 ;;
912 (p0) or GR_P_3 = GR_P_3,GR_TEMP3
913 (p8) or GR_P_0 = GR_P_0,GR_Temp
916 (p0) setf.sig FR_p_1 = GR_P_1 ;;
917 (p0) setf.sig FR_p_2 = GR_P_2
918 (p0) or GR_P_4 = GR_P_4,GR_TEMP4 ;;
922 (p0) setf.sig FR_p_3 = GR_P_3
923 (p0) pmpy2.r GR_M = GR_P_0,GR_x_lo
926 (p0) setf.sig FR_p_4 = GR_P_4
927 // P_1, P_2, P_3, P_4 are integers. They should be
929 // 2^(L-63-64) * P_2;
930 // 2^(L-63-128) * P_3;
931 // 2^(L-63-192) * P_4;
932 // Since each of them need to be multiplied to x, we would scale
933 // both x and the P_j's by some convenient factors: scale each
934 // of P_j's up by 2^(63-L), and scale x down by 2^(L-63).
935 // p_1 := fcvt.xf ( P_1 )
936 // p_2 := fcvt.xf ( P_2 ) * 2^(-64)
937 // p_3 := fcvt.xf ( P_3 ) * 2^(-128)
938 // p_4 := fcvt.xf ( P_4 ) * 2^(-192)
939 // x= Set x's exp to -1 because 2^m*1.x...x *2^(L-63)=2^(-1)*1.x...xxx
940 // --------- --------- ---------
941 // | P_1 | | P_2 | | P_3 |
942 // --------- --------- ---------
946 // ----------------------------------------------------
947 // --------- ---------
949 // --------- ---------
950 // --------- ---------
952 // --------- ---------
953 // --------- ---------
955 // --------- ---------
956 // ====================================================
957 // ----------- --------- --------- ---------
958 // | S_0 | | S_1 | | S_2 | | S_3 |
959 // ----------- --------- --------- ---------
960 // | |___ binary point
961 // |___ possibly one more bit
963 // Let FPSR3 be set to round towards zero with widest precision
964 // and exponent range. Unless an explicit FPSR is given,
965 // round-to-nearest with widest precision and exponent range is
967 (p0) movl GR_TEMP1 = 0x000000000000FFBF
971 (p0) setf.exp FR_ScaleP2 = GR_TEMP1
976 (p0) movl GR_TEMP4 = 0x000000000001003E
980 (p0) setf.exp FR_sigma_C = GR_TEMP4
985 (p0) movl GR_TEMP2 = 0x000000000000FF7F ;;
989 (p0) setf.exp FR_ScaleP3 = GR_TEMP2
990 (p0) fcvt.xuf.s1 FR_p_1 = FR_p_1 ;;
994 (p0) fcvt.xuf.s1 FR_p_2 = FR_p_2
999 (p0) movl GR_Temp = 0x000000000000FFDE ;;
1003 (p0) setf.exp FR_TWOM33 = GR_Temp
1004 (p0) fcvt.xuf.s1 FR_p_3 = FR_p_3 ;;
1008 (p0) fcvt.xuf.s1 FR_p_4 = FR_p_4
1013 // Tmp_C := fmpy.fpsr3( x, p_1 );
1014 // Tmp_B := fmpy.fpsr3( x, p_2 );
1015 // Tmp_A := fmpy.fpsr3( x, p_3 );
1016 // If Tmp_C >= sigma_C then
1018 // C_lo := x*p_1 - C_hi ...fma, exact
1020 // C_hi := fadd.fpsr3(sigma_C, Tmp_C) - sigma_C
1021 // C_lo := x*p_1 - C_hi ...fma, exact
1023 // If Tmp_B >= sigma_B then
1025 // B_lo := x*p_2 - B_hi ...fma, exact
1027 // B_hi := fadd.fpsr3(sigma_B, Tmp_B) - sigma_B
1028 // B_lo := x*p_2 - B_hi ...fma, exact
1030 // If Tmp_A >= sigma_A then
1032 // A_lo := x*p_3 - A_hi ...fma, exact
1034 // A_hi := fadd.fpsr3(sigma_A, Tmp_A) - sigma_A
1035 // Exact, regardless ...of rounding direction
1036 // A_lo := x*p_3 - A_hi ...fma, exact
1038 (p0) fmpy.s3 FR_Tmp_C = FR_X,FR_p_1
1043 (p0) fmpy.s1 FR_p_2 = FR_p_2,FR_ScaleP2
1048 (p0) movl GR_Temp = 0x0000000000000400
1052 (p0) movl GR_TEMP3 = 0x000000000000FF3F ;;
1056 (p0) setf.exp FR_ScaleP4 = GR_TEMP3
1057 (p0) fmpy.s1 FR_p_3 = FR_p_3,FR_ScaleP3 ;;
1061 (p0) movl GR_TEMP4 = 0x0000000000010045 ;;
1065 (p0) setf.exp FR_Tmp2_C = GR_TEMP4
1066 (p0) fmpy.s3 FR_Tmp_B = FR_X,FR_p_2 ;;
1070 (p0) fcmp.ge.unc.s1 p12, p9 = FR_Tmp_C,FR_sigma_C
1075 (p0) fmpy.s3 FR_Tmp_A = FR_X,FR_p_3
1080 (p12) mov FR_C_hi = FR_Tmp_C
1084 (p0) addl GR_BASE = @ltoff(Constants_Bits_of_pi_by_2#), gp
1085 (p9) fadd.s3 FR_C_hi = FR_sigma_C,FR_Tmp_C
1093 // Step 3. Get reduced argument
1094 // If sgn_x == 0 (that is original x is positive)
1095 // D_hi := Pi_by_2_hi
1096 // D_lo := Pi_by_2_lo
1099 // D_hi := neg_Pi_by_2_hi
1100 // D_lo := neg_Pi_by_2_lo
1106 ld8 GR_BASE = [GR_BASE]
1114 (p0) ldfe FR_D_hi = [GR_BASE],16
1115 (p0) fmpy.s1 FR_p_4 = FR_p_4,FR_ScaleP4
1119 (p0) ldfe FR_D_lo = [GR_BASE],0
1120 (p0) fcmp.ge.unc.s1 p13, p10 = FR_Tmp_B,FR_sigma_B
1125 (p13) mov FR_B_hi = FR_Tmp_B
1130 (p12) fms.s1 FR_C_lo = FR_X,FR_p_1,FR_C_hi
1135 (p10) fadd.s3 FR_B_hi = FR_sigma_B,FR_Tmp_B
1140 (p9) fsub.s1 FR_C_hi = FR_C_hi,FR_sigma_C
1145 (p0) fcmp.ge.unc.s1 p14, p11 = FR_Tmp_A,FR_sigma_A
1150 (p14) mov FR_A_hi = FR_Tmp_A
1155 (p11) fadd.s3 FR_A_hi = FR_sigma_A,FR_Tmp_A
1160 (p9) fms.s1 FR_C_lo = FR_X,FR_p_1,FR_C_hi
1161 (p0) cmp.eq.unc p12,p9 = 0x1,GR_sgn_x
1165 (p13) fms.s1 FR_B_lo = FR_X,FR_p_2,FR_B_hi
1170 (p10) fsub.s1 FR_B_hi = FR_B_hi,FR_sigma_B
1175 // Note that C_hi is of integer value. We need only the
1176 // last few bits. Thus we can ensure C_hi is never a big
1177 // integer, freeing us from overflow worry.
1178 // Tmp_C := fadd.fpsr3( C_hi, 2^(70) ) - 2^(70);
1179 // Tmp_C is the upper portion of C_hi
1180 (p0) fadd.s3 FR_Tmp_C = FR_C_hi,FR_Tmp2_C
1185 (p14) fms.s1 FR_A_lo = FR_X,FR_p_3,FR_A_hi
1190 (p11) fsub.s1 FR_A_hi = FR_A_hi,FR_sigma_A
1195 // *******************
1196 // Step 2. Get N and f
1197 // *******************
1198 // We have all the components to obtain
1199 // S_0, S_1, S_2, S_3 and thus N and f. We start by adding
1200 // C_lo and B_hi. This sum together with C_hi estimates
1202 // A := fadd.fpsr3( B_hi, C_lo )
1203 // B := max( B_hi, C_lo )
1204 // b := min( B_hi, C_lo )
1205 (p0) fadd.s3 FR_A = FR_B_hi,FR_C_lo
1210 (p10) fms.s1 FR_B_lo = FR_X,FR_p_2,FR_B_hi
1215 (p0) fsub.s1 FR_Tmp_C = FR_Tmp_C,FR_Tmp2_C
1220 (p0) fmax.s1 FR_B = FR_B_hi,FR_C_lo
1225 (p0) fmin.s1 FR_b = FR_B_hi,FR_C_lo
1230 (p11) fms.s1 FR_A_lo = FR_X,FR_p_3,FR_A_hi
1235 // N := round_to_nearest_integer_value( A );
1236 (p0) fcvt.fx.s1 FR_N = FR_A
1241 // C_hi := C_hi - Tmp_C ...0 <= C_hi < 2^7
1242 (p0) fsub.s1 FR_C_hi = FR_C_hi,FR_Tmp_C
1247 // a := (B - A) + b: Exact - note that a is either 0 or 2^(-64).
1248 (p0) fsub.s1 FR_a = FR_B,FR_A
1253 // f := A - N; Exact because lsb(A) >= 2^(-64) and |f| <= 1/2.
1254 (p0) fnorm.s1 FR_N = FR_N
1259 (p0) fadd.s1 FR_a = FR_a,FR_b
1264 (p0) fsub.s1 FR_f = FR_A,FR_N
1269 // N := convert to integer format( C_hi + N );
1272 (p0) fadd.s1 FR_N = FR_N,FR_C_hi
1277 // f = f + a Exact because a is 0 or 2^(-64);
1278 // the msb of the sum is <= 1/2 and lsb >= 2^(-64).
1279 (p0) fadd.s1 FR_f = FR_f,FR_a
1287 (p0) fcvt.fx.s1 FR_N = FR_N
1292 (p0) fabs FR_f_abs = FR_f
1296 (p0) getf.sig GR_N = FR_N
1303 (p0) add GR_N = GR_N,GR_M ;;
1305 // If sgn_x == 1 (that is original x was negative)
1307 // this maintains N to be non-negative, but still
1308 // equivalent to the (negated N) mod 4.
1311 (p12) sub GR_N = GR_Temp,GR_N
1312 (p0) cmp.eq.unc p12,p9 = 0x0,GR_sgn_x ;;
1317 (p0) fcmp.ge.unc.s1 p13, p10 = FR_f_abs,FR_TWOM33
1322 (p9) fsub.s1 FR_D_hi = f0, FR_D_hi
1327 (p10) fadd.s3 FR_A = FR_A_hi,FR_B_lo
1332 (p13) fadd.s1 FR_g = FR_A_hi,FR_B_lo
1337 (p10) fmax.s1 FR_B = FR_A_hi,FR_B_lo
1342 (p9) fsub.s1 FR_D_lo = f0, FR_D_lo
1347 (p10) fmin.s1 FR_b = FR_A_hi,FR_B_lo
1352 (p0) fsetc.s3 0x7F,0x40
1357 (p10) movl GR_Temp = 0x000000000000FFCD ;;
1361 (p10) setf.exp FR_TWOM50 = GR_Temp
1362 (p10) fadd.s1 FR_f_hi = FR_A,FR_f ;;
1366 // a := (B - A) + b Exact.
1367 // Note that a is either 0 or 2^(-128).
1369 // f_lo := (f - f_hi) + A
1370 // f_lo=f-f_hi is exact because either |f| >= |A|, in which
1371 // case f-f_hi is clearly exact; or otherwise, 0<|f|<|A|
1372 // means msb(f) <= msb(A) = 2^(-64) => |f| = 2^(-64).
1373 // If f = 2^(-64), f-f_hi involves cancellation and is
1374 // exact. If f = -2^(-64), then A + f is exact. Hence
1375 // f-f_hi is -A exactly, giving f_lo = 0.
1376 // f_lo := f_lo + a;
1377 (p10) fsub.s1 FR_a = FR_B,FR_A
1382 (p13) fadd.s1 FR_s_hi = FR_f,FR_g
1387 // If |f| >= 2^(-33)
1390 // g := A_hi + B_lo;
1392 // s_lo := (f - s_hi) + g;
1393 (p13) movl GR_CASE = 0x1 ;;
1400 // A := fadd.fpsr3( A_hi, B_lo )
1401 // B := max( A_hi, B_lo )
1402 // b := min( A_hi, B_lo )
1403 (p10) movl GR_CASE = 0x2
1407 (p10) fsub.s1 FR_f_lo = FR_f,FR_f_hi
1412 (p10) fadd.s1 FR_a = FR_a,FR_b
1417 (p13) fsub.s1 FR_s_lo = FR_f,FR_s_hi
1422 (p13) fadd.s1 FR_s_lo = FR_s_lo,FR_g
1427 (p10) fcmp.ge.unc.s1 p14, p11 = FR_f_abs,FR_TWOM50
1434 (p10) fadd.s1 FR_f_lo = FR_f_lo,FR_A
1439 // If |f| >= 2^(-50) then
1443 // f_lo := (f_lo + A_lo) + x*p_4
1444 // s_hi := f_hi + f_lo
1445 // s_lo := (f_hi - s_hi) + f_lo
1447 (p14) mov FR_s_hi = FR_f_hi
1452 (p10) fadd.s1 FR_f_lo = FR_f_lo,FR_a
1457 (p14) mov FR_s_lo = FR_f_lo
1462 (p11) fadd.s1 FR_f_lo = FR_f_lo,FR_A_lo
1467 (p11) fma.s1 FR_f_lo = FR_X,FR_p_4,FR_f_lo
1472 (p11) fadd.s1 FR_s_hi = FR_f_hi,FR_f_lo
1477 // r_hi := s_hi*D_hi
1478 // r_lo := s_hi*D_hi - r_hi with fma
1479 // r_lo := (s_hi*D_lo + r_lo) + s_lo*D_hi
1480 (p0) fmpy.s1 FR_r_hi = FR_s_hi,FR_D_hi
1485 (p11) fsub.s1 FR_s_lo = FR_f_hi,FR_s_hi
1490 (p0) fms.s1 FR_r_lo = FR_s_hi,FR_D_hi,FR_r_hi
1495 (p11) fadd.s1 FR_s_lo = FR_s_lo,FR_f_lo
1500 // Return N, r_hi, r_lo
1501 // We do not return CASE
1502 (p0) stfe [GR_Address_of_Outputs] = FR_r_hi,16
1507 (p0) fma.s1 FR_r_lo = FR_s_hi,FR_D_lo,FR_r_lo
1512 (p0) fma.s1 FR_r_lo = FR_s_lo,FR_D_hi,FR_r_lo
1517 (p0) stfe [GR_Address_of_Outputs] = FR_r_lo,-16
1523 (p0) br.ret.sptk b0 ;;
1526 .endp __libm_pi_by_2_reduce
1527 ASM_SIZE_DIRECTIVE(__libm_pi_by_2_reduce)