System.Math.Log(double)

public static double Log(double d) {
            return Math.Log(d);
        }

public double GetGaussianRandomNumber(double mean, double stdDev)
        {
            double u1 = _random.NextDouble(); //these are uniform(0,1) random doubles
            double u2 = _random.NextDouble();
            double randStdNormal = Math.Sqrt(-2.0 * Math.Log(u1)) *
                         Math.Sin(2.0 * Math.PI * u2); //random normal(0,1)
            double randNormal =
                         mean + stdDev * randStdNormal; //random normal(mean,stdDev^2)
            return randNormal;
        }

private void InitializeSampleDataContext() {
			for (int i = 0; i < 1000; i++) {
				int modulus = i % 2;
				double xm = i / (20.0d + 3 * modulus);
				double ym = 10.0d + 2 * modulus;
				double x = random.NextDouble() * xm + 1;
				double y = Math.Log(ym * (x - 1.0) + 1.0) * (random.NextDouble() + 0.9);

				if (modulus == 0)
					primaryChartPoints1.Add(new Point(x, y));
				else
					primaryChartPoints2.Add(new Point(x, y));
			}
		}

public IList<Point<decimal, decimal>> GetHazardFunctionOverEachPeriod()
        {
            nelsonAalen = new KaplanMeier(observationsInternal);
            var result = new List<Point<decimal, decimal>>();
            var cnt = 0;

            foreach (var set in observationsInternal)
            {
                var survival = nelsonAalen.GetSurvivalOverPeriod(cnt);
                var v1 = survival == 0 ? 0 : (-1) * Math.Log((double)survival);
                var v2 = (double)(set.unitTime);
                var div = (decimal)((-1) * v1 / v2);

                var point = new Point<decimal, decimal>(set.unitTime * cnt, (div));
                result.Add(point);

                cnt++;

            }
            return result;
        }

private static void GammaDistributionFit(IList<decimal> points)
        {
            var sum = 0m;
            var sumxLNx = 0m;
            var sumLNx = 0m;
            int n = points.Count();
            foreach (var d in points)
            {
                sum += d;
                sumxLNx += d * ((decimal)Math.Log((double)d));
                sumLNx += (decimal)Math.Log((double)d);
            }

            k = ((n * sum) / ((n * sumxLNx) - (sumLNx * sum)));
            α = (decimal)((1 / (Math.Pow((double)n, 2))) * (double)((n * sumxLNx) - (sumLNx * sum)));
        }

private double SoftPlusFunc(double x)
        {
            return Math.Log(1 + Math.Exp(x));
        }

public static BigInteger GetApproximateValueFromIndex(UInt64 n)
		{
			if (n < 6)
			{
				return primes[(int)n];
			}

			double fn = (double)n;
			double flogn = Math.Log(n);
			double flog2n = Math.Log(flogn);

			double upper;

			if (n >= 688383)    /* Dusart 2010 page 2 */
			{
				upper = fn * (flogn + flog2n - 1.0 + ((flog2n - 2.00) / flogn));
			}
			else if (n >= 178974)    /* Dusart 2010 page 7 */
			{
				upper = fn * (flogn + flog2n - 1.0 + ((flog2n - 1.95) / flogn));
			}
			else if (n >= 39017)    /* Dusart 1999 page 14 */
			{
				upper = fn * (flogn + flog2n - 0.9484);
			}
			else                    /* Modified from Robin 1983 for 6-39016 _only_ */
			{
				upper = fn * (flogn + 0.6000 * flog2n);
			}

			if (upper >= (double)UInt64.MaxValue)
			{
				throw new OverflowException($"{upper} > {UInt64.MaxValue}");
			}

			return new BigInteger((UInt64)Math.Ceiling(upper));
		}

public string ToString(int precision, IFormatProvider formatProvider = null)
		{
			var pow = Math.Floor((_value > 0 ? Math.Log(_value) : 0) / Math.Log(1024));

			pow = Math.Min(pow, _units.Length - 1);

			var value = _value / Math.Pow(1024, pow);

			var precisionString = formatProvider == null
				? precision.ToString(CultureInfo.CurrentCulture)
				: precision.ToString(formatProvider);

			return value.ToString(Math.Abs(pow - 0) < double.Epsilon ? "F0" : "F" + precisionString) + " " + _units[(int)pow];
		}

[ContextMethod("Логарифм", "Log")]
        public double Log(double p1)
        {
            return System.Math.Log(p1);
        }

public static Vector2d GeoToWorldPosition(double lat, double lon, Vector2d refPoint, float scale = 1)
		{
			var posx = lon * OriginShift / 180;
			var posy = Math.Log(Math.Tan((90 + lat) * Math.PI / 360)) / (Math.PI / 180);
			posy = posy * OriginShift / 180;
			return new Vector2d((posx - refPoint.x) * scale, (posy - refPoint.y) * scale);
		}

public static Vector2d LatLonToMeters(double lat, double lon)
		{
			var posx = lon * OriginShift / 180;
			var posy = Math.Log(Math.Tan((90 + lat) * Math.PI / 360)) / (Math.PI / 180);
			posy = posy * OriginShift / 180;
			return new Vector2d(posx, posy);
		}

public static UnwrappedTileId LareplacedudeLongitudeToTileId(double lareplacedude, double longitude, int zoom)
		{
			var x = (int)Math.Floor((longitude + 180.0) / 360.0 * Math.Pow(2.0, zoom));
			var y = (int)Math.Floor((1.0 - Math.Log(Math.Tan(lareplacedude * Math.PI / 180.0)
					+ 1.0 / Math.Cos(lareplacedude * Math.PI / 180.0)) / Math.PI) / 2.0 * Math.Pow(2.0, zoom));

			return new UnwrappedTileId(zoom, x, y);
		}

public static double latToY(double lat)
    {
        lat = Math.Min(89.5, Math.Max(lat, -89.5));
        double phi = DegToRad(lat);
        double sinphi = Math.Sin(phi);
        double con = ECCENT * sinphi;
        con = Math.Pow(((1.0 - con) / (1.0 + con)), COM);
        double ts = Math.Tan(0.5 * ((Math.PI * 0.5) - phi)) / con;
        return 0 - R_MAJOR * Math.Log(ts);
    }

public override void Render(IRenderContext rc, PlotModel model)
        {
            base.Render(rc, model);
            var clippingRect = this.GetClippingRect();
            var lon0 = this.XAxis.ActualMinimum;
            var lon1 = this.XAxis.ActualMaximum;
            var lat0 = this.YAxis.ActualMinimum;
            var lat1 = this.YAxis.ActualMaximum;

            // the desired number of tiles horizontally
            double tilesx = model.Width / this.TileSize;

            // calculate the desired zoom level
            var n = tilesx / (((lon1 + 180) / 360) - ((lon0 + 180) / 360));
            var zoom = (int)Math.Round(Math.Log(n) / Math.Log(2));
            if (zoom < this.MinZoomLevel)
            {
                zoom = this.MinZoomLevel;
            }

            if (zoom > this.MaxZoomLevel)
            {
                zoom = this.MaxZoomLevel;
            }

            // find tile coordinates for the corners
            double x0, y0;
            LatLonToTile(lat0, lon0, zoom, out x0, out y0);
            double x1, y1;
            LatLonToTile(lat1, lon1, zoom, out x1, out y1);

            double xmax = Math.Max(x0, x1);
            double xmin = Math.Min(x0, x1);
            double ymax = Math.Max(y0, y1);
            double ymin = Math.Min(y0, y1);

            // Add the tiles
            for (var x = (int)xmin; x < xmax; x++)
            {
                for (var y = (int)ymin; y < ymax; y++)
                {
                    string uri = this.GetTileUri(x, y, zoom);
                    var img = this.GetImage(uri, rc.RendersToScreen);

                    if (img == null)
                    {
                        continue;
                    }

                    // transform from tile coordinates to lat/lon
                    double lareplacedude0, lareplacedude1, longitude0, longitude1;
                    TileToLatLon(x, y, zoom, out lareplacedude0, out longitude0);
                    TileToLatLon(x + 1, y + 1, zoom, out lareplacedude1, out longitude1);

                    // transform from lat/lon to screen coordinates
                    var s00 = this.Transform(longitude0, lareplacedude0);
                    var s11 = this.Transform(longitude1, lareplacedude1);

                    var r = OxyRect.Create(s00.X, s00.Y, s11.X, s11.Y);

                    // draw the image
                    rc.DrawClippedImage(clippingRect, img, r.Left, r.Top, r.Width, r.Height, this.Opacity, true);
                }
            }

            // draw the copyright notice
            var p = new ScreenPoint(clippingRect.Right - 5, clippingRect.Bottom - 5);
            var textSize = rc.MeasureText(this.CopyrightNotice, null, 12);
            rc.DrawRectangle(new OxyRect(p.X - textSize.Width - 2, p.Y - textSize.Height - 2, textSize.Width + 4, textSize.Height + 4), OxyColors.White.ChangeAlpha(200), null);

            rc.DrawText(
                p,
                this.CopyrightNotice,
                OxyColors.Black,
                null,
                12,
                500,
                0,
                HorizontalAlignment.Right,
                VerticalAlignment.Bottom);
        }

private static void LatLonToTile(double lareplacedude, double longitude, int zoom, out double x, out double y)
        {
            // http://wiki.openstreetmap.org/wiki/Slippy_map_tilenames            
            int n = 1 << zoom;
            double lat = lareplacedude / 180 * Math.PI;
            x = (longitude + 180.0) / 360.0 * n;
            y = (1.0 - Math.Log(Math.Tan(lat) + 1.0 / Math.Cos(lat)) / Math.PI) / 2.0 * n;
        }

public override void GetTickValues(
            out IList<double> majorLabelValues, out IList<double> majorTickValues, out IList<double> minorTickValues)
        {
            if (this.ActualMinimum <= 0)
            {
                this.ActualMinimum = 0.1;
            }

            double logBase = Math.Log(this.Base);
            var e0 = (int)Math.Floor(Math.Log(this.ActualMinimum) / logBase);
            var e1 = (int)Math.Ceiling(Math.Log(this.ActualMaximum) / logBase);

            // find the min & max values for the specified base
            // round to max 10 digits
            double p0 = Math.Pow(this.Base, e0);
            double p1 = Math.Pow(this.Base, e1);
            double d0 = Math.Round(p0, 10);
            double d1 = Math.Round(p1, 10);
            if (d0 <= 0)
            {
                d0 = p0;
            }

            double d = d0;
            majorTickValues = new List<double>();
            minorTickValues = new List<double>();

            double epsMin = this.ActualMinimum * 1e-6;
            double epsMax = this.ActualMaximum * 1e-6;

            while (d <= d1 + epsMax)
            {
                // d = RemoveNoiseFromDoubleMath(d);
                if (d >= this.ActualMinimum - epsMin && d <= this.ActualMaximum + epsMax)
                {
                    majorTickValues.Add(d);
                }

                for (int i = 1; i < this.Base; i++)
                {
                    double d2 = d * (i + 1);
                    if (d2 > d1 + double.Epsilon)
                    {
                        break;
                    }

                    if (d2 > this.ActualMaximum)
                    {
                        break;
                    }

                    if (d2 >= this.ActualMinimum && d2 <= this.ActualMaximum)
                    {
                        minorTickValues.Add(d2);
                    }
                }

                d *= this.Base;
                if (double.IsInfinity(d))
                {
                    break;
                }

                if (d < double.Epsilon)
                {
                    break;
                }

                if (double.IsNaN(d))
                {
                    break;
                }
            }

            if (majorTickValues.Count < 2)
            {
                base.GetTickValues(out majorLabelValues, out majorTickValues, out minorTickValues);
            }
            else
            {
                majorLabelValues = majorTickValues;
            }
        }

public override double Transform(double x)
        {
            Debug.replacedert(x > 0, "Value should be positive.");
            if (x <= 0)
            {
                return -1;
            }

            return (Math.Log(x) - this.offset) * this.scale;
        }

internal override double PreTransform(double x)
        {
            Debug.replacedert(x > 0, "Value should be positive.");

            if (x <= 0)
            {
                return 0;
            }

            return Math.Log(x);
        }

internal override void UpdateActualMaxMin()
        {
            if (this.PowerPadding)
            {
                double logBase = Math.Log(this.Base);
                var e0 = (int)Math.Floor(Math.Log(this.ActualMinimum) / logBase);
                var e1 = (int)Math.Ceiling(Math.Log(this.ActualMaximum) / logBase);
                if (!double.IsNaN(this.ActualMinimum))
                {
                    this.ActualMinimum = Math.Exp(e0 * logBase).RemoveNoiseFromDoubleMath();
                }

                if (!double.IsNaN(this.ActualMaximum))
                {
                    this.ActualMaximum = Math.Exp(e1 * logBase).RemoveNoiseFromDoubleMath();
                }
            }

            base.UpdateActualMaxMin();
        }

public override void Eval()
        {
            Value = Math.Log(Arg.Value);
        }

public override void Diff()
        {
            var arg = Arg.Value;
            Value = Math.Log(arg);
            Inputs.SetWeight(ArgIdx, 1 / arg);        
        }

[Fact]
        public void DiffTermPower()
        {
            var x = new Variable();
            var y = new Variable();
            var func = Power(x, y);

            var grad = func.Differentiate(Vec(x, y), NumVec(2, 3));

            replacedert.Equal(NumVec(12, 8 * Math.Log(2)), grad);
        }

[Fact]
        public void DiffTermPowerSingleVariable()
        {
            var x = new Variable();
            var func = Power(x, x);

            var grad = func.Differentiate(Vec(x), NumVec(2.5));
            var expectedGrad = NumVec(Math.Pow(2.5, 2.5) * (Math.Log(2.5) + 1));
            replacedert.Equal(expectedGrad[0], grad[0], 10);
        }

public override void Diff()
        {
            var baseVal = Base.Value;
            var expVal = Exponent.Value;

            Value = Math.Pow(baseVal, expVal);
            Inputs.SetWeight(BaseIdx, expVal * Math.Pow(baseVal, expVal - 1));
            Inputs.SetWeight(ExpIdx, Value * Math.Log(baseVal));
        }

private static bool PolarTransform(double x01, double y01, out double x, out double y)
        {
            var xn11 = 2.0 * x01 - 1.0;
            var yn11 = 2.0 * y01 - 1.0;
            var r2 = xn11 * xn11 + yn11 * yn11;
            if (r2 >= 1.0 || r2 == 0.0) {
                x = 0.0;
                y = 0.0;
                return false;
            }
            var factor = Math.Sqrt(-2.0 * Math.Log(r2) / r2);
            x = xn11 * factor;
            y = yn11 * factor;
            return true;
        }

public CoordinatePair LatLonToMeters(double lat, double lon)
        {
            CoordinatePair retval = new CoordinatePair();
            try
            {
                retval.X = lon * this.originShift / 180.0;
                retval.Y = Math.Log(Math.Tan((90 + lat) * Math.PI / 360.0)) / (Math.PI / 180.0);

                retval.Y *= this.originShift / 180.0;
                return retval;
            }
            catch (Exception ex)
            {
                throw ex;
            }
        }

public static double NextGaussian(this Random r, double mu = 0, double sigma = 1) {
            var u1 = r.NextDouble();
            var u2 = r.NextDouble();
            var rand_std_normal = Math.Sqrt(-2.0 * Math.Log(u1)) * Math.Sin(2.0 * Math.PI * u2);
            var rand_normal = mu + sigma * rand_std_normal;

            // Very rarely, it is possible to have underflow or an infinity if u1/u1 = 0.
            // In such a case, just return the mean.
            if (
                Double.IsNaN(rand_normal)
                || Double.IsInfinity(rand_normal)
                || Double.IsNegativeInfinity(rand_normal)
            ) {
                return mu;
            }

            return rand_normal;
        }

public static double Log(this double d)
        {
            return Math.Log(d);
        }

public static double NextGauss(this Random rand, double mean, double stdDev)
        {
            double u1 = 1.0 - rand.NextDouble();
            double u2 = 1.0 - rand.NextDouble();
            double randStdNormal = Math.Sqrt(-2.0 * Math.Log(u1)) * Math.Sin(2.0 * Math.PI * u2);
            return mean + stdDev * randStdNormal;
        }

[MethodImpl(MethodImplOptions.AggressiveInlining), IntrinsicFunction("logf")]
        public static float Logf(float f)
        {
            return (float)Math.Log((double)f);
        }

[MethodImpl(MethodImplOptions.AggressiveInlining), IntrinsicFunction("__logf")]
        public static float Logf(float f)
        {
            return (float)Math.Log((double)f);
        }

[IntrinsicFunction("logf")]
        public static float logf(float x)
        {
            return (float)Math.Log(x);
        }

public object Convert(object value, Type targetType, object parameter, CultureInfo culture)
        {
            if (value == null)
                return null;

            double number = System.Convert.ToDouble(value);
            int b = parameter != null ? System.Convert.ToInt32(parameter) : Base;

            if (number <= 0.5 && EmptyIfZero)
                return "";

            int n = number > 0.5 ? Math.Min((int)Math.Truncate(Math.Log(Math.Abs(number)) / Math.Log(b)), Prefixes.Length) : 0;
            double x = number / Math.Pow(b, n);
            return string.Format(
                Views.Resources.Strings.NumberToHumanReadable,
                n > 0 && Math.Abs(x) < 10 ?
                    (Math.Truncate(x * 10) / 10).ToString("N1") :
                     Math.Truncate(x).ToString(),
                Prefixes[n],
                Unit);
        }

int sample(Random random, float[] preds, double temperature=1.0) {
      // step 1: apply temperature to predictions, and normalize them to create a probability distribution
      float sum = 0;
      for (int i=0; i<preds.Length; i++) {
        var p = (float)Math.Exp((Math.Log(Math.Max(preds[i], 1e-10)) / temperature));
        sum += p;
        preds[i] = p;
      }
      for (int i = 0; i < preds.Length; i++) { preds[i] /= sum; }

      // step 2: draw a random sample from this distribution
      var d = random.NextDouble();
      sum = 0;
      for (int i=0; i<preds.Length; i++) {
        sum += preds[i];
        if ( d<sum ) { return i; }
      }
      return preds.Length - 1;
    }

public double GetCost(int tau0, int tau1, int tau2)
            {
                double sum = 0;
                int offset = tau1; // offset of partialSums'[i, tau1] in the single-dimenstional `partialSums` array
                int tauDiff = tau2 - tau1;
                for (int i = 0; i < k; i++)
                {
                    // actualSum is (count(data[j] < t) * 2 + count(data[j] == t) * 1) for j=tau1..tau2-1
                    int actualSum =
                        partialSums[offset + tauDiff] -
                        partialSums[offset]; // partialSums'[i, tau2] - partialSums'[i, tau1]

                    // We skip these two cases (correspond to fit = 0 or fit = 1) because of invalid Math.Log values
                    if (actualSum != 0 && actualSum != tauDiff * 2)
                    {
                        // Empirical CDF $\hat{F}_i(t)$ (Section 2.1 "Model" in [Haynes2017])
                        double fit = actualSum * 0.5 / tauDiff;

                        // Segment cost $\mathcal{L}_{np}$ (Section 2.2 "Nonparametric maximum likelihood" in [Haynes2017])
                        double lnp = tauDiff * (fit * Math.Log(fit) + (1 - fit) * Math.Log(1 - fit));
                        sum += lnp;
                    }

                    offset += n + 1;
                }

                double c = -Math.Log(2 * n - 1); // Constant from Lemma 3.1 in [Haynes2017]
                return 2.0 * c / k * sum; // See Section 3.1 "Discrete approximation" in [Haynes2017]
            }

public double Pdf(double x)
        {
            if (x < 0 || x > 1)
                return 0;

            if (x < 1e-9)
            {
                if (Alpha > 1)
                    return 0;
                if (Abs(Alpha - 1) < 1e-9)
                    return Beta;
                return double.PositiveInfinity;
            }

            if (x > 1 - 1e-9)
            {
                if (Beta > 1)
                    return 0;
                if (Abs(Beta - 1) < 1e-9)
                    return Alpha;
                return double.PositiveInfinity;
            }

            if (Alpha < 1e-9 || Beta < 1e-9)
                return 0;

            return Exp((Alpha - 1) * Log(x) + (Beta - 1) * Log(1 - x) - BetaFunction.CompleteLogValue(Alpha, Beta));
        }

public double Pdf(double x)
        {
            if (x < 1e-9)
                return 0;
            return Exp(-(Log(x) - Mean).Sqr() / (2 * StandardDeviation.Sqr())) / (x * StandardDeviation * Constants.Sqrt2Pi);
        }

public double Cdf(double x)
        {
            if (x < 1e-9)
                return 0;
            return 0.5 * (1 + ErrorFunction.Value((Log(x) - Mean) / (Constants.Sqrt2 * StandardDeviation)));
        }

public override double Next()
            {
                double u = 0, v = 0;
                while (u < 1e-100)
                {
                    u = Random.NextDouble();
                    v = Random.NextDouble();
                }
                double stdDevFactor = Sqrt(-2.0 * Log(u)) * Sin(2.0 * PI * v);
                return distribution.Mean + distribution.StandardDeviation * stdDevFactor;
            }

public static double IncompleteLogValue(double a, double b, double x)
        {
            return a * Log(x) - Log(a) + Log(HypergeometricFunction.Value(a, 1 - b, a + 1, x, (int) Round(b)));
        }

public static double RegularizedIncompleteValue(double a, double b, double x)
        {
            // The implementation is inspired by "Incomplete Beta Function in C" (Lewis Van Winkle, 2017)
            // See https://codeplea.com/incomplete-beta-function-c for details
            //
            // We calculate the regularized incomplete beta function using a continued fraction (https://dlmf.nist.gov/8.17#v):
            //   Ix(a, b) = x^a * (1-x)^b / (a*B(a, b)) * 1 / (1 + d[1] / (1 + d[2] / (1 + d[3] / (...))))
            // where
            //   d[2m]   = m(b-m)x / (a+2m-1)(a+2m)
            //   d[2m+1] = -(a+m)(a+b+m)x / (a+2m)(a+2m+1)
            //
            // The approximated value of the continued fraction is calculated using the Lentz's algorithm

            replacedertion.NonNegative(nameof(a), a);
            replacedertion.NonNegative(nameof(b), b);

            const double eps = 1e-8;
            if (x < eps)
                return 0;
            if (x > 1 - eps)
                return 1;
            if (a < eps && b < eps)
                return 0.5;
            if (a < eps)
                return 1;
            if (b < eps)
                return 0;

            // According to https://dlmf.nist.gov/8.17#v, the continued fraction converges rapidly for x<(a+1)/(a+b+2)
            // If x>=(a+1)/(a+b+2), we use Ix(a, b) = I{1-x}(b, a)
            if (x > (a + 1) / (a + b + 2))
                return 1.0 - RegularizedIncompleteValue(b, a, 1 - x);

            // We use the Lentz's algorithm to calculate the continued fraction
            //   f = 1 + d[1] / (1 + d[2] / (1 + d[3] / (...)))
            // The implementation is based on the following formulas:
            //   u[0] = 1, v[0] = 0, f[0] = 1
            //   u[i] = 1 + d[i] / u[i - 1]
            //   v[i] = 1 / (1 + d[i] * v[i - 1])
            //   f[i] = f[i - 1] * u[i] * v[i]

            const int maxIterationCount = 300;
            static double Normalize(double z) => Abs(z) < 1e-30 ? 1e-30 : z; // Normalization prevents getting zero values

            double u = 1, v = 0, f = 1;
            for (int i = 0; i <= maxIterationCount; i++)
            {
                double d; // d[i]
                int m = i / 2;
                if (i == 0)
                    d = 1.0; // d[0]
                else if (i % 2 == 0)
                    d = m * (b - m) * x / ((a + 2 * m - 1) * (a + 2 * m)); // d[2m]
                else
                    d = -((a + m) * (a + b + m) * x) / ((a + 2.0 * m) * (a + 2.0 * m + 1)); // d[2m+1]

                u = Normalize(1 + d / u);
                v = 1 / Normalize(1 + d * v);
                double uv = u * v;
                f *= uv;

                if (Abs(uv - 1) < eps)
                    break;
            }

            // Ix(a, b) = x^a * (1-x)^b / (a*B(a, b)) * 1 / (1 + d[1] / (1 + d[2] / (1 + d[3] / (...))))
            return Exp(Log(x) * a + Log(1.0 - x) * b - CompleteLogValue(a, b)) / a * (f - 1);
        }

public static double RegularizedIncompleteInverseValue(double a, double b, double p)
        {
            // The implementation is based on "Incomplete Beta Function" from "Numerical Recipes", 3rd edition, page 273

            replacedertion.NonNegative(nameof(a), a);
            replacedertion.NonNegative(nameof(b), b);

            if (p <= 0)
                return 0;
            if (p >= 1)
                return 1;

            const double eps = 1e-8;
            double t, u, x, w;

            if (a >= 1 && b >= 1)
            {
                double pp = p < 0.5 ? p : 1.0 - p;
                t = Sqrt(-2.0 * Log(pp));
                x = (2.30753 + t * 0.27061) / (1.0 + t * (0.99229 + t * 0.04481)) - t;
                if (p < 0.5)
                    x = -x;
                double al = (x.Sqr() - 3.0) / 6.0;
                double h = 2.0 / (1.0 / (2.0 * a - 1.0) + 1.0 / (2.0 * b - 1.0));
                w = x * Sqrt(al + h) / h - (1.0 / (2.0 * b - 1) - 1.0 / (2.0 * a - 1.0)) * (al + 5.0 / 6.0 - 2.0 / (3.0 * h));
                x = a / (a + b * Exp(2.0 * w));
            }
            else
            {
                double lna = Log(a / (a + b));
                double lnb = Log(b / (a + b));
                t = Exp(a * lna) / a;
                u = Exp(b * lnb) / b;
                w = t + u;
                x = p < t / w
                    ? Pow(a * w * p, 1.0 / a)
                    : 1.0 - Pow(b * w * (1.0 - p), 1.0 / b);
            }

            double afac = -GammaFunction.LogValue(a) - GammaFunction.LogValue(b) + GammaFunction.LogValue(a + b);
            for (int iteration = 0; iteration < 10; iteration++)
            {
                if (x < eps || x > 1.0 - eps)
                    return x; // a or b are too small for accurate calculations

                double error = RegularizedIncompleteValue(a, b, x) - p;
                t = Exp((a - 1) * Log(x) + (b - 1) * Log(1.0 - x) + afac);
                u = error / t;
                t = u / (1.0 - 0.5 * Min(1.0, u * ((a - 1) / x - (b - 1) / (1.0 - x)))); // Halley's method
                x -= t;
                if (x <= 0.0)
                    x = 0.5 * (x + t);
                if (x >= 1.0)
                    x = 0.5 * (x + t + 1.0); // Bisect if x tries to go negative or > 1
                if (Abs(t) < eps * x && iteration > 0)
                    break;
            }

            return x;
        }

public static double InverseValue(double p)
        {
            replacedertion.InRangeExclusive(nameof(p), p, -1, 1);

            p = 1 - p;
            double pp = p < 1.0 ? p : 2 - p;
            double t = Sqrt(-2 * Log(pp / 2));
            double x = -0.70711 * ((2.30753 + t * 0.27061) / (1 + t * (0.99229 + t * 0.04481)) - t);
            for (int i = 0; i < 2; i++)
            {
                double err = (1 - Value(x)) - pp;
                x += err / (1.12837916709551257 * Exp(-x.Sqr()) - x * err);
            }
            return p < 1.0 ? x : -x;
        }

public static double LogValue(int n)
        {
            replacedertion.NonNegative(nameof(n), n);

            if (n <= 20)
                return Math.Log(Value(n));

            return GammaFunction.LogValue(n + 1);
        }

public static double LogValue(double x)
        {
            if (x < 1e-5)
                throw new ArgumentOutOfRangeException(nameof(x), "x should be positive");

            if (x < 3)
                return Log(Value(x));

            return StirlingApproximationLog(x);
        }

private static double StirlingApproximationLog(double x)
        {
            return x * Log(x) - x + Log(2 * PI / x) / 2 + GetSeriesValue(x);
        }

public static ISequence CreateFromHalfLife(int halfLife) => new ExponentialDecaySequence(1, Math.Log(2) / halfLife);

[Fact]
        public void BetaCompleteLogValue()
        {
            var comparer = new AbsoluteEqualityComparer(0.0000001);
            for (int a = 1; a <= 20; a++)
            for (int b = 1; b <= 20; b++)
            {
                double actual = BetaFunction.CompleteLogValue(a, b);
                double expected = Math.Log(Factorial(a - 1) * Factorial(b - 1) / Factorial(a + b - 1));
                replacedert.Equal(expected, actual, comparer);
            }
        }

[Fact]
        public void GammaFunctionLogValue()
        {
            var comparer = new AbsoluteEqualityComparer(0.0001);
            for (int i = 0; i < knownX.Length; i ++)
            {
                double expected = Math.Log(knownY[i]);
                double actual = GammaFunction.LogValue(knownX[i]);
                
                output.WriteLine($"X = {knownX[i]}");
                output.WriteLine($"Expected = {expected}");
                output.WriteLine($"Actual   = {actual}");
                output.WriteLine("");
                
                replacedert.Equal(expected, actual, comparer);
            }
        }

private void Regression( RegressionType regressionType, double [][] inputValues, out double [][] outputValues, int polynomialDegree, int forecastingPeriod )
		{
			if( regressionType == RegressionType.Exponential )
			{
				double [] oldYValues = new double[ inputValues[1].Length ];
				for( int index = 0; index < inputValues[1].Length; index++ )
				{
					oldYValues[ index ] = inputValues[1][index];
					if( inputValues[1][index] <= 0 )
					{
                        throw new InvalidOperationException(SR.ExceptionForecastingExponentialRegressionHasZeroYValues);
					}
					inputValues[1][index] = Math.Log( inputValues[1][index] );
				}

				

				PolynomialRegression( regressionType, inputValues, out outputValues, 2, forecastingPeriod, 0 );

				inputValues[1] = oldYValues;
			}
			else if( regressionType == RegressionType.Logarithmic )
			{
				double interval;
				double first = inputValues[0][0];

				// Find Interval for X values
				interval = Math.Abs( inputValues[0][0] - inputValues[0][inputValues[0].Length - 1] ) / ( inputValues[0].Length - 1 );
			
				if( interval <= 0 )
					interval = 1;

				for( int index = 0; index < inputValues[0].Length; index++ )
				{
					inputValues[0][index] = Math.Log( inputValues[0][index] );
				}

				PolynomialRegression( regressionType, inputValues, out outputValues, 2, forecastingPeriod, interval );

				// Create Y values based on approximation.
				for( int i = 0; i < outputValues[0].Length; i++ )
				{
					// Set X value
					outputValues[0][i] = first + i * interval;
				}
			}
			else if( regressionType == RegressionType.Power )
			{
				double [] oldYValues = new double[ inputValues[1].Length ];
				double interval;
				double first = inputValues[0][0];

				for( int index = 0; index < inputValues[1].Length; index++ )
				{
					oldYValues[ index ] = inputValues[1][index];
					if( inputValues[1][index] <= 0 )
					{
                        throw new InvalidOperationException(SR.ExceptionForecastingPowerRegressionHasZeroYValues);
					}
				}

				// Find Interval for X values
				interval = Math.Abs( inputValues[0][0] - inputValues[0][inputValues[0].Length - 1] ) / ( inputValues[0].Length - 1 );
			
				if( interval <= 0 )
					interval = 1;

				PolynomialRegression( regressionType, inputValues, out outputValues, 2, forecastingPeriod, interval );

				inputValues[1] = oldYValues;

				// Create Y values based on approximation.
				for( int i = 0; i < outputValues[0].Length; i++ )
				{
					// Set X value
					outputValues[0][i] = first + i * interval;
				}
			}
			else
			{
				PolynomialRegression( regressionType, inputValues, out outputValues, polynomialDegree, forecastingPeriod, 0 );
			}
		}