Return-Path: X-Original-To: apmail-commons-notifications-archive@minotaur.apache.org Delivered-To: apmail-commons-notifications-archive@minotaur.apache.org Received: from mail.apache.org (hermes.apache.org [140.211.11.3]) by minotaur.apache.org (Postfix) with SMTP id 7D5D1187C5 for ; Mon, 27 Jul 2015 19:40:08 +0000 (UTC) Received: (qmail 27673 invoked by uid 500); 27 Jul 2015 19:40:08 -0000 Delivered-To: apmail-commons-notifications-archive@commons.apache.org Received: (qmail 27629 invoked by uid 500); 27 Jul 2015 19:40:08 -0000 Mailing-List: contact notifications-help@commons.apache.org; run by ezmlm Precedence: bulk List-Help: List-Unsubscribe: List-Post: List-Id: Reply-To: dev@commons.apache.org Delivered-To: mailing list notifications@commons.apache.org Received: (qmail 27515 invoked by uid 99); 27 Jul 2015 19:40:08 -0000 Received: from eris.apache.org (HELO hades.apache.org) (140.211.11.105) by apache.org (qpsmtpd/0.29) with ESMTP; Mon, 27 Jul 2015 19:40:08 +0000 Received: from hades.apache.org (localhost [127.0.0.1]) by hades.apache.org (ASF Mail Server at hades.apache.org) with ESMTP id 2F740AC186D for ; Mon, 27 Jul 2015 19:40:08 +0000 (UTC) Content-Type: text/plain; charset="utf-8" MIME-Version: 1.0 Content-Transfer-Encoding: 7bit Subject: svn commit: r959801 [8/10] - in /websites/production/commons/content/proper/commons-math/apidocs: ./ org/apache/commons/math3/exception/class-use/ org/apache/commons/math3/fraction/ org/apache/commons/math3/geometry/euclidean/oned/class-use/ org/apache... Date: Mon, 27 Jul 2015 19:40:06 -0000 To: notifications@commons.apache.org From: luc@apache.org X-Mailer: svnmailer-1.0.9 Message-Id: <20150727194008.2F740AC186D@hades.apache.org> Modified: websites/production/commons/content/proper/commons-math/apidocs/src-html/org/apache/commons/math3/ode/events/EventState.html ============================================================================== --- websites/production/commons/content/proper/commons-math/apidocs/src-html/org/apache/commons/math3/ode/events/EventState.html (original) +++ websites/production/commons/content/proper/commons-math/apidocs/src-html/org/apache/commons/math3/ode/events/EventState.html Mon Jul 27 19:40:06 2015 @@ -304,127 +304,138 @@ 296 ta = forward ? ta + convergence : ta - convergence; 297 ga = f.value(ta); 298 } while ((g0Positive ^ (ga >= 0)) && (forward ^ (ta >= tb))); -299 --i; -300 } else if (Double.isNaN(previousEventTime) || -301 (FastMath.abs(previousEventTime - root) > convergence)) { -302 pendingEventTime = root; -303 pendingEvent = true; -304 return true; -305 } else { -306 // no sign change: there is no event for now -307 ta = tb; -308 ga = gb; -309 } -310 -311 } else { -312 // no sign change: there is no event for now -313 ta = tb; -314 ga = gb; -315 } -316 -317 } -318 -319 // no event during the whole step -320 pendingEvent = false; -321 pendingEventTime = Double.NaN; -322 return false; -323 -324 } catch (LocalMaxCountExceededException lmcee) { -325 throw lmcee.getException(); -326 } +299 +300 if (forward ^ (ta >= tb)) { +301 // we were able to skip this spurious root +302 --i; +303 } else { +304 // we can't avoid this root before the end of the step, +305 // we have to handle it despite it is close to the former one +306 // maybe we have two very close roots +307 pendingEventTime = root; +308 pendingEvent = true; +309 return true; +310 } +311 } else if (Double.isNaN(previousEventTime) || +312 (FastMath.abs(previousEventTime - root) > convergence)) { +313 pendingEventTime = root; +314 pendingEvent = true; +315 return true; +316 } else { +317 // no sign change: there is no event for now +318 ta = tb; +319 ga = gb; +320 } +321 +322 } else { +323 // no sign change: there is no event for now +324 ta = tb; +325 ga = gb; +326 } 327 -328 } +328 } 329 -330 /** Get the occurrence time of the event triggered in the current step. -331 * @return occurrence time of the event triggered in the current -332 * step or infinity if no events are triggered -333 */ -334 public double getEventTime() { -335 return pendingEvent ? -336 pendingEventTime : -337 (forward ? Double.POSITIVE_INFINITY : Double.NEGATIVE_INFINITY); -338 } -339 -340 /** Acknowledge the fact the step has been accepted by the integrator. -341 * @param t value of the independent <i>time</i> variable at the -342 * end of the step -343 * @param y array containing the current value of the state vector -344 * at the end of the step -345 */ -346 public void stepAccepted(final double t, final double[] y) { -347 -348 t0 = t; -349 g0 = handler.g(t, y); +330 // no event during the whole step +331 pendingEvent = false; +332 pendingEventTime = Double.NaN; +333 return false; +334 +335 } catch (LocalMaxCountExceededException lmcee) { +336 throw lmcee.getException(); +337 } +338 +339 } +340 +341 /** Get the occurrence time of the event triggered in the current step. +342 * @return occurrence time of the event triggered in the current +343 * step or infinity if no events are triggered +344 */ +345 public double getEventTime() { +346 return pendingEvent ? +347 pendingEventTime : +348 (forward ? Double.POSITIVE_INFINITY : Double.NEGATIVE_INFINITY); +349 } 350 -351 if (pendingEvent && (FastMath.abs(pendingEventTime - t) <= convergence)) { -352 // force the sign to its value "just after the event" -353 previousEventTime = t; -354 g0Positive = increasing; -355 nextAction = handler.eventOccurred(t, y, !(increasing ^ forward)); -356 } else { -357 g0Positive = g0 >= 0; -358 nextAction = EventHandler.Action.CONTINUE; -359 } -360 } +351 /** Acknowledge the fact the step has been accepted by the integrator. +352 * @param t value of the independent <i>time</i> variable at the +353 * end of the step +354 * @param y array containing the current value of the state vector +355 * at the end of the step +356 */ +357 public void stepAccepted(final double t, final double[] y) { +358 +359 t0 = t; +360 g0 = handler.g(t, y); 361 -362 /** Check if the integration should be stopped at the end of the -363 * current step. -364 * @return true if the integration should be stopped -365 */ -366 public boolean stop() { -367 return nextAction == EventHandler.Action.STOP; -368 } -369 -370 /** Let the event handler reset the state if it wants. -371 * @param t value of the independent <i>time</i> variable at the -372 * beginning of the next step -373 * @param y array were to put the desired state vector at the beginning -374 * of the next step -375 * @return true if the integrator should reset the derivatives too +362 if (pendingEvent && (FastMath.abs(pendingEventTime - t) <= convergence)) { +363 // force the sign to its value "just after the event" +364 previousEventTime = t; +365 g0Positive = increasing; +366 nextAction = handler.eventOccurred(t, y, !(increasing ^ forward)); +367 } else { +368 g0Positive = g0 >= 0; +369 nextAction = EventHandler.Action.CONTINUE; +370 } +371 } +372 +373 /** Check if the integration should be stopped at the end of the +374 * current step. +375 * @return true if the integration should be stopped 376 */ -377 public boolean reset(final double t, final double[] y) { -378 -379 if (!(pendingEvent && (FastMath.abs(pendingEventTime - t) <= convergence))) { -380 return false; -381 } -382 -383 if (nextAction == EventHandler.Action.RESET_STATE) { -384 handler.resetState(t, y); -385 } -386 pendingEvent = false; -387 pendingEventTime = Double.NaN; -388 -389 return (nextAction == EventHandler.Action.RESET_STATE) || -390 (nextAction == EventHandler.Action.RESET_DERIVATIVES); -391 -392 } +377 public boolean stop() { +378 return nextAction == EventHandler.Action.STOP; +379 } +380 +381 /** Let the event handler reset the state if it wants. +382 * @param t value of the independent <i>time</i> variable at the +383 * beginning of the next step +384 * @param y array were to put the desired state vector at the beginning +385 * of the next step +386 * @return true if the integrator should reset the derivatives too +387 */ +388 public boolean reset(final double t, final double[] y) { +389 +390 if (!(pendingEvent && (FastMath.abs(pendingEventTime - t) <= convergence))) { +391 return false; +392 } 393 -394 /** Local wrapper to propagate exceptions. */ -395 private static class LocalMaxCountExceededException extends RuntimeException { -396 -397 /** Serializable UID. */ -398 private static final long serialVersionUID = 20120901L; +394 if (nextAction == EventHandler.Action.RESET_STATE) { +395 handler.resetState(t, y); +396 } +397 pendingEvent = false; +398 pendingEventTime = Double.NaN; 399 -400 /** Wrapped exception. */ -401 private final MaxCountExceededException wrapped; +400 return (nextAction == EventHandler.Action.RESET_STATE) || +401 (nextAction == EventHandler.Action.RESET_DERIVATIVES); 402 -403 /** Simple constructor. -404 * @param exception exception to wrap -405 */ -406 public LocalMaxCountExceededException(final MaxCountExceededException exception) { -407 wrapped = exception; -408 } -409 -410 /** Get the wrapped exception. -411 * @return wrapped exception -412 */ -413 public MaxCountExceededException getException() { -414 return wrapped; -415 } -416 -417 } -418 -419} +403 } +404 +405 /** Local wrapper to propagate exceptions. */ +406 private static class LocalMaxCountExceededException extends RuntimeException { +407 +408 /** Serializable UID. */ +409 private static final long serialVersionUID = 20120901L; +410 +411 /** Wrapped exception. */ +412 private final MaxCountExceededException wrapped; +413 +414 /** Simple constructor. +415 * @param exception exception to wrap +416 */ +417 public LocalMaxCountExceededException(final MaxCountExceededException exception) { +418 wrapped = exception; +419 } +420 +421 /** Get the wrapped exception. +422 * @return wrapped exception +423 */ +424 public MaxCountExceededException getException() { +425 return wrapped; +426 } +427 +428 } +429 +430} Modified: websites/production/commons/content/proper/commons-math/apidocs/src-html/org/apache/commons/math3/special/Gamma.html ============================================================================== --- websites/production/commons/content/proper/commons-math/apidocs/src-html/org/apache/commons/math3/special/Gamma.html (original) +++ websites/production/commons/content/proper/commons-math/apidocs/src-html/org/apache/commons/math3/special/Gamma.html Mon Jul 27 19:40:06 2015 @@ -450,272 +450,280 @@ 442 * @since 2.0 443 */ 444 public static double digamma(double x) { -445 if (x > 0 && x <= S_LIMIT) { -446 // use method 5 from Bernardo AS103 -447 // accurate to O(x) -448 return -GAMMA - 1 / x; -449 } -450 -451 if (x >= C_LIMIT) { -452 // use method 4 (accurate to O(1/x^8) -453 double inv = 1 / (x * x); -454 // 1 1 1 1 -455 // log(x) - --- - ------ + ------- - ------- -456 // 2 x 12 x^2 120 x^4 252 x^6 -457 return FastMath.log(x) - 0.5 / x - inv * ((1.0 / 12) + inv * (1.0 / 120 - inv / 252)); -458 } -459 -460 return digamma(x + 1) - 1 / x; -461 } -462 -463 /** -464 * Computes the trigamma function of x. -465 * This function is derived by taking the derivative of the implementation -466 * of digamma. -467 * -468 * @param x Argument. -469 * @return trigamma(x) to within 10-8 relative or absolute error whichever is smaller -470 * @see <a href="http://en.wikipedia.org/wiki/Trigamma_function">Trigamma</a> -471 * @see Gamma#digamma(double) -472 * @since 2.0 -473 */ -474 public static double trigamma(double x) { -475 if (x > 0 && x <= S_LIMIT) { -476 return 1 / (x * x); -477 } -478 -479 if (x >= C_LIMIT) { -480 double inv = 1 / (x * x); -481 // 1 1 1 1 1 -482 // - + ---- + ---- - ----- + ----- -483 // x 2 3 5 7 -484 // 2 x 6 x 30 x 42 x -485 return 1 / x + inv / 2 + inv / x * (1.0 / 6 - inv * (1.0 / 30 + inv / 42)); -486 } -487 -488 return trigamma(x + 1) + 1 / (x * x); -489 } -490 -491 /** -492 * <p> -493 * Returns the Lanczos approximation used to compute the gamma function. -494 * The Lanczos approximation is related to the Gamma function by the -495 * following equation -496 * <center> -497 * {@code gamma(x) = sqrt(2 * pi) / x * (x + g + 0.5) ^ (x + 0.5) -498 * * exp(-x - g - 0.5) * lanczos(x)}, -499 * </center> -500 * where {@code g} is the Lanczos constant. -501 * </p> -502 * -503 * @param x Argument. -504 * @return The Lanczos approximation. -505 * @see <a href="http://mathworld.wolfram.com/LanczosApproximation.html">Lanczos Approximation</a> -506 * equations (1) through (5), and Paul Godfrey's -507 * <a href="http://my.fit.edu/~gabdo/gamma.txt">Note on the computation -508 * of the convergent Lanczos complex Gamma approximation</a> -509 * @since 3.1 -510 */ -511 public static double lanczos(final double x) { -512 double sum = 0.0; -513 for (int i = LANCZOS.length - 1; i > 0; --i) { -514 sum += LANCZOS[i] / (x + i); -515 } -516 return sum + LANCZOS[0]; -517 } -518 -519 /** -520 * Returns the value of 1 / &Gamma;(1 + x) - 1 for -0&#46;5 &le; x &le; -521 * 1&#46;5. This implementation is based on the double precision -522 * implementation in the <em>NSWC Library of Mathematics Subroutines</em>, -523 * {@code DGAM1}. -524 * -525 * @param x Argument. -526 * @return The value of {@code 1.0 / Gamma(1.0 + x) - 1.0}. -527 * @throws NumberIsTooSmallException if {@code x < -0.5} -528 * @throws NumberIsTooLargeException if {@code x > 1.5} -529 * @since 3.1 -530 */ -531 public static double invGamma1pm1(final double x) { -532 -533 if (x < -0.5) { -534 throw new NumberIsTooSmallException(x, -0.5, true); -535 } -536 if (x > 1.5) { -537 throw new NumberIsTooLargeException(x, 1.5, true); -538 } -539 -540 final double ret; -541 final double t = x <= 0.5 ? x : (x - 0.5) - 0.5; -542 if (t < 0.0) { -543 final double a = INV_GAMMA1P_M1_A0 + t * INV_GAMMA1P_M1_A1; -544 double b = INV_GAMMA1P_M1_B8; -545 b = INV_GAMMA1P_M1_B7 + t * b; -546 b = INV_GAMMA1P_M1_B6 + t * b; -547 b = INV_GAMMA1P_M1_B5 + t * b; -548 b = INV_GAMMA1P_M1_B4 + t * b; -549 b = INV_GAMMA1P_M1_B3 + t * b; -550 b = INV_GAMMA1P_M1_B2 + t * b; -551 b = INV_GAMMA1P_M1_B1 + t * b; -552 b = 1.0 + t * b; -553 -554 double c = INV_GAMMA1P_M1_C13 + t * (a / b); -555 c = INV_GAMMA1P_M1_C12 + t * c; -556 c = INV_GAMMA1P_M1_C11 + t * c; -557 c = INV_GAMMA1P_M1_C10 + t * c; -558 c = INV_GAMMA1P_M1_C9 + t * c; -559 c = INV_GAMMA1P_M1_C8 + t * c; -560 c = INV_GAMMA1P_M1_C7 + t * c; -561 c = INV_GAMMA1P_M1_C6 + t * c; -562 c = INV_GAMMA1P_M1_C5 + t * c; -563 c = INV_GAMMA1P_M1_C4 + t * c; -564 c = INV_GAMMA1P_M1_C3 + t * c; -565 c = INV_GAMMA1P_M1_C2 + t * c; -566 c = INV_GAMMA1P_M1_C1 + t * c; -567 c = INV_GAMMA1P_M1_C + t * c; -568 if (x > 0.5) { -569 ret = t * c / x; -570 } else { -571 ret = x * ((c + 0.5) + 0.5); -572 } -573 } else { -574 double p = INV_GAMMA1P_M1_P6; -575 p = INV_GAMMA1P_M1_P5 + t * p; -576 p = INV_GAMMA1P_M1_P4 + t * p; -577 p = INV_GAMMA1P_M1_P3 + t * p; -578 p = INV_GAMMA1P_M1_P2 + t * p; -579 p = INV_GAMMA1P_M1_P1 + t * p; -580 p = INV_GAMMA1P_M1_P0 + t * p; -581 -582 double q = INV_GAMMA1P_M1_Q4; -583 q = INV_GAMMA1P_M1_Q3 + t * q; -584 q = INV_GAMMA1P_M1_Q2 + t * q; -585 q = INV_GAMMA1P_M1_Q1 + t * q; -586 q = 1.0 + t * q; -587 -588 double c = INV_GAMMA1P_M1_C13 + (p / q) * t; -589 c = INV_GAMMA1P_M1_C12 + t * c; -590 c = INV_GAMMA1P_M1_C11 + t * c; -591 c = INV_GAMMA1P_M1_C10 + t * c; -592 c = INV_GAMMA1P_M1_C9 + t * c; -593 c = INV_GAMMA1P_M1_C8 + t * c; -594 c = INV_GAMMA1P_M1_C7 + t * c; -595 c = INV_GAMMA1P_M1_C6 + t * c; -596 c = INV_GAMMA1P_M1_C5 + t * c; -597 c = INV_GAMMA1P_M1_C4 + t * c; -598 c = INV_GAMMA1P_M1_C3 + t * c; -599 c = INV_GAMMA1P_M1_C2 + t * c; -600 c = INV_GAMMA1P_M1_C1 + t * c; -601 c = INV_GAMMA1P_M1_C0 + t * c; -602 -603 if (x > 0.5) { -604 ret = (t / x) * ((c - 0.5) - 0.5); -605 } else { -606 ret = x * c; -607 } -608 } -609 -610 return ret; -611 } -612 -613 /** -614 * Returns the value of log &Gamma;(1 + x) for -0&#46;5 &le; x &le; 1&#46;5. -615 * This implementation is based on the double precision implementation in -616 * the <em>NSWC Library of Mathematics Subroutines</em>, {@code DGMLN1}. -617 * -618 * @param x Argument. -619 * @return The value of {@code log(Gamma(1 + x))}. -620 * @throws NumberIsTooSmallException if {@code x < -0.5}. -621 * @throws NumberIsTooLargeException if {@code x > 1.5}. -622 * @since 3.1 -623 */ -624 public static double logGamma1p(final double x) -625 throws NumberIsTooSmallException, NumberIsTooLargeException { -626 -627 if (x < -0.5) { -628 throw new NumberIsTooSmallException(x, -0.5, true); -629 } -630 if (x > 1.5) { -631 throw new NumberIsTooLargeException(x, 1.5, true); -632 } -633 -634 return -FastMath.log1p(invGamma1pm1(x)); -635 } -636 -637 -638 /** -639 * Returns the value of ?(x). Based on the <em>NSWC Library of -640 * Mathematics Subroutines</em> double precision implementation, -641 * {@code DGAMMA}. -642 * -643 * @param x Argument. -644 * @return the value of {@code Gamma(x)}. -645 * @since 3.1 -646 */ -647 public static double gamma(final double x) { -648 -649 if ((x == FastMath.rint(x)) && (x <= 0.0)) { -650 return Double.NaN; -651 } -652 -653 final double ret; -654 final double absX = FastMath.abs(x); -655 if (absX <= 20.0) { -656 if (x >= 1.0) { -657 /* -658 * From the recurrence relation -659 * Gamma(x) = (x - 1) * ... * (x - n) * Gamma(x - n), -660 * then -661 * Gamma(t) = 1 / [1 + invGamma1pm1(t - 1)], -662 * where t = x - n. This means that t must satisfy -663 * -0.5 <= t - 1 <= 1.5. -664 */ -665 double prod = 1.0; -666 double t = x; -667 while (t > 2.5) { -668 t -= 1.0; -669 prod *= t; -670 } -671 ret = prod / (1.0 + invGamma1pm1(t - 1.0)); -672 } else { -673 /* -674 * From the recurrence relation -675 * Gamma(x) = Gamma(x + n + 1) / [x * (x + 1) * ... * (x + n)] -676 * then -677 * Gamma(x + n + 1) = 1 / [1 + invGamma1pm1(x + n)], -678 * which requires -0.5 <= x + n <= 1.5. -679 */ -680 double prod = x; -681 double t = x; -682 while (t < -0.5) { -683 t += 1.0; -684 prod *= t; -685 } -686 ret = 1.0 / (prod * (1.0 + invGamma1pm1(t))); -687 } -688 } else { -689 final double y = absX + LANCZOS_G + 0.5; -690 final double gammaAbs = SQRT_TWO_PI / x * -691 FastMath.pow(y, absX + 0.5) * -692 FastMath.exp(-y) * lanczos(absX); -693 if (x > 0.0) { -694 ret = gammaAbs; -695 } else { -696 /* -697 * From the reflection formula -698 * Gamma(x) * Gamma(1 - x) * sin(pi * x) = pi, -699 * and the recurrence relation -700 * Gamma(1 - x) = -x * Gamma(-x), -701 * it is found -702 * Gamma(x) = -pi / [x * sin(pi * x) * Gamma(-x)]. -703 */ -704 ret = -FastMath.PI / -705 (x * FastMath.sin(FastMath.PI * x) * gammaAbs); -706 } -707 } -708 return ret; -709 } -710} +445 if (Double.isNaN(x) || Double.isInfinite(x)) { +446 return x; +447 } +448 +449 if (x > 0 && x <= S_LIMIT) { +450 // use method 5 from Bernardo AS103 +451 // accurate to O(x) +452 return -GAMMA - 1 / x; +453 } +454 +455 if (x >= C_LIMIT) { +456 // use method 4 (accurate to O(1/x^8) +457 double inv = 1 / (x * x); +458 // 1 1 1 1 +459 // log(x) - --- - ------ + ------- - ------- +460 // 2 x 12 x^2 120 x^4 252 x^6 +461 return FastMath.log(x) - 0.5 / x - inv * ((1.0 / 12) + inv * (1.0 / 120 - inv / 252)); +462 } +463 +464 return digamma(x + 1) - 1 / x; +465 } +466 +467 /** +468 * Computes the trigamma function of x. +469 * This function is derived by taking the derivative of the implementation +470 * of digamma. +471 * +472 * @param x Argument. +473 * @return trigamma(x) to within 10-8 relative or absolute error whichever is smaller +474 * @see <a href="http://en.wikipedia.org/wiki/Trigamma_function">Trigamma</a> +475 * @see Gamma#digamma(double) +476 * @since 2.0 +477 */ +478 public static double trigamma(double x) { +479 if (Double.isNaN(x) || Double.isInfinite(x)) { +480 return x; +481 } +482 +483 if (x > 0 && x <= S_LIMIT) { +484 return 1 / (x * x); +485 } +486 +487 if (x >= C_LIMIT) { +488 double inv = 1 / (x * x); +489 // 1 1 1 1 1 +490 // - + ---- + ---- - ----- + ----- +491 // x 2 3 5 7 +492 // 2 x 6 x 30 x 42 x +493 return 1 / x + inv / 2 + inv / x * (1.0 / 6 - inv * (1.0 / 30 + inv / 42)); +494 } +495 +496 return trigamma(x + 1) + 1 / (x * x); +497 } +498 +499 /** +500 * <p> +501 * Returns the Lanczos approximation used to compute the gamma function. +502 * The Lanczos approximation is related to the Gamma function by the +503 * following equation +504 * <center> +505 * {@code gamma(x) = sqrt(2 * pi) / x * (x + g + 0.5) ^ (x + 0.5) +506 * * exp(-x - g - 0.5) * lanczos(x)}, +507 * </center> +508 * where {@code g} is the Lanczos constant. +509 * </p> +510 * +511 * @param x Argument. +512 * @return The Lanczos approximation. +513 * @see <a href="http://mathworld.wolfram.com/LanczosApproximation.html">Lanczos Approximation</a> +514 * equations (1) through (5), and Paul Godfrey's +515 * <a href="http://my.fit.edu/~gabdo/gamma.txt">Note on the computation +516 * of the convergent Lanczos complex Gamma approximation</a> +517 * @since 3.1 +518 */ +519 public static double lanczos(final double x) { +520 double sum = 0.0; +521 for (int i = LANCZOS.length - 1; i > 0; --i) { +522 sum += LANCZOS[i] / (x + i); +523 } +524 return sum + LANCZOS[0]; +525 } +526 +527 /** +528 * Returns the value of 1 / &Gamma;(1 + x) - 1 for -0&#46;5 &le; x &le; +529 * 1&#46;5. This implementation is based on the double precision +530 * implementation in the <em>NSWC Library of Mathematics Subroutines</em>, +531 * {@code DGAM1}. +532 * +533 * @param x Argument. +534 * @return The value of {@code 1.0 / Gamma(1.0 + x) - 1.0}. +535 * @throws NumberIsTooSmallException if {@code x < -0.5} +536 * @throws NumberIsTooLargeException if {@code x > 1.5} +537 * @since 3.1 +538 */ +539 public static double invGamma1pm1(final double x) { +540 +541 if (x < -0.5) { +542 throw new NumberIsTooSmallException(x, -0.5, true); +543 } +544 if (x > 1.5) { +545 throw new NumberIsTooLargeException(x, 1.5, true); +546 } +547 +548 final double ret; +549 final double t = x <= 0.5 ? x : (x - 0.5) - 0.5; +550 if (t < 0.0) { +551 final double a = INV_GAMMA1P_M1_A0 + t * INV_GAMMA1P_M1_A1; +552 double b = INV_GAMMA1P_M1_B8; +553 b = INV_GAMMA1P_M1_B7 + t * b; +554 b = INV_GAMMA1P_M1_B6 + t * b; +555 b = INV_GAMMA1P_M1_B5 + t * b; +556 b = INV_GAMMA1P_M1_B4 + t * b; +557 b = INV_GAMMA1P_M1_B3 + t * b; +558 b = INV_GAMMA1P_M1_B2 + t * b; +559 b = INV_GAMMA1P_M1_B1 + t * b; +560 b = 1.0 + t * b; +561 +562 double c = INV_GAMMA1P_M1_C13 + t * (a / b); +563 c = INV_GAMMA1P_M1_C12 + t * c; +564 c = INV_GAMMA1P_M1_C11 + t * c; +565 c = INV_GAMMA1P_M1_C10 + t * c; +566 c = INV_GAMMA1P_M1_C9 + t * c; +567 c = INV_GAMMA1P_M1_C8 + t * c; +568 c = INV_GAMMA1P_M1_C7 + t * c; +569 c = INV_GAMMA1P_M1_C6 + t * c; +570 c = INV_GAMMA1P_M1_C5 + t * c; +571 c = INV_GAMMA1P_M1_C4 + t * c; +572 c = INV_GAMMA1P_M1_C3 + t * c; +573 c = INV_GAMMA1P_M1_C2 + t * c; +574 c = INV_GAMMA1P_M1_C1 + t * c; +575 c = INV_GAMMA1P_M1_C + t * c; +576 if (x > 0.5) { +577 ret = t * c / x; +578 } else { +579 ret = x * ((c + 0.5) + 0.5); +580 } +581 } else { +582 double p = INV_GAMMA1P_M1_P6; +583 p = INV_GAMMA1P_M1_P5 + t * p; +584 p = INV_GAMMA1P_M1_P4 + t * p; +585 p = INV_GAMMA1P_M1_P3 + t * p; +586 p = INV_GAMMA1P_M1_P2 + t * p; +587 p = INV_GAMMA1P_M1_P1 + t * p; +588 p = INV_GAMMA1P_M1_P0 + t * p; +589 +590 double q = INV_GAMMA1P_M1_Q4; +591 q = INV_GAMMA1P_M1_Q3 + t * q; +592 q = INV_GAMMA1P_M1_Q2 + t * q; +593 q = INV_GAMMA1P_M1_Q1 + t * q; +594 q = 1.0 + t * q; +595 +596 double c = INV_GAMMA1P_M1_C13 + (p / q) * t; +597 c = INV_GAMMA1P_M1_C12 + t * c; +598 c = INV_GAMMA1P_M1_C11 + t * c; +599 c = INV_GAMMA1P_M1_C10 + t * c; +600 c = INV_GAMMA1P_M1_C9 + t * c; +601 c = INV_GAMMA1P_M1_C8 + t * c; +602 c = INV_GAMMA1P_M1_C7 + t * c; +603 c = INV_GAMMA1P_M1_C6 + t * c; +604 c = INV_GAMMA1P_M1_C5 + t * c; +605 c = INV_GAMMA1P_M1_C4 + t * c; +606 c = INV_GAMMA1P_M1_C3 + t * c; +607 c = INV_GAMMA1P_M1_C2 + t * c; +608 c = INV_GAMMA1P_M1_C1 + t * c; +609 c = INV_GAMMA1P_M1_C0 + t * c; +610 +611 if (x > 0.5) { +612 ret = (t / x) * ((c - 0.5) - 0.5); +613 } else { +614 ret = x * c; +615 } +616 } +617 +618 return ret; +619 } +620 +621 /** +622 * Returns the value of log &Gamma;(1 + x) for -0&#46;5 &le; x &le; 1&#46;5. +623 * This implementation is based on the double precision implementation in +624 * the <em>NSWC Library of Mathematics Subroutines</em>, {@code DGMLN1}. +625 * +626 * @param x Argument. +627 * @return The value of {@code log(Gamma(1 + x))}. +628 * @throws NumberIsTooSmallException if {@code x < -0.5}. +629 * @throws NumberIsTooLargeException if {@code x > 1.5}. +630 * @since 3.1 +631 */ +632 public static double logGamma1p(final double x) +633 throws NumberIsTooSmallException, NumberIsTooLargeException { +634 +635 if (x < -0.5) { +636 throw new NumberIsTooSmallException(x, -0.5, true); +637 } +638 if (x > 1.5) { +639 throw new NumberIsTooLargeException(x, 1.5, true); +640 } +641 +642 return -FastMath.log1p(invGamma1pm1(x)); +643 } +644 +645 +646 /** +647 * Returns the value of ?(x). Based on the <em>NSWC Library of +648 * Mathematics Subroutines</em> double precision implementation, +649 * {@code DGAMMA}. +650 * +651 * @param x Argument. +652 * @return the value of {@code Gamma(x)}. +653 * @since 3.1 +654 */ +655 public static double gamma(final double x) { +656 +657 if ((x == FastMath.rint(x)) && (x <= 0.0)) { +658 return Double.NaN; +659 } +660 +661 final double ret; +662 final double absX = FastMath.abs(x); +663 if (absX <= 20.0) { +664 if (x >= 1.0) { +665 /* +666 * From the recurrence relation +667 * Gamma(x) = (x - 1) * ... * (x - n) * Gamma(x - n), +668 * then +669 * Gamma(t) = 1 / [1 + invGamma1pm1(t - 1)], +670 * where t = x - n. This means that t must satisfy +671 * -0.5 <= t - 1 <= 1.5. +672 */ +673 double prod = 1.0; +674 double t = x; +675 while (t > 2.5) { +676 t -= 1.0; +677 prod *= t; +678 } +679 ret = prod / (1.0 + invGamma1pm1(t - 1.0)); +680 } else { +681 /* +682 * From the recurrence relation +683 * Gamma(x) = Gamma(x + n + 1) / [x * (x + 1) * ... * (x + n)] +684 * then +685 * Gamma(x + n + 1) = 1 / [1 + invGamma1pm1(x + n)], +686 * which requires -0.5 <= x + n <= 1.5. +687 */ +688 double prod = x; +689 double t = x; +690 while (t < -0.5) { +691 t += 1.0; +692 prod *= t; +693 } +694 ret = 1.0 / (prod * (1.0 + invGamma1pm1(t))); +695 } +696 } else { +697 final double y = absX + LANCZOS_G + 0.5; +698 final double gammaAbs = SQRT_TWO_PI / x * +699 FastMath.pow(y, absX + 0.5) * +700 FastMath.exp(-y) * lanczos(absX); +701 if (x > 0.0) { +702 ret = gammaAbs; +703 } else { +704 /* +705 * From the reflection formula +706 * Gamma(x) * Gamma(1 - x) * sin(pi * x) = pi, +707 * and the recurrence relation +708 * Gamma(1 - x) = -x * Gamma(-x), +709 * it is found +710 * Gamma(x) = -pi / [x * sin(pi * x) * Gamma(-x)]. +711 */ +712 ret = -FastMath.PI / +713 (x * FastMath.sin(FastMath.PI * x) * gammaAbs); +714 } +715 } +716 return ret; +717 } +718}