/* Square a number.  Optimization of  mul25519(x2, x, x)
517002
private static final long10 sqr(long10 x2, long10 x) {
		long
			x_0=x._0,x_1=x._1,x_2=x._2,x_3=x._3,x_4=x._4,
			x_5=x._5,x_6=x._6,x_7=x._7,x_8=x._8,x_9=x._9;
		long t;
		t = (x_4*x_4) + 2 * ((x_0*x_8) + (x_2*x_6)) + 38 *
			(x_9*x_9) + 4 * ((x_1*x_7) + (x_3*x_5));
		x2._8 = (t & ((1 << 26) - 1));
		t = (t >> 26) + 2 * ((x_0*x_9) + (x_1*x_8) + (x_2*x_7) +
				(x_3*x_6) + (x_4*x_5));
		x2._9 = (t & ((1 << 25) - 1));
		t = 19 * (t >> 25) + (x_0*x_0) + 38 * ((x_2*x_8) +
				(x_4*x_6) + (x_5*x_5)) + 76 * ((x_1*x_9)
				+ (x_3*x_7));
		x2._0 = (t & ((1 << 26) - 1));
		t = (t >> 26) + 2 * (x_0*x_1) + 38 * ((x_2*x_9) +
				(x_3*x_8) + (x_4*x_7) + (x_5*x_6));
		x2._1 = (t & ((1 << 25) - 1));
		t = (t >> 25) + 19 * (x_6*x_6) + 2 * ((x_0*x_2) +
				(x_1*x_1)) + 38 * (x_4*x_8) + 76 *
				((x_3*x_9) + (x_5*x_7));
		x2._2 = (t & ((1 << 26) - 1));
		t = (t >> 26) + 2 * ((x_0*x_3) + (x_1*x_2)) + 38 *
			((x_4*x_9) + (x_5*x_8) + (x_6*x_7));
		x2._3 = (t & ((1 << 25) - 1));
		t = (t >> 25) + (x_2*x_2) + 2 * (x_0*x_4) + 38 *
			((x_6*x_8) + (x_7*x_7)) + 4 * (x_1*x_3) + 76 *
			(x_5*x_9);
		x2._4 = (t & ((1 << 26) - 1));
		t = (t >> 26) + 2 * ((x_0*x_5) + (x_1*x_4) + (x_2*x_3))
			+ 38 * ((x_6*x_9) + (x_7*x_8));
		x2._5 = (t & ((1 << 25) - 1));
		t = (t >> 25) + 19 * (x_8*x_8) + 2 * ((x_0*x_6) +
				(x_2*x_4) + (x_3*x_3)) + 4 * (x_1*x_5) +
				76 * (x_7*x_9);
		x2._6 = (t & ((1 << 26) - 1));
		t = (t >> 26) + 2 * ((x_0*x_7) + (x_1*x_6) + (x_2*x_5) +
				(x_3*x_4)) + 38 * (x_8*x_9);
		x2._7 = (t & ((1 << 25) - 1));
		t = (t >> 25) + x2._8;
		x2._8 = (t & ((1 << 26) - 1));
		x2._9 += (t >> 26);
		return x2;
	}
1172695
public static void square( DMatrix1Row mat )
    {
        int index = 1;
        int indexEnd = mat.numCols;
        for( int i = 0; i < mat.numRows;
             i++ , index += i+1 , indexEnd += mat.numCols ) {
            int indexOther = (i+1)*mat.numCols + i;
            for( ; index < indexEnd; index++, indexOther += mat.numCols) {
                double val = mat.data[ index ];
                mat.data[ index ] = mat.data[ indexOther ];
                mat.data[indexOther] = val;
            }
        }
    }
1012259
function choose(n, k) {
      var ret, j;

      k = Math.round(k);
      if (utils.hasNaN(n, k)) { return NaN; }
      if (k < kSmallMax) {
         if (n - k < k && n >= 0 && n === Math.round(n)) { k = n - k; }
         if (k < 0) { return 0; }
         if (k === 0) { return 1; }
         ret = n;
         for (j = 2; j <= k; j += 1) {
            ret *= (n - j + 1) / j;
         }

         return n === Math.round(n) ? Math.round(ret) : ret;
      }
      if (n < 0) {
         ret = choose(-n + k - 1, k);

         return k % 2 === 1 ? -ret : ret;
      }
      if (n === Math.round(n)) {
         if (n < k) { return 0; }
         if (n - k < kSmallMax) { return choose(n, n - k); }

         return Math.round(Math.exp(lfastchoose(n, k)));
      }
      if (n < k - 1) {
         return signGamma(n - k + 1) * Math.exp(lfastchoose2(n, k));
      }

      return Math.exp(lfastchoose(n, k));
   }
208523
private double optimize(int k0, int kmax, double sumx, double sumx2, double xp, double yp, double sumxy, double sumy) {
    int k = kmax - k0 + 1;
    return (sumxy - xp * sumy - yp * sumx + k * xp * yp) / (sumx2 - 2 * sumx * xp + k * xp * xp);
    // return (-xp * yp * + yp * sumx - sumxy + xp * sumy) / (-xp * xp *
    // kmax - sumx2 + 2 * xp * sumx);
  }
1294034
public void mul(IMatrix b, IMatrix result) {
        if ((b == null) || (columns != b.rows)) return;

        if ((result.rows != rows) || (result.columns != b.columns)) result.reshape(rows, b.columns);
        int i, j, k;
        double realsum, imagsum;
        for (i = 0; i < rows; i++)
            for (k = 0; k < b.columns; k++) {
                realsum = 0;
                imagsum = 0;
                for (j = 0; j < columns; j++) {
                    realsum += realmatrix[i][j] * b.realmatrix[j][k] - imagmatrix[i][j] * b.imagmatrix[j][k];
                    imagsum += realmatrix[i][j] * b.imagmatrix[j][k] + imagmatrix[i][j] * b.realmatrix[j][k];
                }
                result.realmatrix[i][k] = realsum;
                result.imagmatrix[i][k] = imagsum;
            }
    }
510727
def _round_to_next_multiple(n, m):
    
    return n if n % m == 0 else n + m - n % m
601217
void pass_2(int fi, 
	      float in[],  int in0, int istride,
	      float out[], int out0, int ostride,
	      int sign, int product) {
    int k, k1;

    int factor = 2;
    int m = n / factor;
    int q = n / product;
    int product_1 = product / factor;

    int di = istride * m;
    int dj = ostride * product_1;
    int i = in0, j = out0;
    float x_real, x_imag;
    for (k = 0; k < q; k++) {
      float twids[] = twiddle[fi][k];
      float w_real =       twids[0];
      float w_imag = -sign*twids[1];

      for (k1 = 0; k1 < product_1; k1++) {
	float z0_real = in[i];
	float z0_imag = in[i+1];
	float z1_real = in[i+di];
	float z1_imag = in[i+di+1];
	i += istride;

	/* compute x = W(2) z */

	/* apply twiddle factors */
	  
	/* out0 = 1 * (z0 + z1) */
	out[j]   = z0_real + z1_real;
	out[j+1] = z0_imag + z1_imag;
	  
	/* out1 = w * (z0 - z1) */
	x_real = z0_real - z1_real;
	x_imag = z0_imag - z1_imag;
	out[j+dj]   = w_real * x_real - w_imag * x_imag;
	out[j+dj+1] = w_real * x_imag + w_imag * x_real;
	  
	j += ostride;
      }
      j += (factor-1)*dj;
    }}
241112
public void mulRow(int i, double[] arr, int offset)
    {
        int n = cols;
        for (int j=0;j<n;j++)
        {
            mul(i, j, arr[j+offset]);
        }
    }
866810
private int mulsubLong(int[] q, int dh, int dl, int x, int offset) {
        long xLong = x & LONG_MASK;
        offset += 2;
        long product = (dl & LONG_MASK) * xLong;
        long difference = q[offset] - product;
        q[offset--] = (int)difference;
        long carry = (product >>> 32)
                 + (((difference & LONG_MASK) >
                     (((~(int)product) & LONG_MASK))) ? 1:0);
        product = (dh & LONG_MASK) * xLong + carry;
        difference = q[offset] - product;
        q[offset--] = (int)difference;
        carry = (product >>> 32)
                 + (((difference & LONG_MASK) >
                     (((~(int)product) & LONG_MASK))) ? 1:0);
        return (int)carry;
    }
901111
function mpn_cmp(n1, n2, l) {
	while (l > 0) {
		var c = n1[--l] - n2[l];
		if (c != 0) {
			return c;
		}
	}
	return 0;
}
203970
function addNum (num1, num2) {
  let sq1, sq2
  try {
    sq1 = _toString(num1).split('.')[1].length
  } catch (e) {
    sq1 = 0
  }
  try {
    sq2 = _toString(num2).split('.')[1].length
  } catch (e) {
    sq2 = 0
  }
  const m = Math.pow(10, Math.max(sq1, sq2))
  return (Math.round(num1 * m) + Math.round(num2 * m)) / m
}
274340
function fp_sqr(r, n, p, l) {
	return fp_mul(r, n, n, p, l);
}
984774
function blockMult(a, b, rows, cols) {
                            // For small matrices, resort to naive multiplication.
                            if (rows <= 512 || cols <= 512) {
                                return a.mmul(b); // a is equivalent to this
                            }

                            // Apply dynamic padding.
                            if ((rows % 2 === 1) && (cols % 2 === 1)) {
                                a = embed(a, rows + 1, cols + 1);
                                b = embed(b, rows + 1, cols + 1);
                            } else if (rows % 2 === 1) {
                                a = embed(a, rows + 1, cols);
                                b = embed(b, rows + 1, cols);
                            } else if (cols % 2 === 1) {
                                a = embed(a, rows, cols + 1);
                                b = embed(b, rows, cols + 1);
                            }

                            var halfRows = parseInt(a.rows / 2);
                            var halfCols = parseInt(a.columns / 2);
                            // Subdivide input matrices.
                            var a11 = a.subMatrix(0, halfRows - 1, 0, halfCols - 1);
                            var b11 = b.subMatrix(0, halfRows - 1, 0, halfCols - 1);

                            var a12 = a.subMatrix(0, halfRows - 1, halfCols, a.columns - 1);
                            var b12 = b.subMatrix(0, halfRows - 1, halfCols, b.columns - 1);

                            var a21 = a.subMatrix(halfRows, a.rows - 1, 0, halfCols - 1);
                            var b21 = b.subMatrix(halfRows, b.rows - 1, 0, halfCols - 1);

                            var a22 = a.subMatrix(halfRows, a.rows - 1, halfCols, a.columns - 1);
                            var b22 = b.subMatrix(halfRows, b.rows - 1, halfCols, b.columns - 1);

                            // Compute intermediate values.
                            var m1 = blockMult(Matrix.add(a11, a22), Matrix.add(b11, b22), halfRows, halfCols);
                            var m2 = blockMult(Matrix.add(a21, a22), b11, halfRows, halfCols);
                            var m3 = blockMult(a11, Matrix.sub(b12, b22), halfRows, halfCols);
                            var m4 = blockMult(a22, Matrix.sub(b21, b11), halfRows, halfCols);
                            var m5 = blockMult(Matrix.add(a11, a12), b22, halfRows, halfCols);
                            var m6 = blockMult(Matrix.sub(a21, a11), Matrix.add(b11, b12), halfRows, halfCols);
                            var m7 = blockMult(Matrix.sub(a12, a22), Matrix.add(b21, b22), halfRows, halfCols);

                            // Combine intermediate values into the output.
                            var c11 = Matrix.add(m1, m4);
                            c11.sub(m5);
                            c11.add(m7);
                            var c12 = Matrix.add(m3, m5);
                            var c21 = Matrix.add(m2, m4);
                            var c22 = Matrix.sub(m1, m2);
                            c22.add(m3);
                            c22.add(m6);

                            //Crop output to the desired size (undo dynamic padding).
                            var resultat = Matrix.zeros(2 * c11.rows, 2 * c11.columns);
                            resultat = resultat.setSubMatrix(c11, 0, 0);
                            resultat = resultat.setSubMatrix(c12, c11.rows, 0);
                            resultat = resultat.setSubMatrix(c21, 0, c11.columns);
                            resultat = resultat.setSubMatrix(c22, c11.rows, c11.columns);
                            return resultat.subMatrix(0, rows - 1, 0, cols - 1);
                        }
213782
public static double SquaredEuclidean(double[] x, double[] y) {
        double d = 0.0, u;

        for (int i = 0; i < x.length; i++) {
            u = x[i] - y[i];
            d += u * u;
        }

        return d;
    }
25987
private static void mult(int[] src, int srcLen, int v0, int v1, int[] dst) {
        long v = v0 & LONG_MASK;
        long carry = 0;
        for (int j = 0; j < srcLen; j++) {
            long product = v * (src[j] & LONG_MASK) + carry;
            dst[j] = (int) product;
            carry = product >>> 32;
        }
        dst[srcLen] = (int) carry;
        v = v1 & LONG_MASK;
        carry = 0;
        for (int j = 0; j < srcLen; j++) {
            long product = (dst[j + 1] & LONG_MASK) + v * (src[j] & LONG_MASK) + carry;
            dst[j + 1] = (int) product;
            carry = product >>> 32;
        }
        dst[srcLen + 1] = (int) carry;
    }
209379
function mpn_add(r, n1, n2, l) {
	var c = 0;
	for (var i = 0; i < l; ++i) {
		r[i] = (c = (c >>> 31) + n1[i] + n2[i]) & 0x7FFFFFFF;
	}
	return c >>> 31;
}
330734
private static final int mula_small(byte[] p,byte[] q,int m,byte[] x,int n,int z) {
		int v=0;
		for (int i=0;i<n;++i) {
			v+=(q[i+m] & 0xFF)+z*(x[i] & 0xFF);
			p[i+m]=(byte)v;
			v>>=8;
		}
		return v;		
	}