Add factorial and factorialLn tests, codebase cleanup

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
-
-
-
+            Gamma.gammaLn ((float x) + 1.0)

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

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
 
+    ]
+
+[<Tests>]
+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"
+        )
+    ]
+
+[<Tests>]
+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"
+        )
     ]