Calculate kronecker symbol a|p.  Generalization of legendre and jacobi.

@param a an integer
@param p the modulus
@return a|p
19718
public static int kronecker(BigInteger a, BigInteger p) {
    return INSTANCE.get().kroneckerImpl(a, p);
  }
177530
function(p) {

            var refPointDelta = 30;
            var x0 = p.x;
            var y0 = p.y;
            var a = this.a;
            var b = this.b;
            var center = this.bbox().center();
            var m = center.x;
            var n = center.y;

            var q1 = x0 > center.x + a / 2;
            var q3 = x0 < center.x - a / 2;

            var y, x;
            if (q1 || q3) {
                y = x0 > center.x ? y0 - refPointDelta : y0 + refPointDelta;
                x = (a * a / (x0 - m)) - (a * a * (y0 - n) * (y - n)) / (b * b * (x0 - m)) + m;

            } else {
                x = y0 > center.y ? x0 + refPointDelta : x0 - refPointDelta;
                y = ( b * b / (y0 - n)) - (b * b * (x0 - m) * (x - m)) / (a * a * (y0 - n)) + n;
            }

            return (new Point(x, y)).theta(p);

        }
591041
public static double klDivergence(double[] p, double[] q) {
        if (p.length != q.length) {
            throw new IllegalStateException("The length of p and q must be the same.");
        }
        double delta = 1e-8;
        if (!Multinomials.isMultinomial(p, delta)) {
            throw new IllegalStateException("p is not a multinomial");
        }
        if (!Multinomials.isMultinomial(q, delta)) {
            throw new IllegalStateException("q is not a multinomial");
        }
        
        double kl = 0.0;
        for (int i=0; i<p.length; i++) {
            if (p[i] == 0.0 || q[i] == 0.0) {
                continue;
            }
            kl += p[i] * FastMath.log(p[i] / q[i]);
        }
        return kl;
    }
71420
static double expint(int p, final double result[]) {
        //double x = M_E;
        final double xs[] = new double[2];
        final double as[] = new double[2];
        final double ys[] = new double[2];
        //split(x, xs);
        //xs[1] = (double)(2.7182818284590452353602874713526625L - xs[0]);
        //xs[0] = 2.71827697753906250000;
        //xs[1] = 4.85091998273542816811e-06;
        //xs[0] = Double.longBitsToDouble(0x4005bf0800000000L);
        //xs[1] = Double.longBitsToDouble(0x3ed458a2bb4a9b00L);

        /* E */
        xs[0] = 2.718281828459045;
        xs[1] = 1.4456468917292502E-16;

        split(1.0, ys);

        while (p > 0) {
            if ((p & 1) != 0) {
                quadMult(ys, xs, as);
                ys[0] = as[0]; ys[1] = as[1];
            }

            quadMult(xs, xs, as);
            xs[0] = as[0]; xs[1] = as[1];

            p >>= 1;
        }

        if (result != null) {
            result[0] = ys[0];
            result[1] = ys[1];

            resplit(result);
        }

        return ys[0] + ys[1];
    }
1050741
public static double inverseRegularizedIncompleteGamma(double a, double p) {
        if (a <= 0.0) {
            throw new IllegalArgumentException("a must be pos in invgammap");
        }

        final double EPS = 1.0E-8;

        double x, err, t, u, pp;
        double lna1 = 0.0;
        double afac = 0.0;
        double a1 = a - 1;
        double gln = lgamma(a);
        if (p >= 1.) {
            return Math.max(100., a + 100. * Math.sqrt(a));
        }

        if (p <= 0.0) {
            return 0.0;
        }

        if (a > 1.0) {
            lna1 = Math.log(a1);
            afac = Math.exp(a1 * (lna1 - 1.) - gln);
            pp = (p < 0.5) ? p : 1. - p;
            t = Math.sqrt(-2. * Math.log(pp));
            x = (2.30753 + t * 0.27061) / (1. + t * (0.99229 + t * 0.04481)) - t;
            if (p < 0.5) {
                x = -x;
            }
            x = Math.max(1.e-3, a * Math.pow(1. - 1. / (9. * a) - x / (3. * Math.sqrt(a)), 3));
        } else {
            t = 1.0 - a * (0.253 + a * 0.12);
            if (p < t) {
                x = Math.pow(p / t, 1. / a);
            } else {
                x = 1. - Math.log(1. - (p - t) / (1. - t));
            }
        }
        for (int j = 0; j < 12; j++) {
            if (x <= 0.0) {
                return 0.0;
            }
            err = regularizedIncompleteGamma(a, x) - p;
            if (a > 1.) {
                t = afac * Math.exp(-(x - a1) + a1 * (Math.log(x) - lna1));
            } else {
                t = Math.exp(-x + a1 * Math.log(x) - gln);
            }
            u = err / t;
            x -= (t = u / (1. - 0.5 * Math.min(1., u * ((a - 1.) / x - 1))));
            if (x <= 0.) {
                x = 0.5 * (x + t);
            }
            if (Math.abs(t) < EPS * x) {
                break;
            }
        }
        return x;
    }
397403
private static long specialNPrimeB(final long numSamplesI, final double p, final double delta) {
    final double q, numSamplesF, oneMinusDelta;
    double tot, curTerm;
    long m;
    assertTrue(numSamplesI >= 1);
    assertTrue((0.0 < p) && (p < 1.0));
    assertTrue((0.0 < delta) && (delta < 1.0));
    q = 1.0 - p;
    oneMinusDelta = 1.0 - delta;
    numSamplesF = numSamplesI;
    curTerm = Math.pow(p, numSamplesF);  // curTerm = posteriorProbability (k, k, p)
    assertTrue(curTerm > 1e-100); // sanity check for non-use of logarithms
    tot = curTerm;
    m = numSamplesI;
    while (tot < oneMinusDelta) {
      curTerm = (curTerm * q * (m)) / ((m + 1) - numSamplesI);
      tot += curTerm;
      m += 1;
    }
    return (m); // don't need to back up
  }
622140
function (i, p) {
  var k = Math.floor(Math.log(i) / Math.log(2))
    , r = i - Math.pow(2, k);
  if (r === 0) {
    return Math.pow(2, p - k);
  }
  var nr = (i - 2*r).toString(2).split('').reverse().join('');
  return (parseInt(nr, 2) << p - nr.length) + Math.pow(2, p - k - 1);
}
1149861
def _JRStaeckelIntegrandSquared(u,E,Lz,I3U,delta,u0,sinh2u0,v0,sin2v0,
                                potu0v0,pot):
    #potu0v0= potentialStaeckel(u0,v0,pot,delta)
    """"""
    sinh2u= nu.sinh(u)**2.
    dU= (sinh2u+sin2v0)*potentialStaeckel(u,v0,
                                          pot,delta)\
        -(sinh2u0+sin2v0)*potu0v0
    return E*sinh2u-I3U-dU-Lz**2./2./delta**2./sinh2u
194336
public static double getPvalue(int n11, int ndot1, int ndot2, int n1dot, int n2dot) {
        int originalN11=n11;

        int n=ndot1+ndot2;

        int min=Math.max(0,n1dot+ndot1-n);
        int max=Math.min(n1dot,ndot1);

        Set<Integer> combinations = new HashSet<>();
        combinations.add(originalN11);
        for(int r=min;r<=max;++r) {
            combinations.add(r);
        }

        double standardPart = Arithmetics.factorial(n1dot)*Arithmetics.factorial(n2dot)*Arithmetics.factorial(ndot1)*Arithmetics.factorial(ndot2)/Arithmetics.factorial(n);

        double originalPvalue=0.0;
        double Psum=0;
        for(Integer r : combinations) {
            n11=r;

            int n21=ndot1-n11;
            int n12=n1dot-n11;
            int n22=n2dot-n21;

            double pvalue=standardPart/(Arithmetics.factorial(n11)*Arithmetics.factorial(n12)*Arithmetics.factorial(n21)*Arithmetics.factorial(n22));

            if(r==originalN11) {
                originalPvalue=pvalue;
            }

            if(pvalue<=originalPvalue) {
                Psum+=pvalue;
            }
        }
        //combinations = null;    

        return Psum;
    }
1246992
function(p) {
            var m = this.m;
            return {
                x: m[0] * p.x + m[2] * p.y + m[4],
                y: m[1] * p.x + m[3] * p.y + m[5]
            };
        }
540473
public static void extractColumn(DMatrixRMaj P, int col, GeoTuple3D_F64 a) {
		a.x = P.unsafe_get(0,col);
		a.y = P.unsafe_get(1,col);
		a.z = P.unsafe_get(2,col);
	}
854952
public static void ge_double_scalarmult_vartime(ge_p2 r, byte[] a, ge_p3 A, byte[] b) {
        byte[] aslide = new byte[256];
        byte[] bslide = new byte[256];
        ge_cached Ai[] = new ge_cached[8]; /* A,3A,5A,7A,9A,11A,13A,15A */
        for (int count = 0; count < 8; count++)
            Ai[count] = new ge_cached();
        ge_p1p1 t = new ge_p1p1();
        ge_p3 u = new ge_p3();
        ge_p3 A2 = new ge_p3();
        int i;

        slide(aslide, a);
        slide(bslide, b);

        ge_p3_to_cached.ge_p3_to_cached(Ai[0], A);
        ge_p3_dbl.ge_p3_dbl(t, A);
        ge_p1p1_to_p3.ge_p1p1_to_p3(A2, t);
        ge_add.ge_add(t, A2, Ai[0]);
        ge_p1p1_to_p3.ge_p1p1_to_p3(u, t);
        ge_p3_to_cached.ge_p3_to_cached(Ai[1], u);
        ge_add.ge_add(t, A2, Ai[1]);
        ge_p1p1_to_p3.ge_p1p1_to_p3(u, t);
        ge_p3_to_cached.ge_p3_to_cached(Ai[2], u);
        ge_add.ge_add(t, A2, Ai[2]);
        ge_p1p1_to_p3.ge_p1p1_to_p3(u, t);
        ge_p3_to_cached.ge_p3_to_cached(Ai[3], u);
        ge_add.ge_add(t, A2, Ai[3]);
        ge_p1p1_to_p3.ge_p1p1_to_p3(u, t);
        ge_p3_to_cached.ge_p3_to_cached(Ai[4], u);
        ge_add.ge_add(t, A2, Ai[4]);
        ge_p1p1_to_p3.ge_p1p1_to_p3(u, t);
        ge_p3_to_cached.ge_p3_to_cached(Ai[5], u);
        ge_add.ge_add(t, A2, Ai[5]);
        ge_p1p1_to_p3.ge_p1p1_to_p3(u, t);
        ge_p3_to_cached.ge_p3_to_cached(Ai[6], u);
        ge_add.ge_add(t, A2, Ai[6]);
        ge_p1p1_to_p3.ge_p1p1_to_p3(u, t);
        ge_p3_to_cached.ge_p3_to_cached(Ai[7], u);

        ge_p2_0.ge_p2_0(r);

        for (i = 255; i >= 0; --i) {
            if (aslide[i] != 0 || bslide[i] != 0) break;
        }

        for (; i >= 0; --i) {
            ge_p2_dbl.ge_p2_dbl(t, r);

            if (aslide[i] > 0) {
                ge_p1p1_to_p3.ge_p1p1_to_p3(u, t);
                ge_add.ge_add(t, u, Ai[aslide[i] / 2]);
            } else if (aslide[i] < 0) {
                ge_p1p1_to_p3.ge_p1p1_to_p3(u, t);
                ge_sub.ge_sub(t, u, Ai[(-aslide[i]) / 2]);
            }

            if (bslide[i] > 0) {
                ge_p1p1_to_p3.ge_p1p1_to_p3(u, t);
                ge_madd.ge_madd(t, u, Bi[bslide[i] / 2]);
            } else if (bslide[i] < 0) {
                ge_p1p1_to_p3.ge_p1p1_to_p3(u, t);
                ge_msub.ge_msub(t, u, Bi[(-bslide[i]) / 2]);
            }

            ge_p1p1_to_p2.ge_p1p1_to_p2(r, t);
        }
    }
1038308
function(p,phi,theta) {
        var r = 0;
        r += Math.pow(mathUtils.sin(this.m[0] * theta), this.m[1]);
        r += Math.pow(mathUtils.cos(this.m[2] * theta), this.m[3]);
        r += Math.pow(mathUtils.sin(this.m[4] * phi), this.m[5]);
        r += Math.pow(mathUtils.cos(this.m[6] * phi), this.m[7]);

        var sinTheta = mathUtils.sin(theta);
        p.x = r * sinTheta * mathUtils.cos(phi);
        p.y = r * mathUtils.cos(theta);
        p.z = r * sinTheta * mathUtils.sin(phi);
        return p;
    }
520442
function findx0BigA(a, p, q) {
      var w, s, D, B, u, z, zbar, logpg;

      // snTerm -- see formula 34
      logpg = Math.log(p * gamma(a + 1));
      function snTerm(x) {
         return function(n, v) {
            return n === 0 ? 1 : v * x / (a + n);
         };
      }
      // See formula 34, where `n` is off by 1
      function f(x, n) {
         return Math.exp((logpg + x - Math.log(series(snTerm(x), n + 1))) / a);
      }
      B = q * gamma(a);
      D = Math.max(2, a * (a - 1));
      s = funS(p, q);
      w = a + s * Math.sqrt(a) + (s * s - 1) / 3 +
         (s * s * s - 7 * s) / (36 * Math.sqrt(a)) -
         (3 * s * s * s * s + 7 * s * s - 16) / (810 * a) +
         (9 * s * s * s * s * s + 256 * s * s * s - 433 * s) /
            (38880 * a * Math.sqrt(a));
      if (a >= 500 && Math.abs(1 - w / a) < 1e-6) { return w; }
      if (p > 0.5) {
         if (w < 3 * a) { return w; }
         if (B < Math.pow(10, -D)) { return findx0TinyB(a, -Math.log(B)); }
         u = -Math.log(B) + (a - 1) * Math.log(w) -
               Math.log(1 + (1 - a) / (1 + w));
         // formula 33
         return -Math.log(B) + (a - 1) * Math.log(u) -
               Math.log(1 + (1 - a) / (1 + u));
      }
      // handle p <= 0.5
      z = w > 0.15 * (a + 1) ? w : f(f(f(f(w, 1), 3), 3), 4);
      if (z <= 0.01 * (a + 1) || z > 0.7 * (a + 1)) { return z; }
      /* eslint-disable no-extra-parens */
      return (function() {
         // find N for formula 36
         var N, term;

         N = 0;
         term = 1;
         do {
            N += 1;
            term = snTerm(z)(N, term);
         } while (term > 1e-4);
         zbar = f(z, N + 1);
         return zbar * (1 - (a * Math.log(zbar) - zbar - logpg +
                  Math.log(series(snTerm(z), N + 1))) / (a - zbar));
      }());
      /* eslint-enable no-extra-parens */
   }
256028
public static float normalizeAngle (float a) {
        while (a < -FloatMath.PI) {
            a += TWO_PI;
        }
        while (a >= FloatMath.PI) {
            a -= TWO_PI;
        }
        return a;
    }
599583
public static final Point mul(double k, Point p)
    {
        return new AbstractPoint(k*p.getX(), k*p.getY());
    }
408273
function calcBezierAtT(p, t) {
    var x = (1-t)*(1-t)*(1-t)*p[0].x + 3*(1-t)*(1-t)*t*p[1].x + 3*(1-t)*t*t*p[2].x + t*t*t*p[3].x;
    var y = (1-t)*(1-t)*(1-t)*p[0].y + 3*(1-t)*(1-t)*t*p[1].y + 3*(1-t)*t*t*p[2].y + t*t*t*p[3].y;
    return new Point(x, y);
}