问:
Fastest semiprime factorization
Write a program to factorize a semi-prime number in the shortest amount of time.

For testing purposes, use this:
\$38!+1\$ (\$523022617466601111760007224100074291200000001\$)

It is equal to:
\$14029308060317546154181 × 37280713718589679646221\$

答:
Python (w/ pypy2 v7.3.1) ~0.9s

Using a Multiple Polynomial Quadratic Sieve. I took this to be a code challenge, so I opted not to use any external libraries (other than the standard log function, I suppose). When timing, the PyPy JIT should be used, as it results in timings 4-5 times faster than that of cPython.

Update (2013-07-29):
Since originally posting, I've made several minor, but significant changes which increase the overall speed by a factor of about 2.5x.

Update (2014-08-27):
As this post is still receiving attention, I've updated my_math.py correcting two errors, for anyone who may be using it:

isqrt was faulty, sometimes producing incorrect output for values very close to a perfect square. This has been corrected, and the performance increased by using a much better seed.

is_prime has been updated. My previous attempt to remove perfect square 2-sprps was half-hearted, at best. I've added a 3-sprp check - a technique used by Mathmatica - to ensure that the tested value is square-free.

Update (2014-11-24):
If at the end of the calculation no non-trivial congruencies are found, the progam now sieves additional polynomials. This was previously marked in the code as TODO.

Update (2020-08-29):
A few improvements resulting in minor performance gains, including:

Calculating the average contributation of the skipped primes, rather than roughly approximating.

Actually sieving powers of the prime factor base. The modular square root of each power can be computed efficiently through use of Hensel Lifting.

Reusing the sieve array in subsequent iterations, rather than creating a new one.

mpqs.py

import math
import my_math
import time

# Multiple Polynomial Quadratic Sieve
# assumes n is composite
def mpqs(n, verbose=False):
  if verbose:
    time1 = time.time()

  root_n = my_math.isqrt(n)
  root_2n = my_math.isqrt(n+n)

  # formula chosen by experimentation
  # seems to be close to optimal for n < 10^50
  bound = int(7.5 * math.log(n, 10)**2)

  prime = []
  mod_root = []
  log_p = []
  num_prime = 0

  # size of the sieve
  x_max = bound * 5
  x_max_2 = x_max+x_max

  # maximum value on the sieved range
  m_val = (x_max * root_2n) >> 1

  # find a number of small primes for which n is a quadratic residue
  pg = my_math.primes()
  p = pg.next()
  while p < bound or num_prime < 3:

    # legendre (n|p) is only defined for odd p
    if p > 2:
      leg = my_math.legendre(n, p)
    else:
      leg = n & 1

    if leg == 1:
      prime += [p]
      log_p += [math.log(p, 10)]
      r = int(my_math.mod_sqrt(n, p))
      roots = [r]
      q = p
      while q < x_max:
        # find all square roots mod p^n via Hensel Lifting
        r = int((r + (n - r*r)*my_math.mod_inv(r+r, q))%q)
        #assert r*r%q == n%q
        roots += [r]
        q *= p
      mod_root += [roots]
      num_prime += 1
    elif leg == 0:
      if verbose:
        print 'trial division found factors:'
        print p, 'x', n/p
      return p

    p = pg.next()

  # fudging the threshold down a bit makes it easier to find powers of partial-partial
  # relationships, but it also makes the smoothness check slower. reducing by twice the log
  # of the largest prime in the factor base results in cofactors less than that value squared
  thresh = math.log(m_val, 10) - log_p[-1]*2

  # skip small primes. they contribute very little to the log sum
  # and add a lot of unnecessary entries to the table
  # instead, fudge the threshold down a bit, according to expected number of factors
  min_prime = int(thresh*2)
  sp_idx = my_math.binary_search(prime, min_prime)
  sieve_primes = prime[sp_idx:]

  fudge = sum([log_p[i]/(prime[i]-1) for i in xrange(sp_idx)])

  sums = [fudge]*x_max_2

  if verbose:
    print 'smoothness bound:', bound
    print 'sieve size:', x_max
    print 'log threshold:', thresh
    print 'skipping primes less than:', min_prime

  smooth = []
  used_prime = set()
  partial = {}
  num_smooth = 0
  prev_num_smooth = 0
  num_used_prime = 0
  num_partial = 0
  num_poly = 0
  root_A = my_math.isqrt(root_2n / x_max)

  if verbose:
    print 'sieving for smooths...'
  while True:
    # find an integer value A such that:
    # A is =~ sqrt(2*n) / x_max
    # A is a perfect square
    # sqrt(A) is prime, and n is a quadratic residue mod sqrt(A)
    while True:
      root_A = my_math.next_prime(root_A)
      leg = my_math.legendre(n, root_A)
      if leg == 1:
        break
      elif leg == 0:
        if verbose:
          print 'dumb luck found factors:'
          print root_A, 'x', n/root_A
        return root_A

    A = root_A * root_A

    # solve for an adequate B
    # B*B is a quadratic residue mod n, such that B*B-A*C = n
    # this is unsolvable if n is not a quadratic residue mod sqrt(A)
    b = my_math.mod_sqrt(n, root_A)
    B = (b + (n - b*b) * my_math.mod_inv(b + b, root_A))%A

    # B*B-A*C = n <=> C = (B*B-n)/A
    C = (B*B - n) / A

    num_poly += 1

    # sieve for prime factors
    i = sp_idx
    for p in sieve_primes:
      logp = log_p[i]

      e = 0
      q = p
      while q < x_max:
        inv_A = my_math.mod_inv(A, q)
        # modular root of the quadratic
        a = int(((mod_root[i][e] - B) * inv_A)%q)
        b = int(((q - mod_root[i][e] - B) * inv_A)%q)

        amx = a+x_max
        bmx = b+x_max

        apx = amx-q
        bpx = bmx-q

        k = q
        while k < x_max:
          sums[apx+k] += logp
          sums[bpx+k] += logp
          sums[amx-k] += logp
          sums[bmx-k] += logp
          k += q

        q *= p
        e += 1

      i += 1

    # check for smooths
    x = -x_max
    i = 0
    while i < x_max_2:
      v = sums[i]
      if v > thresh:
        vec = set()
        sqr = []
        # because B*B-n = A*C
        # (A*x+B)^2 - n = A*A*x*x+2*A*B*x + B*B - n
        #               = A*(A*x*x+2*B*x+C)
        # gives the congruency
        # (A*x+B)^2 = A*(A*x*x+2*B*x+C) (mod n)
        # because A is chosen to be square, it doesn't need to be sieved
        sieve_val = (A*x + B+B)*x + C

        if sieve_val < 0:
          vec = {-1}
          sieve_val = -sieve_val

        for p in prime:
          while sieve_val%p == 0:
            if p in vec:
              # keep track of perfect square factors
              # to avoid taking the sqrt of a gigantic number at the end
              sqr += [p]
            vec ^= {p}
            sieve_val = int(sieve_val / p)

        if sieve_val == 1:
          # smooth
          smooth += [(vec, (sqr, (A*x+B), root_A))]
          used_prime |= vec
        elif sieve_val in partial:
          # combine two partials to make a (xor) smooth
          # that is, every prime factor with an odd power is in our factor base
          pair_vec, pair_vals = partial[sieve_val]
          sqr += list(vec & pair_vec) + [sieve_val]
          vec ^= pair_vec
          smooth += [(vec, (sqr + pair_vals[0], (A*x+B)*pair_vals[1], root_A*pair_vals[2]))]
          used_prime |= vec
          num_partial += 1
        else:
          # save partial for later pairing
          partial[sieve_val] = (vec, (sqr, A*x+B, root_A))
      x += 1

      # reset the value for the next go
      sums[i] = fudge
      i += 1

    prev_num_smooth = num_smooth
    num_smooth = len(smooth)
    num_used_prime = len(used_prime)

    if verbose:
      print 100 * num_smooth / num_prime, 'percent complete\r',

    if num_smooth > num_used_prime and num_smooth > prev_num_smooth:
      if verbose:
        print '%d polynomials sieved (%d values)'%(num_poly, num_poly*x_max_2)
        print 'found %d smooths (%d from partials) in %.3f seconds'%(num_smooth, num_partial, time.time()-time1)
        print 'solving for non-trivial congruencies...'

      # set up bit fields for gaussian elimination
      masks = []
      mask = 1
      bit_fields = [0]*num_used_prime
      for vec, vals in smooth:
        masks += [mask]
        i = 0
        for p in used_prime:
          if p in vec: bit_fields[i] |= mask
          i += 1
        mask += mask

      # row echelon form
      col_offset = 0
      null_cols = []
      for col in xrange(num_smooth):
        pivot = col-col_offset == num_used_prime or bit_fields[col-col_offset] & masks[col] == 0
        for row in xrange(col+1-col_offset, num_used_prime):
          if bit_fields[row] & masks[col]:
            if pivot:
              bit_fields[col-col_offset], bit_fields[row] = bit_fields[row], bit_fields[col-col_offset]
              pivot = False
            else:
              bit_fields[row] ^= bit_fields[col-col_offset]
        if pivot:
          null_cols += [col]
          col_offset += 1

      # reduced row echelon form
      for row in xrange(num_used_prime):
        # lowest set bit
        mask = bit_fields[row] & -bit_fields[row]
        for up_row in xrange(row):
          if bit_fields[up_row] & mask:
            bit_fields[up_row] ^= bit_fields[row]

      # check for non-trivial congruencies
      for col in null_cols:
        all_vec, (lh, rh, rA) = smooth[col]
        lhs = lh   # sieved values (left hand side)
        rhs = [rh] # sieved values - n (right hand side)
        rAs = [rA] # root_As (cofactor of lhs)
        i = 0
        for field in bit_fields:
          if field & masks[col]:
            vec, (lh, rh, rA) = smooth[i]
            lhs += list(all_vec & vec) + lh
            all_vec ^= vec
            rhs += [rh]
            rAs += [rA]
          i += 1

        factor = my_math.gcd(my_math.list_prod(rAs)*my_math.list_prod(lhs) - my_math.list_prod(rhs), n)
        if 1 < factor < n:
          break
      else:
        if verbose:
          print 'none found.'
        continue
      break

  if verbose:
    print 'factors found:'
    print factor, 'x', n/factor
    print 'time elapsed: %.3f seconds'%(time.time()-time1)
  return factor

if __name__ == "__main__":
  import argparse
  parser = argparse.ArgumentParser(description='Uses a MPQS to factor a composite number')
  parser.add_argument('composite', metavar='number_to_factor', type=long, help='the composite number to factor')
  parser.add_argument('--verbose', dest='verbose', action='store_true', help="enable verbose output")
  args = parser.parse_args()

  if args.verbose:
    mpqs(args.composite, args.verbose)
  else:
    time1 = time.time()
    print mpqs(args.composite)
    print 'time elapsed: %.3f seconds'%(time.time()-time1)
```
<p><strong>my_math.py</strong></p>
<pre class="lang-python prettyprint-override"><code>import math
import mpqs

# primes less than 212
small_primes = [
    2,  3,  5,  7, 11, 13, 17, 19, 23, 29, 31, 37,
   41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89,
   97,101,103,107,109,113,127,131,137,139,149,151,
  157,163,167,173,179,181,191,193,197,199,211]

# pre-calced sieve of eratosthenes for n = 2, 3, 5, 7
indices = [
    1, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47,
   53, 59, 61, 67, 71, 73, 79, 83, 89, 97,101,103,
  107,109,113,121,127,131,137,139,143,149,151,157,
  163,167,169,173,179,181,187,191,193,197,199,209]

# distances between sieve values
offsets = [
  10, 2, 4, 2, 4, 6, 2, 6, 4, 2, 4, 6,
   6, 2, 6, 4, 2, 6, 4, 6, 8, 4, 2, 4,
   2, 4, 8, 6, 4, 6, 2, 4, 6, 2, 6, 6,
   4, 2, 4, 6, 2, 6, 4, 2, 4, 2,10, 2]

# tabulated, mod 105
dindices =[
  0,10, 2, 0, 4, 0, 0, 0, 8, 0, 0, 2, 0, 4, 0,
  0, 6, 2, 0, 4, 0, 0, 4, 6, 0, 0, 6, 0, 0, 2,
  0, 6, 2, 0, 4, 0, 0, 4, 6, 0, 0, 2, 0, 4, 2,
  0, 6, 6, 0, 0, 0, 0, 6, 6, 0, 0, 0, 0, 4, 2,
  0, 6, 2, 0, 4, 0, 0, 4, 6, 0, 0, 2, 0, 6, 2,
  0, 6, 0, 0, 4, 0, 0, 4, 6, 0, 0, 2, 0, 4, 8,
  0, 0, 2, 0,10, 0, 0, 4, 0, 0, 0, 2, 0, 4, 2]

max_int = 2147483647

# returns the index of x in a sorted list a
# or the index of the next larger item if x is not present
# i.e. the proper insertion point for x in a
def binary_search(a, x):
  s = 0
  e = len(a)
  m = e >> 1
  while m != e:
    if a[m] < x:
      s = m
      m = (s + e + 1) >> 1
    else:
      e = m
      m = (s + e) >> 1
  return m

# divide and conquer list product
def list_prod(a):
  while len(a) > 1:
    a = [m*n for m,n in zip(a[::2], a[1::2] + [1])]
  return a[0]

# greatest common divisor of a and b
def gcd(a, b):
  while b:
    a, b = b, a%b
  return a

# extended gcd
def ext_gcd(a, m):
  a = int(a%m)
  x, u = 0, 1
  while a:
    x, u = u, x - (m/a)*u
    m, a = a, m%a
  return (m, x, u)

# legendre symbol (a|m)
# note: returns m-1 if a is a non-residue, instead of -1
def legendre(a, m):
  return pow(a, (m-1) >> 1, m)

# modular inverse of a mod m
def mod_inv(a, m):
  return ext_gcd(a, m)[1]

# modular sqrt(n) mod p
# p must be prime
def mod_sqrt(n, p):
  a = n%p
  if p%4 == 3:
    return pow(a, (p+1) >> 2, p)
  elif p%8 == 5:
    v = pow(a << 1, (p-5) >> 3, p)
    i = ((a*v*v << 1) % p) - 1
    return (a*v*i)%p
  elif p%8 == 1:
    # Shanks' method
    q = p-1
    e = 0
    while q&1 == 0:
      e += 1
      q >>= 1

    n = 2
    while legendre(n, p) != p-1:
      n += 1

    w = pow(a, q, p)
    x = pow(a, (q+1) >> 1, p)
    y = pow(n, q, p)
    r = e
    while True:
      if w == 1:
        return x

      v = w
      k = 0
      while v != 1 and k+1 < r:
        v = (v*v)%p
        k += 1

      if k == 0:
        return x

      d = pow(y, 1 << (r-k-1), p)
      x = (x*d)%p
      y = (d*d)%p
      w = (w*y)%p
      r = k
  else: # p == 2
    return a

#integer sqrt of n
def isqrt(n):
  c = (n << 2)/3
  d = c.bit_length()

  a = d>>1
  if d&1:
    x = 1 << a
    y = (x + (n >> a)) >> 1
  else:
    x = (3 << a) >> 2
    y = (x + (c >> a)) >> 1

  if x != y:
    x = y
    y = (x + n/x) >> 1
    while y < x:
      x = y
      y = (x + n/x) >> 1
  return x

# integer cbrt of n
def icbrt(n):
  d = n.bit_length()

  if d%3 == 2:
    x = 3 << d/3-1
  else:
    x = 1 << d/3

  y = (2*x + n/(x*x))/3
  if x != y:
    x = y
    y = (2*x + n/(x*x))/3
    while y < x:
      x = y
      y = (2*x + n/(x*x))/3
  return x

# strong probable prime
def is_sprp(n, b=2):
  if n < 2: return False
  d = n-1
  s = 0
  while d&1 == 0:
    s += 1
    d >>= 1

  x = pow(b, d, n)
  if x == 1 or x == n-1:
    return True

  for r in xrange(1, s):
    x = (x * x)%n
    if x == 1:
      return False
    elif x == n-1:
      return True

  return False

# lucas probable prime
# assumes D = 1 (mod 4), (D|n) = -1
def is_lucas_prp(n, D):
  P = 1
  Q = (1-D) >> 2

  # n+1 = 2**r*s where s is odd
  s = n+1
  r = 0
  while s&1 == 0:
    r += 1
    s >>= 1

  # calculate the bit reversal of (odd) s
  # e.g. 19 (10011) <=> 25 (11001)
  t = 0
  while s:
    if s&1:
      t += 1
      s -= 1
    else:
      t <<= 1
      s >>= 1

  # use the same bit reversal process to calculate the sth Lucas number
  # keep track of q = Q**n as we go
  U = 0
  V = 2
  q = 1
  # mod_inv(2, n)
  inv_2 = (n+1) >> 1
  while t:
    if t&1:
      # U, V of n+1
      U, V = ((U + V) * inv_2)%n, ((D*U + V) * inv_2)%n
      q = (q * Q)%n
      t -= 1
    else:
      # U, V of n*2
      U, V = (U * V)%n, (V * V - 2 * q)%n
      q = (q * q)%n
      t >>= 1

  # double s until we have the 2**r*sth Lucas number
  while r:
    U, V = (U * V)%n, (V * V - 2 * q)%n
    q = (q * q)%n
    r -= 1

  # primality check
  # if n is prime, n divides the n+1st Lucas number, given the assumptions
  return U == 0

## Baillie-PSW ##
# this is technically a probabalistic test, but there are no known pseudoprimes
def is_bpsw(n):
  if not is_sprp(n, 2): return False

  # idea shamelessly stolen from Mathmatica's PrimeQ
  # if n is a 2-sprp and a 3-sprp, n is necessarily square-free
  if not is_sprp(n, 3): return False

  a = 5
  s = 2
  # if n is a perfect square, this will never terminate
  while legendre(a, n) != n-1:
    s = -s
    a = s-a
  return is_lucas_prp(n, a)

# an 'almost certain' primality check
def is_prime(n):
  if n < 212:
    m = binary_search(small_primes, n)
    return n == small_primes[m]

  for p in small_primes:
    if n%p == 0:
      return False

  # if n is a 32-bit integer, perform full trial division
  if n <= max_int:
    p = 211
    while p*p < n:
      for o in offsets:
        p += o
        if n%p == 0:
          return False
    return True

  return is_bpsw(n)

# next prime strictly larger than n
def next_prime(n):
  if n < 2:
    return 2

  n += 1
  if n < 212:
    m = binary_search(small_primes, n)
    return small_primes[m]

  # find our position in the sieve rotation via binary search
  x = int(n%210)
  m = binary_search(indices, x)
  i = int(n + (indices[m] - x))

  # adjust offsets
  offs = offsets[m:] + offsets[:m]
  while True:
    for o in offs:
      if is_prime(i):
        return i
      i += o

# an infinite prime number generator
def primes(start = 0):
  for n in small_primes[start:]: yield n
  pg = primes(6)
  p = pg.next()
  q = p*p
  sieve = {221: 13, 253: 11}
  n = 211
  while True:
    for o in offsets:
      n += o
      stp = sieve.pop(n, 0)
      if stp:
        nxt = n/stp
        nxt += dindices[nxt%105]
        while nxt*stp in sieve:
          nxt += dindices[nxt%105]
        sieve[nxt*stp] = stp
      elif n < q:
        yield n
      else:
        sieve[q + dindices[p%105]*p] = p
        p = pg.next()
        q = p*p

# true if n is a prime power > 0
def is_prime_power(n):
  if n > 1:
    for p in small_primes:
      if n%p == 0:
        n /= p
        while n%p == 0: n /= p
        return n == 1

    r = isqrt(n)
    if r*r == n:
      return is_prime_power(r)

    s = icbrt(n)
    if s*s*s == n:
      return is_prime_power(s)

    p = 211
    while p*p < r:
      for o in offsets:
        p += o
        if n%p == 0:
          n /= p
          while n%p == 0: n /= p
          return n == 1

    if n <= max_int:
      while p*p < n:
        for o in offsets:
          p += o
          if n%p == 0:
            return False
      return True

    return is_bpsw(n)
  return False

max_trial = 1e10
max_pollard = 1e22

def factor(n):
  if n < max_trial:
    return factor_trial(n)
  for p in small_primes:
    if n%p == 0:
      return [p] + factor(n/p)
  if is_prime(n):
    return [n]
  if n < max_pollard:
    p = pollard_rho(n)
  else:
    p = lenstra_ecf(n) or mpqs.mpqs(n)
  return factor(p) + factor(n/p)

def factor_trial(n):
  a = []
  for p in small_primes:
    while n%p == 0:
      a += [p]
      n /= p
  i = 211
  while i*i < n:
    for o in offsets:
      i += o
      while n%i == 0:
        a += [i]
        n /= i
  if n > 1:
    a += [n]
  return a

def pollard_rho(n):
  # Brent's variant
  y, r, q = 0, 1, 1
  c, m = 9, 40
  g = 1
  while g == 1:
    x = y
    for i in range(r):
      y = (y*y + c) % n
    k = 0
    while k < r and g == 1:
      ys = y
      for j in range(min(m, r-k)):
        y = (y*y + c) % n
        q = q*abs(x-y) % n
      g = gcd(q, n)
      k += m
    r *= 2
  if g == n:
    ys = (ys*ys + c) % n
    g = gcd(n, abs(x-ys))
    while g == 1:
      ys = (ys*ys + c) % n
      g = gcd(n, abs(x-ys))
  return g

def ec_add((x1, z1), (x2, z2), (x0, z0), n):
  t1, t2 = (x1-z1)*(x2+z2), (x1+z1)*(x2-z2)
  x, z = t1+t2, t1-t2
  return (z0*x*x % n, x0*z*z % n)

def ec_double((x, z), (a, b), n):
  t1 = x+z; t1 *= t1
  t2 = x-z; t2 *= t2
  t3 = t1 - t2
  t4 = 4*b*t2
  return (t1*t4 % n, t3*(t4 + a*t3) % n)

def ec_multiply(k, p, C, n):
  # Montgomery ladder algorithm
  p0 = p
  q, p = p, ec_double(p, C, n)
  b = k >> 1
  while b > (b & -b):
    b ^= b & -b
  while b:
    if k&b:
      q, p = ec_add(p, q, p0, n), ec_double(p, C, n)
    else:
      q, p = ec_double(q, C, n), ec_add(p, q, p0, n),
    b >>= 1
  return q

def lenstra_ecf(n, m = 5):
  # Montgomery curves w/ Suyama parameterization.
  # Based on pseudocode found in:
  # "Implementing the Elliptic Curve Method of Factoring in Reconfigurable Hardware"
  # Gaj, Kris et. al
  # http://www.hyperelliptic.org/tanja/SHARCS/talks06/Gaj.pdf
  # Phase 2 is not implemented.
  B1, B2 = 8, 13
  for i in range(m):
    pg = primes()
    p = pg.next()
    k = 1
    while p < B1:
      k *= p**int(math.log(B1, p))
      p = pg.next()
    for s in range(B1, B2):
      u, v = s*s-5, 4*s
      x = u*u*u
      z = v*v*v
      t = pow(v-u, 3, n)
      P = (x, z)
      C = (t*(3*u+v) % n, 4*x*v % n)
      Q = ec_multiply(k, P, C, n)
      g = gcd(Q[1], n)
      if 1 < g < n: return g
    B1, B2 = B2, B1 + B2

if __name__ == "__main__":
  import argparse
  import time
  parser = argparse.ArgumentParser(description='Uses a various methods to factor a composite number')
  parser.add_argument('composite', metavar='number_to_factor', type=long, help='the composite number to factor')
  args = parser.parse_args()

  time1 = time.time()
  print factor(args.composite)
  print 'time elapsed: %.3f seconds'%(time.time()-time1)
</code>```

Sample I/O:

$ pypy mpqs.py --verbose 523022617466601111760007224100074291200000001
smoothness bound: 14998
sieve size: 74990
log threshold: 18.7325510316
skipping primes less than: 37
sieving for smooths...
230 polynomials sieved (34495400 values)
found 849 smooths (242 from partials) in 0.786 seconds
solving for non-trivial congruencies...
factors found:
37280713718589679646221 x 14029308060317546154181
time elapsed: 1.011 seconds
```
<p>Note: not using the <code>--verbose</code> option will give slightly better timings:</p>
<pre class="lang-python prettyprint-override"><code>$ pypy mpqs.py 523022617466601111760007224100074291200000001
37280713718589679646221
time elapsed: 0.902 seconds
</code>```

Basic Concepts

In general, a quadratic sieve is based on the following observation: any odd composite \$n\$ may be represented as:

\$n=(x+d)(x-d)=x^2-d^2\Rightarrow d^2=x^2-n\$

This is not very difficult to confirm. Since \$n\$ is odd, the distance between any two cofactors of \$n\$ must be even \$2d\$, where \$x\$ is the mid point between them. Moreover, the same relation holds for any multiple of \$n\$

\$abn=(ax+ad)(bx-bd)=abx^2-abd^2\Rightarrow abd^2=abx^2-abn\$

Note that if any such \$x\$ and \$d\$ can be found, it will immediately result in a (not necessarily prime) factor of \$n\$, since \$x+d\$ and \$x-d\$ both divide \$n\$ by definition. This relation can be further weakened - at the consequence of allowing potential trivial congruencies - to the following form:

\$d^2\equiv x^2\mod n\$

So in general, if we can find two perfect squares which are equivalent mod \$n\$, then it's fairly likely that we can directly produce a factor of \$n\$ a la \$\gcd(x±d,n)\$. Seems pretty simple, right?

Except it's not. If we intended to conduct an exhaustive search over all possible \$x\$, we would need to search the entire range from \$\sqrt{n}\$ to \$\sqrt{2n}\$, which is marginally smaller than full trial division, but also requires an expensive is_square operation each iteration to confirm the value of \$d\$. Unless it is known beforehand that \$n\$ has factors very near \$\sqrt{n}\$, trial division is likely to be faster.

Perhaps we can weaken this relation even more. Suppose we chose an \$x\$, such that for

\$y\equiv x^2\mod n\$

a full prime factorization of \$y\$ is readily known. If we had enough such relations, we should be able to construct an adequate \$d\$, if we choose a number of \$y\$ such that their product is a perfect square; that is, all prime factors are used an even number of times. In fact, if we have more such \$y\$ than the total number of unique prime factors they contain, a solution is guaranteed to exist; It becomes a system of linear equations. The question now becomes, how do we chose such \$x\$? That's where sieving comes into play.

The Sieve

Consider the polynomial:

\$y(x)=x^2-n\$

Then for any prime \$p\$ and integer \$k\$, the following is true:

\$y(x+kp)=(x+kp)^2-n\\y(x+kp)=x^2+2xkp+(kp)^2-n\\y(x+kp)=y(x)+2xkp+(kp)^2\equiv y(x)\mod p\$

This means that after solving for the roots of the polynomial mod \$p\$ - that is, you've found an \$x\$ such that \$y(x)\equiv 0\mod p\$, ergo \$y\$ is divisible by \$p\$ - then you have found an infinite number of such \$x\$. In this way, you can sieve over a range of \$x\$, identifying small prime factors of \$y\$, hopefully finding some for which all prime factors are small. Such numbers known as \$k-smooth\$, where \$k\$ is the largest prime factor used.

There's a few problems with this approach, though. Not all values of \$x\$ are adequate, in fact, there's only very few of them, centered around \$\sqrt{n}\$. Smaller values will become largely negative (due to the \$-n\$ term), and larger values will become too large, such that it is unlikely that their prime factorization consists only of small primes. There will be a number of such \$x\$, but unless the composite you're factoring is very small, it's highly unlikely that you'll find enough smooths to result in a factorization. And so, for larger \$n\$, it becomes necessary to sieve over multiple polynomials of a given form.

Multiple Polynomials

So we need more polynomials to sieve? How about this:

\$y(x)=(Ax+B)^2-n\$

That'll work. Note that \$A\$ and \$B\$ could literally be any integer value, and the math still holds. All we need to do is choose a few random values, solve for the root of the polynomial, and sieve the values close to zero. At this point we could just call it good enough: if you throw enough stones in random directions, you're bound to break a window sooner or later.

Except, there's a problem with that too. If the slope of the polynomial is large at the x-intercept, there'll only be a few suitable values to sieve per polynomial. It'll work, but you'll end up sieving a whole lot of polynomials before you get what you need. Can we do better?

We can do better. An observation, as a result of Montgomery is as follows: if \$A\$ and \$B\$ are chosen such that there exists some \$C\$ satisfying

\$B^2-n=AC\$

then the entire polynomial can be rewritten as

\$y(x)=(Ax+B)^2-n=(Ax)^2+2ABx+B^2-n=A(Ax^2+2Bx+C)\$

Furthermore, if \$A\$ is chosen to be a perfect square, the leading \$A\$ term can be neglected while sieving, resulting in much smaller values, and a much flatter curve. For such a solution to exist, \$n\$ must be a quadratic residue mod \$\sqrt{A}\$, which can be known immediately by computing the Legendre symbol:
\$(n|\sqrt{A})=1\$. Note that in order to solve for \$B\$, a complete prime factorization of \$\sqrt{A}\$ needs to be known (in order to take the modular square root \$\sqrt{n}\mod\sqrt{A}\$), which is why \$\sqrt{A}\$ is typically chosen to be prime.

It can then be shown that if \$A\approx\frac{\sqrt{2n}}{M}\$, then for all values of \$x\in[-M,M]\$:

\$|y(x)|\le\frac{M\sqrt{2n}}{2}\$

And now, finally, we have all the components necessary to implement our sieve. Or do we?

Powers of Primes as Factors

Our sieve, as described above, has one major flaw. It can identify which values of \$x\$ will result in a \$y\$ divisible by \$p\$, but it cannot identify whether or not this \$y\$ is divisible by a power of \$p\$. In order to determine that, we would need to perform trial division on the value to be sieved, until it is no longer divisible by \$p\$.

Edit: This is incorrect. The roots of \$n\$ mod \$p^k\$ can be computed directly from the roots mod \$p\$ through the use of Hensel Lifting. This current implementation does precisely this.

We seemed to have reached an impasse: the whole point of the sieve was so that we didn't have to do that. Time to check the playbook.

\$\ln(a\cdot b\cdot c\cdot d\cdot\ldots)=\ln(a)+\ln(b)+\ln(c)+\ln(d)+\ln(\ldots)\$

That looks pretty useful. If the sum of the \$\ln\$ of all of the small prime factors of \$y\$ is close to the expected value of \$\ln(y)\$, then it's almost a given that \$y\$ has no other factors. In addition, if we adjust the expected value down a little bit, we can also identify values as smooth which have several powers of primes as factors. In this way, we can use the sieve as a 'pre-screening' process, and only factor those values which are likely to be smooth.

This has a few other advantages as well. Note that small primes contribute very little to the \$\ln\$ sum, but yet they require the most sieve time. Sieving the value 3 requires more time than 11, 13, 17, 19, and 23 combined. Instead, we can just skip the first few primes, and adjust the threshold down accordingly, assuming a certain percentage of them would have passed.

Another result, is that a number of values will be allowed to 'slip through', which are mostly smooth, but contain a single large cofactor. We could just discard these values, but suppose we found another mostly smooth value, with exactly the same cofactor. We can then use these two values to construct a usable \$y\$; since their product will contain this large cofactor squared, it no longer needs to be considered.

Putting it all together

The last thing we need to do is to use these values of \$y\$ construct an adequate \$x\$ and \$d\$. Suppose we only consider the non-square factors of \$y\$, that is, the prime factors of an odd power. Then, each \$y\$ can be expressed in the following manner:

\$y_0=p_0^0\cdot p_1^1\cdot p_2^1\cdots p_n^0\\y_1=p_0^1\cdot p_1^0\cdot p_2^1\cdots p_n^1\\y_2=p_0^0\cdot p_1^0\cdot p_2^0\cdots p_n^1\\y_3=p_0^1\cdot p_1^1\cdot p_2^0\cdots p_n^0\\\vdots\$

which can be expressed in the matrix form:

\$M=\begin{bmatrix}0&1&1&\cdots&0\\1&0&1&\cdots&1\\0&0&0&\cdots&1\\1&1&0&\cdots&0\\\vdots\end{bmatrix}\$

The problem then becomes to find a vector \$v\$ such that \$vM=\vec{0}\mod 2\$, where \$\vec{0}\$ is the null vector. That is, to solve for the left null space of \$M\$. This can be done in a number of ways, the simplest of which is to perform Gaussian Elimination on \$M^T\$, replacing the row addition operation with a row xor. This will result in a number of null space basis vectors, any combination of which will produce a valid solution.

The construction of \$x\$ is fairly straight-forward. It is simply the product of \$Ax+B\$ for each of the \$y\$ used. The construction of \$d\$ is slightly more complicated. If we were to take the product of all \$y\$, we will end up with a value with 10s of thousands, if not 100s of thousands of digits, for which we need to find the square root. This calcuation is impractically expensive. Instead, we can keep track of the even powers of primes during the sieving process, and then use and and xor operations on the vectors of non-square factors to reconstruct the square root.

<s>I seem to have reached the 30000 character limit. Ahh well, I suppose that's good enough.</s> Saved a bunch of bytes by switching to \$MathJax\$.

Title: Fastest semiprime factorization

Body:
Python (w/ pypy2 v7.3.1) ~0.9s

Using a Multiple Polynomial Quadratic Sieve. I took this to be a code challenge, so I opted not to use any external libraries (other than the standard log function, I suppose). When timing, the PyPy JIT should be used, as it results in timings 4-5 times faster than that of cPython.

Update (2013-07-29):
Since originally posting, I've made several minor, but significant changes which increase the overall speed by a factor of about 2.5x.

Update (2014-08-27):
As this post is still receiving attention, I've updated my_math.py correcting two errors, for anyone who may be using it:

isqrt was faulty, sometimes producing incorrect output for values very close to a perfect square. This has been corrected, and the performance increased by using a much better seed.

is_prime has been updated. My previous attempt to remove perfect square 2-sprps was half-hearted, at best. I've added a 3-sprp check - a technique used by Mathmatica - to ensure that the tested value is square-free.

Update (2014-11-24):
If at the end of the calculation no non-trivial congruencies are found, the progam now sieves additional polynomials. This was previously marked in the code as TODO.

Update (2020-08-29):
A few improvements resulting in minor performance gains, including:

Calculating the average contributation of the skipped primes, rather than roughly approximating.

Actually sieving powers of the prime factor base. The modular square root of each power can be computed efficiently through use of Hensel Lifting.

Reusing the sieve array in subsequent iterations, rather than creating a new one.

mpqs.py

import math
import my_math
import time

# Multiple Polynomial Quadratic Sieve
# assumes n is composite
def mpqs(n, verbose=False):
  if verbose:
    time1 = time.time()

  root_n = my_math.isqrt(n)
  root_2n = my_math.isqrt(n+n)

  # formula chosen by experimentation
  # seems to be close to optimal for n < 10^50
  bound = int(7.5 * math.log(n, 10)**2)

  prime = []
  mod_root = []
  log_p = []
  num_prime = 0

  # size of the sieve
  x_max = bound * 5
  x_max_2 = x_max+x_max

  # maximum value on the sieved range
  m_val = (x_max * root_2n) >> 1

  # find a number of small primes for which n is a quadratic residue
  pg = my_math.primes()
  p = pg.next()
  while p < bound or num_prime < 3:

    # legendre (n|p) is only defined for odd p
    if p > 2:
      leg = my_math.legendre(n, p)
    else:
      leg = n & 1

    if leg == 1:
      prime += [p]
      log_p += [math.log(p, 10)]
      r = int(my_math.mod_sqrt(n, p))
      roots = [r]
      q = p
      while q < x_max:
        # find all square roots mod p^n via Hensel Lifting
        r = int((r + (n - r*r)*my_math.mod_inv(r+r, q))%q)
        #assert r*r%q == n%q
        roots += [r]
        q *= p
      mod_root += [roots]
      num_prime += 1
    elif leg == 0:
      if verbose:
        print 'trial division found factors:'
        print p, 'x', n/p
      return p

    p = pg.next()

  # fudging the threshold down a bit makes it easier to find powers of partial-partial
  # relationships, but it also makes the smoothness check slower. reducing by twice the log
  # of the largest prime in the factor base results in cofactors less than that value squared
  thresh = math.log(m_val, 10) - log_p[-1]*2

  # skip small primes. they contribute very little to the log sum
  # and add a lot of unnecessary entries to the table
  # instead, fudge the threshold down a bit, according to expected number of factors
  min_prime = int(thresh*2)
  sp_idx = my_math.binary_search(prime, min_prime)
  sieve_primes = prime[sp_idx:]

  fudge = sum([log_p[i]/(prime[i]-1) for i in xrange(sp_idx)])

  sums = [fudge]*x_max_2

  if verbose:
    print 'smoothness bound:', bound
    print 'sieve size:', x_max
    print 'log threshold:', thresh
    print 'skipping primes less than:', min_prime

  smooth = []
  used_prime = set()
  partial = {}
  num_smooth = 0
  prev_num_smooth = 0
  num_used_prime = 0
  num_partial = 0
  num_poly = 0
  root_A = my_math.isqrt(root_2n / x_max)

  if verbose:
    print 'sieving for smooths...'
  while True:
    # find an integer value A such that:
    # A is =~ sqrt(2*n) / x_max
    # A is a perfect square
    # sqrt(A) is prime, and n is a quadratic residue mod sqrt(A)
    while True:
      root_A = my_math.next_prime(root_A)
      leg = my_math.legendre(n, root_A)
      if leg == 1:
        break
      elif leg == 0:
        if verbose:
          print 'dumb luck found factors:'
          print root_A, 'x', n/root_A
        return root_A

    A = root_A * root_A

    # solve for an adequate B
    # B*B is a quadratic residue mod n, such that B*B-A*C = n
    # this is unsolvable if n is not a quadratic residue mod sqrt(A)
    b = my_math.mod_sqrt(n, root_A)
    B = (b + (n - b*b) * my_math.mod_inv(b + b, root_A))%A

    # B*B-A*C = n <=> C = (B*B-n)/A
    C = (B*B - n) / A

    num_poly += 1

    # sieve for prime factors
    i = sp_idx
    for p in sieve_primes:
      logp = log_p[i]

      e = 0
      q = p
      while q < x_max:
        inv_A = my_math.mod_inv(A, q)
        # modular root of the quadratic
        a = int(((mod_root[i][e] - B) * inv_A)%q)
        b = int(((q - mod_root[i][e] - B) * inv_A)%q)

        amx = a+x_max
        bmx = b+x_max

        apx = amx-q
        bpx = bmx-q

        k = q
        while k < x_max:
          sums[apx+k] += logp
          sums[bpx+k] += logp
          sums[amx-k] += logp
          sums[bmx-k] += logp
          k += q

        q *= p
        e += 1

      i += 1

    # check for smooths
    x = -x_max
    i = 0
    while i < x_max_2:
      v = sums[i]
      if v > thresh:
        vec = set()
        sqr = []
        # because B*B-n = A*C
        # (A*x+B)^2 - n = A*A*x*x+2*A*B*x + B*B - n
        #               = A*(A*x*x+2*B*x+C)
        # gives the congruency
        # (A*x+B)^2 = A*(A*x*x+2*B*x+C) (mod n)
        # because A is chosen to be square, it doesn't need to be sieved
        sieve_val = (A*x + B+B)*x + C

        if sieve_val < 0:
          vec = {-1}
          sieve_val = -sieve_val

        for p in prime:
          while sieve_val%p == 0:
            if p in vec:
              # keep track of perfect square factors
              # to avoid taking the sqrt of a gigantic number at the end
              sqr += [p]
            vec ^= {p}
            sieve_val = int(sieve_val / p)

        if sieve_val == 1:
          # smooth
          smooth += [(vec, (sqr, (A*x+B), root_A))]
          used_prime |= vec
        elif sieve_val in partial:
          # combine two partials to make a (xor) smooth
          # that is, every prime factor with an odd power is in our factor base
          pair_vec, pair_vals = partial[sieve_val]
          sqr += list(vec & pair_vec) + [sieve_val]
          vec ^= pair_vec
          smooth += [(vec, (sqr + pair_vals[0], (A*x+B)*pair_vals[1], root_A*pair_vals[2]))]
          used_prime |= vec
          num_partial += 1
        else:
          # save partial for later pairing
          partial[sieve_val] = (vec, (sqr, A*x+B, root_A))
      x += 1

      # reset the value for the next go
      sums[i] = fudge
      i += 1

    prev_num_smooth = num_smooth
    num_smooth = len(smooth)
    num_used_prime = len(used_prime)

    if verbose:
      print 100 * num_smooth / num_prime, 'percent complete\r',

    if num_smooth > num_used_prime and num_smooth > prev_num_smooth:
      if verbose:
        print '%d polynomials sieved (%d values)'%(num_poly, num_poly*x_max_2)
        print 'found %d smooths (%d from partials) in %.3f seconds'%(num_smooth, num_partial, time.time()-time1)
        print 'solving for non-trivial congruencies...'

      # set up bit fields for gaussian elimination
      masks = []
      mask = 1
      bit_fields = [0]*num_used_prime
      for vec, vals in smooth:
        masks += [mask]
        i = 0
        for p in used_prime:
          if p in vec: bit_fields[i] |= mask
          i += 1
        mask += mask

      # row echelon form
      col_offset = 0
      null_cols = []
      for col in xrange(num_smooth):
        pivot = col-col_offset == num_used_prime or bit_fields[col-col_offset] & masks[col] == 0
        for row in xrange(col+1-col_offset, num_used_prime):
          if bit_fields[row] & masks[col]:
            if pivot:
              bit_fields[col-col_offset], bit_fields[row] = bit_fields[row], bit_fields[col-col_offset]
              pivot = False
            else:
              bit_fields[row] ^= bit_fields[col-col_offset]
        if pivot:
          null_cols += [col]
          col_offset += 1

      # reduced row echelon form
      for row in xrange(num_used_prime):
        # lowest set bit
        mask = bit_fields[row] & -bit_fields[row]
        for up_row in xrange(row):
          if bit_fields[up_row] & mask:
            bit_fields[up_row] ^= bit_fields[row]

      # check for non-trivial congruencies
      for col in null_cols:
        all_vec, (lh, rh, rA) = smooth[col]
        lhs = lh   # sieved values (left hand side)
        rhs = [rh] # sieved values - n (right hand side)
        rAs = [rA] # root_As (cofactor of lhs)
        i = 0
        for field in bit_fields:
          if field & masks[col]:
            vec, (lh, rh, rA) = smooth[i]
            lhs += list(all_vec & vec) + lh
            all_vec ^= vec
            rhs += [rh]
            rAs += [rA]
          i += 1

        factor = my_math.gcd(my_math.list_prod(rAs)*my_math.list_prod(lhs) - my_math.list_prod(rhs), n)
        if 1 < factor < n:
          break
      else:
        if verbose:
          print 'none found.'
        continue
      break

  if verbose:
    print 'factors found:'
    print factor, 'x', n/factor
    print 'time elapsed: %.3f seconds'%(time.time()-time1)
  return factor

if __name__ == "__main__":
  import argparse
  parser = argparse.ArgumentParser(description='Uses a MPQS to factor a composite number')
  parser.add_argument('composite', metavar='number_to_factor', type=long, help='the composite number to factor')
  parser.add_argument('--verbose', dest='verbose', action='store_true', help="enable verbose output")
  args = parser.parse_args()

  if args.verbose:
    mpqs(args.composite, args.verbose)
  else:
    time1 = time.time()
    print mpqs(args.composite)
    print 'time elapsed: %.3f seconds'%(time.time()-time1)

my_math.py

import math
import mpqs

# primes less than 212
small_primes = [
    2,  3,  5,  7, 11, 13, 17, 19, 23, 29, 31, 37,
   41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89,
   97,101,103,107,109,113,127,131,137,139,149,151,
  157,163,167,173,179,181,191,193,197,199,211]

# pre-calced sieve of eratosthenes for n = 2, 3, 5, 7
indices = [
    1, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47,
   53, 59, 61, 67, 71, 73, 79, 83, 89, 97,101,103,
  107,109,113,121,127,131,137,139,143,149,151,157,
  163,167,169,173,179,181,187,191,193,197,199,209]

# distances between sieve values
offsets = [
  10, 2, 4, 2, 4, 6, 2, 6, 4, 2, 4, 6,
   6, 2, 6, 4, 2, 6, 4, 6, 8, 4, 2, 4,
   2, 4, 8, 6, 4, 6, 2, 4, 6, 2, 6, 6,
   4, 2, 4, 6, 2, 6, 4, 2, 4, 2,10, 2]

# tabulated, mod 105
dindices =[
  0,10, 2, 0, 4, 0, 0, 0, 8, 0, 0, 2, 0, 4, 0,
  0, 6, 2, 0, 4, 0, 0, 4, 6, 0, 0, 6, 0, 0, 2,
  0, 6, 2, 0, 4, 0, 0, 4, 6, 0, 0, 2, 0, 4, 2,
  0, 6, 6, 0, 0, 0, 0, 6, 6, 0, 0, 0, 0, 4, 2,
  0, 6, 2, 0, 4, 0, 0, 4, 6, 0, 0, 2, 0, 6, 2,
  0, 6, 0, 0, 4, 0, 0, 4, 6, 0, 0, 2, 0, 4, 8,
  0, 0, 2, 0,10, 0, 0, 4, 0, 0, 0, 2, 0, 4, 2]

max_int = 2147483647

# returns the index of x in a sorted list a
# or the index of the next larger item if x is not present
# i.e. the proper insertion point for x in a
def binary_search(a, x):
  s = 0
  e = len(a)
  m = e >> 1
  while m != e:
    if a[m] < x:
      s = m
      m = (s + e + 1) >> 1
    else:
      e = m
      m = (s + e) >> 1
  return m

# divide and conquer list product
def list_prod(a):
  while len(a) > 1:
    a = [m*n for m,n in zip(a[::2], a[1::2] + [1])]
  return a[0]

# greatest common divisor of a and b
def gcd(a, b):
  while b:
    a, b = b, a%b
  return a

# extended gcd
def ext_gcd(a, m):
  a = int(a%m)
  x, u = 0, 1
  while a:
    x, u = u, x - (m/a)*u
    m, a = a, m%a
  return (m, x, u)

# legendre symbol (a|m)
# note: returns m-1 if a is a non-residue, instead of -1
def legendre(a, m):
  return pow(a, (m-1) >> 1, m)

# modular inverse of a mod m
def mod_inv(a, m):
  return ext_gcd(a, m)[1]

# modular sqrt(n) mod p
# p must be prime
def mod_sqrt(n, p):
  a = n%p
  if p%4 == 3:
    return pow(a, (p+1) >> 2, p)
  elif p%8 == 5:
    v = pow(a << 1, (p-5) >> 3, p)
    i = ((a*v*v << 1) % p) - 1
    return (a*v*i)%p
  elif p%8 == 1:
    # Shanks' method
    q = p-1
    e = 0
    while q&1 == 0:
      e += 1
      q >>= 1

    n = 2
    while legendre(n, p) != p-1:
      n += 1

    w = pow(a, q, p)
    x = pow(a, (q+1) >> 1, p)
    y = pow(n, q, p)
    r = e
    while True:
      if w == 1:
        return x

      v = w
      k = 0
      while v != 1 and k+1 < r:
        v = (v*v)%p
        k += 1

      if k == 0:
        return x

      d = pow(y, 1 << (r-k-1), p)
      x = (x*d)%p
      y = (d*d)%p
      w = (w*y)%p
      r = k
  else: # p == 2
    return a

#integer sqrt of n
def isqrt(n):
  c = (n << 2)/3
  d = c.bit_length()

  a = d>>1
  if d&1:
    x = 1 << a
    y = (x + (n >> a)) >> 1
  else:
    x = (3 << a) >> 2
    y = (x + (c >> a)) >> 1

  if x != y:
    x = y
    y = (x + n/x) >> 1
    while y < x:
      x = y
      y = (x + n/x) >> 1
  return x

# integer cbrt of n
def icbrt(n):
  d = n.bit_length()

  if d%3 == 2:
    x = 3 << d/3-1
  else:
    x = 1 << d/3

  y = (2*x + n/(x*x))/3
  if x != y:
    x = y
    y = (2*x + n/(x*x))/3
    while y < x:
      x = y
      y = (2*x + n/(x*x))/3
  return x

# strong probable prime
def is_sprp(n, b=2):
  if n < 2: return False
  d = n-1
  s = 0
  while d&1 == 0:
    s += 1
    d >>= 1

  x = pow(b, d, n)
  if x == 1 or x == n-1:
    return True

  for r in xrange(1, s):
    x = (x * x)%n
    if x == 1:
      return False
    elif x == n-1:
      return True

  return False

# lucas probable prime
# assumes D = 1 (mod 4), (D|n) = -1
def is_lucas_prp(n, D):
  P = 1
  Q = (1-D) >> 2

  # n+1 = 2**r*s where s is odd
  s = n+1
  r = 0
  while s&1 == 0:
    r += 1
    s >>= 1

  # calculate the bit reversal of (odd) s
  # e.g. 19 (10011) <=> 25 (11001)
  t = 0
  while s:
    if s&1:
      t += 1
      s -= 1
    else:
      t <<= 1
      s >>= 1

  # use the same bit reversal process to calculate the sth Lucas number
  # keep track of q = Q**n as we go
  U = 0
  V = 2
  q = 1
  # mod_inv(2, n)
  inv_2 = (n+1) >> 1
  while t:
    if t&1:
      # U, V of n+1
      U, V = ((U + V) * inv_2)%n, ((D*U + V) * inv_2)%n
      q = (q * Q)%n
      t -= 1
    else:
      # U, V of n*2
      U, V = (U * V)%n, (V * V - 2 * q)%n
      q = (q * q)%n
      t >>= 1

  # double s until we have the 2**r*sth Lucas number
  while r:
    U, V = (U * V)%n, (V * V - 2 * q)%n
    q = (q * q)%n
    r -= 1

  # primality check
  # if n is prime, n divides the n+1st Lucas number, given the assumptions
  return U == 0

## Baillie-PSW ##
# this is technically a probabalistic test, but there are no known pseudoprimes
def is_bpsw(n):
  if not is_sprp(n, 2): return False

  # idea shamelessly stolen from Mathmatica's PrimeQ
  # if n is a 2-sprp and a 3-sprp, n is necessarily square-free
  if not is_sprp(n, 3): return False

  a = 5
  s = 2
  # if n is a perfect square, this will never terminate
  while legendre(a, n) != n-1:
    s = -s
    a = s-a
  return is_lucas_prp(n, a)

# an 'almost certain' primality check
def is_prime(n):
  if n < 212:
    m = binary_search(small_primes, n)
    return n == small_primes[m]

  for p in small_primes:
    if n%p == 0:
      return False

  # if n is a 32-bit integer, perform full trial division
  if n <= max_int:
    p = 211
    while p*p < n:
      for o in offsets:
        p += o
        if n%p == 0:
          return False
    return True

  return is_bpsw(n)

# next prime strictly larger than n
def next_prime(n):
  if n < 2:
    return 2

  n += 1
  if n < 212:
    m = binary_search(small_primes, n)
    return small_primes[m]

  # find our position in the sieve rotation via binary search
  x = int(n%210)
  m = binary_search(indices, x)
  i = int(n + (indices[m] - x))

  # adjust offsets
  offs = offsets[m:] + offsets[:m]
  while True:
    for o in offs:
      if is_prime(i):
        return i
      i += o

# an infinite prime number generator
def primes(start = 0):
  for n in small_primes[start:]: yield n
  pg = primes(6)
  p = pg.next()
  q = p*p
  sieve = {221: 13, 253: 11}
  n = 211
  while True:
    for o in offsets:
      n += o
      stp = sieve.pop(n, 0)
      if stp:
        nxt = n/stp
        nxt += dindices[nxt%105]
        while nxt*stp in sieve:
          nxt += dindices[nxt%105]
        sieve[nxt*stp] = stp
      elif n < q:
        yield n
      else:
        sieve[q + dindices[p%105]*p] = p
        p = pg.next()
        q = p*p

# true if n is a prime power > 0
def is_prime_power(n):
  if n > 1:
    for p in small_primes:
      if n%p == 0:
        n /= p
        while n%p == 0: n /= p
        return n == 1

    r = isqrt(n)
    if r*r == n:
      return is_prime_power(r)

    s = icbrt(n)
    if s*s*s == n:
      return is_prime_power(s)

    p = 211
    while p*p < r:
      for o in offsets:
        p += o
        if n%p == 0:
          n /= p
          while n%p == 0: n /= p
          return n == 1

    if n <= max_int:
      while p*p < n:
        for o in offsets:
          p += o
          if n%p == 0:
            return False
      return True

    return is_bpsw(n)
  return False

max_trial = 1e10
max_pollard = 1e22

def factor(n):
  if n < max_trial:
    return factor_trial(n)
  for p in small_primes:
    if n%p == 0:
      return [p] + factor(n/p)
  if is_prime(n):
    return [n]
  if n < max_pollard:
    p = pollard_rho(n)
  else:
    p = lenstra_ecf(n) or mpqs.mpqs(n)
  return factor(p) + factor(n/p)

def factor_trial(n):
  a = []
  for p in small_primes:
    while n%p == 0:
      a += [p]
      n /= p
  i = 211
  while i*i < n:
    for o in offsets:
      i += o
      while n%i == 0:
        a += [i]
        n /= i
  if n > 1:
    a += [n]
  return a

def pollard_rho(n):
  # Brent's variant
  y, r, q = 0, 1, 1
  c, m = 9, 40
  g = 1
  while g == 1:
    x = y
    for i in range(r):
      y = (y*y + c) % n
    k = 0
    while k < r and g == 1:
      ys = y
      for j in range(min(m, r-k)):
        y = (y*y + c) % n
        q = q*abs(x-y) % n
      g = gcd(q, n)
      k += m
    r *= 2
  if g == n:
    ys = (ys*ys + c) % n
    g = gcd(n, abs(x-ys))
    while g == 1:
      ys = (ys*ys + c) % n
      g = gcd(n, abs(x-ys))
  return g

def ec_add((x1, z1), (x2, z2), (x0, z0), n):
  t1, t2 = (x1-z1)*(x2+z2), (x1+z1)*(x2-z2)
  x, z = t1+t2, t1-t2
  return (z0*x*x % n, x0*z*z % n)

def ec_double((x, z), (a, b), n):
  t1 = x+z; t1 *= t1
  t2 = x-z; t2 *= t2
  t3 = t1 - t2
  t4 = 4*b*t2
  return (t1*t4 % n, t3*(t4 + a*t3) % n)

def ec_multiply(k, p, C, n):
  # Montgomery ladder algorithm
  p0 = p
  q, p = p, ec_double(p, C, n)
  b = k >> 1
  while b > (b & -b):
    b ^= b & -b
  while b:
    if k&b:
      q, p = ec_add(p, q, p0, n), ec_double(p, C, n)
    else:
      q, p = ec_double(q, C, n), ec_add(p, q, p0, n),
    b >>= 1
  return q

def lenstra_ecf(n, m = 5):
  # Montgomery curves w/ Suyama parameterization.
  # Based on pseudocode found in:
  # "Implementing the Elliptic Curve Method of Factoring in Reconfigurable Hardware"
  # Gaj, Kris et. al
  # http://www.hyperelliptic.org/tanja/SHARCS/talks06/Gaj.pdf
  # Phase 2 is not implemented.
  B1, B2 = 8, 13
  for i in range(m):
    pg = primes()
    p = pg.next()
    k = 1
    while p < B1:
      k *= p**int(math.log(B1, p))
      p = pg.next()
    for s in range(B1, B2):
      u, v = s*s-5, 4*s
      x = u*u*u
      z = v*v*v
      t = pow(v-u, 3, n)
      P = (x, z)
      C = (t*(3*u+v) % n, 4*x*v % n)
      Q = ec_multiply(k, P, C, n)
      g = gcd(Q[1], n)
      if 1 < g < n: return g
    B1, B2 = B2, B1 + B2

if __name__ == "__main__":
  import argparse
  import time
  parser = argparse.ArgumentParser(description='Uses a various methods to factor a composite number')
  parser.add_argument('composite', metavar='number_to_factor', type=long, help='the composite number to factor')
  args = parser.parse_args()

  time1 = time.time()
  print factor(args.composite)
  print 'time elapsed: %.3f seconds'%(time.time()-time1)

Sample I/O:

$ pypy mpqs.py --verbose 523022617466601111760007224100074291200000001
smoothness bound: 14998
sieve size: 74990
log threshold: 18.7325510316
skipping primes less than: 37
sieving for smooths...
230 polynomials sieved (34495400 values)
found 849 smooths (242 from partials) in 0.786 seconds
solving for non-trivial congruencies...
factors found:
37280713718589679646221 x 14029308060317546154181
time elapsed: 1.011 seconds

Note: not using the --verbose option will give slightly better timings:

$ pypy mpqs.py 523022617466601111760007224100074291200000001
37280713718589679646221
time elapsed: 0.902 seconds

Basic Concepts

In general, a quadratic sieve is based on the following observation: any odd composite \$n\$ may be represented as:

\$n=(x+d)(x-d)=x^2-d^2\Rightarrow d^2=x^2-n\$

This is not very difficult to confirm. Since \$n\$ is odd, the distance between any two cofactors of \$n\$ must be even \$2d\$, where \$x\$ is the mid point between them. Moreover, the same relation holds for any multiple of \$n\$

\$abn=(ax+ad)(bx-bd)=abx^2-abd^2\Rightarrow abd^2=abx^2-abn\$

Note that if any such \$x\$ and \$d\$ can be found, it will immediately result in a (not necessarily prime) factor of \$n\$, since \$x+d\$ and \$x-d\$ both divide \$n\$ by definition. This relation can be further weakened - at the consequence of allowing potential trivial congruencies - to the following form:

\$d^2\equiv x^2\mod n\$

So in general, if we can find two perfect squares which are equivalent mod \$n\$, then it's fairly likely that we can directly produce a factor of \$n\$ a la \$\gcd(x±d,n)\$. Seems pretty simple, right?

Except it's not. If we intended to conduct an exhaustive search over all possible \$x\$, we would need to search the entire range from \$\sqrt{n}\$ to \$\sqrt{2n}\$, which is marginally smaller than full trial division, but also requires an expensive is_square operation each iteration to confirm the value of \$d\$. Unless it is known beforehand that \$n\$ has factors very near \$\sqrt{n}\$, trial division is likely to be faster.

Perhaps we can weaken this relation even more. Suppose we chose an \$x\$, such that for

\$y\equiv x^2\mod n\$

a full prime factorization of \$y\$ is readily known. If we had enough such relations, we should be able to construct an adequate \$d\$, if we choose a number of \$y\$ such that their product is a perfect square; that is, all prime factors are used an even number of times. In fact, if we have more such \$y\$ than the total number of unique prime factors they contain, a solution is guaranteed to exist; It becomes a system of linear equations. The question now becomes, how do we chose such \$x\$? That's where sieving comes into play.

The Sieve

Consider the polynomial:

\$y(x)=x^2-n\$

Then for any prime \$p\$ and integer \$k\$, the following is true:

\$y(x+kp)=(x+kp)^2-n\\y(x+kp)=x^2+2xkp+(kp)^2-n\\y(x+kp)=y(x)+2xkp+(kp)^2\equiv y(x)\mod p\$

This means that after solving for the roots of the polynomial mod \$p\$ - that is, you've found an \$x\$ such that \$y(x)\equiv 0\mod p\$, ergo \$y\$ is divisible by \$p\$ - then you have found an infinite number of such \$x\$. In this way, you can sieve over a range of \$x\$, identifying small prime factors of \$y\$, hopefully finding some for which all prime factors are small. Such numbers known as \$k-smooth\$, where \$k\$ is the largest prime factor used.

There's a few problems with this approach, though. Not all values of \$x\$ are adequate, in fact, there's only very few of them, centered around \$\sqrt{n}\$. Smaller values will become largely negative (due to the \$-n\$ term), and larger values will become too large, such that it is unlikely that their prime factorization consists only of small primes. There will be a number of such \$x\$, but unless the composite you're factoring is very small, it's highly unlikely that you'll find enough smooths to result in a factorization. And so, for larger \$n\$, it becomes necessary to sieve over multiple polynomials of a given form.

Multiple Polynomials

So we need more polynomials to sieve? How about this:

\$y(x)=(Ax+B)^2-n\$

That'll work. Note that \$A\$ and \$B\$ could literally be any integer value, and the math still holds. All we need to do is choose a few random values, solve for the root of the polynomial, and sieve the values close to zero. At this point we could just call it good enough: if you throw enough stones in random directions, you're bound to break a window sooner or later.

Except, there's a problem with that too. If the slope of the polynomial is large at the x-intercept, there'll only be a few suitable values to sieve per polynomial. It'll work, but you'll end up sieving a whole lot of polynomials before you get what you need. Can we do better?

We can do better. An observation, as a result of Montgomery is as follows: if \$A\$ and \$B\$ are chosen such that there exists some \$C\$ satisfying

\$B^2-n=AC\$

then the entire polynomial can be rewritten as

\$y(x)=(Ax+B)^2-n=(Ax)^2+2ABx+B^2-n=A(Ax^2+2Bx+C)\$

Furthermore, if \$A\$ is chosen to be a perfect square, the leading \$A\$ term can be neglected while sieving, resulting in much smaller values, and a much flatter curve. For such a solution to exist, \$n\$ must be a quadratic residue mod \$\sqrt{A}\$, which can be known immediately by computing the Legendre symbol:
\$(n|\sqrt{A})=1\$. Note that in order to solve for \$B\$, a complete prime factorization of \$\sqrt{A}\$ needs to be known (in order to take the modular square root \$\sqrt{n}\mod\sqrt{A}\$), which is why \$\sqrt{A}\$ is typically chosen to be prime.

It can then be shown that if \$A\approx\frac{\sqrt{2n}}{M}\$, then for all values of \$x\in[-M,M]\$:

\$|y(x)|\le\frac{M\sqrt{2n}}{2}\$

And now, finally, we have all the components necessary to implement our sieve. Or do we?

Powers of Primes as Factors

Our sieve, as described above, has one major flaw. It can identify which values of \$x\$ will result in a \$y\$ divisible by \$p\$, but it cannot identify whether or not this \$y\$ is divisible by a power of \$p\$. In order to determine that, we would need to perform trial division on the value to be sieved, until it is no longer divisible by \$p\$.

Edit: This is incorrect. The roots of \$n\$ mod \$p^k\$ can be computed directly from the roots mod \$p\$ through the use of Hensel Lifting. This current implementation does precisely this.

We seemed to have reached an impasse: the whole point of the sieve was so that we didn't have to do that. Time to check the playbook.

\$\ln(a\cdot b\cdot c\cdot d\cdot\ldots)=\ln(a)+\ln(b)+\ln(c)+\ln(d)+\ln(\ldots)\$

That looks pretty useful. If the sum of the \$\ln\$ of all of the small prime factors of \$y\$ is close to the expected value of \$\ln(y)\$, then it's almost a given that \$y\$ has no other factors. In addition, if we adjust the expected value down a little bit, we can also identify values as smooth which have several powers of primes as factors. In this way, we can use the sieve as a 'pre-screening' process, and only factor those values which are likely to be smooth.

This has a few other advantages as well. Note that small primes contribute very little to the \$\ln\$ sum, but yet they require the most sieve time. Sieving the value 3 requires more time than 11, 13, 17, 19, and 23 combined. Instead, we can just skip the first few primes, and adjust the threshold down accordingly, assuming a certain percentage of them would have passed.

Another result, is that a number of values will be allowed to 'slip through', which are mostly smooth, but contain a single large cofactor. We could just discard these values, but suppose we found another mostly smooth value, with exactly the same cofactor. We can then use these two values to construct a usable \$y\$; since their product will contain this large cofactor squared, it no longer needs to be considered.

Putting it all together

The last thing we need to do is to use these values of \$y\$ construct an adequate \$x\$ and \$d\$. Suppose we only consider the non-square factors of \$y\$, that is, the prime factors of an odd power. Then, each \$y\$ can be expressed in the following manner:

\$y_0=p_0^0\cdot p_1^1\cdot p_2^1\cdots p_n^0\\y_1=p_0^1\cdot p_1^0\cdot p_2^1\cdots p_n^1\\y_2=p_0^0\cdot p_1^0\cdot p_2^0\cdots p_n^1\\y_3=p_0^1\cdot p_1^1\cdot p_2^0\cdots p_n^0\\\vdots\$

which can be expressed in the matrix form:

\$M=\begin{bmatrix}0&1&1&\cdots&0\\1&0&1&\cdots&1\\0&0&0&\cdots&1\\1&1&0&\cdots&0\\\vdots\end{bmatrix}\$

The problem then becomes to find a vector \$v\$ such that \$vM=\vec{0}\mod 2\$, where \$\vec{0}\$ is the null vector. That is, to solve for the left null space of \$M\$. This can be done in a number of ways, the simplest of which is to perform Gaussian Elimination on \$M^T\$, replacing the row addition operation with a row xor. This will result in a number of null space basis vectors, any combination of which will produce a valid solution.

The construction of \$x\$ is fairly straight-forward. It is simply the product of \$Ax+B\$ for each of the \$y\$ used. The construction of \$d\$ is slightly more complicated. If we were to take the product of all \$y\$, we will end up with a value with 10s of thousands, if not 100s of thousands of digits, for which we need to find the square root. This calcuation is impractically expensive. Instead, we can keep track of the even powers of primes during the sieving process, and then use and and xor operations on the vectors of non-square factors to reconstruct the square root.

<s>I seem to have reached the 30000 character limit. Ahh well, I suppose that's good enough.</s> Saved a bunch of bytes by switching to \$MathJax\$.


问:
Pokemon GO Evolution Cost Calculator
Tired of always wondering how many more Pokémon you need to catch to get those high-tier evolutions? Wonder no more! Now you will write a complete program or function to calculate it for you!

The Challenge:

As input, your program will receive a list of the costs in candy to evolve a Pokémon to the next tier. (This list may be separated by any delimiter of your choosing, or as function arguments). Your program will then return or print the number of Pokémon that must be caught, including the one that will be evolved, to evolve through all the tiers given.

How do you calculate this? Like so:

1. Add up all the candy costs: 12 + 50 = 62

2. Subtract 3 candies from the total, this being from the one Pokémon you keep for evolving: 62 - 3 = 59

3. Divide this number by 4 (3 for catching, 1 for giving it to the Professor), always taking the ceil() of the result: ceil(59/4) = 15

4. Finally, add 1 to this total to get the total number of Pokémon you must catch, 16!

Example Input -> Output:

```<code>[4] -> 2
[50] -> 13
[12, 50] -> 16
[25, 100] -> 32
[19, 35, 5, 200] -> 65
</code>```

Winning:

The app has already taken up most of the space on your phone, so your program needs to be as short as possible. The complete program or function with the smallest byte count will be accepted in two weeks! (with any ties being settled by the earliest submitted entry!)

答:
Haskell, 17 bytes

```<code>(+1).(`div`4).sum
</code>```

ContentLicense: CC BY-SA 3.0

Title: Pokemon GO Evolution Cost Calculator

Body:
Haskell, 17 bytes

(+1).(`div`4).sum


问:
Find the minimum and maximum integers in an array, without using builtins

Challenge

Given an array of integers, received from stdin, function arguments, program arguments, or some other method:

Output only the minimum and maximum numbers in the array, through a return value, stdout, or other fitting methods.

Example session

```<code>> minmax( {0, 15, 2, 3, 7, 18, -2, 9, 6, -5, 3, 8, 9, -14} )
-14 18
</code>```

Reference implementation

```<code>// C++14

void minmax(std::vector<int> v) {
    int min = v[0]; int max = v[0];
    for(auto it : v) {
        if (*it < min)
            min = *it;
        if (*it > max)
            max = *it;
    }
    std::cout << min << ' ' << max << std::endl;
}
</code>```

Rules

You may not use a built-in function to calculate the values.

Standard loopholes disallowed.

Creative implementations encouraged.

This is code-golf, shortest answer wins but will not be selected.

Clarifications

If the array contains 1 element you need to output it twice.

If the minimum and maximum values are the same, you need to output them both.

答:
Python, 35 34 bytes

lambda s:sorted(s+s[:1])[::len(s)]
```

<p>Alternative version:</p>

<pre class="lang-python prettyprint-override"><code>lambda s:sorted(s+s)[::len(s)*2-1]
</code>```

Old version, 35 bytes.

lambda s:sorted(s+[s[0]])[::len(s)]
```

Fairly simple: take the input list, append the first element, sort it, then take the first and (length)th element of the resulting list. As the length of the input after appending an element is length + 1, this ends up taking the first and last element of said list, which are the minimum and maximum elements.

Body:
Python, 35 34 bytes

lambda s:sorted(s+s[:1])[::len(s)]

Alternative version:

lambda s:sorted(s+s)[::len(s)*2-1]

Old version, 35 bytes.

lambda s:sorted(s+[s[0]])[::len(s)]

Fairly simple: take the input list, append the first element, sort it, then take the first and (length)th element of the resulting list. As the length of the input after appending an element is length + 1, this ends up taking the first and last element of said list, which are the minimum and maximum elements.


问:
Golf a Brain-Flak Integer
Integers are tedious to represent in Brain-Flak. There are 8 operators:

```<code>()      Evaluates to 1, but does not push anything on any stack
[]      Evaluates to an indeterminate value for the purposes of this question
{}      Removes the top of the stack and evaluates to it
<>      Switches to or back from the alternate stack and evaluates to zero
(foo)   Pushes the value of the expression foo to the stack and evaluates to it
[foo]   Evaluates to the negation of foo
{foo}   Evaluates the expression foo until the top of the stack is zero
<foo>   Evaluates to zero but executes foo anyway
</code>```

foo may consist of multiple operators, in which case they are evaluated and summed. For example (()()) pushes 2 to the stack (and evaluates to 2 too).

Obviously, the (()...()) mechanism is not useful in Code Golf as large numbers would take n*2+2 bytes to represent. Your challenge is therefore to write a program or function that will output in as few bytes as possible a Brain-Flak program that will push a given positive integer n to the active stack. This program must not make any assumptions about the existing contents of the stacks, so it must not leave the stacks exchanged or add or remove extra values from the stacks.

Although your program or function must be capable of returning a working Brain-Flak program for all inputs from 1 to 1,000,000 the winner will be the program or function that generates the smallest set of appropriate Brain-Flak programs for all 1061 prime numbers between 1,000 and 10,000. You should note the total size of your outputs for those 1061 inputs as part of your submission. Your program or function may accept the integer and return the (string) Brain-Flak program in any of the usual acceptable I/O formats. Ties will be broken using the size of your program or function.

答:
Lua, 64664 bytes

Program prints the total length of the programs and the program for the 203rd prime (theres a line you can change to change which one is printed, or uncomment a line to print all the programs)

Right now the only optimization is x = 2*n + 1

Hopefully I will have time to add some more optimizations to lower the score.

local primeS = [[<INSERT PRIMES HERE>]]

local primes = {}

for num in primeS:gmatch("%d+") do
    table.insert(primes, num+0)
end

local progs = {}
progs[0] = ""
progs[1] = "()"
progs[2] = "()()"

local function half(n)
    if progs[n] then return progs[n] end
    local p = ""
    local div = math.floor(n/2)
    local rem = n%2 == 1 and "()" or ""
    return "("..progs[div].."){}"..rem
end

for i = 3, 10000 do

    local bin = half(i)

    progs[i] = progs[i-1] .. "()"

    if #bin < #progs[i] then
        progs[i] = bin
    end

    if i % 1000 == 0 then
        print(i)
    end

end

local n = 203 -- This is the program it outputs
print(n..", ("..progs[203]..")")

local len = 0
for i,v in ipairs(primes) do
    len = len + #progs[v] + 2
    --print(v.." ("..progs[v]..")\n")
end
print("Total len: "..len)
```

Title: Golf a Brain-Flak Integer


答:
Python, <s>59136</s> 58676 characters

Brainflak number golfing function:

m=11111
R=range(0,m)
R[1]="()"
R[2]="()()"
l=2
def a(v,r):
 if v>0 and v<m:
  if isinstance(R[v],int) or len(r)<len(R[v]):
   R[v]=r
   if v<R[0]:
    R[0]=v
def s(v,k):
 S=0
 while v>0:
  S+=v
  v-=k
 return S
p=lambda r:"("+r+")"
w=lambda r:"{({}["+r+"])}{}"
def q(r,v):
 for i in range(1,v):
  r="("+r+")"
 for i in range(1,v):
  r+="{}"
 return r
def e(r,v,k):
 for i in range(0,k):
  r=q(r,v)
 return r
while l<m:
 R[0]=l+1
 a(l*2,q(R[l],2))
 a(l*3,q(R[l],3))
 a(l*5,q(R[l],5))
 a(l*7,q(R[l],7))
 for i in range(1,l):
  a(l+i,R[l]+R[i])
  a(l-i,R[l]+"["+R[i]+"]")
  if l%i==0:
   t=s(l-i,i)
   a(s(l,i),p(R[l])+w(R[i]))
   a(l+2*t,p(R[l])+q(w(R[i]),2))
   a(l+4*t,p(R[l])+e(w(R[i]),2,2))
   a(l+8*t,p(R[l])+e(w(R[i]),2,3))
   a(l+16*t,p(R[l])+e(w(R[i]),2,4))
   a(l+32*t,p(R[l])+e(w(R[i]),2,5))
   a(l+64*t,p(R[l])+e(w(R[i]),2,6))
   a(l+128*t,p(R[l])+e(w(R[i]),2,7))
   a(l+3*t,p(R[l])+q(w(R[i]),3))
   a(l+9*t,p(R[l])+e(w(R[i]),3,2))
   a(l+27*t,p(R[l])+e(w(R[i]),3,3))
   a(l+5*t,p(R[l])+q(w(R[i]),5))
   a(l+6*t,p(R[l])+q(q(w(R[i]),3),2))
   a(l+10*t,p(R[l])+q(q(w(R[i]),5),2))
   a(l+15*t,p(R[l])+q(q(w(R[i]),5),3))
   a(l+12*t,p(R[l])+q(q(q(w(R[i]),3),2),2))
   a(l+18*t,p(R[l])+q(q(q(w(R[i]),3),3),2))
   a(l+20*t,p(R[l])+q(q(q(w(R[i]),5),2),2))
   a(l+24*t,p(R[l])+q(q(q(q(w(R[i]),3),2),2),2))
   a(l+36*t,p(R[l])+q(q(q(q(w(R[i]),3),3),2),2))
   a(l+40*t,p(R[l])+q(q(q(q(w(R[i]),5),2),2),2))
 l=R[0]
f=lambda v:p(R[v])
```

<p><strong>Prime number iteration:</strong></p>

<pre class="lang-python prettyprint-override"><code>def isPrime(v):
 i=2
 while i*i<=v:
  if v%i==0:
   return False
  i+=1
 return True

for i in range(1000,10000):
 if isPrime(i):
  print f(i)
</code>```

Output:

Pastebin

Explanation:

We pre-populate a list R of Brain-flak representation evaluating to individual integers over a larger than necessary range [1,m-1] to define our function f. The representations are formed by taking the lowest unused representation (indexed by l) and forming many new representations from it, keeping only the shortest. The lowest unused representation assumes that all number 1 to l have been assigned a representation, and that these representations have already been used to produce new numbers. If a value less than l gets a shorter representation we must go back and reproduce numbers beginning form that point. The function f produces a program saving the number to the stack by adding parenthesis.

I didn't know any Brainflak when I began this, and greatly appreciate Eamon Olive's answer for pointing out the formula for triangle numbers. Mostly I've generalized the summation and been relentless about checking sums and differences. Adding lots of multiples of sums has had a great effect.

For those who care, here is the scratch code I used to see which formulas were worth it.

Representation Formulas:

Multiplication by small primes:

(X){}

((X)){}{}

((((X)))){}{}{}{}

((((((X)))))){}{}{}{}{}{}

Addition X+Y:

XY

Subtraction X-Y:

X[Y]

Summation to and including X of increment Y:

(X){({}[Y])}{}

Multiples of summations to X of increment Y, plus X:

(X)({({}[Y])}{}){}

(X)(({({}[Y])}{})){}{}

(X)(({({}[Y])}{}){}){}

etc...

Title: Golf a Brain-Flak Integer

Body:
Integers are tedious to represent in Brain-Flak. There are 8 operators:

```<code>()      Evaluates to 1, but does not push anything on any stack
[]      Evaluates to an indeterminate value for the purposes of this question
{}      Removes the top of the stack and evaluates to it
<>      Switches to or back from the alternate stack and evaluates to zero
(foo)   Pushes the value of the expression foo to the stack and evaluates to it
[foo]   Evaluates to the negation of foo
{foo}   Evaluates the expression foo until the top of the stack is zero
<foo>   Evaluates to zero but executes foo anyway
</code>```

foo may consist of multiple operators, in which case they are evaluated and summed. For example (()()) pushes 2 to the stack (and evaluates to 2 too).

Obviously, the (()...()) mechanism is not useful in Code Golf as large numbers would take n*2+2 bytes to represent. Your challenge is therefore to write a program or function that will output in as few bytes as possible a Brain-Flak program that will push a given positive integer n to the active stack. This program must not make any assumptions about the existing contents of the stacks, so it must not leave the stacks exchanged or add or remove extra values from the stacks.

Although your program or function must be capable of returning a working Brain-Flak program for all inputs from 1 to 1,000,000 the winner will be the program or function that generates the smallest set of appropriate Brain-Flak programs for all 1061 prime numbers between 1,000 and 10,000. You should note the total size of your outputs for those 1061 inputs as part of your submission. Your program or function may accept the integer and return the (string) Brain-Flak program in any of the usual acceptable I/O formats. Ties will be broken using the size of your program or function.

Body:
Python, <s>59136</s> 58676 characters

Brainflak number golfing function:

m=11111
R=range(0,m)
R[1]="()"
R[2]="()()"
l=2
def a(v,r):
 if v>0 and v<m:
  if isinstance(R[v],int) or len(r)<len(R[v]):
   R[v]=r
   if v<R[0]:
    R[0]=v
def s(v,k):
 S=0
 while v>0:
  S+=v
  v-=k
 return S
p=lambda r:"("+r+")"
w=lambda r:"{({}["+r+"])}{}"
def q(r,v):
 for i in range(1,v):
  r="("+r+")"
 for i in range(1,v):
  r+="{}"
 return r
def e(r,v,k):
 for i in range(0,k):
  r=q(r,v)
 return r
while l<m:
 R[0]=l+1
 a(l*2,q(R[l],2))
 a(l*3,q(R[l],3))
 a(l*5,q(R[l],5))
 a(l*7,q(R[l],7))
 for i in range(1,l):
  a(l+i,R[l]+R[i])
  a(l-i,R[l]+"["+R[i]+"]")
  if l%i==0:
   t=s(l-i,i)
   a(s(l,i),p(R[l])+w(R[i]))
   a(l+2*t,p(R[l])+q(w(R[i]),2))
   a(l+4*t,p(R[l])+e(w(R[i]),2,2))
   a(l+8*t,p(R[l])+e(w(R[i]),2,3))
   a(l+16*t,p(R[l])+e(w(R[i]),2,4))
   a(l+32*t,p(R[l])+e(w(R[i]),2,5))
   a(l+64*t,p(R[l])+e(w(R[i]),2,6))
   a(l+128*t,p(R[l])+e(w(R[i]),2,7))
   a(l+3*t,p(R[l])+q(w(R[i]),3))
   a(l+9*t,p(R[l])+e(w(R[i]),3,2))
   a(l+27*t,p(R[l])+e(w(R[i]),3,3))
   a(l+5*t,p(R[l])+q(w(R[i]),5))
   a(l+6*t,p(R[l])+q(q(w(R[i]),3),2))
   a(l+10*t,p(R[l])+q(q(w(R[i]),5),2))
   a(l+15*t,p(R[l])+q(q(w(R[i]),5),3))
   a(l+12*t,p(R[l])+q(q(q(w(R[i]),3),2),2))
   a(l+18*t,p(R[l])+q(q(q(w(R[i]),3),3),2))
   a(l+20*t,p(R[l])+q(q(q(w(R[i]),5),2),2))
   a(l+24*t,p(R[l])+q(q(q(q(w(R[i]),3),2),2),2))
   a(l+36*t,p(R[l])+q(q(q(q(w(R[i]),3),3),2),2))
   a(l+40*t,p(R[l])+q(q(q(q(w(R[i]),5),2),2),2))
 l=R[0]
f=lambda v:p(R[v])

Prime number iteration:

def isPrime(v):
 i=2
 while i*i<=v:
  if v%i==0:
   return False
  i+=1
 return True

for i in range(1000,10000):
 if isPrime(i):
  print f(i)

Output:

Pastebin

Explanation:

We pre-populate a list R of Brain-flak representation evaluating to individual integers over a larger than necessary range [1,m-1] to define our function f. The representations are formed by taking the lowest unused representation (indexed by l) and forming many new representations from it, keeping only the shortest. The lowest unused representation assumes that all number 1 to l have been assigned a representation, and that these representations have already been used to produce new numbers. If a value less than l gets a shorter representation we must go back and reproduce numbers beginning form that point. The function f produces a program saving the number to the stack by adding parenthesis.

I didn't know any Brainflak when I began this, and greatly appreciate Eamon Olive's answer for pointing out the formula for triangle numbers. Mostly I've generalized the summation and been relentless about checking sums and differences. Adding lots of multiples of sums has had a great effect.

For those who care, here is the scratch code I used to see which formulas were worth it.

Representation Formulas:

Multiplication by small primes:

(X){}

((X)){}{}

((((X)))){}{}{}{}

((((((X)))))){}{}{}{}{}{}

Addition X+Y:

XY

Subtraction X-Y:

X[Y]

Summation to and including X of increment Y:

(X){({}[Y])}{}

Multiples of summations to X of increment Y, plus X:

(X)({({}[Y])}{}){}

(X)(({({}[Y])}{})){}{}

(X)(({({}[Y])}{}){}){}

etc...


答:
Actually, 4 bytes

```<code>Σ¼≈u
</code>```

Try it online!

Explanation:

```<code>Σ¼≈u
Σ     sum
 ¼    divide by 4
  ≈   floor
   u  add 1
</code>```

Title: Pokemon GO Evolution Cost Calculator

Body:
Tired of always wondering how many more Pokémon you need to catch to get those high-tier evolutions? Wonder no more! Now you will write a complete program or function to calculate it for you!

The Challenge:

As input, your program will receive a list of the costs in candy to evolve a Pokémon to the next tier. (This list may be separated by any delimiter of your choosing, or as function arguments). Your program will then return or print the number of Pokémon that must be caught, including the one that will be evolved, to evolve through all the tiers given.

How do you calculate this? Like so:

1. Add up all the candy costs: 12 + 50 = 62

2. Subtract 3 candies from the total, this being from the one Pokémon you keep for evolving: 62 - 3 = 59

3. Divide this number by 4 (3 for catching, 1 for giving it to the Professor), always taking the ceil() of the result: ceil(59/4) = 15

4. Finally, add 1 to this total to get the total number of Pokémon you must catch, 16!

Example Input -> Output:

```<code>[4] -> 2
[50] -> 13
[12, 50] -> 16
[25, 100] -> 32
[19, 35, 5, 200] -> 65
</code>```

Winning:

The app has already taken up most of the space on your phone, so your program needs to be as short as possible. The complete program or function with the smallest byte count will be accepted in two weeks! (with any ties being settled by the earliest submitted entry!)

Body:
Actually, 4 bytes

Σ¼≈u

Try it online!

Explanation:

Σ¼≈u
Σ     sum
 ¼    divide by 4
  ≈   floor
   u  add 1


问:
The Letter A without A
Your task is to display the letter "A" alone, without anything else, except any form of trailing newlines if you cannot avoid them, doing so in a program and/or snippet. Code that returns (instead of printing) is allowed.

Both the lowercase and uppercase versions of the letter "A" are acceptable (that is, unicode U+0061 or unicode U+0041. Other character encodings that aren't Unicode are allowed, but either way, the resulting output of your code must be the latin letter "A", and not any lookalikes or homoglyphs)

You must not use any of the below characters in your code, regardless of the character encoding that you pick:

"A", whether uppercase or lowercase.

"U", whether lowercase or uppercase.

X, whether uppercase or lowercase.

+

&

#

0

1

4

5

6

7

9

Cheating, loopholes, etc, are not allowed.

Since this is code-golf, the shortest solution, in bytes, that follows all the rules, is the winner.

Validity Checker

This Stack Snippet checks to make sure your code doesn't use the restricted characters. It might not work properly for some character encodings.

var t = prompt("Input your code.");

if (/[AaUuXx+&#0145679]/.test(t)) {
  alert("Contains a disallowed character!");
} else {
  alert("No disallowed characters");
}```
</div>
</div>
</p>

<p><em>This Stack Snippet that makes sure you don't have a disallowed character is also available on JSFiddle.</em></p>

<h3>Leaderboard</h3>

<p><div class="snippet" data-lang="js" data-hide="true" data-console="false" data-babel="false">
<div class="snippet-code snippet-currently-hidden">
<pre class="snippet-code-js lang-js prettyprint-override"><code>var QUESTION_ID=90349,OVERRIDE_USER=58717;function answersUrl(e){return"https://api.stackexchange.com/2.2/questions/"+QUESTION_ID+"/answers?page="+e+"&pagesize=100&order=desc&sort=creation&site=codegolf&filter="+ANSWER_FILTER}function commentUrl(e,s){return"https://api.stackexchange.com/2.2/answers/"+s.join(";")+"/comments?page="+e+"&pagesize=100&order=desc&sort=creation&site=codegolf&filter="+COMMENT_FILTER}function getAnswers(){jQuery.ajax({url:answersUrl(answer_page++),method:"get",dataType:"jsonp",crossDomain:!0,success:function(e){answers.push.apply(answers,e.items),answers_hash=[],answer_ids=[],e.items.forEach(function(e){e.comments=[];var s=+e.share_link.match(/\d+/);answer_ids.push(s),answers_hash[s]=e}),e.has_more||(more_answers=!1),comment_page=1,getComments()}})}function getComments(){jQuery.ajax({url:commentUrl(comment_page++,answer_ids),method:"get",dataType:"jsonp",crossDomain:!0,success:function(e){e.items.forEach(function(e){e.owner.user_id===OVERRIDE_USER&&answers_hash[e.post_id].comments.push(e)}),e.has_more?getComments():more_answers?getAnswers():process()}})}function getAuthorName(e){return e.owner.display_name}function process(){var e=[];answers.forEach(function(s){var r=s.body;s.comments.forEach(function(e){OVERRIDE_REG.test(e.body)&&(r="<h1>"+e.body.replace(OVERRIDE_REG,"")+"</h1>")});var a=r.match(SCORE_REG);a&&e.push({user:getAuthorName(s),size:+a[2],language:a[1],link:s.share_link})}),e.sort(function(e,s){var r=e.size,a=s.size;return r-a});var s={},r=1,a=null,n=1;e.forEach(function(e){e.size!=a&&(n=r),a=e.size,++r;var t=jQuery("#answer-template").html();t=t.replace("{{PLACE}}",n+".").replace("{{NAME}}",e.user).replace("{{LANGUAGE}}",e.language).replace("{{SIZE}}",e.size).replace("{{LINK}}",e.link),t=jQuery(t),jQuery("#answers").append(t);var o=e.language;/<a/.test(o)&&(o=jQuery(o).text()),s[o]=s[o]||{lang:e.language,user:e.user,size:e.size,link:e.link}});var t=[];for(var o in s)s.hasOwnProperty(o)&&t.push(s[o]);t.sort(function(e,s){return e.lang>s.lang?1:e.lang<s.lang?-1:0});for(var c=0;c<t.length;++c){var i=jQuery("#language-template").html(),o=t[c];i=i.replace("{{LANGUAGE}}",o.lang).replace("{{NAME}}",o.user).replace("{{SIZE}}",o.size).replace("{{LINK}}",o.link),i=jQuery(i),jQuery("#languages").append(i)}}var ANSWER_FILTER="!t)IWYnsLAZle2tQ3KqrVveCRJfxcRLe",COMMENT_FILTER="!)Q2B_A2kjfAiU78X(md6BoYk",answers=[],answers_hash,answer_ids,answer_page=1,more_answers=!0,comment_page;getAnswers();var SCORE_REG=/<h\d>\s*([^\n,]*[^\s,]),.*?(\d+)(?=[^\n\d<>]*(?:<(?:s>[^\n<>]*<\/s>|[^\n<>]+>)[^\n\d<>]*)*<\/h\d>)/,OVERRIDE_REG=/^Override\s*header:\s*/i;</code>```

body{text-align:left!important}#answer-list,#language-list{padding:10px;width:290px;float:left}table thead{font-weight:700}table td{padding:5px}```
<pre class="snippet-code-html lang-html prettyprint-override"><code><script src="https://ajax.googleapis.com/ajax/libs/jquery/2.1.1/jquery.min.js"></script> <link rel="stylesheet" type="text/css" href="//cdn.sstatic.net/codegolf/all.css?v=83c949450c8b"> <div id="answer-list"> <h2>Leaderboard</h2> <table class="answer-list"> <thead> <tr><td></td><td>Author</td><td>Language</td><td>Size</td></tr></thead> <tbody id="answers"> </tbody> </table> </div><div id="language-list"> <h2>Winners by Language</h2> <table class="language-list"> <thead> <tr><td>Language</td><td>User</td><td>Score</td></tr></thead> <tbody id="languages"> </tbody> </table> </div><table style="display: none"> <tbody id="answer-template"> <tr><td>{{PLACE}}</td><td>{{NAME}}</td><td>{{LANGUAGE}}</td><td>{{SIZE}}</td><td>Link</td></tr></tbody> </table> <table style="display: none"> <tbody id="language-template"> <tr><td>{{LANGUAGE}}</td><td>{{NAME}}</td><td>{{SIZE}}</td><td>Link</td></tr></tbody> </table></code>```

答:
Lua, 36 bytes

This one took me a while since MANY of the standard Lua functions are taken away from the rules (all of math, string.char, string.match, pairs, even next)

This takes advantage of the fact that Lua has a global _VERSION that usually starts with "Lua" (e.g., Lua 5.1, or similar for other versions), so unless this is run on a non-mainstream interpreter, the third character will be an 'a'

```<code>print(({_VERSION:find("..(.)")})[3])
</code>```

The {...}[3] is to group the results of find which also includes the indices where it matched, and then return the third item which is the matched character (the a)

Title: The Letter A without A

Body:
Lua, 36 bytes

This one took me a while since MANY of the standard Lua functions are taken away from the rules (all of math, string.char, string.match, pairs, even next)

This takes advantage of the fact that Lua has a global _VERSION that usually starts with "Lua" (e.g., Lua 5.1, or similar for other versions), so unless this is run on a non-mainstream interpreter, the third character will be an 'a'

print(({_VERSION:find("..(.)")})[3])

The {...}[3] is to group the results of find which also includes the indices where it matched, and then return the third item which is the matched character (the a)


问:
The Permutation Pigeon-hole Principle
In the game of sudoku, many players like to "pencil in" possible numbers that can go in each square:

The above row can be represented as an array:

```<code>[[1,2,9], [6], [5], [7], [1,2,9], [1,2,9], [3], [1,2,4], [8]]
</code>```

Now, notice that there is only 1 place where a 4 can go.  This effectively lets us simplify the above list to:

```<code>[[1,2,9], [6], [5], [7], [1,2,9], [1,2,9], [3], [4], [8]]
</code>```

The goal of this challenge is to take a list of possible numbers in a permutation, and deduce which possibilities can be eliminated.

As another example, lets say you have the following array of possibilities:

```<code>[[0,1,3], [0,2,3], [1,2], [1,2]]
</code>```

The last two places must be filled with 1 and 2.  Therefore, we can remove those possibilities from the first two elements in the array:

```<code>[[0,3], [0,3], [1,2], [1,2]]
</code>```

As another example:

```<code>[[0,1,2,3], [0,2], [0,2], [0,2]]
</code>```

Its impossible to construct a permutation from the above possibilities, as there's only 1 location for both 1 and 3, and you would want to return an empty array.

You need to input a list of possibilities and output the remaining possibilities after the maximum number of possibilities have been eliminated.

If a particular array is impossible, you either need to return an empty array, or an array where one of the subarrays is empty.

You may assume that the array will be well-formed, and have at least 1 element.

Given an array of size N, you can assume the numbers in the subarray will always be in the range [0:N), and that N <= 10

You may not assume that every number from 0 to N-1 will be present

You may assume that numbers within a single subarray are unique.

If a subarray contains only a single possibility, you can either represent the possibility in an array or by itself.  [[1],[2],[0]], [1,2,0], [[1,2],0,[1,2]] are all valid.

You may accept the array either in a reasonable string format or in list/array format.

Subarrays can be in any order.

Instead of dealing with ragged arrays, you can pad empty places with -1.

Test cases

```<code>[[0]]                                         -> [[0]]
[[1],[0]]                                     -> [[1],[0]]
[[1],[1]]                                     -> []
[[1],[0,1]]                                   -> [[1],[0]]
[[0,1,2],[1,2],[1,2]]                         -> [[0],[1,2],[1,2]]
[[0,1],[1,2],[0,2]]                           -> [[0,1],[1,2],[0,2]]
[[2,1],[1,2],[1,2]]                           -> []
[[0,3],[2,1],[3,0],[3,2]]                     -> [[0,3],[1],[0,3],[2]]
[[0,1],[0,1],[2,3],[2,3,0]]                   -> [[0,1],[0,1],[2,3],[2,3]]
[[0,1],[0,3],[3,2],[0]]                       -> [[1],[3],[2],[0]]
[[3,5,2],[0,2,4],[4,0],[0,1,3,5],[2,1],[2,4]] -> [[3,5],[0,2,4],[4,0],[3,5],[1],[2,4]]
[[6,9,8,4],[4,5],[5,3,6],[3,8,6,1,4],[3,1,9,6],[3,7,0,2,4,5],[9,5,6,8],[6,5,8,1,3,7],[8],[8,0,6,2,5,6,3]] -> [[6,9,4],[4,5],[5,3,6],[3,6,1,4],[3,1,9,6],[0,2],[9,5,6],[7],[8],[0,2]]
[[3,5,0],[5,7],[5,1,2],[1,3,0],[5,3],[5,0],[5,3,7,8,0,6],[7,5,0,1,8],[1,0,8],[0,6]] -> []
[[9,0,2,3,7],[0,7,6,5],[6,9,4,7],[9,1,2,3,0,5],[2,8,5,7,4,6],[6,5,7,1],[5,9,4],[5,9,3,8,1],[5,0,6,4],[0,7,2,1,3,4,8]] -> [[9,0,2,3,7],[0,7,6,5],[6,9,4,7],[9,1,2,3,0,5],[2,8,5,7,4,6],[6,5,7,1],[5,9,4],[5,9,3,8,1],[5,0,6,4],[0,7,2,1,3,4,8]]
[[2,6,0],[0,4,3],[0,6,2],[0,7],[0,9,2,3,6,1,4],[1,7,2],[2,7,8],[8,6,7],[6,5,2,8,0],[5,8,1,4]] -> [[2,6,0],[3],[0,6,2],[0,7],[9],[1],[2,7,8],[8,6,7],[5],[4]]
[[8],[8,0,6,5,7,2,4,1],[8,6,9,3,5,0,7],[3,9,1,0],[9],[9,2,6],[2,8,3],[3,1,6,8,2],[6],[6,4,5,3,0,7]] -> [[8],[5,7,4],[5,7],[0],[9],[2],[3],[1],[6],[4,5,7]]
[[8,1,0],[5,8,7,6,2,0],[6,8,2],[2,4,0,9],[4,1,7,3,6,8],[8,1],[8,0,3],[0,8,2],[0,8,3],[1,8,0]] -> []
</code>```

This is a code-golf so make your answers as short as possible!

答:
Mathematica, 46 bytes

```<code>Union/@Thread@Select[Tuples@#,DuplicateFreeQ]&
</code>```

Title: The Permutation Pigeon-hole Principle

Body:
Mathematica, 46 bytes

Union/@Thread@Select[Tuples@#,DuplicateFreeQ]&


问:
Thwart Lepton compression
Dropbox recently released Lepton (GitHub), a method that losslessly compresses JPEG images round-trip, saving an average of 22%.

Because of the pigeonhole principle, any general compression algorithm cannot be guaranteed to result in a smaller file (general because it doesn't apply to inputs constrained to a specific format). Lepton exploits common characteristics about JPEGs, which if subverted, might pigeonhole it to produce a file larger than the source.

Requirements

Write a program that generates:

A valid JPEG/JFIF image,

with a size between 0.5 MB and 1 MB,

no smaller than 256 × 256 px,

no larger than 4096 × 4096 px,

recognizable by Lepton (it can successfully "compress" to a .lep image), and

decompresses to an identical .jpg (as the input).

APPx, COM, and other metadata, non-graphical marker sections are restricted in the JPEG (injecting arbitrary amounts of random bytes into the image to asymptotically approach 1:1 compression is lame.)

an APP0 JFIF marker is permitted but no thumbnail is allowed (should be exactly 16 bytes)

tl;dr If you're not intentionally shoving metadata into an EXIF segment and you disable any sort of thumbnail your language library of choice wants to put into the image, that should be OK.

Post the code and image.

If you want to write a program that produces a Lepton image that when converted yields a JPEG meeting the criteria, that's fine. It must remain identical across arbitrarily many JPEG → Lepton → JPEG → ... cycles.

Scoring

The byte size of the Lepton image divided by the source JPEG image. Higher (worse Lepton compression) is better. Run Lepton with default flags and switches.

Getting Lepton

A 5-second crash-course to build Lepton:

```<code>git clone https://github.com/dropbox/lepton.git
cd lepton
./autogen.sh && ./configure && make

# fish shell: ./autogen.sh ;and ./configure ;and make
</code>```

Then ./lepton --help should tell you things.

答:
Python 3 + mozjpeg + /dev/urandom, 720×720: avg. score 102%

Depends on the mozjpeg package, the code assumes it's installed into /usr/local/opt/mozjpeg. (on OS X it's trivial to install, just run brew install mozjpeg)

Also it depends on /dev/urandom special file, it's used for generating random data.

The code simply feeds random data to mozjpeg compressor (in TGA format, because cjpeg understands it and it has very simple header), and lets it create an optimized jpeg file. Quality is set to the maximum because it makes DCT coefficients the least compressable, and it doesn't matter much what algorithm is used to compress uncompressable data.

I checked that the jpeg->lepton->jpeg cycle is lossless — it is true.

```<code>import subprocess
from subprocess import PIPE

c_mozjpeg_path = '/usr/local/opt/mozjpeg/bin/cjpeg'
cjpeg_params = '-quality 100 -dc-scan-opt 2 -dct float -targa'
image_size = 720

def write_random_tga_image(w, h, of, rf):
    def wb(value, size):
        of.write(int.to_bytes(value, size, 'little'))

    wb(0, 2)
    wb(3, 1)
    wb(0, 9)
    wb(w, 2)
    wb(h, 2)
    wb(8, 1)
    wb(0, 1)

    data_size = w * h
    while data_size > 0:
        data_size -= of.write(rf.read(data_size))

def main():
    with open('/dev/urandom', 'rb') as rf:
        with open('oops.jpg', 'wb') as f:
            p = subprocess.Popen((c_mozjpeg_path,) + tuple(cjpeg_params.split(' ')), stdin=PIPE, stdout=f)
            write_random_tga_image(image_size, image_size, p.stdin, rf)
            p.communicate()

if __name__ == '__main__':
    main()
</code>```

Code is not golfed, obviously.

Example image:

Fun fact: generated JPEG file is larger than the source uncompressed TGA image, even though JPEG uses lossy compression.

Fun fact 2: Imgur (the default image hosting for SO) does very bad job at this file — for some reason it recompresses it to lower quality, even though it's less than 1MB. So I used Github to upload the example image.

Fun fact 3: In general, mozjpeg indeed does better JPEG compression while staying compatible with existing JPEG decoders. And it also has tool to losslessly optimize JPEG files, too — jpegtran.

Title: Thwart Lepton compression

Body:
Python 3 + mozjpeg + /dev/urandom, 720×720: avg. score 102%

Depends on the mozjpeg package, the code assumes it's installed into /usr/local/opt/mozjpeg. (on OS X it's trivial to install, just run brew install mozjpeg)

Also it depends on /dev/urandom special file, it's used for generating random data.

The code simply feeds random data to mozjpeg compressor (in TGA format, because cjpeg understands it and it has very simple header), and lets it create an optimized jpeg file. Quality is set to the maximum because it makes DCT coefficients the least compressable, and it doesn't matter much what algorithm is used to compress uncompressable data.

I checked that the jpeg->lepton->jpeg cycle is lossless — it is true.

import subprocess
from subprocess import PIPE

c_mozjpeg_path = '/usr/local/opt/mozjpeg/bin/cjpeg'
cjpeg_params = '-quality 100 -dc-scan-opt 2 -dct float -targa'
image_size = 720

def write_random_tga_image(w, h, of, rf):
    def wb(value, size):
        of.write(int.to_bytes(value, size, 'little'))

    wb(0, 2)
    wb(3, 1)
    wb(0, 9)
    wb(w, 2)
    wb(h, 2)
    wb(8, 1)
    wb(0, 1)

    data_size = w * h
    while data_size > 0:
        data_size -= of.write(rf.read(data_size))

def main():
    with open('/dev/urandom', 'rb') as rf:
        with open('oops.jpg', 'wb') as f:
            p = subprocess.Popen((c_mozjpeg_path,) + tuple(cjpeg_params.split(' ')), stdin=PIPE, stdout=f)
            write_random_tga_image(image_size, image_size, p.stdin, rf)
            p.communicate()

if __name__ == '__main__':
    main()

Code is not golfed, obviously.

Example image:

Fun fact: generated JPEG file is larger than the source uncompressed TGA image, even though JPEG uses lossy compression.

Fun fact 2: Imgur (the default image hosting for SO) does very bad job at this file — for some reason it recompresses it to lower quality, even though it's less than 1MB. So I used Github to upload the example image.

Fun fact 3: In general, mozjpeg indeed does better JPEG compression while staying compatible with existing JPEG decoders. And it also has tool to losslessly optimize JPEG files, too — jpegtran.


问:
Text animation format Code Golf
Working with old files from DOS, I found an ancient text animation format, that specify how much time each frame must be shown in the screen. It's very simple:

F1W5F2W3

This example above indicates that Frame 1 must be shown for 5 units of time, and frame 2 for 3 units more. Now I want to convert it to:

1, 1, 1, 1, 1, 2, 2, 2

In the converted example below, frame 1 is repeated 5 times, and frame 2 three times.

Rules: The only allowed letters are F and W (and the lowercase counterparts, f and w). Numbers can go up to 255. Total input String with no more than 255 chars.

Other examples:

Input: F10W3F11W3F10W3F12W3

Output: 10, 10, 10, 11, 11, 11, 10, 10, 10, 12, 12, 12

Input: f1w0f2w10

Output: 2, 2, 2, 2, 2, 2, 2, 2, 2, 2

Short code wins.

Edited - Accepted formats:
Any answer containing the output and ignoring spaces. Ex: 1, 1, 2, 2 or [1,1,2,2] or (1,1,2,2). But NOT 1 1 2 2 or [1; 1; 2; 2] or [1 1 2 2].

答:
Pyth, 13 bytes

```<code>smMcr7Xrz0Gd2
</code>```

Test suite

This program uses one of my favorite tricks, mM, which maps the map function over a list. It is also uses a relic of Pyth's distant past, the r7 function, which splits a string on whitespace and evaluates the contents, which just so happens to be absolutely perfect for this challenge. Here's how it works, and each intermediate state:

```<code>smMcr7Xrz0Gd2     f1W5F2w3

        z        "f1W5F2w3"       Implicit: z is the input.
       r 0       "f1w5f2w3"       Lowercase
      X   Gd     " 1 5 2 3"       Replace all lowercase letters (G) with spaces (d)
    r7           [1, 5, 2, 3]     Split on whitespace, eval each string.
   c        2    [[1, 5], [2, 3]] Split into groups of length 2
 mM                               Apply the map function to each pair
                 [[1, 1, 1, 1, 1], [2, 2, 2]]
s                                 Flatten
                 [1, 1, 1, 1, 1, 2, 2, 2]
</code>```

Title: Text animation format Code Golf

Body:
Pyth, 13 bytes

smMcr7Xrz0Gd2

Test suite

This program uses one of my favorite tricks, mM, which maps the map function over a list. It is also uses a relic of Pyth's distant past, the r7 function, which splits a string on whitespace and evaluates the contents, which just so happens to be absolutely perfect for this challenge. Here's how it works, and each intermediate state:

smMcr7Xrz0Gd2     f1W5F2w3

        z        "f1W5F2w3"       Implicit: z is the input.
       r 0       "f1w5f2w3"       Lowercase
      X   Gd     " 1 5 2 3"       Replace all lowercase letters (G) with spaces (d)
    r7           [1, 5, 2, 3]     Split on whitespace, eval each string.
   c        2    [[1, 5], [2, 3]] Split into groups of length 2
 mM                               Apply the map function to each pair
                 [[1, 1, 1, 1, 1], [2, 2, 2]]
s                                 Flatten
                 [1, 1, 1, 1, 1, 2, 2, 2]


问:
Russian Roulette
Write a program which plays Russian Roulette!

If the program is started,

there should be a 5 in 6 chance of it ending normally after printing "I survived!"

there should be a 1 in 6 chance of the program crashing. (segmentation fault, etc.)

No input, and no other outputs are allowed.

The randomness must be fair: it must have a uniform probability distribution. This means an uninitialized variable (or a RNG without seed) MOD 6 will not be sufficient.

If the solution works with only one dedicated operating system / platform, you will receive a 6 byte penalty to the score.

Shortest code wins, not sooner than 10 days after first valid answer.

答:
vba, 27

```<code>?1/int(6*rnd),"I Survived!"
</code>```

used in immediate window.

On failure, an error window stating:

appears

Title: Russian Roulette

Body:
vba, 27

?1/int(6*rnd),"I Survived!"

used in immediate window.

On failure, an error window stating:

appears