Mercurial > hg > truffle
annotate src/share/vm/runtime/sharedRuntimeTrig.cpp @ 1721:413ad0331a0c
6977924: Changes for 6975078 produce build error with certain gcc versions
Summary: The changes introduced for 6975078 assign badHeapOopVal to the _allocation field in the ResourceObj class. In 32 bit linux builds with certain versions of gcc this assignment will be flagged as an error while compiling allocation.cpp. In 32 bit builds the constant value badHeapOopVal (which is cast to an intptr_t) is negative. The _allocation field is typed as an unsigned intptr_t and gcc catches this as an error.
Reviewed-by: jcoomes, ysr, phh
| author | johnc |
|---|---|
| date | Wed, 18 Aug 2010 10:59:06 -0700 |
| parents | c18cbe5936b8 |
| children | f95d63e2154a |
| rev | line source |
|---|---|
| 0 | 1 /* |
|
1552
c18cbe5936b8
6941466: Oracle rebranding changes for Hotspot repositories
trims
parents:
1485
diff
changeset
|
2 * Copyright (c) 2001, 2005, Oracle and/or its affiliates. All rights reserved. |
| 0 | 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 * | |
|
1552
c18cbe5936b8
6941466: Oracle rebranding changes for Hotspot repositories
trims
parents:
1485
diff
changeset
|
19 * Please contact Oracle, 500 Oracle Parkway, Redwood Shores, CA 94065 USA |
|
c18cbe5936b8
6941466: Oracle rebranding changes for Hotspot repositories
trims
parents:
1485
diff
changeset
|
20 * or visit www.oracle.com if you need additional information or have any |
|
c18cbe5936b8
6941466: Oracle rebranding changes for Hotspot repositories
trims
parents:
1485
diff
changeset
|
21 * questions. |
| 0 | 22 * |
| 23 */ | |
| 24 | |
| 25 #include "incls/_precompiled.incl" | |
| 26 #include "incls/_sharedRuntimeTrig.cpp.incl" | |
| 27 | |
| 28 // This file contains copies of the fdlibm routines used by | |
| 29 // StrictMath. It turns out that it is almost always required to use | |
| 30 // these runtime routines; the Intel CPU doesn't meet the Java | |
| 31 // specification for sin/cos outside a certain limited argument range, | |
| 32 // and the SPARC CPU doesn't appear to have sin/cos instructions. It | |
| 33 // also turns out that avoiding the indirect call through function | |
| 34 // pointer out to libjava.so in SharedRuntime speeds these routines up | |
| 35 // by roughly 15% on both Win32/x86 and Solaris/SPARC. | |
| 36 | |
| 37 // Enabling optimizations in this file causes incorrect code to be | |
| 38 // generated; can not figure out how to turn down optimization for one | |
| 39 // file in the IDE on Windows | |
| 40 #ifdef WIN32 | |
| 41 # pragma optimize ( "", off ) | |
| 42 #endif | |
| 43 | |
| 1485 | 44 /* The above workaround now causes more problems with the latest MS compiler. |
| 45 * Visual Studio 2010's /GS option tries to guard against buffer overruns. | |
| 46 * /GS is on by default if you specify optimizations, which we do globally | |
| 47 * via /W3 /O2. However the above selective turning off of optimizations means | |
| 48 * that /GS issues a warning "4748". And since we treat warnings as errors (/WX) | |
| 49 * then the compilation fails. There are several possible solutions | |
| 50 * (1) Remove that pragma above as obsolete with VS2010 - requires testing. | |
| 51 * (2) Stop treating warnings as errors - would be a backward step | |
| 52 * (3) Disable /GS - may help performance but you lose the security checks | |
| 53 * (4) Disable the warning with "#pragma warning( disable : 4748 )" | |
| 54 * (5) Disable planting the code with __declspec(safebuffers) | |
| 55 * I've opted for (5) although we should investigate the local performance | |
| 56 * benefits of (1) and global performance benefit of (3). | |
| 57 */ | |
| 58 #if defined(WIN32) && (defined(_MSC_VER) && (_MSC_VER >= 1600)) | |
| 59 #define SAFEBUF __declspec(safebuffers) | |
| 60 #else | |
| 61 #define SAFEBUF | |
| 62 #endif | |
| 63 | |
| 0 | 64 #include <math.h> |
| 65 | |
| 66 // VM_LITTLE_ENDIAN is #defined appropriately in the Makefiles | |
| 67 // [jk] this is not 100% correct because the float word order may different | |
| 68 // from the byte order (e.g. on ARM) | |
| 69 #ifdef VM_LITTLE_ENDIAN | |
| 70 # define __HI(x) *(1+(int*)&x) | |
| 71 # define __LO(x) *(int*)&x | |
| 72 #else | |
| 73 # define __HI(x) *(int*)&x | |
| 74 # define __LO(x) *(1+(int*)&x) | |
| 75 #endif | |
| 76 | |
| 77 static double copysignA(double x, double y) { | |
| 78 __HI(x) = (__HI(x)&0x7fffffff)|(__HI(y)&0x80000000); | |
| 79 return x; | |
| 80 } | |
| 81 | |
| 82 /* | |
| 83 * scalbn (double x, int n) | |
| 84 * scalbn(x,n) returns x* 2**n computed by exponent | |
| 85 * manipulation rather than by actually performing an | |
| 86 * exponentiation or a multiplication. | |
| 87 */ | |
| 88 | |
| 89 static const double | |
| 90 two54 = 1.80143985094819840000e+16, /* 0x43500000, 0x00000000 */ | |
| 91 twom54 = 5.55111512312578270212e-17, /* 0x3C900000, 0x00000000 */ | |
| 92 hugeX = 1.0e+300, | |
| 93 tiny = 1.0e-300; | |
| 94 | |
| 95 static double scalbnA (double x, int n) { | |
| 96 int k,hx,lx; | |
| 97 hx = __HI(x); | |
| 98 lx = __LO(x); | |
| 99 k = (hx&0x7ff00000)>>20; /* extract exponent */ | |
| 100 if (k==0) { /* 0 or subnormal x */ | |
| 101 if ((lx|(hx&0x7fffffff))==0) return x; /* +-0 */ | |
| 102 x *= two54; | |
| 103 hx = __HI(x); | |
| 104 k = ((hx&0x7ff00000)>>20) - 54; | |
| 105 if (n< -50000) return tiny*x; /*underflow*/ | |
| 106 } | |
| 107 if (k==0x7ff) return x+x; /* NaN or Inf */ | |
| 108 k = k+n; | |
| 109 if (k > 0x7fe) return hugeX*copysignA(hugeX,x); /* overflow */ | |
| 110 if (k > 0) /* normal result */ | |
| 111 {__HI(x) = (hx&0x800fffff)|(k<<20); return x;} | |
| 112 if (k <= -54) { | |
| 113 if (n > 50000) /* in case integer overflow in n+k */ | |
| 114 return hugeX*copysignA(hugeX,x); /*overflow*/ | |
| 115 else return tiny*copysignA(tiny,x); /*underflow*/ | |
| 116 } | |
| 117 k += 54; /* subnormal result */ | |
| 118 __HI(x) = (hx&0x800fffff)|(k<<20); | |
| 119 return x*twom54; | |
| 120 } | |
| 121 | |
| 122 /* | |
| 123 * __kernel_rem_pio2(x,y,e0,nx,prec,ipio2) | |
| 124 * double x[],y[]; int e0,nx,prec; int ipio2[]; | |
| 125 * | |
| 126 * __kernel_rem_pio2 return the last three digits of N with | |
| 127 * y = x - N*pi/2 | |
| 128 * so that |y| < pi/2. | |
| 129 * | |
| 130 * The method is to compute the integer (mod 8) and fraction parts of | |
| 131 * (2/pi)*x without doing the full multiplication. In general we | |
| 132 * skip the part of the product that are known to be a huge integer ( | |
| 133 * more accurately, = 0 mod 8 ). Thus the number of operations are | |
| 134 * independent of the exponent of the input. | |
| 135 * | |
| 136 * (2/pi) is represented by an array of 24-bit integers in ipio2[]. | |
| 137 * | |
| 138 * Input parameters: | |
| 139 * x[] The input value (must be positive) is broken into nx | |
| 140 * pieces of 24-bit integers in double precision format. | |
| 141 * x[i] will be the i-th 24 bit of x. The scaled exponent | |
| 142 * of x[0] is given in input parameter e0 (i.e., x[0]*2^e0 | |
| 143 * match x's up to 24 bits. | |
| 144 * | |
| 145 * Example of breaking a double positive z into x[0]+x[1]+x[2]: | |
| 146 * e0 = ilogb(z)-23 | |
| 147 * z = scalbn(z,-e0) | |
| 148 * for i = 0,1,2 | |
| 149 * x[i] = floor(z) | |
| 150 * z = (z-x[i])*2**24 | |
| 151 * | |
| 152 * | |
| 153 * y[] ouput result in an array of double precision numbers. | |
| 154 * The dimension of y[] is: | |
| 155 * 24-bit precision 1 | |
| 156 * 53-bit precision 2 | |
| 157 * 64-bit precision 2 | |
| 158 * 113-bit precision 3 | |
| 159 * The actual value is the sum of them. Thus for 113-bit | |
| 160 * precsion, one may have to do something like: | |
| 161 * | |
| 162 * long double t,w,r_head, r_tail; | |
| 163 * t = (long double)y[2] + (long double)y[1]; | |
| 164 * w = (long double)y[0]; | |
| 165 * r_head = t+w; | |
| 166 * r_tail = w - (r_head - t); | |
| 167 * | |
| 168 * e0 The exponent of x[0] | |
| 169 * | |
| 170 * nx dimension of x[] | |
| 171 * | |
| 172 * prec an interger indicating the precision: | |
| 173 * 0 24 bits (single) | |
| 174 * 1 53 bits (double) | |
| 175 * 2 64 bits (extended) | |
| 176 * 3 113 bits (quad) | |
| 177 * | |
| 178 * ipio2[] | |
| 179 * integer array, contains the (24*i)-th to (24*i+23)-th | |
| 180 * bit of 2/pi after binary point. The corresponding | |
| 181 * floating value is | |
| 182 * | |
| 183 * ipio2[i] * 2^(-24(i+1)). | |
| 184 * | |
| 185 * External function: | |
| 186 * double scalbn(), floor(); | |
| 187 * | |
| 188 * | |
| 189 * Here is the description of some local variables: | |
| 190 * | |
| 191 * jk jk+1 is the initial number of terms of ipio2[] needed | |
| 192 * in the computation. The recommended value is 2,3,4, | |
| 193 * 6 for single, double, extended,and quad. | |
| 194 * | |
| 195 * jz local integer variable indicating the number of | |
| 196 * terms of ipio2[] used. | |
| 197 * | |
| 198 * jx nx - 1 | |
| 199 * | |
| 200 * jv index for pointing to the suitable ipio2[] for the | |
| 201 * computation. In general, we want | |
| 202 * ( 2^e0*x[0] * ipio2[jv-1]*2^(-24jv) )/8 | |
| 203 * is an integer. Thus | |
| 204 * e0-3-24*jv >= 0 or (e0-3)/24 >= jv | |
| 205 * Hence jv = max(0,(e0-3)/24). | |
| 206 * | |
| 207 * jp jp+1 is the number of terms in PIo2[] needed, jp = jk. | |
| 208 * | |
| 209 * q[] double array with integral value, representing the | |
| 210 * 24-bits chunk of the product of x and 2/pi. | |
| 211 * | |
| 212 * q0 the corresponding exponent of q[0]. Note that the | |
| 213 * exponent for q[i] would be q0-24*i. | |
| 214 * | |
| 215 * PIo2[] double precision array, obtained by cutting pi/2 | |
| 216 * into 24 bits chunks. | |
| 217 * | |
| 218 * f[] ipio2[] in floating point | |
| 219 * | |
| 220 * iq[] integer array by breaking up q[] in 24-bits chunk. | |
| 221 * | |
| 222 * fq[] final product of x*(2/pi) in fq[0],..,fq[jk] | |
| 223 * | |
| 224 * ih integer. If >0 it indicats q[] is >= 0.5, hence | |
| 225 * it also indicates the *sign* of the result. | |
| 226 * | |
| 227 */ | |
| 228 | |
| 229 | |
| 230 /* | |
| 231 * Constants: | |
| 232 * The hexadecimal values are the intended ones for the following | |
| 233 * constants. The decimal values may be used, provided that the | |
| 234 * compiler will convert from decimal to binary accurately enough | |
| 235 * to produce the hexadecimal values shown. | |
| 236 */ | |
| 237 | |
| 238 | |
| 239 static const int init_jk[] = {2,3,4,6}; /* initial value for jk */ | |
| 240 | |
| 241 static const double PIo2[] = { | |
| 242 1.57079625129699707031e+00, /* 0x3FF921FB, 0x40000000 */ | |
| 243 7.54978941586159635335e-08, /* 0x3E74442D, 0x00000000 */ | |
| 244 5.39030252995776476554e-15, /* 0x3CF84698, 0x80000000 */ | |
| 245 3.28200341580791294123e-22, /* 0x3B78CC51, 0x60000000 */ | |
| 246 1.27065575308067607349e-29, /* 0x39F01B83, 0x80000000 */ | |
| 247 1.22933308981111328932e-36, /* 0x387A2520, 0x40000000 */ | |
| 248 2.73370053816464559624e-44, /* 0x36E38222, 0x80000000 */ | |
| 249 2.16741683877804819444e-51, /* 0x3569F31D, 0x00000000 */ | |
| 250 }; | |
| 251 | |
| 252 static const double | |
| 253 zeroB = 0.0, | |
| 254 one = 1.0, | |
| 255 two24B = 1.67772160000000000000e+07, /* 0x41700000, 0x00000000 */ | |
| 256 twon24 = 5.96046447753906250000e-08; /* 0x3E700000, 0x00000000 */ | |
| 257 | |
| 1485 | 258 static SAFEBUF int __kernel_rem_pio2(double *x, double *y, int e0, int nx, int prec, const int *ipio2) { |
| 0 | 259 int jz,jx,jv,jp,jk,carry,n,iq[20],i,j,k,m,q0,ih; |
| 260 double z,fw,f[20],fq[20],q[20]; | |
| 261 | |
| 262 /* initialize jk*/ | |
| 263 jk = init_jk[prec]; | |
| 264 jp = jk; | |
| 265 | |
| 266 /* determine jx,jv,q0, note that 3>q0 */ | |
| 267 jx = nx-1; | |
| 268 jv = (e0-3)/24; if(jv<0) jv=0; | |
| 269 q0 = e0-24*(jv+1); | |
| 270 | |
| 271 /* set up f[0] to f[jx+jk] where f[jx+jk] = ipio2[jv+jk] */ | |
| 272 j = jv-jx; m = jx+jk; | |
| 273 for(i=0;i<=m;i++,j++) f[i] = (j<0)? zeroB : (double) ipio2[j]; | |
| 274 | |
| 275 /* compute q[0],q[1],...q[jk] */ | |
| 276 for (i=0;i<=jk;i++) { | |
| 277 for(j=0,fw=0.0;j<=jx;j++) fw += x[j]*f[jx+i-j]; q[i] = fw; | |
| 278 } | |
| 279 | |
| 280 jz = jk; | |
| 281 recompute: | |
| 282 /* distill q[] into iq[] reversingly */ | |
| 283 for(i=0,j=jz,z=q[jz];j>0;i++,j--) { | |
| 284 fw = (double)((int)(twon24* z)); | |
| 285 iq[i] = (int)(z-two24B*fw); | |
| 286 z = q[j-1]+fw; | |
| 287 } | |
| 288 | |
| 289 /* compute n */ | |
| 290 z = scalbnA(z,q0); /* actual value of z */ | |
| 291 z -= 8.0*floor(z*0.125); /* trim off integer >= 8 */ | |
| 292 n = (int) z; | |
| 293 z -= (double)n; | |
| 294 ih = 0; | |
| 295 if(q0>0) { /* need iq[jz-1] to determine n */ | |
| 296 i = (iq[jz-1]>>(24-q0)); n += i; | |
| 297 iq[jz-1] -= i<<(24-q0); | |
| 298 ih = iq[jz-1]>>(23-q0); | |
| 299 } | |
| 300 else if(q0==0) ih = iq[jz-1]>>23; | |
| 301 else if(z>=0.5) ih=2; | |
| 302 | |
| 303 if(ih>0) { /* q > 0.5 */ | |
| 304 n += 1; carry = 0; | |
| 305 for(i=0;i<jz ;i++) { /* compute 1-q */ | |
| 306 j = iq[i]; | |
| 307 if(carry==0) { | |
| 308 if(j!=0) { | |
| 309 carry = 1; iq[i] = 0x1000000- j; | |
| 310 } | |
| 311 } else iq[i] = 0xffffff - j; | |
| 312 } | |
| 313 if(q0>0) { /* rare case: chance is 1 in 12 */ | |
| 314 switch(q0) { | |
| 315 case 1: | |
| 316 iq[jz-1] &= 0x7fffff; break; | |
| 317 case 2: | |
| 318 iq[jz-1] &= 0x3fffff; break; | |
| 319 } | |
| 320 } | |
| 321 if(ih==2) { | |
| 322 z = one - z; | |
| 323 if(carry!=0) z -= scalbnA(one,q0); | |
| 324 } | |
| 325 } | |
| 326 | |
| 327 /* check if recomputation is needed */ | |
| 328 if(z==zeroB) { | |
| 329 j = 0; | |
| 330 for (i=jz-1;i>=jk;i--) j |= iq[i]; | |
| 331 if(j==0) { /* need recomputation */ | |
| 332 for(k=1;iq[jk-k]==0;k++); /* k = no. of terms needed */ | |
| 333 | |
| 334 for(i=jz+1;i<=jz+k;i++) { /* add q[jz+1] to q[jz+k] */ | |
| 335 f[jx+i] = (double) ipio2[jv+i]; | |
| 336 for(j=0,fw=0.0;j<=jx;j++) fw += x[j]*f[jx+i-j]; | |
| 337 q[i] = fw; | |
| 338 } | |
| 339 jz += k; | |
| 340 goto recompute; | |
| 341 } | |
| 342 } | |
| 343 | |
| 344 /* chop off zero terms */ | |
| 345 if(z==0.0) { | |
| 346 jz -= 1; q0 -= 24; | |
| 347 while(iq[jz]==0) { jz--; q0-=24;} | |
| 348 } else { /* break z into 24-bit if neccessary */ | |
| 349 z = scalbnA(z,-q0); | |
| 350 if(z>=two24B) { | |
| 351 fw = (double)((int)(twon24*z)); | |
| 352 iq[jz] = (int)(z-two24B*fw); | |
| 353 jz += 1; q0 += 24; | |
| 354 iq[jz] = (int) fw; | |
| 355 } else iq[jz] = (int) z ; | |
| 356 } | |
| 357 | |
| 358 /* convert integer "bit" chunk to floating-point value */ | |
| 359 fw = scalbnA(one,q0); | |
| 360 for(i=jz;i>=0;i--) { | |
| 361 q[i] = fw*(double)iq[i]; fw*=twon24; | |
| 362 } | |
| 363 | |
| 364 /* compute PIo2[0,...,jp]*q[jz,...,0] */ | |
| 365 for(i=jz;i>=0;i--) { | |
| 366 for(fw=0.0,k=0;k<=jp&&k<=jz-i;k++) fw += PIo2[k]*q[i+k]; | |
| 367 fq[jz-i] = fw; | |
| 368 } | |
| 369 | |
| 370 /* compress fq[] into y[] */ | |
| 371 switch(prec) { | |
| 372 case 0: | |
| 373 fw = 0.0; | |
| 374 for (i=jz;i>=0;i--) fw += fq[i]; | |
| 375 y[0] = (ih==0)? fw: -fw; | |
| 376 break; | |
| 377 case 1: | |
| 378 case 2: | |
| 379 fw = 0.0; | |
| 380 for (i=jz;i>=0;i--) fw += fq[i]; | |
| 381 y[0] = (ih==0)? fw: -fw; | |
| 382 fw = fq[0]-fw; | |
| 383 for (i=1;i<=jz;i++) fw += fq[i]; | |
| 384 y[1] = (ih==0)? fw: -fw; | |
| 385 break; | |
| 386 case 3: /* painful */ | |
| 387 for (i=jz;i>0;i--) { | |
| 388 fw = fq[i-1]+fq[i]; | |
| 389 fq[i] += fq[i-1]-fw; | |
| 390 fq[i-1] = fw; | |
| 391 } | |
| 392 for (i=jz;i>1;i--) { | |
| 393 fw = fq[i-1]+fq[i]; | |
| 394 fq[i] += fq[i-1]-fw; | |
| 395 fq[i-1] = fw; | |
| 396 } | |
| 397 for (fw=0.0,i=jz;i>=2;i--) fw += fq[i]; | |
| 398 if(ih==0) { | |
| 399 y[0] = fq[0]; y[1] = fq[1]; y[2] = fw; | |
| 400 } else { | |
| 401 y[0] = -fq[0]; y[1] = -fq[1]; y[2] = -fw; | |
| 402 } | |
| 403 } | |
| 404 return n&7; | |
| 405 } | |
| 406 | |
| 407 | |
| 408 /* | |
| 409 * ==================================================== | |
|
1552
c18cbe5936b8
6941466: Oracle rebranding changes for Hotspot repositories
trims
parents:
1485
diff
changeset
|
410 * Copyright (c) 1993 Oracle and/or its affilates. All rights reserved. |
| 0 | 411 * |
| 412 * Developed at SunPro, a Sun Microsystems, Inc. business. | |
| 413 * Permission to use, copy, modify, and distribute this | |
| 414 * software is freely granted, provided that this notice | |
| 415 * is preserved. | |
| 416 * ==================================================== | |
| 417 * | |
| 418 */ | |
| 419 | |
| 420 /* __ieee754_rem_pio2(x,y) | |
| 421 * | |
| 422 * return the remainder of x rem pi/2 in y[0]+y[1] | |
| 423 * use __kernel_rem_pio2() | |
| 424 */ | |
| 425 | |
| 426 /* | |
| 427 * Table of constants for 2/pi, 396 Hex digits (476 decimal) of 2/pi | |
| 428 */ | |
| 429 static const int two_over_pi[] = { | |
| 430 0xA2F983, 0x6E4E44, 0x1529FC, 0x2757D1, 0xF534DD, 0xC0DB62, | |
| 431 0x95993C, 0x439041, 0xFE5163, 0xABDEBB, 0xC561B7, 0x246E3A, | |
| 432 0x424DD2, 0xE00649, 0x2EEA09, 0xD1921C, 0xFE1DEB, 0x1CB129, | |
| 433 0xA73EE8, 0x8235F5, 0x2EBB44, 0x84E99C, 0x7026B4, 0x5F7E41, | |
| 434 0x3991D6, 0x398353, 0x39F49C, 0x845F8B, 0xBDF928, 0x3B1FF8, | |
| 435 0x97FFDE, 0x05980F, 0xEF2F11, 0x8B5A0A, 0x6D1F6D, 0x367ECF, | |
| 436 0x27CB09, 0xB74F46, 0x3F669E, 0x5FEA2D, 0x7527BA, 0xC7EBE5, | |
| 437 0xF17B3D, 0x0739F7, 0x8A5292, 0xEA6BFB, 0x5FB11F, 0x8D5D08, | |
| 438 0x560330, 0x46FC7B, 0x6BABF0, 0xCFBC20, 0x9AF436, 0x1DA9E3, | |
| 439 0x91615E, 0xE61B08, 0x659985, 0x5F14A0, 0x68408D, 0xFFD880, | |
| 440 0x4D7327, 0x310606, 0x1556CA, 0x73A8C9, 0x60E27B, 0xC08C6B, | |
| 441 }; | |
| 442 | |
| 443 static const int npio2_hw[] = { | |
| 444 0x3FF921FB, 0x400921FB, 0x4012D97C, 0x401921FB, 0x401F6A7A, 0x4022D97C, | |
| 445 0x4025FDBB, 0x402921FB, 0x402C463A, 0x402F6A7A, 0x4031475C, 0x4032D97C, | |
| 446 0x40346B9C, 0x4035FDBB, 0x40378FDB, 0x403921FB, 0x403AB41B, 0x403C463A, | |
| 447 0x403DD85A, 0x403F6A7A, 0x40407E4C, 0x4041475C, 0x4042106C, 0x4042D97C, | |
| 448 0x4043A28C, 0x40446B9C, 0x404534AC, 0x4045FDBB, 0x4046C6CB, 0x40478FDB, | |
| 449 0x404858EB, 0x404921FB, | |
| 450 }; | |
| 451 | |
| 452 /* | |
| 453 * invpio2: 53 bits of 2/pi | |
| 454 * pio2_1: first 33 bit of pi/2 | |
| 455 * pio2_1t: pi/2 - pio2_1 | |
| 456 * pio2_2: second 33 bit of pi/2 | |
| 457 * pio2_2t: pi/2 - (pio2_1+pio2_2) | |
| 458 * pio2_3: third 33 bit of pi/2 | |
| 459 * pio2_3t: pi/2 - (pio2_1+pio2_2+pio2_3) | |
| 460 */ | |
| 461 | |
| 462 static const double | |
| 463 zeroA = 0.00000000000000000000e+00, /* 0x00000000, 0x00000000 */ | |
| 464 half = 5.00000000000000000000e-01, /* 0x3FE00000, 0x00000000 */ | |
| 465 two24A = 1.67772160000000000000e+07, /* 0x41700000, 0x00000000 */ | |
| 466 invpio2 = 6.36619772367581382433e-01, /* 0x3FE45F30, 0x6DC9C883 */ | |
| 467 pio2_1 = 1.57079632673412561417e+00, /* 0x3FF921FB, 0x54400000 */ | |
| 468 pio2_1t = 6.07710050650619224932e-11, /* 0x3DD0B461, 0x1A626331 */ | |
| 469 pio2_2 = 6.07710050630396597660e-11, /* 0x3DD0B461, 0x1A600000 */ | |
| 470 pio2_2t = 2.02226624879595063154e-21, /* 0x3BA3198A, 0x2E037073 */ | |
| 471 pio2_3 = 2.02226624871116645580e-21, /* 0x3BA3198A, 0x2E000000 */ | |
| 472 pio2_3t = 8.47842766036889956997e-32; /* 0x397B839A, 0x252049C1 */ | |
| 473 | |
| 1485 | 474 static SAFEBUF int __ieee754_rem_pio2(double x, double *y) { |
| 0 | 475 double z,w,t,r,fn; |
| 476 double tx[3]; | |
| 477 int e0,i,j,nx,n,ix,hx,i0; | |
| 478 | |
| 479 i0 = ((*(int*)&two24A)>>30)^1; /* high word index */ | |
| 480 hx = *(i0+(int*)&x); /* high word of x */ | |
| 481 ix = hx&0x7fffffff; | |
| 482 if(ix<=0x3fe921fb) /* |x| ~<= pi/4 , no need for reduction */ | |
| 483 {y[0] = x; y[1] = 0; return 0;} | |
| 484 if(ix<0x4002d97c) { /* |x| < 3pi/4, special case with n=+-1 */ | |
| 485 if(hx>0) { | |
| 486 z = x - pio2_1; | |
| 487 if(ix!=0x3ff921fb) { /* 33+53 bit pi is good enough */ | |
| 488 y[0] = z - pio2_1t; | |
| 489 y[1] = (z-y[0])-pio2_1t; | |
| 490 } else { /* near pi/2, use 33+33+53 bit pi */ | |
| 491 z -= pio2_2; | |
| 492 y[0] = z - pio2_2t; | |
| 493 y[1] = (z-y[0])-pio2_2t; | |
| 494 } | |
| 495 return 1; | |
| 496 } else { /* negative x */ | |
| 497 z = x + pio2_1; | |
| 498 if(ix!=0x3ff921fb) { /* 33+53 bit pi is good enough */ | |
| 499 y[0] = z + pio2_1t; | |
| 500 y[1] = (z-y[0])+pio2_1t; | |
| 501 } else { /* near pi/2, use 33+33+53 bit pi */ | |
| 502 z += pio2_2; | |
| 503 y[0] = z + pio2_2t; | |
| 504 y[1] = (z-y[0])+pio2_2t; | |
| 505 } | |
| 506 return -1; | |
| 507 } | |
| 508 } | |
| 509 if(ix<=0x413921fb) { /* |x| ~<= 2^19*(pi/2), medium size */ | |
| 510 t = fabsd(x); | |
| 511 n = (int) (t*invpio2+half); | |
| 512 fn = (double)n; | |
| 513 r = t-fn*pio2_1; | |
| 514 w = fn*pio2_1t; /* 1st round good to 85 bit */ | |
| 515 if(n<32&&ix!=npio2_hw[n-1]) { | |
| 516 y[0] = r-w; /* quick check no cancellation */ | |
| 517 } else { | |
| 518 j = ix>>20; | |
| 519 y[0] = r-w; | |
| 520 i = j-(((*(i0+(int*)&y[0]))>>20)&0x7ff); | |
| 521 if(i>16) { /* 2nd iteration needed, good to 118 */ | |
| 522 t = r; | |
| 523 w = fn*pio2_2; | |
| 524 r = t-w; | |
| 525 w = fn*pio2_2t-((t-r)-w); | |
| 526 y[0] = r-w; | |
| 527 i = j-(((*(i0+(int*)&y[0]))>>20)&0x7ff); | |
| 528 if(i>49) { /* 3rd iteration need, 151 bits acc */ | |
| 529 t = r; /* will cover all possible cases */ | |
| 530 w = fn*pio2_3; | |
| 531 r = t-w; | |
| 532 w = fn*pio2_3t-((t-r)-w); | |
| 533 y[0] = r-w; | |
| 534 } | |
| 535 } | |
| 536 } | |
| 537 y[1] = (r-y[0])-w; | |
| 538 if(hx<0) {y[0] = -y[0]; y[1] = -y[1]; return -n;} | |
| 539 else return n; | |
| 540 } | |
| 541 /* | |
| 542 * all other (large) arguments | |
| 543 */ | |
| 544 if(ix>=0x7ff00000) { /* x is inf or NaN */ | |
| 545 y[0]=y[1]=x-x; return 0; | |
| 546 } | |
| 547 /* set z = scalbn(|x|,ilogb(x)-23) */ | |
| 548 *(1-i0+(int*)&z) = *(1-i0+(int*)&x); | |
| 549 e0 = (ix>>20)-1046; /* e0 = ilogb(z)-23; */ | |
| 550 *(i0+(int*)&z) = ix - (e0<<20); | |
| 551 for(i=0;i<2;i++) { | |
| 552 tx[i] = (double)((int)(z)); | |
| 553 z = (z-tx[i])*two24A; | |
| 554 } | |
| 555 tx[2] = z; | |
| 556 nx = 3; | |
| 557 while(tx[nx-1]==zeroA) nx--; /* skip zero term */ | |
| 558 n = __kernel_rem_pio2(tx,y,e0,nx,2,two_over_pi); | |
| 559 if(hx<0) {y[0] = -y[0]; y[1] = -y[1]; return -n;} | |
| 560 return n; | |
| 561 } | |
| 562 | |
| 563 | |
| 564 /* __kernel_sin( x, y, iy) | |
| 565 * kernel sin function on [-pi/4, pi/4], pi/4 ~ 0.7854 | |
| 566 * Input x is assumed to be bounded by ~pi/4 in magnitude. | |
| 567 * Input y is the tail of x. | |
| 568 * Input iy indicates whether y is 0. (if iy=0, y assume to be 0). | |
| 569 * | |
| 570 * Algorithm | |
| 571 * 1. Since sin(-x) = -sin(x), we need only to consider positive x. | |
| 572 * 2. if x < 2^-27 (hx<0x3e400000 0), return x with inexact if x!=0. | |
| 573 * 3. sin(x) is approximated by a polynomial of degree 13 on | |
| 574 * [0,pi/4] | |
| 575 * 3 13 | |
| 576 * sin(x) ~ x + S1*x + ... + S6*x | |
| 577 * where | |
| 578 * | |
| 579 * |sin(x) 2 4 6 8 10 12 | -58 | |
| 580 * |----- - (1+S1*x +S2*x +S3*x +S4*x +S5*x +S6*x )| <= 2 | |
| 581 * | x | | |
| 582 * | |
| 583 * 4. sin(x+y) = sin(x) + sin'(x')*y | |
| 584 * ~ sin(x) + (1-x*x/2)*y | |
| 585 * For better accuracy, let | |
| 586 * 3 2 2 2 2 | |
| 587 * r = x *(S2+x *(S3+x *(S4+x *(S5+x *S6)))) | |
| 588 * then 3 2 | |
| 589 * sin(x) = x + (S1*x + (x *(r-y/2)+y)) | |
| 590 */ | |
| 591 | |
| 592 static const double | |
| 593 S1 = -1.66666666666666324348e-01, /* 0xBFC55555, 0x55555549 */ | |
| 594 S2 = 8.33333333332248946124e-03, /* 0x3F811111, 0x1110F8A6 */ | |
| 595 S3 = -1.98412698298579493134e-04, /* 0xBF2A01A0, 0x19C161D5 */ | |
| 596 S4 = 2.75573137070700676789e-06, /* 0x3EC71DE3, 0x57B1FE7D */ | |
| 597 S5 = -2.50507602534068634195e-08, /* 0xBE5AE5E6, 0x8A2B9CEB */ | |
| 598 S6 = 1.58969099521155010221e-10; /* 0x3DE5D93A, 0x5ACFD57C */ | |
| 599 | |
| 600 static double __kernel_sin(double x, double y, int iy) | |
| 601 { | |
| 602 double z,r,v; | |
| 603 int ix; | |
| 604 ix = __HI(x)&0x7fffffff; /* high word of x */ | |
| 605 if(ix<0x3e400000) /* |x| < 2**-27 */ | |
| 606 {if((int)x==0) return x;} /* generate inexact */ | |
| 607 z = x*x; | |
| 608 v = z*x; | |
| 609 r = S2+z*(S3+z*(S4+z*(S5+z*S6))); | |
| 610 if(iy==0) return x+v*(S1+z*r); | |
| 611 else return x-((z*(half*y-v*r)-y)-v*S1); | |
| 612 } | |
| 613 | |
| 614 /* | |
| 615 * __kernel_cos( x, y ) | |
| 616 * kernel cos function on [-pi/4, pi/4], pi/4 ~ 0.785398164 | |
| 617 * Input x is assumed to be bounded by ~pi/4 in magnitude. | |
| 618 * Input y is the tail of x. | |
| 619 * | |
| 620 * Algorithm | |
| 621 * 1. Since cos(-x) = cos(x), we need only to consider positive x. | |
| 622 * 2. if x < 2^-27 (hx<0x3e400000 0), return 1 with inexact if x!=0. | |
| 623 * 3. cos(x) is approximated by a polynomial of degree 14 on | |
| 624 * [0,pi/4] | |
| 625 * 4 14 | |
| 626 * cos(x) ~ 1 - x*x/2 + C1*x + ... + C6*x | |
| 627 * where the remez error is | |
| 628 * | |
| 629 * | 2 4 6 8 10 12 14 | -58 | |
| 630 * |cos(x)-(1-.5*x +C1*x +C2*x +C3*x +C4*x +C5*x +C6*x )| <= 2 | |
| 631 * | | | |
| 632 * | |
| 633 * 4 6 8 10 12 14 | |
| 634 * 4. let r = C1*x +C2*x +C3*x +C4*x +C5*x +C6*x , then | |
| 635 * cos(x) = 1 - x*x/2 + r | |
| 636 * since cos(x+y) ~ cos(x) - sin(x)*y | |
| 637 * ~ cos(x) - x*y, | |
| 638 * a correction term is necessary in cos(x) and hence | |
| 639 * cos(x+y) = 1 - (x*x/2 - (r - x*y)) | |
| 640 * For better accuracy when x > 0.3, let qx = |x|/4 with | |
| 641 * the last 32 bits mask off, and if x > 0.78125, let qx = 0.28125. | |
| 642 * Then | |
| 643 * cos(x+y) = (1-qx) - ((x*x/2-qx) - (r-x*y)). | |
| 644 * Note that 1-qx and (x*x/2-qx) is EXACT here, and the | |
| 645 * magnitude of the latter is at least a quarter of x*x/2, | |
| 646 * thus, reducing the rounding error in the subtraction. | |
| 647 */ | |
| 648 | |
| 649 static const double | |
| 650 C1 = 4.16666666666666019037e-02, /* 0x3FA55555, 0x5555554C */ | |
| 651 C2 = -1.38888888888741095749e-03, /* 0xBF56C16C, 0x16C15177 */ | |
| 652 C3 = 2.48015872894767294178e-05, /* 0x3EFA01A0, 0x19CB1590 */ | |
| 653 C4 = -2.75573143513906633035e-07, /* 0xBE927E4F, 0x809C52AD */ | |
| 654 C5 = 2.08757232129817482790e-09, /* 0x3E21EE9E, 0xBDB4B1C4 */ | |
| 655 C6 = -1.13596475577881948265e-11; /* 0xBDA8FAE9, 0xBE8838D4 */ | |
| 656 | |
| 657 static double __kernel_cos(double x, double y) | |
| 658 { | |
| 659 double a,hz,z,r,qx; | |
| 660 int ix; | |
| 661 ix = __HI(x)&0x7fffffff; /* ix = |x|'s high word*/ | |
| 662 if(ix<0x3e400000) { /* if x < 2**27 */ | |
| 663 if(((int)x)==0) return one; /* generate inexact */ | |
| 664 } | |
| 665 z = x*x; | |
| 666 r = z*(C1+z*(C2+z*(C3+z*(C4+z*(C5+z*C6))))); | |
| 667 if(ix < 0x3FD33333) /* if |x| < 0.3 */ | |
| 668 return one - (0.5*z - (z*r - x*y)); | |
| 669 else { | |
| 670 if(ix > 0x3fe90000) { /* x > 0.78125 */ | |
| 671 qx = 0.28125; | |
| 672 } else { | |
| 673 __HI(qx) = ix-0x00200000; /* x/4 */ | |
| 674 __LO(qx) = 0; | |
| 675 } | |
| 676 hz = 0.5*z-qx; | |
| 677 a = one-qx; | |
| 678 return a - (hz - (z*r-x*y)); | |
| 679 } | |
| 680 } | |
| 681 | |
| 682 /* __kernel_tan( x, y, k ) | |
| 683 * kernel tan function on [-pi/4, pi/4], pi/4 ~ 0.7854 | |
| 684 * Input x is assumed to be bounded by ~pi/4 in magnitude. | |
| 685 * Input y is the tail of x. | |
| 686 * Input k indicates whether tan (if k=1) or | |
| 687 * -1/tan (if k= -1) is returned. | |
| 688 * | |
| 689 * Algorithm | |
| 690 * 1. Since tan(-x) = -tan(x), we need only to consider positive x. | |
| 691 * 2. if x < 2^-28 (hx<0x3e300000 0), return x with inexact if x!=0. | |
| 692 * 3. tan(x) is approximated by a odd polynomial of degree 27 on | |
| 693 * [0,0.67434] | |
| 694 * 3 27 | |
| 695 * tan(x) ~ x + T1*x + ... + T13*x | |
| 696 * where | |
| 697 * | |
| 698 * |tan(x) 2 4 26 | -59.2 | |
| 699 * |----- - (1+T1*x +T2*x +.... +T13*x )| <= 2 | |
| 700 * | x | | |
| 701 * | |
| 702 * Note: tan(x+y) = tan(x) + tan'(x)*y | |
| 703 * ~ tan(x) + (1+x*x)*y | |
| 704 * Therefore, for better accuracy in computing tan(x+y), let | |
| 705 * 3 2 2 2 2 | |
| 706 * r = x *(T2+x *(T3+x *(...+x *(T12+x *T13)))) | |
| 707 * then | |
| 708 * 3 2 | |
| 709 * tan(x+y) = x + (T1*x + (x *(r+y)+y)) | |
| 710 * | |
| 711 * 4. For x in [0.67434,pi/4], let y = pi/4 - x, then | |
| 712 * tan(x) = tan(pi/4-y) = (1-tan(y))/(1+tan(y)) | |
| 713 * = 1 - 2*(tan(y) - (tan(y)^2)/(1+tan(y))) | |
| 714 */ | |
| 715 | |
| 716 static const double | |
| 717 pio4 = 7.85398163397448278999e-01, /* 0x3FE921FB, 0x54442D18 */ | |
| 718 pio4lo= 3.06161699786838301793e-17, /* 0x3C81A626, 0x33145C07 */ | |
| 719 T[] = { | |
| 720 3.33333333333334091986e-01, /* 0x3FD55555, 0x55555563 */ | |
| 721 1.33333333333201242699e-01, /* 0x3FC11111, 0x1110FE7A */ | |
| 722 5.39682539762260521377e-02, /* 0x3FABA1BA, 0x1BB341FE */ | |
| 723 2.18694882948595424599e-02, /* 0x3F9664F4, 0x8406D637 */ | |
| 724 8.86323982359930005737e-03, /* 0x3F8226E3, 0xE96E8493 */ | |
| 725 3.59207910759131235356e-03, /* 0x3F6D6D22, 0xC9560328 */ | |
| 726 1.45620945432529025516e-03, /* 0x3F57DBC8, 0xFEE08315 */ | |
| 727 5.88041240820264096874e-04, /* 0x3F4344D8, 0xF2F26501 */ | |
| 728 2.46463134818469906812e-04, /* 0x3F3026F7, 0x1A8D1068 */ | |
| 729 7.81794442939557092300e-05, /* 0x3F147E88, 0xA03792A6 */ | |
| 730 7.14072491382608190305e-05, /* 0x3F12B80F, 0x32F0A7E9 */ | |
| 731 -1.85586374855275456654e-05, /* 0xBEF375CB, 0xDB605373 */ | |
| 732 2.59073051863633712884e-05, /* 0x3EFB2A70, 0x74BF7AD4 */ | |
| 733 }; | |
| 734 | |
| 735 static double __kernel_tan(double x, double y, int iy) | |
| 736 { | |
| 737 double z,r,v,w,s; | |
| 738 int ix,hx; | |
| 739 hx = __HI(x); /* high word of x */ | |
| 740 ix = hx&0x7fffffff; /* high word of |x| */ | |
| 741 if(ix<0x3e300000) { /* x < 2**-28 */ | |
| 742 if((int)x==0) { /* generate inexact */ | |
| 743 if (((ix | __LO(x)) | (iy + 1)) == 0) | |
| 744 return one / fabsd(x); | |
| 745 else { | |
| 746 if (iy == 1) | |
| 747 return x; | |
| 748 else { /* compute -1 / (x+y) carefully */ | |
| 749 double a, t; | |
| 750 | |
| 751 z = w = x + y; | |
| 752 __LO(z) = 0; | |
| 753 v = y - (z - x); | |
| 754 t = a = -one / w; | |
| 755 __LO(t) = 0; | |
| 756 s = one + t * z; | |
| 757 return t + a * (s + t * v); | |
| 758 } | |
| 759 } | |
| 760 } | |
| 761 } | |
| 762 if(ix>=0x3FE59428) { /* |x|>=0.6744 */ | |
| 763 if(hx<0) {x = -x; y = -y;} | |
| 764 z = pio4-x; | |
| 765 w = pio4lo-y; | |
| 766 x = z+w; y = 0.0; | |
| 767 } | |
| 768 z = x*x; | |
| 769 w = z*z; | |
| 770 /* Break x^5*(T[1]+x^2*T[2]+...) into | |
| 771 * x^5(T[1]+x^4*T[3]+...+x^20*T[11]) + | |
| 772 * x^5(x^2*(T[2]+x^4*T[4]+...+x^22*[T12])) | |
| 773 */ | |
| 774 r = T[1]+w*(T[3]+w*(T[5]+w*(T[7]+w*(T[9]+w*T[11])))); | |
| 775 v = z*(T[2]+w*(T[4]+w*(T[6]+w*(T[8]+w*(T[10]+w*T[12]))))); | |
| 776 s = z*x; | |
| 777 r = y + z*(s*(r+v)+y); | |
| 778 r += T[0]*s; | |
| 779 w = x+r; | |
| 780 if(ix>=0x3FE59428) { | |
| 781 v = (double)iy; | |
| 782 return (double)(1-((hx>>30)&2))*(v-2.0*(x-(w*w/(w+v)-r))); | |
| 783 } | |
| 784 if(iy==1) return w; | |
| 785 else { /* if allow error up to 2 ulp, | |
| 786 simply return -1.0/(x+r) here */ | |
| 787 /* compute -1.0/(x+r) accurately */ | |
| 788 double a,t; | |
| 789 z = w; | |
| 790 __LO(z) = 0; | |
| 791 v = r-(z - x); /* z+v = r+x */ | |
| 792 t = a = -1.0/w; /* a = -1.0/w */ | |
| 793 __LO(t) = 0; | |
| 794 s = 1.0+t*z; | |
| 795 return t+a*(s+t*v); | |
| 796 } | |
| 797 } | |
| 798 | |
| 799 | |
| 800 //---------------------------------------------------------------------- | |
| 801 // | |
| 802 // Routines for new sin/cos implementation | |
| 803 // | |
| 804 //---------------------------------------------------------------------- | |
| 805 | |
| 806 /* sin(x) | |
| 807 * Return sine function of x. | |
| 808 * | |
| 809 * kernel function: | |
| 810 * __kernel_sin ... sine function on [-pi/4,pi/4] | |
| 811 * __kernel_cos ... cose function on [-pi/4,pi/4] | |
| 812 * __ieee754_rem_pio2 ... argument reduction routine | |
| 813 * | |
| 814 * Method. | |
| 815 * Let S,C and T denote the sin, cos and tan respectively on | |
| 816 * [-PI/4, +PI/4]. Reduce the argument x to y1+y2 = x-k*pi/2 | |
| 817 * in [-pi/4 , +pi/4], and let n = k mod 4. | |
| 818 * We have | |
| 819 * | |
| 820 * n sin(x) cos(x) tan(x) | |
| 821 * ---------------------------------------------------------- | |
| 822 * 0 S C T | |
| 823 * 1 C -S -1/T | |
| 824 * 2 -S -C T | |
| 825 * 3 -C S -1/T | |
| 826 * ---------------------------------------------------------- | |
| 827 * | |
| 828 * Special cases: | |
| 829 * Let trig be any of sin, cos, or tan. | |
| 830 * trig(+-INF) is NaN, with signals; | |
| 831 * trig(NaN) is that NaN; | |
| 832 * | |
| 833 * Accuracy: | |
| 834 * TRIG(x) returns trig(x) nearly rounded | |
| 835 */ | |
| 836 | |
| 837 JRT_LEAF(jdouble, SharedRuntime::dsin(jdouble x)) | |
| 838 double y[2],z=0.0; | |
| 839 int n, ix; | |
| 840 | |
| 841 /* High word of x. */ | |
| 842 ix = __HI(x); | |
| 843 | |
| 844 /* |x| ~< pi/4 */ | |
| 845 ix &= 0x7fffffff; | |
| 846 if(ix <= 0x3fe921fb) return __kernel_sin(x,z,0); | |
| 847 | |
| 848 /* sin(Inf or NaN) is NaN */ | |
| 849 else if (ix>=0x7ff00000) return x-x; | |
| 850 | |
| 851 /* argument reduction needed */ | |
| 852 else { | |
| 853 n = __ieee754_rem_pio2(x,y); | |
| 854 switch(n&3) { | |
| 855 case 0: return __kernel_sin(y[0],y[1],1); | |
| 856 case 1: return __kernel_cos(y[0],y[1]); | |
| 857 case 2: return -__kernel_sin(y[0],y[1],1); | |
| 858 default: | |
| 859 return -__kernel_cos(y[0],y[1]); | |
| 860 } | |
| 861 } | |
| 862 JRT_END | |
| 863 | |
| 864 /* cos(x) | |
| 865 * Return cosine function of x. | |
| 866 * | |
| 867 * kernel function: | |
| 868 * __kernel_sin ... sine function on [-pi/4,pi/4] | |
| 869 * __kernel_cos ... cosine function on [-pi/4,pi/4] | |
| 870 * __ieee754_rem_pio2 ... argument reduction routine | |
| 871 * | |
| 872 * Method. | |
| 873 * Let S,C and T denote the sin, cos and tan respectively on | |
| 874 * [-PI/4, +PI/4]. Reduce the argument x to y1+y2 = x-k*pi/2 | |
| 875 * in [-pi/4 , +pi/4], and let n = k mod 4. | |
| 876 * We have | |
| 877 * | |
| 878 * n sin(x) cos(x) tan(x) | |
| 879 * ---------------------------------------------------------- | |
| 880 * 0 S C T | |
| 881 * 1 C -S -1/T | |
| 882 * 2 -S -C T | |
| 883 * 3 -C S -1/T | |
| 884 * ---------------------------------------------------------- | |
| 885 * | |
| 886 * Special cases: | |
| 887 * Let trig be any of sin, cos, or tan. | |
| 888 * trig(+-INF) is NaN, with signals; | |
| 889 * trig(NaN) is that NaN; | |
| 890 * | |
| 891 * Accuracy: | |
| 892 * TRIG(x) returns trig(x) nearly rounded | |
| 893 */ | |
| 894 | |
| 895 JRT_LEAF(jdouble, SharedRuntime::dcos(jdouble x)) | |
| 896 double y[2],z=0.0; | |
| 897 int n, ix; | |
| 898 | |
| 899 /* High word of x. */ | |
| 900 ix = __HI(x); | |
| 901 | |
| 902 /* |x| ~< pi/4 */ | |
| 903 ix &= 0x7fffffff; | |
| 904 if(ix <= 0x3fe921fb) return __kernel_cos(x,z); | |
| 905 | |
| 906 /* cos(Inf or NaN) is NaN */ | |
| 907 else if (ix>=0x7ff00000) return x-x; | |
| 908 | |
| 909 /* argument reduction needed */ | |
| 910 else { | |
| 911 n = __ieee754_rem_pio2(x,y); | |
| 912 switch(n&3) { | |
| 913 case 0: return __kernel_cos(y[0],y[1]); | |
| 914 case 1: return -__kernel_sin(y[0],y[1],1); | |
| 915 case 2: return -__kernel_cos(y[0],y[1]); | |
| 916 default: | |
| 917 return __kernel_sin(y[0],y[1],1); | |
| 918 } | |
| 919 } | |
| 920 JRT_END | |
| 921 | |
| 922 /* tan(x) | |
| 923 * Return tangent function of x. | |
| 924 * | |
| 925 * kernel function: | |
| 926 * __kernel_tan ... tangent function on [-pi/4,pi/4] | |
| 927 * __ieee754_rem_pio2 ... argument reduction routine | |
| 928 * | |
| 929 * Method. | |
| 930 * Let S,C and T denote the sin, cos and tan respectively on | |
| 931 * [-PI/4, +PI/4]. Reduce the argument x to y1+y2 = x-k*pi/2 | |
| 932 * in [-pi/4 , +pi/4], and let n = k mod 4. | |
| 933 * We have | |
| 934 * | |
| 935 * n sin(x) cos(x) tan(x) | |
| 936 * ---------------------------------------------------------- | |
| 937 * 0 S C T | |
| 938 * 1 C -S -1/T | |
| 939 * 2 -S -C T | |
| 940 * 3 -C S -1/T | |
| 941 * ---------------------------------------------------------- | |
| 942 * | |
| 943 * Special cases: | |
| 944 * Let trig be any of sin, cos, or tan. | |
| 945 * trig(+-INF) is NaN, with signals; | |
| 946 * trig(NaN) is that NaN; | |
| 947 * | |
| 948 * Accuracy: | |
| 949 * TRIG(x) returns trig(x) nearly rounded | |
| 950 */ | |
| 951 | |
| 952 JRT_LEAF(jdouble, SharedRuntime::dtan(jdouble x)) | |
| 953 double y[2],z=0.0; | |
| 954 int n, ix; | |
| 955 | |
| 956 /* High word of x. */ | |
| 957 ix = __HI(x); | |
| 958 | |
| 959 /* |x| ~< pi/4 */ | |
| 960 ix &= 0x7fffffff; | |
| 961 if(ix <= 0x3fe921fb) return __kernel_tan(x,z,1); | |
| 962 | |
| 963 /* tan(Inf or NaN) is NaN */ | |
| 964 else if (ix>=0x7ff00000) return x-x; /* NaN */ | |
| 965 | |
| 966 /* argument reduction needed */ | |
| 967 else { | |
| 968 n = __ieee754_rem_pio2(x,y); | |
| 969 return __kernel_tan(y[0],y[1],1-((n&1)<<1)); /* 1 -- n even | |
| 970 -1 -- n odd */ | |
| 971 } | |
| 972 JRT_END | |
| 973 | |
| 974 | |
| 975 #ifdef WIN32 | |
| 976 # pragma optimize ( "", on ) | |
| 977 #endif |