diff --git a/packages/repositories.config b/packages/repositories.config index 2cfb26c2..f8ef82c6 100644 --- a/packages/repositories.config +++ b/packages/repositories.config @@ -1,5 +1,6 @@  + \ No newline at end of file diff --git a/src/FSharp/DistributionMonad.fs b/src/FSharp/DistributionMonad.fs new file mode 100644 index 00000000..eb32c91b --- /dev/null +++ b/src/FSharp/DistributionMonad.fs @@ -0,0 +1,98 @@ +namespace MathNet.Numerics + +#nowarn "40" + +open System +open System.Collections +open System.Collections.Generic + +module Distribution = + + type 'a Outcome = { + Value: 'a + Probability : BigRational } + + type 'a Distribution = 'a Outcome seq + + // P(A AND B) = P(A | B) * P(B) + let bind (f: 'a -> 'b Distribution) (dist:'a Distribution) = + dist + |> Seq.map (fun p1 -> + f p1.Value + |> Seq.map (fun p2 -> + { Value = p2.Value; + Probability = + p1.Probability * p2.Probability})) + |> Seq.concat : 'b Distribution + + /// Sequentially compose two actions, passing any value produced by the first as an argument to the second. + let inline (>>=) dist f = bind f dist + /// Flipped >>= + let inline (=<<) f dist = bind f dist + + /// Inject a value into the Distribution type + let returnM (value:'a) = + Seq.singleton { Value = value ; Probability = 1N/1N } + : 'a Distribution + + type DistributionMonadBuilder() = + member this.Bind (r, f) = bind f r + member this.Return x = returnM x + member this.ReturnFrom x = x + + let distribution = DistributionMonadBuilder() + + // Create some helpers + let toUniformDistribution seq : 'a Distribution = + let l = Seq.length seq + seq + |> Seq.map (fun e -> + { Value = e; + Probability = 1N / bignum.FromInt l }) + + let probability (dist:'a Distribution) = + dist + |> Seq.map (fun o -> o.Probability) + |> Seq.sum + + let certainly = returnM + let impossible<'a> :'a Distribution = toUniformDistribution [] + + let fairDice sides = toUniformDistribution [1..sides] + + type CoinSide = + | Heads + | Tails + + let fairCoin = toUniformDistribution [Heads; Tails] + + let filter predicate (dist:'a Distribution) : 'a Distribution = + dist |> Seq.filter (fun o -> predicate o.Value) + + let filterInAnyOrder items dist = + items + |> Seq.fold (fun d item -> filter (Seq.exists ((=) (item))) d) dist + + /// Transforms a Distribution value by using a specified mapping function. + let map f (dist:'a Distribution) : 'b Distribution = + dist + |> Seq.map (fun o -> { Value = f o.Value; Probability = o.Probability }) + + let selectOne values = + [for e in values -> e,values |> Seq.filter ((<>) e)] + |> toUniformDistribution + + let rec selectMany n values = + match n with + | 0 -> certainly ([],values) + | _ -> + distribution { + let! (x,c1) = selectOne values + let! (xs,c2) = selectMany (n-1) c1 + return x::xs,c2} + + let select n values = + selectMany n values + |> map (fst >> List.rev) + + let remove items = Seq.filter (fun v -> Seq.forall ((<>) v) items) diff --git a/src/FSharp/FSharp.fsproj b/src/FSharp/FSharp.fsproj index c8092380..26d2f50d 100644 --- a/src/FSharp/FSharp.fsproj +++ b/src/FSharp/FSharp.fsproj @@ -57,6 +57,11 @@ + + + + + diff --git a/src/FSharp/complex.fs b/src/FSharp/complex.fs new file mode 100644 index 00000000..65707075 --- /dev/null +++ b/src/FSharp/complex.fs @@ -0,0 +1,133 @@ +// First version copied from the F# Power Pack +// https://raw.github.com/fsharp/powerpack/master/src/FSharp.PowerPack/math/complex.fs + +// (c) Microsoft Corporation 2005-2009. + +#nowarn "52" // defensive copy of structs warning + +namespace MathNet.Numerics + + open Microsoft.FSharp.Math + open System + open System.Globalization + + [] + [] + type Complex(real: float, imaginary: float) = + //new() = new Complex(0.0,0.0) + member x.r = real + member x.i = imaginary + override x.ToString() = x.ToString("g") + member x.ToString(fmt) = x.ToString(fmt,CultureInfo.InvariantCulture) + member x.ToString(fmt,fmtprovider:IFormatProvider) = + x.r.ToString(fmt,fmtprovider)+"r"+(if x.i < 0.0 then "-" else "+")+(System.Math.Abs x.i).ToString(fmt,fmtprovider)+"i" + interface IComparable with + member x.CompareTo(obj) = + match obj with + | :? Complex as y -> + let c = compare x.r y.r + if c <> 0 then c else compare x.i y.i + | _ -> invalidArg "obj" "not a Complex number" + override x.Equals(obj) = + match obj with + | :? Complex as y -> x.r = y.r && x.i = y.i + | _ -> false + override x.GetHashCode() = + (hash x.r >>> 5) ^^^ (hash x.r <<< 3) ^^^ (((hash x.i >>> 4) ^^^ (hash x.i <<< 4)) + 0x9e3779b9) + + + type complex = Complex + + [] + module Complex = + let mkRect(a,b) = new Complex(a,b) + let conjugate (c:complex) = mkRect (c.r, -c.i) + let mkPolar(a,b) = mkRect (a * Math.Cos(b), a * Math.Sin(b)) + let cis b = mkPolar(1.0,b) + let zero = mkRect(0.,0.) + let one = mkRect(1.,0.) + let onei = mkRect(0.,1.) + let magnitude (c:complex) = sqrt(c.r*c.r + c.i*c.i) + let phase (c:complex) = Math.Atan2(c.i,c.r) + let realPart (c:complex) = c.r + let imagPart (c:complex) = c.i + let abs (a:complex) = sqrt (a.r**2.0 + a.i**2.0) + let add (a:complex) (b:complex) = mkRect(a.r + b.r, a.i+b.i) + let sub (a:complex) (b:complex) = mkRect(a.r - b.r, a.i-b.i) + let mul (a:complex) (b:complex) = mkRect(a.r * b.r - a.i * b.i, a.i*b.r + b.i*a.r) + let div (x:complex) (y:complex) = + let a = x.r in let b = x.i in + let c = y.r in let d = y.i in + //(a+ib)/(c+id)=(ac+bd+i(bc-ad))/(c2+d2) + let q = c*c + d*d in + mkRect((a*c+b*d)/q, (b*c - a*d)/q) + let neg (a:complex) = mkRect(-a.r,-a.i) + let smul (a:float)(b:complex) = mkRect(a * b.r, a*b.i) + let muls (a:complex) (b:float) = mkRect(a.r *b, a.i*b) + let fmt_of_string numstyle fmtprovider (s:string) = + mkRect (System.Double.Parse(s,numstyle,fmtprovider),0.0) + let of_string s = fmt_of_string NumberStyles.Any CultureInfo.InvariantCulture s + + // ik.(r + i.th) = -k.th + i.k.r + let iscale k (x:complex) = mkRect (-k * x.i , k * x.r) + + // LogN : 'a * 'a -> 'a + // Asin : 'a -> 'a + // Acos : 'a -> 'a + // Atan : 'a -> 'a + // Atan2 : 'a * 'a -> 'a + // Sinh : 'a -> 'a + // Cosh : 'a -> 'a + // Tanh : 'a -> 'a + + let pi = mkRect (Math.PI,0.0) + + // exp(r+it) = exp(r).(cos(t)+i.sin(t)) - De Moivre Theorem + let exp (x:complex) = smul (exp(x.r)) (mkRect(cos(x.i), sin(x.i))) + // x = mag.e^(i.th) = e^ln(mag).e^(i.th) = e^(ln(mag) + i.th) + let log x = mkRect (log(magnitude(x)),phase(x)) + + let sqrt x = mkPolar (sqrt(magnitude x),phase x / 2.0) + + // cos(x) = (exp(i.x) + exp(-i.x))/2 + let cos x = smul 0.5 (add (exp(iscale 1.0 x)) (exp(iscale -1.0 x))) + // sin(x) = (exp(i.x) - exp(-i.x))/2 . (-i) + let sin x = smul 0.5 (sub (exp(iscale 1.0 x)) (exp(iscale -1.0 x))) |> iscale (-1.0) + // tan(x) = (exp(i.x) - exp(-i.x)) . (-i) / (exp(i.x) + exp(-i.x)) + // = (exp(2i.x) - 1.0) . (-i) / (exp(2i.x) + 1.0) + let tan x = let exp2ix = exp(iscale 2.0 x) in + (div (sub exp2ix one) (add exp2ix one)) |> iscale -1.0 + + + type Complex with + static member Create(a,b) = Complex.mkRect (a,b) + static member CreatePolar(a,b) = Complex.mkPolar (a,b) + member x.Magnitude = Complex.magnitude x + member x.Phase = Complex.phase x + member x.RealPart = x.r + member x.ImaginaryPart = x.i + member x.Conjugate = Complex.conjugate x + + static member Sin(x) = Complex.sin(x) + static member Cos(x) = Complex.cos(x) + static member Abs(x) = Complex.abs(x) + static member Tan(x) = Complex.tan(x) + static member Log(x) = Complex.log(x) + static member Exp(x) = Complex.exp(x) + static member Sqrt(x) = Complex.sqrt(x) + + static member Zero = Complex.zero + static member One = Complex.one + static member OneI = Complex.onei + static member ( + ) (a,b) = Complex.add a b + static member ( - ) (a,b) = Complex.sub a b + static member ( * ) (a,b) = Complex.mul a b + static member ( / ) (a,b) = Complex.div a b + static member ( ~- ) a = Complex.neg a + static member ( * ) (a,b) = Complex.smul a b + static member ( * ) (a,b) = Complex.muls a b + + + module ComplexTopLevelOperators = + let complex x y = Complex.mkRect (x,y) + diff --git a/src/FSharp/complex.fsi b/src/FSharp/complex.fsi new file mode 100644 index 00000000..185d9c5e --- /dev/null +++ b/src/FSharp/complex.fsi @@ -0,0 +1,139 @@ +// First version copied from the F# Power Pack +// https://raw.github.com/fsharp/powerpack/master/src/FSharp.PowerPack/math/complex.fsi + +// (c) Microsoft Corporation 2005-2009. + +namespace MathNet.Numerics + + open System + + /// The type of complex numbers stored as pairs of 64-bit floating point numbers in rectangular coordinates + [] + [] + type Complex = + /// The real part of a complex number + member r: float + /// The imaginary part of a complex number + member i: float + /// The polar-coordinate magnitude of a complex number + member Magnitude: float + /// The polar-coordinate phase of a complex number + member Phase: float + /// The real part of a complex number + member RealPart: float + /// The imaginary part of a complex number + member ImaginaryPart: float + /// The conjugate of a complex number, i.e. x-yi + member Conjugate: Complex + /// Create a complex number x+ij using rectangular coordinates + static member Create : float * float -> Complex + /// Create a complex number using magnitude/phase polar coordinates + static member CreatePolar : float * float -> Complex + /// The complex number 0+0i + static member Zero : Complex + /// The complex number 1+0i + static member One : Complex + /// The complex number 0+1i + static member OneI : Complex + /// Add two complex numbers + static member ( + ) : Complex * Complex -> Complex + /// Subtract one complex number from another + static member ( - ) : Complex * Complex -> Complex + /// Multiply two complex numbers + static member ( * ) : Complex * Complex -> Complex + /// Complex division of two complex numbers + static member ( / ) : Complex * Complex -> Complex + /// Unary negation of a complex number + static member ( ~- ) : Complex -> Complex + /// Multiply a scalar by a complex number + static member ( * ) : float * Complex -> Complex + /// Multiply a complex number by a scalar + static member ( * ) : Complex * float -> Complex + + static member Sin : Complex -> Complex + static member Cos : Complex -> Complex + + /// Computes the absolute value of a complex number: e.g. Abs x+iy = sqrt(x**2.0 + y**2.0.) + /// Note: Complex.Abs(z) is the same as z.Magnitude + static member Abs : Complex -> float + static member Tan : Complex -> Complex + static member Log : Complex -> Complex + static member Exp : Complex -> Complex + static member Sqrt : Complex -> Complex + + override ToString : unit -> string + override Equals : obj -> bool + interface System.IComparable + member ToString : format:string -> string + member ToString : format:string * provider:System.IFormatProvider -> string + + /// The type of complex numbers + type complex = Complex + + + [] + [] + module Complex = + + val mkRect: float * float -> complex + + /// The polar-coordinate magnitude of a complex number + val magnitude: complex -> float + /// The polar-coordinate phase of a complex number + val phase : complex -> float + /// The real part of a complex number + val realPart : complex -> float + /// The imaginary part of a complex number + val imagPart : complex -> float + /// Create a complex number using magnitude/phase polar coordinates + val mkPolar : float * float -> complex + /// A complex of magnitude 1 and the given phase and , i.e. cis x = mkPolar 1.0 x + val cis : float -> complex + + /// The conjugate of a complex number, i.e. x-yi + val conjugate : complex -> complex + + /// The complex number 0+0i + val zero : complex + /// The complex number 1+0i + val one : complex + /// The complex number 0+1i + val onei : complex + /// Add two complex numbers + val add : complex -> complex -> complex + /// Subtract one complex number from another + val sub : complex -> complex -> complex + /// Multiply two complex numbers + val mul : complex -> complex -> complex + /// Complex division of two complex numbers + val div : complex -> complex -> complex + /// Unary negation of a complex number + val neg : complex -> complex + /// Multiply a scalar by a complex number + val smul : float -> complex -> complex + /// Multiply a complex number by a scalar + val muls : complex -> float -> complex + + /// pi + val pi : Complex + /// exp(x) = e^x + val exp : Complex -> Complex + /// log(x) is natural log (base e) + val log : Complex -> Complex + /// sqrt(x) and 0 <= phase(x) < pi + val sqrt : Complex -> Complex + /// Sine + val sin : Complex -> Complex + /// Cosine + val cos : Complex -> Complex + /// Tagent + val tan : Complex -> Complex + + + [] + module ComplexTopLevelOperators = + /// Constructs a complex number from both the real and imaginary part. + val complex : float -> float -> complex + + + diff --git a/src/FSharp/q.fs b/src/FSharp/q.fs new file mode 100644 index 00000000..2cb24129 --- /dev/null +++ b/src/FSharp/q.fs @@ -0,0 +1,309 @@ +// First version copied from the F# Power Pack +// https://raw.github.com/fsharp/powerpack/master/src/FSharp.PowerPack/math/q.fs + + +// (c) Microsoft Corporation. All rights reserved + +#nowarn "44" // OK to use the "compiler only" function RangeGeneric +#nowarn "52" // The value has been copied to ensure the original is not mutated by this operation + +namespace MathNet.Numerics + + open System + open System.Numerics + open System.Globalization + + module BigRationalLargeImpl = + let ZeroI = new BigInteger(0) + let OneI = new BigInteger(1) + let bigint (x:int) = new BigInteger(x) + let ToDoubleI (x:BigInteger) = double x + let ToInt32I (x:BigInteger) = int32 x + + open BigRationalLargeImpl + + [] + type BigRationalLarge = + | Q of BigInteger * BigInteger // invariants: (p,q) in lowest form, q >= 0 + + override n.ToString() = + let (Q(p,q)) = n + if q.IsOne then p.ToString() + else p.ToString() + "/" + q.ToString() + + + static member Hash (Q(ap,aq)) = + // This hash code must be identical to the hash for BigInteger when the numbers coincide. + if aq.IsOne then ap.GetHashCode() else (ap.GetHashCode() <<< 3) + aq.GetHashCode() + + + override x.GetHashCode() = BigRationalLarge.Hash(x) + + static member Equals(Q(ap,aq), Q(bp,bq)) = + BigInteger.(=) (ap,bp) && BigInteger.(=) (aq,bq) // normal form, so structural equality + + static member LessThan(Q(ap,aq), Q(bp,bq)) = + BigInteger.(<) (ap * bq,bp * aq) + + // note: performance improvement possible here + static member Compare(p,q) = + if BigRationalLarge.LessThan(p,q) then -1 + elif BigRationalLarge.LessThan(q,p)then 1 + else 0 + + interface System.IComparable with + member this.CompareTo(obj:obj) = + match obj with + | :? BigRationalLarge as that -> BigRationalLarge.Compare(this,that) + | _ -> invalidArg "obj" "the object does not have the correct type" + + override this.Equals(that:obj) = + match that with + | :? BigRationalLarge as that -> BigRationalLarge.Equals(this,that) + | _ -> false + + member x.IsNegative = let (Q(ap,_)) = x in sign ap < 0 + member x.IsPositive = let (Q(ap,_)) = x in sign ap > 0 + + member x.Numerator = let (Q(p,_)) = x in p + member x.Denominator = let (Q(_,q)) = x in q + member x.Sign = (let (Q(p,_)) = x in sign p) + + static member ToDouble (Q(p,q)) = + ToDoubleI p / ToDoubleI q + + static member Normalize (p:BigInteger,q:BigInteger) = + if q.IsZero then + raise (System.DivideByZeroException()) (* throw for any x/0 *) + elif q.IsOne then + Q(p,q) + else + let k = BigInteger.GreatestCommonDivisor(p,q) + let p = p / k + let q = q / k + if sign q < 0 then Q(-p,-q) else Q(p,q) + + static member Rational (p:int,q:int) = BigRationalLarge.Normalize (bigint p,bigint q) + static member RationalZ (p,q) = BigRationalLarge.Normalize (p,q) + + static member Parse (str:string) = + let len = str.Length + if len=0 then invalidArg "str" "empty string"; + let j = str.IndexOf '/' + if j >= 0 then + let p = BigInteger.Parse (str.Substring(0,j)) + let q = BigInteger.Parse (str.Substring(j+1,len-j-1)) + BigRationalLarge.RationalZ (p,q) + else + let p = BigInteger.Parse str + BigRationalLarge.RationalZ (p,OneI) + + static member (~-) (Q(bp,bq)) = Q(-bp,bq) // still coprime, bq >= 0 + static member (+) (Q(ap,aq),Q(bp,bq)) = BigRationalLarge.Normalize ((ap * bq) + (bp * aq),aq * bq) + static member (-) (Q(ap,aq),Q(bp,bq)) = BigRationalLarge.Normalize ((ap * bq) - (bp * aq),aq * bq) + static member (*) (Q(ap,aq),Q(bp,bq)) = BigRationalLarge.Normalize (ap * bp,aq * bq) + static member (/) (Q(ap,aq),Q(bp,bq)) = BigRationalLarge.Normalize (ap * bq,aq * bp) + static member ( ~+ )(n1:BigRationalLarge) = n1 + + + [] + module BigRationalLarge = + open System.Numerics + + let inv (Q(ap,aq)) = BigRationalLarge.Normalize(aq,ap) + + let pown (Q(p,q)) (n:int) = Q(BigInteger.Pow(p,n),BigInteger.Pow (q,n)) // p,q powers still coprime + + let equal (Q(ap,aq)) (Q(bp,bq)) = ap=bp && aq=bq // normal form, so structural equality + let lt a b = BigRationalLarge.LessThan(a,b) + let gt a b = BigRationalLarge.LessThan(b,a) + let lte (Q(ap,aq)) (Q(bp,bq)) = BigInteger.(<=) (ap * bq,bp * aq) + let gte (Q(ap,aq)) (Q(bp,bq)) = BigInteger.(>=) (ap * bq,bp * aq) + + let of_bigint z = BigRationalLarge.RationalZ(z,OneI ) + let of_int n = BigRationalLarge.Rational(n,1) + + // integer part + let integer (Q(p,q)) = + let mutable r = BigInteger(0) + let d = BigInteger.DivRem (p,q,&r) // have p = d.q + r, |r| < |q| + if r < ZeroI + then d - OneI // p = (d-1).q + (r+q) + else d // p = d.q + r + + + //---------------------------------------------------------------------------- + // BigRational + //-------------------------------------------------------------------------- + + [] + [] + type BigRational = + | Z of BigInteger + | Q of BigRationalLarge + + static member ( + )(n1,n2) = + match n1,n2 with + | Z z ,Z zz -> Z (z + zz) + | Q q ,Q qq -> Q (q + qq) + | Z z ,Q qq -> Q (BigRationalLarge.of_bigint z + qq) + | Q q ,Z zz -> Q (q + BigRationalLarge.of_bigint zz) + + static member ( * )(n1,n2) = + match n1,n2 with + | Z z ,Z zz -> Z (z * zz) + | Q q ,Q qq -> Q (q * qq) + | Z z ,Q qq -> Q (BigRationalLarge.of_bigint z * qq) + | Q q ,Z zz -> Q (q * BigRationalLarge.of_bigint zz) + + static member ( - )(n1,n2) = + match n1,n2 with + | Z z ,Z zz -> Z (z - zz) + | Q q ,Q qq -> Q (q - qq) + | Z z ,Q qq -> Q (BigRationalLarge.of_bigint z - qq) + | Q q ,Z zz -> Q (q - BigRationalLarge.of_bigint zz) + + static member ( / )(n1,n2) = + match n1,n2 with + | Z z ,Z zz -> Q (BigRationalLarge.RationalZ(z,zz)) + | Q q ,Q qq -> Q (q / qq) + | Z z ,Q qq -> Q (BigRationalLarge.of_bigint z / qq) + | Q q ,Z zz -> Q (q / BigRationalLarge.of_bigint zz) + + static member ( ~- )(n1) = + match n1 with + | Z z -> Z (-z) + | Q q -> Q (-q) + + static member ( ~+ )(n1:BigRational) = n1 + + // nb. Q and Z hash codes must match up - see notes above + override n.GetHashCode() = + match n with + | Z z -> z.GetHashCode() + | Q q -> q.GetHashCode() + + override this.Equals(obj:obj) = + match obj with + | :? BigRational as that -> BigRational.(=)(this, that) + | _ -> false + + interface System.IComparable with + member n1.CompareTo(obj:obj) = + match obj with + | :? BigRational as n2 -> + if BigRational.(<)(n1, n2) then -1 elif BigRational.(=)(n1, n2) then 0 else 1 + | _ -> invalidArg "obj" "the objects are not comparable" + + static member FromInt (x:int) = Z (bigint x) + static member FromBigInt x = Z x + + static member Zero = BigRational.FromInt(0) + static member One = BigRational.FromInt(1) + + + static member PowN (n,i:int) = + match n with + | Z z -> Z (BigInteger.Pow (z,i)) + | Q q -> Q (BigRationalLarge.pown q i) + + static member op_Equality (n,nn) = + match n,nn with + | Z z ,Z zz -> BigInteger.(=) (z,zz) + | Q q ,Q qq -> (BigRationalLarge.equal q qq) + | Z z ,Q qq -> (BigRationalLarge.equal (BigRationalLarge.of_bigint z) qq) + | Q q ,Z zz -> (BigRationalLarge.equal q (BigRationalLarge.of_bigint zz)) + static member op_Inequality (n,nn) = not (BigRational.op_Equality(n,nn)) + + static member op_LessThan (n,nn) = + match n,nn with + | Z z ,Z zz -> BigInteger.(<) (z,zz) + | Q q ,Q qq -> (BigRationalLarge.lt q qq) + | Z z ,Q qq -> (BigRationalLarge.lt (BigRationalLarge.of_bigint z) qq) + | Q q ,Z zz -> (BigRationalLarge.lt q (BigRationalLarge.of_bigint zz)) + static member op_GreaterThan (n,nn) = + match n,nn with + | Z z ,Z zz -> BigInteger.(>) (z,zz) + | Q q ,Q qq -> (BigRationalLarge.gt q qq) + | Z z ,Q qq -> (BigRationalLarge.gt (BigRationalLarge.of_bigint z) qq) + | Q q ,Z zz -> (BigRationalLarge.gt q (BigRationalLarge.of_bigint zz)) + static member op_LessThanOrEqual (n,nn) = + match n,nn with + | Z z ,Z zz -> BigInteger.(<=) (z,zz) + | Q q ,Q qq -> (BigRationalLarge.lte q qq) + | Z z ,Q qq -> (BigRationalLarge.lte (BigRationalLarge.of_bigint z) qq) + | Q q ,Z zz -> (BigRationalLarge.lte q (BigRationalLarge.of_bigint zz)) + static member op_GreaterThanOrEqual (n,nn) = + match n,nn with + | Z z ,Z zz -> BigInteger.(>=) (z,zz) + | Q q ,Q qq -> (BigRationalLarge.gte q qq) + | Z z ,Q qq -> (BigRationalLarge.gte (BigRationalLarge.of_bigint z) qq) + | Q q ,Z zz -> (BigRationalLarge.gte q (BigRationalLarge.of_bigint zz)) + + + member n.IsNegative = + match n with + | Z z -> sign z < 0 + | Q q -> q.IsNegative + + member n.IsPositive = + match n with + | Z z -> sign z > 0 + | Q q -> q.IsPositive + + member n.Numerator = + match n with + | Z z -> z + | Q q -> q.Numerator + + member n.Denominator = + match n with + | Z _ -> OneI + | Q q -> q.Denominator + + member n.Sign = + if n.IsNegative then -1 + elif n.IsPositive then 1 + else 0 + + static member Abs(n:BigRational) = + if n.IsNegative then -n else n + + static member ToDouble(n:BigRational) = + match n with + | Z z -> ToDoubleI z + | Q q -> BigRationalLarge.ToDouble q + + static member ToBigInt(n:BigRational) = + match n with + | Z z -> z + | Q q -> BigRationalLarge.integer q + + static member ToInt32(n:BigRational) = + match n with + | Z z -> ToInt32I(z) + | Q q -> ToInt32I(BigRationalLarge.integer q ) + + static member op_Explicit (n:BigRational) = BigRational.ToInt32 n + static member op_Explicit (n:BigRational) = BigRational.ToDouble n + static member op_Explicit (n:BigRational) = BigRational.ToBigInt n + + + override n.ToString() = + match n with + | Z z -> z.ToString() + | Q q -> q.ToString() + + member x.StructuredDisplayString = x.ToString() + + static member Parse(s:string) = Q (BigRationalLarge.Parse s) + + type BigNum = BigRational + type bignum = BigNum + + module NumericLiteralN = + let FromZero () = BigRational.Zero + let FromOne () = BigRational.One + let FromInt32 i = BigRational.FromInt i + let FromInt64 (i64:int64) = BigRational.FromBigInt (new BigInteger(i64)) + let FromString s = BigRational.Parse s diff --git a/src/FSharp/q.fsi b/src/FSharp/q.fsi new file mode 100644 index 00000000..ffcf84e2 --- /dev/null +++ b/src/FSharp/q.fsi @@ -0,0 +1,96 @@ +// First version copied from the F# Power Pack +// https://raw.github.com/fsharp/powerpack/master/src/FSharp.PowerPack/math/q.fsi + + +// (c) Microsoft Corporation 2005-2009. + +namespace MathNet.Numerics + + open System + open System.Numerics + + /// The type of arbitrary-sized rational numbers + [] + type BigRational = + /// Return the sum of two rational numbers + static member ( + ) : BigRational * BigRational -> BigRational + /// Return the product of two rational numbers + static member ( * ) : BigRational * BigRational -> BigRational + /// Return the difference of two rational numbers + static member ( - ) : BigRational * BigRational -> BigRational + /// Return the ratio of two rational numbers + static member ( / ) : BigRational * BigRational -> BigRational + /// Return the negation of a rational number + static member ( ~- ): BigRational -> BigRational + /// Return the given rational number + static member ( ~+ ): BigRational -> BigRational + + override ToString: unit -> string + override GetHashCode: unit -> int + interface System.IComparable + + /// Get zero as a rational number + static member Zero : BigRational + /// Get one as a rational number + static member One : BigRational + /// This operator is for use from other .NET languages + static member op_Equality : BigRational * BigRational -> bool + /// This operator is for use from other .NET languages + static member op_Inequality : BigRational * BigRational -> bool + /// This operator is for use from other .NET languages + static member op_LessThan: BigRational * BigRational -> bool + /// This operator is for use from other .NET languages + static member op_GreaterThan: BigRational * BigRational -> bool + /// This operator is for use from other .NET languages + static member op_LessThanOrEqual: BigRational * BigRational -> bool + /// This operator is for use from other .NET languages + static member op_GreaterThanOrEqual: BigRational * BigRational -> bool + + /// Return a boolean indicating if this rational number is strictly negative + member IsNegative: bool + /// Return a boolean indicating if this rational number is strictly positive + member IsPositive: bool + + /// Return the numerator of the normalized rational number + member Numerator: BigInteger + /// Return the denominator of the normalized rational number + member Denominator: BigInteger + + member StructuredDisplayString : string + + /// Return the absolute value of a rational number + static member Abs : BigRational -> BigRational + /// Return the sign of a rational number; 0, +1 or -1 + member Sign : int + /// Return the result of raising the given rational number to the given power + static member PowN : BigRational * int -> BigRational + /// Return the result of converting the given integer to a rational number + static member FromInt : int -> BigRational + /// Return the result of converting the given big integer to a rational number + static member FromBigInt : BigInteger -> BigRational + /// Return the result of converting the given rational number to a floating point number + static member ToDouble: BigRational -> float + /// Return the result of converting the given rational number to a big integer + static member ToBigInt: BigRational -> BigInteger + /// Return the result of converting the given rational number to an integer + static member ToInt32 : BigRational -> int + /// Return the result of converting the given rational number to a floating point number + static member op_Explicit : BigRational -> float + /// Return the result of converting the given rational number to a big integer + static member op_Explicit : BigRational -> BigInteger + /// Return the result of converting the given rational number to an integer + static member op_Explicit : BigRational -> int + /// Return the result of converting the string to a rational number + static member Parse: string -> BigRational + + type BigNum = BigRational + + type bignum = BigRational + + [] + module NumericLiteralN = + val FromZero : unit -> BigRational + val FromOne : unit -> BigRational + val FromInt32 : int32 -> BigRational + val FromInt64 : int64 -> BigRational + val FromString : string -> BigRational \ No newline at end of file diff --git a/src/FSharpUnitTests/BigRationalTests.fs b/src/FSharpUnitTests/BigRationalTests.fs new file mode 100644 index 00000000..1173e196 --- /dev/null +++ b/src/FSharpUnitTests/BigRationalTests.fs @@ -0,0 +1,614 @@ +// First version copied from the F# Power Pack +// https://raw.github.com/fsharp/powerpack/master/src/FSharp.PowerPack.Unittests/BigRationalTests.fs + +namespace MathNet.Numerics.Tests + +open MathNet.Numerics +open NUnit.Framework +open System +open System.Collections +open System.Collections.Generic +open System.Numerics + + +[] +type public BigRationalTests() = + + // BigRational Tests + // ================= + + // Notes: What cases to consider? + // For (p,q) cases q=0, q=1, q<>1. [UPDATE: remove (x,0)] + // For (p,q) when q=1 there could be 2 internal representations, either Z or Q. + // For (p,0) this value can be signed, corresponds to +/- infinity point. [Update: remove it] + // Hashes on (p,1) for both representations must agree. + // For binary operators, try for result with and without HCF (normalisation). + // Also: 0/0 is an acceptable representation. See normalisation code. [Update: remove it]. + + // Overrides to test: + // .ToString() + // .GetHashCode() + // .Equals() + // IComparable.CompareTo() + + // Misc construction. + let natA n = BigRational.FromInt n // internally Z + let natB n = (natA n / natA 7) * natA 7 // internally Q + let ratio p q = BigRational.FromInt p / BigRational.FromInt q + let (/%) b c = BigRational.FromBigInt b / BigRational.FromBigInt c + + // Misc test values + let q0 = natA 0 + let q1 = natA 1 + let q2 = natA 2 + let q3 = natA 3 + let q4 = natA 4 + let q5 = natA 5 + let minIntI = bigint System.Int32.MinValue + let maxIntI = bigint System.Int32.MaxValue + let ran = System.Random() + let nextZ n = bigint (ran.Next(n)) + + // A selection of test points. + let points = + // A selection of integer and reciprical points + let points = + [for i in -13I .. 13I -> i,1I] @ + [for i in -13I .. 13I -> 1I,i] + // Exclude x/0 + let points = [for p,q in points do if q <> 0I then yield p,q ] // PROPOSE: (q,0) never a valid Q value, filter them out of tests... + // Scale by various values, including into BigInt range + let scale (kp,kq) (p,q) = (p*kp,q*kq) + let scales k pqs = List.map (scale k) pqs + let points = List.concat [points; + scales (10000I,1I) points; + scales (1I,10000I) points; + scales (maxIntI,1I) points; + scales (1I,maxIntI) points; + ] + points + let pointsNonZero = [for p,q in points do if p<>0I then yield p,q] // non zero points + + let makeQs p q = + if q = 1I && minIntI <= p && p <= maxIntI then + // (p,1) where p is int32 + let p32 = int32 p + [natA p32;natB p32;BigRational.FromBigInt p] // two reprs for int32 + else + [BigRational.FromBigInt p / BigRational.FromBigInt q] + + let miscQs = [for p,q in points do yield! makeQs p q] + + let product xs ys = [for x in xs do for y in ys do yield x,y] + let vector1s = [for z in points -> z] + let vector2s = product points points + + [] + member this.BasicTests1() = + check "generic format h" "1N" (sprintf "%A" 1N) + check "generic format q" "-1N" (sprintf "%A" (-1N)) + + test "vliwe98" (id -2N = - 2N) + test "d3oc002" (LanguagePrimitives.GenericZero = 0N) + test "d3oc112w" (LanguagePrimitives.GenericOne = 1N) + + check "weioj3h" (sprintf "%O" 3N) "3" + check "weioj3k" (sprintf "%O" (3N / 4N)) "3/4" + check "weioj3k" (sprintf "%O" (3N / 400000000N)) "3/400000000" + check "weioj3l" (sprintf "%O" (3N / 3N)) "1" + check "weioj3q" (sprintf "%O" (-3N)) "-3" + //check "weioj3w" (sprintf "%O" -3N) "-3" + check "weioj3e" (sprintf "%O" (-3N / -3N)) "1" + + // The reason why we do not use hardcoded values is the the representation may change based on the NetFx we are targeting. + // For example, when targeting NetFx4.0, the result is "-3E+61" instead of "-3000....0N" + let v = -30000000000000000000000000000000000000000000000000000000000000N + check "weioj3r" (sprintf "%O" v) ((box v).ToString()) + + + [] + member this.BasicTests2() = + + + // Test arithmetic ops: tests + let test2One name f check ((p,q),(pp,qq)) = + // There may be several ways to construct the test rationals + let zs = makeQs p q + let zzs = makeQs pp qq + let results = [for z in zs do for zz in zzs do yield f (z,zz)] + let refP,refQ = check (p,q) (pp,qq) + let refResult = BigRational.FromBigInt refP / BigRational.FromBigInt refQ + let resOK (result:BigRational) = + result.Numerator * refQ = refP * result.Denominator && + BigRational.Equals(refResult,result) + match List.tryFind (fun result -> not (resOK result)) results with + | None -> () // ok + | Some result -> printf "Test failed. %s (%A,%A) (%A,%A). Expected %A. Observed %A\n" name p q pp qq refResult result + reportFailure "cejkew09" + + let test2All name f check vectors = List.iter (test2One name f check) vectors + + // Test arithmetic ops: call + test2All "add" (BigRational.(+)) (fun (p,q) (pp,qq) -> (p*qq + q*pp,q*qq)) vector2s + test2All "sub" (BigRational.(-)) (fun (p,q) (pp,qq) -> (p*qq - q*pp,q*qq)) vector2s + test2All "mul" (BigRational.(*)) (fun (p,q) (pp,qq) -> (p*pp,q*qq)) vector2s // *) <-- for EMACS + test2All "div" (BigRational.(/)) (fun (p,q) (pp,qq) -> (p*qq,q*pp)) (product points pointsNonZero) + + + + [] + member this.RangeTests() = + // Test x0 .. dx .. x1 + let checkRange3 (x0:BigRational) dx x1 k = + let f (x:BigRational) = x * BigRational.FromBigInt k |> BigRational.ToBigInt + let rangeA = {x0 .. dx .. x1} |> Seq.map f + let rangeB = {f x0 .. f dx .. f x1} + //printf "Length=%d\n" (Seq.length rangeA) + let same = Seq.forall2 (=) rangeA rangeB + check (sprintf "Range3 %A .. %A .. %A scaled to %A" x0 dx x1 k) same true + + checkRange3 (0I /% 1I) (1I /% 7I) (100I /% 1I) (7I*1I) + checkRange3 (0I /% 1I) (1I /% 7I) (100I /% 11I) (7I*11I) + checkRange3 (1I /% 13I) (1I /% 7I) (100I /% 11I) (7I*11I*13I) + for i = 0 to 1000 do + let m = 1000 // max steps is -m to m in steps of 1/m i.e. 2.m^2 + let p0,q0 = nextZ m ,nextZ m + 1I + let p1,q1 = nextZ m ,nextZ m + 1I + let pd,qd = nextZ m + 1I,nextZ m + 1I + checkRange3 (p0 /% q0) (pd /% qd) (p1 /% q1) (q0 * q1 * qd) + + + // Test x0 .. x1 + let checkRange2 (x0:BigRational) x1 = + let z0 = BigRational.ToBigInt x0 + let z01 = BigRational.ToBigInt (x1 - x0) + let f (x:BigRational) = x |> BigRational.ToBigInt + let rangeA = [x0 .. x1] |> List.map f // range with each item rounded down + let rangeB = [z0 .. z0 + z01] // range of same length from the round down start point + check (sprintf "Range2: %A .. %A" x0 x1) rangeA rangeB + + checkRange2 (0I /% 1I) (100I /% 1I) + checkRange2 (0I /% 1I) (100I /% 11I) + checkRange2 (1I /% 13I) (100I /% 11I) + for i = 0 to 1000 do + let m = 10000 // max steps is -m to m in steps of 1 i.e. 2.m + let p0,q0 = nextZ m ,nextZ m + 1I + let p1,q1 = nextZ m ,nextZ m + 1I + checkRange2 (p0 /% q0) (p1 /% q1) //(q0 * q1 * qd) + + // ToString() + // Cases: integer, computed integer, rational<1, rational>1, +/-infinity, nan + (natA 1).ToString() |> check "ToString" "1" + (natA 0).ToString() |> check "ToString" "0" + (natA (-12)).ToString() |> check "ToString" "-12" + (natB 1).ToString() |> check "ToString" "1" + (natB 0).ToString() |> check "ToString" "0" + (natB (-12)).ToString() |> check "ToString" "-12" + (1I /% 3I).ToString() |> check "ToString" "1/3" + (12I /% 5I).ToString() |> check "ToString" "12/5" + //(13I /% 0I).ToString() |> check "ToString" "1/0" // + 1/0. Plan to make this invalid value + //(-13I /% 0I).ToString() |> check "ToString" "1/0" // - 1/0. Plan to make this invalid value + //(0I /% 0I).ToString() |> check "ToString" "0/0" // 0/0. Plan to make this invalid value + + // GetHashCode + // Cases: zero, integer, computed integer, computed by multiple routes. + let checkSameHashGeneric a b = check (sprintf "GenericHash %A %A" a b) (a.GetHashCode()) (b.GetHashCode()) + let checkSameHash (a:BigRational) (b:BigRational) = check (sprintf "BigRationalHash %A %A" a b) (a.GetHashCode()) (b.GetHashCode()); checkSameHashGeneric a b + + List.iter (fun n -> checkSameHash (natA n) (natB n)) [-10 .. 10] + List.iter (fun n -> checkSameHash n ((n * q3 + n * q2) / q5)) miscQs + + // bug 3488: should non-finite values be supported? + //let x = BigRational.FromBigInt (-1I) / BigRational.FromBigInt 0I + //let q2,q3,q5 = BigRational.FromInt 2,BigRational.FromInt 3,BigRational.FromInt 5 + //let x2 = (x * q2 + x * q3) / q5 + //x,x2,x = x2 + + // Test: Zero,One? + check "ZeroA" BigRational.Zero (natA 0) + check "ZeroA" BigRational.Zero (natA 0) + check "OneA" BigRational.One (natB 1) + check "OneB" BigRational.One (natB 1) + + [] + member this.BinaryAndUnaryOperators() = + // Test: generic bop + let testR2One name f check ((p,q),(pp,qq)) = + // There may be several ways to construct the test rationals + let zs = makeQs p q + let zzs = makeQs pp qq + let resultRef = check (p,q) (pp,qq) // : bool + let args = [for z in zs do for zz in zzs do yield (z,zz)] + match List.tryFind (fun (z,zz) -> resultRef <> f (z,zz)) args with + | None -> () // ok + | Some (z,zz) -> printf "Test failed. %s (%A,%A) (%A,%A) = %s %A %A. Expected %A.\n" name p q pp qq name z zz resultRef + reportFailure "cknwe9" + + // Test: generic uop + let testR1One name f check (p,q) = + // There may be several ways to construct the test rationals + let zs = makeQs p q + let resultRef = check (p,q) //: bool + match List.tryFind (fun z -> resultRef <> f z) zs with + | None -> () // ok + | Some z -> printf "Test failed. %s (%A,%A) = %s %A. Expected %A.\n" name p q name z resultRef + reportFailure "vekjkrejvre0" + + let testR2All name f check vectors = List.iter (testR2One name f check) vectors + let testR1All name f check vectors = List.iter (testR1One name f check) vectors + + // Test: relations + let sign (i:BigInteger) = BigInteger(i.Sign) + testR2All "=" BigRational.(=) (fun (p,q) (pp,qq) -> (p*qq = q*pp)) vector2s + testR2All "=" BigRational.op_Equality (fun (p,q) (pp,qq) -> (p*qq = q*pp)) vector2s + testR2All "!=" BigRational.op_Inequality (fun (p,q) (pp,qq) -> (p*qq <> q*pp)) vector2s + // p/q < pp/qq + // iff (p * sign q) / (q * sign q) < (pp * sign qq) / (qq * sign qq) + // iff (p * sign q) * (qq * sign qq) < (pp * sign qq) * (q * sign q) since q*sign q is always +ve. + testR2All "<" BigRational.(<) (fun (p,q) (pp,qq) -> (p * sign q) * (qq * sign qq) < (pp * sign qq) * (q * sign q)) vector2s + testR2All ">" BigRational.(>) (fun (p,q) (pp,qq) -> (p * sign q) * (qq * sign qq) > (pp * sign qq) * (q * sign q)) vector2s + testR2All "<=" BigRational.(<=) (fun (p,q) (pp,qq) -> (p * sign q) * (qq * sign qq) <= (pp * sign qq) * (q * sign q)) vector2s + testR2All ">=" BigRational.(>=) (fun (p,q) (pp,qq) -> (p * sign q) * (qq * sign qq) >= (pp * sign qq) * (q * sign q)) vector2s + + // System.IComparable tests + let BigRationalCompareTo (p:BigRational,q:BigRational) = (p :> System.IComparable).CompareTo(q) + testR2All "IComparable.CompareTo" BigRationalCompareTo (fun (p,q) (pp,qq) -> compare ((p * sign q) * (qq * sign qq)) ((pp * sign qq) * (q * sign q))) vector2s + + // Test: is negative, is positive + testR1All "IsNegative" (fun (x:BigRational) -> x.IsNegative) (fun (p,q) -> sign p * sign q = -1I) vector1s + testR1All "IsPositive" (fun (x:BigRational) -> x.IsPositive) (fun (p,q) -> sign p * sign q = 1I) vector1s + testR1All "IsZero" (fun (x:BigRational) -> x = q0) (fun (p,q) -> sign p = 0I) vector1s + + + let test1One name f check (p,q) = + // There may be several ways to construct the test rationals + let zs = makeQs p q + let results = [for z in zs -> f z] + let refP,refQ = check (p,q) + let refResult = BigRational.FromBigInt refP / BigRational.FromBigInt refQ + let resOK (result:BigRational) = + result.Numerator * refQ = refP * result.Denominator && + BigRational.Equals(refResult,result) + match List.tryFind (fun result -> not (resOK result)) results with + | None -> () // ok + | Some result -> printf "Test failed. %s (%A,%A). Expected %A. Observed %A\n" name p q refResult result + reportFailure "klcwe09wek" + + let test1All name f check vectors = List.iter (test1One name f check) vectors + + test1All "neg" (BigRational.(~-)) (fun (p,q) -> (-p,q)) vector1s + test1All "pos" (BigRational.(~+)) (fun (p,q) -> (p,q)) vector1s // why have ~+ ??? + + // Test: Abs,Sign + test1All "Abs" (BigRational.Abs) (fun (p,q) -> (abs p,abs q)) vector1s + testR1All "Sign" (fun (x:BigRational) -> x.Sign) (fun (p,q) -> check "NonZeroDenom" (sign q <> 0I) true; (sign p * sign q) |> int32) vector1s + + // Test: PowN + test1All "PowN(x,2)" (fun x -> BigRational.PowN(x,2)) (fun (p,q) -> (p*p,q*q)) vector1s + test1All "PowN(x,1)" (fun x -> BigRational.PowN(x,1)) (fun (p,q) -> (p,q)) vector1s + test1All "PowN(x,0)" (fun x -> BigRational.PowN(x,0)) (fun (p,q) -> (1I,1I)) vector1s + + // MatteoT: moved to numbersVS2008\test.ml + //test1All "PowN(x,200)" (fun x -> BigRational.PowN(x,200)) (fun (p,q) -> (BigInteger.Pow(p,200I),BigInteger.Pow(q,200I))) vector1s + + // MatteoT: moved to numbersVS2008\test.ml + //let powers = [0I .. 100I] + //powers |> List.iter (fun i -> test1All "PowN(x,i)" (fun x -> BigRational.PowN(x,int i)) (fun (p,q) -> (BigInteger.Pow(p,i),BigInteger.Pow(q,i))) vector1s) + + // Test: PowN with negative powers - expect exception + testR1All "PowN(x,-1)" (fun x -> throws (fun () -> BigRational.PowN(x,-1))) (fun (p,q) -> true) vector1s + testR1All "PowN(x,-4)" (fun x -> throws (fun () -> BigRational.PowN(x,-4))) (fun (p,q) -> true) vector1s + + + +[] +type BigNumType() = + let g_positive1 = 1000000000000000000000000000000000018N + let g_positive2 = 1000000000000000000000000000000000000N + let g_negative1 = -1000000000000000000000000000000000018N + let g_negative2 = -1000000000000000000000000000000000000N + let g_negative3 = -1000000000000000000000000000000000036N + let g_zero = 0N + let g_normal = 88N + let g_bigintpositive = 1000000000000000000000000000000000018I + let g_bigintnegative = -1000000000000000000000000000000000018I + + // Interfaces + [] + member this.IComparable() = + // Legit IC + let ic = g_positive1 :> IComparable + Assert.AreEqual(ic.CompareTo(g_positive1),0) + checkThrowsArgumentException( fun () -> ic.CompareTo(g_bigintpositive) |> ignore) + + // Base class methods + [] + member this.ObjectToString() = + + // Currently the CLR 4.0 and CLR 2.0 behavior of BigInt.ToString is different, causing this test to fail. + + Assert.AreEqual(g_positive1.ToString(), + "1000000000000000000000000000000000018") + Assert.AreEqual(g_zero.ToString(),"0") + Assert.AreEqual(g_normal.ToString(),"88") + + + [] + member this.System_Object_GetHashCode() = + Assert.AreEqual(g_negative1.GetHashCode(),1210897093) + Assert.AreEqual(g_normal.GetHashCode(),89) + Assert.AreEqual(g_zero.GetHashCode(),1) + () + + // Static methods + [] + member this.Abs() = + Assert.AreEqual(bignum.Abs(g_negative1), g_positive1) + Assert.AreEqual(bignum.Abs(g_negative2), g_positive2) + Assert.AreEqual(bignum.Abs(g_positive1), g_positive1) + Assert.AreEqual(bignum.Abs(g_normal), g_normal) + Assert.AreEqual(bignum.Abs(g_zero), g_zero) + () + + [] + member this.FromBigInt() = + Assert.AreEqual(bignum.FromBigInt(g_bigintpositive), + g_positive1) + Assert.AreEqual(bignum.FromBigInt(g_bigintnegative), + g_negative1) + Assert.AreEqual(bignum.FromBigInt(0I),g_zero) + Assert.AreEqual(bignum.FromBigInt(88I),g_normal) + () + + [] + member this.FromInt() = + Assert.AreEqual(bignum.FromInt(2147483647), 2147483647N) + Assert.AreEqual(bignum.FromInt(-2147483648), -2147483648N) + Assert.AreEqual(bignum.FromInt(0), 0N) + Assert.AreEqual(bignum.FromInt(88), 88N) + () + + [] + member this.One() = + Assert.AreEqual(bignum.One,1N) + () + + [] + member this.Parse() = + Assert.AreEqual(bignum.Parse("100"), 100N) + Assert.AreEqual(bignum.Parse("-100"), -100N) + Assert.AreEqual(bignum.Parse("0"), g_zero) + Assert.AreEqual(bignum.Parse("88"), g_normal) + () + + [] + member this.PowN() = + Assert.AreEqual(bignum.PowN(100N,2), 10000N) + Assert.AreEqual(bignum.PowN(-3N,3), -27N) + Assert.AreEqual(bignum.PowN(g_zero,2147483647), 0N) + Assert.AreEqual(bignum.PowN(g_normal,0), 1N) + () + + + [] + member this.Sign() = + Assert.AreEqual(g_positive1.Sign, 1) + Assert.AreEqual(g_negative1.Sign, -1) + Assert.AreEqual(g_zero.Sign, 0) + Assert.AreEqual(g_normal.Sign, 1) + () + + + + [] + member this.ToBigInt() = + Assert.AreEqual(bignum.ToBigInt(g_positive1), g_bigintpositive) + Assert.AreEqual(bignum.ToBigInt(g_negative1), g_bigintnegative) + Assert.AreEqual(bignum.ToBigInt(g_zero), 0I) + Assert.AreEqual(bignum.ToBigInt(g_normal), 88I) + () + + + + [] + member this.ToDouble() = + Assert.AreEqual(bignum.ToDouble(179769N*1000000000000000N), 1.79769E+20) + Assert.AreEqual(bignum.ToDouble(-179769N*1000000000000000N), -1.79769E+20) + Assert.AreEqual(bignum.ToDouble(0N),0.0) + Assert.AreEqual(bignum.ToDouble(88N),88.0) + Assert.AreEqual(double(179769N*1000000000000000N), 1.79769E+20) + Assert.AreEqual(double(-179769N*1000000000000000N), -1.79769E+20) + Assert.AreEqual(double(0N),0.0) + Assert.AreEqual(double(88N),88.0) + () + + + [] + member this.ToInt32() = + Assert.AreEqual(bignum.ToInt32(2147483647N), 2147483647) + Assert.AreEqual(bignum.ToInt32(-2147483648N), -2147483648) + Assert.AreEqual(bignum.ToInt32(0N), 0) + Assert.AreEqual(bignum.ToInt32(88N), 88) + Assert.AreEqual(int32(2147483647N), 2147483647) + Assert.AreEqual(int32(-2147483648N), -2147483648) + Assert.AreEqual(int32(0N), 0) + Assert.AreEqual(int32(88N), 88) + + + + [] + member this.Zero() = + Assert.AreEqual(bignum.Zero,0N) + () + + // operator methods + [] + member this.test_op_Addition() = + + Assert.AreEqual(100N + 200N, 300N) + Assert.AreEqual((-100N) + (-200N), -300N) + Assert.AreEqual(g_positive1 + g_negative1, 0N) + Assert.AreEqual(g_zero + g_zero,0N) + Assert.AreEqual(g_normal + g_normal, 176N) + Assert.AreEqual(g_normal + g_normal, 176N) + () + + + + [] + member this.test_op_Division() = + Assert.AreEqual(g_positive1 / g_positive1, 1N) + Assert.AreEqual(-100N / 2N, -50N) + Assert.AreEqual(g_zero / g_positive1, 0N) + () + + [] + member this.test_op_Equality() = + + Assert.IsTrue((g_positive1 = g_positive1)) + Assert.IsTrue((g_negative1 = g_negative1)) + Assert.IsTrue((g_zero = g_zero)) + Assert.IsTrue((g_normal = g_normal)) + () + + [] + member this.test_op_GreaterThan() = + Assert.AreEqual((g_positive1 > g_positive2), true) + Assert.AreEqual((g_negative1 > g_negative2), false) + Assert.AreEqual((g_zero > g_zero), false) + Assert.AreEqual((g_normal > g_normal), false) + + + () + [] + member this.test_op_GreaterThanOrEqual() = + Assert.AreEqual((g_positive1 >= g_positive2), true) + Assert.AreEqual((g_positive2 >= g_positive1), false) + Assert.AreEqual((g_negative1 >= g_negative1), true) + Assert.AreEqual((0N >= g_zero), true) + + () + [] + member this.test_op_LessThan() = + Assert.AreEqual((g_positive1 < g_positive2), false) + Assert.AreEqual((g_negative1 < g_negative3), false) + Assert.AreEqual((0N < g_zero), false) + + () + [] + member this.test_op_LessThanOrEqual() = + Assert.AreEqual((g_positive1 <= g_positive2), false) + Assert.AreEqual((g_positive2 <= g_positive1), true) + Assert.AreEqual((g_negative1 <= g_negative1), true) + Assert.AreEqual((0N <= g_zero), true) + + () + + [] + member this.test_op_Multiply() = + Assert.AreEqual(3N * 5N, 15N) + Assert.AreEqual((-3N) * (-5N), 15N) + Assert.AreEqual((-3N) * 5N, -15N) + Assert.AreEqual(0N * 5N, 0N) + + () + + [] + member this.test_op_Range() = + let resultPos = [0N .. 2N] + let seqPos = [0N; 1N; 2N] + verifySeqsEqual resultPos seqPos + + let resultNeg = [-2N .. 0N] + let seqNeg = [-2N; -1N; 0N] + verifySeqsEqual resultNeg seqNeg + + let resultSmall = [0N ..5N] + let seqSmall = [0N; 1N; 2N; 3N; 4N; 5N] + verifySeqsEqual resultSmall seqSmall + + () + + + [] + member this.test_op_RangeStep() = + let resultPos = [0N .. 3N .. 6N] + let seqPos = [0N; 3N; 6N] + verifySeqsEqual resultPos seqPos + + let resultNeg = [-6N .. 3N .. 0N] + let seqNeg = [-6N; -3N; 0N] + verifySeqsEqual resultNeg seqNeg + + let resultSmall = [0N .. 3N .. 9N] + let seqSmall = [0N; 3N; 6N; 9N] + verifySeqsEqual resultSmall seqSmall + + () + + [] + member this.test_op_Subtraction() = + Assert.AreEqual(g_positive1 - g_positive2,18N) + Assert.AreEqual(g_negative1 - g_negative3,18N) + Assert.AreEqual(0N-g_positive1, g_negative1) + () + + [] + member this.test_op_UnaryNegation() = + Assert.AreEqual(-g_positive1, g_negative1) + Assert.AreEqual(-g_negative1, g_positive1) + Assert.AreEqual(-0N,0N) + + () + + [] + member this.test_op_UnaryPlus() = + Assert.AreEqual(+g_positive1,g_positive1) + Assert.AreEqual(+g_negative1,g_negative1) + Assert.AreEqual(+0N, 0N) + + () + + // instance methods + [] + member this.Denominator() = + Assert.AreEqual(g_positive1.Denominator, 1I) + Assert.AreEqual(g_negative1.Denominator, 1I) + Assert.AreEqual(0N.Denominator, 1I) + + () + + [] + member this.IsNegative() = + Assert.IsFalse(g_positive1.IsNegative) + Assert.IsTrue(g_negative1.IsNegative) + + Assert.IsFalse( 0N.IsNegative) + Assert.IsFalse(-0N.IsNegative) + + () + + + [] + member this.IsPositive() = + + Assert.IsTrue(g_positive1.IsPositive) + Assert.IsFalse(g_negative1.IsPositive) + + Assert.IsFalse( 0N.IsPositive) + Assert.IsFalse(-0N.IsPositive) + + () + + [] + member this.Numerator() = + Assert.AreEqual(g_positive1.Numerator, g_bigintpositive) + Assert.AreEqual(g_negative1.Numerator, g_bigintnegative) + Assert.AreEqual(0N.Numerator, 0I) + + () + + + + + diff --git a/src/FSharpUnitTests/DistributionTest.fs b/src/FSharpUnitTests/DistributionTest.fs new file mode 100644 index 00000000..b511f18d --- /dev/null +++ b/src/FSharpUnitTests/DistributionTest.fs @@ -0,0 +1,79 @@ +module MathNet.Numerics.Tests.DistributionTest + +open MathNet.Numerics +open MathNet.Numerics.Distribution +open NUnit.Framework +open FsUnit + +[] +let ``When creating a empty distribution, then the probability should be 1``() = + let actual = distribution { return () } + probability actual |> should equal (1N/1N) + +let sumOfTwoFairDices = distribution { + let! d1 = fairDice 6 + let! d2 = fairDice 6 + return d1 + d2 } + +[] +let ``When creating two fair dices, then P(Sum of 2 dices = 7) should be 1/6``() = + sumOfTwoFairDices + |> filter ((=) 7) + |> probability + |> should equal (1N/6N) + +let fairCoinAndDice = distribution { + let! d = fairDice 6 + let! c = fairCoin + return d,c } + +[] +let ``When creating a fair coin and a fair dice, then P(Heads) should be 1/2``() = + fairCoinAndDice + |> filter (fun (_,c) -> c = Heads) + |> probability + |> should equal (1N/2N) + +[] +let ``When creating a fair coin and a fair dice, then P(Heads and dice > 3) should be 1/4``() = + fairCoinAndDice + |> filter (fun (d,c) -> c = Heads && d > 3) + |> probability + |> should equal (1N/4N) + +// MontyHall Problem +// See Martin Erwig and Steve Kollmansberger's paper +// "Functional Pearls: Probabilistic functional programming in Haskell" + +type Outcome = +| Car +| Goat + +let firstChoice = toUniformDistribution [Car; Goat; Goat] + +let switch firstCoice = + match firstCoice with + | Car -> + // If you had the car and you switch ==> you lose since there are only goats left + certainly Goat + | Goat -> + // If you had the goat, the host has to take out another goat ==> you win + certainly Car + +[] +let ``When making the first choice in a MontyHall situation, the chances to win should be 1/3``() = + firstChoice + |> filter ((=) Car) + |> probability + |> should equal (1N/3N) + +let montyHallWithSwitch = distribution { + let! firstDoor = firstChoice + return! switch firstDoor } + +[] +let ``When switching in a MontyHall situation, the chances to win should be 2/3``() = + montyHallWithSwitch + |> filter ((=) Car) + |> probability + |> should equal (2N/3N) \ No newline at end of file diff --git a/src/FSharpUnitTests/FSharpUnitTests.fsproj b/src/FSharpUnitTests/FSharpUnitTests.fsproj index d53cf840..4eefb688 100644 --- a/src/FSharpUnitTests/FSharpUnitTests.fsproj +++ b/src/FSharpUnitTests/FSharpUnitTests.fsproj @@ -46,14 +46,23 @@ + + + + Always + + + ..\..\packages\NUnit.2.6.2\lib\nunit.framework.dll + True + 3.5 diff --git a/src/FSharpUnitTests/PokerDistributionTest.fs b/src/FSharpUnitTests/PokerDistributionTest.fs new file mode 100644 index 00000000..b2702a71 --- /dev/null +++ b/src/FSharpUnitTests/PokerDistributionTest.fs @@ -0,0 +1,96 @@ +module MathNet.Numerics.Tests.PokerDistributionTest + +open MathNet.Numerics +open MathNet.Numerics.Distribution +open NUnit.Framework +open FsUnit + +type Rank = int +type Suit = | Spades | Hearts | Diamonds | Clubs +type Card = Rank * Suit + +let value = fst +let suit = snd + +let A,K,Q,J,T = 14,13,12,11,10 +let allRanksInSuit suit = [2..A] |> List.map (fun rank -> rank,suit) +let completeDeck = + [Spades; Hearts ; Diamonds; Clubs] + |> List.map allRanksInSuit + |> List.concat + +let isPair c1 c2 = value c1 = value c2 +let isSuited c1 c2 = suit c1 = suit c2 +let isConnected c1 c2 = + let v1,v2 = value c1,value c2 + (v1 - v2 |> abs |> (=) 1) || + (v1 = A && v2 = 2) || + (v1 = 2 && v2 = A) + +[] +let ``When drawing from a full deck, then the probability for an Ace should equal 4/52``() = + completeDeck + |> selectOne |> map fst + |> filter (fun card -> value card = A) + |> probability + |> should equal (4N/52N) + +[] +let ``When drawing from a full deck, then the probability should equal 1/52``() = + completeDeck + |> selectOne |> map fst + |> filter ((=) (A,Spades)) + |> probability + |> should equal (1N/52N) + +[] +let ``When drawing from a full deck, then the probability for the Ace of Clubs and Ace of Spaces (in order) should equal 1/52 * 1/51``() = + completeDeck + |> select 2 + |> filter ((=) [A,Clubs; A,Spades]) + |> probability + |> should equal (1N/52N * 1N/51N) + +[] +let ``When drawing from a full deck, then the probability for the Ace of Clubs and Ace of Spaces (in any order) should equal (1/52 * 1/51) * 2``() = + completeDeck + |> select 2 + |> filterInAnyOrder [A,Clubs; A,Spades] + |> probability + |> should equal ((1N/52N * 1N/51N) * 2N) + +[] +let ``When drawing the Ace of Spades and the Ace of Clubs, then the probability for drawing another Ace should equal 2/50``() = + completeDeck + |> remove [A,Clubs; A,Spades] + |> toUniformDistribution + |> filter (fun card -> value card = A) + |> probability + |> should equal (2N/50N) + + +[] +let ``When drawing from the full deck, then the probability for drawing a Pair preflop should equal 1/17``() = + completeDeck + |> select 2 + |> filter (fun (c1::c2::_) -> isPair c1 c2) + |> probability + |> should equal (1N/17N) + +[] +let ``When drawing from the full deck, then the probability for drawing Suited Connectors should equal 1/25``() = + completeDeck + |> select 2 + |> filter (fun (c1::c2::_) -> isSuited c1 c2 && isConnected c1 c2) + |> probability + |> should equal (2N/51N) + +[] +let ``When holding 3 Spades after the flop, than the probability for drawing a flush should equal 10/47*9/46``() = + completeDeck + |> remove [A,Clubs; A,Spades] // preflop + |> remove [2,Clubs; 3,Spades; 7,Spades] // flop + |> select 2 + |> filter (fun (c1::c2::_) -> suit c1 = Spades && suit c2 = Spades) + |> probability + |> should equal (10N/47N*9N/46N) \ No newline at end of file diff --git a/src/FSharpUnitTests/Utilities.fs b/src/FSharpUnitTests/Utilities.fs new file mode 100644 index 00000000..77975580 --- /dev/null +++ b/src/FSharpUnitTests/Utilities.fs @@ -0,0 +1,128 @@ +// First version copied from the F# Power Pack +// https://raw.github.com/fsharp/powerpack/master/src/FSharp.PowerPack.Unittests/Utilities.fs + +namespace MathNet.Numerics.Tests +open NUnit.Framework +open System +open System.Collections.Generic + +[] +module Utilities = + let test msg b = Assert.IsTrue(b, "MiniTest '" + msg + "'") + let logMessage msg = + System.Console.WriteLine("LOG:" + msg) +// System.Diagnostics.Trace.WriteLine("LOG:" + msg) + let check msg v1 v2 = test msg (v1 = v2) + let reportFailure msg = Assert.Fail msg + let numActiveEnumerators = ref 0 + let throws f = try f() |> ignore; false with e -> true + + let countEnumeratorsAndCheckedDisposedAtMostOnceAtEnd (seq: seq<'a>) = + let enumerator() = + numActiveEnumerators := !numActiveEnumerators + 1; + let disposed = ref false in + let endReached = ref false in + let ie = seq.GetEnumerator() in + { new System.Collections.Generic.IEnumerator<'a> with + member x.Current = + test "rvlrve0" (not !endReached); + test "rvlrve1" (not !disposed); + ie.Current + member x.Dispose() = + test "rvlrve2" !endReached; + test "rvlrve4" (not !disposed); + numActiveEnumerators := !numActiveEnumerators - 1; + disposed := true; + ie.Dispose() + interface System.Collections.IEnumerator with + member x.MoveNext() = + test "rvlrve0" (not !endReached); + test "rvlrve3" (not !disposed); + endReached := not (ie.MoveNext()); + not !endReached + member x.Current = + test "qrvlrve0" (not !endReached); + test "qrvlrve1" (not !disposed); + box ie.Current + member x.Reset() = + ie.Reset() + } in + + { new seq<'a> with + member x.GetEnumerator() = enumerator() + interface System.Collections.IEnumerable with + member x.GetEnumerator() = (enumerator() :> _) } + + let countEnumeratorsAndCheckedDisposedAtMostOnce (seq: seq<'a>) = + let enumerator() = + let disposed = ref false in + let endReached = ref false in + let ie = seq.GetEnumerator() in + numActiveEnumerators := !numActiveEnumerators + 1; + { new System.Collections.Generic.IEnumerator<'a> with + member x.Current = + test "qrvlrve0" (not !endReached); + test "qrvlrve1" (not !disposed); + ie.Current + member x.Dispose() = + test "qrvlrve4" (not !disposed); + numActiveEnumerators := !numActiveEnumerators - 1; + disposed := true; + ie.Dispose() + interface System.Collections.IEnumerator with + member x.MoveNext() = + test "qrvlrve0" (not !endReached); + test "qrvlrve3" (not !disposed); + endReached := not (ie.MoveNext()); + not !endReached + member x.Current = + test "qrvlrve0" (not !endReached); + test "qrvlrve1" (not !disposed); + box ie.Current + member x.Reset() = + ie.Reset() + } in + + { new seq<'a> with + member x.GetEnumerator() = enumerator() + interface System.Collections.IEnumerable with + member x.GetEnumerator() = (enumerator() :> _) } + + // Verifies two sequences are equal (same length, equiv elements) + let verifySeqsEqual seq1 seq2 = + if Seq.length seq1 <> Seq.length seq2 then Assert.Fail() + + let zippedElements = Seq.zip seq1 seq2 + if zippedElements |> Seq.forall (fun (a, b) -> a = b) + then () + else Assert.Fail() + + /// Check that the lamda throws an exception of the given type. Otherwise + /// calls Assert.Fail() + let private checkThrowsExn<'a when 'a :> exn> (f : unit -> unit) = + let funcThrowsAsExpected = + try + let _ = f () + false // Did not throw! + with + | :? 'a + -> true // Thew null ref, OK + | _ -> false // Did now throw a null ref exception! + if funcThrowsAsExpected + then () + else Assert.Fail() + + // Illegitimate exceptions. Once we've scrubbed the library, we should add an + // attribute to flag these exception's usage as a bug. + let checkThrowsNullRefException f = checkThrowsExn f + let checkThrowsIndexOutRangException f = checkThrowsExn f + + // Legit exceptions + let checkThrowsNotSupportedException f = checkThrowsExn f + let checkThrowsArgumentException f = checkThrowsExn f + let checkThrowsArgumentNullException f = checkThrowsExn f + let checkThrowsKeyNotFoundException f = checkThrowsExn f + let checkThrowsDivideByZeroException f = checkThrowsExn f + let checkThrowsInvalidOperationExn f = checkThrowsExn f + + \ No newline at end of file diff --git a/src/FSharpUnitTests/packages.config b/src/FSharpUnitTests/packages.config new file mode 100644 index 00000000..5c3ca54d --- /dev/null +++ b/src/FSharpUnitTests/packages.config @@ -0,0 +1,4 @@ + + + + \ No newline at end of file