// Licensed to the .NET Foundation under one or more agreements.
// The .NET Foundation licenses this file to you under the Apache 2.0 License.
// See the LICENSE file in the project root for more information.

using System;
using System.Collections;
using System.Collections.Generic;
using System.Diagnostics;
using System.Linq;
using System.Numerics;

using Microsoft.Scripting;
using Microsoft.Scripting.Utils;

using IronPython.Runtime;
using IronPython.Runtime.Operations;
using IronPython.Runtime.Types;

[assembly: PythonModule("math", typeof(IronPython.Modules.PythonMath))]
namespace IronPython.Modules {
    public static partial class PythonMath {
        public const string __doc__ = "Provides common mathematical functions.";

#if NET
        public const double tau = Math.Tau;
#else
        public const double tau = 2 * Math.PI;
#endif
        public const double pi = Math.PI;
        public const double e = Math.E;

        private const double degreesToRadians = Math.PI / 180.0;
        private const double radiansToDegrees = 180.0 / Math.PI;

        private const int Bias = 0x3FE;

        public static double degrees(double radians) {
            return Check(radians, radians * radiansToDegrees);
        }

        public static double radians(double degrees) {
            return Check(degrees, degrees * degreesToRadians);
        }

        public static double fmod(double v, double w) {
            return Check(v, w, v % w);
        }

        private static double sum(List<double> partials) {
            // sum the partials the same was as CPython does
            var n = partials.Count;
            var hi = 0.0;

            if (n == 0) return hi;

            var lo = 0.0;

            // sum exact
            while (n > 0) {
                var x = hi;
                var y = partials[--n];
                hi = x + y;
                lo = y - (hi - x);
                if (lo != 0.0)
                    break;
            }

            if (n == 0) return hi;

            // half-even rounding
            if (lo < 0.0 && partials[n - 1] < 0.0 || lo > 0.0 && partials[n - 1] > 0.0) {
                var y = lo * 2.0;
                var x = hi + y;
                var yr = x - hi;
                if (y == yr)
                    hi = x;
            }
            return hi;
        }

        public static double fsum(IEnumerable e) {
            // msum from https://code.activestate.com/recipes/393090/
            var partials = new List<double>();
            foreach (var v in e.Cast<object>().Select(o => Converter.ConvertToDouble(o))) {
                var x = v;
                var i = 0;
                for (var j = 0; j < partials.Count; j++) {
                    var y = partials[j];
                    if (Math.Abs(x) < Math.Abs(y)) {
                        var t = x;
                        x = y;
                        y = t;
                    }
                    var hi = x + y;
                    var lo = y - (hi - x);
                    if (lo != 0) {
                        partials[i++] = lo;
                    }
                    x = hi;
                }
                partials.RemoveRange(i, partials.Count - i);
                partials.Add(x);
            }

            return sum(partials);
        }

        public static PythonTuple frexp(double v) {
            if (Double.IsInfinity(v) || Double.IsNaN(v)) {
                return PythonTuple.MakeTuple(v, 0.0);
            }
            int exponent = 0;
            double mantissa = 0;

            if (v == 0) {
                mantissa = 0;
                exponent = 0;
            } else {
                byte[] vb = BitConverter.GetBytes(v);
                if (BitConverter.IsLittleEndian) {
                    DecomposeLe(vb, out mantissa, out exponent);
                } else {
                    throw new NotImplementedException();
                }
            }

            return PythonTuple.MakeTuple(mantissa, exponent);
        }

        public static PythonTuple modf(double v) {
            if (double.IsInfinity(v)) {
                return PythonTuple.MakeTuple(0.0, v);
            }
            double w = v % 1.0;
            v -= w;
            return PythonTuple.MakeTuple(w, v);
        }

        public static double ldexp(double v, BigInteger w) {
            if (v == 0.0 || double.IsInfinity(v)) {
                return v;
            }
            return Check(v, v * Math.Pow(2.0, (double)w));
        }

        public static double hypot(double v, double w) {
            if (double.IsInfinity(v) || double.IsInfinity(w)) {
                return double.PositiveInfinity;
            }
            return Check(v, w, MathUtils.Hypot(v, w));
        }

        public static double pow(double v, double exp) {
            if (v == 1.0 || exp == 0.0) {
                return 1.0;
            } else if (double.IsNaN(v) || double.IsNaN(exp)) {
                return double.NaN;
            } else if (v == 0.0) {
                if (exp > 0.0) {
                    return 0.0;
                }
                throw PythonOps.ValueError("math domain error");
            } else if (double.IsPositiveInfinity(exp)) {
                if (v > 1.0 || v < -1.0) {
                    return double.PositiveInfinity;
                } else if (v == -1.0) {
                    return 1.0;
                } else {
                    return 0.0;
                }
            } else if (double.IsNegativeInfinity(exp)) {
                if (v > 1.0 || v < -1.0) {
                    return 0.0;
                } else if (v == -1.0) {
                    return 1.0;
                } else {
                    return double.PositiveInfinity;
                }
            }
            return Check(v, exp, Math.Pow(v, exp));
        }

        public static double log(double v0) {
            if (v0 <= 0.0) {
                throw PythonOps.ValueError("math domain error");
            }
            return Check(v0, Math.Log(v0));
        }

        public static double log(double v0, double v1) {
            if (v0 <= 0.0 || v1 == 0.0) {
                throw PythonOps.ValueError("math domain error");
            } else if (v1 == 1.0) {
                throw PythonOps.ZeroDivisionError("float division");
            } else if (v1 == Double.PositiveInfinity) {
                return 0.0;
            }
            return Check(Math.Log(v0, v1));
        }

        public static double log(BigInteger value) {
            if (value.Sign <= 0) {
                throw PythonOps.ValueError("math domain error");
            }
            return value.Log();
        }

        public static double log(object value) {
            // CPython tries float first, then double, so we need
            // an explicit overload which properly matches the order here
            double val;
            if (Converter.TryConvertToDouble(value, out val)) {
                return log(val);
            } else {
                return log(Converter.ConvertToBigInteger(value));
            }
        }

        public static double log(BigInteger value, double newBase) {
            if (newBase <= 0.0 || value <= 0) {
                throw PythonOps.ValueError("math domain error");
            } else if (newBase == 1.0) {
                throw PythonOps.ZeroDivisionError("float division");
            } else if (newBase == Double.PositiveInfinity) {
                return 0.0;
            }
            return Check(value.Log(newBase));
        }

        public static double log(object value, double newBase) {
            // CPython tries float first, then double, so we need
            // an explicit overload which properly matches the order here
            double val;
            if (Converter.TryConvertToDouble(value, out val)) {
                return log(val, newBase);
            } else {
                return log(Converter.ConvertToBigInteger(value), newBase);
            }
        }

        public static double log2(double x) {
            if (x <= 0) throw PythonOps.ValueError("math domain error");
            if (double.IsPositiveInfinity(x) || double.IsNaN(x)) return x;

            if (!BitConverter.IsLittleEndian) return Math.Log(x, 2);

            int exponent = 0;
            double mantissa = 0;

            byte[] vb = BitConverter.GetBytes(x);
            DecomposeLe(vb, out mantissa, out exponent);

            if (x >= 1)
                return Math.Log(mantissa * 2, 2) + (exponent - 1); // similar to CPython for precision
            else
                return Math.Log(mantissa, 2) + exponent;
        }

        public static double log2(BigInteger x) {
            if (x <= 0) throw PythonOps.ValueError("math domain error");

            // cast to double if we can
            var d = (double)x;
            if (!double.IsPositiveInfinity(d)) return log2(d);

            // bring to into double range and try again
            var y = BigInteger.Log(x, 2);
            var z = (int)Math.Ceiling(y) - 1023;
            x >>= z;

            Debug.Assert(!double.IsPositiveInfinity((double)x));

            return log2((double)x) + z;
        }

        public static double log10(double v0) {
            if (v0 <= 0.0) {
                throw PythonOps.ValueError("math domain error");
            }
            return Check(v0, Math.Log10(v0));
        }

        public static double log10(BigInteger value) {
            if (value.Sign <= 0) {
                throw PythonOps.ValueError("math domain error");
            }
            return value.Log10();
        }

        public static double log10(object value) {
            // CPython tries float first, then double, so we need
            // an explicit overload which properly matches the order here
            double val;
            if (Converter.TryConvertToDouble(value, out val)) {
                return log10(val);
            } else {
                return log10(Converter.ConvertToBigInteger(value));
            }
        }

        public static double log1p(double v0) {
            // Calculate log(1.0 + v0) using William Kahan's algorithm for numerical precision

            if (double.IsPositiveInfinity(v0)) {
                return double.PositiveInfinity;
            }

            double v1 = v0 + 1.0;

            // Linear approximation for very small v0
            if (v1 == 1.0) {
                return v0;
            }

            // Apply correction factor
            return log(v1) * (v0 / (v1 - 1.0));
        }

        public static double log1p(BigInteger value) {
            return log(value + BigInteger.One);
        }

        public static double log1p(object value) {
            // CPython tries float first, then double, so we need
            // an explicit overload which properly matches the order here
            double val;
            if (Converter.TryConvertToDouble(value, out val)) {
                return log1p(val);
            } else {
                return log1p(Converter.ConvertToBigInteger(value));
            }
        }

        public static double expm1(double v0) {
            return Check(v0, Math.Tanh(v0 / 2.0) * (Math.Exp(v0) + 1.0));
        }

        public static double asinh(double v0) {
            if (v0 == 0.0 || double.IsInfinity(v0)) {
                return v0;
            }
            // rewrote ln(v0 + sqrt(v0**2 + 1)) for precision
            if (Math.Abs(v0) > 1.0) {
                return Math.Sign(v0) * (Math.Log(Math.Abs(v0)) + Math.Log(1.0 + MathUtils.Hypot(1.0, 1.0 / v0)));
            } else {
                return Math.Log(v0 + MathUtils.Hypot(1.0, v0));
            }
        }

        public static double asinh(object value) {
            // CPython tries float first, then double, so we need
            // an explicit overload which properly matches the order here
            double val;
            if (Converter.TryConvertToDouble(value, out val)) {
                return asinh(val);
            } else {
                return asinh(Converter.ConvertToBigInteger(value));
            }
        }

        public static double acosh(double v0) {
            if (v0 < 1.0) {
                throw PythonOps.ValueError("math domain error");
            } else if (double.IsPositiveInfinity(v0)) {
                return double.PositiveInfinity;
            }
            // rewrote ln(v0 + sqrt(v0**2 - 1)) for precision
            double c = Math.Sqrt(v0 + 1.0);
            return Math.Log(c) + Math.Log(v0 / c + Math.Sqrt(v0 - 1.0));
        }

        public static double acosh(object value) {
            // CPython tries float first, then double, so we need
            // an explicit overload which properly matches the order here
            double val;
            if (Converter.TryConvertToDouble(value, out val)) {
                return acosh(val);
            } else {
                return acosh(Converter.ConvertToBigInteger(value));
            }
        }

        public static double atanh(double v0) {
            if (v0 >= 1.0 || v0 <= -1.0) {
                throw PythonOps.ValueError("math domain error");
            } else if (v0 == 0.0) {
                // preserve +/-0.0
                return v0;
            }

            return Math.Log((1.0 + v0) / (1.0 - v0)) * 0.5;
        }

        public static double atanh(BigInteger value) {
            if (value == 0) {
                return 0;
            } else {
                throw PythonOps.ValueError("math domain error");
            }
        }

        public static double atanh(object value) {
            // CPython tries float first, then double, so we need
            // an explicit overload which properly matches the order here
            double val;
            if (Converter.TryConvertToDouble(value, out val)) {
                return atanh(val);
            } else {
                return atanh(Converter.ConvertToBigInteger(value));
            }
        }

        public static double atan2(double v0, double v1) {
            if (double.IsNaN(v0) || double.IsNaN(v1)) {
                return double.NaN;
            } else if (double.IsInfinity(v0)) {
                if (double.IsPositiveInfinity(v1)) {
                    return pi * 0.25 * Math.Sign(v0);
                } else if (double.IsNegativeInfinity(v1)) {
                    return pi * 0.75 * Math.Sign(v0);
                } else {
                    return pi * 0.5 * Math.Sign(v0);
                }
            } else if (double.IsInfinity(v1)) {
                return v1 > 0.0 ? 0.0 : pi * DoubleOps.Sign(v0);
            }
            return Math.Atan2(v0, v1);
        }

        public static object ceil(CodeContext context, object x) {
            object val;
            if (PythonTypeOps.TryInvokeUnaryOperator(context, x, "__ceil__", out val)) {
                return val;
            }

            throw PythonOps.TypeError("a float is required");
        }

        public static object ceil(double v0) {
            if (double.IsInfinity(v0)) throw PythonOps.OverflowError("cannot convert float infinity to integer");
            if (double.IsNaN(v0)) throw PythonOps.ValueError("cannot convert float NaN to integer");

            var res = Math.Ceiling(v0);
            if (res < int.MinValue || res > int.MaxValue) {
                return (BigInteger)res;
            }
            return (int)res;
        }

        /// <summary>
        /// Error function on real values
        /// </summary>
        public static double erf(double v0) {
            return MathUtils.Erf(v0);
        }

        /// <summary>
        /// Complementary error function on real values: erfc(x) =  1 - erf(x)
        /// </summary>
        public static double erfc(double v0) {
            return MathUtils.ErfComplement(v0);
        }

        public static object factorial(double v0) {
            if (v0 % 1.0 != 0.0) {
                throw PythonOps.ValueError("factorial() only accepts integral values");
            }
            return factorial((BigInteger)v0);
        }

        public static object factorial(BigInteger value) {
            if (value < 0) {
                throw PythonOps.ValueError("factorial() not defined for negative values");
            }
            if (value > SysModule.maxsize) {
                throw PythonOps.OverflowError("factorial() argument should not exceed {0}", SysModule.maxsize);
            }

            BigInteger val = 1;
            for (BigInteger mul = value; mul > BigInteger.One; mul -= BigInteger.One) {
                val *= mul;
            }

            if (val > int.MaxValue) {
                return val;
            }
            return (int)val;
        }

        public static object factorial(object value) {
            // CPython tries float first, then double, so we need
            // an explicit overload which properly matches the order here
            double val;
            if (Converter.TryConvertToDouble(value, out val)) {
                return factorial(val);
            } else {
                return factorial(Converter.ConvertToBigInteger(value));
            }
        }

        public static object floor(CodeContext context, object x) {
            object val;
            if (PythonTypeOps.TryInvokeUnaryOperator(context, x, "__floor__", out val)) {
                return val;
            }

            throw PythonOps.TypeError("a float is required");
        }

        public static object floor(double v0) {
            if (double.IsInfinity(v0)) throw PythonOps.OverflowError("cannot convert float infinity to integer");
            if (double.IsNaN(v0)) throw PythonOps.ValueError("cannot convert float NaN to integer");

            var res = Math.Floor(v0);
            if (res < int.MinValue || res > int.MaxValue) {
                return (BigInteger)res;
            }
            return (int)res;
        }

        /// <summary>
        /// Gamma function on real values
        /// </summary>
        public static double gamma(double v0) {
            return Check(v0, MathUtils.Gamma(v0));
        }

        /// <summary>
        /// Natural log of absolute value of Gamma function
        /// </summary>
        public static double lgamma(double v0) {
            return Check(v0, MathUtils.LogGamma(v0));
        }

        public static object trunc(CodeContext/*!*/ context, object value) {
            object func;
            if (PythonOps.TryGetBoundAttr(value, "__trunc__", out func)) {
                return PythonOps.CallWithContext(context, func);
            } else {
                throw PythonOps.TypeError("type {0} doesn't define __trunc__ method", PythonOps.GetPythonTypeName(value));
            }
        }

        public static bool isfinite(double x) {
            return !double.IsInfinity(x) && !double.IsNaN(x);
        }

        public static bool isinf(double v0) {
            return double.IsInfinity(v0);
        }

        public static bool isinf(BigInteger value) {
            return false;
        }

        public static bool isinf(object value) {
            // CPython tries float first, then double, so we need
            // an explicit overload which properly matches the order here
            double val;
            if (Converter.TryConvertToDouble(value, out val)) {
                return isinf(val);
            }
            return false;
        }

        public static bool isnan(double v0) {
            return double.IsNaN(v0);
        }

        public static bool isnan(BigInteger value) {
            return false;
        }

        public static bool isnan(object value) {
            // CPython tries float first, then double, so we need
            // an explicit overload which properly matches the order here
            double val;
            if (Converter.TryConvertToDouble(value, out val)) {
                return isnan(val);
            }
            return false;
        }

        public static double copysign(double x, double y) {
            return DoubleOps.CopySign(x, y);
        }

        public static double copysign(object x, object y) {
            double val, sign;
            if (!Converter.TryConvertToDouble(x, out val) ||
                !Converter.TryConvertToDouble(y, out sign)) {
                throw PythonOps.TypeError("TypeError: a float is required");
            }
            return DoubleOps.CopySign(val, sign);
        }

        // new in CPython 3.5
        public static object gcd(BigInteger x, BigInteger y) {
            var res = BigInteger.GreatestCommonDivisor(x, y);
            if (res.AsInt32(out var i))
                return i;
            return res;
        }

        public static object gcd(object x, object y) {
            return gcd(ObjectToBigInteger(x), ObjectToBigInteger(y));

            static BigInteger ObjectToBigInteger(object x) {
                BigInteger a;
                switch (PythonOps.Index(x)) {
                    case int i:
                        a = i;
                        break;
                    case BigInteger bi:
                        a = bi;
                        break;
                    default:
                        throw new InvalidOperationException();
                }
                return a;
            }
        }

        public static readonly double nan = double.NaN; // new in CPython 3.5

        public static readonly double inf = double.PositiveInfinity; // new in CPython 3.5

        // new in CPython 3.5
        public static bool isclose(double a, double b, double rel_tol = 1e-09, double abs_tol = 0.0) {
            if (rel_tol < 0 || abs_tol < 0) throw PythonOps.ValueError("tolerances must be non-negative");
            //if (double.IsNaN(a) || double.IsNaN(b)) return false;
            if (a == b) return true;
            if (double.IsInfinity(a) || double.IsInfinity(b)) return false;
            return Math.Abs(a - b) <= Math.Max(rel_tol * Math.Max(Math.Abs(a), Math.Abs(b)), abs_tol);
        }

        #region Private Implementation Details

        private static void SetExponentLe(byte[] v, int exp) {
            exp += Bias;
            ushort oldExp = LdExponentLe(v);
            ushort newExp = (ushort)(oldExp & 0x800f | (exp << 4));
            StExponentLe(v, newExp);
        }

        private static int IntExponentLe(byte[] v) {
            ushort exp = LdExponentLe(v);
            return ((int)((exp & 0x7FF0) >> 4) - Bias);
        }

        private static ushort LdExponentLe(byte[] v) {
            return (ushort)(v[6] | ((ushort)v[7] << 8));
        }

        private static long LdMantissaLe(byte[] v) {
            int i1 = (v[0] | (v[1] << 8) | (v[2] << 16) | (v[3] << 24));
            int i2 = (v[4] | (v[5] << 8) | ((v[6] & 0xF) << 16));

            return i1 | (i2 << 32);
        }

        private static void StExponentLe(byte[] v, ushort e) {
            v[6] = (byte)e;
            v[7] = (byte)(e >> 8);
        }

        private static bool IsDenormalizedLe(byte[] v) {
            ushort exp = LdExponentLe(v);
            long man = LdMantissaLe(v);

            return ((exp & 0x7FF0) == 0 && (man != 0));
        }

        private static void DecomposeLe(byte[] v, out double m, out int e) {
            if (IsDenormalizedLe(v)) {
                m = BitConverter.ToDouble(v, 0);
                m *= Math.Pow(2.0, 1022);
                v = BitConverter.GetBytes(m);
                e = IntExponentLe(v) - 1022;
            } else {
                e = IntExponentLe(v);
            }

            SetExponentLe(v, 0);
            m = BitConverter.ToDouble(v, 0);
        }

        private static double Check(double v) {
            return PythonOps.CheckMath(v);
        }

        private static double Check(double input, double output) {
            return PythonOps.CheckMath(input, output);
        }

        private static double Check(double in0, double in1, double output) {
            return PythonOps.CheckMath(in0, in1, output);
        }

        #endregion
    }
}
