diff --git a/src/FSharp.Stats/SpecialFunctions/Factorial.fs b/src/FSharp.Stats/SpecialFunctions/Factorial.fs index fffe8ff1..139220e6 100644 --- a/src/FSharp.Stats/SpecialFunctions/Factorial.fs +++ b/src/FSharp.Stats/SpecialFunctions/Factorial.fs @@ -8,13 +8,7 @@ module Factorial = // This is the largest integer value for which the factorial function doesn't overflow the floating point format. let private FactorialMax = 170 - //let inline facti n = - // let rec loop acc n = - // if n<=LanguagePrimitives.GenericOne - // then acc - // else loop (n*acc) (n-LanguagePrimitives.GenericOne) - // loop LanguagePrimitives.GenericOne n - + //cache of all factorials in [0..170] let private FactorialCache = //let cache = [| 0 .. FactorialMax |] |> Array.map (fun a -> float a) let cache = Array.zeroCreate (FactorialMax + 1) @@ -29,61 +23,22 @@ module Factorial = let private FactorialLnCache = let cache = Array.zeroCreate (FactorialLnNTop + 1) for i=0 to FactorialLnNTop do - cache.[i] <-Gamma.gammaLn ((float i) + 1.0) + cache.[i] <- Gamma.gammaLn ((float i) + 1.0) cache - /// The factorial functions takes an int x and computes x!. This function will not overflow - /// the floating point format as long as x is at most 170. + /// The factorial functions takes an int x and returns x!. This function will not overflow + /// the floating point format as long as x is at most 170, and will return +infinity for all values > 170 let factorial (x:int) = if x <= FactorialMax then FactorialCache.[x] else System.Double.PositiveInfinity - /// Computes the natural logarithm of the factorial function. + /// Computes the natural logarithm of the factorial function, let factorialLn (x: int) : float = if x < 0 then failwith "Log factorial not defined for n < 0" //if x <= 1 then 0.0 else Gamma.gammaLn ((float x) + 1.0) if x <= FactorialLnNTop then FactorialLnCache.[x] else - Gamma.gammaLn ((float x) + 1.0) - - - -// let cacheFactorial size f = -// let cache = Array.zeroCreate size -// (fun x -> -// match cache.[x] with -// | 0 -> -// let v = f(x) -// cache.[x] <- v -// v -// | y -> y ) -// -// -// let rec factorial = -// cacheFactorial FactorialMax ( -// fun x -> -// if (x = 0) then -// 1 -// else -// x * factorial(x - 1)) -// -// let rec factorialLn = -// cacheFactorial FactorialMax ( -// fun x -> -// if (x = 0) then -// 1 -// else -// x * factorial(x - 1)) -// -// let inline facti n = -// let rec loop acc n = -// if n<=LanguagePrimitives.GenericOne -// then acc -// else loop (n*acc) (n-LanguagePrimitives.GenericOne) -// loop LanguagePrimitives.GenericOne n - - - \ No newline at end of file + Gamma.gammaLn ((float x) + 1.0) \ No newline at end of file diff --git a/tests/FSharp.Stats.Tests/Main.fs b/tests/FSharp.Stats.Tests/Main.fs index 8861812f..e5acde62 100644 --- a/tests/FSharp.Stats.Tests/Main.fs +++ b/tests/FSharp.Stats.Tests/Main.fs @@ -16,7 +16,7 @@ let main argv = //=========================== Special Functions ========================================================= Tests.runTestsWithCLIArgs [] argv SpecialFunctionsTests.gammaFunctionsTests |> ignore Tests.runTestsWithCLIArgs [] argv SpecialFunctionsTests.betaFunctionsTests |> ignore - + Tests.runTestsWithCLIArgs [] argv SpecialFunctionsTests.factorialTests |> ignore //================================ Algebra ============================================================== Tests.runTestsWithCLIArgs [] argv LinearAlgebraTests.managedSVDTests |> ignore Tests.runTestsWithCLIArgs [] argv LinearAlgebraTests.nullspace |> ignore diff --git a/tests/FSharp.Stats.Tests/SpecialFunctions.fs b/tests/FSharp.Stats.Tests/SpecialFunctions.fs index 1456263d..cadedf82 100644 --- a/tests/FSharp.Stats.Tests/SpecialFunctions.fs +++ b/tests/FSharp.Stats.Tests/SpecialFunctions.fs @@ -68,4 +68,38 @@ let betaFunctionsTests = let result = Beta.betaLn 0.0342 170. Expect.floatClose Accuracy.veryHigh result 3.1811881124242447 "Should be equal (double precision)" //rtol=1e-14, atol=0 + ] + +[] +let factorialTests = + testList "SpecialFunctions.Factorial" [ + testCase "Double overflow" (fun _ -> + Expect.equal (Factorial.factorial 171) infinity "Expected factorial of a number larger than 170 to result in infinity (171! is larger than max double)" + ) + testCase "0! equals 1" (fun _ -> + Expect.equal (Factorial.factorial 0) 1. "Expected factorial of 0 to be 1." + ) + testCase "69!" (fun _ -> + Expect.floatClose Accuracy.low (Factorial.factorial 69) 1.7112245e+98 "Expected factorial of 69 to be 1.7112245e+98" + ) + testCase "factorial not defined for negative numbers" (fun _ -> + Expect.throws (fun _ -> Factorial.factorial -69421337 |> ignore) "Expected factorial to fail for negative values" + ) + ] + +[] +let FactorialLnTests = + testList "SpecialFunctions.lnFactorial" [ + testCase "Large value" (fun _ -> + Expect.floatClose Accuracy.low (Factorial.factorialLn 6942) 54467.727976695301612523565124699078303834231913072759124392135342 "Expected factorial of a number larger than 170 to result in infinity (171! is larger than max double)" + ) + testCase "ln(0!) equals 0" (fun _ -> + Expect.equal (Factorial.factorialLn 0) 0. "Expected factorial of 0 to be 1." + ) + testCase "ln(69!)" (fun _ -> + Expect.floatClose Accuracy.low 226.19054832372759333227016852232261788323276357495863628461257077 (Factorial.factorialLn 69) "Expected factorialLn of 69 to be 226.19054832372759333227016852232261788323276357495863628461257077" + ) + testCase "factorialLn not defined for negative numbers" (fun _ -> + Expect.throws (fun _ -> Factorial.factorialLn -69421337 |> ignore) "Expected factorialLn to fail for negative values" + ) ] \ No newline at end of file