File: Utilities\MathUtils.cs
Web Access
Project: src\src\Microsoft.ML.Core\Microsoft.ML.Core.csproj (Microsoft.ML.Core)
// Licensed to the .NET Foundation under one or more agreements.
// The .NET Foundation licenses this file to you under the MIT license.
// See the LICENSE file in the project root for more information.
 
using System;
using System.Collections.Generic;
using Microsoft.ML.Runtime;
 
namespace Microsoft.ML.Internal.Utilities
{
    /// <summary>
    /// Some useful math methods.
    /// </summary>
    [BestFriend]
    internal static class MathUtils
    {
        public static float ToFloat(this Double dbl)
        {
            return (Single)dbl;
        }
 
        // The purpose of this is to catch (at compile time) invocations of ToFloat
        // that are not appropriate. Note that the return type is void.
        public static void ToFloat(this Single dbl)
        {
            Contracts.Assert(false, "Bad use of ToFloat");
            throw Contracts.Except();
        }
 
        public static float Sqrt(float x)
        {
            return Math.Sqrt(x).ToFloat();
        }
 
        public static float Log(float x)
        {
            return Math.Log(x).ToFloat();
        }
 
        public static float Log(float a, float newBase)
        {
            return Math.Log(a, newBase).ToFloat();
        }
 
        public static float Pow(float x, float y)
        {
            return Math.Pow(x, y).ToFloat();
        }
 
        /// <summary>
        /// Finds the best least-squares fit of y = ax + b
        /// </summary>
        /// <param name="x">The x values.</param>
        /// <param name="y">The y values.</param>
        /// <param name="a">The coefficent a.</param>
        /// <param name="b">The intercept b.</param>
        public static void SimpleLinearRegression(float[] x, float[] y, out float a, out float b)
        {
            Contracts.CheckValue(x, nameof(x));
            Contracts.CheckValue(y, nameof(y));
            Contracts.Check(x.Length == y.Length, "Input and output must be same length.");
 
            int m = x.Length;
 
            float sumSqX = 0;
            float sumX = 0;
            float sumXY = 0;
            float sumY = 0;
            float sumSqY = 0;
            for (int i = 0; i < m; i++)
            {
                float xVal = x[i];
                float yVal = y[i];
                sumSqX += xVal * xVal;
                sumX += xVal;
                sumXY += xVal * yVal;
                sumY += yVal;
                sumSqY += yVal * yVal;
            }
 
            float denom = sumSqX * m - sumX * sumX;
            a = (sumXY * m - sumY * sumX) / denom;
            b = (sumSqX * sumY - sumXY * sumX) / denom;
        }
 
        /************************
         * MATH ARRAY FUNCTIONS *
         ************************/
        // REVIEW: Move the vector methods into VectorUtils.
 
        /// <summary>
        /// The product of elements in a
        /// </summary>
        /// <param name="a">an array</param>
        /// <returns>the product of a's elements</returns>
        public static int Product(int[] a)
        {
            Contracts.AssertValue(a);
            int result = 1;
            foreach (var x in a)
                result = checked(result * x);
            return result;
        }
 
        /// <summary>
        /// Find the max element of a
        /// </summary>
        /// <param name="a">an array</param>
        /// <returns>the max element</returns>
        public static float Max(float[] a)
        {
            Contracts.AssertValue(a);
            float result = float.NegativeInfinity;
            foreach (var x in a)
                result = Math.Max(result, x);
            return result;
        }
 
        /// <summary>
        /// Find the minimum element of a
        /// </summary>
        /// <param name="a">an array</param>
        /// <returns>the minimum element</returns>
        public static float Min(float[] a)
        {
            Contracts.AssertValue(a);
            float result = float.PositiveInfinity;
            foreach (var x in a)
                result = Math.Min(result, x);
            return result;
        }
 
        /// <summary>
        /// Finds the first index of the max element of the span.
        /// NaNs are ignored. If all the elements to consider are NaNs, -1 is
        /// returned. The caller should distinguish in this case between two
        /// possibilities:
        /// 1) The number of the element to consider is zero.
        /// 2) All the elements to consider are NaNs.
        /// </summary>
        /// <param name="a">The span of floats.</param>
        /// <returns>the first index of the max element</returns>
        public static int ArgMax(ReadOnlySpan<float> a)
        {
            if (a.IsEmpty)
                return -1;
 
            int amax = -1;
            float max = float.NegativeInfinity;
            for (int i = a.Length - 1; i >= 0; i--)
            {
                if (max <= a[i])
                {
                    amax = i;
                    max = a[i];
                }
            }
 
            return amax;
        }
 
        /// <summary>
        /// Finds the first index of the minimum element of the span.
        /// NaNs are ignored. If all the elements to consider are NaNs, -1 is
        /// returned. The caller should distinguish in this case between two
        /// possibilities:
        /// 1) The number of the element to consider is zero.
        /// 2) All the elements to consider are NaNs.
        /// </summary>
        /// <param name="a">The span of floats.</param>
        /// <returns>the first index of the minimum element</returns>
        public static int ArgMin(ReadOnlySpan<float> a)
        {
            if (a.IsEmpty)
                return -1;
 
            int amin = -1;
            float min = float.PositiveInfinity;
            for (int i = a.Length - 1; i >= 0; i--)
            {
                if (min >= a[i])
                {
                    amin = i;
                    min = a[i];
                }
            }
 
            return amin;
        }
 
        /*****************
         * LOG FUNCTIONS *
         *****************/
 
        private const float LogTolerance = 30;
 
        /// <summary>
        /// computes the "softmax" function: log sum_i exp x_i
        /// </summary>
        /// <param name="inputs">Span of numbers to softmax</param>
        /// <returns>the softmax of the numbers</returns>
        /// <remarks>may have slightly lower roundoff error if inputs are sorted, smallest first</remarks>
        public static float SoftMax(ReadOnlySpan<float> inputs)
        {
            int maxIdx = 0;
            float max = float.NegativeInfinity;
            for (int i = 0; i < inputs.Length; i++)
            {
                if (inputs[i] > max)
                {
                    maxIdx = i;
                    max = inputs[i];
                }
            }
 
            if (float.IsNegativeInfinity(max))
                return float.NegativeInfinity;
 
            if (inputs.Length == 1)
                return max;
 
            double intermediate = 0.0;
            float cutoff = max - LogTolerance;
 
            for (int i = 0; i < inputs.Length; i++)
            {
                if (i == maxIdx)
                    continue;
                if (inputs[i] > cutoff)
                    intermediate += Math.Exp(inputs[i] - max);
            }
 
            if (intermediate > 0.0)
                return (float)(max + Math.Log(1.0 + intermediate));
            return max;
        }
 
        /// <summary>
        /// computes "softmax" function of two arguments: log (exp x + exp y)
        /// </summary>
        public static float SoftMax(float lx, float ly)
        {
            float max;
            float negDiff;
            if (lx > ly)
            {
                max = lx;
                negDiff = ly - lx;
            }
            else
            {
                max = ly;
                negDiff = lx - ly;
            }
            if (float.IsNegativeInfinity(max) || negDiff < -LogTolerance)
            {
                return max;
            }
            else
            {
                return (float)(max + Math.Log(1.0 + Math.Exp(negDiff)));
            }
        }
 
        /*******************
         * OTHER FUNCTIONS *
         *******************/
 
        public const float DefaultMaxRelativeErr = (float)1e-8;
        public const float DefaultMaxAbsErr = (float)1e-12;
 
        /// <summary>
        /// true if two float values are close (using relative comparison)
        /// </summary>
        /// <param name="a"></param>
        /// <param name="b"></param>
        /// <returns></returns>
        public static bool AlmostEqual(float a, float b)
        {
            return AlmostEqual(a, b, DefaultMaxRelativeErr, DefaultMaxAbsErr);
        }
 
        public static bool AlmostEqual(float a, float b, float maxRelErr, float maxAbsError)
        {
            Contracts.Assert(FloatUtils.IsFinite(maxRelErr));
            Contracts.Assert(FloatUtils.IsFinite(maxAbsError));
 
            float absDiff = Math.Abs(a - b);
            if (absDiff < maxAbsError)
                return true;
            float maxAbs = Math.Max(Math.Abs(a), Math.Abs(b));
            return (absDiff / maxAbs) <= maxRelErr;
        }
 
        private static readonly int[] _possiblePrimeMod30 = new int[] { 1, 7, 11, 13, 17, 19, 23, 29 };
        private static readonly double _constantForLogGamma = 0.5 * Math.Log(2 * Math.PI);
        private static readonly double[] _coeffsForLogGamma = { 12.0, -360.0, 1260.0, -1680.0, 1188.0 };
 
        /// <summary>
        /// Returns the log of the gamma function, using the Stirling approximation
        /// </summary>
        /// <param name="x">Argument of function</param>
        /// <returns>Log Gamma(x)</returns>
        /// <remarks>Accurate to eight digits for all x.</remarks>
        public static double LogGamma(double x)
        {
            Contracts.CheckParam(x > 0, nameof(x), "LogGamma invalid for x <= 0");
 
            double res = 0;
            if (x < 6)
            {
                int toAdd = (int)Math.Floor(7 - x);
                double v2 = 1;
                for (int i = 0; i < toAdd; i++)
                    v2 *= (x + i);
                res = -Math.Log(v2);
                x += toAdd;
            }
            x = x - 1;
 
            res += _constantForLogGamma + (x + 0.5) * Math.Log(x) - x;
 
            // correction terms
            double xSquared = x * x;
            double pow = x;
            foreach (double coeff in _coeffsForLogGamma)
            {
                double newRes = res + 1.0 / (coeff * pow);
                if (newRes == res)
                {
                    return res;
                }
                res = newRes;
                pow *= xSquared;
            }
 
            return res;
        }
 
        private static List<double> _logFactorialCache;
        private const int LogFactorialCacheSize = 1000;
 
        /// <summary>
        /// Computes the log factorial of n, using fast methods
        /// </summary>
        /// <param name="n">The number to compute the factorial of</param>
        /// <returns>The log factorial of n</returns>
        public static double LogFactorial(int n)
        {
            Contracts.CheckParam(n >= 0, nameof(n), "LogFactorial is invalid for n < 0.");
 
            if (n >= LogFactorialCacheSize)
                return LogGamma(n + 1);
 
            if (_logFactorialCache == null)
            {
                _logFactorialCache = new List<double>(LogFactorialCacheSize);
                _logFactorialCache.Add(0);
                _logFactorialCache.Add(0);
            }
 
            for (int i = _logFactorialCache.Count; i <= n; i++)
            {
                _logFactorialCache.Add(_logFactorialCache[i - 1] + Math.Log(i));
            }
 
            return _logFactorialCache[n];
        }
 
        /// <summary>
        /// Returns the two-tailed p-value given a t statistic from a distribution
        /// parameterized by the provided number of degrees of freedom.
        /// </summary>
        /// <param name="t">The t-statistic</param>
        /// <param name="df">The degrees of freedom</param>
        /// <returns>The corresponding two-tailed p-value</returns>
        public static Double TStatisticToPValue(Double t, Double df)
        {
            Contracts.CheckParam(df > 0, nameof(df), "Degrees of freedom must be positive");
 
            // REVIEW: Of some interest is calculating p-values for infinite
            // df, but this code will not handle that. In fact I have strong concerns
            // for the numeric stability of this method for larger df, since we have
            // a form of catastrophic cancellation for t. It may be worthwhile folding
            // the terms for incomplete beta into this evaluation, if this becomes
            // problematic.
            Double result = IncompleteBeta(df / (df + t * t), df / 2, 0.5);
            // Clamp the output to 0 to 1, in case for numerical reasons it strayed
            // outside these bounds.
            return Math.Max(0, Math.Min(result, 1));
        }
 
        private delegate Double Sequence(int i);
 
        /// <summary>
        /// Lentz's algorithm for evaluating the continued fraction
        /// b0 + a1 / (b1 + a2 / (b2 + a3 / (b3 + a4 / ...) ) )
        /// </summary>
        /// <param name="a">The <c>a</c> function mapping positive integers to a sequence term</param>
        /// <param name="b">The <c>b</c> function mapping non-negative integers to a sequence term</param>
        /// <param name="tol">Calculate the continued fraction to this tolerance</param>
        /// <returns>The evaluation of the continued fraction</returns>
        private static Double Lentz(Sequence a, Sequence b, Double tol = 1e-15)
        {
            Contracts.AssertValue(a);
            Contracts.AssertValue(b);
            Contracts.Assert(tol > 0);
 
            Double f = Unclamp(b(0));
            Double c = f;
            Double d = 0;
            const int iterMax = 100000;
            for (int i = 1; i < iterMax; ++i)
            {
                Double bi = b(i);
                Double ai = a(i);
                d = 1.0 / Unclamp(bi + ai * d);
                c = Unclamp(bi + ai / c);
                Double ratio = c * d;
                f *= ratio;
                if (Math.Abs(ratio - 1.0) < tol)
                    break;
            }
            return f;
        }
 
        private static Double Beta(Double a, Double b)
        {
            // REVIEW: LogGamma implementation precision is a concern, but
            // is minor, perhaps, compared to the instability of TtoP.
            return Math.Exp(LogGamma(a) + LogGamma(b) - LogGamma(a + b));
        }
 
        private static Double IncompleteBeta(Double x, Double a, Double b)
        {
            Contracts.Assert(0 <= x && x <= 1);
            Contracts.Assert(0 < a);
            Contracts.Assert(0 < b);
            if (x == 0 || x == 1)
                return x;
 
            // This implementation of the incomplete beta will converge quickly
            // only if the minimum of a and b is relatively small, that is, less than
            // few thousand or so -- otherwise it converges slowly. If this method
            // is made public, then this case of larger a and b should be handled
            // better.
 
            // If x > (a+1)/(a+b+2), exploit I_x(a,b) = 1 - I_{x-1}(b,a)...
            // REVIEW: Can we be sure some numerical imprecision "weirdness"
            // won't lead both the above test and (1-x) > (b+1)/(a+b+2) being true?
            // Should we build in some sort of protections?
            if (x * (a + b + 2) > a + 1)
                return 1 - IncompleteBeta(1 - x, b, a);
 
            // This is a fairly naive implementation of I_x(a,b), using a generic
            // non-purpose optimized application of continued fractions.
            Sequence adel = i =>
            {
                if (i == 1)
                    return 1;
                int m = (i - 1) >> 1;
                Double denom = ((a + i - 2) * (a + i - 1));
                if ((i & 1) == 0)
                    return -(a + m) * (a + b + m) * x / denom;
                return m * (b - m) * x / denom;
            };
            Sequence bdel = i => i == 0 ? 0 : 1;
            return Math.Pow(x, a) * Math.Pow(1 - x, b) / (a * Beta(a, b)) * Lentz(adel, bdel);
        }
 
        private static Double Unclamp(Double val)
        {
            const Double bound = 1e-30;
            if (!(-bound <= val && val <= bound))
                return val;
            return val < 0 ? -bound : bound;
        }
 
        /// <summary>
        /// The logistic sigmoid function: 1 / (1 + e^(-x)).
        /// </summary>
        public static float Sigmoid(float x)
        {
#if SLOW_EXP
            return SigmoidSlow(x);
#else
            return SigmoidFast(x);
#endif
        }
 
        /// <summary>
        /// Hyperbolic tangent.
        /// </summary>
        public static float Tanh(float x)
        {
#if SLOW_EXP
            return TanhSlow(x);
#else
            return TanhFast(x);
#endif
        }
 
        /// <summary>
        /// The logistic sigmoid function: 1 / (1 + e^(-x)).
        /// </summary>
        public static float SigmoidSlow(float x)
        {
            // The following two expressions are mathematically equivalent. Due to the potential of getting overflow we should
            // not call exp(x) for large positive x: instead, we modify the expression to compute exp(-x).
            if (x > 0)
                return 1 / (1 + ExpSlow(-x));
            else
            {
                var ex = ExpSlow(x);
                return ex / (1 + ex);
            }
        }
 
        /// <summary>
        /// Hyperbolic tangent.
        /// </summary>
        public static float TanhSlow(float x)
        {
            return Math.Tanh(x).ToFloat();
        }
 
        /// <summary>
        /// The exponential function: e^(x).
        /// </summary>
        public static float ExpSlow(float x)
        {
            return Math.Exp(x).ToFloat();
        }
 
        private const int ExpInf = 128;
        private const float Coef1 = (float)0.013555747234814917704030793;
        private const float Coef2 = (float)0.065588116243247810171479524;
        private const float Coef3 = (float)0.3069678791803394491901401;
 
        // 1 / ln(2).
        private const float RecipLn2 = (float)1.44269504088896340735992468100;
 
        private static float PowerOfTwo(int exp)
        {
            Contracts.Assert(0 <= exp && exp < ExpInf);
            return FloatUtils.GetPowerOfTwoSingle(exp);
        }
 
        /// <summary>
        /// The logistic sigmoid function: 1 / (1 + e^(-x)).
        /// </summary>
        public static float SigmoidFast(float x)
        {
            // This is a loose translation from SSE code
            if (float.IsNaN(x))
                return x;
 
            bool neg = false;
            if (x < 0)
            {
                x = -x;
                neg = true;
            }
 
            // Multiply by 1/ln(2).
            x *= RecipLn2;
            if (x >= ExpInf)
                return neg ? (float)0 : (float)1;
 
            // Get the floor and fractional part.
            int n = (int)x;
            Contracts.Assert(0 <= n && n < ExpInf);
            float f = x - n;
            Contracts.Assert(0 <= f && f < 1);
 
            // Get the integer power of two part.
            float r = PowerOfTwo(n);
            Contracts.Assert(1 <= r && r < float.PositiveInfinity);
 
            // This approximates 2^f for 0 <= f <= 1. Note that it is exact at the endpoints.
            float res = 1 + f + (f - 1) * f * ((Coef1 * f + Coef2) * f + Coef3);
 
            res = 1 / (1 + r * res);
            if (!neg)
                res = 1 - res;
            return res;
        }
 
        /// <summary>
        /// The hyperbolic tangent function.
        /// </summary>
        public static float TanhFast(float x)
        {
            if (float.IsNaN(x))
                return x;
 
            bool neg = false;
            if (x < 0)
            {
                x = -x;
                neg = true;
            }
 
            // Multiply by 2/ln(2).
            x *= 2 * RecipLn2;
            if (x >= ExpInf)
                return neg ? (float)(-1) : (float)1;
 
            // Get the floor and fractional part.
            int n = (int)x;
            Contracts.Assert(0 <= n && n < ExpInf);
            float f = x - n;
            Contracts.Assert(0 <= f && f < 1);
 
            // Get the integer power of two part.
            float r = PowerOfTwo(n);
            Contracts.Assert(1 <= r && r < float.PositiveInfinity);
 
            // This approximates 2^f - 1 for 0 <= f <= 1. Note that it is exact at the endpoints.
            float res = f + (f - 1) * f * ((Coef1 * f + Coef2) * f + Coef3);
 
            res *= r;
            res = (res + (r - 1)) / (res + (r + 1));
            if (neg)
                res = -res;
            return res;
        }
 
        /// <summary>
        /// The exponential function: e^(x).
        /// </summary>
        public static float ExpFast(float x)
        {
            if (float.IsNaN(x))
                return x;
 
            bool neg = false;
            if (x < 0)
            {
                x = -x;
                neg = true;
            }
 
            // Multiply by 1/ln(2). Then we need to calculate 2^x.
            x *= RecipLn2;
            if (x >= ExpInf)
                return neg ? (float)0 : float.PositiveInfinity;
 
            // Get the floor and fractional part.
            int n = (int)x;
            Contracts.Assert(0 <= n && n < ExpInf);
            float f = x - n;
            Contracts.Assert(0 <= f && f < 1);
 
            // Get the integer power of two part.
            float r = PowerOfTwo(n);
            Contracts.Assert(1 <= r && r < float.PositiveInfinity);
 
            // This approximates 2^f for 0 <= f <= 1. Note that it is exact at the endpoints.
            float res = 1 + f + (f - 1) * f * ((Coef1 * f + Coef2) * f + Coef3);
 
            res *= r;
            if (neg)
                res = 1 / res;
            return res;
        }
 
        /// <summary>
        /// Apply a soft max on an array of floats. Note that src and dst may be the same array.
        /// </summary>
        public static void ApplySoftMax(float[] src, float[] dst)
        {
            Contracts.Assert(src.Length == dst.Length);
            ApplySoftMax(src, dst, 0, src.Length);
        }
 
        /// <summary>
        /// Apply a soft max on a range within an array of floats. Note that src and dst may be the same array.
        /// </summary>
        public static void ApplySoftMax(float[] src, float[] dst, int start, int end)
        {
            Contracts.Assert(src.Length == dst.Length);
            Contracts.Assert(0 <= start && start <= end && end <= src.Length);
 
            // Compute max output.
            float maxOut = float.NegativeInfinity;
            for (int i = start; i < end; i++)
                maxOut = Math.Max(maxOut, src[i]);
 
            // Compute exp and sum.
            float sum = 0;
            for (int i = start; i < end; i++)
            {
                dst[i] = ExpFast(src[i] - maxOut);
                sum += dst[i];
            }
 
            // Normalize.
            for (int i = start; i < end; i++)
                dst[i] /= sum;
        }
 
        public static float GetMedianInPlace(float[] src, int count)
        {
            Contracts.Assert(count >= 0);
            Contracts.Assert(Utils.Size(src) >= count);
 
            if (count == 0)
                return float.NaN;
 
            Array.Sort(src, 0, count);
 
            // Skip any NaNs. They sort to the low end.
            int ivMin = 0;
            int ivLim = count;
            while (ivMin < ivLim && float.IsNaN(src[ivMin]))
                ivMin++;
            Contracts.Assert(ivMin <= ivLim);
 
            if (ivMin >= ivLim)
                return float.NaN;
 
            // This assert will fire if Array.Sort changes to put NaNs at the high end.
            Contracts.Assert(!float.IsNaN(src[ivLim - 1]));
 
            // If we're dealing with an odd number of things, just grab the middel item; otherwise,
            // average the two middle items.
            uint cv = (uint)ivMin + (uint)ivLim;
            int iv = (int)(cv / 2);
            if ((cv & 1) != 0)
                return src[iv];
            return (src[iv - 1] + src[iv]) / 2;
        }
 
        public static Double CosineSimilarity(ReadOnlySpan<float> a, ReadOnlySpan<float> b, int aIdx, int bIdx, int len)
        {
            const Double epsilon = 1e-12f;
            Contracts.Assert(len > 0);
            Contracts.Assert(aIdx >= 0 && aIdx <= a.Length - len);
            Contracts.Assert(bIdx >= 0 && bIdx <= b.Length - len);
 
            Double ab = 0;
            Double a2 = 0;
            Double b2 = 0;
 
            for (int lim = aIdx + len; aIdx < lim; aIdx++, bIdx++)
            {
                ab += (Double)a[aIdx] * b[bIdx];
                a2 += (Double)a[aIdx] * a[aIdx];
                b2 += (Double)b[bIdx] * b[bIdx];
            }
 
            Double similarity = ab / (Math.Sqrt(a2 * b2) + epsilon);
            Contracts.Assert(-1 - epsilon <= similarity && similarity <= 1 + epsilon);
            if (Math.Abs(similarity) > 1)
                return similarity > 1 ? 1 : -1;
 
            return similarity;
        }
 
        /// <summary>
        /// Entropy of a given probability
        /// </summary>
        public static Double Entropy(Double prob, bool useLnNotLog2 = false)
        {
            if (prob < 0 || prob > 1)
                return Double.NaN;
            if (prob == 0.0 || prob == 1.0)
                return 0.0;
            return
                useLnNotLog2
                ? -prob * Math.Log(prob) - (1 - prob) * Math.Log(1 - prob)
                : -prob * Math.Log(prob, 2) - (1 - prob) * Math.Log(1 - prob, 2);
        }
 
        /// <summary>
        /// Cross-entropy of two distributions
        /// </summary>
        public static Double CrossEntropy(Double probTrue, Double probPredicted, bool useLnNotLog2 = false)
        {
            if (probTrue < 0 || probTrue > 1 || probPredicted < 0 || probPredicted > 1)
                return Double.NaN;
            if ((probPredicted == 0.0 || probPredicted == 1.0) && (probPredicted == probTrue))
                return 0.0;
            return
                useLnNotLog2
                ? -probTrue * Math.Log(probPredicted) - (1 - probTrue) * Math.Log(1 - probPredicted)
                : -probTrue * Math.Log(probPredicted, 2) - (1 - probTrue) * Math.Log(1 - probPredicted, 2);
        }
 
        /// <summary>
        /// Given a set of values <c>Ln(a1), Ln(a2), ... Ln(an)</c>,
        /// return <c>Ln(a1+a2+...+an)</c>. This is especially useful
        /// when working with log probabilities and likelihoods.
        /// </summary>
        /// <param name="terms"></param>
        public static float LnSum(IEnumerable<float> terms)
        {
            // Two passes to find the overall max is a *lot* simpler,
            // but potentially more computationally intensive.
            float max = float.NegativeInfinity;
            Double soFar = 0;
 
            foreach (float term in terms)
            {
                // At this point, all *prior* terms, Math.Exp(x - max).
                if (float.IsNegativeInfinity(term))
                    continue;
                if (!(term > max))
                    soFar += Math.Exp(term - max);
                else
                {
                    soFar = Math.Exp(max - term) * soFar + 1;
                    max = term;
                }
            }
            return (float)Math.Log(soFar) + max;
        }
 
        /// <summary>
        /// Math.Sin returns the input value for inputs with large magnitude. We return NaN instead, for consistency
        /// with Math.Sin(infinity).
        /// </summary>
        public static double Sin(double a)
        {
            var res = Math.Sin(a);
            return Math.Abs(res) > 1 ? double.NaN : res;
        }
 
        /// <summary>
        /// Math.Cos returns the input value for inputs with large magnitude. We return NaN instead, for consistency
        /// with Math.Cos(infinity).
        /// </summary>
        public static double Cos(double a)
        {
            var res = Math.Cos(a);
            return Math.Abs(res) > 1 ? double.NaN : res;
        }
 
        /// <summary>
        /// Returns the smallest integral value that is greater than or equal to the result of the division.
        /// </summary>
        /// <param name="numerator">Number to be divided.</param>
        /// <param name="denomenator">Number with which to divide the numerator.</param>
        /// <returns></returns>
        public static long DivisionCeiling(long numerator, long denomenator)
        {
            return (checked(numerator + denomenator) - 1) / denomenator;
        }
    }
}