$= z^{-1} mod q$, for z != 0
963223
def inv(z: int) -> int:
    """"""
    # Adapted from curve25519_athlon.c in djb's Curve25519.
    z2 = z * z % q  # 2
    z9 = pow2(z2, 2) * z % q  # 9
    z11 = z9 * z2 % q  # 11
    z2_5_0 = (z11 * z11) % q * z9 % q  # 31 == 2^5 - 2^0
    z2_10_0 = pow2(z2_5_0, 5) * z2_5_0 % q  # 2^10 - 2^0
    z2_20_0 = pow2(z2_10_0, 10) * z2_10_0 % q  # ...
    z2_40_0 = pow2(z2_20_0, 20) * z2_20_0 % q
    z2_50_0 = pow2(z2_40_0, 10) * z2_10_0 % q
    z2_100_0 = pow2(z2_50_0, 50) * z2_50_0 % q
    z2_200_0 = pow2(z2_100_0, 100) * z2_100_0 % q
    z2_250_0 = pow2(z2_200_0, 50) * z2_50_0 % q  # 2^250 - 2^0
    return pow2(z2_250_0, 5) * z11 % q
555844
def  dBinaryRochedz(r, D, q, F):
    
    return -r[2]*(r[0]*r[0]+r[1]*r[1]+r[2]*r[2])**-1.5 -q*r[2]*((r[0]-D)*(r[0]-D)+r[1]*r[1]+r[2]*r[2])**-1.5
694894
def modInv(a, m):
    """"""
    # notice that a and m need to co-prime to each other.
    if coPrime([a, m]):
        linearCombination = extendedEuclid(a, m)
        return linearCombination[1] % m
    else:
        return 0
1263281
def reduce_angle2(a,k0=nil) # divisor of pi or nil for pi*2
        # we could reduce first to pi*2 to avoid the mod k0 operation
        k,r,divisor = num_class.context do
          num_class.context.precision *= 3
          m = k0.nil? ? one_cycle : half_cycle/k0
          a.divmod(m)+[m]
        end
        [r, k.modulo(k0*2).to_i, divisor]
      end
905695
function mod(n, a, b) {
	n = intOf(n);
	a = intOf(a);
	b = intOf(b);
	var rem;
	if (a < 0 || b < 0) {
		var places = (b - a);
		rem = (n - a) % (places + 1);
		rem = rem < 0 ? (rem + (places + 1)) : rem === 0 ? 0 : rem;
		return rem - (places - b);
	}
	if (n === b) { return n; }
	if (n === b + 1) { return a; }
	if (n === a - 1) { return b; }
	rem = n % (b || 1);
	rem = rem < a ? (rem + b) : rem === 0 ? 0 : rem;
	return rem;
}
219428
private static function powi($z, $n) {
        if ($n < 0) return static::div(1,static::powi($z, -$n));

        if ($n == 0) return new Complex(1,0);

        $y = new Complex(1,0);
        while ($n > 1) {
            if ($n % 2 == 0) {
                $n = $n / 2;
            } else {
                $y = static::mul($z, $y);
                $n = ($n-1)/2;
            }
            $z = static::mul($z, $z);
        }

        return static::mul($z, $y);
    }
186183
def factorial_recur(n, mod=None):
    """"""
    if not (isinstance(n, int) and n >= 0):
        raise ValueError("'n' must be a non-negative integer.")
    if mod is not None and not (isinstance(mod, int) and mod > 0):
        raise ValueError("'mod' must be a positive integer")
    if n == 0:
        return 1
    result = n * factorial(n - 1, mod)
    if mod:
        result %= mod
    return result
1340169
function divmod (q, r, n, d, t) {
        n = n | 0;
        t = t | 0;

        var rn = 0;
        var dt = (d[t - 1] & 0xFF) << 8;
        if (t > 1)
            dt |= (d[t - 2] & 0xFF);

        while (n-- >= t) {
            var z = (rn << 16) | ((r[n] & 0xFF) << 8);
            if (n > 0)
                z |= (r[n - 1] & 0xFF);

            var i = n - t + 1;
            z /= dt;
            rn += mula_small(r, r, i, d, t, -z);
            q[i] = (z + rn) & 0xFF;
            /* rn is 0 or -1 (underflow) */
            mula_small(r, r, i, d, t, -rn);
            rn = r[n] & 0xFF;
            r[n] = 0;
        }

        r[t-1] = rn & 0xFF;
    }
29481
def smallest_prime_factor(Q):
    """"""

    A = Q;
    while(A != 1):
        if(np.mod(A, 2) == 0):
            A = A / 2
        elif(np.mod(A, 3) == 0):
            A = A / 3
        elif(np.mod(A, 5) == 0):
            A = A / 5
        elif(np.mod(A, 7) == 0):
            A = A / 7;
        else:
            A = Q + 1;
            Q = A;

    return Q
1006809
function QuadSD_ak1(NN, u, v, p, q, iPar){
                   // Divides p by the quadratic 1, u, v placing the quotient in q and the remainder in a, b
                   // iPar is a dummy variable for passing in the two parameters--a and b--by reference
                   q[0] = iPar.b = p[0];
                   q[1] = iPar.a = -(u*iPar.b) + p[1];

                   for (var i = 2; i < NN; i++){
                       q[i] = -(u*iPar.a + v*iPar.b) + p[i];
                       iPar.b = iPar.a;
                       iPar.a = q[i];
                   } 
                   return;
               }
917514
def mod(ctx, number, divisor):
    
    number = conversions.to_decimal(number, ctx)
    divisor = conversions.to_decimal(divisor, ctx)
    return number - divisor * _int(ctx, number / divisor)
515943
def _modulo(self, decimal_argument):
        
        _times, remainder = self._context.divmod(decimal_argument, 100)

        # match the builtin % behavior by adding the N to the result if negative
        return remainder if remainder >= 0 else remainder + 100
403044
function(a, b) {

      var z = new Complex(a, b);

      a = this['re'];
      b = this['im'];

      if (z['isZero']()) {
        return Complex['ONE'];
      }

      // If the exponent is real
      if (z['im'] === 0) {

        if (b === 0 && a >= 0) {

          return new Complex(Math.pow(a, z['re']), 0);

        } else if (a === 0) { // If base is fully imaginary

          switch ((z['re'] % 4 + 4) % 4) {
            case 0:
              return new Complex(Math.pow(b, z['re']), 0);
            case 1:
              return new Complex(0, Math.pow(b, z['re']));
            case 2:
              return new Complex(-Math.pow(b, z['re']), 0);
            case 3:
              return new Complex(0, -Math.pow(b, z['re']));
          }
        }
      }

      /* I couldn't find a good formula, so here is a derivation and optimization
       *
       * z_1^z_2 = (a + bi)^(c + di)
       *         = exp((c + di) * log(a + bi)
       *         = pow(a^2 + b^2, (c + di) / 2) * exp(i(c + di)atan2(b, a))
       * =>...
       * Re = (pow(a^2 + b^2, c / 2) * exp(-d * atan2(b, a))) * cos(d * log(a^2 + b^2) / 2 + c * atan2(b, a))
       * Im = (pow(a^2 + b^2, c / 2) * exp(-d * atan2(b, a))) * sin(d * log(a^2 + b^2) / 2 + c * atan2(b, a))
       *
       * =>...
       * Re = exp(c * log(sqrt(a^2 + b^2)) - d * atan2(b, a)) * cos(d * log(sqrt(a^2 + b^2)) + c * atan2(b, a))
       * Im = exp(c * log(sqrt(a^2 + b^2)) - d * atan2(b, a)) * sin(d * log(sqrt(a^2 + b^2)) + c * atan2(b, a))
       *
       * =>
       * Re = exp(c * logsq2 - d * arg(z_1)) * cos(d * logsq2 + c * arg(z_1))
       * Im = exp(c * logsq2 - d * arg(z_1)) * sin(d * logsq2 + c * arg(z_1))
       *
       */

      if (a === 0 && b === 0 && z['re'] > 0 && z['im'] >= 0) {
        return Complex['ZERO'];
      }

      var arg = Math.atan2(b, a);
      var loh = logHypot(a, b);

      a = Math.exp(z['re'] * loh - z['im'] * arg);
      b = z['im'] * loh + z['re'] * arg;
      return new Complex(
              a * Math.cos(b),
              a * Math.sin(b));
    }
748807
def ceiling(value, mod=1):
    
    if value == None:
        return None
    mod = int(mod)

    v = int(math_floor(value + mod))
    return v - (v % mod)
441123
def symbolic_part(z, s):
    """
    if s == 1: return symre(z)
    elif s == -1: return symim(z)
    elif s == 0: return z
938507
private function GCD($e, $m)
    {
        $y = $e;
        $x = $m;

        while (bccomp($y, 0) != 0) {
            // modulus function
            $w = bcsub($x, bcmul($y, bcdiv($x, $y, 0)));
            $x = $y;
            $y = $w;
        }

        return $x;
    }
866529
def dBinaryRochedy (r, D, q, F):
    
    return -r[1]*(r[0]*r[0]+r[1]*r[1]+r[2]*r[2])**-1.5 -q*r[1]*((r[0]-D)*(r[0]-D)+r[1]*r[1]+r[2]*r[2])**-1.5 + F*F*(1+q)*r[1]