// Copyright (c) Microsoft Corporation.
// Licensed under the MIT License.

namespace Microsoft.Quantum.Math {
    open Microsoft.Quantum.Convert;
    open Microsoft.Quantum.Diagnostics;

    //
    // Constants PI, E, LogOf2.
    //

    /// # Summary
    /// Represents the ratio of the circumference of a circle to its diameter.
    ///
    /// # Output
    /// A double-precision approximation of the circumference of a circle
    /// to its diameter, π ≈ 3.14159265358979323846.
    ///
    /// # See Also
    /// - Microsoft.Quantum.Math.E
    function PI() : Double {
        3.14159265358979323846
    }

    /// # Summary
    /// Returns the natural logarithmic base to double-precision.
    ///
    /// # Output
    /// A double-precision approximation of the natural logarithmic base,
    /// e ≈ 2.7182818284590452354.
    ///
    /// # See Also
    /// - Microsoft.Quantum.Math.PI
    function E() : Double {
        2.7182818284590452354
    }

    /// # Summary
    /// Returns the natural logarithm of 2.
    ///
    /// # Output
    /// Returns a `Double` equal to 0.6931471805599453.
    function LogOf2() : Double {
        0.6931471805599453
    }

    //
    // Special numbers in IEEE floating-point representation
    //

    /// # Summary
    /// Returns whether a given floating-point value is not a number (i.e. is
    /// NaN).
    ///
    /// # Input
    /// ## d
    /// A floating-point value to be checked.
    ///
    /// # Output
    /// `true` if and only if `d` is not a number.
    function IsNaN(d : Double) : Bool {
        return d != d;
    }

    /// # Summary
    /// Returns whether a given floating-point value is either positive or
    /// negative infinity.
    ///
    /// # Input
    /// ## d
    /// The floating-point value to be checked.
    ///
    /// # Output
    /// `true` if and only if `d` is either positive or negative infinity.
    ///
    /// # Remarks
    /// `NaN` is not a number, and is thus neither a finite number nor
    /// is it infinite. As such, `IsInfinite(0.0 / 0.0)` returns `false`.
    /// To check if a value is `NaN`, use `IsNaN(d)`.
    ///
    /// Note that even though this function returns `true` for both
    /// positive and negative infinities, these values can still be
    /// discriminated by checking `d > 0.0` and `d < 0.0`.
    ///
    /// # Example
    /// ```qsharp
    /// Message($"{IsInfinite(42.0)}"); // false
    /// Message($"{IsInfinite(0.0 / 0.0)}"); // false
    /// Message($"{IsInfinite(-1.0 / 0.0}"); // true
    /// ```
    ///
    /// # See Also
    /// - Microsoft.Quantum.Math.IsNaN
    function IsInfinite(d : Double) : Bool {
        return d == 1.0 / 0.0 or d == -1.0 / 0.0;
    }

    //
    // Sign, Abs, Min, Max, etc.
    //

    /// # Summary
    /// Returns -1, 0 or +1 that indicates the sign of a number.
    function SignI(a : Int) : Int {
        if (a < 0) {
            -1
        } elif (a > 0) {
            + 1
        } else {
            0
        }
    }

    /// # Summary
    /// Returns -1, 0 or +1 that indicates the sign of a number.
    function SignD(a : Double) : Int {
        if (a < 0.0) {
            -1
        } elif (a > 0.0) {
            + 1
        } else {
            0
        }
    }

    /// # Summary
    /// Returns -1, 0 or +1 that indicates the sign of a number.
    function SignL(a : BigInt) : Int {
        if (a < 0L) { 
            -1
        } elif (a > 0L) {
            + 1
        } else {
            0
        }
    }

    /// # Summary
    /// Returns the absolute value of an integer.
    function AbsI(a : Int) : Int {
        a < 0 ? -a | a
    }

    /// # Summary
    /// Returns the absolute value of a double-precision floating-point number.
    function AbsD(a : Double) : Double {
        a < 0.0 ? -a | a
    }

    /// # Summary
    function AbsL(a : BigInt) : BigInt {
        a < 0L ? -a | a
    }

    /// # Summary
    /// Returns the larger of two specified numbers.
    function MaxI(a : Int, b : Int) : Int {
        a > b ? a | b
    }

    /// # Summary
    /// Returns the larger of two specified numbers.
    function MaxD(a : Double, b : Double) : Double {
        a > b ? a | b
    }

    /// # Summary
    /// Returns the larger of two specified numbers.
    function MaxL(a : BigInt, b : BigInt) : BigInt {
        a > b ? a | b
    }

    /// # Summary
    /// Returns the smaller of two specified numbers.
    function MinI(a : Int, b : Int) : Int {
        a < b ? a | b
    }

    /// # Summary
    /// Returns the smaller of two specified numbers.
    function MinD(a : Double, b : Double) : Double {
        a < b ? a | b
    }

    /// # Summary
    /// Returns the smaller of two specified numbers.
    function MinL(a : BigInt, b : BigInt) : BigInt {
        a < b ? a | b
    }

    /// # Summary
    /// Given an array of integers, returns the largest element.
    ///
    /// # Input
    /// ## values
    /// An array to take the maximum of.
    ///
    /// # Output
    /// The largest element of `values`.
    function Max(values : Int[]) : Int {
        Fact(Length(values) > 0, "Array must contain at least one element.");
        mutable max = values[0];
        for element in values[1...] {
            if element > max {
                set max = element;
            }
        }

        max
    }

    /// # Summary
    /// Given an array of integers, returns the smallest element.
    ///
    /// # Input
    /// ## values
    /// An array to take the minimum of.
    ///
    /// # Output
    /// The smallest element of `values`.
    function Min(values : Int[]) : Int {
        Fact(Length(values) > 0, "Array must contain at least one element.");
        mutable min = values[0];
        for element in values[1...] {
            if element < min {
                set min = element;
            }
        }

        min
    }

    //
    // Trigonometric functions
    //

    /// # Summary
    /// Returns the angle whose cosine is the specified number.
    function ArcCos(x : Double) : Double {
        body intrinsic;
    }

    /// # Summary
    /// Returns the angle whose sine is the specified number.
    function ArcSin(y : Double) : Double {
        body intrinsic;
    }

    /// # Summary
    /// Returns the angle whose tangent is the specified number.
    function ArcTan(d : Double) : Double {
        body intrinsic;
    }

    /// # Summary
    /// Returns the angle whose tangent is the quotient of two specified numbers.
    function ArcTan2(y : Double, x : Double) : Double {
        body intrinsic;
    }

    /// # Summary
    /// Returns the cosine of the specified angle.
    function Cos(theta : Double) : Double {
        body intrinsic;
    }

    /// # Summary
    /// Returns the hyperbolic cosine of the specified angle.
    function Cosh(d : Double) : Double {
        body intrinsic;
    }

    /// # Summary
    /// Returns the sine of the specified angle.
    function Sin(theta : Double) : Double {
        body intrinsic;
    }

    /// # Summary
    /// Returns the hyperbolic sine of the specified angle.
    function Sinh(d : Double) : Double {
        body intrinsic;
    }

    /// # Summary
    /// Returns the tangent of the specified angle.
    function Tan(d : Double) : Double {
        body intrinsic;
    }

    /// # Summary
    /// Returns the hyperbolic tangent of the specified angle.
    function Tanh(d : Double) : Double {
        body intrinsic;
    }

    /// # Summary
    /// Computes the inverse hyperbolic cosine of a number.
    function ArcCosh(x : Double) : Double {
        Log(x + Sqrt(x * x - 1.0))
    }

    /// # Summary
    /// Computes the inverse hyperbolic sine of a number.
    function ArcSinh(x : Double) : Double {
        Log(x + Sqrt(x * x + 1.0))
    }


    /// # Summary
    /// Computes the inverse hyperbolic tangent of a number.
    function ArcTanh(x : Double) : Double {
        Log((1.0 + x) / (1.0 - x)) * 0.5
    }

    //
    // Sqrt, Log, exp, etc.
    //

    /// # Summary
    /// Returns the square root of a specified number.
    function Sqrt(d : Double) : Double {
        body intrinsic;
    }

    /// # Summary
    /// Returns the natural (base _e_) logarithm of a specified number.
    function Log(input : Double) : Double {
        body intrinsic;
    }

    /// # Summary
    /// Returns the base-10 logarithm of a specified number.
    function Log10(input : Double) : Double {
        Log(input) / Log(10.0)
    }

    /// # Summary
    /// Computes the base-2 logarithm of a number.
    function Lg(input : Double) : Double {
        Log(input) / Log(2.0)
    }

    //
    // Truncation and Rounding
    //

    /// # Summary
    /// Returns the integral part of a number.
    /// For example: Truncate(3.7) = 3; Truncate(-3.7) = -3
    function Truncate(value : Double) : Int {
        body intrinsic;
    }

    internal function ExtendedTruncation(value : Double) : (Int, Double, Bool) {
        let truncated = Truncate(value);
        (truncated, IntAsDouble(truncated) - value, value >= 0.0)
    }

    /// # Summary
    /// Returns the smallest integer greater than or equal to the specified number.
    /// For example: Ceiling(3.1) = 4; Ceiling(-3.7) = -3
    function Ceiling(value : Double) : Int {
        let (truncated, remainder, isPositive) = ExtendedTruncation(value);
        if AbsD(remainder) <= 1e-15 {
            truncated
        } else {
            isPositive ? truncated + 1 | truncated
        }
    }

    /// # Summary
    /// Returns the largest integer less than or equal to the specified number.
    /// For example: Floor(3.7) = 3; Floor(-3.1) = -4
    function Floor(value : Double) : Int {
        let (truncated, remainder, isPositive) = ExtendedTruncation(value);
        if AbsD(remainder) <= 1e-15 {
            truncated
        } else {
            isPositive ? truncated | truncated - 1
        }
    }

    /// # Summary
    /// Returns the nearest integer to the specified number.
    /// For example: Floor(3.7) = 4; Floor(-3.7) = -4
    function Round(value : Double) : Int {
        let (truncated, remainder, isPositive) = ExtendedTruncation(value);
        if AbsD(remainder) <= 1e-15 {
            truncated
        } else {
            let abs = AbsD(remainder);
            truncated + (abs <= 0.5 ? 0 | (isPositive ? 1 | -1))
        }
    }

    //
    // Modular arithmetic
    //

    /// # Summary
    /// Divides one Integer value by another, returns the result and the remainder as a tuple.
    function DivRemI(dividend : Int, divisor : Int) : (Int, Int) {
        (dividend / divisor, dividend % divisor)
    }

    /// # Summary
    /// Divides one BigInteger value by another, returns the result and the remainder as a tuple.
    function DivRemL(dividend : BigInt, divisor : BigInt) : (BigInt, BigInt) {
        (dividend / divisor, dividend % divisor)
    }

    /// # Summary
    /// Computes the canonical residue of `value` modulo `modulus`.
    /// The result is always in the range 0..modulus-1 even for negative numbers.
    function ModulusI(value : Int, modulus : Int) : Int {
        Fact(modulus > 0, "`modulus` must be positive");
        let r = value % modulus;
        (r < 0) ? (r + modulus) | r
    }

    /// # Summary
    /// Computes the canonical residue of `value` modulo `modulus`.
    /// The result is always in the range 0..modulus-1 even for negative numbers.
    function ModulusL(value : BigInt, modulus : BigInt) : BigInt {
        Fact(modulus > 0L, "`modulus` must be positive");
        let r = value % modulus;
        (r < 0L) ? (r + modulus) | r
    }

    /// # Summary
    /// Returns an integer raised to a given power, with respect to a given
    /// modulus. I.e. (expBase^power) % modulus.
    function ExpModI(expBase : Int, power : Int, modulus : Int) : Int {
        Fact(power >= 0, "`power` must be non-negative");
        Fact(modulus > 0, "`modulus` must be positive");
        Fact(expBase > 0, "`expBase` must be positive");

        // shortcut when modulus is 1
        if modulus == 1 {
            return 0;
        }

        mutable res = 1;
        mutable expPow2mod = expBase % modulus;
        mutable powerBits = power;

        while powerBits > 0 {
            if (powerBits &&& 1) != 0 {
                // if bit pₖ is 1, multiply res by expBase^(2^k) (mod `modulus`)
                set res = (res * expPow2mod) % modulus;
            }

            // update value of expBase^(2^k) (mod `modulus`)
            set expPow2mod = (expPow2mod * expPow2mod) % modulus;
            set powerBits >>>= 1;
        }

        res
    }

    /// # Summary
    /// Returns an integer raised to a given power, with respect to a given
    /// modulus. I.e. (expBase^power) % modulus.
    function ExpModL(expBase : BigInt, power : BigInt, modulus : BigInt) : BigInt {
        Fact(power >= 0L, "`power` must be non-negative");
        Fact(modulus > 0L, "`modulus` must be positive");
        Fact(expBase > 0L, "`expBase` must be positive");

        // shortcut when modulus is 1
        if modulus == 1L {
            return 0L;
        }

        mutable res = 1L;
        mutable expPow2mod = expBase % modulus;
        mutable powerBits = power;

        while powerBits > 0L {
            if (powerBits &&& 1L) != 0L {
                // if bit pₖ is 1, multiply res by expBase^(2ᵏ) (mod `modulus`)
                set res = (res * expPow2mod) % modulus;
            }

            // update value of expBase^(2ᵏ) (mod `modulus`)
            set expPow2mod = (expPow2mod * expPow2mod) % modulus;
            set powerBits >>>= 1;
        }

        res
    }

    /// # Summary
    /// Returns the multiplicative inverse of a modular integer.
    ///
    /// # Description
    /// This will calculate the multiplicative inverse of a
    /// modular integer `b` such that `a • b = 1 (mod modulus)`.
    function InverseModI(a : Int, modulus : Int) : Int {
        let (u, v) = ExtendedGreatestCommonDivisorI(a, modulus);
        let gcd = u * a + v * modulus;
        Fact(gcd == 1, "`a` and `modulus` must be co-prime");
        ModulusI(u, modulus)
    }

    /// # Summary
    /// Returns the multiplicative inverse of a modular integer.
    ///
    /// # Description
    /// This will calculate the multiplicative inverse of a
    /// modular integer `b` such that `a • b = 1 (mod modulus)`.
    function InverseModL(a : BigInt, modulus : BigInt) : BigInt {
        let (u, v) = ExtendedGreatestCommonDivisorL(a, modulus);
        let gcd = u * a + v * modulus;
        Fact(gcd == 1L, "`a` and `modulus` must be co-prime");
        ModulusL(u, modulus)
    }

    //
    // GCD, etc.
    //

    /// # Summary
    /// Computes the greatest common divisor of two integers.
    /// Note: GCD is always positive except that GCD(0,0)=0.
    function GreatestCommonDivisorI(a : Int, b : Int) : Int {
        mutable aa = AbsI(a);
        mutable bb = AbsI(b);
        while bb != 0 {
            let cc = aa % bb;
            set aa = bb;
            set bb = cc;
        }
        aa
    }

    /// # Summary
    /// Computes the greatest common divisor of two integers.
    /// Note: GCD is always positive except that GCD(0,0)=0.
    function GreatestCommonDivisorL(a : BigInt, b : BigInt) : BigInt {
        mutable aa = AbsL(a);
        mutable bb = AbsL(b);
        while bb != 0L {
            let cc = aa % bb;
            set aa = bb;
            set bb = cc;
        }
        aa
    }

    /// # Summary
    /// Returns a tuple (u,v) such that u*a+v*b=GCD(a,b)
    /// Note: GCD is always positive except that GCD(0,0)=0.
    function ExtendedGreatestCommonDivisorI(a : Int, b : Int) : (Int, Int) {
        let signA = SignI(a);
        let signB = SignI(b);
        mutable (s1, s2) = (1, 0);
        mutable (t1, t2) = (0, 1);
        mutable (r1, r2) = (a * signA, b * signB);

        while r2 != 0 {
            let quotient = r1 / r2;
            set (r1, r2) = (r2, r1 - quotient * r2);
            set (s1, s2) = (s2, s1 - quotient * s2);
            set (t1, t2) = (t2, t1 - quotient * t2);
        }

        (s1 * signA, t1 * signB)
    }

    /// # Summary
    /// Returns a tuple (u,v) such that u*a+v*b=GCD(a,b)
    /// Note: GCD is always positive except that GCD(0,0)=0.
    function ExtendedGreatestCommonDivisorL(a : BigInt, b : BigInt) : (BigInt, BigInt) {
        let signA = IntAsBigInt(SignL(a));
        let signB = IntAsBigInt(SignL(b));
        mutable (s1, s2) = (1L, 0L);
        mutable (t1, t2) = (0L, 1L);
        mutable (r1, r2) = (a * signA, b * signB);

        while r2 != 0L {
            let quotient = r1 / r2;
            set (r1, r2) = (r2, r1 - quotient * r2);
            set (s1, s2) = (s2, s1 - quotient * s2);
            set (t1, t2) = (t2, t1 - quotient * t2);
        }

        (s1 * signA, t1 * signB)
    }

    /// # Summary
    /// Returns if two integers are co-prime.
    ///
    /// # Description
    /// Returns true if a and b are co-prime and false otherwise.
    ///
    /// # Input
    /// ## a
    /// the first number of which co-primality is being tested
    /// ## b
    /// the second number of which co-primality is being tested
    ///
    /// # Output
    /// True, if a and b are co-prime (e.g. their greatest common divisor is 1),
    /// and false otherwise
    function IsCoprimeI(a : Int, b : Int) : Bool {
        GreatestCommonDivisorI(a, b) == 1
    }

    /// # Summary
    /// Returns if two integers are co-prime.
    ///
    /// # Description
    /// Returns true if a and b are co-prime and false otherwise.
    ///
    /// # Input
    /// ## a
    /// the first number of which co-primality is being tested
    /// ## b
    /// the second number of which co-primality is being tested
    ///
    /// # Output
    /// True, if a and b are co-prime (e.g. their greatest common divisor is 1),
    /// and false otherwise
    function IsCoprimeL(a : BigInt, b : BigInt) : Bool {
        GreatestCommonDivisorL(a, b) == 1L
    }

    /// # Summary
    /// Finds the continued fraction convergent closest to `fraction`
    /// with the denominator less or equal to `denominatorBound`
    /// Using process similar to this: https://nrich.maths.org/1397
    function ContinuedFractionConvergentI(
        fraction : (Int, Int),
        denominatorBound : Int
    ) : (Int, Int) {
        Fact(denominatorBound > 0, "Denominator bound must be positive");

        let (a, b) = fraction;
        let signA = SignI(a);
        let signB = SignI(b);
        mutable (s1, s2) = (1, 0);
        mutable (t1, t2) = (0, 1);
        mutable (r1, r2) = (a * signA, b * signB);

        while r2 != 0 and AbsI(s2) <= denominatorBound {
            let quotient = r1 / r2;
            set (r1, r2) = (r2, r1 - quotient * r2);
            set (s1, s2) = (s2, s1 - quotient * s2);
            set (t1, t2) = (t2, t1 - quotient * t2);
        }

        if r2 == 0 and AbsI(s2) <= denominatorBound {
            (-t2 * signB, s2 * signA)
        } else {
            (-t1 * signB, s1 * signA)
        }
    }

    /// # Summary
    /// Finds the continued fraction convergent closest to `fraction`
    /// with the denominator less or equal to `denominatorBound`
    /// Using process similar to this: https://nrich.maths.org/1397
    function ContinuedFractionConvergentL(
        fraction : (BigInt, BigInt),
        denominatorBound : BigInt
    ) : (BigInt, BigInt) {
        Fact(denominatorBound > 0L, "Denominator bound must be positive");

        let (a, b) = fraction;
        let signA = IntAsBigInt(SignL(a));
        let signB = IntAsBigInt(SignL(b));
        mutable (s1, s2) = (1L, 0L);
        mutable (t1, t2) = (0L, 1L);
        mutable (r1, r2) = (a * signA, b * signB);

        while r2 != 0L and AbsL(s2) <= denominatorBound {
            let quotient = r1 / r2;
            set (r1, r2) = (r2, r1 - quotient * r2);
            set (s1, s2) = (s2, s1 - quotient * s2);
            set (t1, t2) = (t2, t1 - quotient * t2);
        }

        if r2 == 0L and AbsL(s2) <= denominatorBound {
            (-t2 * signB, s2 * signA)
        } else {
            (-t1 * signB, s1 * signA)
        }
    }

    /// # Summary
    /// Computes the modulus between two real numbers.
    ///
    /// # Input
    /// ## value
    /// A real number x to take the modulus of.
    /// ## modulo
    /// A real number to take the modulus of x with respect to.
    /// ## minValue
    /// The smallest value to be returned by this function.
    ///
    /// # Example
    /// ```qsharp
    ///     // Returns 3 π / 2.
    ///     let y = RealMod(5.5 * PI(), 2.0 * PI(), 0.0);
    ///     // Returns -1.2, since +3.6 and -1.2 are 4.8 apart on the real line,
    ///     // which is a multiple of 2.4.
    ///     let z = RealMod(3.6, 2.4, -1.2);
    /// ```
    function RealMod(value : Double, modulo : Double, minValue : Double) : Double {
        let timesModuloInSegment = (value - minValue) / modulo;
        let fractionalPart = timesModuloInSegment - IntAsDouble(Truncate(timesModuloInSegment));
        modulo * fractionalPart + minValue
    }

    //
    // Binary, bits, etc.
    //

    /// # Summary
    /// For a non-negative integer `a`, returns the number of bits required to represent `a`.
    /// NOTE: This function returns the smallest n such that a < 2^n.
    function BitSizeI(a : Int) : Int {
        Fact(a >= 0, "`a` must be non-negative.");
        mutable number = a;
        mutable size = 0;
        while (number != 0) {
            set size = size + 1;
            set number = number >>> 1;
        }

        size
    }

    /// # Summary
    /// For a non-negative integer `a`, returns the number of bits required to represent `a`.
    /// NOTE: This function returns the smallest n such that a < 2^n.
    function BitSizeL(a : BigInt) : Int {
        Fact(a >= 0L, "`a` must be non-negative.");
        mutable number = a;
        mutable size = 0;
        while (number != 0L) {
            set size = size + 1;
            set number = number >>> 1;
        }

        size
    }

    /// # Summary
    /// For a non-zero integer `a`, returns the number of trailing zero bits
    /// in the binary representation of `a`.
    function TrailingZeroCountI(a : Int) : Int {
        Fact(a != 0, "TrailingZeroCountI: `a` cannot be 0.");

        mutable count = 0;
        mutable n = a;
        while n &&& 1 == 0 {
            set count += 1;
            set n >>>= 1;
        }

        count
    }

    /// # Summary
    /// For a non-zero integer `a`, returns the number of trailing zero bits
    /// in the binary representation of `a`.
    function TrailingZeroCountL(a : BigInt) : Int {
        Fact(a != 0L, "TrailingZeroCountL: `a` cannot be 0.");

        mutable count = 0;
        mutable n = a;
        while n &&& 1L == 0L {
            set count += 1;
            set n >>>= 1;
        }

        count
    }

    /// # Summary
    /// Returns the number of 1 bits in the binary representation of integer `n`.
    function HammingWeightI(n : Int) : Int {
        let i1 = n - ((n >>> 1) &&& 0x5555555555555555);
        let i2 = (i1 &&& 0x3333333333333333) + ((i1 >>> 2) &&& 0x3333333333333333);
        // Multiplication may overflow. See https://github.com/microsoft/qsharp/issues/828
        (((i2 + (i2 >>> 4)) &&& 0xF0F0F0F0F0F0F0F) * 0x101010101010101) >>> 56
    }

    //
    // Combinatorics
    //

    /// # Summary
    /// Returns the factorial of a given number.
    ///
    /// # Description
    /// Returns the factorial of a given nonnegative integer n, where 0 ≤ n ≤ 20.
    ///
    /// # Input
    /// ## n
    /// The number to take the factorial of.
    ///
    /// # Output
    /// The factorial of `n`.
    ///
    /// # Remarks
    /// For inputs greater than 20, please use `Microsoft.Quantum.Math.FactorialL`.
    ///
    /// # See Also
    /// - Microsoft.Quantum.Math.FactorialL
    /// - Microsoft.Quantum.Math.ApproximateFactorial
    function FactorialI(n : Int) : Int {
        Fact(n >= 0, "The factorial is not defined for negative inputs.");
        Fact(n <= 20, "The largest factorial that can be stored as an Int is 20!. Use FactorialL or ApproximateFactorial.");

        [
            1,
            1,
            2,
            6,
            24,
            120,
            720,
            5040,
            40320,
            362880,
            3628800,
            39916800,
            479001600,
            6227020800,
            87178291200,
            1307674368000,
            20922789888000,
            355687428096000,
            6402373705728000,
            121645100408832000,
            2432902008176640000
        ][n]
    }

    /// # Summary
    /// Returns the factorial of a given number.
    ///
    /// # Input
    /// ## n
    /// The number to take the factorial of.
    ///
    /// # Output
    /// The factorial of `n`.
    ///
    /// # See Also
    /// - Microsoft.Quantum.Math.FactorialI
    /// - Microsoft.Quantum.Math.ApproximateFactorial
    function FactorialL(n : Int) : BigInt {
        Fact(n >= 0, "The factorial is not defined for negative inputs.");

        mutable result = 1L;
        for i in 1..n {
            set result *= IntAsBigInt(i);
        }
        result
    }

    /// # Summary
    /// Returns an approximate factorial of a given number.
    ///
    /// # Description
    /// Returns the factorial as `Double`, given an input `n`.
    /// The domain of inputs for this function is `n <= 169`.
    ///
    /// # Remarks
    /// For n > 10, this function uses the Ramanujan approximation with a
    /// relative error of the order of 1 / n⁵.
    ///
    /// # Input
    /// ## n
    /// The number to take the approximate factorial of. Must not be negative.
    ///
    /// # Output
    /// The approximate factorial of `n`.
    ///
    /// # See Also
    /// - Microsoft.Quantum.Math.FactorialI
    /// - Microsoft.Quantum.Math.FactorialL
    function ApproximateFactorial(n : Int) : Double {
        Fact(n >= 0, "The factorial is not defined for negative inputs.");
        Fact(n <= 169, "The largest approximate factorial that can be stored as a Double is 169!. Use FactorialL.");

        // For small enough n, use the exact factorial instead.
        if n <= 20 {
            return IntAsDouble(FactorialI(n));
        }

        let absN = IntAsDouble(n);
        let a = Sqrt(2.0 * PI() * absN);
        let b = (absN / E())^absN;
        let c = E()^(1.0 / (12.0 * absN) - (1.0 / (360.0 * (absN^3.0))));

        a * b * c
    }

    /// # Summary
    /// Returns the natural logarithm of the gamma function (aka the log-gamma
    /// function).
    ///
    /// # Description
    /// The gamma function Γ(x) generalizes the factorial function
    /// to the positive real numbers and is defined as
    /// integral from 0 to ∞ of t¹⁻ˣ⋅e⁻ᵗ𝑑t
    ///
    /// The gamma function has the property that for all positive real numbers
    /// x, Γ(x + 1) = x⋅Γ(x), such that the factorial function
    /// is a special case of Γ, n! = Γ(n + 1) for all natural numbers n.
    ///
    /// # Input
    /// ## x
    /// The point x at which the log-gamma function is to be evaluated.
    ///
    /// # Output
    /// The value ㏑(Γ(x)).
    function LogGammaD(x : Double) : Double {
        // Here, we use the approximation described in Numerical Recipes in C.
        let coefficients = [
            57.1562356658629235,
            -59.5979603554754912,
            14.1360979747417471,
            -0.491913816097620199,
            0.339946499848118887e-4,
            0.465236289270485756e-4,
            -0.983744753048795646e-4,
            0.158088703224912494e-3,
            -0.210264441724104883e-3,
            0.217439618115212643e-3,
            -0.164318106536763890e-3,
            0.844182239838527433e-4,
            -0.261908384015814087e-4,
            0.368991826595316234e-5
        ];

        Fact(x > 0.0, "Γ(x) not defined for x <= 0.");

        mutable y = x;
        let tmp = x + 5.2421875000000000;

        mutable acc = 0.99999999999999709;
        for coeff in coefficients {
            set y += 1.0;
            set acc += coeff / y;
        }

        Log(2.506628274631000 * acc / x) + ((x + 0.5) * Log(tmp) - tmp)
    }

    /// # Summary
    /// Returns the approximate natural logarithm of the factorial of a given
    /// integer.
    ///
    /// # Input
    /// ## n
    /// The number to take the log-factorial of.
    ///
    /// # Output
    /// The natural logarithm of the factorial of the provided input.
    ///
    /// # See Also
    /// - Microsoft.Quantum.Math.ApproximateFactorial
    /// - Microsoft.Quantum.Math.FactorialI
    /// - Microsoft.Quantum.Math.FactorialL
    function LogFactorialD(n : Int) : Double {
        LogGammaD(IntAsDouble(n) + 1.0)
    }

    /// # Summary
    /// Returns the approximate binomial coefficient of two integers.
    ///
    /// # Description
    /// Given two integers n and k, returns the binomial coefficient
    /// binom(n, k), also known as n-choose-k. Computed approximately.
    ///
    /// # Input
    /// ## n
    /// The first of the two integers to compute the binomial coefficient of.
    /// ## k
    /// The second of the two integers to compute the binomial coefficient of.
    ///
    /// # Output
    /// The binomial coefficient n-choose-k.
    function Binom(n : Int, k : Int) : Int {
        // Here, we use the approximation described in Numerical Recipes in C.
        if n < 171 {
            Floor(0.5 + ApproximateFactorial(n) / (ApproximateFactorial(k) * ApproximateFactorial(n - k)))
        } else {
            Floor(0.5 + E()^(LogFactorialD(n) - LogFactorialD(k) - LogFactorialD(n - k)))
        }
    }

    //
    // Norms
    //

    /// # Summary
    /// Returns the squared 2-norm of a vector.
    ///
    /// # Description
    /// Returns the squared 2-norm of a vector; that is, given an input
    /// x̄, returns ∑xᵢ.
    ///
    /// # Input
    /// ## array
    /// The vector whose squared 2-norm is to be returned.
    ///
    /// # Output
    /// The squared 2-norm of `array`.
    function SquaredNorm(array : Double[]) : Double {
        mutable sum = 0.0;
        for element in array {
            set sum += element * element;
        }

        sum
    }

    /// # Summary
    /// Returns the `L(p)` norm of a vector of `Double`s.
    ///
    /// That is, given an array x of type `Double[]`, this returns the p-norm
    /// |x̄|ₚ= (∑(xᵢ)ᵖ)¹ᐟᵖ.
    ///
    /// # Input
    /// ## p
    /// The exponent p in the p-norm.
    ///
    /// # Output
    /// The p-norm |x̄|ₚ.
    function PNorm(p : Double, array : Double[]) : Double {
        if p < 1.0 {
            fail "p must be >= 1.0";
        }

        mutable sum = 0.0;
        for element in array {
            set sum += AbsD(element)^p;
        }

        sum^(1.0 / p)
    }

    /// # Summary
    /// Normalizes a vector of `Double`s in the `L(p)` norm.
    ///
    /// That is, given an array x of type `Double[]`, this returns an array where
    /// all elements are divided by the p-norm |x̄|ₚ.
    /// Function leaves array with norm 0 unchanged.
    ///
    /// # Input
    /// ## p
    /// The exponent p in the p-norm.
    ///
    /// # Output
    /// The array x normalized by the p-norm |x̄|ₚ.
    ///
    /// # See Also
    /// - PNorm
    function PNormalized(p : Double, array : Double[]) : Double[] {
        let norm = PNorm(p, array);
        if (norm == 0.0) {
            return array;
        }

        mutable result = [];
        for element in array {
            set result += [element / norm];
        }

        result
    }

    //
    // Complex numbers
    //

    /// # Summary
    /// Represents a complex number by its real and imaginary components.
    /// The first element of the tuple is the real component,
    /// the second one - the imaginary component.
    ///
    /// # Example
    /// The following snippet defines the imaginary unit 𝑖 = 0 + 1𝑖:
    /// ```qsharp
    /// let imagUnit = Complex(0.0, 1.0);
    /// ```
    newtype Complex = (Real : Double, Imag : Double);

    /// # Summary
    /// Represents a complex number in polar form.
    /// The polar representation of a complex number is c = r⋅𝑒^(t𝑖).
    ///
    /// # Named Items
    /// ## Magnitude
    /// The absolute value r>0 of c.
    /// ## Argument
    /// The phase t ∈ ℝ of c.
    newtype ComplexPolar = (Magnitude : Double, Argument : Double);

    /// # Summary
    /// Returns the squared absolute value of a complex number of type
    /// `Complex`.
    ///
    /// # Input
    /// ## input
    /// Complex number c = x + y𝑖.
    ///
    /// # Output
    /// Squared absolute value |c|² = x² + y².
    function AbsSquaredComplex(input : Complex) : Double {
        input::Real * input::Real + input::Imag * input::Imag
    }

    /// # Summary
    /// Returns the absolute value of a complex number of type
    /// `Complex`.
    ///
    /// # Input
    /// ## input
    /// Complex number c = x + y𝑖.
    ///
    /// # Output
    /// Absolute value |c| = √(x² + y²).
    function AbsComplex(input : Complex) : Double {
        Sqrt(AbsSquaredComplex(input))
    }

    /// # Summary
    /// Returns the phase of a complex number of type
    /// `Complex`.
    ///
    /// # Input
    /// ## input
    /// Complex number c = x + y𝑖.
    ///
    /// # Output
    /// Phase Arg(c) = ArcTan(y,x) ∈ (-𝜋,𝜋].
    function ArgComplex(input : Complex) : Double {
        ArcTan2(input::Imag, input::Real)
    }

    /// # Summary
    /// Returns the squared absolute value of a complex number of type
    /// `ComplexPolar`.
    ///
    /// # Input
    /// ## input
    /// Complex number c = r⋅𝑒^(t𝑖).
    ///
    /// # Output
    /// Squared absolute value |c|² = r².
    function AbsSquaredComplexPolar(input : ComplexPolar) : Double {
        input::Magnitude * input::Magnitude
    }

    /// # Summary
    /// Returns the absolute value of a complex number of type
    /// `ComplexPolar`.
    ///
    /// # Input
    /// ## input
    /// Complex number c = r⋅𝑒^(t𝑖).
    ///
    /// # Output
    /// Absolute value |c| = r.
    function AbsComplexPolar(input : ComplexPolar) : Double {
        input::Magnitude
    }

    /// # Summary
    /// Returns the phase of a complex number of type `ComplexPolar`.
    ///
    /// # Input
    /// ## input
    /// Complex number c = r⋅𝑒^(t𝑖).
    ///
    /// # Output
    /// Phase Arg(c) = t.
    function ArgComplexPolar(input : ComplexPolar) : Double {
        input::Argument
    }

    /// # Summary
    /// Returns the unary negation of an input of type `Complex`.
    ///
    /// # Input
    /// ## input
    /// A value whose negation is to be returned.
    ///
    /// # Output
    /// The unary negation of `input`.
    function NegationC(input : Complex) : Complex {
        Complex(-input::Real, -input::Imag)
    }

    /// # Summary
    /// Returns the unary negation of an input of type `ComplexPolar`
    ///
    /// # Input
    /// ## input
    /// A value whose negation is to be returned.
    ///
    /// # Output
    /// The unary negation of `input`.
    function NegationCP(input : ComplexPolar) : ComplexPolar {
        ComplexPolar(input::Magnitude, input::Argument + PI())
    }

    /// # Summary
    /// Returns the sum of two inputs of type `Complex`.
    ///
    /// # Input
    /// ## a
    /// The first input a to be summed.
    /// ## b
    /// The second input b to be summed.
    ///
    /// # Output
    /// The sum a + b.
    function PlusC(a : Complex, b : Complex) : Complex {
        Complex(a::Real + b::Real, a::Imag + b::Imag)
    }

    /// # Summary
    /// Returns the sum of two inputs of type `ComplexPolar`.
    ///
    /// # Input
    /// ## a
    /// The first input a to be summed.
    /// ## b
    /// The second input b to be summed.
    ///
    /// # Output
    /// The sum a + b.
    function PlusCP(a : ComplexPolar, b : ComplexPolar) : ComplexPolar {
        ComplexAsComplexPolar(
            PlusC(
                ComplexPolarAsComplex(a),
                ComplexPolarAsComplex(b)
            )
        )
    }

    /// # Summary
    /// Returns the difference between two inputs of type `Complex`.
    ///
    /// # Input
    /// ## a
    /// The first input a to be subtracted.
    /// ## b
    /// The second input b to be subtracted.
    ///
    /// # Output
    /// The difference a - b.
    function MinusC(a : Complex, b : Complex) : Complex {
        Complex(a::Real - b::Real, a::Imag - b::Imag)
    }

    /// # Summary
    /// Returns the difference between two inputs of type `ComplexPolar`.
    ///
    /// # Input
    /// ## a
    /// The first input a to be subtracted.
    /// ## b
    /// The second input b to be subtracted.
    ///
    /// # Output
    /// The difference a - b.
    function MinusCP(a : ComplexPolar, b : ComplexPolar) : ComplexPolar {
        PlusCP(a, NegationCP(b))
    }

    /// # Summary
    /// Returns the product of two inputs of type `Complex`.
    ///
    /// # Input
    /// ## a
    /// The first input a to be multiplied.
    /// ## b
    /// The second input b to be multiplied.
    ///
    /// # Output
    /// The product a⋅b.
    function TimesC(a : Complex, b : Complex) : Complex {
        Complex(
            a::Real * b::Real - a::Imag * b::Imag,
            a::Real * b::Imag + a::Imag * b::Real
        )
    }

    /// # Summary
    /// Returns the product of two inputs of type `ComplexPolar`.
    ///
    /// # Input
    /// ## a
    /// The first input a to be multiplied.
    /// ## b
    /// The second input b to be multiplied.
    ///
    /// # Output
    /// The product a⋅b.
    function TimesCP(a : ComplexPolar, b : ComplexPolar) : ComplexPolar {
        ComplexPolar(
            a::Magnitude * b::Magnitude,
            a::Argument + b::Argument
        )
    }

    /// # Summary
    /// Internal. Since it is easiest to define the power of two complex numbers
    /// in Cartesian form as returning in polar form, we define that here, then
    /// convert as needed.
    /// Note that this is a multi-valued function, but only one value is returned.
    internal function PowCAsCP(base : Complex, power : Complex) : ComplexPolar {
        let ((a, b), (c, d)) = (base!, power!);
        let baseSqNorm = a * a + b * b;
        let baseNorm = Sqrt(baseSqNorm);
        let baseArg = ArgComplex(base);

        // We pick the principal value of the multi-valued complex function ㏑ as
        // ㏑(a+b𝑖) = ln(|a+b𝑖|) + 𝑖⋅arg(a+b𝑖) = ln(baseNorm) + 𝑖⋅baseArg
        // Therefore
        // base^power = (a+b𝑖)^(c+d𝑖) = 𝑒^( (c+d𝑖)⋅㏑(a+b𝑖) ) =
        // = 𝑒^( (c+d𝑖)⋅(ln(baseNorm)+𝑖⋅baseArg) ) =
        // = 𝑒^( (c⋅ln(baseNorm) - d⋅baseArg) + 𝑖⋅(c⋅baseArg + d⋅ln(baseNorm)) )
        // magnitude = 𝑒^((c⋅ln(baseNorm) - d⋅baseArg)) = baseNorm^c / 𝑒^(d⋅baseArg)
        // angle = d⋅ln(baseNorm) + c⋅baseArg

        let magnitude = baseNorm^c / E()^(d * baseArg);
        let angle = d * Log(baseNorm) + c * baseArg;

        ComplexPolar(magnitude, angle)
    }

    /// # Summary
    /// Returns a number raised to a given power of type `Complex`.
    /// Note that this is a multi-valued function, but only one value is returned.
    ///
    /// # Input
    /// ## a
    /// The number a that is to be raised.
    /// ## power
    /// The power b to which a should be raised.
    ///
    /// # Output
    /// The power a^b
    function PowC(a : Complex, power : Complex) : Complex {
        ComplexPolarAsComplex(PowCAsCP(a, power))
    }

    /// # Summary
    /// Returns a number raised to a given power of type `ComplexPolar`.
    /// Note that this is a multi-valued function, but only one value is returned.
    ///
    /// # Input
    /// ## a
    /// The number a that is to be raised.
    /// ## power
    /// The power b to which a should be raised.
    ///
    /// # Output
    /// The power a^b
    function PowCP(a : ComplexPolar, power : ComplexPolar) : ComplexPolar {
        PowCAsCP(ComplexPolarAsComplex(a), ComplexPolarAsComplex(power))
    }

    /// # Summary
    /// Returns the quotient of two inputs of type `Complex`.
    ///
    /// # Input
    /// ## a
    /// The first input a to be divided.
    /// ## b
    /// The second input b to be divided.
    ///
    /// # Output
    /// The quotient a / b.
    function DividedByC(a : Complex, b : Complex) : Complex {
        let sqNorm = b::Real * b::Real + b::Imag * b::Imag;
        Complex(
            (a::Real * b::Real + a::Imag * b::Imag) / sqNorm,
            (a::Imag * b::Real - a::Real * b::Imag) / sqNorm
        )
    }

    /// # Summary
    /// Returns the quotient of two inputs of type `ComplexPolar`.
    ///
    /// # Input
    /// ## a
    /// The first input a to be divided.
    /// ## b
    /// The second input b to be divided.
    ///
    /// # Output
    /// The quotient a / b.
    function DividedByCP(a : ComplexPolar, b : ComplexPolar) : ComplexPolar {
        ComplexPolar(a::Magnitude / b::Magnitude, a::Argument - b::Argument)
    }

    //
    // Fixed point
    //

    /// # Summary
    /// Returns the smallest representable number for specific fixed point dimensions.
    ///
    /// # Input
    /// ## integerBits
    /// Number of integer bits (including the sign bit).
    /// ## fractionalBits
    /// Number of fractional bits.
    ///
    /// # Remark
    /// The value can be computed as -2^(p-1), where p is the number of integer bits.
    function SmallestFixedPoint(integerBits : Int, fractionalBits : Int) : Double {
        -(2.0^IntAsDouble(integerBits - 1))
    }

    /// # Summary
    /// Returns the largest representable number for specific fixed point dimensions.
    ///
    /// # Input
    /// ## integerBits
    /// Number of integer bits (including the sign bit).
    /// ## fractionalBits
    /// Number of fractional bits.
    ///
    /// # Remark
    /// The value can be computed as 2^(p-1) - 2^(-q), where p
    /// is the number of integer bits and q is the number of fractional bits.
    function LargestFixedPoint(integerBits : Int, fractionalBits : Int) : Double {
        2.0^IntAsDouble(integerBits - 1) - 2.0^(-IntAsDouble(fractionalBits))
    }

}
