2 people like it.
Like the snippet!

Inverse Gamma function and Inverse Factorial by Approximation

Inverse Gamma function and Inverse Factorial by Approximation based on these references: http://mathforum.org/kb/message.jspa?messageID=342551&tstart=0 https://github.com/DarkoVeberic/LambertW

 1: 
 2: 
 3: 
 4: 
 5: 
 6: 
 7: 
 8: 
 9: 
10: 
11: 
12: 
13: 
14: 
15: 
16: 
17: 
18: 
19: 
20: 
21: 
22: 
23: 
24: 
25: 
26: 
27: 
28: 
29: 
30: 
31: 
32: 
33: 
34: 
35: 
36: 
37: 
38: 
39: 
40: 
41: 
42: 
43: 
44: 
45: 
46: 
47: 
48: 
49: 
50: 
51: 
52: 
53: 
54: 
55: 
56: 
57: 
58:

//More info and native code at:
//https://github.com/fabiogaluppo/samples/tree/master/fragments/InverseGamma

//References:
//David W. Cantrell's Inverse gamma function (and Inverse factorial): http://mathforum.org/kb/message.jspa?messageID=342551&tstart=0
//DarkoVeberic's C++ implementation of the Lambert W(x) function: https://github.com/DarkoVeberic/LambertW

open System
open System.Runtime.InteropServices

[<DllImport("bin\\LambertW.dll", CallingConvention=CallingConvention.Cdecl)>]
extern float LambertW_0(float x)

let c = 0.036534
let ln = Math.Log
let pi = Math.PI
let L x = ln((x + c) / Math.Sqrt(2. * pi))
let W x = LambertW_0 x
let e = Math.E

let AIG x =
  //Approximated Inverse Gamma
  L(x) / (W (L(x) / e)) + 1. / 2.

let InvFact x =
  //Inverse Factorial in terms of rounded AIG
  Math.Round(AIG x) - 1.

//Tests:
let showAIG x =
  printfn "AIG(%25.1f) = %9.6f" x (AIG x)

let showInvFact x =  
  printfn "InvFact(%21.1f) = %9.6f" x (InvFact x)

printfn "Inverse Gamma function and Inverse Factorial by Approximation"
showAIG 1.
showAIG 24.
showAIG 362880.
showAIG 1.216451e+17
showInvFact 6.
showInvFact 24.
showInvFact 3628800.
showInvFact 2.432902e+18

(*

Inverse Gamma function and Inverse Factorial by Approximation
AIG(                      1.0) =  2.021203
AIG(                     24.0) =  4.994871
AIG(                 362880.0) =  9.998053
AIG(     121645100000000000.0) = 19.999281
InvFact(                  6.0) =  3.000000
InvFact(                 24.0) =  4.000000
InvFact(            3628800.0) = 10.000000
InvFact(2432902000000000000.0) = 20.000000

*)

namespace System

namespace System.Runtime

namespace System.Runtime.InteropServices

Multiple items
type DllImportAttribute =
  inherit Attribute
  new : dllName:string -> DllImportAttribute
  val EntryPoint : string
  val CharSet : CharSet
  val SetLastError : bool
  val ExactSpelling : bool
  val PreserveSig : bool
  val CallingConvention : CallingConvention
  val BestFitMapping : bool
  val ThrowOnUnmappableChar : bool
  member Value : string

Full name: System.Runtime.InteropServices.DllImportAttribute

--------------------
DllImportAttribute(dllName: string) : unit

field CallingConvention.Cdecl = 2

Multiple items
val float : value:'T -> float (requires member op_Explicit)

Full name: Microsoft.FSharp.Core.Operators.float

--------------------
type float = Double

Full name: Microsoft.FSharp.Core.float

--------------------
type float<'Measure> = float

Full name: Microsoft.FSharp.Core.float<_>

val LambertW_0 : x:float -> float

Full name: Script.LambertW_0

val x : float

val c : float

Full name: Script.c

val ln : arg00:float -> float

Full name: Script.ln

type Math =
  static val PI : float
  static val E : float
  static member Abs : value:sbyte -> sbyte + 6 overloads
  static member Acos : d:float -> float
  static member Asin : d:float -> float
  static member Atan : d:float -> float
  static member Atan2 : y:float * x:float -> float
  static member BigMul : a:int * b:int -> int64
  static member Ceiling : d:decimal -> decimal + 1 overload
  static member Cos : d:float -> float
  ...

Full name: System.Math

Math.Log(d: float) : float
Math.Log(a: float, newBase: float) : float

val pi : float

Full name: Script.pi

field Math.PI = 3.14159265359

val L : x:float -> float

Full name: Script.L

Math.Sqrt(d: float) : float

val W : x:float -> float

Full name: Script.W

val e : float

Full name: Script.e

field Math.E = 2.71828182846

val AIG : x:float -> float

Full name: Script.AIG

val InvFact : x:float -> float

Full name: Script.InvFact

Math.Round(d: decimal) : decimal
Math.Round(a: float) : float
Math.Round(d: decimal, mode: MidpointRounding) : decimal
Math.Round(d: decimal, decimals: int) : decimal
Math.Round(value: float, mode: MidpointRounding) : float
Math.Round(value: float, digits: int) : float
Math.Round(d: decimal, decimals: int, mode: MidpointRounding) : decimal
Math.Round(value: float, digits: int, mode: MidpointRounding) : float

val showAIG : x:float -> unit

Full name: Script.showAIG

val printfn : format:Printf.TextWriterFormat<'T> -> 'T

Full name: Microsoft.FSharp.Core.ExtraTopLevelOperators.printfn

val showInvFact : x:float -> unit

Full name: Script.showInvFact

Raw view Test code New version

More information

Link:	http://fssnip.net/rw
Posted:	9 years ago
Author:	Fabio Galuppo
Tags:	gamma function , factorial , interop