1.file "libm_sincosl.s" 2 3 4// Copyright (c) 2000 - 2004, Intel Corporation 5// All rights reserved. 6// 7// 8// Redistribution and use in source and binary forms, with or without 9// modification, are permitted provided that the following conditions are 10// met: 11// 12// * Redistributions of source code must retain the above copyright 13// notice, this list of conditions and the following disclaimer. 14// 15// * Redistributions in binary form must reproduce the above copyright 16// notice, this list of conditions and the following disclaimer in the 17// documentation and/or other materials provided with the distribution. 18// 19// * The name of Intel Corporation may not be used to endorse or promote 20// products derived from this software without specific prior written 21// permission. 22 23// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS 24// "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT 25// LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR 26// A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL INTEL OR ITS 27// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, 28// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, 29// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR 30// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY 31// OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY OR TORT (INCLUDING 32// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS 33// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 34// 35// Intel Corporation is the author of this code, and requests that all 36// problem reports or change requests be submitted to it directly at 37// http://www.intel.com/software/products/opensource/libraries/num.htm. 38// 39//********************************************************************* 40// 41// History: 42// 05/13/02 Initial version of sincosl (based on libm's sinl and cosl) 43// 02/10/03 Reordered header: .section, .global, .proc, .align; 44// used data8 for long double table values 45// 10/13/03 Corrected .file name 46// 02/11/04 cisl is moved to the separate file. 47// 10/26/04 Avoided using r14-31 as scratch so not clobbered by dynamic loader 48// 49//********************************************************************* 50// 51// Function: Combined sincosl routine with 3 different API's 52// 53// API's 54//============================================================== 55// 1) void sincosl(long double, long double*s, long double*c) 56// 2) __libm_sincosl - internal LIBM function, that accepts 57// argument in f8 and returns cosine through f8, sine through f9 58// 59// 60//********************************************************************* 61// 62// Resources Used: 63// 64// Floating-Point Registers: f8 (Input x and cosl return value), 65// f9 (sinl returned) 66// f32-f121 67// 68// General Purpose Registers: 69// r32-r61 70// 71// Predicate Registers: p6-p15 72// 73//********************************************************************* 74// 75// IEEE Special Conditions: 76// 77// Denormal fault raised on denormal inputs 78// Overflow exceptions do not occur 79// Underflow exceptions raised when appropriate for sincosl 80// (No specialized error handling for this routine) 81// Inexact raised when appropriate by algorithm 82// 83// sincosl(SNaN) = QNaN, QNaN 84// sincosl(QNaN) = QNaN, QNaN 85// sincosl(inf) = QNaN, QNaN 86// sincosl(+/-0) = +/-0, 1 87// 88//********************************************************************* 89// 90// Mathematical Description 91// ======================== 92// 93// The computation of FSIN and FCOS performed in parallel. 94// 95// Arg = N pi/2 + alpha, |alpha| <= pi/4. 96// 97// cosl( Arg ) = sinl( (N+1) pi/2 + alpha ), 98// 99// therefore, the code for computing sine will produce cosine as long 100// as 1 is added to N immediately after the argument reduction 101// process. 102// 103// Let M = N if sine 104// N+1 if cosine. 105// 106// Now, given 107// 108// Arg = M pi/2 + alpha, |alpha| <= pi/4, 109// 110// let I = M mod 4, or I be the two lsb of M when M is represented 111// as 2's complement. I = [i_0 i_1]. Then 112// 113// sinl( Arg ) = (-1)^i_0 sinl( alpha ) if i_1 = 0, 114// = (-1)^i_0 cosl( alpha ) if i_1 = 1. 115// 116// For example: 117// if M = -1, I = 11 118// sin ((-pi/2 + alpha) = (-1) cos (alpha) 119// if M = 0, I = 00 120// sin (alpha) = sin (alpha) 121// if M = 1, I = 01 122// sin (pi/2 + alpha) = cos (alpha) 123// if M = 2, I = 10 124// sin (pi + alpha) = (-1) sin (alpha) 125// if M = 3, I = 11 126// sin ((3/2)pi + alpha) = (-1) cos (alpha) 127// 128// The value of alpha is obtained by argument reduction and 129// represented by two working precision numbers r and c where 130// 131// alpha = r + c accurately. 132// 133// The reduction method is described in a previous write up. 134// The argument reduction scheme identifies 4 cases. For Cases 2 135// and 4, because |alpha| is small, sinl(r+c) and cosl(r+c) can be 136// computed very easily by 2 or 3 terms of the Taylor series 137// expansion as follows: 138// 139// Case 2: 140// ------- 141// 142// sinl(r + c) = r + c - r^3/6 accurately 143// cosl(r + c) = 1 - 2^(-67) accurately 144// 145// Case 4: 146// ------- 147// 148// sinl(r + c) = r + c - r^3/6 + r^5/120 accurately 149// cosl(r + c) = 1 - r^2/2 + r^4/24 accurately 150// 151// The only cases left are Cases 1 and 3 of the argument reduction 152// procedure. These two cases will be merged since after the 153// argument is reduced in either cases, we have the reduced argument 154// represented as r + c and that the magnitude |r + c| is not small 155// enough to allow the usage of a very short approximation. 156// 157// The required calculation is either 158// 159// sinl(r + c) = sinl(r) + correction, or 160// cosl(r + c) = cosl(r) + correction. 161// 162// Specifically, 163// 164// sinl(r + c) = sinl(r) + c sin'(r) + O(c^2) 165// = sinl(r) + c cos (r) + O(c^2) 166// = sinl(r) + c(1 - r^2/2) accurately. 167// Similarly, 168// 169// cosl(r + c) = cosl(r) - c sinl(r) + O(c^2) 170// = cosl(r) - c(r - r^3/6) accurately. 171// 172// We therefore concentrate on accurately calculating sinl(r) and 173// cosl(r) for a working-precision number r, |r| <= pi/4 to within 174// 0.1% or so. 175// 176// The greatest challenge of this task is that the second terms of 177// the Taylor series 178// 179// r - r^3/3! + r^r/5! - ... 180// 181// and 182// 183// 1 - r^2/2! + r^4/4! - ... 184// 185// are not very small when |r| is close to pi/4 and the rounding 186// errors will be a concern if simple polynomial accumulation is 187// used. When |r| < 2^-3, however, the second terms will be small 188// enough (6 bits or so of right shift) that a normal Horner 189// recurrence suffices. Hence there are two cases that we consider 190// in the accurate computation of sinl(r) and cosl(r), |r| <= pi/4. 191// 192// Case small_r: |r| < 2^(-3) 193// -------------------------- 194// 195// Since Arg = M pi/4 + r + c accurately, and M mod 4 is [i_0 i_1], 196// we have 197// 198// sinl(Arg) = (-1)^i_0 * sinl(r + c) if i_1 = 0 199// = (-1)^i_0 * cosl(r + c) if i_1 = 1 200// 201// can be accurately approximated by 202// 203// sinl(Arg) = (-1)^i_0 * [sinl(r) + c] if i_1 = 0 204// = (-1)^i_0 * [cosl(r) - c*r] if i_1 = 1 205// 206// because |r| is small and thus the second terms in the correction 207// are unnecessary. 208// 209// Finally, sinl(r) and cosl(r) are approximated by polynomials of 210// moderate lengths. 211// 212// sinl(r) = r + S_1 r^3 + S_2 r^5 + ... + S_5 r^11 213// cosl(r) = 1 + C_1 r^2 + C_2 r^4 + ... + C_5 r^10 214// 215// We can make use of predicates to selectively calculate 216// sinl(r) or cosl(r) based on i_1. 217// 218// Case normal_r: 2^(-3) <= |r| <= pi/4 219// ------------------------------------ 220// 221// This case is more likely than the previous one if one considers 222// r to be uniformly distributed in [-pi/4 pi/4]. Again, 223// 224// sinl(Arg) = (-1)^i_0 * sinl(r + c) if i_1 = 0 225// = (-1)^i_0 * cosl(r + c) if i_1 = 1. 226// 227// Because |r| is now larger, we need one extra term in the 228// correction. sinl(Arg) can be accurately approximated by 229// 230// sinl(Arg) = (-1)^i_0 * [sinl(r) + c(1-r^2/2)] if i_1 = 0 231// = (-1)^i_0 * [cosl(r) - c*r*(1 - r^2/6)] i_1 = 1. 232// 233// Finally, sinl(r) and cosl(r) are approximated by polynomials of 234// moderate lengths. 235// 236// sinl(r) = r + PP_1_hi r^3 + PP_1_lo r^3 + 237// PP_2 r^5 + ... + PP_8 r^17 238// 239// cosl(r) = 1 + QQ_1 r^2 + QQ_2 r^4 + ... + QQ_8 r^16 240// 241// where PP_1_hi is only about 16 bits long and QQ_1 is -1/2. 242// The crux in accurate computation is to calculate 243// 244// r + PP_1_hi r^3 or 1 + QQ_1 r^2 245// 246// accurately as two pieces: U_hi and U_lo. The way to achieve this 247// is to obtain r_hi as a 10 sig. bit number that approximates r to 248// roughly 8 bits or so of accuracy. (One convenient way is 249// 250// r_hi := frcpa( frcpa( r ) ).) 251// 252// This way, 253// 254// r + PP_1_hi r^3 = r + PP_1_hi r_hi^3 + 255// PP_1_hi (r^3 - r_hi^3) 256// = [r + PP_1_hi r_hi^3] + 257// [PP_1_hi (r - r_hi) 258// (r^2 + r_hi r + r_hi^2) ] 259// = U_hi + U_lo 260// 261// Since r_hi is only 10 bit long and PP_1_hi is only 16 bit long, 262// PP_1_hi * r_hi^3 is only at most 46 bit long and thus computed 263// exactly. Furthermore, r and PP_1_hi r_hi^3 are of opposite sign 264// and that there is no more than 8 bit shift off between r and 265// PP_1_hi * r_hi^3. Hence the sum, U_hi, is representable and thus 266// calculated without any error. Finally, the fact that 267// 268// |U_lo| <= 2^(-8) |U_hi| 269// 270// says that U_hi + U_lo is approximating r + PP_1_hi r^3 to roughly 271// 8 extra bits of accuracy. 272// 273// Similarly, 274// 275// 1 + QQ_1 r^2 = [1 + QQ_1 r_hi^2] + 276// [QQ_1 (r - r_hi)(r + r_hi)] 277// = U_hi + U_lo. 278// 279// Summarizing, we calculate r_hi = frcpa( frcpa( r ) ). 280// 281// If i_1 = 0, then 282// 283// U_hi := r + PP_1_hi * r_hi^3 284// U_lo := PP_1_hi * (r - r_hi) * (r^2 + r*r_hi + r_hi^2) 285// poly := PP_1_lo r^3 + PP_2 r^5 + ... + PP_8 r^17 286// correction := c * ( 1 + C_1 r^2 ) 287// 288// Else ...i_1 = 1 289// 290// U_hi := 1 + QQ_1 * r_hi * r_hi 291// U_lo := QQ_1 * (r - r_hi) * (r + r_hi) 292// poly := QQ_2 * r^4 + QQ_3 * r^6 + ... + QQ_8 r^16 293// correction := -c * r * (1 + S_1 * r^2) 294// 295// End 296// 297// Finally, 298// 299// V := poly + ( U_lo + correction ) 300// 301// / U_hi + V if i_0 = 0 302// result := | 303// \ (-U_hi) - V if i_0 = 1 304// 305// It is important that in the last step, negation of U_hi is 306// performed prior to the subtraction which is to be performed in 307// the user-set rounding mode. 308// 309// 310// Algorithmic Description 311// ======================= 312// 313// The argument reduction algorithm shares the same code between FSIN and FCOS. 314// The argument reduction description given 315// previously is repeated below. 316// 317// 318// Step 0. Initialization. 319// 320// Step 1. Check for exceptional and special cases. 321// 322// * If Arg is +-0, +-inf, NaN, NaT, go to Step 10 for special 323// handling. 324// * If |Arg| < 2^24, go to Step 2 for reduction of moderate 325// arguments. This is the most likely case. 326// * If |Arg| < 2^63, go to Step 8 for pre-reduction of large 327// arguments. 328// * If |Arg| >= 2^63, go to Step 10 for special handling. 329// 330// Step 2. Reduction of moderate arguments. 331// 332// If |Arg| < pi/4 ...quick branch 333// N_fix := N_inc (integer) 334// r := Arg 335// c := 0.0 336// Branch to Step 4, Case_1_complete 337// Else ...cf. argument reduction 338// N := Arg * two_by_PI (fp) 339// N_fix := fcvt.fx( N ) (int) 340// N := fcvt.xf( N_fix ) 341// N_fix := N_fix + N_inc 342// s := Arg - N * P_1 (first piece of pi/2) 343// w := -N * P_2 (second piece of pi/2) 344// 345// If |s| >= 2^(-33) 346// go to Step 3, Case_1_reduce 347// Else 348// go to Step 7, Case_2_reduce 349// Endif 350// Endif 351// 352// Step 3. Case_1_reduce. 353// 354// r := s + w 355// c := (s - r) + w ...observe order 356// 357// Step 4. Case_1_complete 358// 359// ...At this point, the reduced argument alpha is 360// ...accurately represented as r + c. 361// If |r| < 2^(-3), go to Step 6, small_r. 362// 363// Step 5. Normal_r. 364// 365// Let [i_0 i_1] by the 2 lsb of N_fix. 366// FR_rsq := r * r 367// r_hi := frcpa( frcpa( r ) ) 368// r_lo := r - r_hi 369// 370// If i_1 = 0, then 371// poly := r*FR_rsq*(PP_1_lo + FR_rsq*(PP_2 + ... FR_rsq*PP_8)) 372// U_hi := r + PP_1_hi*r_hi*r_hi*r_hi ...any order 373// U_lo := PP_1_hi*r_lo*(r*r + r*r_hi + r_hi*r_hi) 374// correction := c + c*C_1*FR_rsq ...any order 375// Else 376// poly := FR_rsq*FR_rsq*(QQ_2 + FR_rsq*(QQ_3 + ... + FR_rsq*QQ_8)) 377// U_hi := 1 + QQ_1 * r_hi * r_hi ...any order 378// U_lo := QQ_1 * r_lo * (r + r_hi) 379// correction := -c*(r + S_1*FR_rsq*r) ...any order 380// Endif 381// 382// V := poly + (U_lo + correction) ...observe order 383// 384// result := (i_0 == 0? 1.0 : -1.0) 385// 386// Last instruction in user-set rounding mode 387// 388// result := (i_0 == 0? result*U_hi + V : 389// result*U_hi - V) 390// 391// Return 392// 393// Step 6. Small_r. 394// 395// ...Use flush to zero mode without causing exception 396// Let [i_0 i_1] be the two lsb of N_fix. 397// 398// FR_rsq := r * r 399// 400// If i_1 = 0 then 401// z := FR_rsq*FR_rsq; z := FR_rsq*z *r 402// poly_lo := S_3 + FR_rsq*(S_4 + FR_rsq*S_5) 403// poly_hi := r*FR_rsq*(S_1 + FR_rsq*S_2) 404// correction := c 405// result := r 406// Else 407// z := FR_rsq*FR_rsq; z := FR_rsq*z 408// poly_lo := C_3 + FR_rsq*(C_4 + FR_rsq*C_5) 409// poly_hi := FR_rsq*(C_1 + FR_rsq*C_2) 410// correction := -c*r 411// result := 1 412// Endif 413// 414// poly := poly_hi + (z * poly_lo + correction) 415// 416// If i_0 = 1, result := -result 417// 418// Last operation. Perform in user-set rounding mode 419// 420// result := (i_0 == 0? result + poly : 421// result - poly ) 422// Return 423// 424// Step 7. Case_2_reduce. 425// 426// ...Refer to the write up for argument reduction for 427// ...rationale. The reduction algorithm below is taken from 428// ...argument reduction description and integrated this. 429// 430// w := N*P_3 431// U_1 := N*P_2 + w ...FMA 432// U_2 := (N*P_2 - U_1) + w ...2 FMA 433// ...U_1 + U_2 is N*(P_2+P_3) accurately 434// 435// r := s - U_1 436// c := ( (s - r) - U_1 ) - U_2 437// 438// ...The mathematical sum r + c approximates the reduced 439// ...argument accurately. Note that although compared to 440// ...Case 1, this case requires much more work to reduce 441// ...the argument, the subsequent calculation needed for 442// ...any of the trigonometric function is very little because 443// ...|alpha| < 1.01*2^(-33) and thus two terms of the 444// ...Taylor series expansion suffices. 445// 446// If i_1 = 0 then 447// poly := c + S_1 * r * r * r ...any order 448// result := r 449// Else 450// poly := -2^(-67) 451// result := 1.0 452// Endif 453// 454// If i_0 = 1, result := -result 455// 456// Last operation. Perform in user-set rounding mode 457// 458// result := (i_0 == 0? result + poly : 459// result - poly ) 460// 461// Return 462// 463// 464// Step 8. Pre-reduction of large arguments. 465// 466// ...Again, the following reduction procedure was described 467// ...in the separate write up for argument reduction, which 468// ...is tightly integrated here. 469 470// N_0 := Arg * Inv_P_0 471// N_0_fix := fcvt.fx( N_0 ) 472// N_0 := fcvt.xf( N_0_fix) 473 474// Arg' := Arg - N_0 * P_0 475// w := N_0 * d_1 476// N := Arg' * two_by_PI 477// N_fix := fcvt.fx( N ) 478// N := fcvt.xf( N_fix ) 479// N_fix := N_fix + N_inc 480// 481// s := Arg' - N * P_1 482// w := w - N * P_2 483// 484// If |s| >= 2^(-14) 485// go to Step 3 486// Else 487// go to Step 9 488// Endif 489// 490// Step 9. Case_4_reduce. 491// 492// ...first obtain N_0*d_1 and -N*P_2 accurately 493// U_hi := N_0 * d_1 V_hi := -N*P_2 494// U_lo := N_0 * d_1 - U_hi V_lo := -N*P_2 - U_hi ...FMAs 495// 496// ...compute the contribution from N_0*d_1 and -N*P_3 497// w := -N*P_3 498// w := w + N_0*d_2 499// t := U_lo + V_lo + w ...any order 500// 501// ...at this point, the mathematical value 502// ...s + U_hi + V_hi + t approximates the true reduced argument 503// ...accurately. Just need to compute this accurately. 504// 505// ...Calculate U_hi + V_hi accurately: 506// A := U_hi + V_hi 507// if |U_hi| >= |V_hi| then 508// a := (U_hi - A) + V_hi 509// else 510// a := (V_hi - A) + U_hi 511// endif 512// ...order in computing "a" must be observed. This branch is 513// ...best implemented by predicates. 514// ...A + a is U_hi + V_hi accurately. Moreover, "a" is 515// ...much smaller than A: |a| <= (1/2)ulp(A). 516// 517// ...Just need to calculate s + A + a + t 518// C_hi := s + A t := t + a 519// C_lo := (s - C_hi) + A 520// C_lo := C_lo + t 521// 522// ...Final steps for reduction 523// r := C_hi + C_lo 524// c := (C_hi - r) + C_lo 525// 526// ...At this point, we have r and c 527// ...And all we need is a couple of terms of the corresponding 528// ...Taylor series. 529// 530// If i_1 = 0 531// poly := c + r*FR_rsq*(S_1 + FR_rsq*S_2) 532// result := r 533// Else 534// poly := FR_rsq*(C_1 + FR_rsq*C_2) 535// result := 1 536// Endif 537// 538// If i_0 = 1, result := -result 539// 540// Last operation. Perform in user-set rounding mode 541// 542// result := (i_0 == 0? result + poly : 543// result - poly ) 544// Return 545// 546// Large Arguments: For arguments above 2**63, a Payne-Hanek 547// style argument reduction is used and pi_by_2 reduce is called. 548// 549 550 551RODATA 552.align 64 553 554LOCAL_OBJECT_START(FSINCOSL_CONSTANTS) 555 556sincosl_table_p: 557//data4 0x4E44152A, 0xA2F9836E, 0x00003FFE,0x00000000 // Inv_pi_by_2 558//data4 0xCE81B9F1, 0xC84D32B0, 0x00004016,0x00000000 // P_0 559//data4 0x2168C235, 0xC90FDAA2, 0x00003FFF,0x00000000 // P_1 560//data4 0xFC8F8CBB, 0xECE675D1, 0x0000BFBD,0x00000000 // P_2 561//data4 0xACC19C60, 0xB7ED8FBB, 0x0000BF7C,0x00000000 // P_3 562//data4 0xDBD171A1, 0x8D848E89, 0x0000BFBF,0x00000000 // d_1 563//data4 0x18A66F8E, 0xD5394C36, 0x0000BF7C,0x00000000 // d_2 564data8 0xA2F9836E4E44152A, 0x00003FFE // Inv_pi_by_2 565data8 0xC84D32B0CE81B9F1, 0x00004016 // P_0 566data8 0xC90FDAA22168C235, 0x00003FFF // P_1 567data8 0xECE675D1FC8F8CBB, 0x0000BFBD // P_2 568data8 0xB7ED8FBBACC19C60, 0x0000BF7C // P_3 569data8 0x8D848E89DBD171A1, 0x0000BFBF // d_1 570data8 0xD5394C3618A66F8E, 0x0000BF7C // d_2 571LOCAL_OBJECT_END(FSINCOSL_CONSTANTS) 572 573LOCAL_OBJECT_START(sincosl_table_d) 574//data4 0x2168C234, 0xC90FDAA2, 0x00003FFE,0x00000000 // pi_by_4 575//data4 0x6EC6B45A, 0xA397E504, 0x00003FE7,0x00000000 // Inv_P_0 576data8 0xC90FDAA22168C234, 0x00003FFE // pi_by_4 577data8 0xA397E5046EC6B45A, 0x00003FE7 // Inv_P_0 578data4 0x3E000000, 0xBE000000 // 2^-3 and -2^-3 579data4 0x2F000000, 0xAF000000 // 2^-33 and -2^-33 580data4 0x9E000000, 0x00000000 // -2^-67 581data4 0x00000000, 0x00000000 // pad 582LOCAL_OBJECT_END(sincosl_table_d) 583 584LOCAL_OBJECT_START(sincosl_table_pp) 585//data4 0xA21C0BC9, 0xCC8ABEBC, 0x00003FCE,0x00000000 // PP_8 586//data4 0x720221DA, 0xD7468A05, 0x0000BFD6,0x00000000 // PP_7 587//data4 0x640AD517, 0xB092382F, 0x00003FDE,0x00000000 // PP_6 588//data4 0xD1EB75A4, 0xD7322B47, 0x0000BFE5,0x00000000 // PP_5 589//data4 0xFFFFFFFE, 0xFFFFFFFF, 0x0000BFFD,0x00000000 // C_1 590//data4 0x00000000, 0xAAAA0000, 0x0000BFFC,0x00000000 // PP_1_hi 591//data4 0xBAF69EEA, 0xB8EF1D2A, 0x00003FEC,0x00000000 // PP_4 592//data4 0x0D03BB69, 0xD00D00D0, 0x0000BFF2,0x00000000 // PP_3 593//data4 0x88888962, 0x88888888, 0x00003FF8,0x00000000 // PP_2 594//data4 0xAAAB0000, 0xAAAAAAAA, 0x0000BFEC,0x00000000 // PP_1_lo 595data8 0xCC8ABEBCA21C0BC9, 0x00003FCE // PP_8 596data8 0xD7468A05720221DA, 0x0000BFD6 // PP_7 597data8 0xB092382F640AD517, 0x00003FDE // PP_6 598data8 0xD7322B47D1EB75A4, 0x0000BFE5 // PP_5 599data8 0xFFFFFFFFFFFFFFFE, 0x0000BFFD // C_1 600data8 0xAAAA000000000000, 0x0000BFFC // PP_1_hi 601data8 0xB8EF1D2ABAF69EEA, 0x00003FEC // PP_4 602data8 0xD00D00D00D03BB69, 0x0000BFF2 // PP_3 603data8 0x8888888888888962, 0x00003FF8 // PP_2 604data8 0xAAAAAAAAAAAB0000, 0x0000BFEC // PP_1_lo 605LOCAL_OBJECT_END(sincosl_table_pp) 606 607LOCAL_OBJECT_START(sincosl_table_qq) 608//data4 0xC2B0FE52, 0xD56232EF, 0x00003FD2 // QQ_8 609//data4 0x2B48DCA6, 0xC9C99ABA, 0x0000BFDA // QQ_7 610//data4 0x9C716658, 0x8F76C650, 0x00003FE2 // QQ_6 611//data4 0xFDA8D0FC, 0x93F27DBA, 0x0000BFE9 // QQ_5 612//data4 0xAAAAAAAA, 0xAAAAAAAA, 0x0000BFFC // S_1 613//data4 0x00000000, 0x80000000, 0x0000BFFE,0x00000000 // QQ_1 614//data4 0x0C6E5041, 0xD00D00D0, 0x00003FEF,0x00000000 // QQ_4 615//data4 0x0B607F60, 0xB60B60B6, 0x0000BFF5,0x00000000 // QQ_3 616//data4 0xAAAAAA9B, 0xAAAAAAAA, 0x00003FFA,0x00000000 // QQ_2 617data8 0xD56232EFC2B0FE52, 0x00003FD2 // QQ_8 618data8 0xC9C99ABA2B48DCA6, 0x0000BFDA // QQ_7 619data8 0x8F76C6509C716658, 0x00003FE2 // QQ_6 620data8 0x93F27DBAFDA8D0FC, 0x0000BFE9 // QQ_5 621data8 0xAAAAAAAAAAAAAAAA, 0x0000BFFC // S_1 622data8 0x8000000000000000, 0x0000BFFE // QQ_1 623data8 0xD00D00D00C6E5041, 0x00003FEF // QQ_4 624data8 0xB60B60B60B607F60, 0x0000BFF5 // QQ_3 625data8 0xAAAAAAAAAAAAAA9B, 0x00003FFA // QQ_2 626LOCAL_OBJECT_END(sincosl_table_qq) 627 628LOCAL_OBJECT_START(sincosl_table_c) 629//data4 0xFFFFFFFE, 0xFFFFFFFF, 0x0000BFFD,0x00000000 // C_1 630//data4 0xAAAA719F, 0xAAAAAAAA, 0x00003FFA,0x00000000 // C_2 631//data4 0x0356F994, 0xB60B60B6, 0x0000BFF5,0x00000000 // C_3 632//data4 0xB2385EA9, 0xD00CFFD5, 0x00003FEF,0x00000000 // C_4 633//data4 0x292A14CD, 0x93E4BD18, 0x0000BFE9,0x00000000 // C_5 634data8 0xFFFFFFFFFFFFFFFE, 0x0000BFFD // C_1 635data8 0xAAAAAAAAAAAA719F, 0x00003FFA // C_2 636data8 0xB60B60B60356F994, 0x0000BFF5 // C_3 637data8 0xD00CFFD5B2385EA9, 0x00003FEF // C_4 638data8 0x93E4BD18292A14CD, 0x0000BFE9 // C_5 639LOCAL_OBJECT_END(sincosl_table_c) 640 641LOCAL_OBJECT_START(sincosl_table_s) 642//data4 0xAAAAAAAA, 0xAAAAAAAA, 0x0000BFFC,0x00000000 // S_1 643//data4 0x888868DB, 0x88888888, 0x00003FF8,0x00000000 // S_2 644//data4 0x055EFD4B, 0xD00D00D0, 0x0000BFF2,0x00000000 // S_3 645//data4 0x839730B9, 0xB8EF1C5D, 0x00003FEC,0x00000000 // S_4 646//data4 0xE5B3F492, 0xD71EA3A4, 0x0000BFE5,0x00000000 // S_5 647data8 0xAAAAAAAAAAAAAAAA, 0x0000BFFC // S_1 648data8 0x88888888888868DB, 0x00003FF8 // S_2 649data8 0xD00D00D0055EFD4B, 0x0000BFF2 // S_3 650data8 0xB8EF1C5D839730B9, 0x00003FEC // S_4 651data8 0xD71EA3A4E5B3F492, 0x0000BFE5 // S_5 652data4 0x38800000, 0xB8800000 // two**-14 and -two**-14 653LOCAL_OBJECT_END(sincosl_table_s) 654 655FR_Input_X = f8 656FR_Result = f8 657FR_ResultS = f9 658FR_ResultC = f8 659FR_r = f8 660FR_c = f9 661 662FR_norm_x = f9 663FR_inv_pi_2to63 = f10 664FR_rshf_2to64 = f11 665FR_2tom64 = f12 666FR_rshf = f13 667FR_N_float_signif = f14 668FR_abs_x = f15 669 670FR_r6 = f32 671FR_r7 = f33 672FR_Pi_by_4 = f34 673FR_Two_to_M14 = f35 674FR_Neg_Two_to_M14 = f36 675FR_Two_to_M33 = f37 676FR_Neg_Two_to_M33 = f38 677FR_Neg_Two_to_M67 = f39 678FR_Inv_pi_by_2 = f40 679FR_N_float = f41 680FR_N_fix = f42 681FR_P_1 = f43 682FR_P_2 = f44 683FR_P_3 = f45 684FR_s = f46 685FR_w = f47 686FR_Z = f50 687FR_A = f51 688FR_a = f52 689FR_t = f53 690FR_U_1 = f54 691FR_U_2 = f55 692FR_C_1 = f56 693FR_C_2 = f57 694FR_C_3 = f58 695FR_C_4 = f59 696FR_C_5 = f60 697FR_S_1 = f61 698FR_S_2 = f62 699FR_S_3 = f63 700FR_S_4 = f64 701FR_S_5 = f65 702FR_r_hi = f68 703FR_r_lo = f69 704FR_rsq = f70 705FR_r_cubed = f71 706FR_C_hi = f72 707FR_N_0 = f73 708FR_d_1 = f74 709FR_V_hi = f75 710FR_V_lo = f76 711FR_U_hi = f77 712FR_U_lo = f78 713FR_U_hiabs = f79 714FR_V_hiabs = f80 715FR_PP_8 = f81 716FR_QQ_8 = f101 717FR_PP_7 = f82 718FR_QQ_7 = f102 719FR_PP_6 = f83 720FR_QQ_6 = f103 721FR_PP_5 = f84 722FR_QQ_5 = f104 723FR_PP_4 = f85 724FR_QQ_4 = f105 725FR_PP_3 = f86 726FR_QQ_3 = f106 727FR_PP_2 = f87 728FR_QQ_2 = f107 729FR_QQ_1 = f108 730FR_r_hi_sq = f88 731FR_N_0_fix = f89 732FR_Inv_P_0 = f90 733FR_d_2 = f93 734FR_P_0 = f95 735FR_C_lo = f96 736FR_PP_1 = f97 737FR_PP_1_lo = f98 738FR_ArgPrime = f99 739FR_inexact = f100 740 741FR_Neg_Two_to_M3 = f109 742FR_Two_to_M3 = f110 743 744FR_poly_hiS = f66 745FR_poly_hiC = f112 746 747FR_poly_loS = f67 748FR_poly_loC = f113 749 750FR_polyS = f92 751FR_polyC = f114 752 753FR_cS = FR_c 754FR_cC = f115 755 756FR_corrS = f91 757FR_corrC = f116 758 759FR_U_hiC = f117 760FR_U_loC = f118 761 762FR_VS = f75 763FR_VC = f119 764 765FR_FirstS = f120 766FR_FirstC = f121 767 768FR_U_hiS = FR_U_hi 769FR_U_loS = FR_U_lo 770 771FR_Tmp = f94 772 773 774 775 776sincos_pResSin = r34 777sincos_pResCos = r35 778 779GR_exp_m2_to_m3= r36 780GR_N_Inc = r37 781GR_Cis = r38 782GR_signexp_x = r40 783GR_exp_x = r40 784GR_exp_mask = r41 785GR_exp_2_to_63 = r42 786GR_exp_2_to_m3 = r43 787GR_exp_2_to_24 = r44 788 789GR_N_SignS = r45 790GR_N_SignC = r46 791GR_N_SinCos = r47 792 793GR_sig_inv_pi = r48 794GR_rshf_2to64 = r49 795GR_exp_2tom64 = r50 796GR_rshf = r51 797GR_ad_p = r52 798GR_ad_d = r53 799GR_ad_pp = r54 800GR_ad_qq = r55 801GR_ad_c = r56 802GR_ad_s = r57 803GR_ad_ce = r58 804GR_ad_se = r59 805GR_ad_m14 = r60 806GR_ad_s1 = r61 807 808// For unwind support 809GR_SAVE_B0 = r39 810GR_SAVE_GP = r40 811GR_SAVE_PFS = r41 812 813 814.section .text 815 816GLOBAL_IEEE754_ENTRY(sincosl) 817{ .mlx ///////////////////////////// 1 ///////////////// 818 alloc r32 = ar.pfs,3,27,2,0 819 movl GR_sig_inv_pi = 0xa2f9836e4e44152a // significand of 1/pi 820} 821{ .mlx 822 mov GR_N_Inc = 0x0 823 movl GR_rshf_2to64 = 0x47e8000000000000 // 1.1000 2^(63+64) 824};; 825 826{ .mfi ///////////////////////////// 2 ///////////////// 827 addl GR_ad_p = @ltoff(FSINCOSL_CONSTANTS#), gp 828 fclass.m p6, p0 = FR_Input_X, 0x1E3 // Test x natval, nan, inf 829 mov GR_exp_2_to_m3 = 0xffff - 3 // Exponent of 2^-3 830} 831{ .mfb 832 mov GR_Cis = 0x0 833 fnorm.s1 FR_norm_x = FR_Input_X // Normalize x 834 br.cond.sptk _COMMON_SINCOSL 835};; 836GLOBAL_IEEE754_END(sincosl) 837libm_alias_ldouble_other (__sincos, sincos) 838 839GLOBAL_LIBM_ENTRY(__libm_sincosl) 840{ .mlx ///////////////////////////// 1 ///////////////// 841 alloc r32 = ar.pfs,3,27,2,0 842 movl GR_sig_inv_pi = 0xa2f9836e4e44152a // significand of 1/pi 843} 844{ .mlx 845 mov GR_N_Inc = 0x0 846 movl GR_rshf_2to64 = 0x47e8000000000000 // 1.1000 2^(63+64) 847};; 848 849{ .mfi ///////////////////////////// 2 ///////////////// 850 addl GR_ad_p = @ltoff(FSINCOSL_CONSTANTS#), gp 851 fclass.m p6, p0 = FR_Input_X, 0x1E3 // Test x natval, nan, inf 852 mov GR_exp_2_to_m3 = 0xffff - 3 // Exponent of 2^-3 853} 854{ .mfb 855 mov GR_Cis = 0x1 856 fnorm.s1 FR_norm_x = FR_Input_X // Normalize x 857 nop.b 0 858};; 859 860_COMMON_SINCOSL: 861{ .mfi ///////////////////////////// 3 ///////////////// 862 setf.sig FR_inv_pi_2to63 = GR_sig_inv_pi // Form 1/pi * 2^63 863 nop.f 0 864 mov GR_exp_2tom64 = 0xffff - 64 // Scaling constant to compute N 865} 866{ .mlx 867 setf.d FR_rshf_2to64 = GR_rshf_2to64 // Form const 1.1000 * 2^(63+64) 868 movl GR_rshf = 0x43e8000000000000 // Form const 1.1000 * 2^63 869};; 870 871{ .mfi ///////////////////////////// 4 ///////////////// 872 ld8 GR_ad_p = [GR_ad_p] // Point to Inv_pi_by_2 873 fclass.m p7, p0 = FR_Input_X, 0x0b // Test x denormal 874 nop.i 0 875};; 876 877{ .mfi ///////////////////////////// 5 ///////////////// 878 getf.exp GR_signexp_x = FR_Input_X // Get sign and exponent of x 879 fclass.m p10, p0 = FR_Input_X, 0x007 // Test x zero 880 nop.i 0 881} 882{ .mib 883 mov GR_exp_mask = 0x1ffff // Exponent mask 884 nop.i 0 885(p6) br.cond.spnt SINCOSL_SPECIAL // Branch if x natval, nan, inf 886};; 887 888{ .mfi ///////////////////////////// 6 ///////////////// 889 setf.exp FR_2tom64 = GR_exp_2tom64 // Form 2^-64 for scaling N_float 890 nop.f 0 891 add GR_ad_d = 0x70, GR_ad_p // Point to constant table d 892} 893{ .mib 894 setf.d FR_rshf = GR_rshf // Form right shift const 1.1000 * 2^63 895 mov GR_exp_m2_to_m3 = 0x2fffc // Form -(2^-3) 896(p7) br.cond.spnt SINCOSL_DENORMAL // Branch if x denormal 897};; 898 899SINCOSL_COMMON2: 900{ .mfi ///////////////////////////// 7 ///////////////// 901 and GR_exp_x = GR_exp_mask, GR_signexp_x // Get exponent of x 902 fclass.nm p8, p0 = FR_Input_X, 0x1FF // Test x unsupported type 903 mov GR_exp_2_to_63 = 0xffff + 63 // Exponent of 2^63 904} 905{ .mib 906 add GR_ad_pp = 0x40, GR_ad_d // Point to constant table pp 907 mov GR_exp_2_to_24 = 0xffff + 24 // Exponent of 2^24 908(p10) br.cond.spnt SINCOSL_ZERO // Branch if x zero 909};; 910 911{ .mfi ///////////////////////////// 8 ///////////////// 912 ldfe FR_Inv_pi_by_2 = [GR_ad_p], 16 // Load 2/pi 913 fcmp.eq.s0 p15, p0 = FR_Input_X, f0 // Dummy to set denormal 914 add GR_ad_qq = 0xa0, GR_ad_pp // Point to constant table qq 915} 916{ .mfi 917 ldfe FR_Pi_by_4 = [GR_ad_d], 16 // Load pi/4 for range test 918 nop.f 0 919 cmp.ge p10,p0 = GR_exp_x, GR_exp_2_to_63 // Is |x| >= 2^63 920};; 921 922{ .mfi ///////////////////////////// 9 ///////////////// 923 ldfe FR_P_0 = [GR_ad_p], 16 // Load P_0 for pi/4 <= |x| < 2^63 924 fmerge.s FR_abs_x = f1, FR_norm_x // |x| 925 add GR_ad_c = 0x90, GR_ad_qq // Point to constant table c 926} 927{ .mfi 928 ldfe FR_Inv_P_0 = [GR_ad_d], 16 // Load 1/P_0 for pi/4 <= |x| < 2^63 929 nop.f 0 930 cmp.ge p7,p0 = GR_exp_x, GR_exp_2_to_24 // Is |x| >= 2^24 931};; 932 933{ .mfi ///////////////////////////// 10 ///////////////// 934 ldfe FR_P_1 = [GR_ad_p], 16 // Load P_1 for pi/4 <= |x| < 2^63 935 nop.f 0 936 add GR_ad_s = 0x50, GR_ad_c // Point to constant table s 937} 938{ .mfi 939 ldfe FR_PP_8 = [GR_ad_pp], 16 // Load PP_8 for 2^-3 < |r| < pi/4 940 nop.f 0 941 nop.i 0 942};; 943 944{ .mfi ///////////////////////////// 11 ///////////////// 945 ldfe FR_P_2 = [GR_ad_p], 16 // Load P_2 for pi/4 <= |x| < 2^63 946 nop.f 0 947 add GR_ad_ce = 0x40, GR_ad_c // Point to end of constant table c 948} 949{ .mfi 950 ldfe FR_QQ_8 = [GR_ad_qq], 16 // Load QQ_8 for 2^-3 < |r| < pi/4 951 nop.f 0 952 nop.i 0 953};; 954 955{ .mfi ///////////////////////////// 12 ///////////////// 956 ldfe FR_QQ_7 = [GR_ad_qq], 16 // Load QQ_7 for 2^-3 < |r| < pi/4 957 fma.s1 FR_N_float_signif = FR_Input_X, FR_inv_pi_2to63, FR_rshf_2to64 958 add GR_ad_se = 0x40, GR_ad_s // Point to end of constant table s 959} 960{ .mib 961 ldfe FR_PP_7 = [GR_ad_pp], 16 // Load PP_7 for 2^-3 < |r| < pi/4 962 mov GR_ad_s1 = GR_ad_s // Save pointer to S_1 963(p10) br.cond.spnt SINCOSL_ARG_TOO_LARGE // Branch if |x| >= 2^63 964 // Use Payne-Hanek Reduction 965};; 966 967{ .mfi ///////////////////////////// 13 ///////////////// 968 ldfe FR_P_3 = [GR_ad_p], 16 // Load P_3 for pi/4 <= |x| < 2^63 969 fmerge.se FR_r = FR_norm_x, FR_norm_x // r = x, in case |x| < pi/4 970 add GR_ad_m14 = 0x50, GR_ad_s // Point to constant table m14 971} 972{ .mfb 973 ldfps FR_Two_to_M3, FR_Neg_Two_to_M3 = [GR_ad_d], 8 974 fma.s1 FR_rsq = FR_norm_x, FR_norm_x, f0 // rsq = x*x, in case |x| < pi/4 975(p7) br.cond.spnt SINCOSL_LARGER_ARG // Branch if 2^24 <= |x| < 2^63 976 // Use pre-reduction 977};; 978 979{ .mmf ///////////////////////////// 14 ///////////////// 980 ldfe FR_PP_6 = [GR_ad_pp], 16 // Load PP_6 for normal path 981 ldfe FR_QQ_6 = [GR_ad_qq], 16 // Load QQ_6 for normal path 982 fmerge.se FR_c = f0, f0 // c = 0 in case |x| < pi/4 983};; 984 985{ .mmf ///////////////////////////// 15 ///////////////// 986 ldfe FR_PP_5 = [GR_ad_pp], 16 // Load PP_5 for normal path 987 ldfe FR_QQ_5 = [GR_ad_qq], 16 // Load QQ_5 for normal path 988 nop.f 0 989};; 990 991// Here if 0 < |x| < 2^24 992{ .mfi ///////////////////////////// 17 ///////////////// 993 ldfe FR_S_5 = [GR_ad_se], -16 // Load S_5 if i_1=0 994 fcmp.lt.s1 p6, p7 = FR_abs_x, FR_Pi_by_4 // Test |x| < pi/4 995 nop.i 0 996} 997{ .mfi 998 ldfe FR_C_5 = [GR_ad_ce], -16 // Load C_5 if i_1=1 999 fms.s1 FR_N_float = FR_N_float_signif, FR_2tom64, FR_rshf 1000 nop.i 0 1001};; 1002 1003{ .mmi ///////////////////////////// 18 ///////////////// 1004 ldfe FR_S_4 = [GR_ad_se], -16 // Load S_4 if i_1=0 1005 ldfe FR_C_4 = [GR_ad_ce], -16 // Load C_4 if i_1=1 1006 nop.i 0 1007};; 1008 1009// 1010// N = Arg * 2/pi 1011// Check if Arg < pi/4 1012// 1013// 1014// Case 2: Convert integer N_fix back to normalized floating-point value. 1015// Case 1: p8 is only affected when p6 is set 1016// 1017// 1018// Grab the integer part of N and call it N_fix 1019// 1020{ .mfi ///////////////////////////// 19 ///////////////// 1021(p7) ldfps FR_Two_to_M33, FR_Neg_Two_to_M33 = [GR_ad_d], 8 1022(p6) fma.s1 FR_r_cubed = FR_r, FR_rsq, f0 // r^3 if |x| < pi/4 1023(p6) mov GR_N_Inc = 0x0 // N_IncS if |x| < pi/4 1024};; 1025 1026// If |x| < pi/4, r = x and c = 0 1027// lf |x| < pi/4, is x < 2**(-3). 1028// r = Arg 1029// c = 0 1030{ .mmi ///////////////////////////// 20 ///////////////// 1031(p7) getf.sig GR_N_Inc = FR_N_float_signif 1032 nop.m 0 1033(p6) cmp.lt.unc p8,p0 = GR_exp_x, GR_exp_2_to_m3 // Is |x| < 2^-3 1034};; 1035 1036// 1037// lf |x| < pi/4, is -2**(-3)< x < 2**(-3) - set p8. 1038// If |x| >= pi/4, 1039// Create the right N for |x| < pi/4 and otherwise 1040// Case 2: Place integer part of N in GP register 1041// 1042 1043{ .mbb ///////////////////////////// 21 ///////////////// 1044 nop.m 0 1045(p8) br.cond.spnt SINCOSL_SMALL_R_0 // Branch if 0 < |x| < 2^-3 1046(p6) br.cond.spnt SINCOSL_NORMAL_R_0 // Branch if 2^-3 <= |x| < pi/4 1047};; 1048 1049// Here if pi/4 <= |x| < 2^24 1050{ .mfi 1051 ldfs FR_Neg_Two_to_M67 = [GR_ad_d], 8 // Load -2^-67 1052 fnma.s1 FR_s = FR_N_float, FR_P_1, FR_Input_X // s = -N * P_1 + Arg 1053 nop.i 0 1054} 1055{ .mfi 1056 nop.m 0 1057 fma.s1 FR_w = FR_N_float, FR_P_2, f0 // w = N * P_2 1058 nop.i 0 1059};; 1060 1061{ .mfi 1062 nop.m 0 1063 fms.s1 FR_r = FR_s, f1, FR_w // r = s - w, assume |s| >= 2^-33 1064 nop.i 0 1065};; 1066 1067{ .mfi 1068 nop.m 0 1069 fcmp.lt.s1 p7, p6 = FR_s, FR_Two_to_M33 1070 nop.i 0 1071};; 1072 1073{ .mfi 1074 nop.m 0 1075(p7) fcmp.gt.s1 p7, p6 = FR_s, FR_Neg_Two_to_M33 // p6 if |s| >= 2^-33, else p7 1076 nop.i 0 1077};; 1078 1079{ .mfi 1080 nop.m 0 1081 fms.s1 FR_c = FR_s, f1, FR_r // c = s - r, for |s| >= 2^-33 1082 nop.i 0 1083} 1084{ .mfi 1085 nop.m 0 1086 fma.s1 FR_rsq = FR_r, FR_r, f0 // rsq = r * r, for |s| >= 2^-33 1087 nop.i 0 1088};; 1089 1090{ .mfi 1091 nop.m 0 1092(p7) fma.s1 FR_w = FR_N_float, FR_P_3, f0 1093 nop.i 0 1094};; 1095 1096{ .mmf 1097 ldfe FR_C_1 = [GR_ad_pp], 16 // Load C_1 if i_1=0 1098 ldfe FR_S_1 = [GR_ad_qq], 16 // Load S_1 if i_1=1 1099 frcpa.s1 FR_r_hi, p15 = f1, FR_r // r_hi = frcpa(r) 1100};; 1101 1102{ .mfi 1103 nop.m 0 1104(p6) fcmp.lt.unc.s1 p8, p13 = FR_r, FR_Two_to_M3 // If big s, test r with 2^-3 1105 nop.i 0 1106};; 1107 1108{ .mfi 1109 nop.m 0 1110(p7) fma.s1 FR_U_1 = FR_N_float, FR_P_2, FR_w 1111 nop.i 0 1112};; 1113 1114// 1115// For big s: r = s - w: No futher reduction is necessary 1116// For small s: w = N * P_3 (change sign) More reduction 1117// 1118{ .mfi 1119 nop.m 0 1120(p8) fcmp.gt.s1 p8, p13 = FR_r, FR_Neg_Two_to_M3 // If big s, p8 if |r| < 2^-3 1121 nop.i 0 1122};; 1123 1124{ .mfi 1125 nop.m 0 1126 fma.s1 FR_polyS = FR_rsq, FR_PP_8, FR_PP_7 // poly = rsq*PP_8+PP_7 1127 nop.i 0 1128} 1129{ .mfi 1130 nop.m 0 1131 fma.s1 FR_polyC = FR_rsq, FR_QQ_8, FR_QQ_7 // poly = rsq*QQ_8+QQ_7 1132 nop.i 0 1133};; 1134 1135{ .mfi 1136 nop.m 0 1137(p7) fms.s1 FR_r = FR_s, f1, FR_U_1 1138 nop.i 0 1139};; 1140 1141{ .mfi 1142 nop.m 0 1143(p6) fma.s1 FR_r_cubed = FR_r, FR_rsq, f0 // rcubed = r * rsq 1144 nop.i 0 1145};; 1146 1147{ .mfi 1148// 1149// For big s: Is |r| < 2**(-3)? 1150// For big s: c = S - r 1151// For small s: U_1 = N * P_2 + w 1152// 1153// If p8 is set, prepare to branch to Small_R. 1154// If p9 is set, prepare to branch to Normal_R. 1155// For big s, r is complete here. 1156// 1157// 1158// For big s: c = c + w (w has not been negated.) 1159// For small s: r = S - U_1 1160// 1161 nop.m 0 1162(p6) fms.s1 FR_c = FR_c, f1, FR_w 1163 nop.i 0 1164} 1165{ .mbb 1166 nop.m 0 1167(p8) br.cond.spnt SINCOSL_SMALL_R_1 // Branch if |s|>=2^-33, |r| < 2^-3, 1168 // and pi/4 <= |x| < 2^24 1169(p13) br.cond.sptk SINCOSL_NORMAL_R_1 // Branch if |s|>=2^-33, |r| >= 2^-3, 1170 // and pi/4 <= |x| < 2^24 1171};; 1172 1173SINCOSL_S_TINY: 1174// 1175// Here if |s| < 2^-33, and pi/4 <= |x| < 2^24 1176// 1177{ .mfi 1178 and GR_N_SinCos = 0x1, GR_N_Inc 1179 fms.s1 FR_U_2 = FR_N_float, FR_P_2, FR_U_1 1180 tbit.z p8,p12 = GR_N_Inc, 0 1181};; 1182 1183 1184// 1185// For small s: U_2 = N * P_2 - U_1 1186// S_1 stored constant - grab the one stored with the 1187// coefficients. 1188// 1189{ .mfi 1190 ldfe FR_S_1 = [GR_ad_s1], 16 1191 fma.s1 FR_polyC = f0, f1, FR_Neg_Two_to_M67 1192 sub GR_N_SignS = GR_N_Inc, GR_N_SinCos 1193} 1194{ .mfi 1195 add GR_N_SignC = GR_N_Inc, GR_N_SinCos 1196 nop.f 0 1197 nop.i 0 1198};; 1199 1200{ .mfi 1201 nop.m 0 1202 fms.s1 FR_s = FR_s, f1, FR_r 1203(p8) tbit.z.unc p10,p11 = GR_N_SignC, 1 1204} 1205{ .mfi 1206 nop.m 0 1207 fma.s1 FR_rsq = FR_r, FR_r, f0 1208 nop.i 0 1209};; 1210 1211{ .mfi 1212 nop.m 0 1213 fma.s1 FR_U_2 = FR_U_2, f1, FR_w 1214(p8) tbit.z.unc p8,p9 = GR_N_SignS, 1 1215};; 1216 1217{ .mfi 1218 nop.m 0 1219 fmerge.se FR_FirstS = FR_r, FR_r 1220(p12) tbit.z.unc p14,p15 = GR_N_SignC, 1 1221} 1222{ .mfi 1223 nop.m 0 1224 fma.s1 FR_FirstC = f0, f1, f1 1225 nop.i 0 1226};; 1227 1228{ .mfi 1229 nop.m 0 1230 fms.s1 FR_c = FR_s, f1, FR_U_1 1231(p12) tbit.z.unc p12,p13 = GR_N_SignS, 1 1232};; 1233 1234{ .mfi 1235 nop.m 0 1236 fma.s1 FR_r = FR_S_1, FR_r, f0 1237 nop.i 0 1238};; 1239 1240{ .mfi 1241 nop.m 0 1242 fma.s0 FR_S_1 = FR_S_1, FR_S_1, f0 1243 nop.i 0 1244};; 1245 1246{ .mfi 1247 nop.m 0 1248 fms.s1 FR_c = FR_c, f1, FR_U_2 1249 nop.i 0 1250};; 1251 1252.pred.rel "mutex",p9,p15 1253{ .mfi 1254 nop.m 0 1255(p9) fms.s0 FR_FirstS = f1, f0, FR_FirstS 1256 nop.i 0 1257} 1258{ .mfi 1259 nop.m 0 1260(p15) fms.s0 FR_FirstS = f1, f0, FR_FirstS 1261 nop.i 0 1262};; 1263 1264.pred.rel "mutex",p11,p13 1265{ .mfi 1266 nop.m 0 1267(p11) fms.s0 FR_FirstC = f1, f0, FR_FirstC 1268 nop.i 0 1269} 1270{ .mfi 1271 nop.m 0 1272(p13) fms.s0 FR_FirstC = f1, f0, FR_FirstC 1273 nop.i 0 1274};; 1275 1276{ .mfi 1277 nop.m 0 1278 fma.s1 FR_polyS = FR_r, FR_rsq, FR_c 1279 nop.i 0 1280};; 1281 1282 1283.pred.rel "mutex",p8,p9 1284{ .mfi 1285 nop.m 0 1286(p8) fma.s0 FR_ResultS = FR_FirstS, f1, FR_polyS 1287 nop.i 0 1288} 1289{ .mfi 1290 nop.m 0 1291(p9) fms.s0 FR_ResultS = FR_FirstS, f1, FR_polyS 1292 nop.i 0 1293};; 1294 1295.pred.rel "mutex",p10,p11 1296{ .mfi 1297 nop.m 0 1298(p10) fma.s0 FR_ResultC = FR_FirstC, f1, FR_polyC 1299 nop.i 0 1300} 1301{ .mfi 1302 nop.m 0 1303(p11) fms.s0 FR_ResultC = FR_FirstC, f1, FR_polyC 1304 nop.i 0 1305};; 1306 1307 1308 1309.pred.rel "mutex",p12,p13 1310{ .mfi 1311 nop.m 0 1312(p12) fma.s0 FR_ResultS = FR_FirstC, f1, FR_polyC 1313 nop.i 0 1314} 1315{ .mfi 1316 nop.m 0 1317(p13) fms.s0 FR_ResultS = FR_FirstC, f1, FR_polyC 1318 nop.i 0 1319};; 1320 1321.pred.rel "mutex",p14,p15 1322{ .mfi 1323 nop.m 0 1324(p14) fma.s0 FR_ResultC = FR_FirstS, f1, FR_polyS 1325 nop.i 0 1326} 1327{ .mfb 1328 cmp.eq p10, p0 = 0x1, GR_Cis 1329(p15) fms.s0 FR_ResultC = FR_FirstS, f1, FR_polyS 1330(p10) br.ret.sptk b0 1331};; 1332 1333{ .mmb // exit for sincosl 1334 stfe [sincos_pResSin] = FR_ResultS 1335 stfe [sincos_pResCos] = FR_ResultC 1336 br.ret.sptk b0 1337};; 1338 1339 1340 1341 1342 1343 1344SINCOSL_LARGER_ARG: 1345// 1346// Here if 2^24 <= |x| < 2^63 1347// 1348{ .mfi 1349 ldfe FR_d_1 = [GR_ad_p], 16 // Load d_1 for |x| >= 2^24 path 1350 fma.s1 FR_N_0 = FR_Input_X, FR_Inv_P_0, f0 // N_0 = Arg * Inv_P_0 1351 nop.i 0 1352};; 1353 1354{ .mmi 1355 ldfps FR_Two_to_M14, FR_Neg_Two_to_M14 = [GR_ad_m14] 1356 nop.m 0 1357 nop.i 0 1358};; 1359 1360{ .mfi 1361 ldfe FR_d_2 = [GR_ad_p], 16 // Load d_2 for |x| >= 2^24 path 1362 nop.f 0 1363 nop.i 0 1364};; 1365 1366{ .mfi 1367 nop.m 0 1368 fcvt.fx.s1 FR_N_0_fix = FR_N_0 // N_0_fix = integer part of N_0 1369 nop.i 0 1370};; 1371 1372{ .mfi 1373 nop.m 0 1374 fcvt.xf FR_N_0 = FR_N_0_fix // Make N_0 the integer part 1375 nop.i 0 1376};; 1377 1378{ .mfi 1379 nop.m 0 1380 fnma.s1 FR_ArgPrime = FR_N_0, FR_P_0, FR_Input_X // Arg'=-N_0*P_0+Arg 1381 nop.i 0 1382} 1383{ .mfi 1384 nop.m 0 1385 fma.s1 FR_w = FR_N_0, FR_d_1, f0 // w = N_0 * d_1 1386 nop.i 0 1387};; 1388 1389 1390{ .mfi 1391 nop.m 0 1392 fma.s1 FR_N_float = FR_ArgPrime, FR_Inv_pi_by_2, f0 // N = A' * 2/pi 1393 nop.i 0 1394};; 1395 1396{ .mfi 1397 nop.m 0 1398 fcvt.fx.s1 FR_N_fix = FR_N_float // N_fix is the integer part 1399 nop.i 0 1400};; 1401 1402{ .mfi 1403 nop.m 0 1404 fcvt.xf FR_N_float = FR_N_fix 1405 nop.i 0 1406};; 1407 1408{ .mfi 1409 getf.sig GR_N_Inc = FR_N_fix // N is the integer part of 1410 // the reduced-reduced argument 1411 nop.f 0 1412 nop.i 0 1413};; 1414 1415{ .mfi 1416 nop.m 0 1417 fnma.s1 FR_s = FR_N_float, FR_P_1, FR_ArgPrime // s = -N*P_1 + Arg' 1418 nop.i 0 1419} 1420{ .mfi 1421 nop.m 0 1422 fnma.s1 FR_w = FR_N_float, FR_P_2, FR_w // w = -N*P_2 + w 1423 nop.i 0 1424};; 1425 1426// 1427// For |s| > 2**(-14) r = S + w (r complete) 1428// Else U_hi = N_0 * d_1 1429// 1430{ .mfi 1431 nop.m 0 1432 fcmp.lt.unc.s1 p9, p8 = FR_s, FR_Two_to_M14 1433 nop.i 0 1434};; 1435 1436{ .mfi 1437 nop.m 0 1438(p9) fcmp.gt.s1 p9, p8 = FR_s, FR_Neg_Two_to_M14 // p9 if |s| < 2^-14 1439 nop.i 0 1440};; 1441 1442// 1443// Either S <= -2**(-14) or S >= 2**(-14) 1444// or -2**(-14) < s < 2**(-14) 1445// 1446{ .mfi 1447 nop.m 0 1448(p9) fma.s1 FR_V_hi = FR_N_float, FR_P_2, f0 1449 nop.i 0 1450} 1451{ .mfi 1452 nop.m 0 1453(p9) fma.s1 FR_U_hi = FR_N_0, FR_d_1, f0 1454 nop.i 0 1455};; 1456 1457{ .mfi 1458 nop.m 0 1459(p8) fma.s1 FR_r = FR_s, f1, FR_w 1460 nop.i 0 1461} 1462{ .mfi 1463 nop.m 0 1464(p9) fma.s1 FR_w = FR_N_float, FR_P_3, f0 1465 nop.i 0 1466};; 1467 1468// 1469// We need abs of both U_hi and V_hi - don't 1470// worry about switched sign of V_hi. 1471// 1472// Big s: finish up c = (S - r) + w (c complete) 1473// Case 4: A = U_hi + V_hi 1474// Note: Worry about switched sign of V_hi, so subtract instead of add. 1475// 1476{ .mfi 1477 nop.m 0 1478(p9) fms.s1 FR_A = FR_U_hi, f1, FR_V_hi 1479 nop.i 0 1480} 1481{ .mfi 1482 nop.m 0 1483(p9) fnma.s1 FR_V_lo = FR_N_float, FR_P_2, FR_V_hi 1484 nop.i 0 1485};; 1486 1487{ .mfi 1488 nop.m 0 1489(p9) fmerge.s FR_V_hiabs = f0, FR_V_hi 1490 nop.i 0 1491} 1492{ .mfi 1493 nop.m 0 1494(p9) fms.s1 FR_U_lo = FR_N_0, FR_d_1, FR_U_hi // For small s: U_lo=N_0*d_1-U_hi 1495 nop.i 0 1496};; 1497 1498// 1499// For big s: Is |r| < 2**(-3) 1500// For big s: if p12 set, prepare to branch to Small_R. 1501// For big s: If p13 set, prepare to branch to Normal_R. 1502// 1503{ .mfi 1504 nop.m 0 1505(p9) fmerge.s FR_U_hiabs = f0, FR_U_hi 1506 nop.i 0 1507} 1508{ .mfi 1509 nop.m 0 1510(p8) fms.s1 FR_c = FR_s, f1, FR_r // For big s: c = S - r 1511 nop.i 0 1512};; 1513 1514// 1515// For small S: V_hi = N * P_2 1516// w = N * P_3 1517// Note the product does not include the (-) as in the writeup 1518// so (-) missing for V_hi and w. 1519// 1520{ .mfi 1521 nop.m 0 1522(p8) fcmp.lt.unc.s1 p12, p13 = FR_r, FR_Two_to_M3 1523 nop.i 0 1524};; 1525 1526{ .mfi 1527 nop.m 0 1528(p12) fcmp.gt.s1 p12, p13 = FR_r, FR_Neg_Two_to_M3 1529 nop.i 0 1530};; 1531 1532{ .mfi 1533 nop.m 0 1534(p8) fma.s1 FR_c = FR_c, f1, FR_w 1535 nop.i 0 1536} 1537{ .mfb 1538 nop.m 0 1539(p9) fms.s1 FR_w = FR_N_0, FR_d_2, FR_w 1540(p12) br.cond.spnt SINCOSL_SMALL_R // Branch if |r| < 2^-3 1541 // and 2^24 <= |x| < 2^63 1542};; 1543 1544{ .mib 1545 nop.m 0 1546 nop.i 0 1547(p13) br.cond.sptk SINCOSL_NORMAL_R // Branch if |r| >= 2^-3 1548 // and 2^24 <= |x| < 2^63 1549};; 1550 1551SINCOSL_LARGER_S_TINY: 1552// Here if |s| < 2^-14, and 2^24 <= |x| < 2^63 1553// 1554// Big s: Vector off when |r| < 2**(-3). Recall that p8 will be true. 1555// The remaining stuff is for Case 4. 1556// Small s: V_lo = N * P_2 + U_hi (U_hi is in place of V_hi in writeup) 1557// Note: the (-) is still missing for V_lo. 1558// Small s: w = w + N_0 * d_2 1559// Note: the (-) is now incorporated in w. 1560// 1561{ .mfi 1562 and GR_N_SinCos = 0x1, GR_N_Inc 1563 fcmp.ge.unc.s1 p6, p7 = FR_U_hiabs, FR_V_hiabs 1564 tbit.z p8,p12 = GR_N_Inc, 0 1565} 1566{ .mfi 1567 nop.m 0 1568 fma.s1 FR_t = FR_U_lo, f1, FR_V_lo // C_hi = S + A 1569 nop.i 0 1570};; 1571 1572{ .mfi 1573 sub GR_N_SignS = GR_N_Inc, GR_N_SinCos 1574(p6) fms.s1 FR_a = FR_U_hi, f1, FR_A 1575 add GR_N_SignC = GR_N_Inc, GR_N_SinCos 1576} 1577{ .mfi 1578 nop.m 0 1579(p7) fma.s1 FR_a = FR_V_hi, f1, FR_A 1580 nop.i 0 1581};; 1582 1583{ .mmf 1584 ldfe FR_C_1 = [GR_ad_c], 16 1585 ldfe FR_S_1 = [GR_ad_s], 16 1586 fma.s1 FR_C_hi = FR_s, f1, FR_A 1587};; 1588 1589{ .mmi 1590 ldfe FR_C_2 = [GR_ad_c], 64 1591 ldfe FR_S_2 = [GR_ad_s], 64 1592(p8) tbit.z.unc p10,p11 = GR_N_SignC, 1 1593};; 1594 1595// 1596// r and c have been computed. 1597// Make sure ftz mode is set - should be automatic when using wre 1598// |r| < 2**(-3) 1599// Get [i_0,i_1] - two lsb of N_fix. 1600// 1601// For larger u than v: a = U_hi - A 1602// Else a = V_hi - A (do an add to account for missing (-) on V_hi 1603// 1604{ .mfi 1605 nop.m 0 1606 fma.s1 FR_t = FR_t, f1, FR_w // t = t + w 1607(p8) tbit.z.unc p8,p9 = GR_N_SignS, 1 1608} 1609{ .mfi 1610 nop.m 0 1611(p6) fms.s1 FR_a = FR_a, f1, FR_V_hi 1612 nop.i 0 1613};; 1614 1615// 1616// If u > v: a = (U_hi - A) + V_hi 1617// Else a = (V_hi - A) + U_hi 1618// In each case account for negative missing from V_hi. 1619// 1620{ .mfi 1621 nop.m 0 1622 fms.s1 FR_C_lo = FR_s, f1, FR_C_hi 1623(p12) tbit.z.unc p14,p15 = GR_N_SignC, 1 1624} 1625{ .mfi 1626 nop.m 0 1627(p7) fms.s1 FR_a = FR_U_hi, f1, FR_a 1628 nop.i 0 1629};; 1630 1631{ .mfi 1632 nop.m 0 1633 fma.s1 FR_C_lo = FR_C_lo, f1, FR_A // C_lo = (S - C_hi) + A 1634(p12) tbit.z.unc p12,p13 = GR_N_SignS, 1 1635} 1636{ .mfi 1637 nop.m 0 1638 fma.s1 FR_t = FR_t, f1, FR_a // t = t + a 1639 nop.i 0 1640};; 1641 1642{ .mfi 1643 nop.m 0 1644 fma.s1 FR_r = FR_C_hi, f1, FR_C_lo 1645 nop.i 0 1646};; 1647 1648{ .mfi 1649 nop.m 0 1650 fma.s1 FR_C_lo = FR_C_lo, f1, FR_t // C_lo = C_lo + t 1651 nop.i 0 1652};; 1653 1654 1655{ .mfi 1656 nop.m 0 1657 fma.s1 FR_rsq = FR_r, FR_r, f0 1658 nop.i 0 1659} 1660{ .mfi 1661 nop.m 0 1662 fms.s1 FR_c = FR_C_hi, f1, FR_r 1663 nop.i 0 1664};; 1665 1666{ .mfi 1667 nop.m 0 1668 fma.s1 FR_FirstS = f0, f1, FR_r 1669 nop.i 0 1670} 1671{ .mfi 1672 nop.m 0 1673 fma.s1 FR_FirstC = f0, f1, f1 1674 nop.i 0 1675};; 1676 1677{ .mfi 1678 nop.m 0 1679 fma.s1 FR_polyS = FR_rsq, FR_S_2, FR_S_1 1680 nop.i 0 1681} 1682{ .mfi 1683 nop.m 0 1684 fma.s1 FR_polyC = FR_rsq, FR_C_2, FR_C_1 1685 nop.i 0 1686};; 1687 1688{ .mfi 1689 nop.m 0 1690 fma.s1 FR_r_cubed = FR_rsq, FR_r, f0 1691 nop.i 0 1692} 1693{ .mfi 1694 nop.m 0 1695 fma.s1 FR_c = FR_c, f1, FR_C_lo 1696 nop.i 0 1697};; 1698 1699.pred.rel "mutex",p9,p15 1700{ .mfi 1701 nop.m 0 1702(p9) fms.s0 FR_FirstS = f1, f0, FR_FirstS 1703 nop.i 0 1704} 1705{ .mfi 1706 nop.m 0 1707(p15) fms.s0 FR_FirstS = f1, f0, FR_FirstS 1708 nop.i 0 1709};; 1710 1711.pred.rel "mutex",p11,p13 1712{ .mfi 1713 nop.m 0 1714(p11) fms.s0 FR_FirstC = f1, f0, FR_FirstC 1715 nop.i 0 1716} 1717{ .mfi 1718 nop.m 0 1719(p13) fms.s0 FR_FirstC = f1, f0, FR_FirstC 1720 nop.i 0 1721};; 1722 1723 1724{ .mfi 1725 nop.m 0 1726 fma.s1 FR_polyS = FR_r_cubed, FR_polyS, FR_c 1727 nop.i 0 1728} 1729{ .mfi 1730 nop.m 0 1731 fma.s1 FR_polyC = FR_rsq, FR_polyC, f0 1732 nop.i 0 1733};; 1734 1735 1736 1737.pred.rel "mutex",p8,p9 1738{ .mfi 1739 nop.m 0 1740(p8) fma.s0 FR_ResultS = FR_FirstS, f1, FR_polyS 1741 nop.i 0 1742} 1743{ .mfi 1744 nop.m 0 1745(p9) fms.s0 FR_ResultS = FR_FirstS, f1, FR_polyS 1746 nop.i 0 1747};; 1748 1749.pred.rel "mutex",p10,p11 1750{ .mfi 1751 nop.m 0 1752(p10) fma.s0 FR_ResultC = FR_FirstC, f1, FR_polyC 1753 nop.i 0 1754} 1755{ .mfi 1756 nop.m 0 1757(p11) fms.s0 FR_ResultC = FR_FirstC, f1, FR_polyC 1758 nop.i 0 1759};; 1760 1761 1762 1763.pred.rel "mutex",p12,p13 1764{ .mfi 1765 nop.m 0 1766(p12) fma.s0 FR_ResultS = FR_FirstC, f1, FR_polyC 1767 nop.i 0 1768} 1769{ .mfi 1770 nop.m 0 1771(p13) fms.s0 FR_ResultS = FR_FirstC, f1, FR_polyC 1772 nop.i 0 1773};; 1774 1775.pred.rel "mutex",p14,p15 1776{ .mfi 1777 nop.m 0 1778(p14) fma.s0 FR_ResultC = FR_FirstS, f1, FR_polyS 1779 nop.i 0 1780} 1781{ .mfb 1782 cmp.eq p10, p0 = 0x1, GR_Cis 1783(p15) fms.s0 FR_ResultC = FR_FirstS, f1, FR_polyS 1784(p10) br.ret.sptk b0 1785};; 1786 1787 1788{ .mmb // exit for sincosl 1789 stfe [sincos_pResSin] = FR_ResultS 1790 stfe [sincos_pResCos] = FR_ResultC 1791 br.ret.sptk b0 1792};; 1793 1794 1795 1796SINCOSL_SMALL_R: 1797// 1798// Here if |r| < 2^-3 1799// 1800// Enter with r, c, and N_Inc computed 1801// 1802{ .mfi 1803 nop.m 0 1804 fma.s1 FR_rsq = FR_r, FR_r, f0 // rsq = r * r 1805 nop.i 0 1806};; 1807 1808{ .mmi 1809 ldfe FR_S_5 = [GR_ad_se], -16 // Load S_5 1810 ldfe FR_C_5 = [GR_ad_ce], -16 // Load C_5 1811 nop.i 0 1812};; 1813 1814{ .mmi 1815 ldfe FR_S_4 = [GR_ad_se], -16 // Load S_4 1816 ldfe FR_C_4 = [GR_ad_ce], -16 // Load C_4 1817 nop.i 0 1818};; 1819 1820SINCOSL_SMALL_R_0: 1821// Entry point for 2^-3 < |x| < pi/4 1822SINCOSL_SMALL_R_1: 1823// Entry point for pi/4 < |x| < 2^24 and |r| < 2^-3 1824{ .mfi 1825 ldfe FR_S_3 = [GR_ad_se], -16 // Load S_3 1826 fma.s1 FR_r6 = FR_rsq, FR_rsq, f0 // Z = rsq * rsq 1827 tbit.z p7,p11 = GR_N_Inc, 0 1828} 1829{ .mfi 1830 ldfe FR_C_3 = [GR_ad_ce], -16 // Load C_3 1831 nop.f 0 1832 and GR_N_SinCos = 0x1, GR_N_Inc 1833};; 1834 1835{ .mfi 1836 ldfe FR_S_2 = [GR_ad_se], -16 // Load S_2 1837 fnma.s1 FR_cC = FR_c, FR_r, f0 // c = -c * r 1838 sub GR_N_SignS = GR_N_Inc, GR_N_SinCos 1839} 1840{ .mfi 1841 ldfe FR_C_2 = [GR_ad_ce], -16 // Load C_2 1842 nop.f 0 1843 add GR_N_SignC = GR_N_Inc, GR_N_SinCos 1844};; 1845 1846{ .mmi 1847 ldfe FR_S_1 = [GR_ad_se], -16 // Load S_1 1848 ldfe FR_C_1 = [GR_ad_ce], -16 // Load C_1 1849(p7) tbit.z.unc p9,p10 = GR_N_SignC, 1 1850};; 1851 1852{ .mfi 1853 nop.m 0 1854 fma.s1 FR_r7 = FR_r6, FR_r, f0 // Z = Z * r 1855(p7) tbit.z.unc p7,p8 = GR_N_SignS, 1 1856};; 1857 1858{ .mfi 1859 nop.m 0 1860 fma.s1 FR_poly_loS = FR_rsq, FR_S_5, FR_S_4 // poly_lo=rsq*S_5+S_4 1861(p11) tbit.z.unc p13,p14 = GR_N_SignC, 1 1862} 1863{ .mfi 1864 nop.m 0 1865 fma.s1 FR_poly_loC = FR_rsq, FR_C_5, FR_C_4 // poly_lo=rsq*C_5+C_4 1866 nop.i 0 1867};; 1868 1869{ .mfi 1870 nop.m 0 1871 fma.s1 FR_poly_hiS = FR_rsq, FR_S_2, FR_S_1 // poly_hi=rsq*S_2+S_1 1872(p11) tbit.z.unc p11,p12 = GR_N_SignS, 1 1873} 1874{ .mfi 1875 nop.m 0 1876 fma.s1 FR_poly_hiC = FR_rsq, FR_C_2, FR_C_1 // poly_hi=rsq*C_2+C_1 1877 nop.i 0 1878};; 1879 1880{ .mfi 1881 nop.m 0 1882 fma.s0 FR_FirstS = FR_r, f1, f0 1883 nop.i 0 1884} 1885{ .mfi 1886 nop.m 0 1887 fma.s0 FR_FirstC = f1, f1, f0 1888 nop.i 0 1889};; 1890 1891 1892{ .mfi 1893 nop.m 0 1894 fma.s1 FR_r6 = FR_r6, FR_rsq, f0 1895 nop.i 0 1896} 1897{ .mfi 1898 nop.m 0 1899 fma.s1 FR_r7 = FR_r7, FR_rsq, f0 1900 nop.i 0 1901};; 1902 1903{ .mfi 1904 nop.m 0 1905 fma.s1 FR_poly_loS = FR_rsq, FR_poly_loS, FR_S_3 // p_lo=p_lo*rsq+S_3 1906 nop.i 0 1907} 1908{ .mfi 1909 nop.m 0 1910 fma.s1 FR_poly_loC = FR_rsq, FR_poly_loC, FR_C_3 // p_lo=p_lo*rsq+C_3 1911 nop.i 0 1912};; 1913 1914{ .mfi 1915 nop.m 0 1916 fma.s0 FR_inexact = FR_S_4, FR_S_4, f0 // Dummy op to set inexact 1917 nop.i 0 1918};; 1919 1920{ .mfi 1921 nop.m 0 1922 fma.s1 FR_poly_hiS = FR_poly_hiS, FR_rsq, f0 // p_hi=p_hi*rsq 1923 nop.i 0 1924} 1925{ .mfi 1926 nop.m 0 1927 fma.s1 FR_poly_hiC = FR_poly_hiC, FR_rsq, f0 // p_hi=p_hi*rsq 1928 nop.i 0 1929};; 1930 1931.pred.rel "mutex",p8,p14 1932{ .mfi 1933 nop.m 0 1934(p8) fms.s0 FR_FirstS = f1, f0, FR_FirstS 1935 nop.i 0 1936} 1937{ .mfi 1938 nop.m 0 1939(p14) fms.s0 FR_FirstS = f1, f0, FR_FirstS 1940 nop.i 0 1941};; 1942 1943.pred.rel "mutex",p10,p12 1944{ .mfi 1945 nop.m 0 1946(p10) fms.s0 FR_FirstC = f1, f0, FR_FirstC 1947 nop.i 0 1948} 1949{ .mfi 1950 nop.m 0 1951(p12) fms.s0 FR_FirstC = f1, f0, FR_FirstC 1952 nop.i 0 1953};; 1954 1955{ .mfi 1956 nop.m 0 1957 fma.s1 FR_polyS = FR_r7, FR_poly_loS, FR_cS // poly=Z*poly_lo+c 1958 nop.i 0 1959} 1960{ .mfi 1961 nop.m 0 1962 fma.s1 FR_polyC = FR_r6, FR_poly_loC, FR_cC // poly=Z*poly_lo+c 1963 nop.i 0 1964};; 1965 1966{ .mfi 1967 nop.m 0 1968 fma.s1 FR_poly_hiS = FR_r, FR_poly_hiS, f0 // p_hi=r*p_hi 1969 nop.i 0 1970};; 1971 1972 1973{ .mfi 1974 nop.m 0 1975 fma.s1 FR_polyS = FR_polyS, f1, FR_poly_hiS 1976 nop.i 0 1977} 1978{ .mfi 1979 nop.m 0 1980 fma.s1 FR_polyC = FR_polyC, f1, FR_poly_hiC 1981 nop.i 0 1982};; 1983 1984.pred.rel "mutex",p7,p8 1985{ .mfi 1986 nop.m 0 1987(p7) fma.s0 FR_ResultS = FR_FirstS, f1, FR_polyS 1988 nop.i 0 1989} 1990{ .mfi 1991 nop.m 0 1992(p8) fms.s0 FR_ResultS = FR_FirstS, f1, FR_polyS 1993 nop.i 0 1994};; 1995 1996.pred.rel "mutex",p9,p10 1997{ .mfi 1998 nop.m 0 1999(p9) fma.s0 FR_ResultC = FR_FirstC, f1, FR_polyC 2000 nop.i 0 2001} 2002{ .mfi 2003 nop.m 0 2004(p10) fms.s0 FR_ResultC = FR_FirstC, f1, FR_polyC 2005 nop.i 0 2006};; 2007 2008.pred.rel "mutex",p11,p12 2009{ .mfi 2010 nop.m 0 2011(p11) fma.s0 FR_ResultS = FR_FirstC, f1, FR_polyC 2012 nop.i 0 2013} 2014{ .mfi 2015 nop.m 0 2016(p12) fms.s0 FR_ResultS = FR_FirstC, f1, FR_polyC 2017 nop.i 0 2018};; 2019 2020.pred.rel "mutex",p13,p14 2021{ .mfi 2022 nop.m 0 2023(p13) fma.s0 FR_ResultC = FR_FirstS, f1, FR_polyS 2024 nop.i 0 2025} 2026{ .mfb 2027 cmp.eq p15, p0 = 0x1, GR_Cis 2028(p14) fms.s0 FR_ResultC = FR_FirstS, f1, FR_polyS 2029(p15) br.ret.sptk b0 2030};; 2031 2032 2033{ .mmb // exit for sincosl 2034 stfe [sincos_pResSin] = FR_ResultS 2035 stfe [sincos_pResCos] = FR_ResultC 2036 br.ret.sptk b0 2037};; 2038 2039 2040 2041 2042 2043 2044SINCOSL_NORMAL_R: 2045// 2046// Here if 2^-3 <= |r| < pi/4 2047// THIS IS THE MAIN PATH 2048// 2049// Enter with r, c, and N_Inc having been computed 2050// 2051{ .mfi 2052 ldfe FR_PP_6 = [GR_ad_pp], 16 // Load PP_6 2053 fma.s1 FR_rsq = FR_r, FR_r, f0 // rsq = r * r 2054 nop.i 0 2055} 2056{ .mfi 2057 ldfe FR_QQ_6 = [GR_ad_qq], 16 // Load QQ_6 2058 nop.f 0 2059 nop.i 0 2060};; 2061 2062{ .mmi 2063 ldfe FR_PP_5 = [GR_ad_pp], 16 // Load PP_5 2064 ldfe FR_QQ_5 = [GR_ad_qq], 16 // Load QQ_5 2065 nop.i 0 2066};; 2067 2068 2069 2070SINCOSL_NORMAL_R_0: 2071// Entry for 2^-3 < |x| < pi/4 2072.pred.rel "mutex",p9,p10 2073{ .mmf 2074 ldfe FR_C_1 = [GR_ad_pp], 16 // Load C_1 2075 ldfe FR_S_1 = [GR_ad_qq], 16 // Load S_1 2076 frcpa.s1 FR_r_hi, p6 = f1, FR_r // r_hi = frcpa(r) 2077};; 2078 2079{ .mfi 2080 nop.m 0 2081 fma.s1 FR_polyS = FR_rsq, FR_PP_8, FR_PP_7 // poly = rsq*PP_8+PP_7 2082 nop.i 0 2083} 2084{ .mfi 2085 nop.m 0 2086 fma.s1 FR_polyC = FR_rsq, FR_QQ_8, FR_QQ_7 // poly = rsq*QQ_8+QQ_7 2087 nop.i 0 2088};; 2089 2090{ .mfi 2091 nop.m 0 2092 fma.s1 FR_r_cubed = FR_r, FR_rsq, f0 // rcubed = r * rsq 2093 nop.i 0 2094};; 2095 2096 2097SINCOSL_NORMAL_R_1: 2098// Entry for pi/4 <= |x| < 2^24 2099.pred.rel "mutex",p9,p10 2100{ .mmf 2101 ldfe FR_PP_1 = [GR_ad_pp], 16 // Load PP_1_hi 2102 ldfe FR_QQ_1 = [GR_ad_qq], 16 // Load QQ_1 2103 frcpa.s1 FR_r_hi, p6 = f1, FR_r_hi // r_hi = frpca(frcpa(r)) 2104};; 2105 2106{ .mfi 2107 ldfe FR_PP_4 = [GR_ad_pp], 16 // Load PP_4 2108 fma.s1 FR_polyS = FR_rsq, FR_polyS, FR_PP_6 // poly = rsq*poly+PP_6 2109 and GR_N_SinCos = 0x1, GR_N_Inc 2110} 2111{ .mfi 2112 ldfe FR_QQ_4 = [GR_ad_qq], 16 // Load QQ_4 2113 fma.s1 FR_polyC = FR_rsq, FR_polyC, FR_QQ_6 // poly = rsq*poly+QQ_6 2114 nop.i 0 2115};; 2116 2117{ .mfi 2118 nop.m 0 2119 fma.s1 FR_corrS = FR_C_1, FR_rsq, f0 // corr = C_1 * rsq 2120 sub GR_N_SignS = GR_N_Inc, GR_N_SinCos 2121} 2122{ .mfi 2123 nop.m 0 2124 fma.s1 FR_corrC = FR_S_1, FR_r_cubed, FR_r // corr = S_1 * r^3 + r 2125 add GR_N_SignC = GR_N_Inc, GR_N_SinCos 2126};; 2127 2128{ .mfi 2129 ldfe FR_PP_3 = [GR_ad_pp], 16 // Load PP_3 2130 fma.s1 FR_r_hi_sq = FR_r_hi, FR_r_hi, f0 // r_hi_sq = r_hi * r_hi 2131 tbit.z p7,p11 = GR_N_Inc, 0 2132} 2133{ .mfi 2134 ldfe FR_QQ_3 = [GR_ad_qq], 16 // Load QQ_3 2135 fms.s1 FR_r_lo = FR_r, f1, FR_r_hi // r_lo = r - r_hi 2136 nop.i 0 2137};; 2138 2139{ .mfi 2140 ldfe FR_PP_2 = [GR_ad_pp], 16 // Load PP_2 2141 fma.s1 FR_polyS = FR_rsq, FR_polyS, FR_PP_5 // poly = rsq*poly+PP_5 2142(p7) tbit.z.unc p9,p10 = GR_N_SignC, 1 2143} 2144{ .mfi 2145 ldfe FR_QQ_2 = [GR_ad_qq], 16 // Load QQ_2 2146 fma.s1 FR_polyC = FR_rsq, FR_polyC, FR_QQ_5 // poly = rsq*poly+QQ_5 2147 nop.i 0 2148};; 2149 2150{ .mfi 2151 ldfe FR_PP_1_lo = [GR_ad_pp], 16 // Load PP_1_lo 2152 fma.s1 FR_corrS = FR_corrS, FR_c, FR_c // corr = corr * c + c 2153(p7) tbit.z.unc p7,p8 = GR_N_SignS, 1 2154} 2155{ .mfi 2156 nop.m 0 2157 fnma.s1 FR_corrC = FR_corrC, FR_c, f0 // corr = -corr * c 2158 nop.i 0 2159};; 2160 2161{ .mfi 2162 nop.m 0 2163 fma.s1 FR_U_loS = FR_r, FR_r_hi, FR_r_hi_sq // U_lo = r*r_hi+r_hi_sq 2164(p11) tbit.z.unc p13,p14 = GR_N_SignC, 1 2165} 2166{ .mfi 2167 nop.m 0 2168 fma.s1 FR_U_loC = FR_r_hi, f1, FR_r // U_lo = r_hi + r 2169 nop.i 0 2170};; 2171 2172{ .mfi 2173 nop.m 0 2174 fma.s1 FR_U_hiS = FR_r_hi, FR_r_hi_sq, f0 // U_hi = r_hi*r_hi_sq 2175(p11) tbit.z.unc p11,p12 = GR_N_SignS, 1 2176} 2177{ .mfi 2178 nop.m 0 2179 fma.s1 FR_U_hiC = FR_QQ_1, FR_r_hi_sq, f1 // U_hi = QQ_1*r_hi_sq+1 2180 nop.i 0 2181};; 2182 2183{ .mfi 2184 nop.m 0 2185 fma.s1 FR_polyS = FR_rsq, FR_polyS, FR_PP_4 // poly = poly*rsq+PP_4 2186 nop.i 0 2187} 2188{ .mfi 2189 nop.m 0 2190 fma.s1 FR_polyC = FR_rsq, FR_polyC, FR_QQ_4 // poly = poly*rsq+QQ_4 2191 nop.i 0 2192};; 2193 2194{ .mfi 2195 nop.m 0 2196 fma.s1 FR_U_loS = FR_r, FR_r, FR_U_loS // U_lo = r * r + U_lo 2197 nop.i 0 2198} 2199{ .mfi 2200 nop.m 0 2201 fma.s1 FR_U_loC = FR_r_lo, FR_U_loC, f0 // U_lo = r_lo * U_lo 2202 nop.i 0 2203};; 2204 2205{ .mfi 2206 nop.m 0 2207 fma.s1 FR_U_hiS = FR_PP_1, FR_U_hiS, f0 // U_hi = PP_1 * U_hi 2208 nop.i 0 2209};; 2210 2211{ .mfi 2212 nop.m 0 2213 fma.s1 FR_polyS = FR_rsq, FR_polyS, FR_PP_3 // poly = poly*rsq+PP_3 2214 nop.i 0 2215} 2216{ .mfi 2217 nop.m 0 2218 fma.s1 FR_polyC = FR_rsq, FR_polyC, FR_QQ_3 // poly = poly*rsq+QQ_3 2219 nop.i 0 2220};; 2221 2222{ .mfi 2223 nop.m 0 2224 fma.s1 FR_U_loS = FR_r_lo, FR_U_loS, f0 // U_lo = r_lo * U_lo 2225 nop.i 0 2226} 2227{ .mfi 2228 nop.m 0 2229 fma.s1 FR_U_loC = FR_QQ_1,FR_U_loC, f0 // U_lo = QQ_1 * U_lo 2230 nop.i 0 2231};; 2232 2233{ .mfi 2234 nop.m 0 2235 fma.s1 FR_U_hiS = FR_r, f1, FR_U_hiS // U_hi = r + U_hi 2236 nop.i 0 2237};; 2238 2239{ .mfi 2240 nop.m 0 2241 fma.s1 FR_polyS = FR_rsq, FR_polyS, FR_PP_2 // poly = poly*rsq+PP_2 2242 nop.i 0 2243} 2244{ .mfi 2245 nop.m 0 2246 fma.s1 FR_polyC = FR_rsq, FR_polyC, FR_QQ_2 // poly = poly*rsq+QQ_2 2247 nop.i 0 2248};; 2249 2250{ .mfi 2251 nop.m 0 2252 fma.s1 FR_U_loS = FR_PP_1, FR_U_loS, f0 // U_lo = PP_1 * U_lo 2253 nop.i 0 2254};; 2255 2256{ .mfi 2257 nop.m 0 2258 fma.s1 FR_polyS = FR_rsq, FR_polyS, FR_PP_1_lo // poly =poly*rsq+PP1lo 2259 nop.i 0 2260} 2261{ .mfi 2262 nop.m 0 2263 fma.s1 FR_polyC = FR_rsq, FR_polyC, f0 // poly = poly*rsq 2264 nop.i 0 2265};; 2266 2267 2268.pred.rel "mutex",p8,p14 2269{ .mfi 2270 nop.m 0 2271(p8) fms.s0 FR_U_hiS = f1, f0, FR_U_hiS 2272 nop.i 0 2273} 2274{ .mfi 2275 nop.m 0 2276(p14) fms.s0 FR_U_hiS = f1, f0, FR_U_hiS 2277 nop.i 0 2278};; 2279 2280.pred.rel "mutex",p10,p12 2281{ .mfi 2282 nop.m 0 2283(p10) fms.s0 FR_U_hiC = f1, f0, FR_U_hiC 2284 nop.i 0 2285} 2286{ .mfi 2287 nop.m 0 2288(p12) fms.s0 FR_U_hiC = f1, f0, FR_U_hiC 2289 nop.i 0 2290};; 2291 2292 2293{ .mfi 2294 nop.m 0 2295 fma.s1 FR_VS = FR_U_loS, f1, FR_corrS // V = U_lo + corr 2296 nop.i 0 2297} 2298{ .mfi 2299 nop.m 0 2300 fma.s1 FR_VC = FR_U_loC, f1, FR_corrC // V = U_lo + corr 2301 nop.i 0 2302};; 2303 2304{ .mfi 2305 nop.m 0 2306 fma.s0 FR_inexact = FR_PP_5, FR_PP_4, f0 // Dummy op to set inexact 2307 nop.i 0 2308};; 2309 2310 2311{ .mfi 2312 nop.m 0 2313 fma.s1 FR_polyS = FR_r_cubed, FR_polyS, f0 // poly = poly*r^3 2314 nop.i 0 2315} 2316{ .mfi 2317 nop.m 0 2318 fma.s1 FR_polyC = FR_rsq, FR_polyC, f0 // poly = poly*rsq 2319 nop.i 0 2320};; 2321 2322 2323{ .mfi 2324 nop.m 0 2325 fma.s1 FR_VS = FR_polyS, f1, FR_VS // V = poly + V 2326 nop.i 0 2327} 2328{ .mfi 2329 nop.m 0 2330 fma.s1 FR_VC = FR_polyC, f1, FR_VC // V = poly + V 2331 nop.i 0 2332};; 2333 2334 2335 2336.pred.rel "mutex",p7,p8 2337{ .mfi 2338 nop.m 0 2339(p7) fma.s0 FR_ResultS = FR_U_hiS, f1, FR_VS 2340 nop.i 0 2341} 2342{ .mfi 2343 nop.m 0 2344(p8) fms.s0 FR_ResultS = FR_U_hiS, f1, FR_VS 2345 nop.i 0 2346};; 2347 2348.pred.rel "mutex",p9,p10 2349{ .mfi 2350 nop.m 0 2351(p9) fma.s0 FR_ResultC = FR_U_hiC, f1, FR_VC 2352 nop.i 0 2353} 2354{ .mfi 2355 nop.m 0 2356(p10) fms.s0 FR_ResultC = FR_U_hiC, f1, FR_VC 2357 nop.i 0 2358};; 2359 2360 2361 2362.pred.rel "mutex",p11,p12 2363{ .mfi 2364 nop.m 0 2365(p11) fma.s0 FR_ResultS = FR_U_hiC, f1, FR_VC 2366 nop.i 0 2367} 2368{ .mfi 2369 nop.m 0 2370(p12) fms.s0 FR_ResultS = FR_U_hiC, f1, FR_VC 2371 nop.i 0 2372};; 2373 2374.pred.rel "mutex",p13,p14 2375{ .mfi 2376 nop.m 0 2377(p13) fma.s0 FR_ResultC = FR_U_hiS, f1, FR_VS 2378 nop.i 0 2379} 2380{ .mfb 2381 cmp.eq p15, p0 = 0x1, GR_Cis 2382(p14) fms.s0 FR_ResultC = FR_U_hiS, f1, FR_VS 2383(p15) br.ret.sptk b0 2384};; 2385 2386{ .mmb // exit for sincosl 2387 stfe [sincos_pResSin] = FR_ResultS 2388 stfe [sincos_pResCos] = FR_ResultC 2389 br.ret.sptk b0 2390};; 2391 2392 2393 2394 2395 2396SINCOSL_ZERO: 2397 2398{ .mfi 2399 nop.m 0 2400 fmerge.s FR_ResultS = FR_Input_X, FR_Input_X // If sin, result = input 2401 nop.i 0 2402} 2403{ .mfb 2404 cmp.eq p15, p0 = 0x1, GR_Cis 2405 fma.s0 FR_ResultC = f1, f1, f0 // If cos, result=1.0 2406(p15) br.ret.sptk b0 2407};; 2408 2409{ .mmb // exit for sincosl 2410 stfe [sincos_pResSin] = FR_ResultS 2411 stfe [sincos_pResCos] = FR_ResultC 2412 br.ret.sptk b0 2413};; 2414 2415 2416SINCOSL_DENORMAL: 2417{ .mmb 2418 getf.exp GR_signexp_x = FR_norm_x // Get sign and exponent of x 2419 nop.m 999 2420 br.cond.sptk SINCOSL_COMMON2 // Return to common code 2421} 2422;; 2423 2424 2425SINCOSL_SPECIAL: 2426// 2427// Path for Arg = +/- QNaN, SNaN, Inf 2428// Invalid can be raised. SNaNs 2429// become QNaNs 2430// 2431{ .mfi 2432 cmp.eq p15, p0 = 0x1, GR_Cis 2433 fmpy.s0 FR_ResultS = FR_Input_X, f0 2434 nop.i 0 2435} 2436{ .mfb 2437 nop.m 0 2438 fmpy.s0 FR_ResultC = FR_Input_X, f0 2439(p15) br.ret.sptk b0 2440};; 2441 2442{ .mmb // exit for sincosl 2443 stfe [sincos_pResSin] = FR_ResultS 2444 stfe [sincos_pResCos] = FR_ResultC 2445 br.ret.sptk b0 2446};; 2447 2448GLOBAL_LIBM_END(__libm_sincosl) 2449 2450 2451// ******************************************************************* 2452// ******************************************************************* 2453// ******************************************************************* 2454// 2455// Special Code to handle very large argument case. 2456// Call int __libm_pi_by_2_reduce(x,r,c) for |arguments| >= 2**63 2457// The interface is custom: 2458// On input: 2459// (Arg or x) is in f8 2460// On output: 2461// r is in f8 2462// c is in f9 2463// N is in r8 2464// Be sure to allocate at least 2 GP registers as output registers for 2465// __libm_pi_by_2_reduce. This routine uses r62-63. These are used as 2466// scratch registers within the __libm_pi_by_2_reduce routine (for speed). 2467// 2468// We know also that __libm_pi_by_2_reduce preserves f10-15, f71-127. We 2469// use this to eliminate save/restore of key fp registers in this calling 2470// function. 2471// 2472// ******************************************************************* 2473// ******************************************************************* 2474// ******************************************************************* 2475 2476LOCAL_LIBM_ENTRY(__libm_callout) 2477SINCOSL_ARG_TOO_LARGE: 2478.prologue 2479{ .mfi 2480 nop.f 0 2481.save ar.pfs,GR_SAVE_PFS 2482 mov GR_SAVE_PFS=ar.pfs // Save ar.pfs 2483};; 2484 2485{ .mmi 2486 setf.exp FR_Two_to_M3 = GR_exp_2_to_m3 // Form 2^-3 2487 mov GR_SAVE_GP=gp // Save gp 2488.save b0, GR_SAVE_B0 2489 mov GR_SAVE_B0=b0 // Save b0 2490};; 2491 2492.body 2493// 2494// Call argument reduction with x in f8 2495// Returns with N in r8, r in f8, c in f9 2496// Assumes f71-127 are preserved across the call 2497// 2498{ .mib 2499 setf.exp FR_Neg_Two_to_M3 = GR_exp_m2_to_m3 // Form -(2^-3) 2500 nop.i 0 2501 br.call.sptk b0=__libm_pi_by_2_reduce# 2502};; 2503 2504{ .mfi 2505 mov GR_N_Inc = r8 2506 fcmp.lt.unc.s1 p6, p0 = FR_r, FR_Two_to_M3 2507 mov b0 = GR_SAVE_B0 // Restore return address 2508};; 2509 2510{ .mfi 2511 mov gp = GR_SAVE_GP // Restore gp 2512(p6) fcmp.gt.unc.s1 p6, p0 = FR_r, FR_Neg_Two_to_M3 2513 mov ar.pfs = GR_SAVE_PFS // Restore ar.pfs 2514};; 2515 2516{ .mbb 2517 nop.m 0 2518(p6) br.cond.spnt SINCOSL_SMALL_R // Branch if |r|< 2^-3 for |x| >= 2^63 2519 br.cond.sptk SINCOSL_NORMAL_R // Branch if |r|>=2^-3 for |x| >= 2^63 2520};; 2521 2522LOCAL_LIBM_END(__libm_callout) 2523 2524.type __libm_pi_by_2_reduce#,@function 2525.global __libm_pi_by_2_reduce# 2526