0
votes

in F#, how does one write a generic-math step function?

An (Oliver) Heaviside step function is function that returns zero if x is negative, otherwise it retuns one.

Here is a summary of my attempts so far:

// attempt 1:
let inline stepFct1< ^T when ^T : (static member op_GreaterThan: ^T * float -> bool) 
    >     (x:^T) : ^T = 
    //if (^T : (static member op_GreaterThan) (x 0.0) ) then x  //ouch fails also
    if  (>) x 0.0 then x
    else 0.0

compiler says: error FS0001: A type parameter is missing a constraint 'when ^T : comparison'

// attempt 2:
let inline stepFct2<^T when ^T : (static member (>): ^T * ^T -> bool) > (x:^T) : ^T = 
    match x with 
    | x when x > 0.0 -> 1.0
    | 0.0

FSC says: error FS0010: Unexpected infix operator in pattern

Motivation:

I am trying to rewrite Ian's Cumulative-Normal and Black-Scholes functions here to use Automatic Differentiation (DiffSharp). Ian's Cumulative Normal works on floats, I would like a generic version that works on any numeric type, including AutoDiff.DualG. The cumulative normal function contains a "greater than" statement.

EDIT: Gustavo, thanks, I have accepted your answer - the simple step function now compiles.

But it doesn't seem to help with the Cumulative Normal case. Given this code:

// Cumulative Normal Distribution Function - attempt to write a generic version
let inline CDF(x:^T) : ^T = 
    let (b1,b2,b3)  = (0.319381530, -0.356563782, 1.781477937)
    let (b4,b5)     = (-1.821255978, 1.330274429)
    let (p , c )    = (0.2316419  ,  0.39894228)
    let (zero, one) = (LanguagePrimitives.GenericZero, LanguagePrimitives.GenericOne)
    if x > zero then
        let t = one / (one + p * x) 
        (one - c * exp( -x * x / 2.0)* t * (t*(t*(t*(t*b5+b4)+b3)+b2)+b1)) 
    else
        let t = 1.0 / (one - p * x) 
        (c * exp( -x * x / 2.0)* t * (t*(t*(t*(t*b5+b4)+b3)+b2)+b1))

FSI says:

C:\stdin(116,32): warning FS0064: This construct causes code to be less generic 
than indicated by the type annotations. 
The type variable 'T has been constrained to be type 'float'.

val inline CDF : x:float -> float
> CDF 0.1M;;
CDF 0.1M;;
----^^^^
C:\stdin(122,5): error FS0001: This expression was expected to have type
    float
but here has type
    decimal
>

Does anyone know how to make CDF generic?

1
Your update is really a new question, but you have missed the point of the answer. If you have a floating point literal anywhere, the code is no longer generic. You will need to find generic ways to construct the constants. - John Palmer
@JohnPalmer is right, that's more challenging. The lib I mentioned provides a way by using rational numbers as input and converting them to the destination type. See the update in the answer. - Gus

1 Answers

5
votes

Use LanguagePrimitives.GenericZero / GenericOne and let type inference do the rest

// attempt 1:
let inline stepFct1 x =
    let zero = LanguagePrimitives.GenericZero 
    if x > zero then x
    else zero

I had a look at the link you sent with the function you want to implement. FSharpPlus (F#+) may help you to write generic math code since it contains a dedicated Generic Numbers module. Or at least you can grab some techniques from there.

UPDATE

Regarding your updated question, which takes the complexity to a higher level, here is a solution using the latest version of the F#+ project:

let inline CDF(x:^T) : ^T = 
    let num x = fromRational (x </ratio/> 1000000000I)
    let (b1,b2,b3)  = (num 319381530I   , num -356563782I  , num 1781477937I)
    let (b4,b5)     = (num -1821255978I , num 1330274429I)
    let (p , c )    = (num  0231641900I , num 0398942280I)
    let (zero, one, two) = 0G, 1G, 2G
    if x > zero then
        let t = one / (one + p * x) 
        (one - c * exp( -x * x / two)* t * (t*(t*(t*(t*b5+b4)+b3)+b2)+b1)) 
    else
        let t = one / (one - p * x) 
        (c * exp( -x * x / two)* t * (t*(t*(t*(t*b5+b4)+b3)+b2)+b1))

Unfortunately at this time I realised some functions were marked as internal in the library and therefore were not exposed, but I re-created them in a working example here so you can test your function, which works nicely with float and float32.

A new version will be released before end of this year, but in the mean time you can branch it, remove the internals and compile it, or just re-create the functions as I did in the linked example.

If you are interested in Generic Maths feel free to contribute with code or use cases.