diff --git a/src/FSharpPortable/FSharpPortable.fsproj b/src/FSharpPortable/FSharpPortable.fsproj index 67038619..420908c3 100644 --- a/src/FSharpPortable/FSharpPortable.fsproj +++ b/src/FSharpPortable/FSharpPortable.fsproj @@ -19,7 +19,7 @@ false false ..\..\out\debug\Portable\ - TRACE;DEBUG;PORTABLE + TRACE;DEBUG;PORTABLE;FX_NO_BIGINT 3 @@ -29,7 +29,7 @@ true true ..\..\out\lib\Portable\ - TRACE;PORTABLE + TRACE;PORTABLE;FX_NO_BIGINT 3 ..\..\out\lib\Portable\MathNet.Numerics.FSharp.XML @@ -62,6 +62,25 @@ LinearAlgebra.Double.fs + + + + + + complex.fsi + + + complex.fs + + + q.fsi + + + q.fs + + + DistributionMonad.fs + diff --git a/src/FSharpPortable/n.fs b/src/FSharpPortable/n.fs new file mode 100644 index 00000000..a5c97799 --- /dev/null +++ b/src/FSharpPortable/n.fs @@ -0,0 +1,1611 @@ +// First version copied from the F# compiler sources +// https://raw.github.com/fsharp/fsharp/master/src/fsharp/FSharp.Core/math/n.fs + +//---------------------------------------------------------------------------- +// Copyright (c) 2002-2012 Microsoft Corporation. +// +// This source code is subject to terms and conditions of the Apache License, Version 2.0. A +// copy of the license can be found in the License.html file at the root of this distribution. +// By using this source code in any fashion, you are agreeing to be bound +// by the terms of the Apache License, Version 2.0. +// +// You must not remove this notice, or any other, from this software. +//---------------------------------------------------------------------------- + +namespace Microsoft.FSharp.Math + +#if FX_NO_BIGINT +open System +open System.Diagnostics.CodeAnalysis +open Microsoft.FSharp.Core +open Microsoft.FSharp.Core.LanguagePrimitives.IntrinsicOperators +open Microsoft.FSharp.Core.Operators +open Microsoft.FSharp.Collections +open Microsoft.FSharp.Primitives.Basics + +type ints = int array + +[] +type internal BigNat = + + // Have n = sum (from i=0 to bound) a.[i] * baseN ^ i + // Have 0 <= a.[i] < baseN. + //------ + // Invariant: bound is least such, i.e. bound=0 or (a.[bound-1] is highest coeff). + // Zero is {bound=0,a=...}. + // Naturals are a normal form, + // but not structurally so, + // since arrays may have non-contributing cells at a.[bound] and beyond. + // + { mutable bound : int; // non-zero coeff must be 0...(bound-1) + digits : ints // must have at least elts 0...(bound-1), + // maybe more (which should be zero!). + // Actually, the "zero" condition may be relaxed. + // + } + + +module internal BigNatModule = + + //------------------------------------------------------------------------- + // misc + //----------------------------------------------------------------------- + + #if SELFTEST + let check b = if not b then failwith "assertion failwith" + #endif + + module FFT = + let rec pow32 x n = + if n=0 then 1 + elif n % 2 = 0 then pow32 (x*x) (n / 2) + else x* pow32 (x*x) (n / 2) + + let leastBounding2Power b = + let rec findBounding2Power b tp i = if b<=tp then tp,i else findBounding2Power b (tp*2) (i+1) in + findBounding2Power b 1 0 + + //------------------------------------------------------------------------- + // p = 2^k.m + 1 prime and w primitive 2^k root of 1 mod p + //----------------------------------------------------------------------- + + // Given p = 2^k.m + 1 prime and w a primitive 2^k root of unity (mod p). + // Required to define arithmetic ops for Fp = field modulo p. + // The following are possible choices for p. + + // p, k, m, g, w + // let p,k,m,g,w = 97L, 4, 6, 5, 8 // p is 7 bit + // let p,k,m,g,w = 769L, 8, 3, 7, 7 // p is 10 bit + // let p,k,m,g,w = 7681L, 8, 30, 13, 198 // p is 13 bit + // let p,k,m,g,w = 12289L, 10, 12, 11, 49 // p is 14 bit + // let p,k,m,g,w = 167772161L, 25, 5, 557092, 39162105 // p is 28 bit + // let p,k,m,g,w = 469762049L, 26, 7, 1226571, 288772249 // p is 29 bit + + + let p,k,m,g,w = 2013265921L, 27, 15, 31, 440564289 // p is 31 bit + let primeP = p + + let maxBitsInsideFp = 30 + + + //------------------------------------------------------------------------- + // Fp = finite field mod p - rep is uint32 + //----------------------------------------------------------------------- + + + type fp = uint32 + // operations in Fp (finite field size p) + module Fp = + //module I = UInt32 + let p = 2013265921ul : fp + let p64 = 2013265921UL : uint64 + let toInt (x:fp) : int = int32 x + let ofInt32 (x:int) : fp = uint32 x + + let mzero : fp = 0ul + let mone : fp = 1ul + let mtwo : fp = 2ul + let inline madd (x:fp) (y:fp) : fp = (x + y) % p + let inline msub (x:fp) (y:fp) : fp = (x + p - y) % p + let inline mmul (x:fp) (y:fp) : fp = uint32 ((uint64 x * uint64 y) % p64) + + let rec mpow x n = + if n=0 then mone + elif n % 2=0 then mpow (mmul x x) (n / 2) + else mmul x (mpow (mmul x x) (n / 2)) + + let rec mpowL x n = + if n = 0L then mone + elif n % 2L = 0L then mpowL (mmul x x) (n / 2L) + else mmul x (mpowL (mmul x x) (n / 2L)) + + // Have the w is primitive 2^kth root of 1 in Zp + let m2PowNthRoot n = + // Find x s.t. x is (2^n)th root of unity. + // + // pow w (pow 2 k) = 1 primitively. + // = pow w (pow 2 ((k-n)+n)) + // = pow w (pow 2 (k-n) * pow 2 n) + // = pow (pow w (pow 2 (k-n))) (pow 2 n) + // + // Take wn = pow (pow w (pow 2 (k-n))) + + mpow (uint32 w) (pow32 2 (k-n)) + + let minv x = mpowL x (primeP - 2L) + + + //------------------------------------------------------------------------- + // FFT - in place low garbage + //----------------------------------------------------------------------- + + open Fp + let rec computeFFT lambda mu n w (u: _[]) (res: _[]) offset = + // Given n a 2-power, + // w an nth root of 1 in Fp, and + // lambda, mu and u(x) defining + // poly(lambda,mu,x) = sum(i pow32 2 i) + + let computeFftPaddedPolynomialProduct bigK k u v = + // REQUIRES: bigK = 2^k + // REQUIRES: Array lengths of u and v = bigK. + // REQUIRES: degree(uv) <= bigK-1 + // --- + // Given u,v polynomials. + // Computes the product polynomial by FFT. + // For correctness, + // require the result coeff to be in range [0,p-1], for p defining Fp above. + + #if SELFTEST + check ( k <= maxTwoPower ); + check ( bigK = twoPowerTable.[k] ); + check ( u.Length = bigK ); + check ( v.Length = bigK ); + #endif + // Find 2^k primitive root of 1 + let w = m2PowNthRoot k + // FFT + let n = bigK + let uT = computFftInPlace n w u + let vT = computFftInPlace n w v + // Evaluate + let rT = Array.init n (fun i -> mmul uT.[i] vT.[i]) + // INV FFT + let r = computeInverseFftInPlace n w rT + r + + let padTo n (u: _ array) = + let uBound = u.Length + Array.init n (fun i -> if i mmul uT.[i] vT.[i]) + // INV FFT + let r = computeInverseFftInPlace n w rT + Array.map Fp.toInt r + + + //------------------------------------------------------------------------- + // fp exports + //----------------------------------------------------------------------- + + open Fp + let mzero = mzero + let mone = mone + let maxFp = msub Fp.p mone + + //------------------------------------------------------------------------- + // FFT - reference implementation + //----------------------------------------------------------------------- + + #if SELFTEST + open Fp + let rec computeFftReference n w u = + // Given n a 2-power, + // w an nth root of 1 in Fp, and + // u(x) = sum(i u.[2*i]) + let uodd = Array.init (n/2) (fun i -> u.[2*i+1]) + let uevenFT = computeFftReference (n/2) (mmul w w) ueven + let uoddFT = computeFftReference (n/2) (mmul w w) uodd + Array.init n + (fun j -> + if j < n/2 then + madd + (uevenFT.[j]) + (mmul + (mpow w j) + (uoddFT.[j])) + else + let j = j - (n/2) + msub + (uevenFT.[j]) + (mmul + (mpow w j) + (uoddFT.[j]))) + #endif + + open FFT + + type n = BigNat + + let bound (n: n) = n.bound + let setBound (n: n) (v:int32) = n.bound <- v + let coeff (n:n) i = n.digits.[i] + let coeff64 (n:n) i = int64 (coeff n i) + let setCoeff (n:n) i v = n.digits.[i] <- v + + let rec pow64 x n = + if n=0 then 1L + elif n % 2 = 0 then pow64 (x * x) (n / 2) + else x * (pow64 (x * x) (n / 2)) + + let rec pow32 x n = + if n=0 then 1 + elif n % 2 = 0 then pow32 (x*x) (n / 2) + else x* pow32 (x*x) (n / 2) + + let hash(n) = + let mutable res = 0 + for i = 0 to n.bound - 1 do // could stop soon, it's "hash" + res <- n.digits.[i] + (res <<< 3) + done; + res + + //---------------------------------------------------------------------------- + // misc + //-------------------------------------------------------------------------- + +#if CHECKED + let check b str = if not b then failwith ("check failed: " + str) +#endif + let maxInt a b = if a int32 + let inline div64base (x:int64) = int64 (uint64 x >>> baseBits) + + let divbase x = int32 (uint32 x >>> baseBits) + let modbase x = (x &&& baseMask) + + let inline index z i = if i < z.bound then z.digits.[i] else 0 + + let createN b = { bound = b; + digits = Array.zeroCreate b } + let copyN x = { bound = x.bound; + digits = Array.copy x.digits } // could copy just enough... + + let normN n = + // normalises bound + let rec findLeastBound (na:ints) i = if i = -1 || na.[i]<>0 then i+1 else findLeastBound na (i-1) + let bound = findLeastBound n.digits (n.bound-1) + n.bound <- bound; + n + + let boundInt = 2 // int will fit with bound=2 + let boundInt64 = 3 // int64 will fit with bound=3 + let boundBase = 1 // base will fit with bound=1 - obviously! + +//---------------------------------------------------------------------------- +// base, coefficients, poly +//-------------------------------------------------------------------------- + + let embed x = + let x = if x<0 then 0 else x // no -ve naturals + if x < baseN then + let r = createN 1 + r.digits.[0] <- x; + normN r + else + let r = createN boundInt + for i = 0 to boundInt - 1 do + r.digits.[i] <- (x / pow32 baseN i) % baseN + done; + normN r + + let embed64 x = + let x = if x<0L then 0L else x // no -ve naturals + let r = createN boundInt64 + for i = 0 to boundInt64-1 do + r.digits.[i] <- int32 ( (x / pow64 baseNi64 i) % baseNi64) + done; + normN r + + let eval n = + if n.bound = 1 + then n.digits.[0] + else + let mutable acc = 0 + for i = n.bound-1 downto 0 do + acc <- n.digits.[i] + baseN * acc + done; + acc + + let eval64 n = + if n.bound = 1 + then int64 n.digits.[0] + else + let mutable acc = 0L + for i = n.bound-1 downto 0 do + acc <- int64 (n.digits.[i]) + baseNi64 * acc + done; + acc + + let one = embed 1 + let zero = embed 0 + + let restrictTo d n = + { bound = minInt d n.bound; digits = n.digits} + + let shiftUp d n = + let m = createN (n.bound+d) + for i = 0 to n.bound-1 do + m.digits.[i+d] <- n.digits.[i] + done; + m + + let shiftDown d n = + if n.bound-d<=0 then + zero + else + let m = createN (n.bound-d) + for i = 0 to m.bound-1 do + m.digits.[i] <- n.digits.[i+d] + done; + m + + let degree n = n.bound-1 + + +//---------------------------------------------------------------------------- +// add, sub +//-------------------------------------------------------------------------- + + // addition + let rec addP i n c p q r = // p+q + c + if i0 then + r.digits.[i] <- modbase x; + let c = divbase x + // if p (or q) exhausted and c zero could switch to copying mode + subP (i+1) n c p q r + else + let x = x + baseN // add baseN + r.digits.[i] <- modbase x; + let c = divbase x - 1 // sub baseN + // if p (or q) exhausted and c zero could switch to copying mode + subP (i+1) n c p q r + else + let underflow = c<>0 + underflow + + let sub p q = + // NOTE: x-y=0 when x<=y, it is natural subtraction + let rbound = maxInt p.bound q.bound + let r = createN rbound + let carry = 0 + let underflow = subP 0 rbound carry p q r + if underflow then + embed 0 + else + normN r + + +//---------------------------------------------------------------------------- +// isZero, equal, ordering, sign, min, max +//-------------------------------------------------------------------------- + + let isZero p = p.bound=0 + let IsZero p = isZero p + let isOne p = p.bound=1 && p.digits.[0] = 1 + + let equal p q = + (p.bound = q.bound) && + (let rec check (pa:ints) (qa:ints) i = + // HAVE: pa.[j] = qa.[j] for i < j < p.bound + (i = -1) || (pa.[i]=qa.[i] && check pa qa (i-1)) + + check p.digits q.digits (p.bound-1)) + + let shiftCompare p pn q qn = + if p.bound + pn < q.bound + qn then -1 + elif p.bound + pn > q.bound + pn then 1 + else + let rec check (pa:ints) (qa:ints) i = + // HAVE: pa.[j-pn] = qa.[j-qn] for i < j < p.bound + // Looking for most significant differing coeffs to determine ordering + if i = -1 then + 0 + else + let pai = if i < pn then 0 else pa.[i-pn] + let qai = if i < qn then 0 else qa.[i-qn] + if pai = qai then check pa qa (i-1) + elif pai < qai then -1 + else 1 + + check p.digits q.digits (p.bound + pn - 1) + + let compare p q = + if p.bound < q.bound then -1 + elif p.bound > q.bound then 1 + else + let rec check (pa:ints) (qa:ints) i = + // HAVE: pa.[j] = qa.[j] for i < j < p.bound + // Looking for most significant differing coeffs to determine ordering + if i = -1 then 0 + elif pa.[i]=qa.[i] then check pa qa (i-1) + elif pa.[i]] + let lt p q = compare p q = -1 + [] + let gt p q = compare p q = 1 + [] + let lte p q = compare p q <> 1 + [] + let gte p q = compare p q <> -1 + + [] + let min a b = if lt a b then a else b + [] + let max a b = if lt a b then b else a + + +//---------------------------------------------------------------------------- +// scale +//-------------------------------------------------------------------------- + + // REQUIRE: baseN + baseN.2^32 < Int64.maxInt + let rec contributeArr (a:ints) i (c:int64) = + // Given c and require c < baseN.2^32 + // Compute: r <- r + c . B^i + // via r.digits.[i] <- r.digits.[i] + c and normalised + let x = int64 a.[i] + c + // HAVE: x < baseN + baseN.2^32 + let c = div64base x + let x = mod64base x + // HAVE: c < 1 + 2^32 < baseN.2^32, recursive call ok + // HAVE: x < baseN + a.[i] <- x; // store residue x + if c>0L then + contributeArr a (i+1) c // contribute carry next position + + let inline contribute r i c = contributeArr r.digits i c + + // REQUIRE: maxInt < 2^32 + [] + let rec scale (k:int) (p:n) = + // Given k and p and require k < 2^32 + // Computes "scalar" product k.p + // + let rbound = p.bound + boundInt + let r = createN rbound + let k = int64 k + for i = 0 to p.bound-1 do + let kpi = k * int64 p.digits.[i] + // HAVE: kpi < 2^32 * baseN which meets "contribute" requirement + contribute r i kpi + done; + normN r + + +//---------------------------------------------------------------------------- +// mulSchoolBook +//-------------------------------------------------------------------------- + + // multiplication: naively O(n^2) +(* + let mulSchoolBook' p q = + let rbound = p.bound + q.bound + boundBase*2 + let r = createN rbound + let pa = p.digits + let qa = q.digits + for i = 0 to p.bound-1 do + for j = 0 to q.bound-1 do + contribute r (i+j) (int64 pa.[i] * int64 qa.[j]) + done + done; + normN r +*) + + let mulSchoolBookBothSmall p q = + let r = createN 2 + let rak = int64 p * int64 q + setCoeff r 0 (mod64base rak); + setCoeff r 1 (int32 (div64base rak)) + normN r + + let rec mulSchoolBookCarry r c k = + if ( c > 0L ) then + // ToAdd = c.B^k + let rak = (coeff64 r k) + c + setCoeff r k (mod64base rak); + mulSchoolBookCarry r (div64base rak) (k + 1) + + let mulSchoolBookOneSmall p q = + let bp = bound(p) + let rbound = bp + 1 + let r = createN rbound + let q = int64 q + let mutable c = 0L + for i = 0 to bp-1 do + let rak = c + (coeff64 r i) + (coeff64 p i) * q + setCoeff r i (mod64base rak); + c <- div64base rak; + mulSchoolBookCarry r c bp + normN r + + + // multiplication: naively O(n^2) -- this version - unchecked - is faster + let mulSchoolBookNeitherSmall p q = + let rbound = p.bound + q.bound + let r = createN rbound + let ra = r.digits + let pa = p.digits + let qa = q.digits + // ToAdd p*q + for i = 0 to p.bound-1 do + // ToAdd p.[i] * q * B^i + let pai = int64 pa.[i] + let mutable c = 0L + let mutable k = i // k = i + j + // ToAdd = pi.qj.B^(i+j) for j = 0,j+1... + for j = 0 to q.bound-1 do + // ToAdd = c.B^k + pi.qj.B^(i+j) for j = j,j+1... and k = i+j + let qaj = int64 qa.[j] + let rak = int64 ra.[k] + c + pai * qaj + ra.[k] <- int32 (mod64base rak); + c <- div64base rak; + k <- k + 1; + mulSchoolBookCarry r c k + normN r + + let mulSchoolBook p q = + let pSmall = (bound(p) = 1) + let qSmall = (bound(q) = 1) + if (pSmall && qSmall) then mulSchoolBookBothSmall (coeff p 0) (coeff q 0) + elif pSmall then mulSchoolBookOneSmall q (coeff p 0) + elif qSmall then mulSchoolBookOneSmall p (coeff q 0) + else mulSchoolBookNeitherSmall p q + + +//---------------------------------------------------------------------------- +// quickMulUsingFft +//-------------------------------------------------------------------------- + + // The FFT polynomial multiplier requires the result coeffs fit inside Fp. + // + // OVERVIEW: + // The numbers are recoded as polynomials to be evaluated at (x=2^bigL). + // The polynomials are FFT multiplied, requiring result coeff to fit Fp. + // The result product is recovered by evaluating the poly at (x=2^bigL). + // + // REF: + // QuickMul: Practical FFT-base Integer Multiplication, + // Chee Yap and Chen Yi. + // + // There is choice of how to encode the nats polynomials. + // The choice is the (2^bigL) base to use. + // For bigL=1, the FFT will cater for a product of upto 256M bits. + // Larger bigL have less reach, but compute faster. + // So plan to choose bigL depending on the number of bits product. + // + // DETERMINING THE K,L BOUNDS. + // + // Given representing using K-vectors, K a power of 2, K=2^k, and + // If choosing inputs to have L-bit coefficients. + // + // The result coeff are: + // + // res(i) = sum (j] + type encoding = + { bigL : int; // bits per input coeff + twoToBigL : int; // 2^bigL + k : int; + bigK : int; // bigK = 2^k, number of terms polynomials + bigN : int; // bits result (under-estimate of limit) + split : int; // baseBits / bigL + splits : int array; + } + +#if CHECKED + let _ = check (baseBits=24) "24bit" +#endif + // Requiring baseN mod 2^bigL = 0 gave quick encoding, but... + // also a terrible drop performance when the bigK jumped by more than needed! + // Below, it choose a minimal bigK to hold the product. + + let mkEncoding (bigL,k,bigK,bigN) = +#if CHECKED + check (bigK = pow32 2 k) "bigK"; + check (bigN = bigK * bigL) "bigN"; + check (2 * bigL + k <= maxBitsInsideFp) "constraint"; +#endif + { bigL = bigL; + twoToBigL = pow32 2 bigL; + k = k; + bigK = bigK; + bigN = bigN; + split = baseBits/bigL; // should divide exactly + splits = Array.init (baseBits/bigL) (fun i -> pow32 2 (bigL*i)) + } + + let table = + [| // bigL , k , bigK , bigN // + mkEncoding ( 1 , 28 , 268435456 , 268435456 ) ; + mkEncoding ( 2 , 26 , 67108864 , 134217728 ) ; + mkEncoding ( 3 , 24 , 16777216 , 50331648 ) ; + mkEncoding ( 4 , 22 , 4194304 , 16777216 ) ; + mkEncoding ( 5 , 20 , 1048576 , 5242880 ) ; + mkEncoding ( 6 , 18 , 262144 , 1572864 ) ; + mkEncoding ( 7 , 16 , 65536 , 458752 ) ; + mkEncoding ( 8 , 14 , 16384 , 131072 ) ; + mkEncoding ( 9 , 12 , 4096 , 36864 ) ; + mkEncoding ( 10 , 10 , 1024 , 10240 ) ; + mkEncoding ( 11 , 8 , 256 , 2816 ) ; + mkEncoding ( 12 , 6 , 64 , 768 ) ; + mkEncoding ( 13 , 4 , 16 , 208 ) ; + |] + + let calculateTableTow bigL = + // Given L. + // Have L via "log2 K <= maxBitsInsideFp - 2L". + // Have N via "N = K.L" + // + let k = maxBitsInsideFp - 2*bigL + let bigK = pow64 2L k + let N = bigK * int64 bigL + bigL,k,bigK,N + + let encodingGivenResultBits bitsRes = + // choose maximum bigL s.t. bitsRes < bigN + // EXCEPTION: fails is bitsRes exceeds 2^28 (largest bigN table) + let rec selectFrom i = + if i+1 < table.Length && bitsRes < table.[i+1].bigN then + selectFrom (i+1) + else + table.[i] + + if bitsRes >= table.[0].bigN then + failwith "Product is huge, around 268435456 bits, beyond quickmul" + else + selectFrom 0 + + let bitmask = Array.init baseBits (fun i -> (pow32 2 i - 1)) + let twopowers = Array.init baseBits (fun i -> (pow32 2 i)) + let twopowersI64 = Array.init baseBits (fun i -> (pow64 2L i)) + // bitmask(k) = 2^k - 1 + // twopowers(k) = 2^k // + + let wordBits word = + let rec hi k = + if k=0 then 0 + elif (word &&& twopowers.[k-1]) <> 0 then k + else hi (k-1) + + hi baseBits + + let bits u = + if u.bound=0 then 0 + else degree u * baseBits + wordBits u.digits.[degree u] + + let extractBits n enc bi = + let bj = bi + enc.bigL - 1 // the last bit (inclusive) + let biw = bi / baseBits // first bit is this index pos + let bjw = bj / baseBits // last bit is this index pos + if biw <> bjw then + // two words + let x = index n biw + let y = index n bjw // bjw = biw+1 + let xbit = bi % baseBits // start bit x + let nxbits = baseBits - xbit // number of bitsin x + let x = x >>> xbit // shift down x so bit0 is first + let y = y <<< nxbits // shift up y so it starts where x finished + let x = x ||| y // combine them + let x = x &&& bitmask.[enc.bigL] // mask out (high y bits) to get required bits + x + else + // one word + let x = index n biw + let xbit = bi % baseBits // start bit x + let x = x >>> xbit + let x = x &&& bitmask.[enc.bigL] + x + + let encodePoly enc n = + // Find poly s.t. n = poly evaluated at x=2^bigL + // with 0 <= pi < 2^bigL. + // + let poly = Array.create enc.bigK (Fp.ofInt32 0) + let biMax = n.bound * baseBits + let rec encoder i bi = + // bi = i * bigL + if i=enc.bigK || bi > biMax then + () // done + else + ( let pi = extractBits n enc bi + poly.[i] <- Fp.ofInt32 pi; + let i = i + 1 + let bi = bi + enc.bigL + encoder i bi + ) + + encoder 0 0; + poly + + let decodeResultBits enc (poly : fp array) = + // Decoding evaluates poly(x) (coeff Fp) at X = 2^bigL. + // A bound on number of result bits is "enc.bigN + boundInt", but that takes HUGE STEPS. + // Garbage has a cost, so we minimize it by working out a tight bound. + // + // poly(X) = sum i=0..n coeff_i * X^i where n is highest non-zero coeff. + // <= 2^maxBitsInsideFp * (1 + X + ... X^n) + // <= 2^maxBitsInsideFp * (X^(n+1) - 1) / (X - 1) + // <= 2^maxBitsInsideFp * X^(n+1) / (X - 1) + // <= 2^maxBitsInsideFp * X^(n+1) / (X/2) provided X/2 <= X-1 + // <= 2^maxBitsInsideFp * X^n * 2 + // <= 2^maxBitsInsideFp * (2^bigL)^n * 2 + // <= 2^(maxBitsInsideFp + bigL.n + 1) + // + let mutable n = 0 + for i = 0 to poly.Length-1 do + if poly.[i] <> mzero then n <- i + done; + let rbits = maxBitsInsideFp + enc.bigL * n + 1 + rbits + 1 // +1 since 2^1 requires 2 bits not 1 + + // REQUIRE: bigL <= baseBits + let decodePoly enc poly = + // Find n = poly evaluated at x=2^bigL + // Note, 0 <= pi < maxFp. + // + let rbound = (decodeResultBits enc poly) / baseBits + 1 + let r = createN rbound + let rec evaluate i j d = + // HAVE: bigL.i = j * baseBits + d and d= rbound then +#if CHECKED + check (poly.[i] = mzero) "decodePoly"; +#endif + () + else ( + let x = int64 (Fp.toInt poly.[i]) * twopowersI64.[d] + // HAVE: x < 2^32 . 2^baseBits = 2^32.baseN + contribute r j x + ); + let i = i + 1 + let d = d + enc.bigL + let j,d = if d >= baseBits then j+1 , d-baseBits else j,d + // HAVE: d < baseBits, note: bigL minDigitsKaratsuba then + let k = bmax / 2 + let a0 = restrictTo k p + let a1 = shiftDown k p + let b0 = restrictTo k q + let b1 = shiftDown k q + let q0 = mul a0 b0 + let q1 = mul (add a0 a1) (add b0 b1) + let q2 = mul a1 b1 + let p0 = q0 + let p1 = sub q1 (add q0 q2) + let p2 = q2 + let r = add p0 (shiftUp k (add p1 (shiftUp k p2))) + r + else + mulSchoolBook p q + + let rec mulKaratsuba x y = recMulKaratsuba mulKaratsuba x y + + +//---------------------------------------------------------------------------- +// mul - composite +//-------------------------------------------------------------------------- + + let productDigitsUpperSchoolBook = (64000 / baseBits) + // When is it worth switching away from SchoolBook? + // SchoolBook overhead is low, so although it's O(n^2) it remains competitive. + // + // 28/3/2006: + // The FFT can take over from SchoolBook at around 64000 bits. + // Note, FFT performance is stepwise, according to enc from table. + // The steps are big steps (meaning sudden jumps/drops perf). + // + + let singleDigitForceSchoolBook = (32000 / baseBits) + // If either argument is "small" then stay with SchoolBook. + // + + let productDigitsUpperFft = (table.[0].bigN / baseBits) + // QuickMul is good upto a finite (but huge) limit: + // Limit 268,435,456 bits product. + // + // From the code: + // let bitsRes = bits u + bits v + // fails when bitsRes >= table.[0].bigN + // So, not applicable when: + // P1: table.[0].bigN <= bits(u) + bits(v) + // P2: table.[0].bigN <= .. <= baseBits * (u.bound + v.bound) + // P3: table.[0].bigN <= .. <= baseBits * (u.bound + v.bound) + // P4: table.[0].bigN / baseBits <= u.bound + v.bound + // + + // Summary of mul algorithm choice: + // 0 <= uv_bound < upper_school_book - Schoolbook + // upper_school_book <= uv_bound < upper_fft - QuickMul + // upper_fft <= uv_bound < ... - Karatsuba + // + // NOTE: + // - Karatsuba current implementation has high garbage cost. + // - However, a linear space cost is possible... + // - Meantime, switch to Karatsuba only beyond FFT range. + // + + let rec mul p q = + let pqBound = p.bound + q.bound + if pqBound < productDigitsUpperSchoolBook || + p.bound < singleDigitForceSchoolBook || + q.bound < singleDigitForceSchoolBook + then + // Within school-book initial range: + mulSchoolBook p q + else + if pqBound < productDigitsUpperFft then + // Inside QuickMul FFT range: + quickMulUsingFft p q + else + // Beyond QuickMul FFT range, or maybe between Schoolbook and QuickMul (no!): + // Use karatsuba method, with "mul" as recursive multiplier, + // so will reduce sizes of products on recursive calls, + // and QuickMul will take over if they fall within it's range. + // + recMulKaratsuba mul p q + + +//---------------------------------------------------------------------------- +// division - scaleSubInPlace +//-------------------------------------------------------------------------- + + let scaleSubInPlace x f a n = + // Have x = sumR 0 xd (\i.xi.B^i) where xd = degree x + // a = sumR 0 ad (\i.digitsi.B^i) where ad = degree a + // f < B + // n < xd + // Assumes "f.digits.B^n < x". + // Required to remove f.digits.B^n from x place. + //------ + // Result = x_initial - f.digits.B^n + // = x_initial - f.[sumR 0 ad (\i.digitsi.B^(i+n))] + // State: j = 0 + // z = f * a0 + // Invariant(x,z,j,n): + // P1: x_result = x - z.B^(j+n) - f.[sumR (j+1) ad (\i.digitsi.B^i+n)] + // P2: z < B^2 - 1, and so has form z = zHi.B + zLo for zHi,zLo < B. + // Base: Invariant holds initially. + // Step: (a) Remove zLo from x: + // If zLo <= x_(j+n) then zLo <- 0 + // x_(j+n) <- x_(j+n) - zLo + // else zLo <- 0 + // x_(j+n) <- x_(j+n) + (B - zLo) + // = x_(j+n) - zLo + B + // zHi <- zHi + 1 + // Here, invariant P1 still holds, P2 may break. + // (b) Advance j: + // Have z = zHi.B since zLo = 0. + // j <- j + 1 + // z <- zHi + f * a_(j+1) + // P2 holds: + // Have z <= B + (B-1)*(B-1) = B + B^2 - 2B + 1 = B^2 - B + 1 + // Have z <= B^2 - 1 when B >= 2 which is required for B being a base. + // P1 holds, + // moved f.digits_(j+1).B^(j+1+n) factor over. + // + // Once j+1 exceeds ad, summation is zero and it contributes no more terms (b). + // Continue until z = 0, which happens since z decreases towards 0. + // Done. + // + let invariant (_,_,_) = () + #if CHECKED + let x_initial = copyN x + let x_result = sub x_initial (shiftUp n (scale f a)) + let invariant (z,j,n) = + let P1 = + equal + x_result + (sub x (add (shiftUp (j+n) (embed64 z)) + (mul (embed f) + (shiftUp (j+1+n) (shiftDown (j+1) a))))) + let P2 = z < baseNi64 * baseNi64 - 1L + check P1 "P1"; + check P2 "P2" + + #endif + let xres = x + let x,xd = x.digits,degree x + let a,ad = a.digits,degree a + let f = int64 f + let mutable j = 0 + let mutable z = f * int64 a.[0] + while( z > 0L || j < ad ) do + if j > xd then failwith "scaleSubInPlace: pre-condition did not apply, result would be -ve"; + invariant(z,j,n); // P1,P2 hold + let mutable zLo = mod64base z |> int32 + let mutable zHi = div64base z + if zLo <= x.[j+n] then + x.[j+n] <- x.[j+n] - zLo + else ( + x.[j+n] <- x.[j+n] + (baseN - zLo); + zHi <- zHi + 1L + ); + // P1 holds + if j < ad then + z <- zHi + f * int64 a.[j+1] + else + z <- zHi; + j <- j + 1; + // P1,P2 hold + done; + ignore (normN xres) + + // + let scaleSub x f a n = + let freshx = add x zero + scaleSubInPlace freshx f a n; + normN freshx +(* + + let scaleSub2 x f a n = sub x (shiftUp n (mul (embed f) a)) + + let x = (mul (embed 234234234) (pow (embed 10) (embed 20))) + let f = 2 + let a = (embed 1231231231) + let n = 2 + let res = scaleSub x f a n + let res2 = scaleSub2 x f a n + + let x, xd, f, a, ad, n = freshx.digits, freshx.bound, f, a.digits, a.bound, n + *) + + +//---------------------------------------------------------------------------- +// division - scaleAddInPlace +//-------------------------------------------------------------------------- + + let scaleAddInPlace x f a n = + // Have x = sumR 0 xd (\i.xi.B^i) + // a = sumR 0 ad (\i.digitsi.B^i) + // f < B + // n < xd + // Required to add f.digits.B^n to x place. + // Assumes result will fit with x (0...xd). + //------ + // Result = x_initial + f.digits.B^n + // = x_initial + f.[sumR 0 ad (\i.digitsi.B^i+n)] + // State: j = 0 + // z = f * a0 + // Invariant(x,z,j,n): + // P1: x_result = x + z.B^(j+n) + f.[sumR (j+1) ad (\i.digitsi.B^i+n)] + // P2: z < B^2 - 1, and so has form z = zHi.B + zLo for zHi,zLo < B. + // Base: Invariant holds initially. + // Step: (a) Add zLo to x: + // If zLo < B - x_(j+n) then zLo <- 0 + // x_(j+n) <- x_(j+n) + zLo + // else zLo <- 0 + // x_(j+n) <- zLo - (B - x_(j+n)) + // = x_(j+n) + zLo - B + // zHi <- zHi + 1 + // Here, invariant P1 still holds, P2 may break. + // (b) Advance j: + // Have z = zHi.B since zLo = 0. + // j <- j + 1 + // z <- zHi + f * a_(j+1) + // P2 holds: + // Have z <= B + (B-1)*(B-1) = B + B^2 - 2B + 1 = B^2 - B + 1 + // Have z <= B^2 - 1 when B >= 2 which is required for B being a base. + // P1 holds, + // moved f.digits_(j+1).B^(j+1+n) factor over. + // + // Once j+1 exceeds ad, summation is zero and it contributes no more terms (b). + // Continue until z = 0, which happens since z decreases towards 0. + // Done. + // + let invariant (_,_,_) = () +#if CHECKED + let x_initial = copyN x + let x_result = add x_initial (shiftUp n (scale f a)) + let invariant (z,j,n) = + let P1 = + equal + x_result + (add x (add (shiftUp (j+n) (embed64 z)) + (mul (embed f) + (shiftUp (j+1+n) (shiftDown (j+1) a))))) + let P2 = z < baseNi64 * baseNi64 - 1L + check P1 "P1"; + check P2 "P2" + +#endif + let xres = x + let x,xd = x.digits,degree x + let a,ad = a.digits,degree a + let f = int64 f + let mutable j = 0 + let mutable z = f * int64 a.[0] + while( z > 0L || j < ad ) do + if j > xd then failwith "scaleSubInPlace: pre-condition did not apply, result would be -ve"; + invariant(z,j,n); // P1,P2 hold + let mutable zLo = mod64base z |> int32 + let mutable zHi = div64base z + if zLo < baseN - x.[j+n] then + x.[j+n] <- x.[j+n] + zLo + else ( + x.[j+n] <- zLo - (baseN - x.[j+n]); + zHi <- zHi + 1L + ); + // P1 holds + if j < ad then + z <- zHi + f * int64 a.[j+1] + else + z <- zHi; + j <- j + 1; + // P1,P2 hold + done; + ignore (normN xres) + + // + let scaleAdd x f a n = + let freshx = add x zero + scaleAddInPlace freshx f a n; + normN freshx + +(* + let scaleAdd2 x f a n = add x (shiftUp n (mul (embed f) a)) + + let x = (mul (embed 234234234) (pow (embed 10) (embed 20))) + let f = 2 + let a = (embed 1231231231) + let n = 2 + let res = scaleAdd x f a n + let res2 = scaleAdd2 x f a n + + let x, xd, f, a, ad, n = freshx.digits, freshx.bound, f, a.digits, a.bound, n +*) + +//---------------------------------------------------------------------------- +// division - removeFactor +//-------------------------------------------------------------------------- + + (* + let removeFactorReference x a n = + let ff = div x (shiftUp n a) + toInt ff + *) + + let removeFactor x a n = + // Assumes x < a.B^(n+1) + // Choose f s.t. + // (a) f.digits.B^n <= x + // (b) f=0 iff x < a.B^n + // + let dega,degx = degree a,degree x + if degx < dega + n then 0 else // possible with "normalisation" + let aa,xa = a.digits,x.digits + let f = + if dega = 0 then // a = a0 + if degx = n then + xa.[n] / aa.[0] + else ( +#if CHECKED + check (degx = n+1) "removeFactor degx#1"; +#endif + let f64 = (int64 xa.[degx] * baseNi64 + int64 xa.[degx-1]) / int64 aa.[0] + int32 f64 + ) + else // a = sumR 0 dega (\i.digitsi.B^i) + if degx = dega + n then + xa.[degx] / (aa.[dega] + 1) // +1 to bound above a + else ( +#if CHECKED + check (degx = dega+n+1) "removeFactor degx#2"; +#endif + let f64 = (int64 xa.[degx] * baseNi64 + int64 xa.[degx-1]) + / (int64 aa.[dega] + 1L) // +1 to bound above a + int32 f64 + ) + + if f = 0 then + let lte = (shiftCompare a n x 0) <> 1 + if lte then 1 else 0 + else + f + + +//---------------------------------------------------------------------------- +// division - divmod +//-------------------------------------------------------------------------- + + let divmod b a = + // Returns d,r where b = d.digits + r and r0 then + scaleSubInPlace x f a n; + scaleAddInPlace d f one n; + Invariant(d,x,n,p) + else + finished <- f=0 && n=0; + if not finished then + if p = m+n then + Invariant(d,x,n-1,p); + n <- n-1 + else + Invariant(d,x,n-1,p-1); + n <- n-1; + p <- p-1 + // Have: "b = d.digits + x" return d,x + normN d,normN x + + //---------------------------------------------------------------------------- + // div, mod + //-------------------------------------------------------------------------- + + [] + let div b a = fst (divmod b a) + [] + let rem b a = snd (divmod b a) + // rem b a, for small a can do (base mod a) trick - O(N) + + + //---------------------------------------------------------------------------- + // hcf + //-------------------------------------------------------------------------- + + let hcf a b = + // Have: 0 <= a,b since naturals + let rec hcfloop a b = // Require: 0 <= a <= b + if equal zero a then b + else + // Have: 0 < a <= b + let _,r = divmod b a + // Have: r < a from divmod + hcfloop r a // Have: 0 <= r < a + + if lt a b then hcfloop a b // Have: 0 <= a < b + else hcfloop b a // Have: 0 <= b <= a + + //---------------------------------------------------------------------------- + // pow + //-------------------------------------------------------------------------- + + let two = embed 2 + let powi x n = + let rec power acc x n = + if n=0 then acc + elif n % 2=0 then power acc (mul x x) (n / 2) + else power (mul x acc) (mul x x) (n / 2) + + power one x n + + let pow x n = + let rec power acc x n = + if isZero n then acc + else + let ndiv2,nmod2 = divmod n two // use: intdivmod when available + if isZero nmod2 then power acc (mul x x) ndiv2 + else power (mul x acc) (mul x x) ndiv2 + + power one x n + +//---------------------------------------------------------------------------- +// float n +//-------------------------------------------------------------------------- + + let toFloat n = + let basef = float baseN + let rec evalFloat acc k i = + if i = n.bound then + acc + else + evalFloat (acc + k * float n.digits.[i]) (k * basef) (i+1) + evalFloat 0.0 1.0 0 + +//---------------------------------------------------------------------------- +// n <-> int +//-------------------------------------------------------------------------- + + let ofInt32 n = embed n + let ofInt64 n = embed64 n + + /// Convert BigNat to uint32 otherwise OverflowException. + let toUInt32 n = + match n.bound with + | 0 -> 0u + | 1 -> n.digits.[0] |> uint32 + | 2 -> let xA,xB = n.digits.[0],n.digits.[1] + if xB > baseMask32B then raise (System.OverflowException()) + ( uint32 (xA &&& baseMask32A)) + + ((uint32 (xB &&& baseMask32B)) <<< baseShift32B) + | _ -> raise (System.OverflowException()) + + /// Convert BigNat to uint64 otherwise OverflowException. + let toUInt64 n = + match n.bound with + | 0 -> 0UL + | 1 -> n.digits.[0] |> uint64 + | 2 -> let xA,xB = n.digits.[0],n.digits.[1] + ( uint64 (xA &&& baseMask64A)) + + ((uint64 (xB &&& baseMask64B)) <<< baseShift64B) + | 3 -> let xA,xB,xC = n.digits.[0],n.digits.[1],n.digits.[2] + if xC > baseMask64C then raise (System.OverflowException()) + ( uint64 (xA &&& baseMask64A)) + + ((uint64 (xB &&& baseMask64B)) <<< baseShift64B) + + ((uint64 (xC &&& baseMask64C)) <<< baseShift64C) + | _ -> raise (System.OverflowException()) + + +//---------------------------------------------------------------------------- +// n -> string +//-------------------------------------------------------------------------- + + +#if CHECKED + let checks = false +#endif + let toString n = + // Much better complexity than naive_string_of_z. + // It still does "nDigit" calls to (int)divmod, + // but the degree on which it is called halves (not decrements) each time. + // + let degn = degree n + let rec route prior k ten2k = + if degree ten2k > degn + then (k,ten2k)::prior + else route ((k,ten2k)::prior) (k+1) (mul ten2k ten2k) + let kten2ks = route [] 0 (embed 10) + let rec collect isLeading digits n = function + | [] -> + // Have 0 <= n < 10^1, so collect a single digit (if needed) + let n = eval n +#if CHECKED + if checks then check (0 <= n) "toString: digit0"; + if checks then check (n <= 9) "toString: digit9"; +#endif + if isLeading && n=0 then digits // suppress leading 0 + else string n :: digits + | (_,ten2k)::prior -> +#if CHECKED + if checks then check (lt n (mul ten2k ten2k)) "string_of_int: bound n"; +#endif + // Have 0 <= n < (ten2k)^2 and ten2k = 10^(2^k) + let nH,nL = divmod n ten2k +#if CHECKED + if checks then check (lt nH ten2k) "string_of_int: bound nH"; + if checks then check (lt nL ten2k) "string_of_int: bound nL"; +#endif + // Have 0 <= nH,nL < (ten2k) and ten2k = 10^(2^k) + if isLeading && isZero nH then + // suppress leading 0s + let digits = collect isLeading digits nL prior + digits + else + let digits = collect false digits nL prior + let digits = collect isLeading digits nH prior + digits + + let prior = kten2ks + let digits = collect true [] n prior + match digits with + | [] -> "0" + | _ -> digits |> Array.ofList |> System.String.Concat + +//---------------------------------------------------------------------------- +// n <- string +//-------------------------------------------------------------------------- + + let ofString (str:string) = + // Would it be better to split string half and combine results? + let len = str.Length + if System.String.IsNullOrEmpty str then invalidArg "str" "empty string"; + let ten = embed 10 + let rec build acc i = + if i=len then + acc + else + let c = str.[i] + let d = int c - int '0' + if 0 <= d && d <= 9 then + build (add (mul ten acc) (embed d)) (i+1) + else + raise (new System.FormatException("The value could not be parsed")) + + build (embed 0) 0 + + let isSmall n = (n.bound <= 1) + let getSmall n = index n 0 + + //---------------------------------------------------------------------------- + // factorial + //-------------------------------------------------------------------------- + + let factorial n = + //***** + // Factorial(n) = 1.2.3.....(n-1).n + // + // Factorial is sometimes used as a test for multiplication. + // The QuickMul FFT multiplier takes over only when both operands reach a given size. + // How to compute factorial? + // + // (a) Factorial(n) = factorial(n-1).n + // This is unlikely to make use of the FFT (n never large enough). + // (b) Factorial(n) = (1.2.3.4....k) . (k.[k+1]...(n-1).n) + // Applied recursively QuickMul FFT will take over on large products. + // + //**** + let rec productR a b = + if equal a b then a + else + let m = div (add a b) (ofInt32 2) + mul (productR a m) (productR (add m one) b) + + productR one n + + +#endif diff --git a/src/FSharpPortable/n.fsi b/src/FSharpPortable/n.fsi new file mode 100644 index 00000000..bf01cbba --- /dev/null +++ b/src/FSharpPortable/n.fsi @@ -0,0 +1,71 @@ +// First version copied from the F# compiler sources +// https://raw.github.com/fsharp/fsharp/master/src/fsharp/FSharp.Core/math/n.fsi + +//---------------------------------------------------------------------------- +// Copyright (c) 2002-2012 Microsoft Corporation. +// +// This source code is subject to terms and conditions of the Apache License, Version 2.0. A +// copy of the license can be found in the License.html file at the root of this distribution. +// By using this source code in any fashion, you are agreeing to be bound +// by the terms of the Apache License, Version 2.0. +// +// You must not remove this notice, or any other, from this software. +//---------------------------------------------------------------------------- + +namespace Microsoft.FSharp.Math + +#if FX_NO_BIGINT +open Microsoft.FSharp.Collections +open Microsoft.FSharp.Core + +/// Abstract internal type +[] +type internal BigNat + +module internal BigNatModule = + + val zero : BigNat + val one : BigNat + val two : BigNat + + val add : BigNat -> BigNat -> BigNat + val sub : BigNat -> BigNat -> BigNat + val mul : BigNat -> BigNat -> BigNat + val divmod : BigNat -> BigNat -> BigNat * BigNat + val div : BigNat -> BigNat -> BigNat + val rem : BigNat -> BigNat -> BigNat + val hcf : BigNat -> BigNat -> BigNat + + val min : BigNat -> BigNat -> BigNat + val max : BigNat -> BigNat -> BigNat + val scale : int -> BigNat -> BigNat + val powi : BigNat -> int -> BigNat + val pow : BigNat -> BigNat -> BigNat + + val IsZero : BigNat -> bool + val isZero : BigNat -> bool + val isOne : BigNat -> bool + val equal : BigNat -> BigNat -> bool + val compare : BigNat -> BigNat -> int + val lt : BigNat -> BigNat -> bool + val gt : BigNat -> BigNat -> bool + val lte : BigNat -> BigNat -> bool + val gte : BigNat -> BigNat -> bool + + val hash : BigNat -> int + val toFloat : BigNat -> float + val ofInt32 : int -> BigNat + val ofInt64 : int64 -> BigNat + val toString : BigNat -> string + val ofString : string -> BigNat + + val toUInt32 : BigNat -> uint32 + val toUInt64 : BigNat -> uint64 + + val factorial : BigNat -> BigNat + // val randomBits : int -> BigNat + val bits : BigNat -> int + val isSmall : BigNat -> bool (* will fit in int32 (but not nec all int32) *) + val getSmall : BigNat -> int32 (* get the value, if it satisfies isSmall *) + +#endif diff --git a/src/FSharpPortable/z.fs b/src/FSharpPortable/z.fs new file mode 100644 index 00000000..b2d56594 --- /dev/null +++ b/src/FSharpPortable/z.fs @@ -0,0 +1,411 @@ +// First version copied from the F# compiler sources +// https://raw.github.com/fsharp/fsharp/master/src/fsharp/FSharp.Core/math/z.fs + +//---------------------------------------------------------------------------- +// Copyright (c) 2002-2012 Microsoft Corporation. +// +// This source code is subject to terms and conditions of the Apache License, Version 2.0. A +// copy of the license can be found in the License.html file at the root of this distribution. +// By using this source code in any fashion, you are agreeing to be bound +// by the terms of the Apache License, Version 2.0. +// +// You must not remove this notice, or any other, from this software. +//---------------------------------------------------------------------------- + +#nowarn "44" // This construct is deprecated. This function is for use by compiled F# code and should not be used directly +namespace System.Numerics + +#if FX_NO_BIGINT + open Microsoft.FSharp.Collections + open Microsoft.FSharp.Core + open Microsoft.FSharp.Core.Operators + open Microsoft.FSharp.Core.LanguagePrimitives.IntrinsicOperators + open Microsoft.FSharp.Primitives.Basics + open Microsoft.FSharp.Math + open System + open System.Globalization + + + // INVARIANT: signInt = 1 or -1 + // value(z) = signInt * v + // NOTE: 0 has two repns (+1,0) or (-1,0). + [] + [] + [] + type BigInteger(signInt:int, v : BigNat) = + + static let smallLim = 4096 + static let smallPosTab = Array.init smallLim BigNatModule.ofInt32 + static let one = BigInteger(1) + static let zero = BigInteger(0) + + static member internal nat n = + if BigNatModule.isSmall n && BigNatModule.getSmall n < smallLim + then smallPosTab.[BigNatModule.getSmall n] + else n + static member internal create (s,n) = BigInteger(s,BigInteger.nat n) + static member internal posn n = BigInteger(1,BigInteger.nat n) + static member internal negn n = BigInteger(-1,BigInteger.nat n) + + + member x.Sign = if x.IsZero then 0 else signInt + member x.SignInt = signInt + member internal x.V = v + + static member op_Equality (x:BigInteger, y:BigInteger) = + //System.Console.WriteLine("x = {0}",box x) + //System.Console.WriteLine("y = {0}",box y) + match x.SignInt,y.SignInt with + | 1, 1 -> BigNatModule.equal x.V y.V // +1.xv = +1.yv iff xv = yv + | -1, -1 -> BigNatModule.equal x.V y.V // -1.xv = -1.yv iff xv = yv + | 1,-1 -> BigNatModule.isZero x.V && BigNatModule.isZero y.V // 1.xv = -1.yv iff xv=0 and yv=0 + | -1, 1 -> BigNatModule.isZero x.V && BigNatModule.isZero y.V // -1.xv = 1.yv iff xv=0 and yv=0 + | _ -> invalidArg "x" "signs should be +/- 1" + + static member op_Inequality (x:BigInteger, y:BigInteger) = not (BigInteger.op_Equality(x,y)) // CA2226: OperatorsShouldHaveSymmetricalOverloads + + static member op_LessThan (x:BigInteger, y:BigInteger) = + match x.SignInt,y.SignInt with + | 1, 1 -> BigNatModule.lt x.V y.V // 1.xv < 1.yv iff xv < yv + | -1,-1 -> BigNatModule.lt y.V x.V // -1.xv < -1.yv iff yv < xv + | 1,-1 -> false // 1.xv < -1.yv iff 0 <= 1.xv < -1.yv <= 0 iff false + | -1, 1 -> not (BigNatModule.isZero x.V) || not (BigNatModule.isZero y.V) + // -1.xv < 1.yv + // (a) xv=0 and yv=0, then false + // (b) xv<>0, -1.xv < 0 <= 1.yv, so true + // (c) yv<>0, -1.xv <= 0 < 1.yv, so true + | _ -> invalidArg "x" "signs should be +/- 1" + + static member op_GreaterThan (x:BigInteger, y:BigInteger) = // Follow lt by +/- symmetry + match x.SignInt,y.SignInt with + | 1, 1 -> BigNatModule.gt x.V y.V + | -1,-1 -> BigNatModule.gt y.V x.V + | 1,-1 -> not (BigNatModule.isZero x.V) || not (BigNatModule.isZero y.V) + | -1, 1 -> false + | _ -> invalidArg "x" "signs should be +/- 1" + + static member internal compare(n,nn) = if BigInteger.op_LessThan(n,nn) then -1 elif BigInteger.op_Equality(n,nn) then 0 else 1 + static member internal hash (z:BigInteger) = z.SignInt + BigNatModule.hash(z.V) + + override x.ToString() = + match x.SignInt with + | 1 -> BigNatModule.toString x.V // positive + | -1 -> + if BigNatModule.isZero x.V + then "0" // not negative infact, but zero. + else "-" + BigNatModule.toString x.V // negative + | _ -> invalidOp "signs should be +/- 1" + + member x.StructuredDisplayString = x.ToString() + + interface System.IComparable with + member this.CompareTo(obj:obj) = + match obj with + | :? BigInteger as that -> BigInteger.compare(this,that) + | _ -> invalidArg "obj" "the objects are not comparable" + + override this.Equals(obj) = + match obj with + | :? BigInteger as that -> BigInteger.op_Equality(this, that) + | _ -> false + + override x.GetHashCode() = BigInteger.hash(x) + + + new (n:int) = + if n>=0 + then BigInteger (1,BigInteger.nat(BigNatModule.ofInt32 n)) + elif (n = System.Int32.MinValue) + then BigInteger(-1,BigInteger.nat(BigNatModule.ofInt64 (-(int64 n)))) + else BigInteger(-1,BigInteger.nat(BigNatModule.ofInt32 (-n))) + + + new (n:int64) = + if n>=0L + then BigInteger(1,BigInteger.nat (BigNatModule.ofInt64 n)) + elif (n = System.Int64.MinValue) + then BigInteger(-1,BigInteger.nat (BigNatModule.add (BigNatModule.ofInt64 System.Int64.MaxValue) BigNatModule.one) ) + else BigInteger(-1,BigInteger.nat (BigNatModule.ofInt64 (-n))) + + static member One = one + static member Zero = zero + static member (~-) (z:BigInteger) = BigInteger.create(-1 * z.SignInt,z.V) + static member Scale(k,z:BigInteger) = + if k<0 + then BigInteger.create(-z.SignInt, (BigNatModule.scale (-k) z.V)) // k.zsign.zv = -zsign.(-k.zv) + else BigInteger.create(z.SignInt, (BigNatModule.scale k z.V)) // k.zsign.zv = zsign.k.zv + + // Result: 1.nx - 1.ny (integer subtraction) + static member internal subnn (nx,ny) = + if BigNatModule.gte nx ny + then BigInteger.posn (BigNatModule.sub nx ny) // nx >= ny, result +ve, +1.(nx - ny) + else BigInteger.negn (BigNatModule.sub ny nx) // nx < ny, result -ve, -1.(ny - nx) + + static member internal addnn (nx,ny) = + BigInteger.posn (BigNatModule.add nx ny) // Compute "nx + ny" to be integer + + member x.IsZero = BigNatModule.isZero x.V // signx.xv = 0 iff xv=0, since signx is +1,-1 + member x.IsOne = (x.SignInt = 1) && BigNatModule.isOne x.V // signx.xv = 1 iff signx = +1 and xv = 1 + static member (+) (x:BigInteger,y:BigInteger) = + if y.IsZero then x else + if x.IsZero then y else + match x.SignInt,y.SignInt with + | 1, 1 -> BigInteger.addnn(x.V,y.V) // 1.xv + 1.yv = (xv + yv) + | -1,-1 -> -(BigInteger.addnn(x.V,y.V)) // -1.xv + -1.yv = -(xv + yv) + | 1,-1 -> BigInteger.subnn (x.V,y.V) // 1.xv + -1.yv = (xv - yv) + | -1, 1 -> BigInteger.subnn(y.V,x.V) // -1.xv + 1.yv = (yv - xv) + | _ -> invalidArg "x" "signs should be +/- 1" + + static member (-) (x:BigInteger,y:BigInteger) = + if y.IsZero then x else + match x.SignInt,y.SignInt with + | 1, 1 -> BigInteger.subnn(x.V,y.V) // 1.xv - 1.yv = (xv - yv) + | -1,-1 -> BigInteger.subnn(y.V,x.V) // -1.xv - -1.yv = (yv - xv) + | 1,-1 -> BigInteger.addnn(x.V,y.V) // 1.xv - -1.yv = (xv + yv) + | -1, 1 -> -(BigInteger.addnn(x.V,y.V)) // -1.xv - 1.yv = -(xv + yv) + | _ -> invalidArg "x" "signs should be +/- 1" + + static member ( * ) (x:BigInteger,y:BigInteger) = + if x.IsZero then x + elif y.IsZero then y + elif x.IsOne then y + elif y.IsOne then x + else + let m = (BigNatModule.mul x.V y.V) + BigInteger.create (x.SignInt * y.SignInt,m) // xsign.xv * ysign.yv = (xsign.ysign).(xv.yv) + + static member DivRem (x:BigInteger,y:BigInteger,rem:BigInteger byref) = + let d,r = BigNatModule.divmod x.V y.V + // HAVE: |x| = d.|y| + r and 0 <= r < |y| + // HAVE: xv = d.yv + r and 0 <= r < yv + match x.SignInt,y.SignInt with + | 1, 1 -> rem <- BigInteger.posn r ; BigInteger.posn d // 1.xv = 1.d.( 1.yv) + ( 1.r) + | -1,-1 -> rem <- BigInteger.negn r ; BigInteger.posn d // -1.xv = 1.d.(-1.yv) + (-1.r) + | 1,-1 -> rem <- BigInteger.posn r ; BigInteger.negn d // 1.xv = -1.d.(-1.yv) + ( 1.r) + | -1, 1 -> rem <- BigInteger.negn r ; BigInteger.negn d // -1.xv = -1.d.( 1.yv) + (-1.r) + | _ -> invalidArg "x" "signs should be +/- 1" + + static member (/) (x:BigInteger,y:BigInteger) = + let mutable rem = new BigInteger(0) + BigInteger.DivRem(x,y,&rem) + static member (%) (x:BigInteger,y:BigInteger) = + let mutable rem = new BigInteger(0) + BigInteger.DivRem(x,y,&rem) |> ignore ; rem + static member GreatestCommonDivisor (x:BigInteger,y:BigInteger) = BigInteger.posn (BigNatModule.hcf x.V y.V) // hcf (xsign.xv,ysign.yv) = hcf (xv,yv) + + member x.IsNegative = x.SignInt = -1 && not (x.IsZero) // signx.xv < 0 iff signx = -1 and xv<>0 + member x.IsPositive = x.SignInt = 1 && not (x.IsZero) // signx.xv > 0 iff signx = +1 and xv<>0 + static member Abs (x:BigInteger) = if x.SignInt = -1 then -x else x + + static member op_LessThanOrEqual (x:BigInteger,y:BigInteger) = + match x.SignInt,y.SignInt with + | 1, 1 -> BigNatModule.lte x.V y.V // 1.xv <= 1.yv iff xv <= yv + | -1,-1 -> BigNatModule.lte y.V x.V // -1.xv <= -1.yv iff yv <= xv + | 1,-1 -> BigNatModule.isZero x.V && BigNatModule.isZero y.V // 1.xv <= -1.yv, + // (a) if xv=0 and yv=0 then true + // (b) otherwise false, only meet at zero. + + | -1, 1 -> true // -1.xv <= 1.yv, true + | _ -> invalidArg "x" "signs should be +/- 1" + + static member op_GreaterThanOrEqual (x:BigInteger,y:BigInteger) = // Follow lte by +/- symmetry + match x.SignInt,y.SignInt with + | 1, 1 -> BigNatModule.gte x.V y.V + | -1,-1 -> BigNatModule.gte y.V x.V + | 1,-1 -> true + | -1, 1 -> BigNatModule.isZero x.V && BigNatModule.isZero y.V + | _ -> invalidArg "x" "signs should be +/- 1" + + + static member Pow (x:BigInteger,y:int32) = + if y < 0 then invalidArg "y" "The input must be non-negative." + let yval = BigInteger(y) + BigInteger.create ((if BigNatModule.isZero (BigNatModule.rem yval.V BigNatModule.two) then 1 else x.SignInt), BigNatModule.pow x.V yval.V) + + static member op_Explicit (x:BigInteger) = + let u = BigNatModule.toUInt32 x.V + if u <= uint32 System.Int32.MaxValue then + // Handle range [-MaxValue,MaxValue] + x.SignInt * int32 u + elif x.SignInt = -1 && u = uint32 (System.Int32.MaxValue + 1) then + //assert(System.Int32.MinValue = 0 - System.Int32.MaxValue - 1) + // Handle MinValue = -(MaxValue+1) special case not covered by the above + System.Int32.MinValue + else + raise (System.OverflowException()) + + static member op_Explicit (x:BigInteger) = + let u = BigNatModule.toUInt64 x.V + if u <= uint64 System.Int64.MaxValue then + (* Handle range [-MaxValue,MaxValue] *) + int64 x.SignInt * int64 u + elif x.SignInt = -1 && u = uint64 (System.Int64.MaxValue + 1L) then + //assert(System.Int64.MinValue = 0 - System.Int64.MaxValue - 1L) + (* Handle MinValue = -(MaxValue+1) special case not covered by the above *) + System.Int64.MinValue + else + raise (System.OverflowException()) + + static member op_Explicit (x:BigInteger) = + match x.SignInt with + | 1 -> BigNatModule.toFloat x.V // float (1.xv) = float (xv) + | -1 -> - (BigNatModule.toFloat x.V) // float (-1.xv) = - float (xv) + | _ -> invalidArg "x" "signs should be +/- 1" + + static member Parse(text:string) = + let len = text.Length + if len = 0 then raise (new System.FormatException("The value could not be parsed")) + if text.[0..0] = "-" then + BigInteger.negn (BigNatModule.ofString text.[1..len-1]) + else + BigInteger.posn (BigNatModule.ofString text) + + member internal x.IsSmall = BigNatModule.isSmall (x.V) + static member Factorial (x:BigInteger) = + if x.IsNegative then invalidArg "x" "The input must be non-negative." + if x.IsPositive then BigInteger.posn (BigNatModule.factorial x.V) + else BigInteger.One + + static member ( ~+ )(n1:BigInteger) = n1 + + static member FromInt64(x:int64) = new BigInteger(x) + static member FromInt32(x:int32) = new BigInteger(x) +#endif + +namespace Microsoft.FSharp.Core + + + type bigint = System.Numerics.BigInteger + + open System + open System.Diagnostics.CodeAnalysis + open System.Globalization + open Microsoft.FSharp.Core.Operators + open Microsoft.FSharp.Core.LanguagePrimitives.IntrinsicOperators + open System.Numerics + +#if FX_NO_BIGINT + // FxCop suppressions + [] + [] + [] + [] + [] + [] + [] + [] + [] + [] + [] + do() +#endif + + [] + module NumericLiterals = + + module NumericLiteralI = + + let tab64 = new System.Collections.Generic.Dictionary() + let tabParse = new System.Collections.Generic.Dictionary() + + let FromInt64Dynamic (x64:int64) : obj = + lock tab64 (fun () -> + let mutable res = Unchecked.defaultof<_> + let ok = tab64.TryGetValue(x64,&res) + if ok then res else + res <- BigInteger(x64) + tab64.[x64] <- res + res) + + let inline get32 (x32:int32) = FromInt64Dynamic (int64 x32) + + let inline isOX s = not (System.String.IsNullOrEmpty(s)) && s.Length > 2 && s.[0] = '0' && s.[1] = 'x' + + let FromZero () : 'T = + (get32 0 :?> 'T) + + let FromOne () : 'T = + (get32 1 :?> 'T) + + let FromInt32 (i:int32): 'T = + (get32 i :?> 'T) + + let FromInt64 (i:int64): 'T = + (FromInt64Dynamic i :?> 'T) + + let getParse s = + lock tabParse (fun () -> + let mutable res = Unchecked.defaultof<_> + let ok = tabParse.TryGetValue(s,&res) + if ok then + res + else +#if FSHARP_CORE_PORTABLE + // SL5 (and therefore Portable Profile47) does not have Parse, so make our own simple implementation + let parse(s : string) = + // ws* sign? digits+ ws* + let mutable i = 0 + // leading whitespace + while i < s.Length && System.Char.IsWhiteSpace(s.[i]) do + i <- i + 1 + if i = s.Length then + raise <| new System.ArgumentException() + // optional sign + let mutable isNegative = false + if s.[i] = '+' then + i <- i + 1 + elif s.[i] = '-' then + isNegative <- true + i <- i + 1 + if i = s.Length then + raise <| new System.ArgumentException() + // digits + let startDigits = i + while i < s.Length && System.Char.IsDigit(s.[i]) do + i <- i + 1 + let endDigits = i + let len = endDigits - startDigits + if len = 0 then + raise <| new System.ArgumentException() + // trailing whitespace + while i < s.Length && System.Char.IsWhiteSpace(s.[i]) do + i <- i + 1 + if i <> s.Length then + raise <| new System.ArgumentException() + // text is now valid, parse it + let mutable r = new System.Numerics.BigInteger(int(s.[startDigits]) - int('0')) + let ten = new System.Numerics.BigInteger(10) + for j in startDigits+1 .. endDigits-1 do + r <- r * ten + r <- r + new System.Numerics.BigInteger(int(s.[j]) - int('0')) + if isNegative then + r <- new System.Numerics.BigInteger(0) - r + r + let v = parse s +#else + let v = +#if FX_NO_BIGINT + BigInteger.Parse s +#else +#if FX_NO_BIGINT_CULTURE_PARSE + BigInteger.Parse s +#else + if isOX s then + BigInteger.Parse (s.[2..],NumberStyles.AllowHexSpecifier,CultureInfo.InvariantCulture) + else + BigInteger.Parse (s,NumberStyles.AllowLeadingSign,CultureInfo.InvariantCulture) +#endif +#endif +#endif + res <- v + tabParse.[s] <- res + res) + + let FromStringDynamic (s:string) : obj = + getParse s + + let FromString (s:string) : 'T = + (FromStringDynamic s :?> 'T) + + + diff --git a/src/FSharpPortable/z.fsi b/src/FSharpPortable/z.fsi new file mode 100644 index 00000000..cbc66872 --- /dev/null +++ b/src/FSharpPortable/z.fsi @@ -0,0 +1,124 @@ +// First version copied from the F# compiler sources +// https://raw.github.com/fsharp/fsharp/master/src/fsharp/FSharp.Core/math/z.fsi + +//---------------------------------------------------------------------------- +// Copyright (c) 2002-2012 Microsoft Corporation. +// +// This source code is subject to terms and conditions of the Apache License, Version 2.0. A +// copy of the license can be found in the License.html file at the root of this distribution. +// By using this source code in any fashion, you are agreeing to be bound +// by the terms of the Apache License, Version 2.0. +// +// You must not remove this notice, or any other, from this software. +//---------------------------------------------------------------------------- + +namespace System.Numerics +#if FX_NO_BIGINT + + open System + open Microsoft.FSharp.Collections + open Microsoft.FSharp.Core + + /// The type of arbitrary-sized integers + [] + [] + type BigInteger = + /// Return the sum of two big integers + static member ( + ) : x:BigInteger * y:BigInteger -> BigInteger + /// Return the modulus of big integers + static member ( % ) : x:BigInteger * y:BigInteger -> BigInteger + /// Return the product of big integers + static member ( * ) : x:BigInteger * y:BigInteger -> BigInteger + /// Return the difference of two big integers + static member ( - ) : x:BigInteger * y:BigInteger -> BigInteger + /// Return the ratio of two big integers + static member ( / ) : x:BigInteger * y:BigInteger -> BigInteger + /// Return the negation of a big integer + static member (~-) : x:BigInteger -> BigInteger + /// Return the given big integer + static member (~+) : x:BigInteger -> BigInteger + /// Convert a big integer to a floating point number + static member op_Explicit : x:BigInteger -> float + /// Convert a big integer to a 64-bit signed integer + static member op_Explicit : x:BigInteger -> int64 + /// Convert a big integer to a 32-bit signed integer + static member op_Explicit : x:BigInteger -> int32 + /// Parse a big integer from a string format + static member Parse : text:string -> BigInteger + /// Return the sign of a big integer: 0, +1 or -1 + member Sign : int + /// Compute the ratio and remainder of two big integers + static member DivRem : x:BigInteger * y:BigInteger * rem:BigInteger byref -> BigInteger + + /// This operator is for consistency when this type be used from other CLI languages + static member op_LessThan : x:BigInteger * y:BigInteger -> bool + /// This operator is for consistency when this type be used from other CLI languages + static member op_LessThanOrEqual : x:BigInteger * y:BigInteger -> bool + /// This operator is for consistency when this type be used from other CLI languages + static member op_GreaterThan : x:BigInteger * y:BigInteger -> bool + /// This operator is for consistency when this type be used from other CLI languages + static member op_GreaterThanOrEqual : x:BigInteger * y:BigInteger -> bool + /// This operator is for consistency when this type be used from other CLI languages + static member op_Equality : x:BigInteger * y:BigInteger -> bool + /// This operator is for consistency when this type be used from other CLI languages + static member op_Inequality : x:BigInteger * y:BigInteger -> bool + + /// Return the greatest common divisor of two big integers + static member GreatestCommonDivisor : x:BigInteger * y:BigInteger -> BigInteger + /// Return n^m for two big integers + static member Pow : x:BigInteger * y:int32 -> BigInteger + /// Compute the absolute value of a big integer + static member Abs : x:BigInteger -> BigInteger + /// Get the big integer for zero + static member Zero : BigInteger + /// Get the big integer for one + static member One : BigInteger + + /// Return true if a big integer is 'zero' + member IsZero : bool + /// Return true if a big integer is 'one' + member IsOne : bool + interface System.IComparable + override Equals : obj -> bool + override GetHashCode : unit -> int + override ToString : unit -> string + + /// Construct a BigInteger value for the given integer + new : x:int -> BigInteger + /// Construct a BigInteger value for the given 64-bit integer + new : x:int64 -> BigInteger +#if SILVERLIGHT + /// Provides custom formatting for BigInteger values + member StructuredDisplayString : string // This needs to be accessible in order for formatting code to invoke it via reflection in Silverlight. +#endif +#endif + + +namespace Microsoft.FSharp.Core + + type bigint = System.Numerics.BigInteger + + [] + /// Provides a default implementations of F# numeric literal syntax for literals fo the form 'dddI' + module NumericLiterals = + + /// Provides a default implementations of F# numeric literal syntax for literals fo the form 'dddI' + module NumericLiteralI = + open System.Numerics + + /// Provides a default implementations of F# numeric literal syntax for literals fo the form 'dddI' + val FromZero : value:unit -> 'T + /// Provides a default implementations of F# numeric literal syntax for literals fo the form 'dddI' + val FromOne : value:unit -> 'T + /// Provides a default implementations of F# numeric literal syntax for literals fo the form 'dddI' + val FromInt32 : value:int32 -> 'T + /// Provides a default implementations of F# numeric literal syntax for literals fo the form 'dddI' + val FromInt64 : value:int64 -> 'T + /// Provides a default implementations of F# numeric literal syntax for literals fo the form 'dddI' + val FromString : text:string -> 'T + /// Provides a default implementations of F# numeric literal syntax for literals fo the form 'dddI' + val FromInt64Dynamic : value:int64 -> obj + /// Provides a default implementations of F# numeric literal syntax for literals fo the form 'dddI' + val FromStringDynamic : text:string -> obj + +