Compute turbulence using Perlin noise.

@param x       the x value
@param y       the y value
@param octaves number of octaves of turbulence
@return turbulence value at (x,y)
192234
public float turbulence3(float x, float y, float z, float octaves) {
        float t = 0.0f;

        for (float f = 1.0f; f <= octaves; f *= 2)
            t += Math.abs(noise3(f * x, f * y, f * z)) / f;
        return t;
    }
123880
public Matrix4x3d scaleXY(double x, double y, Matrix4x3d dest) {
        return scale(x, y, 1.0, dest);
    }
1233564
public static double logSubtract(double x, double y) {
        if (FastMath.useLogAddTable) {
            return SmoothedLogAddTable.logSubtract(x,y);
        } else {
            return FastMath.logSubtractExact(x,y);
        }
    }
784758
function(x,y){
        var _lsw = (x&0xFFFF)+(y&0xFFFF),
            _msw = (x>>16)+(y>>16)+(_lsw>>16);
        return (_msw<<16)|(_lsw&0xFFFF);
    }
332469
public static double beta(double x, double y) {
        return Math.exp(Gamma.lgamma(x) + Gamma.lgamma(y) - Gamma.lgamma(x + y));
    }
5393
public static double uncenteredCorrelation(NumberVector x, NumberVector y) {
    final int xdim = x.getDimensionality(), ydim = y.getDimensionality();
    if(xdim != ydim) {
      throw new IllegalArgumentException("Invalid arguments: number vectors differ in dimensionality.");
    }
    double sumXX = 0., sumYY = 0., sumXY = 0.;
    for(int i = 0; i < xdim; i++) {
      final double xv = x.doubleValue(i), yv = y.doubleValue(i);
      sumXX += xv * xv;
      sumYY += yv * yv;
      sumXY += xv * yv;
    }
    // One or both series were constant:
    if(!(sumXX > 0. && sumYY > 0.)) {
      return (sumXX == sumYY) ? 1. : 0.;
    }
    return sumXY / FastMath.sqrt(sumXX * sumYY);
  }
561426
protected static void computeAlphas( InterleavedF64 yf , InterleavedF64 kf , double lambda ,
										 InterleavedF64 alphaf ) {

		for( int y = 0; y < kf.height; y++ ) {

			int index = yf.startIndex + y*yf.stride;

			for( int x = 0; x < kf.width; x++, index += 2 ) {
				double a = yf.data[index];
				double b = yf.data[index+1];

				double c = kf.data[index] + lambda;
				double d = kf.data[index+1];

				double bottom = c*c + d*d;

				alphaf.data[index] = (a*c + b*d)/bottom;
				alphaf.data[index+1] = (b*c - a*d)/bottom;
			}
		}
	}
1313039
public static FloatMatrix pairwiseSquaredDistances(FloatMatrix X, FloatMatrix Y) {
		if (X.rows != Y.rows)
			throw new IllegalArgumentException(
					"Matrices must have same number of rows");
	
		FloatMatrix XX = X.mul(X).columnSums();
		FloatMatrix YY = Y.mul(Y).columnSums();
	
		FloatMatrix Z = X.transpose().mmul(Y);
		Z.muli(-2.0f); //Z.print();
		Z.addiColumnVector(XX);
		Z.addiRowVector(YY);
	
		return Z;
	}
778092
public static double atan2(double y, double x) {
        if (x != x || y != y) {
            return Double.NaN;
        }

        if (y == 0) {
            final double result = x * y;
            final double invx = 1d / x;
            final double invy = 1d / y;

            if (invx == 0) { // X is infinite
                if (x > 0) {
                    return y; // return +/- 0.0
                } else {
                    return copySign(Math.PI, y);
                }
            }

            if (x < 0 || invx < 0) {
                if (y < 0 || invy < 0) {
                    return -Math.PI;
                } else {
                    return Math.PI;
                }
            } else {
                return result;
            }
        }

        // y cannot now be zero

        if (y == Double.POSITIVE_INFINITY) {
            if (x == Double.POSITIVE_INFINITY) {
                return Math.PI * F_1_4;
            }

            if (x == Double.NEGATIVE_INFINITY) {
                return Math.PI * F_3_4;
            }

            return Math.PI * F_1_2;
        }

        if (y == Double.NEGATIVE_INFINITY) {
            if (x == Double.POSITIVE_INFINITY) {
                return -Math.PI * F_1_4;
            }

            if (x == Double.NEGATIVE_INFINITY) {
                return -Math.PI * F_3_4;
            }

            return -Math.PI * F_1_2;
        }

        if (x == Double.POSITIVE_INFINITY) {
            if (y > 0 || 1 / y > 0) {
                return 0d;
            }

            if (y < 0 || 1 / y < 0) {
                return -0d;
            }
        }

        if (x == Double.NEGATIVE_INFINITY)
        {
            if (y > 0.0 || 1 / y > 0.0) {
                return Math.PI;
            }

            if (y < 0 || 1 / y < 0) {
                return -Math.PI;
            }
        }

        // Neither y nor x can be infinite or NAN here

        if (x == 0) {
            if (y > 0 || 1 / y > 0) {
                return Math.PI * F_1_2;
            }

            if (y < 0 || 1 / y < 0) {
                return -Math.PI * F_1_2;
            }
        }

        // Compute ratio r = y/x
        final double r = y / x;
        if (Double.isInfinite(r)) { // bypass calculations that can create NaN
            return atan(r, 0, x < 0);
        }

        double ra = doubleHighPart(r);
        double rb = r - ra;

        // Split x
        final double xa = doubleHighPart(x);
        final double xb = x - xa;

        rb += (y - ra * xa - ra * xb - rb * xa - rb * xb) / x;

        final double temp = ra + rb;
        rb = -(temp - ra - rb);
        ra = temp;

        if (ra == 0) { // Fix up the sign so atan works correctly
            ra = copySign(0d, y);
        }

        // Call atan
        final double result = atan(ra, rb, x < 0);

        return result;
    }
660752
public static void fft(boolean forward, long m, double[] x , double[] y)
    {
        int n, i, i1, j, k, i2, l, l1, l2;
        double c1, c2, tx, ty, t1, t2, u1, u2, z;

        /* Calculate the number of points */
        n = 1;
        for (i = 0; i < m; i++)
        {
            n *= 2;
        }

        /* Do the bit reversal */
        i2 = n >> 1;
        j = 0;
        for (i = 0; i < n - 1; i++)
        {
            if (i < j)
            {
                tx = x[i];
                ty = y[i];
                x[i] = x[j];
                y[i] = y[j];
                x[j] = tx;
                y[j] = ty;
            }
            k = i2;
            while (k <= j)
            {
                j -= k;
                k >>= 1;
            }
            j += k;
        }

        /* Compute the FFT */
        c1 = -1.0;
        c2 = 0.0;
        l2 = 1;
        for (l = 0; l < m; l++)
        {
            l1 = l2;
            l2 <<= 1;
            u1 = 1.0;
            u2 = 0.0;
            for (j = 0; j < l1; j++)
            {
                for (i = j; i < n; i += l2)
                {
                    i1 = i + l1;
                    t1 = u1 * x[i1] - u2 * y[i1];
                    t2 = u1 * y[i1] + u2 * x[i1];
                    x[i1] = x[i] - t1;
                    y[i1] = y[i] - t2;
                    x[i] += t1;
                    y[i] += t2;
                }
                z = u1 * c1 - u2 * c2;
                u2 = u1 * c2 + u2 * c1;
                u1 = z;
            }
            c2 = Math.sqrt((1.0 - c1) / 2.0);
            if (forward)
            {
                c2 = -c2;
            }
            c1 = Math.sqrt((1.0 + c1) / 2.0);
        }

        /* Scaling for forward transform */
        if (forward)
        {
            for (i = 0; i < n; i++)
            {
                x[i] /= n;
                y[i] /= n;
            }
        }
    }
704700
def xcorr(x, y, maxlags):
    
    xlen = len(x)
    ylen = len(y)
    assert xlen == ylen

    c = np.correlate(x, y, mode=2)

    # normalize
    c /= np.sqrt(np.dot(x, x) * np.dot(y, y))

    lags = np.arange(-maxlags, maxlags + 1)
    c = c[xlen - 1 - maxlags:xlen + maxlags]

    return c
421281
private float calculateScale(int x, int y) {
		final float centerX = getWidth() / 2f;
		final float centerY = getHeight() / 2f;
		final float maxDistance = (float) Math.sqrt(centerX * centerX + centerY * centerY);

		final float deltaX = centerX - x;
		final float deltaY = centerY - y;
		final float distance = (float) Math.sqrt(deltaX * deltaX + deltaY * deltaY);
		final float scale = 0.5f + (distance / maxDistance) * 0.5f;
		return scale;
	}
857653
public static double coefficient(double[] x, double[] y) {
    final int xdim = x.length;
    final int ydim = y.length;
    if(xdim != ydim) {
      throw new IllegalArgumentException("Invalid arguments: arrays differ in length.");
    }
    if(xdim == 0) {
      throw new IllegalArgumentException("Empty vector.");
    }
    // Inlined computation of Pearson correlation, to avoid allocating objects!
    // This is a numerically stabilized version, avoiding sum-of-squares.
    double sumXX = 0., sumYY = 0., sumXY = 0.;
    double sumX = x[0], sumY = y[0];
    int i = 1;
    while(i < xdim) {
      final double xv = x[i], yv = y[i];
      // Delta to previous mean
      final double deltaX = xv * i - sumX;
      final double deltaY = yv * i - sumY;
      // Increment count first
      final double oldi = i; // Convert to double!
      ++i;
      final double f = 1. / (i * oldi);
      // Update
      sumXX += f * deltaX * deltaX;
      sumYY += f * deltaY * deltaY;
      // should equal deltaY * neltaX!
      sumXY += f * deltaX * deltaY;
      // Update sums
      sumX += xv;
      sumY += yv;
    }
    // One or both series were constant:
    if(!(sumXX > 0. && sumYY > 0.)) {
      return (sumXX == sumYY) ? 1. : 0.;
    }
    return sumXY / FastMath.sqrt(sumXX * sumYY);
  }
9859
public static double dot(DoubleMatrix x, DoubleMatrix y) {
		//return NativeBlas.ddot(x.length, x.data, 0, 1, y.data, 0, 1);
		return JavaBlas.rdot(x.length, x.data, 0, 1, y.data, 0, 1);
	}
1103538
private int computeOffset(int x, int y) {
    if(y > x) {
      return computeOffset(y, x);
    }
    return ((x * (x + 1)) >> 1) + y;
  }
156706
protected DoubleMatrix2D[] getHessFi(DoubleMatrix1D X){
		DoubleMatrix2D[] ret = new DoubleMatrix2D[request.getFi().length];
		double[] x = X.toArray();
		for(int i=0; i<request.getFi().length; i++){
			double[][] hess = request.getFi()[i].hessian(x);
			if(hess == FunctionsUtils.ZEROES_2D_ARRAY_PLACEHOLDER){
				ret[i] = FunctionsUtils.ZEROES_MATRIX_PLACEHOLDER;
			}else{
				ret[i] = F2.make(hess);
			}
		}
		return ret;
	}
1004396
def _nfw_func(self, x):
        

        c = 0.000001

        if isinstance(x, np.ndarray):
            x[np.where(x<c)] = c
            nfwvals = np.ones_like(x)
            inds1 = np.where(x < 1)
            inds2 = np.where(x > 1)

            nfwvals[inds1] = (1 - x[inds1] ** 2) ** -.5 * np.arctanh((1 - x[inds1] ** 2) ** .5)
            nfwvals[inds2] = (x[inds2] ** 2 - 1) ** -.5 * np.arctan((x[inds2] ** 2 - 1) ** .5)

            return nfwvals

        elif isinstance(x, float) or isinstance(x, int):
            x = max(x, c)
            if x == 1:
                return 1
            if x < 1:
                return (1 - x ** 2) ** -.5 * np.arctanh((1 - x ** 2) ** .5)
            else:
                return (x ** 2 - 1) ** -.5 * np.arctan((x ** 2 - 1) ** .5)