1 /* Copyright (c) 2018, Oracle and/or its affiliates. All rights reserved. 2 * Copyright (c) 2018, Cavium. All rights reserved. (By BELLSOFT) 3 * DO NOT ALTER OR REMOVE COPYRIGHT NOTICES OR THIS FILE HEADER. 4 * 5 * This code is free software; you can redistribute it and/or modify it 6 * under the terms of the GNU General Public License version 2 only, as 7 * published by the Free Software Foundation. 8 * 9 * This code is distributed in the hope that it will be useful, but WITHOUT 10 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or 11 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License 12 * version 2 for more details (a copy is included in the LICENSE file that 13 * accompanied this code). 14 * 15 * You should have received a copy of the GNU General Public License version 16 * 2 along with this work; if not, write to the Free Software Foundation, 17 * Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA. 18 * 19 * Please contact Oracle, 500 Oracle Parkway, Redwood Shores, CA 94065 USA 20 * or visit www.oracle.com if you need additional information or have any 21 * questions. 22 * 23 */ 24 25 #include "precompiled.hpp" 26 #include "asm/assembler.hpp" 27 #include "asm/assembler.inline.hpp" 28 #include "runtime/stubRoutines.hpp" 29 #include "macroAssembler_aarch64.hpp" 30 31 // The following code is a optimized version of fdlibm sin/cos implementation 32 // (C code is in share/runtime/sharedRuntimeTrig.cpp) adapted for AARCH64. 33 34 // Please refer to sin/cos approximation via polynomial and 35 // trigonometric argument reduction techniques to the following literature: 36 // 37 // [1] Muller, Jean-Michel, Nicolas Brisebarre, Florent De Dinechin, 38 // Claude-Pierre Jeannerod, Vincent Lefevre, Guillaume Melquiond, 39 // Nathalie Revol, Damien Stehlé, and Serge Torres: 40 // Handbook of floating-point arithmetic. 41 // Springer Science & Business Media, 2009. 42 // [2] K. C. Ng 43 // Argument Reduction for Huge Arguments: Good to the Last Bit 44 // July 13, 1992, SunPro 45 // 46 // HOW TO READ THIS CODE: 47 // This code consists of several functions. Each function has following header: 48 // 1) Description 49 // 2) C-pseudo code with differences from fdlibm marked by comments starting 50 // with "NOTE". Check unmodified fdlibm code in 51 // share/runtime/SharedRuntimeTrig.cpp 52 // 3) Brief textual description of changes between fdlibm and current 53 // implementation along with optimization notes (if applicable) 54 // 4) Assumptions, input and output 55 // 5) (Optional) additional notes about intrinsic implementation 56 // Each function is separated in blocks which follow the pseudo-code structure 57 // 58 // HIGH-LEVEL ALGORITHM DESCRIPTION: 59 // - entry point: generate_dsin_dcos(...); 60 // - check corner cases: NaN, INF, tiny argument. 61 // - check if |x| < Pi/4. Then approximate sin/cos via polynomial (kernel_sin/kernel_cos) 62 // -- else proceed to argument reduction routine (__ieee754_rem_pio2) and 63 // use reduced argument to get result via kernel_sin/kernel_cos 64 // 65 // HIGH-LEVEL CHANGES BETWEEN INTRINSICS AND FDLIBM: 66 // 1) two_over_pi table fdlibm representation is int[], while intrinsic version 67 // has these int values converted to double representation to load converted 68 // double values directly (see stubRoutines_aarch4::_two_over_pi) 69 // 2) Several loops are unrolled and vectorized: see comments in code after 70 // labels: SKIP_F_LOAD, RECOMP_FOR1_CHECK, RECOMP_FOR2 71 // 3) fdlibm npio2_hw table now has "prefix" with constants used in 72 // calculation. These constants are loaded from npio2_hw table instead of 73 // constructing it in code (see stubRoutines_aarch64.cpp) 74 // 4) Polynomial coefficients for sin and cos are moved to table sin_coef 75 // and cos_coef to use the same optimization as in 3). It allows to load most of 76 // required constants via single instruction 77 // 78 // 79 // 80 ///* __ieee754_rem_pio2(x,y) 81 // * 82 // * returns the remainder of x rem pi/2 in y[0]+y[1] (i.e. like x div pi/2) 83 // * x is input argument, y[] is hi and low parts of reduced argument (x) 84 // * uses __kernel_rem_pio2() 85 // */ 86 // // use tables(see stubRoutines_aarch64.cpp): two_over_pi and modified npio2_hw 87 // 88 // BEGIN __ieee754_rem_pio2 PSEUDO CODE 89 // 90 //static int __ieee754_rem_pio2(double x, double *y) { 91 // double z,w,t,r,fn; 92 // double tx[3]; 93 // int e0,i,j,nx,n,ix,hx,i0; 94 // 95 // i0 = ((*(int*)&two24A)>>30)^1; /* high word index */ 96 // hx = *(i0+(int*)&x); /* high word of x */ 97 // ix = hx&0x7fffffff; 98 // if(ix<0x4002d97c) { /* |x| < 3pi/4, special case with n=+-1 */ 99 // if(hx>0) { 100 // z = x - pio2_1; 101 // if(ix!=0x3ff921fb) { /* 33+53 bit pi is good enough */ 102 // y[0] = z - pio2_1t; 103 // y[1] = (z-y[0])-pio2_1t; 104 // } else { /* near pi/2, use 33+33+53 bit pi */ 105 // z -= pio2_2; 106 // y[0] = z - pio2_2t; 107 // y[1] = (z-y[0])-pio2_2t; 108 // } 109 // return 1; 110 // } else { /* negative x */ 111 // z = x + pio2_1; 112 // if(ix!=0x3ff921fb) { /* 33+53 bit pi is good enough */ 113 // y[0] = z + pio2_1t; 114 // y[1] = (z-y[0])+pio2_1t; 115 // } else { /* near pi/2, use 33+33+53 bit pi */ 116 // z += pio2_2; 117 // y[0] = z + pio2_2t; 118 // y[1] = (z-y[0])+pio2_2t; 119 // } 120 // return -1; 121 // } 122 // } 123 // if(ix<=0x413921fb) { /* |x| ~<= 2^19*(pi/2), medium size */ 124 // t = fabsd(x); 125 // n = (int) (t*invpio2+half); 126 // fn = (double)n; 127 // r = t-fn*pio2_1; 128 // w = fn*pio2_1t; /* 1st round good to 85 bit */ 129 // // NOTE: y[0] = r-w; is moved from if/else below to be before "if" 130 // y[0] = r-w; 131 // if(n<32&&ix!=npio2_hw[n-1]) { 132 // // y[0] = r-w; /* quick check no cancellation */ // NOTE: moved earlier 133 // } else { 134 // j = ix>>20; 135 // // y[0] = r-w; // NOTE: moved earlier 136 // i = j-(((*(i0+(int*)&y[0]))>>20)&0x7ff); 137 // if(i>16) { /* 2nd iteration needed, good to 118 */ 138 // t = r; 139 // w = fn*pio2_2; 140 // r = t-w; 141 // w = fn*pio2_2t-((t-r)-w); 142 // y[0] = r-w; 143 // i = j-(((*(i0+(int*)&y[0]))>>20)&0x7ff); 144 // if(i>49) { /* 3rd iteration need, 151 bits acc */ 145 // t = r; /* will cover all possible cases */ 146 // w = fn*pio2_3; 147 // r = t-w; 148 // w = fn*pio2_3t-((t-r)-w); 149 // y[0] = r-w; 150 // } 151 // } 152 // } 153 // y[1] = (r-y[0])-w; 154 // if(hx<0) {y[0] = -y[0]; y[1] = -y[1]; return -n;} 155 // else return n; 156 // } 157 // /* 158 // * all other (large) arguments 159 // */ 160 // // NOTE: this check is removed, because it was checked in dsin/dcos 161 // // if(ix>=0x7ff00000) { /* x is inf or NaN */ 162 // // y[0]=y[1]=x-x; return 0; 163 // // } 164 // /* set z = scalbn(|x|,ilogb(x)-23) */ 165 // *(1-i0+(int*)&z) = *(1-i0+(int*)&x); 166 // e0 = (ix>>20)-1046; /* e0 = ilogb(z)-23; */ 167 // *(i0+(int*)&z) = ix - (e0<<20); 168 // 169 // // NOTE: "for" loop below in unrolled. See comments in asm code 170 // for(i=0;i<2;i++) { 171 // tx[i] = (double)((int)(z)); 172 // z = (z-tx[i])*two24A; 173 // } 174 // 175 // tx[2] = z; 176 // nx = 3; 177 // 178 // // NOTE: while(tx[nx-1]==zeroA) nx--; is unrolled. See comments in asm code 179 // while(tx[nx-1]==zeroA) nx--; /* skip zero term */ 180 // 181 // n = __kernel_rem_pio2(tx,y,e0,nx,2,two_over_pi); 182 // if(hx<0) {y[0] = -y[0]; y[1] = -y[1]; return -n;} 183 // return n; 184 //} 185 // 186 // END __ieee754_rem_pio2 PSEUDO CODE 187 // 188 // Changes between fdlibm and intrinsic for __ieee754_rem_pio2: 189 // 1. INF/NaN check for huge argument is removed in comparison with fdlibm 190 // code, because this check is already done in dcos/dsin code 191 // 2. Most constants are now loaded from table instead of direct initialization 192 // 3. Two loops are unrolled 193 // Assumptions: 194 // 1. Assume |X| >= PI/4 195 // 2. Assume rscratch1 = 0x3fe921fb00000000 (~ PI/4) 196 // 3. Assume ix = r3 197 // Input and output: 198 // 1. Input: X = r0 199 // 2. Return n in r2, y[0] == y0 == v4, y[1] == y1 == v5 200 // NOTE: general purpose register names match local variable names in C code 201 // NOTE: fpu registers are actively reused. See comments in code about their usage 202 void MacroAssembler::generate__ieee754_rem_pio2(address npio2_hw, 203 address two_over_pi, address pio2) { 204 const int64_t PIO2_1t = 0x3DD0B4611A626331UL; 205 const int64_t PIO2_2 = 0x3DD0B4611A600000UL; 206 const int64_t PIO2_2t = 0x3BA3198A2E037073UL; 207 Label X_IS_NEGATIVE, X_IS_MEDIUM_OR_LARGE, X_IS_POSITIVE_LONG_PI, LARGE_ELSE, 208 REDUCTION_DONE, X_IS_MEDIUM_BRANCH_DONE, X_IS_LARGE, NX_SET, 209 X_IS_NEGATIVE_LONG_PI; 210 Register X = r0, n = r2, ix = r3, jv = r4, tmp5 = r5, jx = r6, 211 tmp3 = r7, iqBase = r10, ih = r11, i = r17; 212 // initializing constants first 213 // rscratch1 = 0x3fe921fb00000000 (see assumptions) 214 movk(rscratch1, 0x3ff9, 48); // was 0x3fe921fb0..0 now it's 0x3ff921fb0..0 215 mov(rscratch2, 0x4002d97c); // 3*PI/4 high word 216 movk(rscratch1, 0x5440, 16); // now rscratch1 == PIO2_1 217 fmovd(v1, rscratch1); // v1 = PIO2_1 218 cmp(rscratch2, ix); 219 br(LE, X_IS_MEDIUM_OR_LARGE); 220 221 block_comment("if(ix<0x4002d97c) {... /* |x| ~< 3pi/4 */ "); { 222 cmp(X, zr); 223 br(LT, X_IS_NEGATIVE); 224 225 block_comment("if(hx>0) {"); { 226 fsubd(v2, v0, v1); // v2 = z = x - pio2_1 227 cmp(ix, rscratch1, LSR, 32); 228 mov(n, 1); 229 br(EQ, X_IS_POSITIVE_LONG_PI); 230 231 block_comment("case: hx > 0 && ix!=0x3ff921fb {"); { /* 33+53 bit pi is good enough */ 232 mov(rscratch2, PIO2_1t); 233 fmovd(v27, rscratch2); 234 fsubd(v4, v2, v27); // v4 = y[0] = z - pio2_1t; 235 fsubd(v5, v2, v4); 236 fsubd(v5, v5, v27); // v5 = y[1] = (z-y[0])-pio2_1t 237 b(REDUCTION_DONE); 238 } 239 240 block_comment("case: hx > 0 &*& ix==0x3ff921fb {"); { /* near pi/2, use 33+33+53 bit pi */ 241 bind(X_IS_POSITIVE_LONG_PI); 242 mov(rscratch1, PIO2_2); 243 mov(rscratch2, PIO2_2t); 244 fmovd(v27, rscratch1); 245 fmovd(v6, rscratch2); 246 fsubd(v2, v2, v27); // z-= pio2_2 247 fsubd(v4, v2, v6); // y[0] = z - pio2_2t 248 fsubd(v5, v2, v4); 249 fsubd(v5, v5, v6); // v5 = (z - y[0]) - pio2_2t 250 b(REDUCTION_DONE); 251 } 252 } 253 254 block_comment("case: hx <= 0)"); { 255 bind(X_IS_NEGATIVE); 256 faddd(v2, v0, v1); // v2 = z = x + pio2_1 257 cmp(ix, rscratch1, LSR, 32); 258 mov(n, -1); 259 br(EQ, X_IS_NEGATIVE_LONG_PI); 260 261 block_comment("case: hx <= 0 && ix!=0x3ff921fb) {"); { /* 33+53 bit pi is good enough */ 262 mov(rscratch2, PIO2_1t); 263 fmovd(v27, rscratch2); 264 faddd(v4, v2, v27); // v4 = y[0] = z + pio2_1t; 265 fsubd(v5, v2, v4); 266 faddd(v5, v5, v27); // v5 = y[1] = (z-y[0]) + pio2_1t 267 b(REDUCTION_DONE); 268 } 269 270 block_comment("case: hx <= 0 && ix==0x3ff921fb"); { /* near pi/2, use 33+33+53 bit pi */ 271 bind(X_IS_NEGATIVE_LONG_PI); 272 mov(rscratch1, PIO2_2); 273 mov(rscratch2, PIO2_2t); 274 fmovd(v27, rscratch1); 275 fmovd(v6, rscratch2); 276 faddd(v2, v2, v27); // z += pio2_2 277 faddd(v4, v2, v6); // y[0] = z + pio2_2t 278 fsubd(v5, v2, v4); 279 faddd(v5, v5, v6); // v5 = (z - y[0]) + pio2_2t 280 b(REDUCTION_DONE); 281 } 282 } 283 } 284 bind(X_IS_MEDIUM_OR_LARGE); 285 mov(rscratch1, 0x413921fb); 286 cmp(ix, rscratch1); // ix < = 0x413921fb ? 287 br(GT, X_IS_LARGE); 288 289 block_comment("|x| ~<= 2^19*(pi/2), medium size"); { 290 lea(ih, ExternalAddress(npio2_hw)); 291 ld1(v4, v5, v6, v7, T1D, ih); 292 fabsd(v31, v0); // v31 = t = |x| 293 add(ih, ih, 64); 294 fmaddd(v2, v31, v5, v4); // v2 = t * invpio2 + half (invpio2 = 53 bits of 2/pi, half = 0.5) 295 fcvtzdw(n, v2); // n = (int) v2 296 frintzd(v2, v2); 297 fmsubd(v3, v2, v6, v31); // v3 = r = t - fn * pio2_1 298 fmuld(v26, v2, v7); // v26 = w = fn * pio2_1t 299 fsubd(v4, v3, v26); // y[0] = r - w. Calculated before branch 300 cmp(n, (u1)32); 301 br(GT, LARGE_ELSE); 302 subw(tmp5, n, 1); // tmp5 = n - 1 303 ldrw(jv, Address(ih, tmp5, Address::lsl(2))); 304 cmp(ix, jv); 305 br(NE, X_IS_MEDIUM_BRANCH_DONE); 306 307 block_comment("else block for if(n<32&&ix!=npio2_hw[n-1])"); { 308 bind(LARGE_ELSE); 309 fmovd(jx, v4); 310 lsr(tmp5, ix, 20); // j = ix >> 20 311 lsl(jx, jx, 1); 312 sub(tmp3, tmp5, jx, LSR, 32 + 20 + 1); // r7 = j-(((*(i0+(int*)&y[0]))>>20)&0x7ff); 313 314 block_comment("if(i>16)"); { 315 cmp(tmp3, (u1)16); 316 br(LE, X_IS_MEDIUM_BRANCH_DONE); 317 // i > 16. 2nd iteration needed 318 ldpd(v6, v7, Address(ih, -32)); 319 fmovd(v28, v3); // t = r 320 fmuld(v29, v2, v6); // w = v29 = fn * pio2_2 321 fsubd(v3, v28, v29); // r = t - w 322 fsubd(v31, v28, v3); // v31 = (t - r) 323 fsubd(v31, v29, v31); // v31 = w - (t - r) = - ((t - r) - w) 324 fmaddd(v26, v2, v7, v31); // v26 = w = fn*pio2_2t - ((t - r) - w) 325 fsubd(v4, v3, v26); // y[0] = r - w 326 fmovd(jx, v4); 327 lsl(jx, jx, 1); 328 sub(tmp3, tmp5, jx, LSR, 32 + 20 + 1); // r7 = j-(((*(i0+(int*)&y[0]))>>20)&0x7ff); 329 330 block_comment("if(i>49)"); { 331 cmp(tmp3, (u1)49); 332 br(LE, X_IS_MEDIUM_BRANCH_DONE); 333 // 3rd iteration need, 151 bits acc 334 ldpd(v6, v7, Address(ih, -16)); 335 fmovd(v28, v3); // save "r" 336 fmuld(v29, v2, v6); // v29 = fn * pio2_3 337 fsubd(v3, v28, v29); // r = r - w 338 fsubd(v31, v28, v3); // v31 = (t - r) 339 fsubd(v31, v29, v31); // v31 = w - (t - r) = - ((t - r) - w) 340 fmaddd(v26, v2, v7, v31); // v26 = w = fn*pio2_3t - ((t - r) - w) 341 fsubd(v4, v3, v26); // y[0] = r - w 342 } 343 } 344 } 345 block_comment("medium x tail"); { 346 bind(X_IS_MEDIUM_BRANCH_DONE); 347 fsubd(v5, v3, v4); // v5 = y[1] = (r - y[0]) 348 fsubd(v5, v5, v26); // v5 = y[1] = (r - y[0]) - w 349 cmp(X, zr); 350 br(GT, REDUCTION_DONE); 351 fnegd(v4, v4); 352 negw(n, n); 353 fnegd(v5, v5); 354 b(REDUCTION_DONE); 355 } 356 } 357 358 block_comment("all other (large) arguments"); { 359 bind(X_IS_LARGE); 360 lsr(rscratch1, ix, 20); // ix >> 20 361 movz(tmp5, 0x4170, 48); 362 subw(rscratch1, rscratch1, 1046); // e0 363 fmovd(v10, tmp5); // init two24A value 364 subw(jv, ix, rscratch1, LSL, 20); // ix - (e0<<20) 365 lsl(jv, jv, 32); 366 subw(rscratch2, rscratch1, 3); 367 bfm(jv, X, 0, 31); // jv = z 368 movw(i, 24); 369 fmovd(v26, jv); // v26 = z 370 371 block_comment("unrolled for(i=0;i<2;i++) {tx[i] = (double)((int)(z));z = (z-tx[i])*two24A;}"); { 372 // tx[0,1,2] = v6,v7,v26 373 frintzd(v6, v26); // v6 = (double)((int)v26) 374 sdivw(jv, rscratch2, i); // jv = (e0 - 3)/24 375 fsubd(v26, v26, v6); 376 sub(sp, sp, 560); 377 fmuld(v26, v26, v10); 378 frintzd(v7, v26); // v7 = (double)((int)v26) 379 movw(jx, 2); // calculate jx as nx - 1, which is initially 2. Not a part of unrolled loop 380 fsubd(v26, v26, v7); 381 } 382 383 block_comment("nx calculation with unrolled while(tx[nx-1]==zeroA) nx--;"); { 384 fcmpd(v26, 0.0); // if NE then jx == 2. else it's 1 or 0 385 add(iqBase, sp, 480); // base of iq[] 386 fmuld(v3, v26, v10); 387 br(NE, NX_SET); 388 fcmpd(v7, 0.0); // v7 == 0 => jx = 0. Else jx = 1 389 csetw(jx, NE); 390 } 391 bind(NX_SET); 392 generate__kernel_rem_pio2(two_over_pi, pio2); 393 // now we have y[0] = v4, y[1] = v5 and n = r2 394 cmp(X, zr); 395 br(GE, REDUCTION_DONE); 396 fnegd(v4, v4); 397 fnegd(v5, v5); 398 negw(n, n); 399 } 400 bind(REDUCTION_DONE); 401 } 402 403 ///* 404 // * __kernel_rem_pio2(x,y,e0,nx,prec,ipio2) 405 // * double x[],y[]; int e0,nx,prec; int ipio2[]; 406 // * 407 // * __kernel_rem_pio2 return the last three digits of N with 408 // * y = x - N*pi/2 409 // * so that |y| < pi/2. 410 // * 411 // * The method is to compute the integer (mod 8) and fraction parts of 412 // * (2/pi)*x without doing the full multiplication. In general we 413 // * skip the part of the product that are known to be a huge integer ( 414 // * more accurately, = 0 mod 8 ). Thus the number of operations are 415 // * independent of the exponent of the input. 416 // * 417 // * NOTE: 2/pi int representation is converted to double 418 // * // (2/pi) is represented by an array of 24-bit integers in ipio2[]. 419 // * 420 // * Input parameters: 421 // * x[] The input value (must be positive) is broken into nx 422 // * pieces of 24-bit integers in double precision format. 423 // * x[i] will be the i-th 24 bit of x. The scaled exponent 424 // * of x[0] is given in input parameter e0 (i.e., x[0]*2^e0 425 // * match x's up to 24 bits. 426 // * 427 // * Example of breaking a double positive z into x[0]+x[1]+x[2]: 428 // * e0 = ilogb(z)-23 429 // * z = scalbn(z,-e0) 430 // * for i = 0,1,2 431 // * x[i] = floor(z) 432 // * z = (z-x[i])*2**24 433 // * 434 // * 435 // * y[] ouput result in an array of double precision numbers. 436 // * The dimension of y[] is: 437 // * 24-bit precision 1 438 // * 53-bit precision 2 439 // * 64-bit precision 2 440 // * 113-bit precision 3 441 // * The actual value is the sum of them. Thus for 113-bit 442 // * precsion, one may have to do something like: 443 // * 444 // * long double t,w,r_head, r_tail; 445 // * t = (long double)y[2] + (long double)y[1]; 446 // * w = (long double)y[0]; 447 // * r_head = t+w; 448 // * r_tail = w - (r_head - t); 449 // * 450 // * e0 The exponent of x[0] 451 // * 452 // * nx dimension of x[] 453 // * 454 // * prec an interger indicating the precision: 455 // * 0 24 bits (single) 456 // * 1 53 bits (double) 457 // * 2 64 bits (extended) 458 // * 3 113 bits (quad) 459 // * 460 // * NOTE: ipio2[] array below is converted to double representation 461 // * //ipio2[] 462 // * // integer array, contains the (24*i)-th to (24*i+23)-th 463 // * // bit of 2/pi after binary point. The corresponding 464 // * // floating value is 465 // * 466 // * ipio2[i] * 2^(-24(i+1)). 467 // * 468 // * Here is the description of some local variables: 469 // * 470 // * jk jk+1 is the initial number of terms of ipio2[] needed 471 // * in the computation. The recommended value is 2,3,4, 472 // * 6 for single, double, extended,and quad. 473 // * 474 // * jz local integer variable indicating the number of 475 // * terms of ipio2[] used. 476 // * 477 // * jx nx - 1 478 // * 479 // * jv index for pointing to the suitable ipio2[] for the 480 // * computation. In general, we want 481 // * ( 2^e0*x[0] * ipio2[jv-1]*2^(-24jv) )/8 482 // * is an integer. Thus 483 // * e0-3-24*jv >= 0 or (e0-3)/24 >= jv 484 // * Hence jv = max(0,(e0-3)/24). 485 // * 486 // * jp jp+1 is the number of terms in PIo2[] needed, jp = jk. 487 // * 488 // * q[] double array with integral value, representing the 489 // * 24-bits chunk of the product of x and 2/pi. 490 // * 491 // * q0 the corresponding exponent of q[0]. Note that the 492 // * exponent for q[i] would be q0-24*i. 493 // * 494 // * PIo2[] double precision array, obtained by cutting pi/2 495 // * into 24 bits chunks. 496 // * 497 // * f[] ipio2[] in floating point 498 // * 499 // * iq[] integer array by breaking up q[] in 24-bits chunk. 500 // * 501 // * fq[] final product of x*(2/pi) in fq[0],..,fq[jk] 502 // * 503 // * ih integer. If >0 it indicates q[] is >= 0.5, hence 504 // * it also indicates the *sign* of the result. 505 // * 506 // */ 507 // 508 // Use PIo2 table(see stubRoutines_aarch64.cpp) 509 // 510 // BEGIN __kernel_rem_pio2 PSEUDO CODE 511 // 512 //static int __kernel_rem_pio2(double *x, double *y, int e0, int nx, int prec, /* NOTE: converted to double */ const double *ipio2 // const int *ipio2) { 513 // int jz,jx,jv,jp,jk,carry,n,iq[20],i,j,k,m,q0,ih; 514 // double z,fw,f[20],fq[20],q[20]; 515 // 516 // /* initialize jk*/ 517 // // jk = init_jk[prec]; // NOTE: prec==2 for double. jk is always 4. 518 // jp = jk; // NOTE: always 4 519 // 520 // /* determine jx,jv,q0, note that 3>q0 */ 521 // jx = nx-1; 522 // jv = (e0-3)/24; if(jv<0) jv=0; 523 // q0 = e0-24*(jv+1); 524 // 525 // /* set up f[0] to f[jx+jk] where f[jx+jk] = ipio2[jv+jk] */ 526 // j = jv-jx; m = jx+jk; 527 // 528 // // NOTE: split into two for-loops: one with zeroB and one with ipio2[j]. It 529 // // allows the use of wider loads/stores 530 // for(i=0;i<=m;i++,j++) f[i] = (j<0)? zeroB : /* NOTE: converted to double */ ipio2[j]; //(double) ipio2[j]; 531 // 532 // // NOTE: unrolled and vectorized "for". See comments in asm code 533 // /* compute q[0],q[1],...q[jk] */ 534 // for (i=0;i<=jk;i++) { 535 // for(j=0,fw=0.0;j<=jx;j++) fw += x[j]*f[jx+i-j]; q[i] = fw; 536 // } 537 // 538 // jz = jk; 539 //recompute: 540 // /* distill q[] into iq[] reversingly */ 541 // for(i=0,j=jz,z=q[jz];j>0;i++,j--) { 542 // fw = (double)((int)(twon24* z)); 543 // iq[i] = (int)(z-two24B*fw); 544 // z = q[j-1]+fw; 545 // } 546 // 547 // /* compute n */ 548 // z = scalbnA(z,q0); /* actual value of z */ 549 // z -= 8.0*floor(z*0.125); /* trim off integer >= 8 */ 550 // n = (int) z; 551 // z -= (double)n; 552 // ih = 0; 553 // if(q0>0) { /* need iq[jz-1] to determine n */ 554 // i = (iq[jz-1]>>(24-q0)); n += i; 555 // iq[jz-1] -= i<<(24-q0); 556 // ih = iq[jz-1]>>(23-q0); 557 // } 558 // else if(q0==0) ih = iq[jz-1]>>23; 559 // else if(z>=0.5) ih=2; 560 // 561 // if(ih>0) { /* q > 0.5 */ 562 // n += 1; carry = 0; 563 // for(i=0;i<jz ;i++) { /* compute 1-q */ 564 // j = iq[i]; 565 // if(carry==0) { 566 // if(j!=0) { 567 // carry = 1; iq[i] = 0x1000000- j; 568 // } 569 // } else iq[i] = 0xffffff - j; 570 // } 571 // if(q0>0) { /* rare case: chance is 1 in 12 */ 572 // switch(q0) { 573 // case 1: 574 // iq[jz-1] &= 0x7fffff; break; 575 // case 2: 576 // iq[jz-1] &= 0x3fffff; break; 577 // } 578 // } 579 // if(ih==2) { 580 // z = one - z; 581 // if(carry!=0) z -= scalbnA(one,q0); 582 // } 583 // } 584 // 585 // /* check if recomputation is needed */ 586 // if(z==zeroB) { 587 // j = 0; 588 // for (i=jz-1;i>=jk;i--) j |= iq[i]; 589 // if(j==0) { /* need recomputation */ 590 // for(k=1;iq[jk-k]==0;k++); /* k = no. of terms needed */ 591 // 592 // for(i=jz+1;i<=jz+k;i++) { /* add q[jz+1] to q[jz+k] */ 593 // f[jx+i] = /* NOTE: converted to double */ ipio2[jv+i]; //(double) ipio2[jv+i]; 594 // for(j=0,fw=0.0;j<=jx;j++) fw += x[j]*f[jx+i-j]; 595 // q[i] = fw; 596 // } 597 // jz += k; 598 // goto recompute; 599 // } 600 // } 601 // 602 // /* chop off zero terms */ 603 // if(z==0.0) { 604 // jz -= 1; q0 -= 24; 605 // while(iq[jz]==0) { jz--; q0-=24;} 606 // } else { /* break z into 24-bit if necessary */ 607 // z = scalbnA(z,-q0); 608 // if(z>=two24B) { 609 // fw = (double)((int)(twon24*z)); 610 // iq[jz] = (int)(z-two24B*fw); 611 // jz += 1; q0 += 24; 612 // iq[jz] = (int) fw; 613 // } else iq[jz] = (int) z ; 614 // } 615 // 616 // /* convert integer "bit" chunk to floating-point value */ 617 // fw = scalbnA(one,q0); 618 // for(i=jz;i>=0;i--) { 619 // q[i] = fw*(double)iq[i]; fw*=twon24; 620 // } 621 // 622 // /* compute PIo2[0,...,jp]*q[jz,...,0] */ 623 // for(i=jz;i>=0;i--) { 624 // for(fw=0.0,k=0;k<=jp&&k<=jz-i;k++) fw += PIo2[k]*q[i+k]; 625 // fq[jz-i] = fw; 626 // } 627 // 628 // // NOTE: switch below is eliminated, because prec is always 2 for doubles 629 // /* compress fq[] into y[] */ 630 // //switch(prec) { 631 // //case 0: 632 // // fw = 0.0; 633 // // for (i=jz;i>=0;i--) fw += fq[i]; 634 // // y[0] = (ih==0)? fw: -fw; 635 // // break; 636 // //case 1: 637 // //case 2: 638 // fw = 0.0; 639 // for (i=jz;i>=0;i--) fw += fq[i]; 640 // y[0] = (ih==0)? fw: -fw; 641 // fw = fq[0]-fw; 642 // for (i=1;i<=jz;i++) fw += fq[i]; 643 // y[1] = (ih==0)? fw: -fw; 644 // // break; 645 // //case 3: /* painful */ 646 // // for (i=jz;i>0;i--) { 647 // // fw = fq[i-1]+fq[i]; 648 // // fq[i] += fq[i-1]-fw; 649 // // fq[i-1] = fw; 650 // // } 651 // // for (i=jz;i>1;i--) { 652 // // fw = fq[i-1]+fq[i]; 653 // // fq[i] += fq[i-1]-fw; 654 // // fq[i-1] = fw; 655 // // } 656 // // for (fw=0.0,i=jz;i>=2;i--) fw += fq[i]; 657 // // if(ih==0) { 658 // // y[0] = fq[0]; y[1] = fq[1]; y[2] = fw; 659 // // } else { 660 // // y[0] = -fq[0]; y[1] = -fq[1]; y[2] = -fw; 661 // // } 662 // //} 663 // return n&7; 664 //} 665 // 666 // END __kernel_rem_pio2 PSEUDO CODE 667 // 668 // Changes between fdlibm and intrinsic: 669 // 1. One loop is unrolled and vectorized (see comments in code) 670 // 2. One loop is split into 2 loops (see comments in code) 671 // 3. Non-double code is removed(last switch). Sevaral variables became 672 // constants because of that (see comments in code) 673 // 4. Use of jx, which is nx-1 instead of nx 674 // Assumptions: 675 // 1. Assume |X| >= PI/4 676 // Input and output: 677 // 1. Input: X = r0, jx == nx - 1 == r6, e0 == rscratch1 678 // 2. Return n in r2, y[0] == y0 == v4, y[1] == y1 == v5 679 // NOTE: general purpose register names match local variable names in C code 680 // NOTE: fpu registers are actively reused. See comments in code about their usage 681 void MacroAssembler::generate__kernel_rem_pio2(address two_over_pi, address pio2) { 682 Label Q_DONE, JX_IS_0, JX_IS_2, COMP_INNER_LOOP, RECOMP_FOR2, Q0_ZERO_CMP_LT, 683 RECOMP_CHECK_DONE_NOT_ZERO, Q0_ZERO_CMP_DONE, COMP_FOR, Q0_ZERO_CMP_EQ, 684 INIT_F_ZERO, RECOMPUTE, IH_FOR_INCREMENT, IH_FOR_STORE, RECOMP_CHECK_DONE, 685 Z_IS_LESS_THAN_TWO24B, Z_IS_ZERO, FW_Y1_NO_NEGATION, 686 RECOMP_FW_UPDATED, Z_ZERO_CHECK_DONE, FW_FOR1, IH_AFTER_SWITCH, IH_HANDLED, 687 CONVERTION_FOR, FW_Y0_NO_NEGATION, FW_FOR1_DONE, FW_FOR2, FW_FOR2_DONE, 688 IH_FOR, SKIP_F_LOAD, RECOMP_FOR1, RECOMP_FIRST_FOR, INIT_F_COPY, 689 RECOMP_FOR1_CHECK; 690 Register tmp2 = r1, n = r2, jv = r4, tmp5 = r5, jx = r6, 691 tmp3 = r7, iqBase = r10, ih = r11, tmp4 = r12, tmp1 = r13, 692 jz = r14, j = r15, twoOverPiBase = r16, i = r17, qBase = r18; 693 // jp = jk == init_jk[prec] = init_jk[2] == {2,3,4,6}[2] == 4 694 // jx = nx - 1 695 lea(twoOverPiBase, ExternalAddress(two_over_pi)); 696 cmpw(jv, zr); 697 addw(tmp4, jx, 4); // tmp4 = m = jx + jk = jx + 4. jx is in {0,1,2} so m is in [4,5,6] 698 cselw(jv, jv, zr, GE); 699 fmovd(v26, 0.0); 700 addw(tmp5, jv, 1); // jv+1 701 subsw(j, jv, jx); 702 add(qBase, sp, 320); // base of q[] 703 msubw(rscratch1, i, tmp5, rscratch1); // q0 = e0-24*(jv+1) 704 // use double f[20], fq[20], q[20], iq[20] on stack, which is 705 // (20 + 20 + 20) x 8 + 20 x 4 = 560 bytes. From lower to upper addresses it 706 // will contain f[20], fq[20], q[20], iq[20] 707 // now initialize f[20] indexes 0..m (inclusive) 708 // for(i=0;i<=m;i++,j++) f[i] = (j<0)? zeroB : /* NOTE: converted to double */ ipio2[j]; // (double) ipio2[j]; 709 mov(tmp5, sp); 710 711 block_comment("for(i=0;i<=m;i++,j++) f[i] = (j<0)? zeroB : /* NOTE: converted to double */ ipio2[j]; // (double) ipio2[j];"); { 712 eorw(i, i, i); 713 br(GE, INIT_F_COPY); 714 bind(INIT_F_ZERO); 715 stpq(v26, v26, Address(post(tmp5, 32))); 716 addw(i, i, 4); 717 addsw(j, j, 4); 718 br(LT, INIT_F_ZERO); 719 subw(i, i, j); 720 movw(j, zr); 721 bind(INIT_F_COPY); 722 add(tmp1, twoOverPiBase, j, LSL, 3); // ipio2[j] start address 723 ld1(v18, v19, v20, v21, T16B, tmp1); 724 add(tmp5, sp, i, ext::uxtx, 3); 725 st1(v18, v19, v20, v21, T16B, tmp5); 726 } 727 // v18..v21 can actually contain f[0..7] 728 cbz(i, SKIP_F_LOAD); // i == 0 => f[i] == f[0] => already loaded 729 ld1(v18, v19, v20, v21, T2D, Address(sp)); // load f[0..7] 730 bind(SKIP_F_LOAD); 731 // calculate 2^q0 and 2^-q0, which we'll need further. 732 // q0 is exponent. So, calculate biased exponent(q0+1023) 733 negw(tmp4, rscratch1); 734 addw(tmp5, rscratch1, 1023); 735 addw(tmp4, tmp4, 1023); 736 // Unroll following for(s) depending on jx in [0,1,2] 737 // for (i=0;i<=jk;i++) { 738 // for(j=0,fw=0.0;j<=jx;j++) fw += x[j]*f[jx+i-j]; q[i] = fw; 739 // } 740 // Unrolling for jx == 0 case: 741 // q[0] = x[0] * f[0] 742 // q[1] = x[0] * f[1] 743 // q[2] = x[0] * f[2] 744 // q[3] = x[0] * f[3] 745 // q[4] = x[0] * f[4] 746 // 747 // Vectorization for unrolled jx == 0 case: 748 // {q[0], q[1]} = {f[0], f[1]} * x[0] 749 // {q[2], q[3]} = {f[2], f[3]} * x[0] 750 // q[4] = f[4] * x[0] 751 // 752 // Unrolling for jx == 1 case: 753 // q[0] = x[0] * f[1] + x[1] * f[0] 754 // q[1] = x[0] * f[2] + x[1] * f[1] 755 // q[2] = x[0] * f[3] + x[1] * f[2] 756 // q[3] = x[0] * f[4] + x[1] * f[3] 757 // q[4] = x[0] * f[5] + x[1] * f[4] 758 // 759 // Vectorization for unrolled jx == 1 case: 760 // {q[0], q[1]} = {f[0], f[1]} * x[1] 761 // {q[2], q[3]} = {f[2], f[3]} * x[1] 762 // q[4] = f[4] * x[1] 763 // {q[0], q[1]} += {f[1], f[2]} * x[0] 764 // {q[2], q[3]} += {f[3], f[4]} * x[0] 765 // q[4] += f[5] * x[0] 766 // 767 // Unrolling for jx == 2 case: 768 // q[0] = x[0] * f[2] + x[1] * f[1] + x[2] * f[0] 769 // q[1] = x[0] * f[3] + x[1] * f[2] + x[2] * f[1] 770 // q[2] = x[0] * f[4] + x[1] * f[3] + x[2] * f[2] 771 // q[3] = x[0] * f[5] + x[1] * f[4] + x[2] * f[3] 772 // q[4] = x[0] * f[6] + x[1] * f[5] + x[2] * f[4] 773 // 774 // Vectorization for unrolled jx == 2 case: 775 // {q[0], q[1]} = {f[0], f[1]} * x[2] 776 // {q[2], q[3]} = {f[2], f[3]} * x[2] 777 // q[4] = f[4] * x[2] 778 // {q[0], q[1]} += {f[1], f[2]} * x[1] 779 // {q[2], q[3]} += {f[3], f[4]} * x[1] 780 // q[4] += f[5] * x[1] 781 // {q[0], q[1]} += {f[2], f[3]} * x[0] 782 // {q[2], q[3]} += {f[4], f[5]} * x[0] 783 // q[4] += f[6] * x[0] 784 block_comment("unrolled and vectorized computation of q[0]..q[jk]"); { 785 cmpw(jx, 1); 786 lsl(tmp5, tmp5, 52); // now it's 2^q0 double value 787 lsl(tmp4, tmp4, 52); // now it's 2^-q0 double value 788 br(LT, JX_IS_0); 789 add(i, sp, 8); 790 ldpq(v26, v27, i); // load f[1..4] 791 br(GT, JX_IS_2); 792 // jx == 1 793 fmulxvs(v28, T2D, v18, v7); // f[0,1] * x[1] 794 fmulxvs(v29, T2D, v19, v7); // f[2,3] * x[1] 795 fmuld(v30, v20, v7); // f[4] * x[1] 796 fmlavs(v28, T2D, v26, v6, 0); 797 fmlavs(v29, T2D, v27, v6, 0); 798 fmlavs(v30, T2D, v6, v20, 1); // v30 += f[5] * x[0] 799 b(Q_DONE); 800 bind(JX_IS_2); 801 fmulxvs(v28, T2D, v18, v3); // f[0,1] * x[2] 802 fmulxvs(v29, T2D, v19, v3); // f[2,3] * x[2] 803 fmuld(v30, v20, v3); // f[4] * x[2] 804 fmlavs(v28, T2D, v26, v7, 0); 805 fmlavs(v29, T2D, v27, v7, 0); 806 fmlavs(v30, T2D, v7, v20, 1); // v30 += f[5] * x[1] 807 fmlavs(v28, T2D, v19, v6, 0); 808 fmlavs(v29, T2D, v20, v6, 0); 809 fmlavs(v30, T2D, v6, v21, 0); // v30 += f[6] * x[0] 810 b(Q_DONE); 811 bind(JX_IS_0); 812 fmulxvs(v28, T2D, v18, v6); // f[0,1] * x[0] 813 fmulxvs(v29, T2D, v19, v6); // f[2,3] * x[0] 814 fmuld(v30, v20, v6); // f[4] * x[0] 815 bind(Q_DONE); 816 st1(v28, v29, v30, T2D, Address(qBase)); // save calculated q[0]...q[jk] 817 } 818 movz(i, 0x3E70, 48); 819 movw(jz, 4); 820 fmovd(v17, i); // v17 = twon24 821 fmovd(v30, tmp5); // 2^q0 822 fmovd(v21, 0.125); 823 fmovd(v20, 8.0); 824 fmovd(v22, tmp4); // 2^-q0 825 826 block_comment("recompute loop"); { 827 bind(RECOMPUTE); 828 // for(i=0,j=jz,z=q[jz];j>0;i++,j--) { 829 // fw = (double)((int)(twon24* z)); 830 // iq[i] = (int)(z-two24A*fw); 831 // z = q[j-1]+fw; 832 // } 833 block_comment("distill q[] into iq[] reversingly"); { 834 eorw(i, i, i); 835 movw(j, jz); 836 add(tmp2, qBase, jz, LSL, 3); // q[jz] address 837 ldrd(v18, post(tmp2, -8)); // z = q[j] and moving address to q[j-1] 838 bind(RECOMP_FIRST_FOR); 839 ldrd(v27, post(tmp2, -8)); 840 fmuld(v29, v17, v18); // twon24*z 841 frintzd(v29, v29); // (double)(int) 842 fmsubd(v28, v10, v29, v18); // v28 = z-two24A*fw 843 fcvtzdw(tmp1, v28); // (int)(z-two24A*fw) 844 strw(tmp1, Address(iqBase, i, Address::lsl(2))); 845 faddd(v18, v27, v29); 846 add(i, i, 1); 847 subs(j, j, 1); 848 br(GT, RECOMP_FIRST_FOR); 849 } 850 // compute n 851 fmuld(v18, v18, v30); 852 fmuld(v2, v18, v21); 853 frintmd(v2, v2); // v2 = floor(v2) == rounding towards -inf 854 fmsubd(v18, v2, v20, v18); // z -= 8.0*floor(z*0.125); 855 movw(ih, 2); 856 frintzd(v2, v18); // v2 = (double)((int)z) 857 fcvtzdw(n, v18); // n = (int) z; 858 fsubd(v18, v18, v2); // z -= (double)n; 859 860 block_comment("q0-dependent initialization"); { 861 cmpw(rscratch1, 0); // if (q0 > 0) 862 br(LT, Q0_ZERO_CMP_LT); 863 subw(j, jz, 1); // j = jz - 1 864 ldrw(tmp2, Address(iqBase, j, Address::lsl(2))); // tmp2 = iq[jz-1] 865 br(EQ, Q0_ZERO_CMP_EQ); 866 movw(tmp4, 24); 867 subw(tmp4, tmp4, rscratch1); // == 24 - q0 868 lsrvw(i, tmp2, tmp4); // i = iq[jz-1] >> (24-q0) 869 lslvw(tmp5, i, tmp4); 870 subw(tmp2, tmp2, tmp5); // iq[jz-1] -= i<<(24-q0); 871 strw(tmp2, Address(iqBase, j, Address::lsl(2))); // store iq[jz-1] 872 subw(rscratch2, tmp4, 1); // == 23 - q0 873 addw(n, n, i); // n+=i 874 lsrvw(ih, tmp2, rscratch2); // ih = iq[jz-1] >> (23-q0) 875 b(Q0_ZERO_CMP_DONE); 876 bind(Q0_ZERO_CMP_EQ); 877 lsr(ih, tmp2, 23); // ih = iq[z-1] >> 23 878 b(Q0_ZERO_CMP_DONE); 879 bind(Q0_ZERO_CMP_LT); 880 fmovd(v4, 0.5); 881 fcmpd(v18, v4); 882 cselw(ih, zr, ih, LT); // if (z<0.5) ih = 0 883 } 884 bind(Q0_ZERO_CMP_DONE); 885 cmpw(ih, zr); 886 br(LE, IH_HANDLED); 887 888 block_comment("if(ih>) {"); { 889 // use rscratch2 as carry 890 891 block_comment("for(i=0;i<jz ;i++) {...}"); { 892 addw(n, n, 1); 893 eorw(i, i, i); 894 eorw(rscratch2, rscratch2, rscratch2); 895 bind(IH_FOR); 896 ldrw(j, Address(iqBase, i, Address::lsl(2))); // j = iq[i] 897 movw(tmp3, 0x1000000); 898 subw(tmp3, tmp3, rscratch2); 899 cbnzw(rscratch2, IH_FOR_STORE); 900 cbzw(j, IH_FOR_INCREMENT); 901 movw(rscratch2, 1); 902 bind(IH_FOR_STORE); 903 subw(tmp3, tmp3, j); 904 strw(tmp3, Address(iqBase, i, Address::lsl(2))); // iq[i] = 0xffffff - j 905 bind(IH_FOR_INCREMENT); 906 addw(i, i, 1); 907 cmpw(i, jz); 908 br(LT, IH_FOR); 909 } 910 911 block_comment("if(q0>0) {"); { 912 cmpw(rscratch1, zr); 913 br(LE, IH_AFTER_SWITCH); 914 // tmp3 still has iq[jz-1] value. no need to reload 915 // now, zero high tmp3 bits (rscratch1 number of bits) 916 movw(j, -1); 917 subw(i, jz, 1); // set i to jz-1 918 lsrv(j, j, rscratch1); 919 andw(tmp3, tmp3, j, LSR, 8); // we have 24-bit-based constants 920 strw(tmp3, Address(iqBase, i, Address::lsl(2))); // save iq[jz-1] 921 } 922 bind(IH_AFTER_SWITCH); 923 cmpw(ih, 2); 924 br(NE, IH_HANDLED); 925 926 block_comment("if(ih==2) {"); { 927 fmovd(v25, 1.0); 928 fsubd(v18, v25, v18); // z = one - z; 929 cbzw(rscratch2, IH_HANDLED); 930 fsubd(v18, v18, v30); // z -= scalbnA(one,q0); 931 } 932 } 933 bind(IH_HANDLED); 934 // check if recomputation is needed 935 fcmpd(v18, 0.0); 936 br(NE, RECOMP_CHECK_DONE_NOT_ZERO); 937 938 block_comment("if(z==zeroB) {"); { 939 940 block_comment("for (i=jz-1;i>=jk;i--) j |= iq[i];"); { 941 subw(i, jz, 1); 942 eorw(j, j, j); 943 b(RECOMP_FOR1_CHECK); 944 bind(RECOMP_FOR1); 945 ldrw(tmp1, Address(iqBase, i, Address::lsl(2))); 946 orrw(j, j, tmp1); 947 subw(i, i, 1); 948 bind(RECOMP_FOR1_CHECK); 949 cmpw(i, 4); 950 br(GE, RECOMP_FOR1); 951 } 952 cbnzw(j, RECOMP_CHECK_DONE); 953 954 block_comment("if(j==0) {"); { 955 // for(k=1;iq[jk-k]==0;k++); // let's unroll it. jk == 4. So, read 956 // iq[3], iq[2], iq[1], iq[0] until non-zero value 957 ldp(tmp1, tmp3, iqBase); // iq[0..3] 958 movw(j, 2); 959 cmp(tmp3, zr); 960 csel(tmp1, tmp1, tmp3, EQ); // set register for further consideration 961 cselw(j, j, zr, EQ); // set initial k. Use j as k 962 cmp(zr, tmp1, LSR, 32); 963 addw(i, jz, 1); 964 csincw(j, j, j, NE); 965 966 block_comment("for(i=jz+1;i<=jz+k;i++) {...}"); { 967 addw(jz, i, j); // i = jz+1, j = k-1. j+i = jz+k (which is a new jz) 968 bind(RECOMP_FOR2); 969 addw(tmp1, jv, i); 970 ldrd(v29, Address(twoOverPiBase, tmp1, Address::lsl(3))); 971 addw(tmp2, jx, i); 972 strd(v29, Address(sp, tmp2, Address::lsl(3))); 973 // f[jx+i] = /* NOTE: converted to double */ ipio2[jv+i]; //(double) ipio2[jv+i]; 974 // since jx = 0, 1 or 2 we can unroll it: 975 // for(j=0,fw=0.0;j<=jx;j++) fw += x[j]*f[jx+i-j]; 976 // f[jx+i-j] == (for first iteration) f[jx+i], which is already v29 977 add(tmp2, sp, tmp2, ext::uxtx, 3); // address of f[jx+i] 978 ldpd(v4, v5, Address(tmp2, -16)); // load f[jx+i-2] and f[jx+i-1] 979 fmuld(v26, v6, v29); // initial fw 980 cbzw(jx, RECOMP_FW_UPDATED); 981 fmaddd(v26, v7, v5, v26); 982 cmpw(jx, 1); 983 br(EQ, RECOMP_FW_UPDATED); 984 fmaddd(v26, v3, v4, v26); 985 bind(RECOMP_FW_UPDATED); 986 strd(v26, Address(qBase, i, Address::lsl(3))); // q[i] = fw; 987 addw(i, i, 1); 988 cmpw(i, jz); // jz here is "old jz" + k 989 br(LE, RECOMP_FOR2); 990 } 991 b(RECOMPUTE); 992 } 993 } 994 } 995 bind(RECOMP_CHECK_DONE); 996 // chop off zero terms 997 fcmpd(v18, 0.0); 998 br(EQ, Z_IS_ZERO); 999 1000 block_comment("else block of if(z==0.0) {"); { 1001 bind(RECOMP_CHECK_DONE_NOT_ZERO); 1002 fmuld(v18, v18, v22); 1003 fcmpd(v18, v10); // v10 is stil two24A 1004 br(LT, Z_IS_LESS_THAN_TWO24B); 1005 fmuld(v1, v18, v17); // twon24*z 1006 frintzd(v1, v1); // v1 = (double)(int)(v1) 1007 fmsubd(v2, v10, v1, v18); 1008 fcvtzdw(tmp3, v1); // (int)fw 1009 fcvtzdw(tmp2, v2); // double to int 1010 strw(tmp2, Address(iqBase, jz, Address::lsl(2))); 1011 addw(rscratch1, rscratch1, 24); 1012 addw(jz, jz, 1); 1013 strw(tmp3, Address(iqBase, jz, Address::lsl(2))); // iq[jz] = (int) fw 1014 b(Z_ZERO_CHECK_DONE); 1015 bind(Z_IS_LESS_THAN_TWO24B); 1016 fcvtzdw(tmp3, v18); // (int)z 1017 strw(tmp3, Address(iqBase, jz, Address::lsl(2))); // iq[jz] = (int) z 1018 b(Z_ZERO_CHECK_DONE); 1019 } 1020 1021 block_comment("if(z==0.0) {"); { 1022 bind(Z_IS_ZERO); 1023 subw(jz, jz, 1); 1024 ldrw(tmp1, Address(iqBase, jz, Address::lsl(2))); 1025 subw(rscratch1, rscratch1, 24); 1026 cbz(tmp1, Z_IS_ZERO); 1027 } 1028 bind(Z_ZERO_CHECK_DONE); 1029 // convert integer "bit" chunk to floating-point value 1030 // v17 = twon24 1031 // update v30, which was scalbnA(1.0, <old q0>); 1032 addw(tmp2, rscratch1, 1023); // biased exponent 1033 lsl(tmp2, tmp2, 52); // put at correct position 1034 mov(i, jz); 1035 fmovd(v30, tmp2); 1036 1037 block_comment("for(i=jz;i>=0;i--) {q[i] = fw*(double)iq[i]; fw*=twon24;}"); { 1038 bind(CONVERTION_FOR); 1039 ldrw(tmp1, Address(iqBase, i, Address::lsl(2))); 1040 scvtfwd(v31, tmp1); 1041 fmuld(v31, v31, v30); 1042 strd(v31, Address(qBase, i, Address::lsl(3))); 1043 fmuld(v30, v30, v17); 1044 subsw(i, i, 1); 1045 br(GE, CONVERTION_FOR); 1046 } 1047 add(rscratch2, sp, 160); // base for fq 1048 // reusing twoOverPiBase 1049 lea(twoOverPiBase, ExternalAddress(pio2)); 1050 1051 block_comment("compute PIo2[0,...,jp]*q[jz,...,0]. for(i=jz;i>=0;i--) {...}"); { 1052 movw(i, jz); 1053 movw(tmp2, zr); // tmp2 will keep jz - i == 0 at start 1054 bind(COMP_FOR); 1055 // for(fw=0.0,k=0;k<=jp&&k<=jz-i;k++) fw += PIo2[k]*q[i+k]; 1056 fmovd(v30, 0.0); 1057 add(tmp5, qBase, i, LSL, 3); // address of q[i+k] for k==0 1058 movw(tmp3, 4); 1059 movw(tmp4, zr); // used as k 1060 cmpw(tmp2, 4); 1061 add(tmp1, qBase, i, LSL, 3); // used as q[i] address 1062 cselw(tmp3, tmp2, tmp3, LE); // min(jz - i, jp) 1063 1064 block_comment("for(fw=0.0,k=0;k<=jp&&k<=jz-i;k++) fw += PIo2[k]*q[i+k];"); { 1065 bind(COMP_INNER_LOOP); 1066 ldrd(v18, Address(tmp1, tmp4, Address::lsl(3))); // q[i+k] 1067 ldrd(v19, Address(twoOverPiBase, tmp4, Address::lsl(3))); // PIo2[k] 1068 fmaddd(v30, v18, v19, v30); // fw += PIo2[k]*q[i+k]; 1069 addw(tmp4, tmp4, 1); // k++ 1070 cmpw(tmp4, tmp3); 1071 br(LE, COMP_INNER_LOOP); 1072 } 1073 strd(v30, Address(rscratch2, tmp2, Address::lsl(3))); // fq[jz-i] 1074 add(tmp2, tmp2, 1); 1075 subsw(i, i, 1); 1076 br(GE, COMP_FOR); 1077 } 1078 1079 block_comment("switch(prec) {...}. case 2:"); { 1080 // compress fq into y[] 1081 // remember prec == 2 1082 1083 block_comment("for (i=jz;i>=0;i--) fw += fq[i];"); { 1084 fmovd(v4, 0.0); 1085 mov(i, jz); 1086 bind(FW_FOR1); 1087 ldrd(v1, Address(rscratch2, i, Address::lsl(3))); 1088 subsw(i, i, 1); 1089 faddd(v4, v4, v1); 1090 br(GE, FW_FOR1); 1091 } 1092 bind(FW_FOR1_DONE); 1093 // v1 contains fq[0]. so, keep it so far 1094 fsubd(v5, v1, v4); // fw = fq[0] - fw 1095 cbzw(ih, FW_Y0_NO_NEGATION); 1096 fnegd(v4, v4); 1097 bind(FW_Y0_NO_NEGATION); 1098 1099 block_comment("for (i=1;i<=jz;i++) fw += fq[i];"); { 1100 movw(i, 1); 1101 cmpw(jz, 1); 1102 br(LT, FW_FOR2_DONE); 1103 bind(FW_FOR2); 1104 ldrd(v1, Address(rscratch2, i, Address::lsl(3))); 1105 addw(i, i, 1); 1106 cmp(i, jz); 1107 faddd(v5, v5, v1); 1108 br(LE, FW_FOR2); 1109 } 1110 bind(FW_FOR2_DONE); 1111 cbz(ih, FW_Y1_NO_NEGATION); 1112 fnegd(v5, v5); 1113 bind(FW_Y1_NO_NEGATION); 1114 add(sp, sp, 560); 1115 } 1116 } 1117 1118 ///* __kernel_sin( x, y, iy) 1119 // * kernel sin function on [-pi/4, pi/4], pi/4 ~ 0.7854 1120 // * Input x is assumed to be bounded by ~pi/4 in magnitude. 1121 // * Input y is the tail of x. 1122 // * Input iy indicates whether y is 0. (if iy=0, y assume to be 0). 1123 // * 1124 // * Algorithm 1125 // * 1. Since sin(-x) = -sin(x), we need only to consider positive x. 1126 // * 2. if x < 2^-27 (hx<0x3e400000 0), return x with inexact if x!=0. 1127 // * 3. sin(x) is approximated by a polynomial of degree 13 on 1128 // * [0,pi/4] 1129 // * 3 13 1130 // * sin(x) ~ x + S1*x + ... + S6*x 1131 // * where 1132 // * 1133 // * |sin(x) 2 4 6 8 10 12 | -58 1134 // * |----- - (1+S1*x +S2*x +S3*x +S4*x +S5*x +S6*x )| <= 2 1135 // * | x | 1136 // * 1137 // * 4. sin(x+y) = sin(x) + sin'(x')*y 1138 // * ~ sin(x) + (1-x*x/2)*y 1139 // * For better accuracy, let 1140 // * 3 2 2 2 2 1141 // * r = x *(S2+x *(S3+x *(S4+x *(S5+x *S6)))) 1142 // * then 3 2 1143 // * sin(x) = x + (S1*x + (x *(r-y/2)+y)) 1144 // */ 1145 //static const double 1146 //S1 = -1.66666666666666324348e-01, /* 0xBFC55555, 0x55555549 */ 1147 //S2 = 8.33333333332248946124e-03, /* 0x3F811111, 0x1110F8A6 */ 1148 //S3 = -1.98412698298579493134e-04, /* 0xBF2A01A0, 0x19C161D5 */ 1149 //S4 = 2.75573137070700676789e-06, /* 0x3EC71DE3, 0x57B1FE7D */ 1150 //S5 = -2.50507602534068634195e-08, /* 0xBE5AE5E6, 0x8A2B9CEB */ 1151 //S6 = 1.58969099521155010221e-10; /* 0x3DE5D93A, 0x5ACFD57C */ 1152 // 1153 // NOTE: S1..S6 were moved into a table: StubRoutines::aarch64::_dsin_coef 1154 // 1155 // BEGIN __kernel_sin PSEUDO CODE 1156 // 1157 //static double __kernel_sin(double x, double y, bool iy) 1158 //{ 1159 // double z,r,v; 1160 // 1161 // // NOTE: not needed. moved to dsin/dcos 1162 // //int ix; 1163 // //ix = high(x)&0x7fffffff; /* high word of x */ 1164 // 1165 // // NOTE: moved to dsin/dcos 1166 // //if(ix<0x3e400000) /* |x| < 2**-27 */ 1167 // // {if((int)x==0) return x;} /* generate inexact */ 1168 // 1169 // z = x*x; 1170 // v = z*x; 1171 // r = S2+z*(S3+z*(S4+z*(S5+z*S6))); 1172 // if(iy==0) return x+v*(S1+z*r); 1173 // else return x-((z*(half*y-v*r)-y)-v*S1); 1174 //} 1175 // 1176 // END __kernel_sin PSEUDO CODE 1177 // 1178 // Changes between fdlibm and intrinsic: 1179 // 1. Removed |x| < 2**-27 check, because if was done earlier in dsin/dcos 1180 // 2. Constants are now loaded from table dsin_coef 1181 // 3. C code parameter "int iy" was modified to "bool iyIsOne", because 1182 // iy is always 0 or 1. Also, iyIsOne branch was moved into 1183 // generation phase instead of taking it during code execution 1184 // Input ans output: 1185 // 1. Input for generated function: X argument = x 1186 // 2. Input for generator: x = register to read argument from, iyIsOne 1187 // = flag to use low argument low part or not, dsin_coef = coefficients 1188 // table address 1189 // 3. Return sin(x) value in v0 1190 void MacroAssembler::generate_kernel_sin(FloatRegister x, bool iyIsOne, 1191 address dsin_coef) { 1192 FloatRegister y = v5, z = v6, v = v7, r = v16, S1 = v17, S2 = v18, 1193 S3 = v19, S4 = v20, S5 = v21, S6 = v22, half = v23; 1194 lea(rscratch2, ExternalAddress(dsin_coef)); 1195 ldpd(S5, S6, Address(rscratch2, 32)); 1196 fmuld(z, x, x); // z = x*x; 1197 ld1(S1, S2, S3, S4, T1D, Address(rscratch2)); 1198 fmuld(v, z, x); // v = z*x; 1199 1200 block_comment("calculate r = S2+z*(S3+z*(S4+z*(S5+z*S6)))"); { 1201 fmaddd(r, z, S6, S5); 1202 // initialize "half" in current block to utilize 2nd FPU. However, it's 1203 // not a part of this block 1204 fmovd(half, 0.5); 1205 fmaddd(r, z, r, S4); 1206 fmaddd(r, z, r, S3); 1207 fmaddd(r, z, r, S2); 1208 } 1209 1210 if (!iyIsOne) { 1211 // return x+v*(S1+z*r); 1212 fmaddd(S1, z, r, S1); 1213 fmaddd(v0, v, S1, x); 1214 } else { 1215 // return x-((z*(half*y-v*r)-y)-v*S1); 1216 fmuld(S6, half, y); // half*y 1217 fmsubd(S6, v, r, S6); // half*y-v*r 1218 fmsubd(S6, z, S6, y); // y - z*(half*y-v*r) = - (z*(half*y-v*r)-y) 1219 fmaddd(S6, v, S1, S6); // - (z*(half*y-v*r)-y) + v*S1 == -((z*(half*y-v*r)-y)-v*S1) 1220 faddd(v0, x, S6); 1221 } 1222 } 1223 1224 ///* 1225 // * __kernel_cos( x, y ) 1226 // * kernel cos function on [-pi/4, pi/4], pi/4 ~ 0.785398164 1227 // * Input x is assumed to be bounded by ~pi/4 in magnitude. 1228 // * Input y is the tail of x. 1229 // * 1230 // * Algorithm 1231 // * 1. Since cos(-x) = cos(x), we need only to consider positive x. 1232 // * 2. if x < 2^-27 (hx<0x3e400000 0), return 1 with inexact if x!=0. 1233 // * 3. cos(x) is approximated by a polynomial of degree 14 on 1234 // * [0,pi/4] 1235 // * 4 14 1236 // * cos(x) ~ 1 - x*x/2 + C1*x + ... + C6*x 1237 // * where the remez error is 1238 // * 1239 // * | 2 4 6 8 10 12 14 | -58 1240 // * |cos(x)-(1-.5*x +C1*x +C2*x +C3*x +C4*x +C5*x +C6*x )| <= 2 1241 // * | | 1242 // * 1243 // * 4 6 8 10 12 14 1244 // * 4. let r = C1*x +C2*x +C3*x +C4*x +C5*x +C6*x , then 1245 // * cos(x) = 1 - x*x/2 + r 1246 // * since cos(x+y) ~ cos(x) - sin(x)*y 1247 // * ~ cos(x) - x*y, 1248 // * a correction term is necessary in cos(x) and hence 1249 // * cos(x+y) = 1 - (x*x/2 - (r - x*y)) 1250 // * For better accuracy when x > 0.3, let qx = |x|/4 with 1251 // * the last 32 bits mask off, and if x > 0.78125, let qx = 0.28125. 1252 // * Then 1253 // * cos(x+y) = (1-qx) - ((x*x/2-qx) - (r-x*y)). 1254 // * Note that 1-qx and (x*x/2-qx) is EXACT here, and the 1255 // * magnitude of the latter is at least a quarter of x*x/2, 1256 // * thus, reducing the rounding error in the subtraction. 1257 // */ 1258 // 1259 //static const double 1260 //C1 = 4.16666666666666019037e-02, /* 0x3FA55555, 0x5555554C */ 1261 //C2 = -1.38888888888741095749e-03, /* 0xBF56C16C, 0x16C15177 */ 1262 //C3 = 2.48015872894767294178e-05, /* 0x3EFA01A0, 0x19CB1590 */ 1263 //C4 = -2.75573143513906633035e-07, /* 0xBE927E4F, 0x809C52AD */ 1264 //C5 = 2.08757232129817482790e-09, /* 0x3E21EE9E, 0xBDB4B1C4 */ 1265 //C6 = -1.13596475577881948265e-11; /* 0xBDA8FAE9, 0xBE8838D4 */ 1266 // 1267 // NOTE: C1..C6 were moved into a table: StubRoutines::aarch64::_dcos_coef 1268 // 1269 // BEGIN __kernel_cos PSEUDO CODE 1270 // 1271 //static double __kernel_cos(double x, double y) 1272 //{ 1273 // double a,h,z,r,qx=0; 1274 // 1275 // // NOTE: ix is already initialized in dsin/dcos. Reuse value from register 1276 // //int ix; 1277 // //ix = high(x)&0x7fffffff; /* ix = |x|'s high word*/ 1278 // 1279 // // NOTE: moved to dsin/dcos 1280 // //if(ix<0x3e400000) { /* if x < 2**27 */ 1281 // // if(((int)x)==0) return one; /* generate inexact */ 1282 // //} 1283 // 1284 // z = x*x; 1285 // r = z*(C1+z*(C2+z*(C3+z*(C4+z*(C5+z*C6))))); 1286 // if(ix < 0x3FD33333) /* if |x| < 0.3 */ 1287 // return one - (0.5*z - (z*r - x*y)); 1288 // else { 1289 // if(ix > 0x3fe90000) { /* x > 0.78125 */ 1290 // qx = 0.28125; 1291 // } else { 1292 // set_high(&qx, ix-0x00200000); /* x/4 */ 1293 // set_low(&qx, 0); 1294 // } 1295 // h = 0.5*z-qx; 1296 // a = one-qx; 1297 // return a - (h - (z*r-x*y)); 1298 // } 1299 //} 1300 // 1301 // END __kernel_cos PSEUDO CODE 1302 // 1303 // Changes between fdlibm and intrinsic: 1304 // 1. Removed |x| < 2**-27 check, because if was done earlier in dsin/dcos 1305 // 2. Constants are now loaded from table dcos_coef 1306 // Input and output: 1307 // 1. Input for generated function: X argument = x 1308 // 2. Input for generator: x = register to read argument from, dcos_coef 1309 // = coefficients table address 1310 // 2. Return cos(x) value in v0 1311 void MacroAssembler::generate_kernel_cos(FloatRegister x, address dcos_coef) { 1312 Register ix = r3; 1313 FloatRegister qx = v1, h = v2, a = v3, y = v5, z = v6, r = v7, C1 = v18, 1314 C2 = v19, C3 = v20, C4 = v21, C5 = v22, C6 = v23, one = v25, half = v26; 1315 Label IX_IS_LARGE, SET_QX_CONST, DONE, QX_SET; 1316 lea(rscratch2, ExternalAddress(dcos_coef)); 1317 ldpd(C5, C6, Address(rscratch2, 32)); // load C5, C6 1318 fmuld(z, x, x); // z=x^2 1319 ld1(C1, C2, C3, C4, T1D, Address(rscratch2)); // load C1..C3\4 1320 block_comment("calculate r = z*(C1+z*(C2+z*(C3+z*(C4+z*(C5+z*C6)))))"); { 1321 fmaddd(r, z, C6, C5); 1322 fmovd(half, 0.5); 1323 fmaddd(r, z, r, C4); 1324 fmuld(y, x, y); 1325 fmaddd(r, z, r, C3); 1326 mov(rscratch1, 0x3FD33333); 1327 fmaddd(r, z, r, C2); 1328 fmuld(x, z, z); // x = z^2 1329 fmaddd(r, z, r, C1); // r = C1+z(C2+z(C4+z(C5+z*C6))) 1330 } 1331 // need to multiply r by z to have "final" r value 1332 fmovd(one, 1.0); 1333 cmp(ix, rscratch1); 1334 br(GT, IX_IS_LARGE); 1335 block_comment("if(ix < 0x3FD33333) return one - (0.5*z - (z*r - x*y))"); { 1336 // return 1.0 - (0.5*z - (z*r - x*y)) = 1.0 - (0.5*z + (x*y - z*r)) 1337 fmsubd(v0, x, r, y); 1338 fmaddd(v0, half, z, v0); 1339 fsubd(v0, one, v0); 1340 b(DONE); 1341 } 1342 block_comment("if(ix >= 0x3FD33333)"); { 1343 bind(IX_IS_LARGE); 1344 movz(rscratch2, 0x3FE9, 16); 1345 cmp(ix, rscratch2); 1346 br(GT, SET_QX_CONST); 1347 block_comment("set_high(&qx, ix-0x00200000); set_low(&qx, 0);"); { 1348 subw(rscratch2, ix, 0x00200000); 1349 lsl(rscratch2, rscratch2, 32); 1350 fmovd(qx, rscratch2); 1351 } 1352 b(QX_SET); 1353 bind(SET_QX_CONST); 1354 block_comment("if(ix > 0x3fe90000) qx = 0.28125;"); { 1355 fmovd(qx, 0.28125); 1356 } 1357 bind(QX_SET); 1358 fnmsub(C6, x, r, y); // z*r - xy 1359 fnmsub(h, half, z, qx); // h = 0.5*z - qx 1360 fsubd(a, one, qx); // a = 1-qx 1361 fsubd(C6, h, C6); // = h - (z*r - x*y) 1362 fsubd(v0, a, C6); 1363 } 1364 bind(DONE); 1365 } 1366 1367 // generate_dsin_dcos creates stub for dsin and dcos 1368 // Generation is done via single call because dsin and dcos code is almost the 1369 // same(see C code below). These functions work as follows: 1370 // 1) handle corner cases: |x| ~< pi/4, x is NaN or INF, |x| < 2**-27 1371 // 2) perform argument reduction if required 1372 // 3) call kernel_sin or kernel_cos which approximate sin/cos via polynomial 1373 // 1374 // BEGIN dsin/dcos PSEUDO CODE 1375 // 1376 //dsin_dcos(jdouble x, bool isCos) { 1377 // double y[2],z=0.0; 1378 // int n, ix; 1379 // 1380 // /* High word of x. */ 1381 // ix = high(x); 1382 // 1383 // /* |x| ~< pi/4 */ 1384 // ix &= 0x7fffffff; 1385 // if(ix <= 0x3fe921fb) return isCos ? __kernel_cos : __kernel_sin(x,z,0); 1386 // 1387 // /* sin/cos(Inf or NaN) is NaN */ 1388 // else if (ix>=0x7ff00000) return x-x; 1389 // else if (ix<0x3e400000) { /* if ix < 2**27 */ 1390 // if(((int)x)==0) return isCos ? one : x; /* generate inexact */ 1391 // } 1392 // /* argument reduction needed */ 1393 // else { 1394 // n = __ieee754_rem_pio2(x,y); 1395 // switch(n&3) { 1396 // case 0: return isCos ? __kernel_cos(y[0],y[1]) : __kernel_sin(y[0],y[1], true); 1397 // case 1: return isCos ? -__kernel_sin(y[0],y[1],true) : __kernel_cos(y[0],y[1]); 1398 // case 2: return isCos ? -__kernel_cos(y[0],y[1]) : -__kernel_sin(y[0],y[1], true); 1399 // default: 1400 // return isCos ? __kernel_sin(y[0],y[1],1) : -__kernel_cos(y[0],y[1]); 1401 // } 1402 // } 1403 //} 1404 // END dsin/dcos PSEUDO CODE 1405 // 1406 // Changes between fdlibm and intrinsic: 1407 // 1. Moved ix < 2**27 from kernel_sin/kernel_cos into dsin/dcos 1408 // 2. Final switch use equivalent bit checks(tbz/tbnz) 1409 // Input ans output: 1410 // 1. Input for generated function: X = r0 1411 // 2. Input for generator: isCos = generate sin or cos, npio2_hw = address 1412 // of npio2_hw table, two_over_pi = address of two_over_pi table, 1413 // pio2 = address if pio2 table, dsin_coef = address if dsin_coef table, 1414 // dcos_coef = address of dcos_coef table 1415 // 3. Return result in v0 1416 // NOTE: general purpose register names match local variable names in C code 1417 void MacroAssembler::generate_dsin_dcos(bool isCos, address npio2_hw, 1418 address two_over_pi, address pio2, address dsin_coef, address dcos_coef) { 1419 const int POSITIVE_INFINITY_OR_NAN_PREFIX = 0x7FF0; 1420 1421 Label DONE, ARG_REDUCTION, TINY_X, RETURN_SIN, EARLY_CASE; 1422 Register X = r0, absX = r1, n = r2, ix = r3; 1423 FloatRegister y0 = v4, y1 = v5; 1424 block_comment("check |x| ~< pi/4, NaN, Inf and |x| < 2**-27 cases"); { 1425 fmovd(X, v0); 1426 mov(rscratch2, 0x3e400000); 1427 mov(rscratch1, 0x3fe921fb00000000); // pi/4. shifted to reuse later 1428 ubfm(absX, X, 0, 62); // absX 1429 movz(r10, POSITIVE_INFINITY_OR_NAN_PREFIX, 48); 1430 cmp(rscratch2, absX, LSR, 32); 1431 lsr(ix, absX, 32); // set ix 1432 br(GT, TINY_X); // handle tiny x (|x| < 2^-27) 1433 cmp(ix, rscratch1, LSR, 32); 1434 br(LE, EARLY_CASE); // if(ix <= 0x3fe921fb) return 1435 cmp(absX, r10); 1436 br(LT, ARG_REDUCTION); 1437 // X is NaN or INF(i.e. 0x7FF* or 0xFFF*). Return NaN (mantissa != 0). 1438 // Set last bit unconditionally to make it NaN 1439 orr(r10, r10, 1); 1440 fmovd(v0, r10); 1441 ret(lr); 1442 } 1443 block_comment("kernel_sin/kernel_cos: if(ix<0x3e400000) {<fast return>}"); { 1444 bind(TINY_X); 1445 if (isCos) { 1446 fmovd(v0, 1.0); 1447 } 1448 ret(lr); 1449 } 1450 bind(ARG_REDUCTION); /* argument reduction needed */ 1451 block_comment("n = __ieee754_rem_pio2(x,y);"); { 1452 generate__ieee754_rem_pio2(npio2_hw, two_over_pi, pio2); 1453 } 1454 block_comment("switch(n&3) {case ... }"); { 1455 if (isCos) { 1456 eorw(absX, n, n, LSR, 1); 1457 tbnz(n, 0, RETURN_SIN); 1458 } else { 1459 tbz(n, 0, RETURN_SIN); 1460 } 1461 generate_kernel_cos(y0, dcos_coef); 1462 if (isCos) { 1463 tbz(absX, 0, DONE); 1464 } else { 1465 tbz(n, 1, DONE); 1466 } 1467 fnegd(v0, v0); 1468 ret(lr); 1469 bind(RETURN_SIN); 1470 generate_kernel_sin(y0, true, dsin_coef); 1471 if (isCos) { 1472 tbz(absX, 0, DONE); 1473 } else { 1474 tbz(n, 1, DONE); 1475 } 1476 fnegd(v0, v0); 1477 ret(lr); 1478 } 1479 bind(EARLY_CASE); 1480 eor(y1, T8B, y1, y1); 1481 if (isCos) { 1482 generate_kernel_cos(v0, dcos_coef); 1483 } else { 1484 generate_kernel_sin(v0, false, dsin_coef); 1485 } 1486 bind(DONE); 1487 ret(lr); 1488 }