/* y^2 = x^3 + 486662 x^2 + x  over GF(2^255-19) /* t1 = ax + az
t2 = ax - az
766368
function mont_prep (t1, t2, ax, az) {
        add(t1, ax, az);
        sub(t2, ax, az);
    }
188910
func (x Int) Mod(y Int) Int {
	var quo, rem big.Int
	quo.QuoRem(x.bigint, y.bigint, &rem)
	if (x.bigint.Sign() < 0) != (y.bigint.Sign() < 0) && rem.Sign() != 0 {
		rem.Add(&rem, y.bigint)
	}
	return Int{&rem}
}
632540
function x_to_y2 (t, y2, x) {
        sqr(t, x);
        mul_small(y2, x, 486662);
        add(t, t, y2);
        add(t, t, C1);
        mul(y2, t, x);
    }
129479
private void computeDandH( DenseMatrix64F param , DenseMatrix64F x , DenseMatrix64F y )
    {
        func.compute(param,x, tempDH);
        subtractEquals(tempDH, y);

        if (jacobianFactory != null)
        {
            jacobianFactory.computeJacobian(param, x, jacobian);
        }
        else
        {
            computeNumericalJacobian(param,x,jacobian);
        }

        int numParam = param.getNumElements();
        int length = y.getNumElements();

        // d = average{ (f(x_i;p) - y_i) * jacobian(:,i) }
        for( int i = 0; i < numParam; i++ ) {
            double total = 0;
            for( int j = 0; j < length; j++ ) {
                total += tempDH.get(j,0)*jacobian.get(i,j);
            }
            d.set(i,0,total/length);
        }

        // compute the approximation of the hessian
        multTransB(jacobian,jacobian,H);
        scale(1.0/length,H);
    }
199158
function int64add(dst, x, y) {
        var w0 = (x.l & 0xffff) + (y.l & 0xffff);
        var w1 = (x.l >>> 16) + (y.l >>> 16) + (w0 >>> 16);
        var w2 = (x.h & 0xffff) + (y.h & 0xffff) + (w1 >>> 16);
        var w3 = (x.h >>> 16) + (y.h >>> 16) + (w2 >>> 16);
        dst.l = (w0 & 0xffff) | (w1 << 16);
        dst.h = (w2 & 0xffff) | (w3 << 16);
      }
1240683
function recip (y, x, sqrtassist) {
        var t0 = createUnpackedArray();
        var t1 = createUnpackedArray();
        var t2 = createUnpackedArray();
        var t3 = createUnpackedArray();
        var t4 = createUnpackedArray();

        /* the chain for x^(2^255-21) is straight from djb's implementation */
        var i;
        sqr(t1, x); /*  2 === 2 * 1	*/
        sqr(t2, t1); /*  4 === 2 * 2	*/
        sqr(t0, t2); /*  8 === 2 * 4	*/
        mul(t2, t0, x); /*  9 === 8 + 1	*/
        mul(t0, t2, t1); /* 11 === 9 + 2	*/
        sqr(t1, t0); /* 22 === 2 * 11	*/
        mul(t3, t1, t2); /* 31 === 22 + 9 === 2^5   - 2^0	*/
        sqr(t1, t3); /* 2^6   - 2^1	*/
        sqr(t2, t1); /* 2^7   - 2^2	*/
        sqr(t1, t2); /* 2^8   - 2^3	*/
        sqr(t2, t1); /* 2^9   - 2^4	*/
        sqr(t1, t2); /* 2^10  - 2^5	*/
        mul(t2, t1, t3); /* 2^10  - 2^0	*/
        sqr(t1, t2); /* 2^11  - 2^1	*/
        sqr(t3, t1); /* 2^12  - 2^2	*/
        for (i = 1; i < 5; i++) {
            sqr(t1, t3);
            sqr(t3, t1);
        } /* t3 */ /* 2^20  - 2^10	*/
        mul(t1, t3, t2); /* 2^20  - 2^0	*/
        sqr(t3, t1); /* 2^21  - 2^1	*/
        sqr(t4, t3); /* 2^22  - 2^2	*/
        for (i = 1; i < 10; i++) {
            sqr(t3, t4);
            sqr(t4, t3);
        } /* t4 */ /* 2^40  - 2^20	*/
        mul(t3, t4, t1); /* 2^40  - 2^0	*/
        for (i = 0; i < 5; i++) {
            sqr(t1, t3);
            sqr(t3, t1);
        } /* t3 */ /* 2^50  - 2^10	*/
        mul(t1, t3, t2); /* 2^50  - 2^0	*/
        sqr(t2, t1); /* 2^51  - 2^1	*/
        sqr(t3, t2); /* 2^52  - 2^2	*/
        for (i = 1; i < 25; i++) {
            sqr(t2, t3);
            sqr(t3, t2);
        } /* t3 */ /* 2^100 - 2^50 */
        mul(t2, t3, t1); /* 2^100 - 2^0	*/
        sqr(t3, t2); /* 2^101 - 2^1	*/
        sqr(t4, t3); /* 2^102 - 2^2	*/
        for (i = 1; i < 50; i++) {
            sqr(t3, t4);
            sqr(t4, t3);
        } /* t4 */ /* 2^200 - 2^100 */
        mul(t3, t4, t2); /* 2^200 - 2^0	*/
        for (i = 0; i < 25; i++) {
            sqr(t4, t3);
            sqr(t3, t4);
        } /* t3 */ /* 2^250 - 2^50	*/
        mul(t2, t3, t1); /* 2^250 - 2^0	*/
        sqr(t1, t2); /* 2^251 - 2^1	*/
        sqr(t2, t1); /* 2^252 - 2^2	*/
        if (sqrtassist !== 0) {
            mul(y, x, t2); /* 2^252 - 3 */
        } else {
            sqr(t1, t2); /* 2^253 - 2^3	*/
            sqr(t2, t1); /* 2^254 - 2^4	*/
            sqr(t1, t2); /* 2^255 - 2^5	*/
            mul(y, t1, t0); /* 2^255 - 21	*/
        }
    }
10617
function f1(d, m, r) {
    const t = m + d;
    const I = (t << r) | (t >>> (32 - r));
    return ((sBox[0][I >>> 24] ^ sBox[1][(I >>> 16) & 255]) - sBox[2][(I >>> 8) & 255]) + sBox[3][I & 255];
  }
306748
func f(r0, r1, r2, r3 *uint32, k0, k1 uint32) {
	k0 ^= *r0
	k1 ^= *r1

	t := sbox4_4404[byte(k0)]
	t ^= sbox3_3033[byte(k0>>8)]
	t ^= sbox2_0222[byte(k0>>16)]
	t ^= sbox1_1110[byte(k0>>24)]
	*r3 ^= (t >> 8) | (t << (32 - 8))

	k0 = t
	k0 ^= sbox1_1110[byte(k1)]
	k0 ^= sbox4_4404[byte(k1>>8)]
	k0 ^= sbox3_3033[byte(k1>>16)]
	k0 ^= sbox2_0222[byte(k1>>24)]

	*r2 ^= k0
	*r3 ^= k0
}
221584
def _g_a (self, x, a, b, s):
        
        return 2*b/(1.0+np.exp(s*(x-a)))
130915
static public double ScafferF6(double x, double y) {
    double temp1 = x * x + y * y;
    double temp2 = Math.sin(Math.sqrt(temp1));
    double temp3 = 1.0 + 0.001 * temp1;
    return (0.5 + ((temp2 * temp2 - 0.5) / (temp3 * temp3)));
  }
409022
function egcd32 (x, y, a, b) {
        var an, bn = 32, qn, i;
        for (i = 0; i < 32; i++)
            x[i] = y[i] = 0;
        x[0] = 1;
        an = numsize(a, 32);
        if (an === 0)
            return y; /* division by zero */
        var temp = new Array(32);
        while (true) {
            qn = bn - an + 1;
            divmod(temp, b, bn, a, an);
            bn = numsize(b, bn);
            if (bn === 0)
                return x;
            mula32(y, x, temp, qn, -1);

            qn = an - bn + 1;
            divmod(temp, a, an, b, bn);
            an = numsize(a, an);
            if (an === 0)
                return y;
            mula32(x, y, temp, qn, -1);
        }
    }
889326
function pchi2(n, y) {
        var w, eps, v, s, qe;

        if (n === 1) {
            w = pnorm(1 - y / 2);
            return w * w;
        } else if (n === 2) {
            return -2.0 * Math.log(y);
        }

        eps = 1.0e-5;
        v = 0.0;
        s = 10.0;
        while (true) {
            v = v + 2;

            if (s <= eps) { break; }
            if ((qe = q_chi2(n, v) - y) === 0.0) { break; }
            if (qe < 0.0) {
                v = v - s;
                s = s / 10.0;
            }
        }
        return v;
    }
1008607
function(d2) {

      var w = this['w'];
      var x = this['x'];
      var y = this['y'];
      var z = this['z'];

      var n = w * w + x * x + y * y + z * z;
      var s = n === 0 ? 0 : 2 / n;
      var wx = s * w * x, wy = s * w * y, wz = s * w * z;
      var xx = s * x * x, xy = s * x * y, xz = s * x * z;
      var yy = s * y * y, yz = s * y * z, zz = s * z * z;

      if (d2) {
        return [
          [1 - (yy + zz), xy - wz, xz + wy],
          [xy + wz, 1 - (xx + zz), yz - wx],
          [xz - wy, yz + wx, 1 - (xx + yy)]];
      }

      return [
        1 - (yy + zz), xy - wz, xz + wy,
        xy + wz, 1 - (xx + zz), yz - wx,
        xz - wy, yz + wx, 1 - (xx + yy)];
    }
23777
function() { 
            var a = this.clone(),
                i = 1,
                b = a.clone().diff(),
                c = a.clone().gcd(b),
                w = a.divide(c)[0];
            var output = Polynomial.fromArray([new Frac(1)], a.variable);
            while(!c.equalsNumber(1)) { 
                var y = w.gcd(c); 
                var z = w.divide(y)[0];
                //one of the factors may have shown up since it's square but smaller than the 
                //one where finding
                if(!z.equalsNumber(1) && i>1) {
                    var t = z.clone();
                    for(var j=1; j<i; j++)
                        t.multiply(z.clone());
                    z = t;
                }
                output = output.multiply(z); 
                i++;
                w = y;
                c = c.divide(y)[0];
            }
            return [output, w, i];
        }
231211
function rsRem(x) {

        var b = x >>> 24 & 0xFF
          , g2 = (b << 1 ^ ((b & 0x80) !== 0 ? RS_GF_FDBK : 0)) & 0xFF
          , g3 = b >>> 1 ^ ((b & 0x01) !== 0 ? RS_GF_FDBK >>> 1 : 0 ) ^ g2
          , result = x << 8 ^ g3 << 24 ^ g2 << 16 ^ g3 << 8 ^ b;
        return result;
      }
122015
func feSquare(h, f *fieldElement) {
	f0 := f[0]
	f1 := f[1]
	f2 := f[2]
	f3 := f[3]
	f4 := f[4]
	f5 := f[5]
	f6 := f[6]
	f7 := f[7]
	f8 := f[8]
	f9 := f[9]
	f0_2 := 2 * f0
	f1_2 := 2 * f1
	f2_2 := 2 * f2
	f3_2 := 2 * f3
	f4_2 := 2 * f4
	f5_2 := 2 * f5
	f6_2 := 2 * f6
	f7_2 := 2 * f7
	f5_38 := 38 * f5 // 1.31*2^30
	f6_19 := 19 * f6 // 1.31*2^30
	f7_38 := 38 * f7 // 1.31*2^30
	f8_19 := 19 * f8 // 1.31*2^30
	f9_38 := 38 * f9 // 1.31*2^30
	f0f0 := int64(f0) * int64(f0)
	f0f1_2 := int64(f0_2) * int64(f1)
	f0f2_2 := int64(f0_2) * int64(f2)
	f0f3_2 := int64(f0_2) * int64(f3)
	f0f4_2 := int64(f0_2) * int64(f4)
	f0f5_2 := int64(f0_2) * int64(f5)
	f0f6_2 := int64(f0_2) * int64(f6)
	f0f7_2 := int64(f0_2) * int64(f7)
	f0f8_2 := int64(f0_2) * int64(f8)
	f0f9_2 := int64(f0_2) * int64(f9)
	f1f1_2 := int64(f1_2) * int64(f1)
	f1f2_2 := int64(f1_2) * int64(f2)
	f1f3_4 := int64(f1_2) * int64(f3_2)
	f1f4_2 := int64(f1_2) * int64(f4)
	f1f5_4 := int64(f1_2) * int64(f5_2)
	f1f6_2 := int64(f1_2) * int64(f6)
	f1f7_4 := int64(f1_2) * int64(f7_2)
	f1f8_2 := int64(f1_2) * int64(f8)
	f1f9_76 := int64(f1_2) * int64(f9_38)
	f2f2 := int64(f2) * int64(f2)
	f2f3_2 := int64(f2_2) * int64(f3)
	f2f4_2 := int64(f2_2) * int64(f4)
	f2f5_2 := int64(f2_2) * int64(f5)
	f2f6_2 := int64(f2_2) * int64(f6)
	f2f7_2 := int64(f2_2) * int64(f7)
	f2f8_38 := int64(f2_2) * int64(f8_19)
	f2f9_38 := int64(f2) * int64(f9_38)
	f3f3_2 := int64(f3_2) * int64(f3)
	f3f4_2 := int64(f3_2) * int64(f4)
	f3f5_4 := int64(f3_2) * int64(f5_2)
	f3f6_2 := int64(f3_2) * int64(f6)
	f3f7_76 := int64(f3_2) * int64(f7_38)
	f3f8_38 := int64(f3_2) * int64(f8_19)
	f3f9_76 := int64(f3_2) * int64(f9_38)
	f4f4 := int64(f4) * int64(f4)
	f4f5_2 := int64(f4_2) * int64(f5)
	f4f6_38 := int64(f4_2) * int64(f6_19)
	f4f7_38 := int64(f4) * int64(f7_38)
	f4f8_38 := int64(f4_2) * int64(f8_19)
	f4f9_38 := int64(f4) * int64(f9_38)
	f5f5_38 := int64(f5) * int64(f5_38)
	f5f6_38 := int64(f5_2) * int64(f6_19)
	f5f7_76 := int64(f5_2) * int64(f7_38)
	f5f8_38 := int64(f5_2) * int64(f8_19)
	f5f9_76 := int64(f5_2) * int64(f9_38)
	f6f6_19 := int64(f6) * int64(f6_19)
	f6f7_38 := int64(f6) * int64(f7_38)
	f6f8_38 := int64(f6_2) * int64(f8_19)
	f6f9_38 := int64(f6) * int64(f9_38)
	f7f7_38 := int64(f7) * int64(f7_38)
	f7f8_38 := int64(f7_2) * int64(f8_19)
	f7f9_76 := int64(f7_2) * int64(f9_38)
	f8f8_19 := int64(f8) * int64(f8_19)
	f8f9_38 := int64(f8) * int64(f9_38)
	f9f9_38 := int64(f9) * int64(f9_38)
	h0 := f0f0 + f1f9_76 + f2f8_38 + f3f7_76 + f4f6_38 + f5f5_38
	h1 := f0f1_2 + f2f9_38 + f3f8_38 + f4f7_38 + f5f6_38
	h2 := f0f2_2 + f1f1_2 + f3f9_76 + f4f8_38 + f5f7_76 + f6f6_19
	h3 := f0f3_2 + f1f2_2 + f4f9_38 + f5f8_38 + f6f7_38
	h4 := f0f4_2 + f1f3_4 + f2f2 + f5f9_76 + f6f8_38 + f7f7_38
	h5 := f0f5_2 + f1f4_2 + f2f3_2 + f6f9_38 + f7f8_38
	h6 := f0f6_2 + f1f5_4 + f2f4_2 + f3f3_2 + f7f9_76 + f8f8_19
	h7 := f0f7_2 + f1f6_2 + f2f5_2 + f3f4_2 + f8f9_38
	h8 := f0f8_2 + f1f7_4 + f2f6_2 + f3f5_4 + f4f4 + f9f9_38
	h9 := f0f9_2 + f1f8_2 + f2f7_2 + f3f6_2 + f4f5_2
	var carry [10]int64

	carry[0] = (h0 + (1 << 25)) >> 26
	h1 += carry[0]
	h0 -= carry[0] << 26
	carry[4] = (h4 + (1 << 25)) >> 26
	h5 += carry[4]
	h4 -= carry[4] << 26

	carry[1] = (h1 + (1 << 24)) >> 25
	h2 += carry[1]
	h1 -= carry[1] << 25
	carry[5] = (h5 + (1 << 24)) >> 25
	h6 += carry[5]
	h5 -= carry[5] << 25

	carry[2] = (h2 + (1 << 25)) >> 26
	h3 += carry[2]
	h2 -= carry[2] << 26
	carry[6] = (h6 + (1 << 25)) >> 26
	h7 += carry[6]
	h6 -= carry[6] << 26

	carry[3] = (h3 + (1 << 24)) >> 25
	h4 += carry[3]
	h3 -= carry[3] << 25
	carry[7] = (h7 + (1 << 24)) >> 25
	h8 += carry[7]
	h7 -= carry[7] << 25

	carry[4] = (h4 + (1 << 25)) >> 26
	h5 += carry[4]
	h4 -= carry[4] << 26
	carry[8] = (h8 + (1 << 25)) >> 26
	h9 += carry[8]
	h8 -= carry[8] << 26

	carry[9] = (h9 + (1 << 24)) >> 25
	h0 += carry[9] * 19
	h9 -= carry[9] << 25

	carry[0] = (h0 + (1 << 25)) >> 26
	h1 += carry[0]
	h0 -= carry[0] << 26

	h[0] = int32(h0)
	h[1] = int32(h1)
	h[2] = int32(h2)
	h[3] = int32(h3)
	h[4] = int32(h4)
	h[5] = int32(h5)
	h[6] = int32(h6)
	h[7] = int32(h7)
	h[8] = int32(h8)
	h[9] = int32(h9)
}
1040434
function findx0TinyB(a, y) {
      var c1;

      c1 = (a - 1) * Math.log(y);
      return y + Polynomial.new([
         // c5
         (a - 1) * Polynomial.new([
            -0.25, (11 * a - 17) / 6, -3 * a * a + 13 * a - 13,
            (2 * a * a * a - 25 * a * a + 72 * a - 61) / 2,
            (25 * a * a * a - 195 * a * a + 477 * a - 379) / 12
         ]).evalAt(c1),
         // c4
         (a - 1) * Polynomial.new([
            1 / 3, (3 * a - 5) / 2, a * a - 6 * a + 7,
            (11 * a * a - 46 * a + 47) / 6
         ]).evalAt(c1),
         // c3
         (a - 1) * (-0.5 * c1 * c1 + (a - 2) * c1 + (3 * a - 5) / 2),
         // c2
         (a - 1) * (1 + c1),
         // c1
         c1
      ]).evalAt(1 / y);
   }