Aleksey Khudyakov avatar Aleksey Khudyakov committed e6d2b88

Add strilingError. it was omitted somehow

Comments (0)

Files changed (1)

Numeric/SpecFunctions.hs

     -- * Factorial
   , factorial
   , logFactorial
+  , stirlingError
     -- * Combinatorics
   , choose
     -- * References
           z = ((-(5.95238095238e-4 * y) + 7.936500793651e-4) * y -
                2.7777777777778e-3) * y + 8.3333333333333e-2
 
+-- | Calculate the error term of the Stirling approximation.  This is
+-- only defined for non-negative values.
+--
+-- > stirlingError @n@ = @log(n!) - log(sqrt(2*pi*n)*(n/e)^n)
+stirlingError :: Double -> Double
+stirlingError n
+  | n <= 15.0   = case properFraction (n+n) of
+                    (i,0) -> sfe `U.unsafeIndex` i
+                    _     -> logGamma (n+1.0) - (n+0.5) * log n + n -
+                             m_ln_sqrt_2_pi
+  | n > 500     = (s0-s1/nn)/n
+  | n > 80      = (s0-(s1-s2/nn)/nn)/n
+  | n > 35      = (s0-(s1-(s2-s3/nn)/nn)/nn)/n
+  | otherwise   = (s0-(s1-(s2-(s3-s4/nn)/nn)/nn)/nn)/n
+  where
+    nn = n*n
+    s0 = 0.083333333333333333333        -- 1/12
+    s1 = 0.00277777777777777777778      -- 1/360
+    s2 = 0.00079365079365079365079365   -- 1/1260
+    s3 = 0.000595238095238095238095238  -- 1/1680
+    s4 = 0.0008417508417508417508417508 -- 1/1188
+    sfe = U.fromList [ 0.0,
+                0.1534264097200273452913848,   0.0810614667953272582196702,
+                0.0548141210519176538961390,   0.0413406959554092940938221,
+                0.03316287351993628748511048,  0.02767792568499833914878929,
+                0.02374616365629749597132920,  0.02079067210376509311152277,
+                0.01848845053267318523077934,  0.01664469118982119216319487,
+                0.01513497322191737887351255,  0.01387612882307074799874573,
+                0.01281046524292022692424986,  0.01189670994589177009505572,
+                0.01110455975820691732662991,  0.010411265261972096497478567,
+                0.009799416126158803298389475, 0.009255462182712732917728637,
+                0.008768700134139385462952823, 0.008330563433362871256469318,
+                0.007934114564314020547248100, 0.007573675487951840794972024,
+                0.007244554301320383179543912, 0.006942840107209529865664152,
+                0.006665247032707682442354394, 0.006408994188004207068439631,
+                0.006171712263039457647532867, 0.005951370112758847735624416,
+                0.005746216513010115682023589, 0.005554733551962801371038690 ]
 
 
 ----------------------------------------------------------------
Tip: Filter by directory path e.g. /media app.js to search for public/media/app.js.
Tip: Use camelCasing e.g. ProjME to search for ProjectModifiedEvent.java.
Tip: Filter by extension type e.g. /repo .js to search for all .js files in the /repo directory.
Tip: Separate your search with spaces e.g. /ssh pom.xml to search for src/ssh/pom.xml.
Tip: Use ↑ and ↓ arrow keys to navigate and return to view the file.
Tip: You can also navigate files with Ctrl+j (next) and Ctrl+k (previous) and view the file with Ctrl+o.
Tip: You can also navigate files with Alt+j (next) and Alt+k (previous) and view the file with Alt+o.