math.Arithmetic Maven / Gradle / Ivy
Go to download
Show more of this group Show more artifacts with this name
Show all versions of finwhale Show documentation
Show all versions of finwhale Show documentation
Statistical distributions library (in statu nascendi)
/*
Copyright ? 1999 CERN - European Organization for Nuclear Research.
Permission to use, copy, modify, distribute and sell this software and its documentation for any purpose
is hereby granted without fee, provided that the above copyright notice appear in all copies and
that both that copyright notice and this permission notice appear in supporting documentation.
CERN makes no representations about the suitability of this software for any purpose.
It is provided "as is" without expressed or implied warranty.
*/
package math;
import java.math.BigDecimal;
import java.math.RoundingMode;
import math.rng.DefaultRng;
/**
* Arithmetic functions.
*/
public final class Arithmetic {
// for method stirlingCorrection(...)
private static final double[] stirlingCorrection = { 0.0,
8.106146679532726e-02, 4.134069595540929e-02,
2.767792568499834e-02, 2.079067210376509e-02,
1.664469118982119e-02, 1.387612882307075e-02,
1.189670994589177e-02, 1.041126526197209e-02,
9.255462182712733e-03, 8.330563433362871e-03,
7.573675487951841e-03, 6.942840107209530e-03,
6.408994188004207e-03, 5.951370112758848e-03,
5.554733551962801e-03, 5.207655919609640e-03,
4.901395948434738e-03, 4.629153749334029e-03,
4.385560249232324e-03, 4.166319691996922e-03,
3.967954218640860e-03, 3.787618068444430e-03,
3.622960224683090e-03, 3.472021382978770e-03,
3.333155636728090e-03, 3.204970228055040e-03,
3.086278682608780e-03, 2.976063983550410e-03,
2.873449362352470e-03, 2.777674929752690e-03, };
// for method logFactorial(...)
// log(k!) for k = 0, ..., 29
private static final double[] logFactorials = { 0.00000000000000000,
0.00000000000000000, 0.69314718055994531, 1.79175946922805500,
3.17805383034794562, 4.78749174278204599, 6.57925121201010100,
8.52516136106541430, 10.60460290274525023, 12.80182748008146961,
15.10441257307551530, 17.50230784587388584, 19.98721449566188615,
22.55216385312342289, 25.19122118273868150, 27.89927138384089157,
30.67186010608067280, 33.50507345013688888, 36.39544520803305358,
39.33988418719949404, 42.33561646075348503, 45.38013889847690803,
48.47118135183522388, 51.60667556776437357, 54.78472939811231919,
58.00360522298051994, 61.26170176100200198, 64.55753862700633106,
67.88974313718153498, 71.25703896716800901 };
// k! for k = 0, ..., 20
private static final long[] longFactorials = { 1L, 1L, 2L, 6L, 24L, 120L,
720L, 5040L, 40320L, 362880L, 3628800L, 39916800L, 479001600L,
6227020800L, 87178291200L, 1307674368000L, 20922789888000L,
355687428096000L, 6402373705728000L, 121645100408832000L,
2432902008176640000L };
// k! for k = 21, ..., 170
private static final double[] doubleFactorials = { 5.109094217170944E19,
1.1240007277776077E21, 2.585201673888498E22, 6.204484017332394E23,
1.5511210043330984E25, 4.032914611266057E26, 1.0888869450418352E28,
3.048883446117138E29, 8.841761993739701E30, 2.652528598121911E32,
8.222838654177924E33, 2.6313083693369355E35, 8.68331761881189E36,
2.952327990396041E38, 1.0333147966386144E40, 3.719933267899013E41,
1.3763753091226346E43, 5.23022617466601E44, 2.0397882081197447E46,
8.15915283247898E47, 3.34525266131638E49, 1.4050061177528801E51,
6.041526306337384E52, 2.6582715747884495E54, 1.196222208654802E56,
5.502622159812089E57, 2.5862324151116827E59, 1.2413915592536068E61,
6.082818640342679E62, 3.0414093201713376E64, 1.5511187532873816E66,
8.06581751709439E67, 4.274883284060024E69, 2.308436973392413E71,
1.2696403353658264E73, 7.109985878048632E74, 4.052691950487723E76,
2.350561331282879E78, 1.386831185456898E80, 8.32098711274139E81,
5.075802138772246E83, 3.146997326038794E85, 1.9826083154044396E87,
1.2688693218588414E89, 8.247650592082472E90, 5.443449390774432E92,
3.6471110918188705E94, 2.48003554243683E96, 1.7112245242814127E98,
1.1978571669969892E100, 8.504785885678624E101,
6.123445837688612E103, 4.470115461512686E105,
3.307885441519387E107, 2.4809140811395404E109,
1.8854947016660506E111, 1.451830920282859E113,
1.1324281178206295E115, 8.94618213078298E116, 7.15694570462638E118,
5.797126020747369E120, 4.7536433370128435E122,
3.94552396972066E124, 3.314240134565354E126,
2.8171041143805494E128, 2.4227095383672744E130,
2.107757298379527E132, 1.854826422573984E134,
1.6507955160908465E136, 1.4857159644817605E138,
1.3520015276784033E140, 1.2438414054641305E142,
1.156772507081641E144, 1.0873661566567426E146,
1.0329978488239061E148, 9.916779348709491E149,
9.619275968248216E151, 9.426890448883248E153,
9.332621544394415E155, 9.332621544394418E157, 9.42594775983836E159,
9.614466715035125E161, 9.902900716486178E163,
1.0299016745145631E166, 1.0813967582402912E168,
1.1462805637347086E170, 1.2265202031961373E172,
1.324641819451829E174, 1.4438595832024942E176,
1.5882455415227423E178, 1.7629525510902457E180,
1.974506857221075E182, 2.2311927486598138E184,
2.543559733472186E186, 2.925093693493014E188,
3.393108684451899E190, 3.96993716080872E192,
4.6845258497542896E194, 5.574585761207606E196,
6.689502913449135E198, 8.094298525273444E200,
9.875044200833601E202, 1.2146304367025332E205,
1.506141741511141E207, 1.882677176888926E209,
2.3721732428800483E211, 3.0126600184576624E213,
3.856204823625808E215, 4.974504222477287E217,
6.466855489220473E219, 8.471580690878813E221,
1.1182486511960037E224, 1.4872707060906847E226,
1.99294274616152E228, 2.690472707318049E230,
3.6590428819525483E232, 5.0128887482749884E234,
6.917786472619482E236, 9.615723196941089E238,
1.3462012475717523E241, 1.8981437590761713E243,
2.6953641378881633E245, 3.8543707171800694E247,
5.550293832739308E249, 8.047926057471989E251,
1.1749972043909107E254, 1.72724589045464E256,
2.5563239178728637E258, 3.8089226376305687E260,
5.7133839564458575E262, 8.627209774233244E264,
1.3113358856834527E267, 2.0063439050956838E269,
3.0897696138473515E271, 4.789142901463393E273,
7.471062926282892E275, 1.1729568794264134E278,
1.8532718694937346E280, 2.946702272495036E282,
4.714723635992061E284, 7.590705053947223E286,
1.2296942187394494E289, 2.0044015765453032E291,
3.287218585534299E293, 5.423910666131583E295,
9.003691705778434E297, 1.5036165148649983E300,
2.5260757449731988E302, 4.2690680090047056E304,
7.257415615308004E306 };
// for method factoPow()
// n!/n^n for n = 0, ..., 500
//@formatter:off
private static final double[] factoPows = {
1.0, 1.0, 0.5, 0.2222222222222222, 0.09375, 0.0384,
0.015432098765432098, 0.006119899021666143, 0.00240325927734375,
9.36656708416885E-4, 3.6288E-4, 1.3990594886818849E-4,
5.372321709247828E-5, 2.055969825338982E-5, 7.845413755460156E-6,
2.9862813725549966E-6, 1.1342267125513672E-6, 4.299687097977245E-7,
1.627181239384263E-7, 6.148599467297143E-8, 2.32019615953125E-8,
8.744575307296545E-9, 3.2920590329493666E-9, 1.2380956566238489E-9,
4.651958899973984E-10, 1.7464069942802745E-10, 6.551066071283736E-11,
2.455619981919136E-11, 9.198477380186612E-12, 3.443474578407548E-12,
1.2883160975104712E-12, 4.817361748848426E-13, 1.8004142466391044E-13,
6.725516941137585E-14, 2.511200704691097E-14, 9.372406058407147E-15,
3.4965933035017217E-15, 1.3039864624082533E-15, 4.861203460210078E-16,
1.8116128928263687E-16, 6.749093037883734E-17, 2.513568929493553E-17,
9.35853822486767E-18, 3.483397861864393E-18, 1.2962285476089634E-18,
4.82223911246876E-19, 1.793533262585642E-19, 6.66911532924369E-20,
2.4793017907155796E-20, 9.215034101804921E-21, 3.424322470251202E-21,
1.2722312750921944E-21, 4.725791134077835E-22, 1.755105280903671E-22,
6.517109097038799E-23, 2.419539903336549E-23, 8.981300465842497E-24,
3.333318559896576E-24, 1.236938082646538E-24, 4.589389441131137E-25,
1.7025498093327752E-25, 6.315165678191623E-26, 2.3421333871715883E-26,
8.685249677894244E-27, 3.2203165380318905E-27, 1.1938838356471498E-27,
4.42562328829495E-28, 1.6403525997799646E-28, 6.079275720609595E-29,
2.2527849431516092E-29, 8.347227202377159E-30, 3.092577738453789E-30,
1.1456610223247468E-30, 4.2437515686630335E-31, 1.5718214119030096E-31,
5.821259634147174E-32, 2.1557197619542006E-32, 7.982339651754881E-33,
2.9555045641404006E-33, 1.094202048651836E-33, 4.050687717855241E-34,
1.4994300566878794E-34, 5.549970801498484E-35, 2.0541067935114653E-35,
7.601931343991518E-36, 2.8131585658572497E-36, 1.0409611999103514E-36,
3.851639429476385E-37, 1.425043510365707E-37, 5.272088466442438E-38,
1.9503382601179696E-38, 7.214570491098713E-39, 2.6686086775143567E-39,
9.870377041123554E-40, 3.65054388316495E-40, 1.3500719515971733E-40,
4.992663296375074E-41, 1.8462230151016932E-41, 6.82673496078692E-42,
2.5241745843582964E-42, 9.332621544394427E-43, 3.4503748253869174E-43,
1.2755799391582677E-43, 4.715505767071247E-44, 1.7431247815788975E-44,
6.443304575882822E-45, 2.3816021326359414E-45, 8.802590402549339E-46,
3.253365591819818E-46, 1.2023659923256258E-46, 4.443470256622595E-47,
1.642063616040405E-47, 6.067923645149806E-48, 2.242192916317133E-48,
8.28493086617883E-49, 3.0611748798826205E-49, 1.1310219160032624E-49,
4.1786674309371387E-50, 1.5437919701070235E-50, 5.703273458254764E-51,
2.1069019326139996E-51, 7.783042630583783E-52, 2.8750122923531045E-52,
1.0619778352758944E-52, 3.92262583825212E-53, 1.448852691202891E-53,
5.35128052444453E-54, 1.976412365871322E-54, 7.299346475035788E-55,
2.6957348629103867E-55, 9.955368869046248E-56, 3.6764166212254784E-56,
1.357623865663977E-56, 5.013277852874676E-57, 1.8511936055520708E-57,
6.835493017610423E-58, 2.5239221838619777E-58, 9.319022620872382E-59,
3.440750878188984E-59, 1.2703538780218175E-59, 4.690131222298997E-60,
1.7315467316230336E-60, 6.392526592528376E-61, 2.3599354585800467E-61,
8.711984891369027E-62, 3.2160563784586906E-62, 1.1871890678637917E-62,
4.3823389020648563E-63, 1.6176404755011315E-63, 5.971014948465491E-64,
2.203964285160501E-64, 8.134883110057956E-65, 3.00253845654203E-65,
1.1081957030698288E-65, 4.090110949889288E-66, 1.5095398404108963E-66,
5.5711525566499164E-67, 2.0560639502671574E-67, 7.587861495197822E-68,
2.8002286266800977E-68, 1.0333775343358301E-68, 3.8134326059141613E-69,
1.4072289965098046E-69, 5.192843304271672E-70, 1.9161852988107498E-70,
7.070688953528846E-71, 2.6090235625867065E-71, 9.626898778288872E-72,
3.552115173719195E-72, 1.310629665354808E-72, 4.835766925193849E-73,
1.7841987392899996E-73, 6.5828458296778185E-74, 2.4287164478649817E-74,
8.960509453689453E-75, 3.3058369695520314E-75, 1.219615920133218E-75,
4.4994324937839855E-76, 1.6599135790494626E-76, 6.123593794681787E-77,
2.2590222381394057E-77, 8.333509759043794E-78, 3.074176330186928E-78,
1.134026117237392E-78, 4.1832215550469095E-79, 1.543093614958859E-79,
5.6920322031473625E-80, 2.0995980329558125E-80, 7.74459599464029E-81,
2.85663823156747E-81, 1.0536725180347681E-81, 3.886422903593632E-82,
1.4334697044076614E-82, 5.287143530529042E-83, 1.950059498431422E-83,
7.192317727065618E-84, 2.65267574061198E-84, 9.783490415232438E-85,
3.6082606004887183E-85, 1.3307498896896168E-85, 4.907829957616485E-86,
1.8099940863785973E-86, 6.675125551376249E-87, 2.4617072711875098E-87,
9.078374983370637E-88, 3.347916493100407E-88, 1.234627532534424E-88,
4.552942751534994E-89, 1.678971645534168E-89, 6.191409411320687E-90,
2.2831304214168704E-90, 8.419125868841565E-91, 3.1045482525179376E-91,
1.1447878254335282E-91, 4.221305768426715E-92, 1.5565529039073294E-92,
5.739529209325878E-93, 2.116333018437095E-93, 7.803458780210424E-94,
2.8773037473193345E-94, 1.0609129759005472E-94, 3.9117341490226224E-95,
1.4422962097694012E-95, 5.317839108309043E-96, 1.9607018779274333E-96,
7.22908971850219E-97, 2.665332418100625E-97, 9.826863041459781E-98,
3.6230490775885395E-98, 1.3357628464105089E-98, 4.9247066662060635E-99,
1.815629691496367E-99, 6.693759762227838E-100, 2.4677941944612486E-100,
9.097954503505493E-101, 3.354089314819583E-101, 1.2365211780507929E-101,
4.558527735610431E-102, 1.6805203797151638E-102, 6.195254293288485E-103,
2.2838663978785864E-103, 8.419348522182607E-104, 3.1037200869608617E-104,
1.1441498379275434E-104, 4.217737893326829E-105, 1.5547931880030164E-105,
5.731417668555093E-106, 2.1127489947052878E-106, 7.788076248501916E-107,
2.8708398098537606E-107, 1.0582401321138131E-107, 3.900821029109967E-108,
1.4378857535340492E-108, 5.300164164498457E-109, 1.9536684993215962E-109,
7.201269517448687E-110, 2.6543850425066388E-110, 9.783978021815971E-111,
3.60631561265367E-111, 1.3292563314363002E-111, 4.899485877771569E-112,
1.805880775156918E-112, 6.656170767839589E-113, 2.4533339961654032E-113,
9.042442649659337E-114, 3.3328190678497415E-114, 1.22838533593063E-114,
4.527458146450243E-115, 1.6686729384675867E-115, 6.150139715423595E-116,
2.2667088686152173E-116, 8.354174273010658E-117, 3.0789906094480836E-117,
1.1347763252739112E-117, 4.182242959542931E-118, 1.5413645345559194E-118,
5.680657394985829E-119, 2.093577227165486E-119, 7.71572156818343E-120,
2.843552836391193E-120, 1.0479565372989326E-120, 3.862090742327071E-121,
1.4233081270792701E-121, 5.245328081886801E-122, 1.9330526255616414E-122,
7.123804760336865E-123, 2.6252923258053868E-123, 9.674770640969486E-124,
3.5653408388611528E-124, 1.3138894170724144E-124, 4.841879516283081E-125,
1.7842946084267174E-125, 6.575315354694613E-126, 2.4230594299528936E-126,
8.929127392607606E-127, 3.2904209581109934E-127, 1.2125269883744376E-127,
4.4681617809198583E-128, 1.6465082305680545E-128, 6.067314120190377E-129,
2.2357674438264404E-129, 8.238617985357862E-130, 3.0358458044850794E-130,
1.1186717595775604E-130, 4.122145123202054E-131, 1.5189431966300763E-131,
5.597027708260043E-132, 2.0623912822454442E-132, 7.59945364258836E-133,
2.8002150755576136E-133, 1.0318062088081041E-133, 3.801916962102057E-134,
1.4008927083586238E-134, 5.161844223336349E-135, 1.9019657812800466E-135,
7.008067952345277E-136, 2.5822109780452196E-136, 9.514434261294118E-137,
3.5056782786114466E-137, 1.2916920915390393E-137, 4.759308006368524E-138,
1.7535837187627455E-138, 6.461109364087714E-139, 2.3805953866166293E-139,
8.771260387335206E-140, 3.231739553480377E-140, 1.1907172619663312E-140,
4.387113846532066E-141, 1.6163936318765144E-141, 5.955432966052453E-142,
2.194206747972747E-142, 8.0842504604286535E-143, 2.9785163867542895E-143,
1.0973830791834613E-143, 4.043100797695488E-144, 1.4895974957821232E-144,
5.4880917749400255E-145, 2.021956814115056E-145, 7.449385951137214E-146,
2.744524909164136E-146, 1.011141627837606E-146, 3.725245743920311E-147,
1.372448315515089E-147, 5.0563277269788585E-148, 1.8628272732268658E-148,
6.862907211033043E-149, 2.5283770749683204E-149, 9.314804590706179E-150,
3.43165683125897E-150, 1.2642476735657576E-150, 4.657562523756934E-151,
1.7158663336325764E-151, 6.32130072853293E-152, 2.328775908734385E-152,
8.579208181428029E-153, 3.160566860525624E-153, 1.1643432673495662E-153,
4.2893884563898343E-154, 1.5801853054648018E-154, 5.821286397060868E-155,
2.1445107642460705E-155, 7.900158792929099E-156, 2.910326716837242E-156,
1.0721264639159536E-156, 3.949559306679236E-157, 1.454955215618703E-157,
5.359804923821673E-158, 1.9744526417782592E-158, 7.273490044990256E-159,
2.679398895963787E-159, 9.870298935204573E-160, 3.6359814031264123E-160,
1.3394034840219444E-160, 4.934005480046516E-161, 1.8175496000226797E-161,
6.695320391881055E-162, 2.4663511743923165E-162, 9.085250989153236E-163,
3.346704892346962E-163, 1.2328106829603568E-163, 4.541233670739288E-164,
1.6728222906579478E-164, 6.162036549324029E-165, 2.2698504590841385E-165,
8.361202532477945E-166, 3.0799152720361095E-166, 1.1345073959739233E-166,
4.179019882986145E-167, 1.539359636673928E-167, 5.670277466730624E-168,
2.088656843839046E-168, 7.693579561657345E-169, 2.8339252196987657E-169,
1.0438711885074012E-169, 3.8450680561815985E-170, 1.4163146698618915E-170,
5.216918899439567E-171, 1.921617891683406E-171, 7.078131150113517E-172,
2.607166741497965E-172, 9.603236989240896E-173, 3.5372449045585916E-173,
1.302900538180764E-173, 4.799058618646969E-174, 1.7676631565212049E-174,
6.510909157444578E-175, 2.3981836574604653E-175, 8.833278555507694E-176,
3.253569611755524E-176, 1.1983866898501176E-176, 4.414002282734271E-177,
1.625798952221923E-177, 5.988248864023971E-178, 2.2056244326339168E-178,
8.123852230162674E-179, 2.9922040084441497E-179, 1.1020952378674422E-179,
4.059249982288677E-180, 1.4951033731900594E-180, 5.506750544095511E-181,
2.028235358386933E-181, 7.470334791908225E-182, 2.751443236970537E-182,
1.01339750609043E-182, 3.7324834618661153E-183, 1.3747215671187737E-183,
5.063261978414284E-184, 1.864854083675686E-184, 6.868440224565288E-185,
2.5297069012654846E-185, 9.317107559677908E-186, 3.431554014119103E-186,
1.2638612579490836E-186, 4.654861897808366E-187, 1.7144034926118035E-187,
6.314197077798479E-188, 2.3255306957391882E-188, 8.564950913649061E-189,
3.154471036011241E-189, 1.1617885149300188E-189, 4.278844341376186E-190,
1.575885921907071E-190, 5.803927090641182E-191, 2.1375585033017118E-191,
7.87250537645065E-192, 2.8993912322160236E-192, 1.0678237620125319E-192,
3.9327039070615877E-193, 1.4483777096250334E-193, 5.334225044897781E-194,
1.964534949885737E-194, 7.235142901749931E-195, 2.664608447241313E-195,
9.813380205359242E-196, 3.614121867783531E-196, 1.3310240650353442E-196,
4.901939809651039E-197, 1.805298400757316E-197, 6.6485812667658745E-198,
2.4485441051546535E-198, 9.017493772133676E-199, 3.320953122904365E-199,
1.2230342481580683E-199, 4.504156533802387E-200, 1.6587743823710399E-200,
6.1088597730320516E-201, 2.2497382338537297E-201, 8.285196687399061E-202,
3.051213952399593E-202, 1.1236770467290546E-202, 4.138179940548775E-203,
1.5239694124595614E-203, 5.612316459435819E-204, 2.066841049510087E-204,
7.611513738834091E-205, 2.803070539473724E-205, 1.0322765377523142E-205,
3.8015186757326975E-206, 1.399965220364452E-206, 5.155566760200242E-207,
1.8986050759965313E-207, 6.9918470363895875E-208, 2.5748282264084767E-208,
9.482081204794577E-209, 3.491870660898607E-209, 1.2859133392450805E-209,
4.735483539242872E-210, 1.743877773704403E-210, 6.421948798953034E-211,
2.364921468008692E-211, 8.708948075247807E-212, 3.207109401492489E-212,
1.181030051687586E-212, 4.349178432038412E-213, 1.6015946520850137E-213,
5.897895930514445E-214, 2.17190444702779E-214, 7.998037608511487E-215,
2.945271200813268E-215, 1.0845916703268429E-215, 3.9939844265475195E-216
};
//@formatter:on
private static final int LF_LENGTH = longFactorials.length;
private static final int DF_LENGTH = doubleFactorials.length;
private static final int LENGTH_SUM = LF_LENGTH + DF_LENGTH;
// Constants for Stirling approximation in logFactorial() and
// stirlingCorrection()
private static final double C0 = 9.18938533204672742e-01;
private static final double C1 = 8.33333333333333333e-02; // +1/12
private static final double C3 = -2.77777777777777778e-03; // -1/360
private static final double C5 = 7.93650793650793651e-04; // +1/1260
private static final double C7 = -5.95238095238095238e-04; // -1/1680
/**
* Efficiently returns the binomial coefficient, often also referred to as
* "n over k" or "n choose k".
*
* The binomial coefficient is defined as
* (n * n-1 * ... * n-k+1 ) / ( 1 * 2 * ... * k ).
*
* - k<0: 0.
*
- k==0: 1.
*
- k==1: n.
*
- else: (n * n-1 * ... * n-k+1 ) / ( 1 * 2 * ... * k ).
*
*
* @return the binomial coefficient.
*/
public static double binomial(final double n, final long k) {
if (k < 0) {
return 0;
}
if (k == 0) {
return 1;
}
if (k == 1) {
return n;
}
// binomial(n,k) = (n * n-1 * ... * n-k+1 ) / ( 1 * 2 * ... * k )
double a = n - k + 1;
double b = 1;
double binomial = 1;
for (long i = k; i-- > 0; /**/) {
binomial *= (a++) / (b++);
}
return binomial;
}
/**
* Efficiently returns the binomial coefficient, often also referred to as
* "n over k" or "n choose k".
*
* The binomial coefficient is defined as
*
* - k<0: 0.
*
- k==0 || k==n: 1.
*
- k==1 || k==n-1: n.
*
- else: (n * n-1 * ... * n-k+1 ) / ( 1 * 2 * ... * k ).
*
*
* @return the binomial coefficient.
*/
public static double binomial(final long n, long k) {
if (k < 0) {
return 0;
}
if (k == 0 || k == n) {
return 1;
}
if (k == 1 || k == n - 1) {
return n;
}
// try quick version and see whether we get numeric overflows.
// factorial(..) is O(1); requires no loop; only a table lookup.
if (n > k) {
int max = longFactorials.length + doubleFactorials.length;
if (n < max) { // if (n! < inf && k! < inf)
final double n_fac = factorial((int) n);
final double k_fac = factorial((int) k);
final double n_minus_k_fac = factorial((int) (n - k));
final double nk = n_minus_k_fac * k_fac;
if (nk != Double.POSITIVE_INFINITY) { // no numeric overflow?
// now this is completely safe and accurate
return n_fac / nk;
}
}
if (k > n / 2) {
k = n - k; // quicker
}
}
// binomial(n,k) = (n * n-1 * ... * n-k+1 ) / ( 1 * 2 * ... * k )
long a = n - k + 1;
long b = 1;
double binomial = 1;
for (long i = k; i-- > 0; /**/) {
binomial *= ((double) (a++)) / (b++);
}
return binomial;
}
/**
* Returns the smallest long >= value.
*
* Examples: 1.0 -> 1, 1.2 -> 2, 1.9 -> 2.
*
* This method is safer than using (long) Math.ceil(value), because of
* possible rounding error.
*/
public static long ceil(final double value) {
return Math.round(Math.ceil(value));
}
/**
* Instantly returns the factorial k!.
* If {@code k >= 171} {@link Double#POSITIVE_INFINITY} is returned.
*
* @param k an integer >= 0
* @return the factorial of k
*/
public static double factorial(final int k) {
if (k < 0) {
throw new IllegalArgumentException("k < 0: (k = " + k + ")");
}
if (k < LF_LENGTH) {
return longFactorials[k];
}
if (k < LENGTH_SUM) {
return doubleFactorials[k - LF_LENGTH];
}
return Double.POSITIVE_INFINITY;
}
/**
* Returns the largest long <= value.
*
* Examples:
*
*
* 1.0 -> 1, 1.2 -> 1, 1.9 -> 1
* 2.0 -> 2, 2.2 -> 2, 2.9 -> 2
*
*
* This method is safer than using (long) Math.floor(value), because of
* possible rounding error.
*/
public static long floor(final double value) {
return Math.round(Math.floor(value));
}
/**
* Returns logbase(x).
*/
public static double log(final double base, final double x) {
return Math.log(x) / Math.log(base);
}
/**
* Returns log2(x).
*/
public static double log2(final double x) {
// 1.0 / Math.log(2) == 1.4426950408889634
return Math.log(x) * 1.4426950408889634;
}
/**
* Returns log(k!).
*
* Tries to avoid overflows. For k < 30 simply looks up a table in
* O(1). For k >= 30 uses Stirlings approximation.
*
* @param k
* must hold k >= 0.
*/
public static double logFactorial(final int k) {
if (k >= 30) {
final double r = 1.0 / (double) k;
final double rr = r * r;
return (k + 0.5) * Math.log(k) - k + C0 + r
* (C1 + rr * (C3 + rr * (C5 + rr * C7)));
} else {
return logFactorials[k];
}
}
/**
* Returns the value of {@code n!/n^n}.
* @param n an integer >= 0
* @return the value of {@code n!/n^n}
*/
public static double factoPow(int n) {
if (n < 0) {
throw new IllegalArgumentException("n < 0 : " + n);
}
if (n < factoPows.length) {
return factoPows[n];
}
double x = 1.0 / n;
for (int i = 2; i <= n; i++) {
x *= (double) i / n;
}
return x;
}
/**
* Instantly returns the factorial k!.
*
* @param k
* must hold k >= 0 && k < 21.
*/
public static long longFactorial(final int k) {
if (k < 0) {
throw new IllegalArgumentException("Negative k = " + k);
}
if (k < longFactorials.length) {
return longFactorials[k];
}
throw new IllegalArgumentException("Overflow for k = " + k);
}
/**
* Returns the StirlingCorrection.
*
* Correction term of the Stirling approximation for log(k!)
* (series in 1/k, or table values for small k) with int parameter k.
*
*
*
* log k! = (k + 1/2)log(k + 1) - (k + 1) + (1/2)log(2Pi) +
* stirlingCorrection(k + 1)
*
* log k! = (k + 1/2)log(k) - k + (1/2)log(2Pi) +
* stirlingCorrection(k)
*
*/
public static double stirlingCorrection(final int k) {
if (k > 30) {
final double r = 1.0 / (double) k;
final double rr = r * r;
return r * (C1 + rr * (C3 + rr * (C5 + rr * C7)));
} else {
return stirlingCorrection[k];
}
}
/**
* Checks whether the passed {@code value} is one of {@link Double#NaN},
* {@link Double#NEGATIVE_INFINITY} or {@link Double#POSITIVE_INFINITY}.
* Returns {@code true} if this is the case, otherwise returns
* {@code false}.
*
* @param value
* the double to check
* @return {@code true} when {@code value} is finite and not
* {@link Double#NaN}, {@code false} in all other cases
*/
public static boolean isBadNum(double value) {
if (Double.isNaN(value) || Double.isInfinite(value)) {
return true;
}
return false;
}
/**
* Checks whether the passed {@code value} qualifies as a probability, i.e.,
* whether it fulfills {@code 0.0 <= value <= 1.0}.
*
* @param value
* the double to check
* @return {@code true} when {@code value} lies within {@code [0, 1]},
* otherwise {@code false}
*/
public static boolean isProbability(double value) {
return value >= 0.0 && value <= 1.0;
}
public static int roundToIntStochastically(double x) {
if (isBadNum(x)) {
return 0;
}
double fraction = x % 1.0;
if (Math.abs(fraction) == 0.5) {
// stochastic rounding if the fractional part is exactly 0.5
if (DefaultRng.getGlobalPseudoRandom().nextDouble() >= 0.5) {
x += 0.5;
} else {
x -= 0.5;
}
}
return (int) Math.rint(x);
}
public static double round(double x, int scale) {
if (isBadNum(x)) {
return x;
}
return BigDecimal.valueOf(x).setScale(scale, RoundingMode.HALF_EVEN).doubleValue();
}
public static double round(double x) {
return round(x, 3);
}
private Arithmetic() {
}
}