You are viewing a plain text version of this content. The canonical link for it is here.
Posted to commits@mahout.apache.org by sr...@apache.org on 2009/11/25 04:31:49 UTC
svn commit: r883972 [1/10] - in /lucene/mahout/trunk:
core/src/main/java/org/apache/mahout/fpm/pfpgrowth/
matrix/src/main/java/org/apache/mahout/jet/math/
matrix/src/main/java/org/apache/mahout/jet/random/
matrix/src/main/java/org/apache/mahout/jet/ran...
Author: srowen
Date: Wed Nov 25 03:31:47 2009
New Revision: 883972
URL: http://svn.apache.org/viewvc?rev=883972&view=rev
Log:
MAHOUT-204 -- Drew's author-cleanup patch (in parts)
Modified:
lucene/mahout/trunk/core/src/main/java/org/apache/mahout/fpm/pfpgrowth/package.html
lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/math/Arithmetic.java
lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/math/Bessel.java
lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/math/Constants.java
lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/math/Functions.java
lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/math/IntFunctions.java
lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/math/Mult.java
lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/math/PlusMult.java
lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/math/Polynomial.java
lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/AbstractContinousDistribution.java
lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/AbstractDiscreteDistribution.java
lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/AbstractDistribution.java
lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/Benchmark.java
lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/Beta.java
lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/Binomial.java
lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/BreitWigner.java
lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/BreitWignerMeanSquare.java
lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/ChiSquare.java
lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/Distributions.java
lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/Empirical.java
lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/EmpiricalWalker.java
lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/Exponential.java
lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/ExponentialPower.java
lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/Fun.java
lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/Gamma.java
lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/HyperGeometric.java
lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/Hyperbolic.java
lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/Logarithmic.java
lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/NegativeBinomial.java
lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/Normal.java
lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/Poisson.java
lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/PoissonSlow.java
lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/Stack.java
lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/StudentT.java
lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/Uniform.java
lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/VonMises.java
lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/Zeta.java
lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/engine/Benchmark.java
lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/engine/DRand.java
lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/engine/MersenneTwister.java
lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/engine/MersenneTwister64.java
lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/engine/RandomEngine.java
lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/engine/RandomGenerator.java
lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/engine/RandomSeedGenerator.java
lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/engine/RandomSeedTable.java
lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/package.html
lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/sampling/RandomSampler.java
lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/sampling/RandomSamplingAssistant.java
lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/random/sampling/WeightedRandomSampler.java
lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/stat/Descriptive.java
lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/stat/Gamma.java
lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/stat/Probability.java
lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/stat/quantile/Buffer.java
lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/stat/quantile/DoubleBuffer.java
lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/stat/quantile/DoubleBufferSet.java
lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/stat/quantile/DoubleQuantileEstimator.java
lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/stat/quantile/EquiDepthHistogram.java
lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/stat/quantile/ExactDoubleQuantileFinder.java
lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/stat/quantile/KnownDoubleQuantileEstimator.java
lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/stat/quantile/Quantile1Test.java
lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/stat/quantile/QuantileCalc.java
lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/stat/quantile/QuantileFinderFactory.java
lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/stat/quantile/QuantileFinderTest.java
lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/stat/quantile/UnknownDoubleQuantileEstimator.java
lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/stat/quantile/Utils.java
lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/stat/quantile/package.html
Modified: lucene/mahout/trunk/core/src/main/java/org/apache/mahout/fpm/pfpgrowth/package.html
URL: http://svn.apache.org/viewvc/lucene/mahout/trunk/core/src/main/java/org/apache/mahout/fpm/pfpgrowth/package.html?rev=883972&r1=883971&r2=883972&view=diff
==============================================================================
--- lucene/mahout/trunk/core/src/main/java/org/apache/mahout/fpm/pfpgrowth/package.html (original)
+++ lucene/mahout/trunk/core/src/main/java/org/apache/mahout/fpm/pfpgrowth/package.html Wed Nov 25 03:31:47 2009
@@ -39,10 +39,10 @@
FPGrowth<String> fp = new FPGrowth<String>();
Set<String> features = new HashSet<String>();
fp.generateTopKStringFrequentPatterns(
- new StringRecordIterator(new FileLineIterable(new File(input), encoding, false), pattern),
+ new StringRecordIterator(new FileLineIterable(new File(input), encoding, false), pattern),
fp.generateFList(
- new StringRecordIterator(new FileLineIterable(new File(input), encoding, false), pattern), minSupport),
- minSupport,
+ new StringRecordIterator(new FileLineIterable(new File(input), encoding, false), pattern), minSupport),
+ minSupport,
maxHeapSize,
features,
new StringOutputConvertor(new SequenceFileOutputCollector<Text, TopKStringPatterns>(writer))
Modified: lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/math/Arithmetic.java
URL: http://svn.apache.org/viewvc/lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/math/Arithmetic.java?rev=883972&r1=883971&r2=883972&view=diff
==============================================================================
--- lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/math/Arithmetic.java (original)
+++ lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/math/Arithmetic.java Wed Nov 25 03:31:47 2009
@@ -16,220 +16,220 @@
*/
@Deprecated
public class Arithmetic extends Constants {
- // 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 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
- protected 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
- };
+ // for method logFactorial(...)
+ // log(k!) for k = 0, ..., 29
+ protected 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
- protected 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 = 0, ..., 20
+ protected 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
- protected 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
- };
-
+ // k! for k = 21, ..., 170
+ protected 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
+ };
+
/**
* Makes this class non instantiable, but still let's others inherit from it.
*/
@@ -246,18 +246,18 @@
* @return the binomial coefficient.
*/
public static double binomial(double n, 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;
+ 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".
@@ -271,35 +271,35 @@
* @return the binomial coefficient.
*/
public static double binomial(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)
- double n_fac = factorial((int)n);
- double k_fac = factorial((int)k);
- double n_minus_k_fac = factorial((int) (n-k));
- 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;
+ 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)
+ double n_fac = factorial((int)n);
+ double k_fac = factorial((int)k);
+ double n_minus_k_fac = factorial((int) (n-k));
+ 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 <code>long >= value</code>.
@@ -307,7 +307,7 @@
* This method is safer than using (long) Math.ceil(value), because of possible rounding error.
*/
public static long ceil(double value) {
- return Math.round(Math.ceil(value));
+ return Math.round(Math.ceil(value));
}
/**
* Evaluates the series of Chebyshev polynomials Ti at argument x/2.
@@ -345,63 +345,63 @@
* @param N the number of coefficients.
*/
public static double chbevl( double x, double coef[], int N ) throws ArithmeticException {
- double b0, b1, b2;
+ double b0, b1, b2;
- int p = 0;
- int i;
+ int p = 0;
+ int i;
- b0 = coef[p++];
- b1 = 0.0;
- i = N - 1;
+ b0 = coef[p++];
+ b1 = 0.0;
+ i = N - 1;
- do {
- b2 = b1;
- b1 = b0;
- b0 = x * b1 - b2 + coef[p++];
- } while( --i > 0);
+ do {
+ b2 = b1;
+ b1 = b0;
+ b0 = x * b1 - b2 + coef[p++];
+ } while( --i > 0);
- return( 0.5*(b0-b2) );
+ return( 0.5*(b0-b2) );
}
/**
* Returns the factorial of the argument.
*/
static private long fac1(int j) {
- long i = j;
- if(j < 0) i = Math.abs(j);
- if (i>longFactorials.length) throw new IllegalArgumentException("Overflow");
-
- long d = 1;
- while (i > 1) d *= i--;
-
- if (j < 0) return -d;
- else return d;
+ long i = j;
+ if(j < 0) i = Math.abs(j);
+ if (i>longFactorials.length) throw new IllegalArgumentException("Overflow");
+
+ long d = 1;
+ while (i > 1) d *= i--;
+
+ if (j < 0) return -d;
+ else return d;
}
/**
* Returns the factorial of the argument.
*/
static private double fac2(int j) {
- long i = j;
- if (j < 0) i = Math.abs(j);
-
- double d = 1.0;
- while (i > 1) d *= i--;
-
- if (j < 0) return -d;
- else return d;
+ long i = j;
+ if (j < 0) i = Math.abs(j);
+
+ double d = 1.0;
+ while (i > 1) d *= i--;
+
+ if (j < 0) return -d;
+ else return d;
}
/**
* Instantly returns the factorial <tt>k!</tt>.
* @param k must hold <tt>k >= 0</tt>.
*/
static public double factorial(int k) {
- if (k<0) throw new IllegalArgumentException();
-
- int length1 = longFactorials.length;
- if (k<length1) return longFactorials[k];
-
- int length2 = doubleFactorials.length;
- if (k<length1+length2) return doubleFactorials[k-length1];
- else return Double.POSITIVE_INFINITY;
+ if (k<0) throw new IllegalArgumentException();
+
+ int length1 = longFactorials.length;
+ if (k<length1) return longFactorials[k];
+
+ int length2 = doubleFactorials.length;
+ if (k<length1+length2) return doubleFactorials[k-length1];
+ else return Double.POSITIVE_INFINITY;
}
/**
* Returns the largest <code>long <= value</code>.
@@ -411,27 +411,27 @@
* This method is safer than using (long) Math.floor(value), because of possible rounding error.
*/
public static long floor(double value) {
- return Math.round(Math.floor(value));
+ return Math.round(Math.floor(value));
}
/**
* Returns <tt>log<sub>base</sub>value</tt>.
*/
public static double log(double base, double value) {
- return Math.log(value) / Math.log(base);
+ return Math.log(value) / Math.log(base);
}
/**
* Returns <tt>log<sub>10</sub>value</tt>.
*/
static public double log10(double value) {
- // 1.0 / Math.log(10) == 0.43429448190325176
- return Math.log(value) * 0.43429448190325176;
+ // 1.0 / Math.log(10) == 0.43429448190325176
+ return Math.log(value) * 0.43429448190325176;
}
/**
* Returns <tt>log<sub>2</sub>value</tt>.
*/
static public double log2(double value) {
- // 1.0 / Math.log(2) == 1.4426950408889634
- return Math.log(value) * 1.4426950408889634;
+ // 1.0 / Math.log(2) == 1.4426950408889634
+ return Math.log(value) * 1.4426950408889634;
}
/**
* Returns <tt>log(k!)</tt>.
@@ -441,30 +441,30 @@
* @param k must hold <tt>k >= 0</tt>.
*/
public static double logFactorial(int k) {
- if (k >= 30) {
- double r, rr;
- final double C0 = 9.18938533204672742e-01;
- final double C1 = 8.33333333333333333e-02;
- final double C3 = -2.77777777777777778e-03;
- final double C5 = 7.93650793650793651e-04;
- final double C7 = -5.95238095238095238e-04;
+ if (k >= 30) {
+ double r, rr;
+ final double C0 = 9.18938533204672742e-01;
+ final double C1 = 8.33333333333333333e-02;
+ final double C3 = -2.77777777777777778e-03;
+ final double C5 = 7.93650793650793651e-04;
+ final double C7 = -5.95238095238095238e-04;
- r = 1.0 / (double) k;
- rr = r * r;
- return (k + 0.5)*Math.log(k) - k + C0 + r*(C1 + rr*(C3 + rr*(C5 + rr*C7)));
- }
- else
- return logFactorials[k];
+ r = 1.0 / (double) k;
+ rr = r * r;
+ return (k + 0.5)*Math.log(k) - k + C0 + r*(C1 + rr*(C3 + rr*(C5 + rr*C7)));
+ }
+ else
+ return logFactorials[k];
}
/**
* Instantly returns the factorial <tt>k!</tt>.
* @param k must hold <tt>k >= 0 && k < 21</tt>.
*/
static public long longFactorial(int k) throws IllegalArgumentException{
- if (k<0) throw new IllegalArgumentException("Negative k");
-
- if (k<longFactorials.length) return longFactorials[k];
- throw new IllegalArgumentException("Overflow");
+ if (k<0) throw new IllegalArgumentException("Negative k");
+
+ if (k<longFactorials.length) return longFactorials[k];
+ throw new IllegalArgumentException("Overflow");
}
/**
* Returns the StirlingCorrection.
@@ -482,24 +482,24 @@
* </tt>
*/
public static double stirlingCorrection(int k) {
- final double C1 = 8.33333333333333333e-02; // +1/12
- final double C3 = -2.77777777777777778e-03; // -1/360
- final double C5 = 7.93650793650793651e-04; // +1/1260
- final double C7 = -5.95238095238095238e-04; // -1/1680
+ final double C1 = 8.33333333333333333e-02; // +1/12
+ final double C3 = -2.77777777777777778e-03; // -1/360
+ final double C5 = 7.93650793650793651e-04; // +1/1260
+ final double C7 = -5.95238095238095238e-04; // -1/1680
- double r, rr;
+ double r, rr;
- if (k > 30) {
- r = 1.0 / (double) k;
- rr = r * r;
- return r*(C1 + rr*(C3 + rr*(C5 + rr*C7)));
- }
- else return stirlingCorrection[k];
+ if (k > 30) {
+ r = 1.0 / (double) k;
+ rr = r * r;
+ return r*(C1 + rr*(C3 + rr*(C5 + rr*C7)));
+ }
+ else return stirlingCorrection[k];
}
/**
* Equivalent to <tt>Math.round(binomial(n,k))</tt>.
*/
private static long xlongBinomial(long n, long k) {
- return Math.round(binomial(n,k));
+ return Math.round(binomial(n,k));
}
}
Modified: lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/math/Bessel.java
URL: http://svn.apache.org/viewvc/lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/math/Bessel.java?rev=883972&r1=883971&r2=883972&view=diff
==============================================================================
--- lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/math/Bessel.java (original)
+++ lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/math/Bessel.java Wed Nov 25 03:31:47 2009
@@ -16,277 +16,277 @@
*/
@Deprecated
public class Bessel extends Constants {
- /****************************************
- * COEFFICIENTS FOR METHODS i0, i0e *
- ****************************************/
-
- /**
- * Chebyshev coefficients for exp(-x) I0(x)
- * in the interval [0,8].
- *
- * lim(x->0){ exp(-x) I0(x) } = 1.
- */
- protected static final double[] A_i0 = {
- -4.41534164647933937950E-18,
- 3.33079451882223809783E-17,
- -2.43127984654795469359E-16,
- 1.71539128555513303061E-15,
- -1.16853328779934516808E-14,
- 7.67618549860493561688E-14,
- -4.85644678311192946090E-13,
- 2.95505266312963983461E-12,
- -1.72682629144155570723E-11,
- 9.67580903537323691224E-11,
- -5.18979560163526290666E-10,
- 2.65982372468238665035E-9,
- -1.30002500998624804212E-8,
- 6.04699502254191894932E-8,
- -2.67079385394061173391E-7,
- 1.11738753912010371815E-6,
- -4.41673835845875056359E-6,
- 1.64484480707288970893E-5,
- -5.75419501008210370398E-5,
- 1.88502885095841655729E-4,
- -5.76375574538582365885E-4,
- 1.63947561694133579842E-3,
- -4.32430999505057594430E-3,
- 1.05464603945949983183E-2,
- -2.37374148058994688156E-2,
- 4.93052842396707084878E-2,
- -9.49010970480476444210E-2,
- 1.71620901522208775349E-1,
- -3.04682672343198398683E-1,
- 6.76795274409476084995E-1
- };
-
-
-
- /**
- * Chebyshev coefficients for exp(-x) sqrt(x) I0(x)
- * in the inverted interval [8,infinity].
- *
- * lim(x->inf){ exp(-x) sqrt(x) I0(x) } = 1/sqrt(2pi).
- */
- protected static final double[] B_i0 = {
- -7.23318048787475395456E-18,
- -4.83050448594418207126E-18,
- 4.46562142029675999901E-17,
- 3.46122286769746109310E-17,
- -2.82762398051658348494E-16,
- -3.42548561967721913462E-16,
- 1.77256013305652638360E-15,
- 3.81168066935262242075E-15,
- -9.55484669882830764870E-15,
- -4.15056934728722208663E-14,
- 1.54008621752140982691E-14,
- 3.85277838274214270114E-13,
- 7.18012445138366623367E-13,
- -1.79417853150680611778E-12,
- -1.32158118404477131188E-11,
- -3.14991652796324136454E-11,
- 1.18891471078464383424E-11,
- 4.94060238822496958910E-10,
- 3.39623202570838634515E-9,
- 2.26666899049817806459E-8,
- 2.04891858946906374183E-7,
- 2.89137052083475648297E-6,
- 6.88975834691682398426E-5,
- 3.36911647825569408990E-3,
- 8.04490411014108831608E-1
- };
-
-
-
- /****************************************
- * COEFFICIENTS FOR METHODS i1, i1e *
- ****************************************/
- /**
- * Chebyshev coefficients for exp(-x) I1(x) / x
- * in the interval [0,8].
- *
- * lim(x->0){ exp(-x) I1(x) / x } = 1/2.
- */
- protected static final double[] A_i1 = {
- 2.77791411276104639959E-18,
- -2.11142121435816608115E-17,
- 1.55363195773620046921E-16,
- -1.10559694773538630805E-15,
- 7.60068429473540693410E-15,
- -5.04218550472791168711E-14,
- 3.22379336594557470981E-13,
- -1.98397439776494371520E-12,
- 1.17361862988909016308E-11,
- -6.66348972350202774223E-11,
- 3.62559028155211703701E-10,
- -1.88724975172282928790E-9,
- 9.38153738649577178388E-9,
- -4.44505912879632808065E-8,
- 2.00329475355213526229E-7,
- -8.56872026469545474066E-7,
- 3.47025130813767847674E-6,
- -1.32731636560394358279E-5,
- 4.78156510755005422638E-5,
- -1.61760815825896745588E-4,
- 5.12285956168575772895E-4,
- -1.51357245063125314899E-3,
- 4.15642294431288815669E-3,
- -1.05640848946261981558E-2,
- 2.47264490306265168283E-2,
- -5.29459812080949914269E-2,
- 1.02643658689847095384E-1,
- -1.76416518357834055153E-1,
- 2.52587186443633654823E-1
- };
-
- /*
- * Chebyshev coefficients for exp(-x) sqrt(x) I1(x)
- * in the inverted interval [8,infinity].
- *
- * lim(x->inf){ exp(-x) sqrt(x) I1(x) } = 1/sqrt(2pi).
- */
- protected static final double[] B_i1 = {
- 7.51729631084210481353E-18,
- 4.41434832307170791151E-18,
- -4.65030536848935832153E-17,
- -3.20952592199342395980E-17,
- 2.96262899764595013876E-16,
- 3.30820231092092828324E-16,
- -1.88035477551078244854E-15,
- -3.81440307243700780478E-15,
- 1.04202769841288027642E-14,
- 4.27244001671195135429E-14,
- -2.10154184277266431302E-14,
- -4.08355111109219731823E-13,
- -7.19855177624590851209E-13,
- 2.03562854414708950722E-12,
- 1.41258074366137813316E-11,
- 3.25260358301548823856E-11,
- -1.89749581235054123450E-11,
- -5.58974346219658380687E-10,
- -3.83538038596423702205E-9,
- -2.63146884688951950684E-8,
- -2.51223623787020892529E-7,
- -3.88256480887769039346E-6,
- -1.10588938762623716291E-4,
- -9.76109749136146840777E-3,
- 7.78576235018280120474E-1
- };
-
-
-
-
- /****************************************
- * COEFFICIENTS FOR METHODS k0, k0e *
- ****************************************/
- /* Chebyshev coefficients for K0(x) + log(x/2) I0(x)
- * in the interval [0,2]. The odd order coefficients are all
- * zero; only the even order coefficients are listed.
- *
- * lim(x->0){ K0(x) + log(x/2) I0(x) } = -EUL.
- */
- protected static final double[] A_k0 = {
- 1.37446543561352307156E-16,
- 4.25981614279661018399E-14,
- 1.03496952576338420167E-11,
- 1.90451637722020886025E-9,
- 2.53479107902614945675E-7,
- 2.28621210311945178607E-5,
- 1.26461541144692592338E-3,
- 3.59799365153615016266E-2,
- 3.44289899924628486886E-1,
- -5.35327393233902768720E-1
- };
-
- /* Chebyshev coefficients for exp(x) sqrt(x) K0(x)
- * in the inverted interval [2,infinity].
- *
- * lim(x->inf){ exp(x) sqrt(x) K0(x) } = sqrt(pi/2).
- */
- protected static final double[] B_k0 = {
- 5.30043377268626276149E-18,
- -1.64758043015242134646E-17,
- 5.21039150503902756861E-17,
- -1.67823109680541210385E-16,
- 5.51205597852431940784E-16,
- -1.84859337734377901440E-15,
- 6.34007647740507060557E-15,
- -2.22751332699166985548E-14,
- 8.03289077536357521100E-14,
- -2.98009692317273043925E-13,
- 1.14034058820847496303E-12,
- -4.51459788337394416547E-12,
- 1.85594911495471785253E-11,
- -7.95748924447710747776E-11,
- 3.57739728140030116597E-10,
- -1.69753450938905987466E-9,
- 8.57403401741422608519E-9,
- -4.66048989768794782956E-8,
- 2.76681363944501510342E-7,
- -1.83175552271911948767E-6,
- 1.39498137188764993662E-5,
- -1.28495495816278026384E-4,
- 1.56988388573005337491E-3,
- -3.14481013119645005427E-2,
- 2.44030308206595545468E0
- };
-
-
-
-
- /****************************************
- * COEFFICIENTS FOR METHODS k1, k1e *
- ****************************************/
- /* Chebyshev coefficients for x(K1(x) - log(x/2) I1(x))
- * in the interval [0,2].
- *
- * lim(x->0){ x(K1(x) - log(x/2) I1(x)) } = 1.
- */
- protected static final double[] A_k1 = {
- -7.02386347938628759343E-18,
- -2.42744985051936593393E-15,
- -6.66690169419932900609E-13,
- -1.41148839263352776110E-10,
- -2.21338763073472585583E-8,
- -2.43340614156596823496E-6,
- -1.73028895751305206302E-4,
- -6.97572385963986435018E-3,
- -1.22611180822657148235E-1,
- -3.53155960776544875667E-1,
- 1.52530022733894777053E0
- };
-
- /* Chebyshev coefficients for exp(x) sqrt(x) K1(x)
- * in the interval [2,infinity].
- *
- * lim(x->inf){ exp(x) sqrt(x) K1(x) } = sqrt(pi/2).
- */
- protected static final double[] B_k1 = {
- -5.75674448366501715755E-18,
- 1.79405087314755922667E-17,
- -5.68946255844285935196E-17,
- 1.83809354436663880070E-16,
- -6.05704724837331885336E-16,
- 2.03870316562433424052E-15,
- -7.01983709041831346144E-15,
- 2.47715442448130437068E-14,
- -8.97670518232499435011E-14,
- 3.34841966607842919884E-13,
- -1.28917396095102890680E-12,
- 5.13963967348173025100E-12,
- -2.12996783842756842877E-11,
- 9.21831518760500529508E-11,
- -4.19035475934189648750E-10,
- 2.01504975519703286596E-9,
- -1.03457624656780970260E-8,
- 5.74108412545004946722E-8,
- -3.50196060308781257119E-7,
- 2.40648494783721712015E-6,
- -1.93619797416608296024E-5,
- 1.95215518471351631108E-4,
- -2.85781685962277938680E-3,
- 1.03923736576817238437E-1,
- 2.72062619048444266945E0
- };
+ /****************************************
+ * COEFFICIENTS FOR METHODS i0, i0e *
+ ****************************************/
+
+ /**
+ * Chebyshev coefficients for exp(-x) I0(x)
+ * in the interval [0,8].
+ *
+ * lim(x->0){ exp(-x) I0(x) } = 1.
+ */
+ protected static final double[] A_i0 = {
+ -4.41534164647933937950E-18,
+ 3.33079451882223809783E-17,
+ -2.43127984654795469359E-16,
+ 1.71539128555513303061E-15,
+ -1.16853328779934516808E-14,
+ 7.67618549860493561688E-14,
+ -4.85644678311192946090E-13,
+ 2.95505266312963983461E-12,
+ -1.72682629144155570723E-11,
+ 9.67580903537323691224E-11,
+ -5.18979560163526290666E-10,
+ 2.65982372468238665035E-9,
+ -1.30002500998624804212E-8,
+ 6.04699502254191894932E-8,
+ -2.67079385394061173391E-7,
+ 1.11738753912010371815E-6,
+ -4.41673835845875056359E-6,
+ 1.64484480707288970893E-5,
+ -5.75419501008210370398E-5,
+ 1.88502885095841655729E-4,
+ -5.76375574538582365885E-4,
+ 1.63947561694133579842E-3,
+ -4.32430999505057594430E-3,
+ 1.05464603945949983183E-2,
+ -2.37374148058994688156E-2,
+ 4.93052842396707084878E-2,
+ -9.49010970480476444210E-2,
+ 1.71620901522208775349E-1,
+ -3.04682672343198398683E-1,
+ 6.76795274409476084995E-1
+ };
+
+
+
+ /**
+ * Chebyshev coefficients for exp(-x) sqrt(x) I0(x)
+ * in the inverted interval [8,infinity].
+ *
+ * lim(x->inf){ exp(-x) sqrt(x) I0(x) } = 1/sqrt(2pi).
+ */
+ protected static final double[] B_i0 = {
+ -7.23318048787475395456E-18,
+ -4.83050448594418207126E-18,
+ 4.46562142029675999901E-17,
+ 3.46122286769746109310E-17,
+ -2.82762398051658348494E-16,
+ -3.42548561967721913462E-16,
+ 1.77256013305652638360E-15,
+ 3.81168066935262242075E-15,
+ -9.55484669882830764870E-15,
+ -4.15056934728722208663E-14,
+ 1.54008621752140982691E-14,
+ 3.85277838274214270114E-13,
+ 7.18012445138366623367E-13,
+ -1.79417853150680611778E-12,
+ -1.32158118404477131188E-11,
+ -3.14991652796324136454E-11,
+ 1.18891471078464383424E-11,
+ 4.94060238822496958910E-10,
+ 3.39623202570838634515E-9,
+ 2.26666899049817806459E-8,
+ 2.04891858946906374183E-7,
+ 2.89137052083475648297E-6,
+ 6.88975834691682398426E-5,
+ 3.36911647825569408990E-3,
+ 8.04490411014108831608E-1
+ };
+
+
+
+ /****************************************
+ * COEFFICIENTS FOR METHODS i1, i1e *
+ ****************************************/
+ /**
+ * Chebyshev coefficients for exp(-x) I1(x) / x
+ * in the interval [0,8].
+ *
+ * lim(x->0){ exp(-x) I1(x) / x } = 1/2.
+ */
+ protected static final double[] A_i1 = {
+ 2.77791411276104639959E-18,
+ -2.11142121435816608115E-17,
+ 1.55363195773620046921E-16,
+ -1.10559694773538630805E-15,
+ 7.60068429473540693410E-15,
+ -5.04218550472791168711E-14,
+ 3.22379336594557470981E-13,
+ -1.98397439776494371520E-12,
+ 1.17361862988909016308E-11,
+ -6.66348972350202774223E-11,
+ 3.62559028155211703701E-10,
+ -1.88724975172282928790E-9,
+ 9.38153738649577178388E-9,
+ -4.44505912879632808065E-8,
+ 2.00329475355213526229E-7,
+ -8.56872026469545474066E-7,
+ 3.47025130813767847674E-6,
+ -1.32731636560394358279E-5,
+ 4.78156510755005422638E-5,
+ -1.61760815825896745588E-4,
+ 5.12285956168575772895E-4,
+ -1.51357245063125314899E-3,
+ 4.15642294431288815669E-3,
+ -1.05640848946261981558E-2,
+ 2.47264490306265168283E-2,
+ -5.29459812080949914269E-2,
+ 1.02643658689847095384E-1,
+ -1.76416518357834055153E-1,
+ 2.52587186443633654823E-1
+ };
+
+ /*
+ * Chebyshev coefficients for exp(-x) sqrt(x) I1(x)
+ * in the inverted interval [8,infinity].
+ *
+ * lim(x->inf){ exp(-x) sqrt(x) I1(x) } = 1/sqrt(2pi).
+ */
+ protected static final double[] B_i1 = {
+ 7.51729631084210481353E-18,
+ 4.41434832307170791151E-18,
+ -4.65030536848935832153E-17,
+ -3.20952592199342395980E-17,
+ 2.96262899764595013876E-16,
+ 3.30820231092092828324E-16,
+ -1.88035477551078244854E-15,
+ -3.81440307243700780478E-15,
+ 1.04202769841288027642E-14,
+ 4.27244001671195135429E-14,
+ -2.10154184277266431302E-14,
+ -4.08355111109219731823E-13,
+ -7.19855177624590851209E-13,
+ 2.03562854414708950722E-12,
+ 1.41258074366137813316E-11,
+ 3.25260358301548823856E-11,
+ -1.89749581235054123450E-11,
+ -5.58974346219658380687E-10,
+ -3.83538038596423702205E-9,
+ -2.63146884688951950684E-8,
+ -2.51223623787020892529E-7,
+ -3.88256480887769039346E-6,
+ -1.10588938762623716291E-4,
+ -9.76109749136146840777E-3,
+ 7.78576235018280120474E-1
+ };
+
+
+
+
+ /****************************************
+ * COEFFICIENTS FOR METHODS k0, k0e *
+ ****************************************/
+ /* Chebyshev coefficients for K0(x) + log(x/2) I0(x)
+ * in the interval [0,2]. The odd order coefficients are all
+ * zero; only the even order coefficients are listed.
+ *
+ * lim(x->0){ K0(x) + log(x/2) I0(x) } = -EUL.
+ */
+ protected static final double[] A_k0 = {
+ 1.37446543561352307156E-16,
+ 4.25981614279661018399E-14,
+ 1.03496952576338420167E-11,
+ 1.90451637722020886025E-9,
+ 2.53479107902614945675E-7,
+ 2.28621210311945178607E-5,
+ 1.26461541144692592338E-3,
+ 3.59799365153615016266E-2,
+ 3.44289899924628486886E-1,
+ -5.35327393233902768720E-1
+ };
+
+ /* Chebyshev coefficients for exp(x) sqrt(x) K0(x)
+ * in the inverted interval [2,infinity].
+ *
+ * lim(x->inf){ exp(x) sqrt(x) K0(x) } = sqrt(pi/2).
+ */
+ protected static final double[] B_k0 = {
+ 5.30043377268626276149E-18,
+ -1.64758043015242134646E-17,
+ 5.21039150503902756861E-17,
+ -1.67823109680541210385E-16,
+ 5.51205597852431940784E-16,
+ -1.84859337734377901440E-15,
+ 6.34007647740507060557E-15,
+ -2.22751332699166985548E-14,
+ 8.03289077536357521100E-14,
+ -2.98009692317273043925E-13,
+ 1.14034058820847496303E-12,
+ -4.51459788337394416547E-12,
+ 1.85594911495471785253E-11,
+ -7.95748924447710747776E-11,
+ 3.57739728140030116597E-10,
+ -1.69753450938905987466E-9,
+ 8.57403401741422608519E-9,
+ -4.66048989768794782956E-8,
+ 2.76681363944501510342E-7,
+ -1.83175552271911948767E-6,
+ 1.39498137188764993662E-5,
+ -1.28495495816278026384E-4,
+ 1.56988388573005337491E-3,
+ -3.14481013119645005427E-2,
+ 2.44030308206595545468E0
+ };
+
+
+
+
+ /****************************************
+ * COEFFICIENTS FOR METHODS k1, k1e *
+ ****************************************/
+ /* Chebyshev coefficients for x(K1(x) - log(x/2) I1(x))
+ * in the interval [0,2].
+ *
+ * lim(x->0){ x(K1(x) - log(x/2) I1(x)) } = 1.
+ */
+ protected static final double[] A_k1 = {
+ -7.02386347938628759343E-18,
+ -2.42744985051936593393E-15,
+ -6.66690169419932900609E-13,
+ -1.41148839263352776110E-10,
+ -2.21338763073472585583E-8,
+ -2.43340614156596823496E-6,
+ -1.73028895751305206302E-4,
+ -6.97572385963986435018E-3,
+ -1.22611180822657148235E-1,
+ -3.53155960776544875667E-1,
+ 1.52530022733894777053E0
+ };
+
+ /* Chebyshev coefficients for exp(x) sqrt(x) K1(x)
+ * in the interval [2,infinity].
+ *
+ * lim(x->inf){ exp(x) sqrt(x) K1(x) } = sqrt(pi/2).
+ */
+ protected static final double[] B_k1 = {
+ -5.75674448366501715755E-18,
+ 1.79405087314755922667E-17,
+ -5.68946255844285935196E-17,
+ 1.83809354436663880070E-16,
+ -6.05704724837331885336E-16,
+ 2.03870316562433424052E-15,
+ -7.01983709041831346144E-15,
+ 2.47715442448130437068E-14,
+ -8.97670518232499435011E-14,
+ 3.34841966607842919884E-13,
+ -1.28917396095102890680E-12,
+ 5.13963967348173025100E-12,
+ -2.12996783842756842877E-11,
+ 9.21831518760500529508E-11,
+ -4.19035475934189648750E-10,
+ 2.01504975519703286596E-9,
+ -1.03457624656780970260E-8,
+ 5.74108412545004946722E-8,
+ -3.50196060308781257119E-7,
+ 2.40648494783721712015E-6,
+ -1.93619797416608296024E-5,
+ 1.95215518471351631108E-4,
+ -2.85781685962277938680E-3,
+ 1.03923736576817238437E-1,
+ 2.72062619048444266945E0
+ };
/**
* Makes this class non instantiable, but still let's others inherit from it.
@@ -305,14 +305,14 @@
* @param x the value to compute the bessel function of.
*/
static public double i0(double x) throws ArithmeticException {
- double y;
- if( x < 0 ) x = -x;
- if( x <= 8.0 ) {
- y = (x/2.0) - 2.0;
- return( Math.exp(x) * Arithmetic.chbevl( y, A_i0, 30 ) );
- }
+ double y;
+ if( x < 0 ) x = -x;
+ if( x <= 8.0 ) {
+ y = (x/2.0) - 2.0;
+ return( Math.exp(x) * Arithmetic.chbevl( y, A_i0, 30 ) );
+ }
- return( Math.exp(x) * Arithmetic.chbevl( 32.0/x - 2.0, B_i0, 25 ) / Math.sqrt(x) );
+ return( Math.exp(x) * Arithmetic.chbevl( 32.0/x - 2.0, B_i0, 25 ) / Math.sqrt(x) );
}
/**
* Returns the exponentially scaled modified Bessel function
@@ -324,15 +324,15 @@
* @param x the value to compute the bessel function of.
*/
static public double i0e(double x) throws ArithmeticException {
- double y;
-
- if( x < 0 ) x = -x;
- if( x <= 8.0 ) {
- y = (x/2.0) - 2.0;
- return( Arithmetic.chbevl( y, A_i0, 30 ) );
- }
+ double y;
+
+ if( x < 0 ) x = -x;
+ if( x <= 8.0 ) {
+ y = (x/2.0) - 2.0;
+ return( Arithmetic.chbevl( y, A_i0, 30 ) );
+ }
- return( Arithmetic.chbevl( 32.0/x - 2.0, B_i0, 25 ) / Math.sqrt(x) );
+ return( Arithmetic.chbevl( 32.0/x - 2.0, B_i0, 25 ) / Math.sqrt(x) );
}
/**
* Returns the modified Bessel function of order 1 of the
@@ -347,21 +347,21 @@
* @param x the value to compute the bessel function of.
*/
static public double i1(double x) throws ArithmeticException {
- double y, z;
+ double y, z;
- z = Math.abs(x);
- if( z <= 8.0 )
- {
- y = (z/2.0) - 2.0;
- z = Arithmetic.chbevl( y, A_i1, 29 ) * z * Math.exp(z);
- }
- else
- {
- z = Math.exp(z) * Arithmetic.chbevl( 32.0/z - 2.0, B_i1, 25 ) / Math.sqrt(z);
- }
- if( x < 0.0 )
- z = -z;
- return( z );
+ z = Math.abs(x);
+ if( z <= 8.0 )
+ {
+ y = (z/2.0) - 2.0;
+ z = Arithmetic.chbevl( y, A_i1, 29 ) * z * Math.exp(z);
+ }
+ else
+ {
+ z = Math.exp(z) * Arithmetic.chbevl( 32.0/z - 2.0, B_i1, 25 ) / Math.sqrt(z);
+ }
+ if( x < 0.0 )
+ z = -z;
+ return( z );
}
/**
* Returns the exponentially scaled modified Bessel function
@@ -376,16 +376,16 @@
z = Math.abs(x);
if( z <= 8.0 )
- {
- y = (z/2.0) - 2.0;
- z = Arithmetic.chbevl( y, A_i1, 29 ) * z;
- }
+ {
+ y = (z/2.0) - 2.0;
+ z = Arithmetic.chbevl( y, A_i1, 29 ) * z;
+ }
else
- {
- z = Arithmetic.chbevl( 32.0/z - 2.0, B_i1, 25 ) / Math.sqrt(z);
- }
+ {
+ z = Arithmetic.chbevl( 32.0/z - 2.0, B_i1, 25 ) / Math.sqrt(z);
+ }
if( x < 0.0 )
- z = -z;
+ z = -z;
return( z );
}
/**
@@ -393,64 +393,64 @@
* @param x the value to compute the bessel function of.
*/
static public double j0(double x) throws ArithmeticException {
- double ax;
+ double ax;
- if( (ax=Math.abs(x)) < 8.0 ) {
- double y=x*x;
- double ans1=57568490574.0+y*(-13362590354.0+y*(651619640.7
- +y*(-11214424.18+y*(77392.33017+y*(-184.9052456)))));
- double ans2=57568490411.0+y*(1029532985.0+y*(9494680.718
- +y*(59272.64853+y*(267.8532712+y*1.0))));
-
- return ans1/ans2;
-
- }
- else {
- double z=8.0/ax;
- double y=z*z;
- double xx=ax-0.785398164;
- double ans1=1.0+y*(-0.1098628627e-2+y*(0.2734510407e-4
- +y*(-0.2073370639e-5+y*0.2093887211e-6)));
- double ans2 = -0.1562499995e-1+y*(0.1430488765e-3
- +y*(-0.6911147651e-5+y*(0.7621095161e-6
- -y*0.934935152e-7)));
-
- return Math.sqrt(0.636619772/ax)*
- (Math.cos(xx)*ans1-z*Math.sin(xx)*ans2);
- }
+ if( (ax=Math.abs(x)) < 8.0 ) {
+ double y=x*x;
+ double ans1=57568490574.0+y*(-13362590354.0+y*(651619640.7
+ +y*(-11214424.18+y*(77392.33017+y*(-184.9052456)))));
+ double ans2=57568490411.0+y*(1029532985.0+y*(9494680.718
+ +y*(59272.64853+y*(267.8532712+y*1.0))));
+
+ return ans1/ans2;
+
+ }
+ else {
+ double z=8.0/ax;
+ double y=z*z;
+ double xx=ax-0.785398164;
+ double ans1=1.0+y*(-0.1098628627e-2+y*(0.2734510407e-4
+ +y*(-0.2073370639e-5+y*0.2093887211e-6)));
+ double ans2 = -0.1562499995e-1+y*(0.1430488765e-3
+ +y*(-0.6911147651e-5+y*(0.7621095161e-6
+ -y*0.934935152e-7)));
+
+ return Math.sqrt(0.636619772/ax)*
+ (Math.cos(xx)*ans1-z*Math.sin(xx)*ans2);
+ }
}
/**
* Returns the Bessel function of the first kind of order 1 of the argument.
* @param x the value to compute the bessel function of.
*/
static public double j1(double x) throws ArithmeticException {
- double ax;
- double y;
- double ans1, ans2;
-
- if ((ax=Math.abs(x)) < 8.0) {
- y=x*x;
- ans1=x*(72362614232.0+y*(-7895059235.0+y*(242396853.1
- +y*(-2972611.439+y*(15704.48260+y*(-30.16036606))))));
- ans2=144725228442.0+y*(2300535178.0+y*(18583304.74
- +y*(99447.43394+y*(376.9991397+y*1.0))));
- return ans1/ans2;
- }
- else {
- double z=8.0/ax;
- double xx=ax-2.356194491;
- y=z*z;
-
- ans1=1.0+y*(0.183105e-2+y*(-0.3516396496e-4
- +y*(0.2457520174e-5+y*(-0.240337019e-6))));
- ans2=0.04687499995+y*(-0.2002690873e-3
- +y*(0.8449199096e-5+y*(-0.88228987e-6
- +y*0.105787412e-6)));
- double ans=Math.sqrt(0.636619772/ax)*
- (Math.cos(xx)*ans1-z*Math.sin(xx)*ans2);
- if (x < 0.0) ans = -ans;
- return ans;
- }
+ double ax;
+ double y;
+ double ans1, ans2;
+
+ if ((ax=Math.abs(x)) < 8.0) {
+ y=x*x;
+ ans1=x*(72362614232.0+y*(-7895059235.0+y*(242396853.1
+ +y*(-2972611.439+y*(15704.48260+y*(-30.16036606))))));
+ ans2=144725228442.0+y*(2300535178.0+y*(18583304.74
+ +y*(99447.43394+y*(376.9991397+y*1.0))));
+ return ans1/ans2;
+ }
+ else {
+ double z=8.0/ax;
+ double xx=ax-2.356194491;
+ y=z*z;
+
+ ans1=1.0+y*(0.183105e-2+y*(-0.3516396496e-4
+ +y*(0.2457520174e-5+y*(-0.240337019e-6))));
+ ans2=0.04687499995+y*(-0.2002690873e-3
+ +y*(0.8449199096e-5+y*(-0.88228987e-6
+ +y*0.105787412e-6)));
+ double ans=Math.sqrt(0.636619772/ax)*
+ (Math.cos(xx)*ans1-z*Math.sin(xx)*ans2);
+ if (x < 0.0) ans = -ans;
+ return ans;
+ }
}
/**
* Returns the Bessel function of the first kind of order <tt>n</tt> of the argument.
@@ -458,55 +458,55 @@
* @param x the value to compute the bessel function of.
*/
static public double jn(int n, double x) throws ArithmeticException {
- int j,m;
- double ax,bj,bjm,bjp,sum,tox,ans;
- boolean jsum;
-
- final double ACC = 40.0;
- final double BIGNO = 1.0e+10;
- final double BIGNI = 1.0e-10;
-
- if(n == 0) return j0(x);
- if(n == 1) return j1(x);
-
- ax=Math.abs(x);
- if(ax == 0.0) return 0.0;
-
- if (ax > (double)n) {
- tox=2.0/ax;
- bjm=j0(ax);
- bj=j1(ax);
- for (j=1;j<n;j++) {
- bjp=j*tox*bj-bjm;
- bjm=bj;
- bj=bjp;
- }
- ans=bj;
- }
- else {
- tox=2.0/ax;
- m=2*((n+(int)Math.sqrt(ACC*n))/2);
- jsum=false;
- bjp=ans=sum=0.0;
- bj=1.0;
- for (j=m;j>0;j--) {
- bjm=j*tox*bj-bjp;
- bjp=bj;
- bj=bjm;
- if (Math.abs(bj) > BIGNO) {
- bj *= BIGNI;
- bjp *= BIGNI;
- ans *= BIGNI;
- sum *= BIGNI;
- }
- if (jsum) sum += bj;
- jsum=!jsum;
- if (j == n) ans=bjp;
- }
- sum=2.0*sum-bj;
- ans /= sum;
- }
- return x < 0.0 && n%2 == 1 ? -ans : ans;
+ int j,m;
+ double ax,bj,bjm,bjp,sum,tox,ans;
+ boolean jsum;
+
+ final double ACC = 40.0;
+ final double BIGNO = 1.0e+10;
+ final double BIGNI = 1.0e-10;
+
+ if(n == 0) return j0(x);
+ if(n == 1) return j1(x);
+
+ ax=Math.abs(x);
+ if(ax == 0.0) return 0.0;
+
+ if (ax > (double)n) {
+ tox=2.0/ax;
+ bjm=j0(ax);
+ bj=j1(ax);
+ for (j=1;j<n;j++) {
+ bjp=j*tox*bj-bjm;
+ bjm=bj;
+ bj=bjp;
+ }
+ ans=bj;
+ }
+ else {
+ tox=2.0/ax;
+ m=2*((n+(int)Math.sqrt(ACC*n))/2);
+ jsum=false;
+ bjp=ans=sum=0.0;
+ bj=1.0;
+ for (j=m;j>0;j--) {
+ bjm=j*tox*bj-bjp;
+ bjp=bj;
+ bj=bjm;
+ if (Math.abs(bj) > BIGNO) {
+ bj *= BIGNI;
+ bjp *= BIGNI;
+ ans *= BIGNI;
+ sum *= BIGNI;
+ }
+ if (jsum) sum += bj;
+ jsum=!jsum;
+ if (j == n) ans=bjp;
+ }
+ sum=2.0*sum-bj;
+ ans /= sum;
+ }
+ return x < 0.0 && n%2 == 1 ? -ans : ans;
}
/**
* Returns the modified Bessel function of the third kind
@@ -519,18 +519,18 @@
* @param x the value to compute the bessel function of.
*/
static public double k0(double x) throws ArithmeticException {
- double y, z;
+ double y, z;
- if( x <= 0.0 ) throw new ArithmeticException();
- if( x <= 2.0 ) {
- y = x * x - 2.0;
- y = Arithmetic.chbevl( y, A_k0, 10 ) - Math.log( 0.5 * x ) * i0(x);
- return( y );
- }
-
- z = 8.0/x - 2.0;
- y = Math.exp(-x) * Arithmetic.chbevl( z, B_k0, 25 ) / Math.sqrt(x);
- return(y);
+ if( x <= 0.0 ) throw new ArithmeticException();
+ if( x <= 2.0 ) {
+ y = x * x - 2.0;
+ y = Arithmetic.chbevl( y, A_k0, 10 ) - Math.log( 0.5 * x ) * i0(x);
+ return( y );
+ }
+
+ z = 8.0/x - 2.0;
+ y = Math.exp(-x) * Arithmetic.chbevl( z, B_k0, 25 ) / Math.sqrt(x);
+ return(y);
}
/**
* Returns the exponentially scaled modified Bessel function
@@ -539,17 +539,17 @@
* @param x the value to compute the bessel function of.
*/
static public double k0e(double x) throws ArithmeticException {
- double y;
+ double y;
- if( x <= 0.0 ) throw new ArithmeticException();
- if( x <= 2.0 ) {
- y = x * x - 2.0;
- y = Arithmetic.chbevl( y, A_k0, 10 ) - Math.log( 0.5 * x ) * i0(x);
- return( y * Math.exp(x) );
- }
+ if( x <= 0.0 ) throw new ArithmeticException();
+ if( x <= 2.0 ) {
+ y = x * x - 2.0;
+ y = Arithmetic.chbevl( y, A_k0, 10 ) - Math.log( 0.5 * x ) * i0(x);
+ return( y * Math.exp(x) );
+ }
- y = Arithmetic.chbevl( 8.0/x - 2.0, B_k0, 25 ) / Math.sqrt(x);
- return(y);
+ y = Arithmetic.chbevl( 8.0/x - 2.0, B_k0, 25 ) / Math.sqrt(x);
+ return(y);
}
/**
* Returns the modified Bessel function of the third kind
@@ -562,17 +562,17 @@
* @param x the value to compute the bessel function of.
*/
static public double k1(double x) throws ArithmeticException {
- double y, z;
+ double y, z;
- z = 0.5 * x;
- if( z <= 0.0 ) throw new ArithmeticException();
- if( x <= 2.0 ) {
- y = x * x - 2.0;
- y = Math.log(z) * i1(x) + Arithmetic.chbevl( y, A_k1, 11 ) / x;
- return( y );
- }
+ z = 0.5 * x;
+ if( z <= 0.0 ) throw new ArithmeticException();
+ if( x <= 2.0 ) {
+ y = x * x - 2.0;
+ y = Math.log(z) * i1(x) + Arithmetic.chbevl( y, A_k1, 11 ) / x;
+ return( y );
+ }
- return( Math.exp(-x) * Arithmetic.chbevl( 8.0/x - 2.0, B_k1, 25 ) / Math.sqrt(x) );
+ return( Math.exp(-x) * Arithmetic.chbevl( 8.0/x - 2.0, B_k1, 25 ) / Math.sqrt(x) );
}
/**
* Returns the exponentially scaled modified Bessel function
@@ -583,16 +583,16 @@
* @param x the value to compute the bessel function of.
*/
static public double k1e(double x) throws ArithmeticException {
- double y;
+ double y;
- if( x <= 0.0 ) throw new ArithmeticException();
- if( x <= 2.0 ) {
- y = x * x - 2.0;
- y = Math.log( 0.5 * x ) * i1(x) + Arithmetic.chbevl( y, A_k1, 11 ) / x;
- return( y * Math.exp(x) );
- }
+ if( x <= 0.0 ) throw new ArithmeticException();
+ if( x <= 2.0 ) {
+ y = x * x - 2.0;
+ y = Math.log( 0.5 * x ) * i1(x) + Arithmetic.chbevl( y, A_k1, 11 ) / x;
+ return( y * Math.exp(x) );
+ }
- return( Arithmetic.chbevl( 8.0/x - 2.0, B_k1, 25 ) / Math.sqrt(x) );
+ return( Arithmetic.chbevl( 8.0/x - 2.0, B_k1, 25 ) / Math.sqrt(x) );
}
/**
* Returns the modified Bessel function of the third kind
@@ -608,220 +608,220 @@
static public double kn(int nn, double x) throws ArithmeticException {
/*
Algorithm for Kn.
- n-1
- -n - (n-k-1)! 2 k
+ n-1
+ -n - (n-k-1)! 2 k
K (x) = 0.5 (x/2) > -------- (-x /4)
n - k!
- k=0
+ k=0
- inf. 2 k
- n n - (x /4)
+ inf. 2 k
+ n n - (x /4)
+ (-1) 0.5(x/2) > {p(k+1) + p(n+k+1) - 2log(x/2)} ---------
- - k! (n+k)!
- k=0
+ - k! (n+k)!
+ k=0
where p(m) is the psi function: p(1) = -EUL and
- m-1
- -
- p(m) = -EUL + > 1/k
- -
- k=1
+ m-1
+ -
+ p(m) = -EUL + > 1/k
+ -
+ k=1
For large x,
- 2 2 2
- u-1 (u-1 )(u-3 )
+ 2 2 2
+ u-1 (u-1 )(u-3 )
K (z) = sqrt(pi/2z) exp(-z) { 1 + ------- + ------------ + ...}
v 1 2
- 1! (8z) 2! (8z)
+ 1! (8z) 2! (8z)
asymptotically, where
- 2
- u = 4 v .
+ 2
+ u = 4 v .
*/
- final double EUL = 5.772156649015328606065e-1;
- final double MAXNUM = Double.MAX_VALUE;
- final int MAXFAC = 31;
-
-
- double k, kf, nk1f, nkf, zn, t, s, z0, z;
- double ans, fn, pn, pk, zmn, tlg, tox;
- int i, n;
-
- if( nn < 0 )
- n = -nn;
- else
- n = nn;
-
- if( n > MAXFAC ) throw new ArithmeticException("Overflow");
- if( x <= 0.0 ) throw new IllegalArgumentException();
-
- if( x <= 9.55 ) {
- ans = 0.0;
- z0 = 0.25 * x * x;
- fn = 1.0;
- pn = 0.0;
- zmn = 1.0;
- tox = 2.0/x;
-
- if( n > 0 ) {
- /* compute factorial of n and psi(n) */
- pn = -EUL;
- k = 1.0;
- for( i=1; i<n; i++ ) {
- pn += 1.0/k;
- k += 1.0;
- fn *= k;
- }
-
- zmn = tox;
-
- if( n == 1 ) {
- ans = 1.0/x;
- }
- else {
- nk1f = fn/n;
- kf = 1.0;
- s = nk1f;
- z = -z0;
- zn = 1.0;
- for( i=1; i<n; i++ ) {
- nk1f = nk1f/(n-i);
- kf = kf * i;
- zn *= z;
- t = nk1f * zn / kf;
- s += t;
- if( (MAXNUM - Math.abs(t)) < Math.abs(s) ) throw new ArithmeticException("Overflow");
- if( (tox > 1.0) && ((MAXNUM/tox) < zmn) ) throw new ArithmeticException("Overflow");
- zmn *= tox;
- }
- s *= 0.5;
- t = Math.abs(s);
- if( (zmn > 1.0) && ((MAXNUM/zmn) < t) ) throw new ArithmeticException("Overflow");
- if( (t > 1.0) && ((MAXNUM/t) < zmn) ) throw new ArithmeticException("Overflow");
- ans = s * zmn;
- }
- }
-
-
- tlg = 2.0 * Math.log( 0.5 * x );
- pk = -EUL;
- if( n == 0 ) {
- pn = pk;
- t = 1.0;
- }
- else {
- pn = pn + 1.0/n;
- t = 1.0/fn;
- }
- s = (pk+pn-tlg)*t;
- k = 1.0;
- do {
- t *= z0 / (k * (k+n));
- pk += 1.0/k;
- pn += 1.0/(k+n);
- s += (pk+pn-tlg)*t;
- k += 1.0;
- }
- while( Math.abs(t/s) > MACHEP );
-
- s = 0.5 * s / zmn;
- if( (n & 1) > 0)
- s = -s;
- ans += s;
-
- return(ans);
- }
-
-
-
- /* Asymptotic expansion for Kn(x) */
- /* Converges to 1.4e-17 for x > 18.4 */
- if( x > MAXLOG ) throw new ArithmeticException("Underflow");
- k = n;
- pn = 4.0 * k * k;
- pk = 1.0;
- z0 = 8.0 * x;
- fn = 1.0;
- t = 1.0;
- s = t;
- nkf = MAXNUM;
- i = 0;
- do {
- z = pn - pk * pk;
- t = t * z /(fn * z0);
- nk1f = Math.abs(t);
- if( (i >= n) && (nk1f > nkf) ) {
- ans = Math.exp(-x) * Math.sqrt( Math.PI/(2.0*x) ) * s;
- return(ans);
- }
- nkf = nk1f;
- s += t;
- fn += 1.0;
- pk += 2.0;
- i += 1;
- } while( Math.abs(t/s) > MACHEP );
+ final double EUL = 5.772156649015328606065e-1;
+ final double MAXNUM = Double.MAX_VALUE;
+ final int MAXFAC = 31;
+
+
+ double k, kf, nk1f, nkf, zn, t, s, z0, z;
+ double ans, fn, pn, pk, zmn, tlg, tox;
+ int i, n;
+
+ if( nn < 0 )
+ n = -nn;
+ else
+ n = nn;
+
+ if( n > MAXFAC ) throw new ArithmeticException("Overflow");
+ if( x <= 0.0 ) throw new IllegalArgumentException();
+
+ if( x <= 9.55 ) {
+ ans = 0.0;
+ z0 = 0.25 * x * x;
+ fn = 1.0;
+ pn = 0.0;
+ zmn = 1.0;
+ tox = 2.0/x;
+
+ if( n > 0 ) {
+ /* compute factorial of n and psi(n) */
+ pn = -EUL;
+ k = 1.0;
+ for( i=1; i<n; i++ ) {
+ pn += 1.0/k;
+ k += 1.0;
+ fn *= k;
+ }
+
+ zmn = tox;
+
+ if( n == 1 ) {
+ ans = 1.0/x;
+ }
+ else {
+ nk1f = fn/n;
+ kf = 1.0;
+ s = nk1f;
+ z = -z0;
+ zn = 1.0;
+ for( i=1; i<n; i++ ) {
+ nk1f = nk1f/(n-i);
+ kf = kf * i;
+ zn *= z;
+ t = nk1f * zn / kf;
+ s += t;
+ if( (MAXNUM - Math.abs(t)) < Math.abs(s) ) throw new ArithmeticException("Overflow");
+ if( (tox > 1.0) && ((MAXNUM/tox) < zmn) ) throw new ArithmeticException("Overflow");
+ zmn *= tox;
+ }
+ s *= 0.5;
+ t = Math.abs(s);
+ if( (zmn > 1.0) && ((MAXNUM/zmn) < t) ) throw new ArithmeticException("Overflow");
+ if( (t > 1.0) && ((MAXNUM/t) < zmn) ) throw new ArithmeticException("Overflow");
+ ans = s * zmn;
+ }
+ }
+
+
+ tlg = 2.0 * Math.log( 0.5 * x );
+ pk = -EUL;
+ if( n == 0 ) {
+ pn = pk;
+ t = 1.0;
+ }
+ else {
+ pn = pn + 1.0/n;
+ t = 1.0/fn;
+ }
+ s = (pk+pn-tlg)*t;
+ k = 1.0;
+ do {
+ t *= z0 / (k * (k+n));
+ pk += 1.0/k;
+ pn += 1.0/(k+n);
+ s += (pk+pn-tlg)*t;
+ k += 1.0;
+ }
+ while( Math.abs(t/s) > MACHEP );
+
+ s = 0.5 * s / zmn;
+ if( (n & 1) > 0)
+ s = -s;
+ ans += s;
+
+ return(ans);
+ }
+
+
+
+ /* Asymptotic expansion for Kn(x) */
+ /* Converges to 1.4e-17 for x > 18.4 */
+ if( x > MAXLOG ) throw new ArithmeticException("Underflow");
+ k = n;
+ pn = 4.0 * k * k;
+ pk = 1.0;
+ z0 = 8.0 * x;
+ fn = 1.0;
+ t = 1.0;
+ s = t;
+ nkf = MAXNUM;
+ i = 0;
+ do {
+ z = pn - pk * pk;
+ t = t * z /(fn * z0);
+ nk1f = Math.abs(t);
+ if( (i >= n) && (nk1f > nkf) ) {
+ ans = Math.exp(-x) * Math.sqrt( Math.PI/(2.0*x) ) * s;
+ return(ans);
+ }
+ nkf = nk1f;
+ s += t;
+ fn += 1.0;
+ pk += 2.0;
+ i += 1;
+ } while( Math.abs(t/s) > MACHEP );
- ans = Math.exp(-x) * Math.sqrt( Math.PI/(2.0*x) ) * s;
- return(ans);
+ ans = Math.exp(-x) * Math.sqrt( Math.PI/(2.0*x) ) * s;
+ return(ans);
}
/**
* Returns the Bessel function of the second kind of order 0 of the argument.
* @param x the value to compute the bessel function of.
*/
static public double y0(double x) throws ArithmeticException {
- if (x < 8.0) {
- double y=x*x;
- double ans1 = -2957821389.0+y*(7062834065.0+y*(-512359803.6
- +y*(10879881.29+y*(-86327.92757+y*228.4622733))));
- double ans2=40076544269.0+y*(745249964.8+y*(7189466.438
- +y*(47447.26470+y*(226.1030244+y*1.0))));
-
- return (ans1/ans2)+0.636619772*j0(x)*Math.log(x);
- }
- else {
- double z=8.0/x;
- double y=z*z;
- double xx=x-0.785398164;
-
- double ans1=1.0+y*(-0.1098628627e-2+y*(0.2734510407e-4
- +y*(-0.2073370639e-5+y*0.2093887211e-6)));
- double ans2 = -0.1562499995e-1+y*(0.1430488765e-3
- +y*(-0.6911147651e-5+y*(0.7621095161e-6
- +y*(-0.934945152e-7))));
- return Math.sqrt(0.636619772/x)*
- (Math.sin(xx)*ans1+z*Math.cos(xx)*ans2);
- }
+ if (x < 8.0) {
+ double y=x*x;
+ double ans1 = -2957821389.0+y*(7062834065.0+y*(-512359803.6
+ +y*(10879881.29+y*(-86327.92757+y*228.4622733))));
+ double ans2=40076544269.0+y*(745249964.8+y*(7189466.438
+ +y*(47447.26470+y*(226.1030244+y*1.0))));
+
+ return (ans1/ans2)+0.636619772*j0(x)*Math.log(x);
+ }
+ else {
+ double z=8.0/x;
+ double y=z*z;
+ double xx=x-0.785398164;
+
+ double ans1=1.0+y*(-0.1098628627e-2+y*(0.2734510407e-4
+ +y*(-0.2073370639e-5+y*0.2093887211e-6)));
+ double ans2 = -0.1562499995e-1+y*(0.1430488765e-3
+ +y*(-0.6911147651e-5+y*(0.7621095161e-6
+ +y*(-0.934945152e-7))));
+ return Math.sqrt(0.636619772/x)*
+ (Math.sin(xx)*ans1+z*Math.cos(xx)*ans2);
+ }
}
/**
* Returns the Bessel function of the second kind of order 1 of the argument.
* @param x the value to compute the bessel function of.
*/
static public double y1(double x) throws ArithmeticException {
- if (x < 8.0) {
- double y=x*x;
- double ans1=x*(-0.4900604943e13+y*(0.1275274390e13
- +y*(-0.5153438139e11+y*(0.7349264551e9
- +y*(-0.4237922726e7+y*0.8511937935e4)))));
- double ans2=0.2499580570e14+y*(0.4244419664e12
- +y*(0.3733650367e10+y*(0.2245904002e8
- +y*(0.1020426050e6+y*(0.3549632885e3+y)))));
- return (ans1/ans2)+0.636619772*(j1(x)*Math.log(x)-1.0/x);
- }
- else {
- double z=8.0/x;
- double y=z*z;
- double xx=x-2.356194491;
- double ans1=1.0+y*(0.183105e-2+y*(-0.3516396496e-4
- +y*(0.2457520174e-5+y*(-0.240337019e-6))));
- double ans2=0.04687499995+y*(-0.2002690873e-3
- +y*(0.8449199096e-5+y*(-0.88228987e-6
- +y*0.105787412e-6)));
- return Math.sqrt(0.636619772/x)*
- (Math.sin(xx)*ans1+z*Math.cos(xx)*ans2);
- }
+ if (x < 8.0) {
+ double y=x*x;
+ double ans1=x*(-0.4900604943e13+y*(0.1275274390e13
+ +y*(-0.5153438139e11+y*(0.7349264551e9
+ +y*(-0.4237922726e7+y*0.8511937935e4)))));
+ double ans2=0.2499580570e14+y*(0.4244419664e12
+ +y*(0.3733650367e10+y*(0.2245904002e8
+ +y*(0.1020426050e6+y*(0.3549632885e3+y)))));
+ return (ans1/ans2)+0.636619772*(j1(x)*Math.log(x)-1.0/x);
+ }
+ else {
+ double z=8.0/x;
+ double y=z*z;
+ double xx=x-2.356194491;
+ double ans1=1.0+y*(0.183105e-2+y*(-0.3516396496e-4
+ +y*(0.2457520174e-5+y*(-0.240337019e-6))));
+ double ans2=0.04687499995+y*(-0.2002690873e-3
+ +y*(0.8449199096e-5+y*(-0.88228987e-6
+ +y*0.105787412e-6)));
+ return Math.sqrt(0.636619772/x)*
+ (Math.sin(xx)*ans1+z*Math.cos(xx)*ans2);
+ }
}
/**
* Returns the Bessel function of the second kind of order <tt>n</tt> of the argument.
@@ -829,19 +829,19 @@
* @param x the value to compute the bessel function of.
*/
static public double yn(int n, double x) throws ArithmeticException {
- double by,bym,byp,tox;
+ double by,bym,byp,tox;
- if(n == 0) return y0(x);
- if(n == 1) return y1(x);
+ if(n == 0) return y0(x);
+ if(n == 1) return y1(x);
- tox=2.0/x;
- by=y1(x);
- bym=y0(x);
- for (int j=1;j<n;j++) {
- byp=j*tox*by-bym;
- bym=by;
- by=byp;
- }
- return by;
+ tox=2.0/x;
+ by=y1(x);
+ bym=y0(x);
+ for (int j=1;j<n;j++) {
+ byp=j*tox*by-bym;
+ bym=by;
+ by=byp;
+ }
+ return by;
}
}
Modified: lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/math/Constants.java
URL: http://svn.apache.org/viewvc/lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/math/Constants.java?rev=883972&r1=883971&r2=883972&view=diff
==============================================================================
--- lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/math/Constants.java (original)
+++ lucene/mahout/trunk/matrix/src/main/java/org/apache/mahout/jet/math/Constants.java Wed Nov 25 03:31:47 2009
@@ -19,19 +19,19 @@
/*
* machine constants
*/
- protected static final double MACHEP = 1.11022302462515654042E-16;
- protected static final double MAXLOG = 7.09782712893383996732E2;
- protected static final double MINLOG = -7.451332191019412076235E2;
- protected static final double MAXGAM = 171.624376956302725;
- protected static final double SQTPI = 2.50662827463100050242E0;
- protected static final double SQRTH = 7.07106781186547524401E-1;
- protected static final double LOGPI = 1.14472988584940017414;
+ protected static final double MACHEP = 1.11022302462515654042E-16;
+ protected static final double MAXLOG = 7.09782712893383996732E2;
+ protected static final double MINLOG = -7.451332191019412076235E2;
+ protected static final double MAXGAM = 171.624376956302725;
+ protected static final double SQTPI = 2.50662827463100050242E0;
+ protected static final double SQRTH = 7.07106781186547524401E-1;
+ protected static final double LOGPI = 1.14472988584940017414;
- protected static final double big = 4.503599627370496e15;
- protected static final double biginv = 2.22044604925031308085e-16;
+ protected static final double big = 4.503599627370496e15;
+ protected static final double biginv = 2.22044604925031308085e-16;
- /*
+ /*
* MACHEP = 1.38777878078144567553E-17 2**-56
* MAXLOG = 8.8029691931113054295988E1 log(2**127)
* MINLOG = -8.872283911167299960540E1 log(2**-128)