1 /* 2 * Licensed to the Apache Software Foundation (ASF) under one or more 3 * contributor license agreements. See the NOTICE file distributed with 4 * this work for additional information regarding copyright ownership. 5 * The ASF licenses this file to You under the Apache License, Version 2.0 6 * (the "License"); you may not use this file except in compliance with 7 * the License. You may obtain a copy of the License at 8 * 9 * http://www.apache.org/licenses/LICENSE-2.0 10 * 11 * Unless required by applicable law or agreed to in writing, software 12 * distributed under the License is distributed on an "AS IS" BASIS, 13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 14 * See the License for the specific language governing permissions and 15 * limitations under the License. 16 */ 17 18 /* 19 * acos, asin, atan, cosh, sinh, tanh, exp, expm1, log, log10, log1p, and cbrt 20 * have been implemented with the following license. 21 * ==================================================== 22 * Copyright (C) 2004 by Sun Microsystems, Inc. All rights reserved. 23 * 24 * Permission to use, copy, modify, and distribute this 25 * software is freely granted, provided that this notice 26 * is preserved. 27 * ==================================================== 28 */ 29 30 package java.lang; 31 32 /** 33 * Class StrictMath provides basic math constants and operations such as 34 * trigonometric functions, hyperbolic functions, exponential, logarithms, etc. 35 * <p> 36 * In contrast to class {@link Math}, the methods in this class return exactly 37 * the same results on all platforms. Algorithms based on these methods thus 38 * behave the same (e.g. regarding numerical convergence) on all platforms, 39 * complying with the slogan "write once, run everywhere". On the other side, 40 * the implementation of class StrictMath may be less efficient than that of 41 * class Math, as class StrictMath cannot utilize platform specific features 42 * such as an extended precision math co-processors. 43 * <p> 44 * The methods in this class are specified using the "Freely Distributable Math 45 * Library" (fdlibm), version 5.3. 46 * <p> 47 * <a href="http://www.netlib.org/fdlibm/">http://www.netlib.org/fdlibm/</a> 48 */ 49 public final class StrictMath { 50 /** 51 * The double value closest to e, the base of the natural logarithm. 52 */ 53 public static final double E = Math.E; 54 55 /** 56 * The double value closest to pi, the ratio of a circle's circumference to 57 * its diameter. 58 */ 59 public static final double PI = Math.PI; 60 61 /** 62 * Prevents this class from being instantiated. 63 */ 64 private StrictMath() { 65 } 66 67 /** 68 * Returns the absolute value of the argument. 69 * <p> 70 * Special cases: 71 * <ul> 72 * <li>{@code abs(-0.0) = +0.0}</li> 73 * <li>{@code abs(+infinity) = +infinity}</li> 74 * <li>{@code abs(-infinity) = +infinity}</li> 75 * <li>{@code abs(NaN) = NaN}</li> 76 * </ul> 77 */ 78 public static double abs(double d) { 79 return Math.abs(d); 80 } 81 82 /** 83 * Returns the absolute value of the argument. 84 * <p> 85 * Special cases: 86 * <ul> 87 * <li>{@code abs(-0.0) = +0.0}</li> 88 * <li>{@code abs(+infinity) = +infinity}</li> 89 * <li>{@code abs(-infinity) = +infinity}</li> 90 * <li>{@code abs(NaN) = NaN}</li> 91 * </ul> 92 */ 93 public static float abs(float f) { 94 return Math.abs(f); 95 } 96 97 /** 98 * Returns the absolute value of the argument. 99 * <p> 100 * If the argument is {@code Integer.MIN_VALUE}, {@code Integer.MIN_VALUE} 101 * is returned. 102 */ 103 public static int abs(int i) { 104 return Math.abs(i); 105 } 106 107 /** 108 * Returns the absolute value of the argument. 109 * <p> 110 * If the argument is {@code Long.MIN_VALUE}, {@code Long.MIN_VALUE} is 111 * returned. 112 */ 113 public static long abs(long l) { 114 return Math.abs(l); 115 } 116 117 private static final double PIO2_HI = 1.57079632679489655800e+00; 118 private static final double PIO2_LO = 6.12323399573676603587e-17; 119 private static final double PS0 = 1.66666666666666657415e-01; 120 private static final double PS1 = -3.25565818622400915405e-01; 121 private static final double PS2 = 2.01212532134862925881e-01; 122 private static final double PS3 = -4.00555345006794114027e-02; 123 private static final double PS4 = 7.91534994289814532176e-04; 124 private static final double PS5 = 3.47933107596021167570e-05; 125 private static final double QS1 = -2.40339491173441421878e+00; 126 private static final double QS2 = 2.02094576023350569471e+00; 127 private static final double QS3 = -6.88283971605453293030e-01; 128 private static final double QS4 = 7.70381505559019352791e-02; 129 private static final double HUGE = 1.000e+300; 130 private static final double PIO4_HI = 7.85398163397448278999e-01; 131 132 /** 133 * Returns the closest double approximation of the arc cosine of the 134 * argument within the range {@code [0..pi]}. 135 * <p> 136 * Special cases: 137 * <ul> 138 * <li>{@code acos((anything > 1) = NaN}</li> 139 * <li>{@code acos((anything < -1) = NaN}</li> 140 * <li>{@code acos(NaN) = NaN}</li> 141 * </ul> 142 * 143 * @param x 144 * the value to compute arc cosine of. 145 * @return the arc cosine of the argument. 146 */ 147 public static double acos(double x) { 148 double z, p, q, r, w, s, c, df; 149 int hx, ix; 150 final long bits = Double.doubleToRawLongBits(x); 151 hx = (int) (bits >>> 32); 152 ix = hx & 0x7fffffff; 153 if (ix >= 0x3ff00000) { /* |x| >= 1 */ 154 if ((((ix - 0x3ff00000) | ((int) bits))) == 0) { /* |x|==1 */ 155 if (hx > 0) { 156 return 0.0; /* ieee_acos(1) = 0 */ 157 } else { 158 return 3.14159265358979311600e+00 + 2.0 * PIO2_LO; /* ieee_acos(-1)= pi */ 159 } 160 } 161 return (x - x) / (x - x); /* ieee_acos(|x|>1) is NaN */ 162 } 163 164 if (ix < 0x3fe00000) { /* |x| < 0.5 */ 165 if (ix <= 0x3c600000) { 166 return PIO2_HI + PIO2_LO;/* if|x|<2**-57 */ 167 } 168 169 z = x * x; 170 p = z * (PS0 + z 171 * (PS1 + z * (PS2 + z * (PS3 + z * (PS4 + z * PS5))))); 172 q = 1.00000000000000000000e+00 + z * (QS1 + z * (QS2 + z * (QS3 + z * QS4))); 173 r = p / q; 174 return PIO2_HI - (x - (PIO2_LO - x * r)); 175 } else if (hx < 0) { /* x < -0.5 */ 176 z = (1.00000000000000000000e+00 + x) * 0.5; 177 p = z * (PS0 + z 178 * (PS1 + z * (PS2 + z * (PS3 + z * (PS4 + z * PS5))))); 179 q = 1.00000000000000000000e+00 + z * (QS1 + z * (QS2 + z * (QS3 + z * QS4))); 180 s = StrictMath.sqrt(z); 181 r = p / q; 182 w = r * s - PIO2_LO; 183 return 3.14159265358979311600e+00 - 2.0 * (s + w); 184 } else { /* x > 0.5 */ 185 z = (1.00000000000000000000e+00 - x) * 0.5; 186 s = StrictMath.sqrt(z); 187 df = s; 188 df = Double.longBitsToDouble( 189 Double.doubleToRawLongBits(df) & 0xffffffffL << 32); 190 c = (z - df * df) / (s + df); 191 p = z * (PS0 + z 192 * (PS1 + z * (PS2 + z * (PS3 + z * (PS4 + z * PS5))))); 193 q = 1.00000000000000000000e+00 + z * (QS1 + z * (QS2 + z * (QS3 + z * QS4))); 194 r = p / q; 195 w = r * s + c; 196 return 2.0 * (df + w); 197 } 198 } 199 200 /** 201 * Returns the closest double approximation of the arc sine of the argument 202 * within the range {@code [-pi/2..pi/2]}. 203 * <p> 204 * Special cases: 205 * <ul> 206 * <li>{@code asin((anything > 1)) = NaN}</li> 207 * <li>{@code asin((anything < -1)) = NaN}</li> 208 * <li>{@code asin(NaN) = NaN}</li> 209 * </ul> 210 * 211 * @param x 212 * the value whose arc sine has to be computed. 213 * @return the arc sine of the argument. 214 */ 215 public static double asin(double x) { 216 double t, w, p, q, c, r, s; 217 int hx, ix; 218 final long bits = Double.doubleToRawLongBits(x); 219 hx = (int) (bits >>> 32); 220 ix = hx & 0x7fffffff; 221 if (ix >= 0x3ff00000) { /* |x|>= 1 */ 222 if ((((ix - 0x3ff00000) | ((int) bits))) == 0) { 223 /* ieee_asin(1)=+-pi/2 with inexact */ 224 return x * PIO2_HI + x * PIO2_LO; 225 } 226 return (x - x) / (x - x); /* ieee_asin(|x|>1) is NaN */ 227 } else if (ix < 0x3fe00000) { /* |x|<0.5 */ 228 if (ix < 0x3e400000) { /* if |x| < 2**-27 */ 229 if (HUGE + x > 1.00000000000000000000e+00) { 230 return x;/* return x with inexact if x!=0 */ 231 } 232 } else { 233 t = x * x; 234 p = t * (PS0 + t 235 * (PS1 + t * (PS2 + t * (PS3 + t * (PS4 + t * PS5))))); 236 q = 1.00000000000000000000e+00 + t * (QS1 + t * (QS2 + t * (QS3 + t * QS4))); 237 w = p / q; 238 return x + x * w; 239 } 240 } 241 /* 1> |x|>= 0.5 */ 242 w = 1.00000000000000000000e+00 - Math.abs(x); 243 t = w * 0.5; 244 p = t * (PS0 + t * (PS1 + t * (PS2 + t * (PS3 + t * (PS4 + t * PS5))))); 245 q = 1.00000000000000000000e+00 + t * (QS1 + t * (QS2 + t * (QS3 + t * QS4))); 246 s = StrictMath.sqrt(t); 247 if (ix >= 0x3FEF3333) { /* if |x| > 0.975 */ 248 w = p / q; 249 t = PIO2_HI - (2.0 * (s + s * w) - PIO2_LO); 250 } else { 251 w = s; 252 w = Double.longBitsToDouble( 253 Double.doubleToRawLongBits(w) & 0xffffffffL << 32); 254 c = (t - w * w) / (s + w); 255 r = p / q; 256 p = 2.0 * s * r - (PIO2_LO - 2.0 * c); 257 q = PIO4_HI - 2.0 * w; 258 t = PIO4_HI - (p - q); 259 } 260 return (hx > 0) ? t : -t; 261 } 262 263 private static final double[] ATANHI = { 4.63647609000806093515e-01, 264 7.85398163397448278999e-01, 9.82793723247329054082e-01, 265 1.57079632679489655800e+00 }; 266 private static final double[] ATANLO = { 2.26987774529616870924e-17, 267 3.06161699786838301793e-17, 1.39033110312309984516e-17, 268 6.12323399573676603587e-17 }; 269 private static final double AT0 = 3.33333333333329318027e-01; 270 private static final double AT1 = -1.99999999998764832476e-01; 271 private static final double AT2 = 1.42857142725034663711e-01; 272 private static final double AT3 = -1.11111104054623557880e-01; 273 private static final double AT4 = 9.09088713343650656196e-02; 274 private static final double AT5 = -7.69187620504482999495e-02; 275 private static final double AT6 = 6.66107313738753120669e-02; 276 private static final double AT7= -5.83357013379057348645e-02; 277 private static final double AT8 = 4.97687799461593236017e-02; 278 private static final double AT9 = -3.65315727442169155270e-02; 279 private static final double AT10 = 1.62858201153657823623e-02; 280 281 /** 282 * Returns the closest double approximation of the arc tangent of the 283 * argument within the range {@code [-pi/2..pi/2]}. 284 * <p> 285 * Special cases: 286 * <ul> 287 * <li>{@code atan(+0.0) = +0.0}</li> 288 * <li>{@code atan(-0.0) = -0.0}</li> 289 * <li>{@code atan(+infinity) = +pi/2}</li> 290 * <li>{@code atan(-infinity) = -pi/2}</li> 291 * <li>{@code atan(NaN) = NaN}</li> 292 * </ul> 293 * 294 * @param x 295 * the value whose arc tangent has to be computed. 296 * @return the arc tangent of the argument. 297 */ 298 public static double atan(double x) { 299 double w, s1, s2, z; 300 int ix, hx, id; 301 302 final long bits = Double.doubleToRawLongBits(x); 303 hx = (int) (bits >>> 32); 304 ix = hx & 0x7fffffff; 305 if (ix >= 0x44100000) { /* if |x| >= 2^66 */ 306 if (ix > 0x7ff00000 || (ix == 0x7ff00000 && (((int) bits) != 0))) { 307 return x + x; /* NaN */ 308 } 309 if (hx > 0) { 310 return ATANHI[3] + ATANLO[3]; 311 } else { 312 return -ATANHI[3] - ATANLO[3]; 313 } 314 } 315 if (ix < 0x3fdc0000) { /* |x| < 0.4375 */ 316 if (ix < 0x3e200000) { /* |x| < 2^-29 */ 317 if (HUGE + x > 1.00000000000000000000e+00) { 318 return x; /* raise inexact */ 319 } 320 } 321 id = -1; 322 } else { 323 x = Math.abs(x); 324 if (ix < 0x3ff30000) { /* |x| < 1.1875 */ 325 if (ix < 0x3fe60000) { /* 7/16 <=|x|<11/16 */ 326 id = 0; 327 x = (2.0 * x - 1.00000000000000000000e+00) / (2.0 + x); 328 } else { /* 11/16<=|x|< 19/16 */ 329 id = 1; 330 x = (x - 1.00000000000000000000e+00) / (x + 1.00000000000000000000e+00); 331 } 332 } else { 333 if (ix < 0x40038000) { /* |x| < 2.4375 */ 334 id = 2; 335 x = (x - 1.5) / (1.00000000000000000000e+00 + 1.5 * x); 336 } else { /* 2.4375 <= |x| < 2^66 */ 337 id = 3; 338 x = -1.0 / x; 339 } 340 } 341 } 342 343 /* end of argument reduction */ 344 z = x * x; 345 w = z * z; 346 /* break sum from i=0 to 10 aT[i]z**(i+1) into odd and even poly */ 347 s1 = z * (AT0 + w * (AT2 + w 348 * (AT4 + w * (AT6 + w * (AT8 + w * AT10))))); 349 s2 = w * (AT1 + w * (AT3 + w * (AT5 + w * (AT7 + w * AT9)))); 350 if (id < 0) { 351 return x - x * (s1 + s2); 352 } else { 353 z = ATANHI[id] - ((x * (s1 + s2) - ATANLO[id]) - x); 354 return (hx < 0) ? -z : z; 355 } 356 } 357 358 private static final double PI_O_4 = 7.8539816339744827900E-01; 359 private static final double PI_O_2 = 1.5707963267948965580E+00; 360 private static final double PI_LO = 1.2246467991473531772E-16; 361 362 /** 363 * Returns the closest double approximation of the arc tangent of 364 * {@code y/x} within the range {@code [-pi..pi]}. This is the angle of the 365 * polar representation of the rectangular coordinates (x,y). 366 * <p> 367 * Special cases: 368 * <ul> 369 * <li>{@code atan2((anything), NaN ) = NaN;}</li> 370 * <li>{@code atan2(NaN , (anything) ) = NaN;}</li> 371 * <li>{@code atan2(+0.0, +(anything but NaN)) = +0.0}</li> 372 * <li>{@code atan2(-0.0, +(anything but NaN)) = -0.0}</li> 373 * <li>{@code atan2(+0.0, -(anything but NaN)) = +pi}</li> 374 * <li>{@code atan2(-0.0, -(anything but NaN)) = -pi}</li> 375 * <li>{@code atan2(+(anything but 0 and NaN), 0) = +pi/2}</li> 376 * <li>{@code atan2(-(anything but 0 and NaN), 0) = -pi/2}</li> 377 * <li>{@code atan2(+(anything but infinity and NaN), +infinity)} {@code =} 378 * {@code +0.0}</li> 379 * <li>{@code atan2(-(anything but infinity and NaN), +infinity)} {@code =} 380 * {@code -0.0}</li> 381 * <li>{@code atan2(+(anything but infinity and NaN), -infinity) = +pi}</li> 382 * <li>{@code atan2(-(anything but infinity and NaN), -infinity) = -pi}</li> 383 * <li>{@code atan2(+infinity, +infinity ) = +pi/4}</li> 384 * <li>{@code atan2(-infinity, +infinity ) = -pi/4}</li> 385 * <li>{@code atan2(+infinity, -infinity ) = +3pi/4}</li> 386 * <li>{@code atan2(-infinity, -infinity ) = -3pi/4}</li> 387 * <li>{@code atan2(+infinity, (anything but,0, NaN, and infinity))} 388 * {@code =} {@code +pi/2}</li> 389 * <li>{@code atan2(-infinity, (anything but,0, NaN, and infinity))} 390 * {@code =} {@code -pi/2}</li> 391 * </ul> 392 * 393 * @param y 394 * the numerator of the value whose atan has to be computed. 395 * @param x 396 * the denominator of the value whose atan has to be computed. 397 * @return the arc tangent of {@code y/x}. 398 */ 399 public static double atan2(double y, double x) { 400 double z; 401 int k, m, hx, hy, ix, iy; 402 int lx, ly; // watch out, should be unsigned 403 404 final long yBits = Double.doubleToRawLongBits(y); 405 final long xBits = Double.doubleToRawLongBits(x); 406 407 hx = (int) (xBits >>> 32); // __HI(x); 408 ix = hx & 0x7fffffff; 409 lx = (int) xBits; // __LO(x); 410 hy = (int) (yBits >>> 32); // __HI(y); 411 iy = hy & 0x7fffffff; 412 ly = (int) yBits; // __LO(y); 413 if (((ix | ((lx | -lx) >> 31)) > 0x7ff00000) 414 || ((iy | ((ly | -ly) >> 31)) > 0x7ff00000)) { /* x or y is NaN */ 415 return x + y; 416 } 417 if ((hx - 0x3ff00000 | lx) == 0) { 418 return StrictMath.atan(y); /* x=1.0 */ 419 } 420 421 m = ((hy >> 31) & 1) | ((hx >> 30) & 2); /* 2*sign(x)+sign(y) */ 422 423 /* when y = 0 */ 424 if ((iy | ly) == 0) { 425 switch (m) { 426 case 0: 427 case 1: 428 return y; /* ieee_atan(+-0,+anything)=+-0 */ 429 case 2: 430 return 3.14159265358979311600e+00 + TINY;/* ieee_atan(+0,-anything) = pi */ 431 case 3: 432 return -3.14159265358979311600e+00 - TINY;/* ieee_atan(-0,-anything) =-pi */ 433 } 434 } 435 /* when x = 0 */ 436 if ((ix | lx) == 0) 437 return (hy < 0) ? -PI_O_2 - TINY : PI_O_2 + TINY; 438 439 /* when x is INF */ 440 if (ix == 0x7ff00000) { 441 if (iy == 0x7ff00000) { 442 switch (m) { 443 case 0: 444 return PI_O_4 + TINY;/* ieee_atan(+INF,+INF) */ 445 case 1: 446 return -PI_O_4 - TINY;/* ieee_atan(-INF,+INF) */ 447 case 2: 448 return 3.0 * PI_O_4 + TINY;/* ieee_atan(+INF,-INF) */ 449 case 3: 450 return -3.0 * PI_O_4 - TINY;/* ieee_atan(-INF,-INF) */ 451 } 452 } else { 453 switch (m) { 454 case 0: 455 return 0.0; /* ieee_atan(+...,+INF) */ 456 case 1: 457 return -0.0; /* ieee_atan(-...,+INF) */ 458 case 2: 459 return 3.14159265358979311600e+00 + TINY; /* ieee_atan(+...,-INF) */ 460 case 3: 461 return -3.14159265358979311600e+00 - TINY; /* ieee_atan(-...,-INF) */ 462 } 463 } 464 } 465 /* when y is INF */ 466 if (iy == 0x7ff00000) 467 return (hy < 0) ? -PI_O_2 - TINY : PI_O_2 + TINY; 468 469 /* compute y/x */ 470 k = (iy - ix) >> 20; 471 if (k > 60) { 472 z = PI_O_2 + 0.5 * PI_LO; /* |y/x| > 2**60 */ 473 } else if (hx < 0 && k < -60) { 474 z = 0.0; /* |y|/x < -2**60 */ 475 } else { 476 z = StrictMath.atan(Math.abs(y / x)); /* safe to do y/x */ 477 } 478 479 switch (m) { 480 case 0: 481 return z; /* ieee_atan(+,+) */ 482 case 1: 483 // __HI(z) ^= 0x80000000; 484 z = Double.longBitsToDouble( 485 Double.doubleToRawLongBits(z) ^ (0x80000000L << 32)); 486 return z; /* ieee_atan(-,+) */ 487 case 2: 488 return 3.14159265358979311600e+00 - (z - PI_LO);/* ieee_atan(+,-) */ 489 default: /* case 3 */ 490 return (z - PI_LO) - 3.14159265358979311600e+00;/* ieee_atan(-,-) */ 491 } 492 } 493 494 private static final int B1 = 715094163; 495 private static final int B2 = 696219795; 496 private static final double C = 5.42857142857142815906e-01; 497 private static final double D = -7.05306122448979611050e-01; 498 private static final double CBRTE = 1.41428571428571436819e+00; 499 private static final double F = 1.60714285714285720630e+00; 500 private static final double G = 3.57142857142857150787e-01; 501 502 /** 503 * Returns the closest double approximation of the cube root of the 504 * argument. 505 * <p> 506 * Special cases: 507 * <ul> 508 * <li>{@code cbrt(+0.0) = +0.0}</li> 509 * <li>{@code cbrt(-0.0) = -0.0}</li> 510 * <li>{@code cbrt(+infinity) = +infinity}</li> 511 * <li>{@code cbrt(-infinity) = -infinity}</li> 512 * <li>{@code cbrt(NaN) = NaN}</li> 513 * </ul> 514 * 515 * @param x 516 * the value whose cube root has to be computed. 517 * @return the cube root of the argument. 518 */ 519 public static double cbrt(double x) { 520 if (x < 0) { 521 return -cbrt(-x); 522 } 523 int hx; 524 double r, s, w; 525 int sign; // caution: should be unsigned 526 long bits = Double.doubleToRawLongBits(x); 527 528 hx = (int) (bits >>> 32); 529 sign = hx & 0x80000000; /* sign= sign(x) */ 530 hx ^= sign; 531 if (hx >= 0x7ff00000) { 532 return (x + x); /* ieee_cbrt(NaN,INF) is itself */ 533 } 534 535 if ((hx | ((int) bits)) == 0) { 536 return x; /* ieee_cbrt(0) is itself */ 537 } 538 539 // __HI(x) = hx; /* x <- |x| */ 540 bits &= 0x00000000ffffffffL; 541 bits |= ((long) hx << 32); 542 543 long tBits = Double.doubleToRawLongBits(0.0) & 0x00000000ffffffffL; 544 double t = 0.0; 545 /* rough cbrt to 5 bits */ 546 if (hx < 0x00100000) { /* subnormal number */ 547 // __HI(t)=0x43500000; /*set t= 2**54*/ 548 tBits |= 0x43500000L << 32; 549 t = Double.longBitsToDouble(tBits); 550 t *= x; 551 552 // __HI(t)=__HI(t)/3+B2; 553 tBits = Double.doubleToRawLongBits(t); 554 long tBitsHigh = tBits >> 32; 555 tBits &= 0x00000000ffffffffL; 556 tBits |= ((tBitsHigh / 3) + B2) << 32; 557 t = Double.longBitsToDouble(tBits); 558 559 } else { 560 // __HI(t)=hx/3+B1; 561 tBits |= ((long) ((hx / 3) + B1)) << 32; 562 t = Double.longBitsToDouble(tBits); 563 } 564 565 /* new cbrt to 23 bits, may be implemented in single precision */ 566 r = t * t / x; 567 s = C + r * t; 568 t *= G + F / (s + CBRTE + D / s); 569 570 /* chopped to 20 bits and make it larger than ieee_cbrt(x) */ 571 tBits = Double.doubleToRawLongBits(t); 572 tBits &= 0xFFFFFFFFL << 32; 573 tBits += 0x00000001L << 32; 574 t = Double.longBitsToDouble(tBits); 575 576 /* one step newton iteration to 53 bits with error less than 0.667 ulps */ 577 s = t * t; /* t*t is exact */ 578 r = x / s; 579 w = t + t; 580 r = (r - t) / (w + r); /* r-s is exact */ 581 t = t + t * r; 582 583 /* retore the sign bit */ 584 tBits = Double.doubleToRawLongBits(t); 585 tBits |= ((long) sign) << 32; 586 return Double.longBitsToDouble(tBits); 587 } 588 589 /** 590 * Returns the double conversion of the most negative (closest to negative 591 * infinity) integer value greater than or equal to the argument. 592 * <p> 593 * Special cases: 594 * <ul> 595 * <li>{@code ceil(+0.0) = +0.0}</li> 596 * <li>{@code ceil(-0.0) = -0.0}</li> 597 * <li>{@code ceil((anything in range (-1,0)) = -0.0}</li> 598 * <li>{@code ceil(+infinity) = +infinity}</li> 599 * <li>{@code ceil(-infinity) = -infinity}</li> 600 * <li>{@code ceil(NaN) = NaN}</li> 601 * </ul> 602 */ 603 public static native double ceil(double d); 604 605 private static final long ONEBITS = Double.doubleToRawLongBits(1.00000000000000000000e+00) 606 & 0x00000000ffffffffL; 607 608 /** 609 * Returns the closest double approximation of the hyperbolic cosine of the 610 * argument. 611 * <p> 612 * Special cases: 613 * <ul> 614 * <li>{@code cosh(+infinity) = +infinity}</li> 615 * <li>{@code cosh(-infinity) = +infinity}</li> 616 * <li>{@code cosh(NaN) = NaN}</li> 617 * </ul> 618 * 619 * @param x 620 * the value whose hyperbolic cosine has to be computed. 621 * @return the hyperbolic cosine of the argument. 622 */ 623 public static double cosh(double x) { 624 double t, w; 625 int ix; 626 final long bits = Double.doubleToRawLongBits(x); 627 ix = (int) (bits >>> 32) & 0x7fffffff; 628 629 /* x is INF or NaN */ 630 if (ix >= 0x7ff00000) { 631 return x * x; 632 } 633 634 /* |x| in [0,0.5*ln2], return 1+ieee_expm1(|x|)^2/(2*ieee_exp(|x|)) */ 635 if (ix < 0x3fd62e43) { 636 t = expm1(Math.abs(x)); 637 w = 1.00000000000000000000e+00 + t; 638 if (ix < 0x3c800000) 639 return w; /* ieee_cosh(tiny) = 1 */ 640 return 1.00000000000000000000e+00 + (t * t) / (w + w); 641 } 642 643 /* |x| in [0.5*ln2,22], return (ieee_exp(|x|)+1/ieee_exp(|x|)/2; */ 644 if (ix < 0x40360000) { 645 t = exp(Math.abs(x)); 646 return 0.5 * t + 0.5 / t; 647 } 648 649 /* |x| in [22, ieee_log(maxdouble)] return half*ieee_exp(|x|) */ 650 if (ix < 0x40862E42) { 651 return 0.5 * exp(Math.abs(x)); 652 } 653 654 /* |x| in [log(maxdouble), overflowthresold] */ 655 final long lx = ((ONEBITS >>> 29) + ((int) bits)) & 0x00000000ffffffffL; 656 // watch out: lx should be an unsigned int 657 // lx = *( (((*(unsigned*)&one)>>29)) + (unsigned*)&x); 658 if (ix < 0x408633CE || (ix == 0x408633ce) && (lx <= 0x8fb9f87dL)) { 659 w = exp(0.5 * Math.abs(x)); 660 t = 0.5 * w; 661 return t * w; 662 } 663 664 /* |x| > overflowthresold, ieee_cosh(x) overflow */ 665 return HUGE * HUGE; 666 } 667 668 /** 669 * Returns the closest double approximation of the cosine of the argument. 670 * <p> 671 * Special cases: 672 * <ul> 673 * <li>{@code cos(+infinity) = NaN}</li> 674 * <li>{@code cos(-infinity) = NaN}</li> 675 * <li>{@code cos(NaN) = NaN}</li> 676 * </ul> 677 * 678 * @param d 679 * the angle whose cosine has to be computed, in radians. 680 * @return the cosine of the argument. 681 */ 682 public static native double cos(double d); 683 684 private static final double TWON24 = 5.96046447753906250000e-08; 685 private static final double TWO54 = 1.80143985094819840000e+16, 686 TWOM54 = 5.55111512312578270212e-17; 687 private static final double TWOM1000 = 9.33263618503218878990e-302; 688 private static final double O_THRESHOLD = 7.09782712893383973096e+02; 689 private static final double U_THRESHOLD = -7.45133219101941108420e+02; 690 private static final double INVLN2 = 1.44269504088896338700e+00; 691 private static final double P1 = 1.66666666666666019037e-01; 692 private static final double P2 = -2.77777777770155933842e-03; 693 private static final double P3 = 6.61375632143793436117e-05; 694 private static final double P4 = -1.65339022054652515390e-06; 695 private static final double P5 = 4.13813679705723846039e-08; 696 697 /** 698 * Returns the closest double approximation of the raising "e" to the power 699 * of the argument. 700 * <p> 701 * Special cases: 702 * <ul> 703 * <li>{@code exp(+infinity) = +infinity}</li> 704 * <li>{@code exp(-infinity) = +0.0}</li> 705 * <li>{@code exp(NaN) = NaN}</li> 706 * </ul> 707 * 708 * @param x 709 * the value whose exponential has to be computed. 710 * @return the exponential of the argument. 711 */ 712 public static double exp(double x) { 713 double y, c, t; 714 double hi = 0, lo = 0; 715 int k = 0, xsb; 716 int hx; // should be unsigned, be careful! 717 final long bits = Double.doubleToRawLongBits(x); 718 int lowBits = (int) bits; 719 int highBits = (int) (bits >>> 32); 720 hx = highBits & 0x7fffffff; 721 xsb = (highBits >>> 31) & 1; 722 723 /* filter out non-finite argument */ 724 if (hx >= 0x40862E42) { /* if |x|>=709.78... */ 725 if (hx >= 0x7ff00000) { 726 if (((hx & 0xfffff) | lowBits) != 0) { 727 return x + x; /* NaN */ 728 } else { 729 return (xsb == 0) ? x : 0.0; /* ieee_exp(+-inf)={inf,0} */ 730 } 731 } 732 733 if (x > O_THRESHOLD) { 734 return HUGE * HUGE; /* overflow */ 735 } 736 737 if (x < U_THRESHOLD) { 738 return TWOM1000 * TWOM1000; /* underflow */ 739 } 740 } 741 742 /* argument reduction */ 743 if (hx > 0x3fd62e42) { /* if |x| > 0.5 ln2 */ 744 if (hx < 0x3FF0A2B2) { /* and |x| < 1.5 ln2 */ 745 hi = x - ((xsb == 0) ? 6.93147180369123816490e-01 : 746 -6.93147180369123816490e-01); // LN2HI[xsb]; 747 lo = (xsb == 0) ? 1.90821492927058770002e-10 : 748 -1.90821492927058770002e-10; // LN2LO[xsb]; 749 k = 1 - xsb - xsb; 750 } else { 751 k = (int) (INVLN2 * x + ((xsb == 0) ? 0.5 : -0.5 ));//halF[xsb]); 752 t = k; 753 hi = x - t * 6.93147180369123816490e-01; //ln2HI[0]; /* t*ln2HI is exact here */ 754 lo = t * 1.90821492927058770002e-10; //ln2LO[0]; 755 } 756 x = hi - lo; 757 } else if (hx < 0x3e300000) { /* when |x|<2**-28 */ 758 if (HUGE + x > 1.00000000000000000000e+00) 759 return 1.00000000000000000000e+00 + x;/* trigger inexact */ 760 } else { 761 k = 0; 762 } 763 764 /* x is now in primary range */ 765 t = x * x; 766 c = x - t * (P1 + t * (P2 + t * (P3 + t * (P4 + t * P5)))); 767 if (k == 0) { 768 return 1.00000000000000000000e+00 - ((x * c) / (c - 2.0) - x); 769 } else { 770 y = 1.00000000000000000000e+00 - ((lo - (x * c) / (2.0 - c)) - hi); 771 } 772 long yBits = Double.doubleToRawLongBits(y); 773 if (k >= -1021) { 774 yBits += ((long) (k << 20)) << 32; /* add k to y's exponent */ 775 return Double.longBitsToDouble(yBits); 776 } else { 777 yBits += ((long) ((k + 1000) << 20)) << 32;/* add k to y's exponent */ 778 return Double.longBitsToDouble(yBits) * TWOM1000; 779 } 780 } 781 782 private static final double TINY = 1.0e-300; 783 private static final double LN2_HI = 6.93147180369123816490e-01; 784 private static final double LN2_LO = 1.90821492927058770002e-10; 785 private static final double Q1 = -3.33333333333331316428e-02; 786 private static final double Q2 = 1.58730158725481460165e-03; 787 private static final double Q3 = -7.93650757867487942473e-05; 788 private static final double Q4 = 4.00821782732936239552e-06; 789 private static final double Q5 = -2.01099218183624371326e-07; 790 791 /** 792 * Returns the closest double approximation of <i>{@code e}</i><sup> 793 * {@code d}</sup>{@code - 1}. If the argument is very close to 0, it is 794 * much more accurate to use {@code expm1(d)+1} than {@code exp(d)} (due to 795 * cancellation of significant digits). 796 * <p> 797 * Special cases: 798 * <ul> 799 * <li>{@code expm1(+0.0) = +0.0}</li> 800 * <li>{@code expm1(-0.0) = -0.0}</li> 801 * <li>{@code expm1(+infinity) = +infinity}</li> 802 * <li>{@code expm1(-infinity) = -1.0}</li> 803 * <li>{@code expm1(NaN) = NaN}</li> 804 * </ul> 805 * 806 * @param x 807 * the value to compute the <i>{@code e}</i><sup>{@code d}</sup> 808 * {@code - 1} of. 809 * @return the <i>{@code e}</i><sup>{@code d}</sup>{@code - 1} value of the 810 * argument. 811 */ 812 public static double expm1(double x) { 813 double y, hi, lo, t, e, hxs, hfx, r1, c = 0.0; 814 int k, xsb; 815 long yBits = 0; 816 final long bits = Double.doubleToRawLongBits(x); 817 int highBits = (int) (bits >>> 32); 818 int lowBits = (int) (bits); 819 int hx = highBits & 0x7fffffff; // caution: should be unsigned! 820 xsb = highBits & 0x80000000; /* sign bit of x */ 821 y = xsb == 0 ? x : -x; /* y = |x| */ 822 823 /* filter out huge and non-finite argument */ 824 if (hx >= 0x4043687A) { /* if |x|>=56*ln2 */ 825 if (hx >= 0x40862E42) { /* if |x|>=709.78... */ 826 if (hx >= 0x7ff00000) { 827 if (((hx & 0xfffff) | lowBits) != 0) { 828 return x + x; /* NaN */ 829 } else { 830 return (xsb == 0) ? x : -1.0;/* ieee_exp(+-inf)={inf,-1} */ 831 } 832 } 833 if (x > O_THRESHOLD) { 834 return HUGE * HUGE; /* overflow */ 835 } 836 } 837 if (xsb != 0) { /* x < -56*ln2, return -1.0 with inexact */ 838 if (x + TINY < 0.0) { /* raise inexact */ 839 return TINY - 1.00000000000000000000e+00; /* return -1 */ 840 } 841 } 842 } 843 /* argument reduction */ 844 if (hx > 0x3fd62e42) { /* if |x| > 0.5 ln2 */ 845 if (hx < 0x3FF0A2B2) { /* and |x| < 1.5 ln2 */ 846 if (xsb == 0) { 847 hi = x - LN2_HI; 848 lo = LN2_LO; 849 k = 1; 850 } else { 851 hi = x + LN2_HI; 852 lo = -LN2_LO; 853 k = -1; 854 } 855 } else { 856 k = (int) (INVLN2 * x + ((xsb == 0) ? 0.5 : -0.5)); 857 t = k; 858 hi = x - t * LN2_HI; /* t*ln2_hi is exact here */ 859 lo = t * LN2_LO; 860 } 861 x = hi - lo; 862 c = (hi - x) - lo; 863 } else if (hx < 0x3c900000) { /* when |x|<2**-54, return x */ 864 // t = huge+x; /* return x with inexact flags when x!=0 */ 865 // return x - (t-(huge+x)); 866 return x; // inexact flag is not set, but Java ignors this flag 867 // anyway 868 } else { 869 k = 0; 870 } 871 872 /* x is now in primary range */ 873 hfx = 0.5 * x; 874 hxs = x * hfx; 875 r1 = 1.00000000000000000000e+00 + hxs * (Q1 + hxs * (Q2 + hxs * (Q3 + hxs * (Q4 + hxs * Q5)))); 876 t = 3.0 - r1 * hfx; 877 e = hxs * ((r1 - t) / (6.0 - x * t)); 878 if (k == 0) { 879 return x - (x * e - hxs); /* c is 0 */ 880 } else { 881 e = (x * (e - c) - c); 882 e -= hxs; 883 if (k == -1) { 884 return 0.5 * (x - e) - 0.5; 885 } 886 887 if (k == 1) { 888 if (x < -0.25) { 889 return -2.0 * (e - (x + 0.5)); 890 } else { 891 return 1.00000000000000000000e+00 + 2.0 * (x - e); 892 } 893 } 894 895 if (k <= -2 || k > 56) { /* suffice to return ieee_exp(x)-1 */ 896 y = 1.00000000000000000000e+00 - (e - x); 897 yBits = Double.doubleToRawLongBits(y); 898 yBits += (((long) k) << 52); /* add k to y's exponent */ 899 return Double.longBitsToDouble(yBits) - 1.00000000000000000000e+00; 900 } 901 902 long tBits = Double.doubleToRawLongBits(1.00000000000000000000e+00) & 0x00000000ffffffffL; 903 904 if (k < 20) { 905 tBits |= (((long) 0x3ff00000) - (0x200000 >> k)) << 32; 906 y = Double.longBitsToDouble(tBits) - (e - x); 907 yBits = Double.doubleToRawLongBits(y); 908 yBits += (((long) k) << 52); /* add k to y's exponent */ 909 return Double.longBitsToDouble(yBits); 910 } else { 911 tBits |= ((((long) 0x3ff) - k) << 52); /* 2^-k */ 912 y = x - (e + Double.longBitsToDouble(tBits)); 913 y += 1.00000000000000000000e+00; 914 yBits = Double.doubleToRawLongBits(y); 915 yBits += (((long) k) << 52); /* add k to y's exponent */ 916 return Double.longBitsToDouble(yBits); 917 } 918 } 919 } 920 921 /** 922 * Returns the double conversion of the most positive (closest to positive 923 * infinity) integer less than or equal to the argument. 924 * <p> 925 * Special cases: 926 * <ul> 927 * <li>{@code floor(+0.0) = +0.0}</li> 928 * <li>{@code floor(-0.0) = -0.0}</li> 929 * <li>{@code floor(+infinity) = +infinity}</li> 930 * <li>{@code floor(-infinity) = -infinity}</li> 931 * <li>{@code floor(NaN) = NaN}</li> 932 * </ul> 933 */ 934 public static native double floor(double d); 935 936 /** 937 * Returns {@code sqrt(}<i>{@code x}</i><sup>{@code 2}</sup>{@code +} <i> 938 * {@code y}</i><sup>{@code 2}</sup>{@code )}. The final result is without 939 * medium underflow or overflow. 940 * <p> 941 * Special cases: 942 * <ul> 943 * <li>{@code hypot(+infinity, (anything including NaN)) = +infinity}</li> 944 * <li>{@code hypot(-infinity, (anything including NaN)) = +infinity}</li> 945 * <li>{@code hypot((anything including NaN), +infinity) = +infinity}</li> 946 * <li>{@code hypot((anything including NaN), -infinity) = +infinity}</li> 947 * <li>{@code hypot(NaN, NaN) = NaN}</li> 948 * </ul> 949 * 950 * @param x 951 * a double number. 952 * @param y 953 * a double number. 954 * @return the {@code sqrt(}<i>{@code x}</i><sup>{@code 2}</sup>{@code +} 955 * <i> {@code y}</i><sup>{@code 2}</sup>{@code )} value of the 956 * arguments. 957 */ 958 public static native double hypot(double x, double y); 959 960 /** 961 * Returns the remainder of dividing {@code x} by {@code y} using the IEEE 962 * 754 rules. The result is {@code x-round(x/p)*p} where {@code round(x/p)} 963 * is the nearest integer (rounded to even), but without numerical 964 * cancellation problems. 965 * <p> 966 * Special cases: 967 * <ul> 968 * <li>{@code IEEEremainder((anything), 0) = NaN}</li> 969 * <li>{@code IEEEremainder(+infinity, (anything)) = NaN}</li> 970 * <li>{@code IEEEremainder(-infinity, (anything)) = NaN}</li> 971 * <li>{@code IEEEremainder(NaN, (anything)) = NaN}</li> 972 * <li>{@code IEEEremainder((anything), NaN) = NaN}</li> 973 * <li>{@code IEEEremainder(x, +infinity) = x } where x is anything but 974 * +/-infinity</li> 975 * <li>{@code IEEEremainder(x, -infinity) = x } where x is anything but 976 * +/-infinity</li> 977 * </ul> 978 * 979 * @param x 980 * the numerator of the operation. 981 * @param y 982 * the denominator of the operation. 983 * @return the IEEE754 floating point reminder of of {@code x/y}. 984 */ 985 public static native double IEEEremainder(double x, double y); 986 987 private static final double LG1 = 6.666666666666735130e-01; 988 private static final double LG2 = 3.999999999940941908e-01; 989 private static final double LG3 = 2.857142874366239149e-01; 990 private static final double LG4 = 2.222219843214978396e-01; 991 private static final double LG5 = 1.818357216161805012e-01; 992 private static final double LG6 = 1.531383769920937332e-01; 993 private static final double LG7 = 1.479819860511658591e-01; 994 995 /** 996 * Returns the closest double approximation of the natural logarithm of the 997 * argument. 998 * <p> 999 * Special cases: 1000 * <ul> 1001 * <li>{@code log(+0.0) = -infinity}</li> 1002 * <li>{@code log(-0.0) = -infinity}</li> 1003 * <li>{@code log((anything < 0) = NaN}</li> 1004 * <li>{@code log(+infinity) = +infinity}</li> 1005 * <li>{@code log(-infinity) = NaN}</li> 1006 * <li>{@code log(NaN) = NaN}</li> 1007 * </ul> 1008 * 1009 * @param x 1010 * the value whose log has to be computed. 1011 * @return the natural logarithm of the argument. 1012 */ 1013 public static double log(double x) { 1014 double hfsq, f, s, z, R, w, t1, t2, dk; 1015 int hx, i, j, k = 0; 1016 int lx; // watch out, should be unsigned 1017 1018 long bits = Double.doubleToRawLongBits(x); 1019 hx = (int) (bits >>> 32); /* high word of x */ 1020 lx = (int) bits; /* low word of x */ 1021 1022 if (hx < 0x00100000) { /* x < 2**-1022 */ 1023 if (((hx & 0x7fffffff) | lx) == 0) { 1024 return -TWO54 / 0.0; /* ieee_log(+-0)=-inf */ 1025 } 1026 1027 if (hx < 0) { 1028 return (x - x) / 0.0; /* ieee_log(-#) = NaN */ 1029 } 1030 1031 k -= 54; 1032 x *= TWO54; /* subnormal number, scale up x */ 1033 bits = Double.doubleToRawLongBits(x); 1034 hx = (int) (bits >>> 32); /* high word of x */ 1035 } 1036 1037 if (hx >= 0x7ff00000) { 1038 return x + x; 1039 } 1040 1041 k += (hx >> 20) - 1023; 1042 hx &= 0x000fffff; 1043 bits &= 0x00000000ffffffffL; 1044 i = (hx + 0x95f64) & 0x100000; 1045 bits |= ((long) hx | (i ^ 0x3ff00000)) << 32; /* normalize x or x/2 */ 1046 x = Double.longBitsToDouble(bits); 1047 k += (i >> 20); 1048 f = x - 1.0; 1049 1050 if ((0x000fffff & (2 + hx)) < 3) { /* |f| < 2**-20 */ 1051 if (f == 0.0) { 1052 if (k == 0) { 1053 return 0.0; 1054 } else { 1055 dk = k; 1056 } 1057 return dk * LN2_HI + dk * LN2_LO; 1058 } 1059 1060 R = f * f * (0.5 - 0.33333333333333333 * f); 1061 if (k == 0) { 1062 return f - R; 1063 } else { 1064 dk = k; 1065 return dk * LN2_HI - ((R - dk * LN2_LO) - f); 1066 } 1067 } 1068 s = f / (2.0 + f); 1069 dk = k; 1070 z = s * s; 1071 i = hx - 0x6147a; 1072 w = z * z; 1073 j = 0x6b851 - hx; 1074 t1 = w * (LG2 + w * (LG4 + w * LG6)); 1075 t2 = z * (LG1 + w * (LG3 + w * (LG5 + w * LG7))); 1076 i |= j; 1077 R = t2 + t1; 1078 if (i > 0) { 1079 hfsq = 0.5 * f * f; 1080 if (k == 0) { 1081 return f - (hfsq - s * (hfsq + R)); 1082 } else { 1083 return dk * LN2_HI 1084 - ((hfsq - (s * (hfsq + R) + dk * LN2_LO)) - f); 1085 } 1086 } else { 1087 if (k == 0) { 1088 return f - s * (f - R); 1089 } else { 1090 return dk * LN2_HI - ((s * (f - R) - dk * LN2_LO) - f); 1091 } 1092 } 1093 } 1094 1095 private static final double IVLN10 = 4.34294481903251816668e-01; 1096 private static final double LOG10_2HI = 3.01029995663611771306e-01; 1097 private static final double LOG10_2LO = 3.69423907715893078616e-13; 1098 1099 /** 1100 * Returns the closest double approximation of the base 10 logarithm of the 1101 * argument. 1102 * <p> 1103 * Special cases: 1104 * <ul> 1105 * <li>{@code log10(+0.0) = -infinity}</li> 1106 * <li>{@code log10(-0.0) = -infinity}</li> 1107 * <li>{@code log10((anything < 0) = NaN}</li> 1108 * <li>{@code log10(+infinity) = +infinity}</li> 1109 * <li>{@code log10(-infinity) = NaN}</li> 1110 * <li>{@code log10(NaN) = NaN}</li> 1111 * </ul> 1112 * 1113 * @param x 1114 * the value whose base 10 log has to be computed. 1115 * @return the the base 10 logarithm of x 1116 */ 1117 public static double log10(double x) { 1118 double y, z; 1119 int i, k = 0, hx; 1120 int lx; // careful: lx should be unsigned! 1121 long bits = Double.doubleToRawLongBits(x); 1122 hx = (int) (bits >> 32); /* high word of x */ 1123 lx = (int) bits; /* low word of x */ 1124 if (hx < 0x00100000) { /* x < 2**-1022 */ 1125 if (((hx & 0x7fffffff) | lx) == 0) { 1126 return -TWO54 / 0.0; /* ieee_log(+-0)=-inf */ 1127 } 1128 1129 if (hx < 0) { 1130 return (x - x) / 0.0; /* ieee_log(-#) = NaN */ 1131 } 1132 1133 k -= 54; 1134 x *= TWO54; /* subnormal number, scale up x */ 1135 bits = Double.doubleToRawLongBits(x); 1136 hx = (int) (bits >> 32); /* high word of x */ 1137 } 1138 1139 if (hx >= 0x7ff00000) { 1140 return x + x; 1141 } 1142 1143 k += (hx >> 20) - 1023; 1144 i = (int) (((k & 0x00000000ffffffffL) & 0x80000000) >>> 31); 1145 hx = (hx & 0x000fffff) | ((0x3ff - i) << 20); 1146 y = k + i; 1147 bits &= 0x00000000ffffffffL; 1148 bits |= ((long) hx) << 32; 1149 x = Double.longBitsToDouble(bits); // __HI(x) = hx; 1150 z = y * LOG10_2LO + IVLN10 * log(x); 1151 return z + y * LOG10_2HI; 1152 } 1153 1154 private static final double LP1 = 6.666666666666735130e-01, 1155 LP2 = 3.999999999940941908e-01, /* 3FD99999 9997FA04 */ 1156 LP3 = 2.857142874366239149e-01, /* 3FD24924 94229359 */ 1157 LP4 = 2.222219843214978396e-01, /* 3FCC71C5 1D8E78AF */ 1158 LP5 = 1.818357216161805012e-01, /* 3FC74664 96CB03DE */ 1159 LP6 = 1.531383769920937332e-01, /* 3FC39A09 D078C69F */ 1160 LP7 = 1.479819860511658591e-01; /* 3FC2F112 DF3E5244 */ 1161 1162 /** 1163 * Returns the closest double approximation of the natural logarithm of the 1164 * sum of the argument and 1. If the argument is very close to 0, it is much 1165 * more accurate to use {@code log1p(d)} than {@code log(1.0+d)} (due to 1166 * numerical cancellation). 1167 * <p> 1168 * Special cases: 1169 * <ul> 1170 * <li>{@code log1p(+0.0) = +0.0}</li> 1171 * <li>{@code log1p(-0.0) = -0.0}</li> 1172 * <li>{@code log1p((anything < 1)) = NaN}</li> 1173 * <li>{@code log1p(-1.0) = -infinity}</li> 1174 * <li>{@code log1p(+infinity) = +infinity}</li> 1175 * <li>{@code log1p(-infinity) = NaN}</li> 1176 * <li>{@code log1p(NaN) = NaN}</li> 1177 * </ul> 1178 * 1179 * @param x 1180 * the value to compute the {@code ln(1+d)} of. 1181 * @return the natural logarithm of the sum of the argument and 1. 1182 */ 1183 1184 public static double log1p(double x) { 1185 double hfsq, f = 0.0, c = 0.0, s, z, R, u = 0.0; 1186 int k, hx, hu = 0, ax; 1187 1188 final long bits = Double.doubleToRawLongBits(x); 1189 hx = (int) (bits >>> 32); /* high word of x */ 1190 ax = hx & 0x7fffffff; 1191 1192 k = 1; 1193 if (hx < 0x3FDA827A) { /* x < 0.41422 */ 1194 if (ax >= 0x3ff00000) { /* x <= -1.0 */ 1195 if (x == -1.0) { 1196 return -TWO54 / 0.0; /* ieee_log1p(-1)=+inf */ 1197 } else { 1198 return (x - x) / (x - x); /* ieee_log1p(x<-1)=NaN */ 1199 } 1200 } 1201 if (ax < 0x3e200000) { 1202 if (TWO54 + x > 0.0 && ax < 0x3c900000) { 1203 return x; 1204 } else { 1205 return x - x * x * 0.5; 1206 } 1207 } 1208 if (hx > 0 || hx <= 0xbfd2bec3) { 1209 k = 0; 1210 f = x; 1211 hu = 1; 1212 } /* -0.2929<x<0.41422 */ 1213 } 1214 1215 if (hx >= 0x7ff00000) { 1216 return x + x; 1217 } 1218 1219 if (k != 0) { 1220 long uBits; 1221 if (hx < 0x43400000) { 1222 u = 1.0 + x; 1223 uBits = Double.doubleToRawLongBits(u); 1224 hu = (int) (uBits >>> 32); 1225 k = (hu >> 20) - 1023; 1226 c = (k > 0) ? 1.0 - (u - x) : x - (u - 1.0);/* correction term */ 1227 c /= u; 1228 } else { 1229 uBits = Double.doubleToRawLongBits(x); 1230 hu = (int) (uBits >>> 32); 1231 k = (hu >> 20) - 1023; 1232 c = 0; 1233 } 1234 hu &= 0x000fffff; 1235 if (hu < 0x6a09e) { 1236 // __HI(u) = hu|0x3ff00000; /* normalize u */ 1237 uBits &= 0x00000000ffffffffL; 1238 uBits |= ((long) hu | 0x3ff00000) << 32; 1239 u = Double.longBitsToDouble(uBits); 1240 } else { 1241 k += 1; 1242 // __HI(u) = hu|0x3fe00000; /* normalize u/2 */ 1243 uBits &= 0xffffffffL; 1244 uBits |= ((long) hu | 0x3fe00000) << 32; 1245 u = Double.longBitsToDouble(uBits); 1246 hu = (0x00100000 - hu) >> 2; 1247 } 1248 f = u - 1.0; 1249 } 1250 hfsq = 0.5 * f * f; 1251 if (hu == 0) { /* |f| < 2**-20 */ 1252 if (f == 0.0) { 1253 if (k == 0) { 1254 return 0.0; 1255 } else { 1256 c += k * LN2_LO; 1257 return k * LN2_HI + c; 1258 } 1259 } 1260 1261 R = hfsq * (1.0 - 0.66666666666666666 * f); 1262 if (k == 0) { 1263 return f - R; 1264 } else { 1265 return k * LN2_HI - ((R - (k * LN2_LO + c)) - f); 1266 } 1267 } 1268 1269 s = f / (2.0 + f); 1270 z = s * s; 1271 R = z * (LP1 + z * (LP2 + z 1272 * (LP3 + z * (LP4 + z * (LP5 + z * (LP6 + z * LP7)))))); 1273 if (k == 0) { 1274 return f - (hfsq - s * (hfsq + R)); 1275 } else { 1276 return k * LN2_HI 1277 - ((hfsq - (s * (hfsq + R) + (k * LN2_LO + c))) - f); 1278 } 1279 } 1280 1281 /** 1282 * Returns the most positive (closest to positive infinity) of the two 1283 * arguments. 1284 * <p> 1285 * Special cases: 1286 * <ul> 1287 * <li>{@code max(NaN, (anything)) = NaN}</li> 1288 * <li>{@code max((anything), NaN) = NaN}</li> 1289 * <li>{@code max(+0.0, -0.0) = +0.0}</li> 1290 * <li>{@code max(-0.0, +0.0) = +0.0}</li> 1291 * </ul> 1292 */ 1293 public static double max(double d1, double d2) { 1294 if (d1 > d2) 1295 return d1; 1296 if (d1 < d2) 1297 return d2; 1298 /* if either arg is NaN, return NaN */ 1299 if (d1 != d2) 1300 return Double.NaN; 1301 /* max( +0.0,-0.0) == +0.0 */ 1302 if (d1 == 0.0 && 1303 ((Double.doubleToLongBits(d1) & Double.doubleToLongBits(d2)) & 0x8000000000000000L) == 0) 1304 return 0.0; 1305 return d1; 1306 } 1307 1308 /** 1309 * Returns the most positive (closest to positive infinity) of the two 1310 * arguments. 1311 * <p> 1312 * Special cases: 1313 * <ul> 1314 * <li>{@code max(NaN, (anything)) = NaN}</li> 1315 * <li>{@code max((anything), NaN) = NaN}</li> 1316 * <li>{@code max(+0.0, -0.0) = +0.0}</li> 1317 * <li>{@code max(-0.0, +0.0) = +0.0}</li> 1318 * </ul> 1319 */ 1320 public static float max(float f1, float f2) { 1321 if (f1 > f2) 1322 return f1; 1323 if (f1 < f2) 1324 return f2; 1325 /* if either arg is NaN, return NaN */ 1326 if (f1 != f2) 1327 return Float.NaN; 1328 /* max( +0.0,-0.0) == +0.0 */ 1329 if (f1 == 0.0f && 1330 ((Float.floatToIntBits(f1) & Float.floatToIntBits(f2)) & 0x80000000) == 0) 1331 return 0.0f; 1332 return f1; 1333 } 1334 1335 /** 1336 * Returns the most positive (closest to positive infinity) of the two 1337 * arguments. 1338 */ 1339 public static int max(int i1, int i2) { 1340 return Math.max(i1, i2); 1341 } 1342 1343 /** 1344 * Returns the most positive (closest to positive infinity) of the two 1345 * arguments. 1346 */ 1347 public static long max(long l1, long l2) { 1348 return l1 > l2 ? l1 : l2; 1349 } 1350 1351 /** 1352 * Returns the most negative (closest to negative infinity) of the two 1353 * arguments. 1354 * <p> 1355 * Special cases: 1356 * <ul> 1357 * <li>{@code min(NaN, (anything)) = NaN}</li> 1358 * <li>{@code min((anything), NaN) = NaN}</li> 1359 * <li>{@code min(+0.0, -0.0) = -0.0}</li> 1360 * <li>{@code min(-0.0, +0.0) = -0.0}</li> 1361 * </ul> 1362 */ 1363 public static double min(double d1, double d2) { 1364 if (d1 > d2) 1365 return d2; 1366 if (d1 < d2) 1367 return d1; 1368 /* if either arg is NaN, return NaN */ 1369 if (d1 != d2) 1370 return Double.NaN; 1371 /* min( +0.0,-0.0) == -0.0 */ 1372 if (d1 == 0.0 && 1373 ((Double.doubleToLongBits(d1) | Double.doubleToLongBits(d2)) & 0x8000000000000000l) != 0) 1374 return 0.0 * (-1.0); 1375 return d1; 1376 } 1377 1378 /** 1379 * Returns the most negative (closest to negative infinity) of the two 1380 * arguments. 1381 * <p> 1382 * Special cases: 1383 * <ul> 1384 * <li>{@code min(NaN, (anything)) = NaN}</li> 1385 * <li>{@code min((anything), NaN) = NaN}</li> 1386 * <li>{@code min(+0.0, -0.0) = -0.0}</li> 1387 * <li>{@code min(-0.0, +0.0) = -0.0}</li> 1388 * </ul> 1389 */ 1390 public static float min(float f1, float f2) { 1391 if (f1 > f2) 1392 return f2; 1393 if (f1 < f2) 1394 return f1; 1395 /* if either arg is NaN, return NaN */ 1396 if (f1 != f2) 1397 return Float.NaN; 1398 /* min( +0.0,-0.0) == -0.0 */ 1399 if (f1 == 0.0f && 1400 ((Float.floatToIntBits(f1) | Float.floatToIntBits(f2)) & 0x80000000) != 0) 1401 return 0.0f * (-1.0f); 1402 return f1; 1403 } 1404 1405 /** 1406 * Returns the most negative (closest to negative infinity) of the two 1407 * arguments. 1408 */ 1409 public static int min(int i1, int i2) { 1410 return Math.min(i1, i2); 1411 } 1412 1413 /** 1414 * Returns the most negative (closest to negative infinity) of the two 1415 * arguments. 1416 */ 1417 public static long min(long l1, long l2) { 1418 return l1 < l2 ? l1 : l2; 1419 } 1420 1421 /** 1422 * Returns the closest double approximation of the result of raising 1423 * {@code x} to the power of {@code y}. 1424 * <p> 1425 * Special cases: 1426 * <ul> 1427 * <li>{@code pow((anything), +0.0) = 1.0}</li> 1428 * <li>{@code pow((anything), -0.0) = 1.0}</li> 1429 * <li>{@code pow(x, 1.0) = x}</li> 1430 * <li>{@code pow((anything), NaN) = NaN}</li> 1431 * <li>{@code pow(NaN, (anything except 0)) = NaN}</li> 1432 * <li>{@code pow(+/-(|x| > 1), +infinity) = +infinity}</li> 1433 * <li>{@code pow(+/-(|x| > 1), -infinity) = +0.0}</li> 1434 * <li>{@code pow(+/-(|x| < 1), +infinity) = +0.0}</li> 1435 * <li>{@code pow(+/-(|x| < 1), -infinity) = +infinity}</li> 1436 * <li>{@code pow(+/-1.0 , +infinity) = NaN}</li> 1437 * <li>{@code pow(+/-1.0 , -infinity) = NaN}</li> 1438 * <li>{@code pow(+0.0, (+anything except 0, NaN)) = +0.0}</li> 1439 * <li>{@code pow(-0.0, (+anything except 0, NaN, odd integer)) = +0.0}</li> 1440 * <li>{@code pow(+0.0, (-anything except 0, NaN)) = +infinity}</li> 1441 * <li>{@code pow(-0.0, (-anything except 0, NAN, odd integer))} {@code =} 1442 * {@code +infinity}</li> 1443 * <li>{@code pow(-0.0, (odd integer)) = -pow( +0 , (odd integer) )}</li> 1444 * <li>{@code pow(+infinity, (+anything except 0, NaN)) = +infinity}</li> 1445 * <li>{@code pow(+infinity, (-anything except 0, NaN)) = +0.0}</li> 1446 * <li>{@code pow(-infinity, (anything)) = -pow(0, (-anything))}</li> 1447 * <li>{@code pow((-anything), (integer))} {@code =} 1448 * {@code pow(-1,(integer))*pow(+anything,integer)}</li> 1449 * <li>{@code pow((-anything except 0 and infinity), (non-integer))} 1450 * {@code =} {@code NAN}</li> 1451 * </ul> 1452 * 1453 * @param x 1454 * the base of the operation. 1455 * @param y 1456 * the exponent of the operation. 1457 * @return {@code x} to the power of {@code y}. 1458 */ 1459 public static native double pow(double x, double y); 1460 1461 /** 1462 * Returns a pseudo-random number between 0.0 (inclusive) and 1.0 1463 * (exclusive). 1464 * 1465 * @return a pseudo-random number. 1466 */ 1467 public static double random() { 1468 return Math.random(); 1469 } 1470 1471 /** 1472 * Returns the double conversion of the result of rounding the argument to 1473 * an integer. Tie breaks are rounded towards even. 1474 * <p> 1475 * Special cases: 1476 * <ul> 1477 * <li>{@code rint(+0.0) = +0.0}</li> 1478 * <li>{@code rint(-0.0) = -0.0}</li> 1479 * <li>{@code rint(+infinity) = +infinity}</li> 1480 * <li>{@code rint(-infinity) = -infinity}</li> 1481 * <li>{@code rint(NaN) = NaN}</li> 1482 * </ul> 1483 * 1484 * @param d 1485 * the value to be rounded. 1486 * @return the closest integer to the argument (as a double). 1487 */ 1488 public static native double rint(double d); 1489 1490 /** 1491 * Returns the result of rounding the argument to an integer. The result is 1492 * equivalent to {@code (long) Math.floor(d+0.5)}. 1493 * <p> 1494 * Special cases: 1495 * <ul> 1496 * <li>{@code round(+0.0) = +0.0}</li> 1497 * <li>{@code round(-0.0) = +0.0}</li> 1498 * <li>{@code round((anything > Long.MAX_VALUE) = Long.MAX_VALUE}</li> 1499 * <li>{@code round((anything < Long.MIN_VALUE) = Long.MIN_VALUE}</li> 1500 * <li>{@code round(+infinity) = Long.MAX_VALUE}</li> 1501 * <li>{@code round(-infinity) = Long.MIN_VALUE}</li> 1502 * <li>{@code round(NaN) = +0.0}</li> 1503 * </ul> 1504 * 1505 * @param d 1506 * the value to be rounded. 1507 * @return the closest integer to the argument. 1508 */ 1509 public static long round(double d) { 1510 return Math.round(d); 1511 } 1512 1513 /** 1514 * Returns the result of rounding the argument to an integer. The result is 1515 * equivalent to {@code (int) Math.floor(f+0.5)}. 1516 * <p> 1517 * Special cases: 1518 * <ul> 1519 * <li>{@code round(+0.0) = +0.0}</li> 1520 * <li>{@code round(-0.0) = +0.0}</li> 1521 * <li>{@code round((anything > Integer.MAX_VALUE) = Integer.MAX_VALUE}</li> 1522 * <li>{@code round((anything < Integer.MIN_VALUE) = Integer.MIN_VALUE}</li> 1523 * <li>{@code round(+infinity) = Integer.MAX_VALUE}</li> 1524 * <li>{@code round(-infinity) = Integer.MIN_VALUE}</li> 1525 * <li>{@code round(NaN) = +0.0}</li> 1526 * </ul> 1527 * 1528 * @param f 1529 * the value to be rounded. 1530 * @return the closest integer to the argument. 1531 */ 1532 public static int round(float f) { 1533 return Math.round(f); 1534 } 1535 1536 /** 1537 * Returns the signum function of the argument. If the argument is less than 1538 * zero, it returns -1.0. If the argument is greater than zero, 1.0 is 1539 * returned. If the argument is either positive or negative zero, the 1540 * argument is returned as result. 1541 * <p> 1542 * Special cases: 1543 * <ul> 1544 * <li>{@code signum(+0.0) = +0.0}</li> 1545 * <li>{@code signum(-0.0) = -0.0}</li> 1546 * <li>{@code signum(+infinity) = +1.0}</li> 1547 * <li>{@code signum(-infinity) = -1.0}</li> 1548 * <li>{@code signum(NaN) = NaN}</li> 1549 * </ul> 1550 * 1551 * @param d 1552 * the value whose signum has to be computed. 1553 * @return the value of the signum function. 1554 */ 1555 public static double signum(double d) { 1556 return Math.signum(d); 1557 } 1558 1559 /** 1560 * Returns the signum function of the argument. If the argument is less than 1561 * zero, it returns -1.0. If the argument is greater than zero, 1.0 is 1562 * returned. If the argument is either positive or negative zero, the 1563 * argument is returned as result. 1564 * <p> 1565 * Special cases: 1566 * <ul> 1567 * <li>{@code signum(+0.0) = +0.0}</li> 1568 * <li>{@code signum(-0.0) = -0.0}</li> 1569 * <li>{@code signum(+infinity) = +1.0}</li> 1570 * <li>{@code signum(-infinity) = -1.0}</li> 1571 * <li>{@code signum(NaN) = NaN}</li> 1572 * </ul> 1573 * 1574 * @param f 1575 * the value whose signum has to be computed. 1576 * @return the value of the signum function. 1577 */ 1578 public static float signum(float f) { 1579 return Math.signum(f); 1580 } 1581 1582 private static final double shuge = 1.0e307; 1583 1584 /** 1585 * Returns the closest double approximation of the hyperbolic sine of the 1586 * argument. 1587 * <p> 1588 * Special cases: 1589 * <ul> 1590 * <li>{@code sinh(+0.0) = +0.0}</li> 1591 * <li>{@code sinh(-0.0) = -0.0}</li> 1592 * <li>{@code sinh(+infinity) = +infinity}</li> 1593 * <li>{@code sinh(-infinity) = -infinity}</li> 1594 * <li>{@code sinh(NaN) = NaN}</li> 1595 * </ul> 1596 * 1597 * @param x 1598 * the value whose hyperbolic sine has to be computed. 1599 * @return the hyperbolic sine of the argument. 1600 */ 1601 public static double sinh(double x) { 1602 double t, w, h; 1603 int ix, jx; 1604 final long bits = Double.doubleToRawLongBits(x); 1605 1606 jx = (int) (bits >>> 32); 1607 ix = jx & 0x7fffffff; 1608 1609 /* x is INF or NaN */ 1610 if (ix >= 0x7ff00000) { 1611 return x + x; 1612 } 1613 1614 h = 0.5; 1615 if (jx < 0) { 1616 h = -h; 1617 } 1618 1619 /* |x| in [0,22], return sign(x)*0.5*(E+E/(E+1))) */ 1620 if (ix < 0x40360000) { /* |x|<22 */ 1621 if (ix < 0x3e300000) /* |x|<2**-28 */ 1622 if (shuge + x > 1.00000000000000000000e+00) { 1623 return x;/* ieee_sinh(tiny) = tiny with inexact */ 1624 } 1625 t = expm1(Math.abs(x)); 1626 if (ix < 0x3ff00000) 1627 return h * (2.0 * t - t * t / (t + 1.00000000000000000000e+00)); 1628 return h * (t + t / (t + 1.00000000000000000000e+00)); 1629 } 1630 1631 /* |x| in [22, ieee_log(maxdouble)] return 0.5*ieee_exp(|x|) */ 1632 if (ix < 0x40862E42) { 1633 return h * exp(Math.abs(x)); 1634 } 1635 1636 /* |x| in [log(maxdouble), overflowthresold] */ 1637 final long lx = ((ONEBITS >>> 29) + ((int) bits)) & 0x00000000ffffffffL; 1638 // lx = *( (((*(unsigned*)&one)>>29)) + (unsigned*)&x); 1639 if (ix < 0x408633CE || (ix == 0x408633ce) && (lx <= 0x8fb9f87dL)) { 1640 w = exp(0.5 * Math.abs(x)); 1641 t = h * w; 1642 return t * w; 1643 } 1644 1645 /* |x| > overflowthresold, ieee_sinh(x) overflow */ 1646 return x * shuge; 1647 } 1648 1649 /** 1650 * Returns the closest double approximation of the sine of the argument. 1651 * <p> 1652 * Special cases: 1653 * <ul> 1654 * <li>{@code sin(+0.0) = +0.0}</li> 1655 * <li>{@code sin(-0.0) = -0.0}</li> 1656 * <li>{@code sin(+infinity) = NaN}</li> 1657 * <li>{@code sin(-infinity) = NaN}</li> 1658 * <li>{@code sin(NaN) = NaN}</li> 1659 * </ul> 1660 * 1661 * @param d 1662 * the angle whose sin has to be computed, in radians. 1663 * @return the sine of the argument. 1664 */ 1665 public static native double sin(double d); 1666 1667 /** 1668 * Returns the closest double approximation of the square root of the 1669 * argument. 1670 * <p> 1671 * Special cases: 1672 * <ul> 1673 * <li>{@code sqrt(+0.0) = +0.0}</li> 1674 * <li>{@code sqrt(-0.0) = -0.0}</li> 1675 * <li>{@code sqrt( (anything < 0) ) = NaN}</li> 1676 * <li>{@code sqrt(+infinity) = +infinity}</li> 1677 * <li>{@code sqrt(NaN) = NaN}</li> 1678 * </ul> 1679 */ 1680 public static native double sqrt(double d); 1681 1682 /** 1683 * Returns the closest double approximation of the tangent of the argument. 1684 * <p> 1685 * Special cases: 1686 * <ul> 1687 * <li>{@code tan(+0.0) = +0.0}</li> 1688 * <li>{@code tan(-0.0) = -0.0}</li> 1689 * <li>{@code tan(+infinity) = NaN}</li> 1690 * <li>{@code tan(-infinity) = NaN}</li> 1691 * <li>{@code tan(NaN) = NaN}</li> 1692 * </ul> 1693 * 1694 * @param d 1695 * the angle whose tangent has to be computed, in radians. 1696 * @return the tangent of the argument. 1697 */ 1698 public static native double tan(double d); 1699 1700 /** 1701 * Returns the closest double approximation of the hyperbolic tangent of the 1702 * argument. The absolute value is always less than 1. 1703 * <p> 1704 * Special cases: 1705 * <ul> 1706 * <li>{@code tanh(+0.0) = +0.0}</li> 1707 * <li>{@code tanh(-0.0) = -0.0}</li> 1708 * <li>{@code tanh(+infinity) = +1.0}</li> 1709 * <li>{@code tanh(-infinity) = -1.0}</li> 1710 * <li>{@code tanh(NaN) = NaN}</li> 1711 * </ul> 1712 * 1713 * @param x 1714 * the value whose hyperbolic tangent has to be computed. 1715 * @return the hyperbolic tangent of the argument 1716 */ 1717 public static double tanh(double x) { 1718 double t, z; 1719 int jx, ix; 1720 1721 final long bits = Double.doubleToRawLongBits(x); 1722 /* High word of |x|. */ 1723 jx = (int) (bits >>> 32); 1724 ix = jx & 0x7fffffff; 1725 1726 /* x is INF or NaN */ 1727 if (ix >= 0x7ff00000) { 1728 if (jx >= 0) { 1729 return 1.00000000000000000000e+00 / x + 1.00000000000000000000e+00; /* ieee_tanh(+-inf)=+-1 */ 1730 } else { 1731 return 1.00000000000000000000e+00 / x - 1.00000000000000000000e+00; /* ieee_tanh(NaN) = NaN */ 1732 } 1733 } 1734 1735 /* |x| < 22 */ 1736 if (ix < 0x40360000) { /* |x|<22 */ 1737 if (ix < 0x3c800000) { /* |x|<2**-55 */ 1738 return x * (1.00000000000000000000e+00 + x);/* ieee_tanh(small) = small */ 1739 } 1740 1741 if (ix >= 0x3ff00000) { /* |x|>=1 */ 1742 t = Math.expm1(2.0 * Math.abs(x)); 1743 z = 1.00000000000000000000e+00 - 2.0 / (t + 2.0); 1744 } else { 1745 t = Math.expm1(-2.0 * Math.abs(x)); 1746 z = -t / (t + 2.0); 1747 } 1748 /* |x| > 22, return +-1 */ 1749 } else { 1750 z = 1.00000000000000000000e+00 - TINY; /* raised inexact flag */ 1751 } 1752 return (jx >= 0) ? z : -z; 1753 } 1754 1755 /** 1756 * Returns the measure in degrees of the supplied radian angle. The result 1757 * is {@code angrad * 180 / pi}. 1758 * <p> 1759 * Special cases: 1760 * <ul> 1761 * <li>{@code toDegrees(+0.0) = +0.0}</li> 1762 * <li>{@code toDegrees(-0.0) = -0.0}</li> 1763 * <li>{@code toDegrees(+infinity) = +infinity}</li> 1764 * <li>{@code toDegrees(-infinity) = -infinity}</li> 1765 * <li>{@code toDegrees(NaN) = NaN}</li> 1766 * </ul> 1767 * 1768 * @param angrad 1769 * an angle in radians. 1770 * @return the degree measure of the angle. 1771 */ 1772 public static double toDegrees(double angrad) { 1773 return Math.toDegrees(angrad); 1774 } 1775 1776 /** 1777 * Returns the measure in radians of the supplied degree angle. The result 1778 * is {@code angdeg / 180 * pi}. 1779 * <p> 1780 * Special cases: 1781 * <ul> 1782 * <li>{@code toRadians(+0.0) = +0.0}</li> 1783 * <li>{@code toRadians(-0.0) = -0.0}</li> 1784 * <li>{@code toRadians(+infinity) = +infinity}</li> 1785 * <li>{@code toRadians(-infinity) = -infinity}</li> 1786 * <li>{@code toRadians(NaN) = NaN}</li> 1787 * </ul> 1788 * 1789 * @param angdeg 1790 * an angle in degrees. 1791 * @return the radian measure of the angle. 1792 */ 1793 public static double toRadians(double angdeg) { 1794 return Math.toRadians(angdeg); 1795 } 1796 1797 /** 1798 * Returns the argument's ulp (unit in the last place). The size of a ulp of 1799 * a double value is the positive distance between this value and the double 1800 * value next larger in magnitude. For non-NaN {@code x}, 1801 * {@code ulp(-x) == ulp(x)}. 1802 * <p> 1803 * Special cases: 1804 * <ul> 1805 * <li>{@code ulp(+0.0) = Double.MIN_VALUE}</li> 1806 * <li>{@code ulp(-0.0) = Double.MIN_VALUE}</li> 1807 * <li>{@code ulp(+infinity) = infinity}</li> 1808 * <li>{@code ulp(-infinity) = infinity}</li> 1809 * <li>{@code ulp(NaN) = NaN}</li> 1810 * </ul> 1811 * 1812 * @param d 1813 * the floating-point value to compute ulp of. 1814 * @return the size of a ulp of the argument. 1815 */ 1816 public static double ulp(double d) { 1817 // special cases 1818 if (Double.isInfinite(d)) { 1819 return Double.POSITIVE_INFINITY; 1820 } else if (d == Double.MAX_VALUE || d == -Double.MAX_VALUE) { 1821 return pow(2, 971); 1822 } 1823 d = Math.abs(d); 1824 return nextafter(d, Double.MAX_VALUE) - d; 1825 } 1826 1827 /** 1828 * Returns the argument's ulp (unit in the last place). The size of a ulp of 1829 * a float value is the positive distance between this value and the float 1830 * value next larger in magnitude. For non-NaN {@code x}, 1831 * {@code ulp(-x) == ulp(x)}. 1832 * <p> 1833 * Special cases: 1834 * <ul> 1835 * <li>{@code ulp(+0.0) = Float.MIN_VALUE}</li> 1836 * <li>{@code ulp(-0.0) = Float.MIN_VALUE}</li> 1837 * <li>{@code ulp(+infinity) = infinity}</li> 1838 * <li>{@code ulp(-infinity) = infinity}</li> 1839 * <li>{@code ulp(NaN) = NaN}</li> 1840 * </ul> 1841 * 1842 * @param f 1843 * the floating-point value to compute ulp of. 1844 * @return the size of a ulp of the argument. 1845 */ 1846 public static float ulp(float f) { 1847 return Math.ulp(f); 1848 } 1849 1850 private static native double nextafter(double x, double y); 1851 1852 /** 1853 * Returns a double with the given magnitude and the sign of {@code sign}. 1854 * If {@code sign} is NaN, the sign of the result is positive. 1855 * 1856 * @since 1.6 1857 */ 1858 public static double copySign(double magnitude, double sign) { 1859 // We manually inline Double.isNaN here because the JIT can't do it yet. 1860 // With Double.isNaN: 236.3ns 1861 // With manual inline: 141.2ns 1862 // With no check (i.e. Math's behavior): 110.0ns 1863 // (Tested on a Nexus One.) 1864 long magnitudeBits = Double.doubleToRawLongBits(magnitude); 1865 long signBits = Double.doubleToRawLongBits((sign != sign) ? 1.0 : sign); 1866 magnitudeBits = (magnitudeBits & ~Double.SIGN_MASK) 1867 | (signBits & Double.SIGN_MASK); 1868 return Double.longBitsToDouble(magnitudeBits); 1869 } 1870 1871 /** 1872 * Returns a float with the given magnitude and the sign of {@code sign}. If 1873 * {@code sign} is NaN, the sign of the result is positive. 1874 * 1875 * @since 1.6 1876 */ 1877 public static float copySign(float magnitude, float sign) { 1878 // We manually inline Float.isNaN here because the JIT can't do it yet. 1879 // With Float.isNaN: 214.7ns 1880 // With manual inline: 112.3ns 1881 // With no check (i.e. Math's behavior): 93.1ns 1882 // (Tested on a Nexus One.) 1883 int magnitudeBits = Float.floatToRawIntBits(magnitude); 1884 int signBits = Float.floatToRawIntBits((sign != sign) ? 1.0f : sign); 1885 magnitudeBits = (magnitudeBits & ~Float.SIGN_MASK) 1886 | (signBits & Float.SIGN_MASK); 1887 return Float.intBitsToFloat(magnitudeBits); 1888 } 1889 1890 /** 1891 * Returns the exponent of float {@code f}. 1892 * 1893 * @since 1.6 1894 */ 1895 public static int getExponent(float f) { 1896 return Math.getExponent(f); 1897 } 1898 1899 /** 1900 * Returns the exponent of double {@code d}. 1901 * 1902 * @since 1.6 1903 */ 1904 public static int getExponent(double d) { 1905 return Math.getExponent(d); 1906 } 1907 1908 /** 1909 * Returns the next double after {@code start} in the given 1910 * {@code direction}. 1911 * 1912 * @since 1.6 1913 */ 1914 public static double nextAfter(double start, double direction) { 1915 if (start == 0 && direction == 0) { 1916 return direction; 1917 } 1918 return nextafter(start, direction); 1919 } 1920 1921 /** 1922 * Returns the next float after {@code start} in the given {@code direction} 1923 * . 1924 * 1925 * @since 1.6 1926 */ 1927 public static float nextAfter(float start, double direction) { 1928 return Math.nextAfter(start, direction); 1929 } 1930 1931 /** 1932 * Returns the next double larger than {@code d}. 1933 * 1934 * @since 1.6 1935 */ 1936 public static double nextUp(double d) { 1937 return Math.nextUp(d); 1938 } 1939 1940 /** 1941 * Returns the next float larger than {@code f}. 1942 * 1943 * @since 1.6 1944 */ 1945 public static float nextUp(float f) { 1946 return Math.nextUp(f); 1947 } 1948 1949 /** 1950 * Returns {@code d} * 2^{@code scaleFactor}. The result may be rounded. 1951 * 1952 * @since 1.6 1953 */ 1954 public static double scalb(double d, int scaleFactor) { 1955 if (Double.isNaN(d) || Double.isInfinite(d) || d == 0) { 1956 return d; 1957 } 1958 // change double to long for calculation 1959 long bits = Double.doubleToLongBits(d); 1960 // the sign of the results must be the same of given d 1961 long sign = bits & Double.SIGN_MASK; 1962 // calculates the factor of the result 1963 long factor = (int) ((bits & Double.EXPONENT_MASK) >> Double.MANTISSA_BITS) 1964 - Double.EXPONENT_BIAS + scaleFactor; 1965 1966 // calculates the factor of sub-normal values 1967 int subNormalFactor = Long.numberOfLeadingZeros(bits & ~Double.SIGN_MASK) 1968 - Double.EXPONENT_BITS; 1969 if (subNormalFactor < 0) { 1970 // not sub-normal values 1971 subNormalFactor = 0; 1972 } 1973 if (Math.abs(d) < Double.MIN_NORMAL) { 1974 factor = factor - subNormalFactor; 1975 } 1976 if (factor > Double.MAX_EXPONENT) { 1977 return (d > 0 ? Double.POSITIVE_INFINITY : Double.NEGATIVE_INFINITY); 1978 } 1979 1980 long result; 1981 // if result is a sub-normal 1982 if (factor < -Double.EXPONENT_BIAS) { 1983 // the number of digits that shifts 1984 long digits = factor + Double.EXPONENT_BIAS + subNormalFactor; 1985 if (Math.abs(d) < Double.MIN_NORMAL) { 1986 // origin d is already sub-normal 1987 result = shiftLongBits(bits & Double.MANTISSA_MASK, digits); 1988 } else { 1989 // origin d is not sub-normal, change mantissa to sub-normal 1990 result = shiftLongBits(bits & Double.MANTISSA_MASK | 0x0010000000000000L, digits - 1); 1991 } 1992 } else { 1993 if (Math.abs(d) >= Double.MIN_NORMAL) { 1994 // common situation 1995 result = ((factor + Double.EXPONENT_BIAS) << Double.MANTISSA_BITS) | (bits & Double.MANTISSA_MASK); 1996 } else { 1997 // origin d is sub-normal, change mantissa to normal style 1998 result = ((factor + Double.EXPONENT_BIAS) << Double.MANTISSA_BITS) | ((bits << (subNormalFactor + 1)) & Double.MANTISSA_MASK); 1999 } 2000 } 2001 return Double.longBitsToDouble(result | sign); 2002 } 2003 2004 /** 2005 * Returns {@code d} * 2^{@code scaleFactor}. The result may be rounded. 2006 * 2007 * @since 1.6 2008 */ 2009 public static float scalb(float d, int scaleFactor) { 2010 if (Float.isNaN(d) || Float.isInfinite(d) || d == 0) { 2011 return d; 2012 } 2013 int bits = Float.floatToIntBits(d); 2014 int sign = bits & Float.SIGN_MASK; 2015 int factor = ((bits & Float.EXPONENT_MASK) >> Float.MANTISSA_BITS) 2016 - Float.EXPONENT_BIAS + scaleFactor; 2017 // calculates the factor of sub-normal values 2018 int subNormalFactor = Integer.numberOfLeadingZeros(bits & ~Float.SIGN_MASK) - Float.EXPONENT_BITS; 2019 if (subNormalFactor < 0) { 2020 // not sub-normal values 2021 subNormalFactor = 0; 2022 } 2023 if (Math.abs(d) < Float.MIN_NORMAL) { 2024 factor = factor - subNormalFactor; 2025 } 2026 if (factor > Float.MAX_EXPONENT) { 2027 return (d > 0 ? Float.POSITIVE_INFINITY : Float.NEGATIVE_INFINITY); 2028 } 2029 2030 int result; 2031 // if result is a sub-normal 2032 if (factor < -Float.EXPONENT_BIAS) { 2033 // the number of digits that shifts 2034 int digits = factor + Float.EXPONENT_BIAS + subNormalFactor; 2035 if (Math.abs(d) < Float.MIN_NORMAL) { 2036 // origin d is already sub-normal 2037 result = shiftIntBits(bits & Float.MANTISSA_MASK, digits); 2038 } else { 2039 // origin d is not sub-normal, change mantissa to sub-normal 2040 result = shiftIntBits(bits & Float.MANTISSA_MASK | 0x00800000, digits - 1); 2041 } 2042 } else { 2043 if (Math.abs(d) >= Float.MIN_NORMAL) { 2044 // common situation 2045 result = ((factor + Float.EXPONENT_BIAS) << Float.MANTISSA_BITS) 2046 | (bits & Float.MANTISSA_MASK); 2047 } else { 2048 // origin d is sub-normal, change mantissa to normal style 2049 result = ((factor + Float.EXPONENT_BIAS) 2050 << Float.MANTISSA_BITS) | ( 2051 (bits << (subNormalFactor + 1)) & Float.MANTISSA_MASK); 2052 } 2053 } 2054 return Float.intBitsToFloat(result | sign); 2055 } 2056 2057 // Shifts integer bits as float, if the digits is positive, left-shift; if 2058 // not, shift to right and calculate its carry. 2059 private static int shiftIntBits(int bits, int digits) { 2060 if (digits > 0) { 2061 return bits << digits; 2062 } 2063 // change it to positive 2064 int absDigits = -digits; 2065 if (Integer.numberOfLeadingZeros(bits & ~Float.SIGN_MASK) 2066 <= (32 - absDigits)) { 2067 // some bits will remain after shifting, calculates its carry 2068 if ((((bits >> (absDigits - 1)) & 0x1) == 0) || Integer.numberOfTrailingZeros(bits) == (absDigits - 1)) { 2069 return bits >> absDigits; 2070 } 2071 return ((bits >> absDigits) + 1); 2072 } 2073 return 0; 2074 } 2075 2076 // Shifts long bits as double, if the digits is positive, left-shift; if 2077 // not, shift to right and calculate its carry. 2078 private static long shiftLongBits(long bits, long digits) { 2079 if (digits > 0) { 2080 return bits << digits; 2081 } 2082 // change it to positive 2083 long absDigits = -digits; 2084 if (Long.numberOfLeadingZeros(bits & ~Double.SIGN_MASK) 2085 <= (64 - absDigits)) { 2086 // some bits will remain after shifting, calculates its carry 2087 if ((((bits >> (absDigits - 1)) & 0x1) == 0) || Long.numberOfTrailingZeros(bits) == (absDigits - 1)) { 2088 return bits >> absDigits; 2089 } 2090 return ((bits >> absDigits) + 1); 2091 } 2092 return 0; 2093 } 2094 } 2095
