/* interpolated 1./sqrt(p) where .5 <= p < 1.
516363
static float invsqlook(float a) {
        double d = a * (2.f * (float) INVSQ_LOOKUP_SZ) - (float) INVSQ_LOOKUP_SZ;
        int i = (int) d;
        return INVSQ_LOOKUP[i] + ((float) (d - i)) * (INVSQ_LOOKUP[i + 1] - INVSQ_LOOKUP[i]);
    }
527296
public static double inverfc(double p) {
        double x, err, t, pp;
        if (p >= 2.0) {
            return -100.;
        }
        if (p <= 0.0) {
            return 100.;
        }
        pp = (p < 1.0) ? p : 2. - p;
        t = Math.sqrt(-2. * Math.log(pp / 2.));
        x = -0.70711 * ((2.30753 + t * 0.27061) / (1. + t * (0.99229 + t * 0.04481)) - t);
        for (int j = 0; j < 2; j++) {
            err = erfc(x) - pp;
            x += err / (1.12837916709551257 * Math.exp(-x * x) - x * err);
        }
        return (p < 1.0 ? x : -x);
    }
926031
function(ratio, opt) {

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

            opt = opt || {};
            var precision = (opt.precision === undefined) ? this.PRECISION : opt.precision;
            var subdivisions = (opt.subdivisions === undefined) ? this.getSubdivisions({ precision: precision }) : opt.subdivisions;
            var localOpt = { precision: precision, subdivisions: subdivisions };

            var curveLength = this.length(localOpt);
            var length = curveLength * ratio;

            return this.tAtLength(length, localOpt);
        }
479584
def sqrt(self, val, flag):
        
        if val.iszero():
            return val
        sw = self.p % 8
        if sw == 3 or sw == 7:
            res = val ** ((self.p + 1) / 4)
        elif sw == 5:
            x = val ** ((self.p + 1) / 4)
            if x == 1:
                res = val ** ((self.p + 3) / 8)
            else:
                res = (4 * val) ** ((self.p - 5) / 8) * 2 * val
        else:
            raise Exception("modsqrt non supported for (p%8)==1")

        if res.value % 2 == flag:
            return res
        else:
            return -res
1202093
private double rdist(double[] x1, double[] x2) {
        double d = 0.0;
        for (int i = 0; i < x1.length; i++) {
            double t = x1[i] - x2[i];
            d += t * t;
        }
        return Math.sqrt(d);
    }
1052151
static int isqrt(long n) {
        int x = (int) ((Math.sqrt(8.0 * n + 1.0) - 1.0) / 2.0);
        if (x * (x + 1L) / 2 > n)
            x--;
        return x;
    }
225817
@Override
    public double inverseCumulativeProbability(final double p) throws OutOfRangeException {
        /*
         * IMPLEMENTATION NOTES
         * --------------------
         * Where applicable, use is made of the one-sided Chebyshev inequality
         * to bracket the root. This inequality states that
         * P(X - mu >= k * sig) <= 1 / (1 + k^2),
         * mu: mean, sig: standard deviation. Equivalently
         * 1 - P(X < mu + k * sig) <= 1 / (1 + k^2),
         * F(mu + k * sig) >= k^2 / (1 + k^2).
         *
         * For k = sqrt(p / (1 - p)), we find
         * F(mu + k * sig) >= p,
         * and (mu + k * sig) is an upper-bound for the root.
         *
         * Then, introducing Y = -X, mean(Y) = -mu, sd(Y) = sig, and
         * P(Y >= -mu + k * sig) <= 1 / (1 + k^2),
         * P(-X >= -mu + k * sig) <= 1 / (1 + k^2),
         * P(X <= mu - k * sig) <= 1 / (1 + k^2),
         * F(mu - k * sig) <= 1 / (1 + k^2).
         *
         * For k = sqrt((1 - p) / p), we find
         * F(mu - k * sig) <= p,
         * and (mu - k * sig) is a lower-bound for the root.
         *
         * In cases where the Chebyshev inequality does not apply, geometric
         * progressions 1, 2, 4, ... and -1, -2, -4, ... are used to bracket
         * the root.
         */
        if (p < 0.0 || p > 1.0) {
            throw new OutOfRangeException(p, 0, 1);
        }

        double lowerBound = getSupportLowerBound();
        if (p == 0.0) {
            return lowerBound;
        }

        double upperBound = getSupportUpperBound();
        if (p == 1.0) {
            return upperBound;
        }

        final double mu = getNumericalMean();
        final double sig = FastMath.sqrt(getNumericalVariance());
        final boolean chebyshevApplies;
        chebyshevApplies = !(Double.isInfinite(mu) || Double.isNaN(mu) || Double.isInfinite(sig) || Double.isNaN(sig));

        if (lowerBound == Double.NEGATIVE_INFINITY) {
            if (chebyshevApplies) {
                lowerBound = mu - sig * FastMath.sqrt((1. - p) / p);
            } else {
                lowerBound = -1.0;
                while (cumulativeProbability(lowerBound) >= p) {
                    lowerBound *= 2.0;
                }
            }
        }

        if (upperBound == Double.POSITIVE_INFINITY) {
            if (chebyshevApplies) {
                upperBound = mu + sig * FastMath.sqrt(p / (1. - p));
            } else {
                upperBound = 1.0;
                while (cumulativeProbability(upperBound) < p) {
                    upperBound *= 2.0;
                }
            }
        }

        final UnivariateFunction toSolve = new UnivariateFunction() {

            public double value(final double x) {
                return cumulativeProbability(x) - p;
            }
        };

        double x = UnivariateSolverUtils.solve(toSolve, lowerBound, upperBound, getSolverAbsoluteAccuracy());

        if (!isSupportConnected()) {
            /* Test for plateau. */
            final double dx = getSolverAbsoluteAccuracy();
            if (x - dx >= getSupportLowerBound()) {
                double px = cumulativeProbability(x);
                if (cumulativeProbability(x - dx) == px) {
                    upperBound = x;
                    while (upperBound - lowerBound > dx) {
                        final double midPoint = 0.5 * (lowerBound + upperBound);
                        if (cumulativeProbability(midPoint) < px) {
                            lowerBound = midPoint;
                        } else {
                            upperBound = midPoint;
                        }
                    }
                    return upperBound;
                }
            }
        }
        return x;
    }
549536
protected function interpolate($x, $y, $nw, $ne, $sw, $se)
    {
        list($r0, $g0, $b0) = $this->getRGB($nw);
        list($r1, $g1, $b1) = $this->getRGB($ne);
        list($r2, $g2, $b2) = $this->getRGB($sw);
        list($r3, $g3, $b3) = $this->getRGB($se);

        $cx = 1.0 - $x;
        $cy = 1.0 - $y;

        $m0 = $cx * $r0 + $x * $r1;
        $m1 = $cx * $r2 + $x * $r3;
        $r  = (int) ($cy * $m0 + $y * $m1);

        $m0 = $cx * $g0 + $x * $g1;
        $m1 = $cx * $g2 + $x * $g3;
        $g  = (int) ($cy * $m0 + $y * $m1);

        $m0 = $cx * $b0 + $x * $b1;
        $m1 = $cx * $b2 + $x * $b3;
        $b  = (int) ($cy * $m0 + $y * $m1);

        return ($r << 16) | ($g << 8) | $b;
    }
673953
@Override
    public double quantile(double p) {
        if (p < 0.0 || p > 1.0) {
            throw new IllegalArgumentException("Invalid p: " + p);
        }

        if (sigma == 0.0) {
            if (p < 1.0) {
                return mu - 1E-10;
            } else {
                return mu;
            }
        }
        return -1.41421356237309505 * sigma * Erf.inverfc(2.0 * p) + mu;
    }
507603
public static double interpolate(Coordinate firstCoordinate, Coordinate lastCoordinate, Coordinate toBeInterpolated) {
        if (Double.isNaN(firstCoordinate.z)) {
            return Double.NaN;
        }
        if (Double.isNaN(lastCoordinate.z)) {
            return Double.NaN;
        }
        return firstCoordinate.z + (lastCoordinate.z - firstCoordinate.z) * firstCoordinate.distance(toBeInterpolated)
                / (firstCoordinate.distance(toBeInterpolated) + toBeInterpolated.distance(lastCoordinate));
    }
732782
function gauss(x, s) {
    // 2.5066282746310005 is sqrt(2*pi)
    return Math.exp(-0.5 * x*x / (s*s)) / (s * 2.5066282746310005);
}
1161757
def rec_sqrt(x, context=None):
    
    return _apply_function_in_current_context(
        BigFloat,
        mpfr.mpfr_rec_sqrt,
        (BigFloat._implicit_convert(x),),
        context,
    )
1226304
def interpolate(p, axis_values, pixelgrid, order=1, mode='constant', cval=0.0):
    
    # convert requested parameter combination into a coordinate
    p_ = np.array([np.searchsorted(av_,val) for av_, val in zip(axis_values,p)])
    lowervals_stepsize = np.array([[av_[p__-1], av_[p__]-av_[p__-1]] \
                         for av_, p__ in zip(axis_values,p_)])
    p_coord = (p-lowervals_stepsize[:,0])/lowervals_stepsize[:,1] + p_-1
    
    # interpolate
    if order > 1:
        prefilter = False
    else:
        prefilter = False

    return [ndimage.map_coordinates(pixelgrid[...,i],p_coord, order=order,
                                    prefilter=prefilter, mode=mode, cval=cval) \
                for i in range(np.shape(pixelgrid)[-1])]
696928
function sqrt( out, x ) {
	var len = x.length,
		i;
	if ( out.length !== len ) {
		throw new Error( 'sqrt()::invalid input arguments. Input and output matrices must be the same length.' );
	}
	for ( i = 0; i < len; i++ ) {
		out.data[ i ] = SQRT( x.data[ i ] );
	}
	return out;
}
403526
def sqrt(x, context=None):
    
    return _apply_function_in_current_context(
        BigFloat,
        mpfr.mpfr_sqrt,
        (BigFloat._implicit_convert(x),),
        context,
    )
209899
private double handleCase1() {

		double stpc = SearchInterpolate.cubic2(fx, gx, stx, fp, gp, stp);
		double stpq = SearchInterpolate.quadratic(fx,gx,stx,fp,stp);

		// If the cubic step is closer to stx than the quadratic step take that
		bracket = true;
		if( Math.abs(stpc-stx) < Math.abs(stpq-stx)) {
			return stpc;
		} else {
			// otherwise go with the average of the twp
			return stpc + (stpq-stpc)/2.0;
		}
	}
1178597
function clamp$1(v, min, max) {
    if (min === void 0) { min = 0; }
    if (max === void 0) { max = 100; }
    return Math.max(min, Math.min(max, v));
}