title: Fastest Home Prime Generator

text:
# Question
Title: Fastest Home Prime Generator
What is a home prime?

For an example, take HP(4). First, find the prime factors. The prime factors of 4 (*in numerical order from least to greatest, always*) are 2, 2. Take those factors as a literal number. 2, 2 becomes 22. This process of factoring continues until you reach a prime number.

```
number    prime factors
4         2, 2
22        2, 11
211       211 is prime

```

Once you reach a prime number, the sequence ends. HP(4)=211. Here's a longer example, with 14:

```
number    prime factors
14        2, 7
27        3, 3, 3
333       3, 3, 37
3337      47, 71
4771      13, 367
13367     13367 is prime

```

Your challenge is to create a program that will calculate HP(x) given x, and do it *as quickly as possible*. You may use whatever resources you want, other than a list of known home primes.

Take note, these numbers become very large very fast. At x=8, HP(x) jumps all the way to 3331113965338635107. HP(49) has yet to be found.

Program speed will be tested on a Raspberry Pi 2, averaging the following inputs:

```
16
20
64
65
80

```

If you have a Raspberry Pi 2, time the program yourself, then average the times.

# Answer
> 6 votes
# Mathematica, HP(80) in ~0.88s

```
NestWhile[
  FromDigits[
    Flatten[IntegerDigits /@
      ConstantArray @@@ FactorInteger[#]]] &, #, CompositeQ] &

```

Anonymous function. Takes a number as input and returns a number as output.

# Answer
> 6 votes
# PyPy 5.4.1 64bit (linux), HP(80) ~ 1.54s

*The 32bit version will time slightly slower.*

I use four different factorization methods, with empirically determined breakpoints:

I tried for a while to find a clean break between ECF and MPQS, but there doesn't seem to be one. However, if *n* contains a small factor, ECF will usually find it almost immediately, so I've opted to test only a few curves, before moving on to MPQS.

Currently, it's only ~2x slower than Mathmatica, which I certainly consider a success.

---

**home-prime.py**

```
import math
import my_math
import mpqs

max_trial = 1e10
max_pollard = 1e22

def factor(n):
  if n < max_trial:
    return factor_trial(n)
  for p in my_math.small_primes:
    if n%p == 0:
      return [p] + factor(n/p)
  if my_math.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 my_math.small_primes:
    while n%p == 0:
      a += [p]
      n /= p
  i = 211
  while i*i < n:
    for o in my_math.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 = my_math.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 = my_math.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 = my_math.gcd(Q[1], n)
      if 1 < g < n: return g
    B1, B2 = B2, B1 + B2

if __name__ == '__main__':
  import time
  import sys
  for n in sys.argv[1:]:
    t0 = time.time()
    i = int(n)
    f = []
    while len(f) != 1:
      f = sorted(factor(i))
      #print i, f
      i = int(''.join(map(str, f)))
    t1 = time.time()-t0
    print n, i
    print '%.3fs'%(t1)
    print

```

**Sample Timings**

```
    $ pypy home-prime.py 8 16 20 64 65 80
8 3331113965338635107
0.005s

16 31636373
0.001s

20 3318308475676071413
0.004s

64 1272505013723
0.000s

65 1381321118321175157763339900357651
0.397s

80 313169138727147145210044974146858220729781791489
1.537s

```

The average of the 5 is about 0.39s.

---

**Dependencies**

`mpqs.py` is taken directly from my answer to Fastest semiprime factorization with a few very minor modifications.

**mpqs.py**

```
import math
import my_math
import time

# Multiple Polynomial Quadratic Sieve
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(5 * math.log(n, 10)**2)

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

  # find a number of small primes for which n is a quadratic residue
  p = 2
  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]
      mod_root += [int(my_math.mod_sqrt(n, p))]
      log_p += [math.log(p, 10)]
      num_prime += 1
    elif leg == 0:
      if verbose:
        print 'trial division found factors:'
        print p, 'x', n/p
      return p

    p = my_math.next_prime(p)

  # size of the sieve
  x_max = bound*8

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

  # fudging the threshold down a bit makes it easier to find powers of primes as factors
  # as well as partial-partial relationships, but it also makes the smoothness check slower.
  # there's a happy medium somewhere, depending on how efficient the smoothness check is
  thresh = math.log(m_val, 10) * 0.735

  # 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, assuming ~1/4 of them pass
  min_prime = int(thresh*3)
  fudge = sum(log_p[i] for i,p in enumerate(prime) if p < min_prime)/4
  thresh -= fudge

  sieve_primes = [p for p in prime if p >= min_prime]
  sp_idx = prime.index(sieve_primes[0])

  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
    sums = [0.0]*(2*x_max)
    i = sp_idx
    for p in sieve_primes:
      logp = log_p[i]

      inv_A = my_math.mod_inv(A, p)
      # modular root of the quadratic
      a = int(((mod_root[i] - B) * inv_A)%p)
      b = int(((p - mod_root[i] - B) * inv_A)%p)

      amx = a+x_max
      bmx = b+x_max

      ax = amx-p
      bx = bmx-p

      k = p
      while k < x_max:
        sums[k+ax] += logp
        sums[k+bx] += logp
        sums[amx-k] += logp
        sums[bmx-k] += logp
        k += p

      if k+ax < x_max:
        sums[k+ax] += logp
      if k+bx < x_max:
        sums[k+bx] += logp
      if amx-k > 0:
        sums[amx-k] += logp
      if bmx-k > 0:
        sums[bmx-k] += logp
      i += 1

    # check for smooths
    x = -x_max
    for v in sums:
      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

    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 %f seconds'%(num_smooth, num_partial, time.time()-time1)
        print 'solving for non-trivial congruencies...'

      used_prime_list = sorted(list(used_prime))

      # 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_list:
          if p in vec: bit_fields[i] |= mask
          i += 1
        mask <<= 1

      # 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: %f 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: %f seconds'%(time.time()-time1)

```

`my_math.py` is taken from the same post as `mpqs.py`, however I've also added in the unbounded prime number generator I used in my answer to Find the largest gap between good primes.

**my\_math.py**

```
# 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):
  size = len(a)
  if size == 1:
    return a[0]
  return list_prod(a[:size>>1]) * list_prod(a[size>>1:])

# 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:
    # Shank's 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*4/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

  # first odd larger than n
  n = (n + 1) | 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

```

# Answer
> 2 votes
# C#

```
using System;
using System.Linq;

public class Program
{
    public static void Main(string[] args) {

        Console.Write("Enter Number: ");

        Int64 n = Convert.ToInt64(Console.ReadLine());

        Console.WriteLine("Start Time: " + DateTime.Now.ToString("HH:mm:ss.ffffff"));
        Console.WriteLine("Number, Factors");

        HomePrime(n);

        Console.WriteLine("End Time: " + DateTime.Now.ToString("HH:mm:ss.ffffff"));
        Console.ReadLine();
    }

    public static void HomePrime(Int64 num) {
        string s = FindFactors(num);
        if (CheckPrime(num,s) == true) {
            Console.WriteLine("{0} is prime", num);
        } else {
            Console.WriteLine("{0}, {1}", num, s);
            HomePrime(Convert.ToInt64(RemSp(s)));
        }
    }

    public static string FindFactors(Int64 num) {
        Int64 n, r, t = num;
        string f = "";
        for (n = num; n >= 2; n--) {
            r = CalcP(n, t);
            if (r != 0) {
                f = f + " " + r.ToString();
                t = n / r;
                n = n / r + 1;
            }
        }
        return f;
    }

    public static Int64 CalcP(Int64 num, Int64 tot) {
        for (Int64 i = 2; i <= tot; i++) {
            if (num % i == 0) {
                return i;
            }
        }
        return 0;
    }

    public static string RemSp(string str) {
        return new string(str.ToCharArray().Where(c => !Char.IsWhiteSpace(c)).ToArray());
    }

    public static bool CheckPrime(Int64 num, string s) {
        if (s == "") {
            return false;
        } else if (num == Convert.ToInt64(RemSp(s))) {
            return true;
        }
        return false;
    }

}

```

This is a more optimized version of the previous code, with some unnecessary redundant parts removed.

Output (on my i7 laptop):

```
Enter Number: 16
Start Time: 18:09:51.636445
Number, Factors
16,  2 2 2 2
2222,  2 11 101
211101,  3 11 6397
3116397,  3 163 6373
31636373 is prime
End Time: 18:09:51.637954

```

Test online

# Answer
> 2 votes
# Python 2 + primefac 1.1

I don't have a Raspberry Pi to test it on.

```
from primefac import primefac

def HP(n):
    factors = list(primefac(n))

    #print n, factors

    if len(factors) == 1 and n in factors:
        return n

    n = ""
    for f in sorted(factors):
        n += str(f)
    return HP(int(n))

```

**Try it online**

The `primefac` function returns a list of all prime factors of `n`. In its definition, it calls `isprime(n)`, which uses a combination of trial division, Euler's method, and the Miller-Rabin primality test. I'd recommend downloading the package and viewing the source.

I tried using `n = n * 10 ** int(floor(log10(f))+1) + f` instead of string concatenation, but it's much slower.

# Answer
> 1 votes
# Perl + ntheory, HP(80) in 0.35s on PC

No Raspberry Pi on hand.

```
use ntheory ":all";
use feature "say";
sub hp {
  my $n = shift;
  while (!is_prime($n)) {
    $n = join "",factor($n);
  }
  $n;
}
say hp($_) for (16,20,64,65,80);

```

The primality test is ES BPSW, plus a single random-base M-R for larger numbers. At this size we could use `is_provable_prime` (n-1 and/or ECPP) with no noticeable difference in speed, but that would change for \>300-digit values with no real benefit. Factoring includes trial, power, Rho-Brent, P-1, SQUFOF, ECM, QS depending on the size.

For these inputs it runs about the same speed as Charles' Pari/GP code on the OEIS site. ntheory has faster factoring for small numbers, and my P-1 and ECM are pretty good, but the QS isn't great so I would expect Pari to be faster at some point.

# Answer
> 0 votes
# pari/gp, HP(80) in 0.279 s on my laptop

Try it online!

```
homeprime(n) = {
    while(!isprime(n),
        fact = factor(n);
        n_str = "";
        for(i=1, matsize(fact)[1],
            for(j=1, fact[i,2],
                n_str = concat(n_str, Str(fact[i, 1]));
            );
        );
        n = eval(n_str);
    );
    return(n);
}

input_list = [16, 20, 64, 65, 80];
{
for(i=1, #input_list,
    n = input_list[i];
    total = 0.0;
    count = 10; \\ set the number of times you want to repeat the function call
    for(j=1, count,
        t0 = gettime();
        print(homeprime(n));
        t1 = gettime();
        taken = (t1 - t0); \\ gettime() returns milliseconds
        total += taken;
    );
    print("\t\t", total / count, "\n");
)
}

```

---
Tags: math, sequence, primes, fastest-code
---


# Question
Title: Find the first run of numbers summing to n
Given as input a positive nonzero integer `n` \>= 10 and a sequence of digits `0-9` (which may be taken as a string or a list), find the first contiguous subsequence of digits in the sequence that sums to `n` and output the start and end indexes. You may use zero- or one-based indexing. If no such subsequence exists, your program may output any *constant* value.

## Examples

These examples use zero-based indexing.

```
Input: 10 123456789
Output: 0 3
Input: 32 444444444
Output: 0 7
Input: 15 123456789
Output: 0 4
Input: 33 444444444
Output: No solutions

```

This is , so shortest program wins!

# Answer
> 6 votes
# Excel (ms365), 155 bytes

Assuming list as input and 0-indexed output:

Formula in `C1`:

```
=LET(x,TOCOL(B:B,1),y,COUNT(x),REDUCE("No solutions",SEQUENCE(y),LAMBDA(a,b,IFERROR(HSTACK(y-b,y-b+MATCH(A1,SCAN(0,DROP(x,y-b),LAMBDA(c,d,c+d)),0)-1),a))))

```

The idea here is to use `SCAN()` to iterate from each element in reversed order in column B and `MATCH()` the value in `A1` against each resulting array. If found, then `HSTACK()` both the current iteration and the position the element is found.

Unfortunately this process is rather verbose in Excel, but fun to attempt this challenge nonetheless.

---

# Answer
> 5 votes
# Vyxal `g`, 11 bytes

```
ẏ:Ẋ'ƒṡ¹İ∑⁰=

```

Try it Online!

A simple little approach.

## Explained

```
ẏ:Ẋ'ƒṡ¹İ∑⁰=
ẏ:              # The range [0, len(input)), twice.
  Ẋ             # cartesian product of the two lists and sort to get the first
   '            # get all items of that where
    ƒṡ          # the list turned into an inclusive range between the two numbers
      ¹İ        # indexed into the sequence of numbers
        ∑⁰=     # summed equals n
# g flag gets the smallest pair

```

# Answer
> 5 votes
# \>\<\> (Fish), 126 bytes

```
i01:}r[l0$>:?v~]}r4[{:}]=?v1+$:@$:@l5-$-)?v20.
/[-2lr@@+@:$@/    ;n+oan:$/-3lr1;?=-2l:+1~/
\}]r1-a0.                 \[}]r20.

```

Try it

## Explanation

Takes the array of characters as the initial value of the stack and the target as the input over STDIN. Outputs nothing if no solution exists, and an exclusive range if it does.

`i01` push the target value, the start of the range, and the length to the stack.

`:}r[l0$>:` Push the top "length" elements of the stack to a new stack. Set it's length, that's the amount of elements we want to compare.

`@$:@+@@rl2-]}]r1-a0.` This adds the top of the stack to the accumulator, rotates the stack, then subtract one from the number of remaining rotations necessary.

`~]}` pop 0, the remaining number of rotations, and combine the stack back with the previous. Now the sum of the top N elements of the input are on top of the stack.

`r4[{:}]=` Copy the target sum to the end of the input.

```
=?v
   \$:nao+n;

```

If the sum equals the target, output.

`1+$:@$:@15-$-)?v` Add one to the length. If it's greater than the remaining length of the list go down.

`~1+:l2-=?;` Add one to the starting value. If it's greater than the length of the input exit.

`1rl3-[}]r` Set the length back to 1. Shift the rest of the stack over one so it starts at the second position.

# Answer
> 5 votes
# 05AB1E, 11 bytes

```
ā<ã.Δ¹sŸèOQ

```

Inputs in the order `sequence, n`.
0-based output; will output `-1` if no result can be found; and `[i,i]` if `n` is a single digit that's present in the input-`sequence`.

Try it online or verify all test cases.

**Explanation:**

```
ā            # Push a list in the range [1, length of the first (implicit) input sequence]
 <           # Decrease each by to make it a 0-based [1,length)
  ã          # Cartesian power of 2 to get all possible pairs of this list
   .Δ        # Find the first pair that's truthy for:
     ¹       #  Push the first input-sequence again
      s      #  Swap so the current pair is at the top
       Ÿ     #  Pop and push a list in the range of this pair
        Ù    #  Uniquify this list (in case the pair was [a,a])
         è   #  Index each into the first input-sequence
          O  #  Sum these digits together
           Q #  Check if this sum is equal to the second (implicit) input
             # (after which the result is output implicitly)

```

# Answer
> 5 votes
# Python 3.8 (pre-release), 75 bytes

```
f=lambda v,l,a=0,b=1:v==(q:=sum(l[a:b]))and(a,b-1)or f(v,l,a+(q>v),b+(q<v))

```

Try it online!

Switched to a list of integers instead of a string for -9 bytes

old:

```
f=lambda v,l,a=0,b=1:v==(q:=sum(map(int,l[a:b])))and(a,b-1)or f(v,l,a+(q>v),b+(q<v))

```

Try it online!

Takes an integer and a string as input. Returns a tuple with start and end index if a solution is found, and an error if not.

# Answer
> 4 votes
# JavaScript (ES6), 64 bytes

*-1 thanks to @l4m2*

Expects `(n)(array)`. Returns either `[start, end]` or *false*.

```
n=>a=>a.some((_,i)=>a.some((v,j)=>j<i||(s-=v)?0:o=[i,j],s=n))&&o

```

Try it online!

### Commented

```
n =>               // n = target sum
a =>               // a[] = array of digits
a.some((_, i) =>   // for each digit at index i in a[]:
  a.some((v, j) => //   for each digit v at index j in a[]:
    j < i ||       //     if j is less than i
    (s -= v)       //     or subtracting v from s
    ?              //     does not result in 0:
      0            //       do nothing
    :              //     else:
      o = [i, j],  //       set o = [i, j] and exit
    s = n          //     start with s = n
  )                //   end of inner some()
)                  // end of outer some()
&& o               // if successful, return o[]

```

# Answer
> 4 votes
# Japt, 15 bytes

```
ðÊï æÈrõ xgU ¶V

```

Try it

```
ðÊï æÈrõ xgU ¶V     :Implicit input of digit array U and integer V
ð                   :0-based indices of elements of U that are truthy under
 Ê                  :  Factorial
  ï                 :Cartesian product with itself
    æ               :First element/pair that returns true when
     È              :Passed through the following function
      r             :  Reduce by
       õ            :    Inclusive range
         x          :  Reduce by addition after
          gU        :    Indexing each into U
             ¶V     :  Is equal to V?

```

# Answer
> 4 votes
# R 4.1, ~~78 bytes~~ ~~71 bytes~~ 65 bytes

## v1.2 65 bytes

-5 bytes thanks to pajonk:

1. -4 for pointing out that `F <- F + 1;k <- m[, F];sum(s[k:k[2]])` can be expressed as `sum(s[(k <- m[, F <- F + 1]):k[2]])`. This also means you can lose the curly braces as the expression is on one line, and you do not need semi-colons.
2. -1 because `while(a!=b)` can actually be expressed as `while(a-b)`. This is safe as negative numbers are also `TRUE` when coerced to `logical`.

(And also -1 byte by pointing out I had miscounted.)

```
\(s,n,m=combn(seq(s),2)){while(sum(s[(k=m[,F<-F+1]):k[2]])-n)0
k}

```

\[Attempt This Online!\]https://ato.pxeger.com/run?1=m72waMGSNAUb3aWlJWm6FjcdYzSKdfJ0cm2T83OT8jSKUws1ijV1jDQ1q8szMnNSNYpLczWKozWybXOjddxsdN20DWM1rbKjjWJjNXXzNA24smuhBrkpK4SkFpcoJCcWpxZzpWkYWlnqKBgaaAKZRakFGiZgjo6CsZEmXNIUQ9JYE2LaggUQGgA)

## v 1.1 71 bytes

```
\(s,n,m=combn(seq(s),2)){while({F=F+1;k<-m[,F];sum(s[k:k[2]])!=n})0
k}

```

-7 bytes thanks to Giuseppe

Two changes:

1. Not replacing `seq()` with the unary `-` as it's only being used once (I was using it twice before and failed to notice that was no longer the case).
2. Replacing `j` with `F` is yet another ludicrous thing that R allows you to do that I would never have tried in the vein of these coercion tips by @rturnbull (link from Giuseppe - thanks). For non-R uses, `T` and `F` are built-in aliases for TRUE and FALSE. This means you can do this:

```
T == TRUE # TRUE
F == FALSE # TRUE
F == TRUE # FALSE
F  <- F+1
F == FALSE # FALSE
F == TRUE # TRUE

```

I have no idea why this is allowed in R but it's saved 5 characters compared with initialising `j`.

## v.1.0 78 bytes

```
\(s,n,`-`=seq,j=0,m=combn(-s,2)){while({j=j+1;k<-m[,j];sum(s[k:k[2]])!=n})0
k}

```

Attempt this online

### How?

This is a function which expects a vector, `s` and a target value `n`, as arguments. Basically it does this:

1. Uses `combn()` to create a `m`, a matrix of all possible length 2 combinations of the indices of `s`. E.g. `combn(1:5,2)` produces:

```
     [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,]    1    1    1    1    2    2    2    3    3     4
[2,]    2    3    4    5    3    4    5    4    5     5

```
2. Loops through the columns of this to subset `s`, using row 1 as the start index and row 2 as the end index.
3. Breaks the loop if the subset adds up to `n`. It exits with an error if there is no solution. I wasn't entirely clear if an error counts as a constant, but this nice Python answer behaves in the same way, which I took as a green light.
4. There are a couple of golfing tricks both from Giuseppe. Taking advantage of how `seq(s)` behaves if `s` is a vector and reassigning the unary `-` with ``-` = seq`, means we can do `m = combn(-s, 2)` instead of the positively verbose `m = combn(seq_along(s), 2)`.

The test cases are included in the footer of ATO link. Only my second Code Golf answer ever so very interested in any suggestions / improvements.

# Answer
> 3 votes
# Perl 5 List::Util, 82 bytes

```
sub{$s=pop;($i,$j)=($_%@_,int$_/@_),$s-sum(@_[$i..$j])||return($i,$j)for 0..@_*@_}

```

Try it online!

# Answer
> 3 votes
# Pyth, 18 bytes

```
KEhfqQs:KhTheT^UK2

```

Try it online!

Takes n, sequence as separate inputs. Outputs \[start, end\] or an index error if no sequence exists.

### Explanation

```
                      # implicitly assign Q = eval(input())
KE                    # assign K = eval(input())
              ^  2    # cartesian product of two copies of
               UK     # range(len(K))
   f                  # filter this on lambda T
      s               #   sum of
       :KhTheT        #   slice K from first element of product to second element plus 1
    qQ                #   equals Q
  h                   # take the first element (this will give us the indices if they exist or an index error if they don't)

```

# Answer
> 3 votes
# C (clang), <sup>~~90~~</sup> 88 bytes

```
i;j;a;f(*s,n){for(a=i=-1;s[++i]*a;)for(a=n,j=i;a*s[j];)a-=s[j++]-48;*s++=!a*i-1;*s=j-1;}

```

Try it online!

*Saved 2 bytes thanks to ceilingcat!!!*

Inputs a pointer to a wide charater string of the sequence and \\$n\\$.
Returns the first and last position in the array referenced by the input pointer or \\$-1\\$ if no sum exists.

# Answer
> 3 votes
# Jelly, 12 bytes

```
Jṗ2r/ị³S=ʋƇḢ

```

Try it online!

Full program, that takes a list of digits and an integer \\$n\\$, and outputs the 1-indexed indices, or `0` if there's no solution.

## How it works

```
Jṗ2r/ị³S=ʋƇḢ - Main link. Takes a list of digits D on the left, and n on the right
J            - Indices; [1, 2, ..., len(D)]
 ṗ2          - Cartesian square
         ʋƇ  - Filter the indices on the dyad f([i, j], n):
   r/        -   Range; [i, i+1, ..., j-1, j]
     ị³      -   Index into D
       S     -   Sum
        =    -   Does that equal n?
           Ḣ - Take the first pair, or zero if there's no solution

```

# Answer
> 2 votes
# Thunno, \\$ 25 \log\_{256}(96) \approx \\$ 20.58 bytes

```
LRDZ!gAu1+:z0sAISz1=kAJ0~

```

Takes a list of digits, then a number. Outputs 0 if it doesn't find the subsequence.

#### Explanation

```
LRDZ!gAu1+:z0sAISz1=kAJ0~  # Implicit input
LR                         # range(len(input))
  DZ!                      # Cartesian product with itself
     g              k      # Filter for items where the following is truthy:
      Au                   #  Dump the pair onto the stack
        1+                 #  Add one to the last one (because of how Thunno's range works)
          :                #  Range between the two numbers
             sAI           #  Indexed into
           z0              #  The first input (this will return a list)
                S          #  The sum of this list
                 z1=       #  Equals the second input?
                     AJ    # After the filter, get the first item of the list
                       0~  # If the list was empty, push 0

```

#### Screenshots

# Answer
> 2 votes
# C (clang), 76 bytes

Takes a list of digits and a number as input and returns \\$\[start, end)\\$.

```
i;j;f(*l,s,n){for(i=j=0;n<0|s/~j*n;)n-=n>0?l[j++]:-l[i++];l[1]=j;*l=n?-1:i;}

```

## How it works

```
i;j;             // start and end of the sliding window, respectively.
f(*l, s, n) {    // input a list of integers, its size and the number.
  for(           // main loop
    i = j = 0;   // start with i = 0 and j = 0 (empty window)
    n < 0 |      // the sum of elements in the window exceeds n, we will remove elements in the next iteration
    s / ~j * n;) // j < l and the sum of elements does not match n
      n -=       // subtract from n:
        n > 0 ?  // if n > 0 ...
          l[j++] // increase the window size.
        :        // else (n < 0) ...
         -l[i++] // decrease the window size.
      ;          // end main loop
  l[1] = j;      // l doubles as return placeholder. l[1] is the end of the window.
  *l =           // l[0] is:
    n ?          // if n != 0 ...
      -1         // there was no solution, set it to -1.
    :            // else (n == 0) ...
      i;         // we found the solution, set it to the start of the window.
}

```

Try it online!

# Answer
> 2 votes
# Desmos, 79 bytes

```
L=[1...l.length]
T=[0^{(l[I...E].total-s)^2}(I,E)forE=L,I=L]
f(s,l)=T[T.x>0][1]

```

\\$f\\$ returns a point `(I,E)` that has the starting index `I` and the ending index `E` as the \\$x\\$ and \\$y\\$ coordinates, respectively, and returns the answer one-indexed.

For no solution, \\$f\\$ returns the point `(undefined,undefined)`.

Note that because of the limitations of Desmos, this will not work for input lists longer than \\$100\\$ elements.

Try It On Desmos!

Try It On Desmos! - Prettified

# Answer
> 2 votes
# Julia, ~~69~~ ~~65~~ 47 bytes

```
~=:;y^x=argmax(k->sum(x[k])==y,(q=keys(x)).~q')

```

Attempt This Online!

Uses 1-based indexing. Prints the range `1:1` if no solution exists.

* -4 bytes thanks to MarcMush: Replace `1:length(x)` with `keys(x)`
* -18 bytes (!!) thanks to MarcMush: Improve range generation

# Answer
> 1 votes
# Charcoal, 25 bytes

```
ＮθＩ⌊ΦＥη⟦κ⌕ＥηΣ✂ηκ⊕μ¹θ⟧⊕↨ι⁰

```

Try it online! Link is to verbose version of code. Outputs `None` if there is no solution. Explanation:

```
Ｎθ                          First input as a number
      η                     Second input
     Ｅ                     Map over digits
        κ                   Current index
           η                Second input
          Ｅ                 Map over digits
              η             Second input
             ✂    ¹         Sliced from
               κ            Outer index to
                 μ          Inner index
                ⊕           Incremented
            Σ               Take the sum
         ⌕                  Find index of
                   θ        First input
       ⟦            ⟧       Make into list
    Φ                       Filtered where
                      ↨ι⁰   Last element
                     ⊕      Is not `-1`
   ⌊                        Take the minimum (first if it exists)
  Ｉ                         Cast to string
                            Implicitly print

```

# Answer
> 1 votes
# Factor \+ `math.combinatorics math.unicode`, 69 bytes

```
[| n s | s length iota 2 [ last2 s subseq Σ n = ] find-combination ]

```

Try it online!

Returns a 0-based start index and 1-based end index. Returns `f` if there is no solution.

# Answer
> 1 votes
## C99:

```
int main() {
    int x = 22;
    int a[] = {1, 2, 3, 4, 5, 6, 7, 8, 9}, n = sizeof(a) / sizeof(int);

    int i = 0, j = 0, s = 0;
    while ((s < x && j < n) || (s > x && i < n))
        s += s < x ? a[j++] : -a[i++];

    printf("found: %d, i: %d, j: %d", s == x, i, j - 1);
    return 0;
}

```

The algorithm actually takes 66 bytes:

```
inti=0,j=0,s=0;while((s<x&&j<n)||(s>x&&i<n))s+=s<x?a[j++]:-a[i++];

```

Try it!

# Answer
> 1 votes
# Haskell, 71 bytes

```
x#y=take 1[[s,e]|i<-[[0..length x-1]],s<-i,e<-i,sum[x!!i|i<-[s..e]]==y]

```

Try it online!

# Haskell, 77 bytes

```
x#y=take 1[[s,e]|let m=length x-1,s<-[0..m],e<-[s..m],sum[x!!i|i<-[s..e]]==y]

```

Try it online!

# Answer
> 1 votes
# Ruby, ~~58~~ 56 bytes

Returns `nil` if no subsequence is found.

-2 bytes thanks to Value Ink

```
->n,a{(b=*0...a.size).product(b).find{a[_1.._2].sum==n}}

```

Attempt This Online!

# Answer
> 0 votes
# Arturo, 79 bytes

```
$[n,a][l:-size a 1loop 0..l'x[loop x..l'y[if n=sum slice a x y->return[x,y]]]0]

```

Try it

Returns `[start end+1]` or `0` if there's no solution.

# Answer
> 0 votes
# Scala, 91 bytes

Golfed version. Try it online!

```
(a,t)=>a.indices.flatMap(s=>(s until a.length).filter(e=>a.slice(s,e+1).sum==t).map((s,_)))

```

Ungolfed version. Try it online!

```
object Main {

  def findSubarrays(arr: Array[Int], target: Int): Seq[(Int, Int)] = {
    for {
      start <- arr.indices
      end <- start until arr.length
      if arr.slice(start, end + 1).sum == target
    } yield (start, end)
  }

  def main(args: Array[String]): Unit = {
    val numbers = "12345".map(_.asDigit).toArray
    val target = 12
    println(findSubarrays(numbers, target))
  }

}

```

---
Tags: code-golf, number, subsequence
---


title: Fastest code to find the next prime

text:
# Question
Title: Fastest code to find the next prime
The problem is as follows.

**Input:** An integer `n`

**Output:** The smallest prime bigger than `n`.

The challenge is to give the fastest code possible to do this. I will test the code on values starting at size  roughly `10^8` `10^200` and doubling in size until it takes more than one minute 10 seconds on my computer.

**The winning code** will find the next prime for the largest input size.

By way of comparison, a simple sieve written in python is able to find the next prime bigger than `10^8` in about `20` seconds.

The requirement that I can test it on my 4GB RAM ubuntu computer is strict. All code must be free (in both senses) and if it uses libraries they must also be free and easily installable. Any false primes reported will immediately disqualify the submission.

I will award separate commendations for the winners in each programming language too if the code is entirely written in that language without the use of external libraries. I will also keep a running table of the fastest times as the competition goes on so people can see how they are doing.

**Table so far**

* **Python.** An astonishing `357` digit prime `343239883006530485749095039954069660863471765007165270469723172959277159169882802606127982033072727748864815569574042901856099399985832190628701414555752857600000000000000000000000000000000000000002872284792758930912601189043411951050852357613658978971208596097634095500808832510259693761982135208603287199546795000697807728609476163156438356035166156820611` was the final number under 10 seconds using the code supplied by `primo`. Will anyone beat this first entry?

# Answer
> 23 votes
## Python ~451 digits

This is part of the library I wrote for a semiprime factorization problem, with unnecessary functions removed. It uses the Baillie-PSW primality test, which is technically a probabilistic test, but to date, there are no known pseudoprimes – and there's even a cash reward if you're able to find one (or for supplying a proof that none exist).

**Edit**: I hadn't realized that Python has built-in modular exponentiation. Replacing my own for the built-in results in a performance boost of about 33%.

my\_math.py

```
# 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)

# strong probable prime
def is_sprp(n, b=2):
  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 range(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 > 0:
    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 > 0:
    if t&1 == 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 > 0:
      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

# primes less than 212
small_primes = set([
    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]

max_int = 2147483647

# an 'almost certain' primality check
def is_prime(n):
  if n < 212:
    return n in small_primes

  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:
    i = 211
    while i*i < n:
      for o in offsets:
        i += o
        if n%i == 0:
          return False
    return True

  # Baillie-PSW
  # this is technically a probabalistic test, but there are no known pseudoprimes
  if not is_sprp(n): return False
  a = 5
  s = 2
  while legendre(a, n) != n-1:
    s = -s
    a = s-a
  return is_lucas_prp(n, a)

# next prime strictly larger than n
def next_prime(n):
  if n < 2:
    return 2
  # first odd larger than n
  n = (n + 1) | 1
  if n < 212:
    while True:
      if n in small_primes:
        return n
      n += 2

  # find our position in the sieve rotation via binary search
  x = int(n%210)
  s = 0
  e = 47
  m = 24
  while m != e:
    if indices[m] < x:
      s = m
      m = (s + e + 1) >> 1
    else:
      e = m
      m = (s + e) >> 1

  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

```

A sample test script:

```
from time import clock
from my_math import *

n = i = 317**79
while True:
  i *= 317
  time1 = clock()
  n, o = next_prime(i), n
  span = clock()-time1
  if span > 10:
    break
  print(len(str(n)), span)
print(o)

```

A factor of 317 was chosen, because it's approximately the square root of `10000`, adding roughly 2.5 digits per iteration (and because doubling was too slow to sit through). Output shows the current number of digits, and the time taken.

Sample results:

```
201 0.13121248650317288
203 0.059535499623555505
206 0.9157767258129175
208 0.2583420518529589
211 0.15367400046653978
213 0.32343915218274955
216 1.3962866788935466
218 0.5986165839513125
221 0.973842206202185
223 2.346910291671148
...
428 0.932809896229827
431 4.345940056627313
433 9.511724255457068
436 6.089835998709333
438 1.3793498894412721
441 4.290633027381972
443 3.5102506044762833
446 3.1629148397352083
448 3.364759208223404
451 7.34668009481652
1551197868099891386459896063244381932060770425565921999885096817830297496627504652115239001983985153119775350914638552307445919773021758654815641382344720913548160379485681746575245251059529720935264144339378936233043585239478807971817857394193701584822359805681429741446927344534491412763713568490429195862973508863067230162660278070962484418979417980291904500349345162151774412157280412235743457342694749679453616265540134456421369622519723266737913

```

All code is now python 3 compatible.

# Answer
> 6 votes
## C++ with GMP: 567 digits

Uses the Miller-Rabin implementation in GMP. It might return a false positive, but good luck actually hitting one with probability 2^-200.

```
#include <gmp.h>
#include <stdlib.h>
#include <stdio.h>
#include <sys/time.h>

double time() {
  struct timeval t;
  gettimeofday(&t, NULL);
  return t.tv_usec  * 1e-6 + t.tv_sec;
}

int main(int argc, char *argv[]) {
  mpz_t n, m;
  mpz_init_set_ui(n, 10);
  mpz_pow_ui(n, n, 200);
  mpz_init(m);
  for (int i = 0; true; i++, mpz_mul_ui(n, n, 2)) {
    double start = time();
    for (mpz_add_ui(m, n, 1); !mpz_millerrabin(m, 100); mpz_add_ui(m, m, 2)) ;
    double t = time() - start;
    gmp_printf("%d %Zd %f\n", i, m, t);
    if (t > 10.0) break;
  }
}

```

Finds the prime `10^200 * 2^1216 + 361` (567 digits) before running over time on my slow laptop.

# Answer
> 5 votes
# Perl with GMP module, 1300 digits

Using my module Math::Prime::Util and its GMP back end. The exact crossover point will depend on your computer and whether you have the latest GMP library. All code is free (the modules are on github and CPAN, and GMP is freely available). I've run them on AWS's Ubuntu as well as a desktop Ubuntu (and Fedora, and AIX, and NetBSD, etc.).

The core code is in C and C+GMP. next\_prime from MPU sees the number is too large and forwards it to the GMP back end (or pure Perl code if the back end isn't installed). That stringiifies and converts to a mpz, and will convert the result back into the input object type or Math::BigInt. next\_prime itself does:

* a mod 30 wheel
* keeps track of the remainder mod 23# so it can do native modulos for primes up to 23
* probable prime test on things that pass these.

The probable prime test is:

* check for tiny divisors using mpz\_gcd\_ui (in 64-bit two of these check up to 101)
* check for small divisors using singly-calculated large primorials. This is either primes up to 10k or 40k depending on the input size.
* for values larger than 2^1600, performs further trial division using a treesieve. This could be done more efficiently.
* finally, ES BPSW is done (Miller-Rabin test with base 2 followed by an extra strong Lucas test).

Everything prior to the ES BPSW is just optimization, which of course we want for next\_prime. next\_prime is also implemented in Perl using the Math::BigInt module (in core with optional Pari and GMP back ends). That does AES BPSW (like Pari) but isn't as optimized.

I have thought about the merits of a partial-sieve-based version, using a range of, for instance, 2 merits. I'm just not sure if this would really be better, as most of the time we'd be doing unnecessary sieving as the gap was small, and sometimes for a large gap we'd have to repeat it multiple times.

The library implements ECPP (including certificates) so we could run a proof on the result, but 1200 digits is really too large for the tiny default set of included polynomials (there is a method for downloading larger sets -- proofs would take a bit under 15 minutes which is a bit faster than Pari's APR-CL but a bit slower than WraithX's mpz\_aprcl). One downside of ECPP vs. APR-CL is that it has more time variance so there is a good chance of it exceeding 10 seconds on some number before the average time gets there. With a proof I think we're limited to something in the 400 digit range unless we allow multi-threaded software.

```
#!/usr/bin/env perl
use warnings;
use strict;
use Math::Prime::Util ":all";
use Math::Prime::Util::GMP;  # Barf if the backend isn't installed
use Time::HiRes qw(gettimeofday tv_interval);
use Math::GMP;

my $n = Math::GMP->new(10) ** 200;
while (1) {
  my $start = [gettimeofday];
  my $np = next_prime($n);
  my $sec = tv_interval($start);
  my $len = length($n);
  die "next_prime $len = +",$np-$n," in $sec seconds\n" if $sec > 10;
  warn "  next_prime $len = +",$np-$n," in $sec seconds\n";
  $n *= 10;
}

```

I decided to try with the same sequence used by primo. It got to 1191 digits, as that is where we hit a gap of 18138. I also tested primo's code using the latest my\_math.py. It gets to 630 digits with the 10^e sequence and 641 with his sequence. Very impressive for compact all-Python code without lots of pretests.

# Answer
> 1 votes
# pari/gp, 703 digits

Try it online!

```
n = i = 317^79;

{
while(1,
  i *= 317;
  time1 = gettime();
  n = nextprime(i);
  o = n;
  span = gettime() - time1;

  if (span > 10*1000, break);

  print("Length: ", #Str(n), ", Time: ", span, "ms");
);
}

print(o);

```

---
Tags: fastest-code, primes
---


title: Compute OEIS A005434

text:
# Question
Title: Compute OEIS A005434
The task is to compute OEIS A005434 as quickly as possible.

Consider a binary string `S` of length `n`. Indexing from `1`, we can determine if `S[1..i+1]` matches `S[n-i..n]` exactly for all `i` in order from `0` to `n-1`. For example,

```
S = 01010

```

gives

```
[Y, N, Y, N, Y].

```

This is because `0` matches `0`, `01` does not match `10`, `010` matches `010`, `0101` does not match `1010` and finally `01010` matches itself.

Define `f(n)` to be the number of distinct arrays of `Y`s and `N`s one gets when iterating over all `2^n` different possible bit strings `S` of length `n`.

The observant will notice this question is a simpler variant of another recent question of mine. However, I expect that clever tricks can make this much faster and easier.

## Task

For increasing `n` starting at `1`, your code should output `n, f(n)`.

## Example answers

For `n = 1..24`, the correct answers are:

```
1, 2, 3, 4, 6, 8, 10, 13, 17, 21, 27, 30, 37, 47, 57, 62, 75, 87, 102, 116, 135, 155, 180, 194

```

## Scoring

Your code should iterate up from `n = 1` giving the answer for each `n` in turn. I will time the entire run, killing it after two minutes.

Your score is the highest `n` you get to in that time.

In the case of a tie, the first answer wins.

## Where will my code be tested?

I will run your code under Virtualbox in a Lubuntu guest VM (on my Windows 7 host).

My laptop has 8GB of RAM and an Intel i7 5600U@2.6 GHz (Broadwell) CPU with 2 cores and 4 threads. The instruction set includes SSE4.2, AVX, AVX2, FMA3 and TSX.

## Leading entries per language

* **n = 599** in **Rust** bu Anders Kaseorg.
* **n= 30** in **C** by Grimy. Parallel version gets to **32** when run natively in cygwin.

# Answer
# Rust, n ≈ 660

```
use std::collections::HashMap;
use std::iter::once;
use std::rc::Rc;

type Memo = HashMap<(u32, u32, Rc<Vec<u32>>), u64>;

fn f(memo: &mut Memo, mut n: u32, p: u32, mut s: Rc<Vec<u32>>) -> u64 {
    debug_assert!(p != 0);
    let d = n / p;
    debug_assert!(d >= 1);
    let r = n - p * if d >= 2 { d - 1 } else { 1 };

    let k = s.binary_search(&(n - r + 1)).unwrap_or_else(|i| i);
    for &i in &s[..k] {
        if i % p != 0 {
            return 0;
        }
    }

    if d >= 3 {
        let o = n - (p + r);
        n = p + r;
        s = Rc::new(s[k..].iter().map(|i| i - o).collect());
    } else if n == p {
        return 1;
    } else if k != 0 {
        s = Rc::new(s[k..].to_vec());
    }

    let query = (n, p, s);
    if let Some(&c) = memo.get(&query) {
        return c;
    }
    let (n, p, s) = query;

    let t = Rc::new(s.iter().map(|i| i - p).collect::<Vec<_>>());
    let c = if d < 2 {
        (1..r + 1).map(|q| f(memo, r, q, t.clone())).sum()
    } else if r == p {
        (1..p + 1)
            .filter(|&q| p % q != 0 || q == p)
            .map(|q| f(memo, r, q, t.clone()))
            .sum()
    } else {
        let t = match t.binary_search(&p) {
            Ok(_) => t,
            Err(k) => {
                Rc::new(t[..k]
                            .iter()
                            .cloned()
                            .chain(once(p))
                            .chain(t[k..].iter().cloned())
                            .collect::<Vec<_>>())
            }
        };
        (1..t.first().unwrap() + 1)
            .filter(|&q| p % q != 0 || q == p)
            .map(|q| f(memo, r, q, t.clone()))
            .sum()
    };
    memo.insert((n, p, s), c);
    c
}

fn main() {
    let mut memo = HashMap::new();
    let s = Rc::new(Vec::new());
    for n in 1.. {
        println!("{} {}",
                 n,
                 (1..n + 1)
                     .map(|p| f(&mut memo, n, p, s.clone()))
                     .sum::<u64>());
    }
}

```

Try it online!

### How it works

This is a memoized implementation of the recursive predicate Ξ given in Leo Guibas, “Periods in strings” (1981). The function `f(memo, n, p, s)` finds the number of correlations of length `n` with smallest period `p` and also each of the periods in the set `s`.
> 4 votes

# Answer
Just a simple brute-force search, to get the challenge started:

```
#include <stdio.h>
#include <stdint.h>
#include <string.h>

typedef uint16_t u16;
typedef uint64_t u64;

static u64 map[1<<16];

int main(void)
{
    for (u64 n = 1;; ++n) {
        u64 result = 1;
        u64 mask = (1ul << n) - 1;
        memset(map, 0, sizeof(map));

        #pragma omp parallel
        #pragma omp for
        for (u64 x = 1ul << (n - 1); x < 1ul << n; ++x) {

            u64 r = 0;
            for (u64 i = 1; i < n; ++i)
                r |= (u64) (x >> i == (x & (mask >> i))) << i;
            if (!r)
                continue;

            u16 h = (u16) (r ^ r >> 13 ^ r >> 27);
            while (map[h] && map[h] != r)
                ++h;

            if (!map[h]) {
                #pragma omp critical
                if (!map[h]) {
                    map[h] = r;
                    ++result;
                }
            }
        }

        printf("%ld\n", result);
    }
}

```

Compile with `clang -fopenmp -Weverything -O3 -march=native`. On my machine it reaches n=34 in 2 minutes.

EDIT: sprinkled some OMP directives for easy parallelism.
> 2 votes

# Answer
# Python, n≈427

Rewrite @Anders Kaseorg's Rust answer in Python

TIO Link

```
import bisect
from collections import defaultdict
import time

memo = defaultdict(int)
s = []

def f(n, p, s):
    assert(p != 0)
    d = n // p
    assert(d >= 1)
    r = n - p * (d - 1 if d >= 2 else 1)

    k = bisect.bisect_left(s, n - r + 1)
    for i in s[:k]:
        if i % p != 0:
            return 0

    if d >= 3:
        o = n - (p + r)
        n = p + r
        s = [i - o for i in s[k:]]
    elif n == p:
        return 1
    elif k != 0:
        s = s[k:]

    query = (n, p, tuple(s))
    if query in memo:
        return memo[query]

    t = [i - p for i in s]
    if d < 2:
        c = sum(f(r, q, t) for q in range(1, r + 1))
    elif r == p:
        c = sum(f(r, q, t) for q in range(1, p + 1) if p % q != 0 or q == p)
    else:
        k = bisect.bisect_left(t, p)
        if k < len(t) and t[k] == p:
            pass
        else:
            t.insert(k, p)
        c = sum(f(r, q, t) for q in range(1, t[0] + 1) if p % q != 0 or q == p)
    memo[query] = c
    return c

start_time = time.time()
n = 1
while time.time() - start_time < 120:  # run for at most 2 minutes
    print(n, sum(f(n, p, s) for p in range(1, n + 1)))
    n += 1

```
> 0 votes

---
Tags: array, combinatorics, fastest-code
---


# Question
Title: Find the maximum determinant for each size Toeplitz matrix
For a fixed n, consider the n by n Toeplitz matrices with entries which are either 0 or 1. The aim is to find maximum determinant over all such Toeplitz matrices.

**Task**

For each `n` from 1 upwards, output the maximum determinant over all n by n Toeplitz matrices with entries which are either 0 or 1. There should be one output per `n` which should have the maximum determinant and also an example matrix that reaches it.

**Score**

Your score is the largest `n` your code gets to in 2 minutes on my computer. To clarify a little, your code can run for 2 minutes in total, this is not 2 minutes per `n`.

**Tie breaker**

If two entries get the same `n` score then the winning entry will be the one that gets to the highest `n` in the shortest time on my machine. If the two best entries are equal on this criterion too then the winner will be the answer submitted first.

**Languages and libraries**

You can use any freely available language and libraries you like. I must be able to run your code so please include a full explanation for how to run/compile your code in linux if at all possible.

**My Machine** The timings will be run on my machine. This is a standard ubuntu install on an AMD FX-8350 Eight-Core Processor. This also means I need to be able to run your code.

**Small answers**

For n = 1..10 the outputs should be 1,1,2,3,5,9,32,56,125,315

This sequence is not in OEIS and so the winning entry also gets to propose a new entry there.

**Entries so far**

* `n=10` `n=11` by Vioz in **Python**
* `n=9` by Tyilo in **C**
* `n=12` by Legendre in **J**
* `n=10` by Tensibai in **R**
* `n=14` by SteelRaven in **C++**
* `n=14` by RetoKoradi in **C++**

# Answer
# C++ with pthreads

This gets to n=14 in just under 1 minute on my machine. But since that's just a 2-core laptop, I hope that the 8-core test machine can finish n=15 in under 2 minutes. It takes about 4:20 minutes on my machine.

I was really hoping to come up with something more efficient. There has **got** to be a way to calculate the determinate of a binary matrix more efficiently. I wanted to come up with some kind of dynamic programming approach that counts the +1 and -1 terms in the determinant calculation. But it just hasn't quite come together so far.

Since the bounty is about to expire, I implemented the standard brute force approach:

* Loop over all possible Toeplitz matrices.
* Skip one of the two in each transposed matrix pair. Since the matrix is described by bitmask values, this is simple to do by skipping all values where the reverse of the bitmask is smaller than the bitmask itself.
* The determinate is calculated with a text book LR decomposition. Except for some minor performance tuning, the main improvement I made to the algorithm from my college numerical methods book is that I use a simpler pivot strategy.
* Parallelization is done with pthreads. Just using regular spacing for the values processed by each thread caused very bad load balancing, so I introduced some swizzling.

I tested this on Mac OS, but I used similar code on Ubuntu before, so I hope this will compile and run without a hitch:

1. Save the code in a file with a `.cpp` extension, e.g. `optim.cpp`.
2. Compile with `gcc -Ofast optim.cpp -lpthread -lstdc++`.
3. Run with `time ./a.out 14 8`. The first argument is the maximum `n`. 14 should finish in under 2 minutes for sure, but it would be great if you could try 15 as well. The second argument is the number of threads. Using the same value as the number of cores of the machine is normally a good start, but trying some variations could potentially improve the times.

Let me know if you have any problems building or running the code.

```
#include <stdint.h>
#include <pthread.h>
#include <cstdlib>
#include <iostream>

static int NMax = 14;
static int ThreadCount = 4;

static pthread_mutex_t ThreadMutex;
static pthread_cond_t ThreadCond;
static int BarrierCount = 0;

static float* MaxDetA;
static uint32_t* MaxDescrA;

static inline float absVal(float val)
{
    return val < 0.0f ? -val : val;
}

static uint32_t reverse(int n, uint32_t descr)
{
    uint32_t descrRev = 0;
    for (int iBit = 0; iBit < 2 * n - 1; ++iBit)
    {
        descrRev <<= 1;
        descrRev |= descr & 1;
        descr >>= 1;
    }

    return descrRev;
}

static void buildMat(int n, float mat[], uint32_t descr)
{
    int iDiag;
    for (iDiag = 1 - n; iDiag < 0; ++iDiag)
    {
        float val = static_cast<float>(descr & 1);
        descr >>= 1;
        for (int iRow = 0; iRow < n + iDiag; ++iRow)
        {
            mat[iRow * (n + 1) - iDiag] = val;
        }
    }

    for ( ; iDiag < n; ++iDiag)
    {
        float val = static_cast<float>(descr & 1);
        descr >>= 1;
        for (int iCol = 0; iCol < n - iDiag; ++iCol)
        {
            mat[iCol * (n + 1) + iDiag * n] = val;
        }
    }
}

static float determinant(int n, float mat[])
{
    float det = 1.0f;
    for (int k = 0; k < n - 1; ++k)
    {
        float maxVal = 0.0f;
        int pk = 0;
        for (int i = k; i < n; ++i)
        {
            float q = absVal(mat[i * n + k]);
            if (q > maxVal)
            {
                maxVal = q;
                pk = i;
            }
        }

        if (pk != k)
        {
            det = -det;
            for (int j = 0; j < n; ++j)
            {
                float t = mat[k * n + j];
                mat[k * n + j] = mat[pk * n + j];
                mat[pk * n + j] = t;
            }
        }

        float s = mat[k * n + k];
        det *= s;

        s = 1.0f / s;
        for (int i = k + 1; i < n; ++i)
        {
            mat[i * n + k] *= s;
            for (int j = k + 1; j < n; ++j)
            {
                mat[i * n + j] -= mat[i * n + k] * mat[k * n + j];
            }
        }
    }

    det *= mat[n * n - 1];

    return det;
}

static void threadBarrier()
{
    pthread_mutex_lock(&ThreadMutex);

    ++BarrierCount;
    if (BarrierCount <= ThreadCount)
    {
        pthread_cond_wait(&ThreadCond, &ThreadMutex);
    }
    else
    {
        pthread_cond_broadcast(&ThreadCond);
        BarrierCount = 0;
    }

    pthread_mutex_unlock(&ThreadMutex);
}

static void* threadFunc(void* pData)
{
    int* pThreadIdx = static_cast<int*>(pData);
    int threadIdx = *pThreadIdx;

    float* mat = new float[NMax * NMax];

    for (int n = 1; n <= NMax; ++n)
    {
        uint32_t descrRange(1u << (2 * n - 1));
        float maxDet = 0.0f;
        uint32_t maxDescr = 0;

        uint32_t descrInc = threadIdx;
        for (uint32_t descrBase = 0;
             descrBase + descrInc < descrRange;
             descrBase += ThreadCount)
        {
            uint32_t descr = descrBase + descrInc;
            descrInc = (descrInc + 1) % ThreadCount;

            if (reverse(n, descr) > descr)
            {
                continue;
            }

            buildMat(n, mat, descr);
            float det = determinant(n, mat);
            if (det > maxDet)
            {
                maxDet = det;
                maxDescr = descr;
            }
        }

        MaxDetA[threadIdx] = maxDet;
        MaxDescrA[threadIdx] = maxDescr;

        threadBarrier();
        // Let main thread output results.
        threadBarrier();
    }

    delete[] mat;

    return 0;
}

static void printMat(int n, float mat[])
{
    for (int iRow = 0; iRow < n; ++iRow)
    {
        for (int iCol = 0; iCol < n; ++iCol)
        {
            std::cout << " " << mat[iRow * n + iCol];
        }
        std::cout << std::endl;
    }

    std::cout << std::endl;
}

int main(int argc, char* argv[])
{
    if (argc > 1)
    {
        NMax = atoi(argv[1]);
        if (NMax > 16)
        {
            NMax = 16;
        }
    }

    if (argc > 2)
    {
        ThreadCount = atoi(argv[2]);
    }

    MaxDetA = new float[ThreadCount];
    MaxDescrA = new uint32_t[ThreadCount];

    pthread_mutex_init(&ThreadMutex, 0);
    pthread_cond_init(&ThreadCond, 0);

    int* threadIdxA = new int[ThreadCount];
    pthread_t* threadA = new pthread_t[ThreadCount];

    for (int iThread = 0; iThread < ThreadCount; ++iThread)
    {
        threadIdxA[iThread] = iThread;
        pthread_create(threadA + iThread, 0, threadFunc, threadIdxA + iThread);
    }

    float* mat = new float[NMax * NMax];

    for (int n = 1; n <= NMax; ++n)
    {
        threadBarrier();

        float maxDet = 0.0f;
        uint32_t maxDescr = 0;

        for (int iThread = 0; iThread < ThreadCount; ++iThread)
        {
            if (MaxDetA[iThread] > maxDet)
            {
                maxDet = MaxDetA[iThread];
                maxDescr = MaxDescrA[iThread];
            }
        }

        std::cout << "n = " << n << " det = " << maxDet << std::endl;
        buildMat(n, mat, maxDescr);
        printMat(n, mat);

        threadBarrier();
    }

    delete[] mat;

    delete[] MaxDetA;
    delete[] MaxDescrA;

    delete[] threadIdxA;
    delete[] threadA;

    return 0;
}

```
> 3 votes

# Answer
# J

**Update: Improved code to search over half the values. Now calculates `n=12` comfortably within 120 seconds (down from 217s to 60s).**

You will need the latest version of  installed.

```
#!/usr/bin/jconsole

dim =: -:@>:@#
take =: i.@dim
rotstack =: |."0 1~ take
toep =: (dim (|."1 @: {."1) rotstack)"1
det =: -/ . * @: toep
ps =: 3 : ',/(0 1 ,"0 1/ ,.y)'
canonical =: #. >: [: #. |. " 1

lss =: 3 : 0
  shape =. (2^y), y
  shape $ ,>{;/(y,2)$0 1
)

ls =: (canonical@:lss) # lss
ans =: >./ @: det @: ls @: <: @: +:

display =: 3 : 0
echo 'n = ';y;'the answer is';ans y
)
display"0 (1 + i.13)
exit''

```

---

Run this and kill when two minutes are up. My results (MBP 2014 - 16GB of RAM):

```
┌────┬─┬─────────────┬─┐
│n = │1│the answer is│1│
└────┴─┴─────────────┴─┘
┌────┬─┬─────────────┬─┐
│n = │2│the answer is│1│
└────┴─┴─────────────┴─┘
┌────┬─┬─────────────┬─┐
│n = │3│the answer is│2│
└────┴─┴─────────────┴─┘
┌────┬─┬─────────────┬─┐
│n = │4│the answer is│3│
└────┴─┴─────────────┴─┘
┌────┬─┬─────────────┬─┐
│n = │5│the answer is│5│
└────┴─┴─────────────┴─┘
┌────┬─┬─────────────┬─┐
│n = │6│the answer is│9│
└────┴─┴─────────────┴─┘
┌────┬─┬─────────────┬──┐
│n = │7│the answer is│32│
└────┴─┴─────────────┴──┘
┌────┬─┬─────────────┬──┐
│n = │8│the answer is│56│
└────┴─┴─────────────┴──┘
┌────┬─┬─────────────┬───┐
│n = │9│the answer is│125│
└────┴─┴─────────────┴───┘
┌────┬──┬─────────────┬───┐
│n = │10│the answer is│315│
└────┴──┴─────────────┴───┘
┌────┬──┬─────────────┬────┐
│n = │11│the answer is│1458│
└────┴──┴─────────────┴────┘
┌────┬──┬─────────────┬────┐
│n = │12│the answer is│2673│
└────┴──┴─────────────┴────┘

```

Total run time = 61.83s.

---

Just for fun

```
┌────┬──┬─────────────┬────┐
│n = │13│the answer is│8118│
└────┴──┴─────────────┴────┘

```

This took approximately 210 seconds on its own.
> 8 votes

# Answer
# Python 2

This is a very straightforward solution, and probably won't win the contest. But hey, it works!

I'll give a quick overview of what exactly is happening.

1. I first generate every possible starting row for `n`. For example, when `n=2`, this will generate an array length 2n+1, where each row is length 2n-1. It would look like this: `[[0,0,0],[0,0,1],[0,1,0],[0,1,1],[1,0,0],[1,0,1],[1,1,0],[1,1,1]]`.
2. Then, for each of those possible starting rows, I rotate it `n` times and slice off the first `n` items to generate the appropriate matrix, and use `scipy` to calculate the determinant, all while keeping track of the maximum value. At the end of this, I simply print out the maximum, increment `n` by 1, and keep going until 10 minutes has passed.

To run this, you will need scipy installed.

*Edit 1: Changed how initial rows were built by using itertools.product instead, thanks Sp3000!*

*Edit 2: Removed storage of possible starting rows for a minimal improvement in speed.*

*Edit 3: Changed to `scipy` to have more control over how `det` worked.*

```
from scipy import linalg
from collections import deque
from time import time
from itertools import product

c=1
t=time()
while 1:
    m=0
    for d in product(range(2),repeat=2*c-1):
        a=deque(d)
        l=[d[0:c]]
        for x in xrange(c-1):
            a.rotate(1)
            l+=[list(a)[0:c]]
        m=max(m,linalg.det(l,overwrite_a=True,check_finite=False))
    print m,'in',time()-t,'s'
    c+=1

```

Here's some sample output on my home machine (i7-4510U, 8GB RAM):

```
1.0 in 0.0460000038147 s
1.0 in 0.0520000457764 s
2.0 in 0.0579998493195 s
3.0 in 0.0659999847412 s
5.0 in 0.0829999446869 s
9.0 in 0.134999990463 s
32.0 in 0.362999916077 s
56.0 in 1.28399991989 s
125.0 in 5.34999990463 s
315.0 in 27.6089999676 s
1458.0 in 117.513000011 s

```
> 4 votes

# Answer
# C++

Bruteforce with use of OpenMP for parallelization and simple optimization to avoid evaluation of determinant for transposed matrices.

```
$ lscpu
...
Model name:            Intel(R) Core(TM) i5-2410M CPU @ 2.30GHz
...
$ g++ -O2 toepl.cpp -fopenmp
$ timeout 2m ./a.out
1 1
2 1
3 2
4 3
5 5
6 9
7 32
8 56
9 125
10 315
11 1458
12 2673
13 8118
14 22386

```

```
#include <cmath>

#include <algorithm>
#include <iostream>
#include <vector>

using namespace std;

void updateReverses(vector < int > & reverses) {
  int reversesCnt = reverses.size();
  for(int i = 0; i < reversesCnt; ++i){
    reverses[i] <<= 1;
    reverses.push_back(reverses[i] | 1);
  }
}

const double eps = 1e-9;

double determinant(vector < vector < double > > & matrix) {
  int n = matrix.size();
  double det = 1;
  if(n == 1) return matrix[0][0];
  for(int i = 0; i < n; ++i){
    int p = i;
    for(int j = i + 1; j < n; ++j)
      if(fabs(matrix[j][i]) > fabs(matrix[p][i]))
        p = j;
    if(fabs(matrix[p][i]) < eps)
      return 0;
    matrix[i].swap(matrix[p]);
    if(i != p) det *= -1;
    det *= matrix[i][i];
    matrix[i][i] = 1. / matrix[i][i];
    for(int j = i + 1; j < n; ++j)
      matrix[i][j] *= matrix[i][i];
    for(int j = i + 1; j < n; ++j){
      if(fabs(matrix[j][i]) < eps) continue;
      for(int k = i + 1; k < n; ++k)
        matrix[j][k] -= matrix[i][k] * matrix[j][i];
    }
  }
  return det;
}

int main() {
  vector < int > reverses(1, 0);
  reverses.reserve(1 << 30);
  updateReverses(reverses);
  for(int n = 1;; ++n){
    double res = 0;
    int topMask = 1 << (2 * n - 1);
    vector < vector < double > > matrix(n, vector < double > (n));
#pragma omp parallel for reduction(max:res) firstprivate(matrix) schedule(dynamic,1<<10)
    for(int mask = 0; mask < topMask; ++mask){
      if(mask < reverses[mask]) continue;
      for(int i = 0; i < n; ++i)
        for(int j = 0; j < n; ++j)
          matrix[i][j] = (mask >> (i - j + n - 1)) & 1;
      res = max(res, determinant(matrix));
    }
    cout << n << ' ' << res << endl;
    updateReverses(reverses);
    updateReverses(reverses);
  }
}

```
> 4 votes

# Answer
# R

You'll have to install R and the packages listed with `install.packages("package_name")`

Didn't get under 2 mins on my machine with this version (I've to try with a parallel modification)

```
library(pracma)
library(stringr)
library(R.utils)
library(microbenchmark)

f <- function(n) {
  #If n is 1, return 1 to avoid code complexity on this special case
  if(n==1) { return(1) }
  # Generate matrices and get their determinants
  dets <- sapply(strsplit(intToBin( 0:(2^n - 1)), ""), function(x) {
              sapply( strsplit( intToBin( 0:(2^(n-1) - 1) ), ""),
                    function(y) {
                      det(Toeplitz(x,c(x[1],y)))
                    })

              })
  #Get the maximum determinant and return it
  res <- max(abs(dets))
  return(res)
}

```

Call and output:

```
> sapply(1:10,f)
 [1]   1   1   2   3   5   9  32  56 125 315

```

Benchmark on my machine:

```
> microbenchmark(sapply(1:10,f),times=1L)
Unit: seconds
            expr      min       lq     mean   median       uq      max neval
 sapply(1:10, f) 66.35315 66.35315 66.35315 66.35315 66.35315 66.35315     1

```

For information, for a 1:11 range, it takes 285 seconds.
> 2 votes

# Answer
# C

Compile with:

```
$ clang -Ofast 52851.c -o 52851

```

Run with:

```
$ ./52851

```

Can output the maximum determinant for `n = 1..10` in ~115 seconds on my computer.

The program is just getting the determinant every possible binary Toeplitz matrix of size `n`, however every determinant of matrices of size `5x5` or smaller will be cached using memoization.

At first I wrongly assumed that every submatrix of a Toeplitz matrix will also be a Toeplitz matrix, so I only needed to memoize `2^(2n-1)` values instead of `2^(n^2)` for each `n`. I made the program before realizing my mistake, so this submission is just a fix of that program.

---

```
#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <limits.h>
#include <string.h>

#define ELEMENTS(x) (sizeof(x) / sizeof(*x))

int *dets[6];

void print_matrix(int n, int c) {
    for(int row = 0; row < n; row++) {
        for(int col = 0; col < n; col++) {
            int j = n - 1 - row + col;
            int val = !!(c & (1 << j));
            printf("%d ", val);
        }
        puts("");
    }
}

int det(int n, uint8_t *m) {
    if(n == 1) {
        return m[0];
    }

    int i = 0;

    if(n < ELEMENTS(dets)) {
        for(int j = 0; j < n * n; j++) {
            i *= 2;
            i += m[j];
        }

        int v = dets[n][i];
        if(v != INT_MIN) {
            return v;
        }
    }

    int v = 0;

    uint8_t *sub = malloc((n - 1) * (n - 1));

    for(int removed = 0; removed < n; removed++) {
        if(m[removed]) {
            uint8_t *p = sub;
            for(int row = 1; row < n; row++) {
                for(int col = 0; col < n; col++) {
                    if(col == removed) {
                        continue;
                    }

                    *p = m[col + row * n];

                    p++;
                }
            }

            v += (removed % 2 == 0? 1: -1) * det(n - 1, sub);
        }
    }

    free(sub);

    if(n < ELEMENTS(dets)) {
        dets[n][i] = v;
    }
    return v;
}

int main(void) {
    for(int i = 2; i < ELEMENTS(dets); i++) {
        int combinations = 1 << (i * i);
        dets[i] = malloc(combinations * sizeof(**dets));
        for(int j = 0; j < combinations; j++) {
            dets[i][j] = INT_MIN;
        }
    }

    puts("1: 1");

    for(int n = 2; n < 65; n++) {
        int vars = 2 * n - 1;
        size_t combinations = 1 << vars;

        int best = -1;
        int max = -1;

        uint8_t *sub = malloc((n - 1) * (n - 1));

        for(int c = 0; c < combinations; c++) {
            int d = 0;
            for(int i = 0; i < n; i++) {
                if(c & (1 << (n - 1 + i))) {
                    uint8_t *p = sub;
                    for(int row = 1; row < n; row++) {
                        for(int col = 0; col < n; col++) {
                            if(col == i) {
                                continue;
                            }

                            int j = n - 1 - row + col;
                            *p = !!(c & (1 << j));

                            p++;
                        }
                    }
                    d += (i % 2 == 0? 1: -1) * det(n - 1, sub);
                }
            }

            if(d > max) {
                max = d;
                best = c;
            }
        }

        free(sub);

        printf("%d: %d\n", n, max);
        //print_matrix(n, best);
    }

    return 0;
}

```
> 1 votes

# Answer
# PARI/GP, n=11

This is brute force but taking advantage of `det(A^T) = det(A)`. I'm only posting it to demonstrate how easy it is to skip transposes. The lowest bit of `b1` holds the top left cell, and the other bits hold the rest of the top row. `b2` holds the rest of the left column. We simply enforce `b2 <= (b1>>1)`.

```
{ for(n=1,11,
    res=0;
    for(b1=0,2^n-1,
      for(b2=0,b1>>1,
        res=max(res,matdet(matrix(n,n,i,j,bittest(if(i>j,b2>>(i-j-1),b1>>(j-i)),0))));
      )
    );
    print(n" "res);
  )
}

```

Regarding computing Toeplitz determinants in `O(n^2)` time: In my limited research, I've kept running into a requirement that all leading principal minors must be nonzero in order for the algorithms to work, which is a major obstacle for this task. Feel free to give me pointers if you know more about this than I do.
> 1 votes

# Answer
# Rust, \\$n \leq 14 \\$ in 2 min on my laptop

### `src/main.rs`

```
use std::time::{Instant, Duration};
use rulinalg::matrix::Matrix;

fn bit_get(n: u32, i: u32) -> f64 {
    ((n >> i) & 1) as f64
}

fn main() {
    let start = Instant::now();
    let max_duration = Duration::from_secs(120); // 2 minutes

    let mut n = 1;

    while start.elapsed() < max_duration {
        let mut res = 0f64;
        for b1 in 0..(1 << n) {
            for b2 in 0..=(b1 / 2) {
                let values: Vec<f64> = (0..n)
                    .flat_map(|r| (0..n)
                        .map(move |c| if r > c {
                            bit_get(b2, r as u32 - c as u32 - 1)
                        } else {
                            bit_get(b1 as u32, c as u32 - r as u32)
                        }))
                    .collect();
                let matrix = Matrix::new(n as usize, n as usize, values);
                let det = matrix.det();
                res = res.max(det);
            }
        }
        println!("{} {}", n, res);
        n += 1;
    }
}

```

### `Cargo.toml`

```
[package]
name = "rust_hello"
version = "0.1.0"
edition = "2021"

# See more keys and their definitions at https://doc.rust-lang.org/cargo/reference/manifest.html

[dependencies]
ndarray = "0.15"
ndarray-linalg = "0.13"
rulinalg = "0.4"

```

### Build and running

```
# win10 environment

$ set Copt-level=3
$ set Clto=fat
$ set Ctarget-cpu=native
$ cargo build --release
$ target\release\rust_hello.exe
1 1
2 1
3 2
4 3
5 5
6 9
7 32
8 56
9 125
10 315
11 1458
12 2673
13 8117.999999999997  # roughly 8118
14 22386.000000000007 # roughly 22386

```
> 1 votes

# Answer
# Mathematica

So slow. It can be optimized more.

Try it online!

```
Clear["Global`*"];
Table[res = 0;
 For[b1 = 0, b1 <= 2^n - 1, b1++,
  For[b2 = 0, b2 <= Floor[b1/2], b2++,
   matrix =
    Array[If[#1 > #2, BitGet[b2, #1 - #2 - 1],
       BitGet[b1, #2 - #1]] &, {n, n}, 0];
   res = Max[res, Det[matrix]];]];
 Print[n, " ", res], {n, 1, 11}]

```

```
1,1,2,3,5,9,32,56,125,315,1458

```
> 0 votes

---
Tags: math, fastest-code
---


title: Carmichael Numbers

text:
# Question
Title: Carmichael Numbers
Write a program to print Carmichael numbers. The challenge is to print the maximum number of them in less than 5 secs. The numbers should be in increasing order. Please do not use built in functions. Also please don't store the carmichael numbers in any form. To keep a common time measurement use Ideone. The program to print the maximum number of Carmichael number wins.

```
Sample Output
561
1105
1729
2465
2821
6601
8911
.
.
.

```

**OEIS A002997**

# Answer
## Python

This generates 1332 Carmichael numbers in under 5 seconds. Remove the print to actually get ideone to run to completion.

```
N=180000

maxprime = 18*N+1
prime = [1]*maxprime
prime[0] = prime[1] = 0
for i in xrange(maxprime):
  if prime[i]:
    for j in xrange(2*i,maxprime,i):prime[j]=0

for i in xrange(N):
  a=6*i+1
  b=12*i+1
  c=18*i+1
  if prime[a] and prime[b] and prime[c]:print a*b*c

```

Note that you could make this a whole lot faster in C.
> 2 votes

# Answer
**unGolf3D dear Java :)** maxCount=1193221 (at my local) and 552721 (at ideone) in 4.8s

```
    public class Main{
    public static void main(String[] args){
          long start=System.currentTimeMillis();
          int N=1200000;
          int list[]=new int[N];
            for (int i=0;i<N;i++) {
                list[i]=(i+1);
            }
            for (int i=1;i<N;i++) {
                if (Math.sqrt(list[i])-Math.floor(Math.sqrt(list[i]))==0) {
                    for (int j=i;j<N;j+=(i+1)){
                        list[j]=-1;
                    }
                }
            }
        boolean Carmichael=false;
        int chemCount=0;
        int product=1;
          for(int i=561  ;(System.currentTimeMillis()-start)/1000.0<=5;i++){
              if(!isPrime(i)&&list[i-1]==i){
                  for(int div=2;div<Math.sqrt(i);div++){
                      if(i%div==0 && isPrime(div)){
                          if((i-1)%(div-1)==0){
                              Carmichael= true;
                              product*=div;
                              chemCount++;
                              if (product==i) break;
                              }
                          else{
                              Carmichael=false;
                              break;
                              }
                      }
                  }
                  if (Carmichael && chemCount>=3 && product==i)
                      System.out.println(i);
                  Carmichael=false;
                   chemCount=0;
                   product=1;
              }
          }
  }

        static boolean isPrime(long number){
            if (number==2)
                return true;
            if (number==1)
                return false;
            for (long i=2;i<=Math.sqrt(number);i++) {
                if (number%i==0)
                    return false;
            }
            return true;
        }
    }

```

My JVM gets this to 1193221 number in 5 sec... I am sure this is slower than any other scripting languages and algo!! :( ... OTOH, http://www.ideone.com/clone/oNJin fails to create this big process.. calculates upto 552721 in 4.8 sec and breaks after that throwing http://en.wikipedia.org/wiki/SIGXCPU :)
> 1 votes

# Answer
## C++

Machine i5-760 @ 2.8Ghz, only one core used

Numbers found in 5 seconds: 129

```
#define WIN32_MEAN_AND_LEAN
#include <windows.h>
#include <iostream>
using namespace std;

bool IsPrime (int v)
{
  bool
    is_prime = true;

  for (int i = 3, i2 = 9 ; i2 < v ; ++i, i2 += i + i + 1)
  {
    if (v % i == 0)
    {
      is_prime = false;
      break;
    }
  }

  return is_prime;
}

int Prime (int i, int &i2)
{
  static int
    primes [10000] = {2, 3, 5},
    primes2 [10000] = {4, 9, 25},
    num_primes = 3;

  while (i >= num_primes)
  {
    int
      v = primes [num_primes - 1];

    do
    {
      v += 2;
    } while (!IsPrime (v));

    primes [num_primes] = v;
    primes2 [num_primes] = v * v;
    ++num_primes;
  }

  i2 = primes2 [i];
  return primes [i];
}

bool IsCarmichael (int n)
{
   if (n < 2)
   {
      return false;
   }

   int
     k = n--,
     i2 = 4;

   for (int i = 1, p = 2 ; i2 <= k ; p = Prime (i++, i2))
   {
      if (k % p == 0)
      {
         k /= p;

         if (k % p == 0 || n % (p - 1) != 0)
         {
            return false;
         }

         i = 0;
      }
   }

   return --k != n && n % k == 0;
}

int main ()
{
  DWORD
    end = GetTickCount () + 5000;

  int
    count = 0,
    n = 561;

  while (GetTickCount () < end)
  {
    for (int i = 0 ; i < 1000 ; ++i)
    {
      if (IsCarmichael (n))
      {
        ++count;
        cout << n << endl;
      }
      ++n;
    }
  }

  cout << count << " Carmichael numbers found out of " << (n - 561) << " values tested." << endl;
}

```
> 1 votes

# Answer
# D Language, 139 numbers, 18900973 max

Implemented Korselt's criterion to extract Carmichael numbers, considering only odd numbers greater than 3 since an even number cannot satisfy the relation.

## Code

```
import std.stdio, std.math;

bool isCarmichael(uint n) pure nothrow @nogc {
   auto maxValue = ceil(sqrt(float(n)));
   uint cn = n;

   for (uint i = 3; i < maxValue; i += 2) {
      if (cn % i == 0) {
         cn /= i;

         if (cn % i == 0)
            return false;

         if ((n-1) % (i-1))
            return false;
      }
   }

   if (cn > 1)
      return false;

   return true;
}

void main() {
   for (uint n = 3; ; n += 2) {
      if (isCarmichael(n))
         writeln(n);
   }
}

```

## Output

```
561
1105
1729
2465
2821
6601
8911
10585
15841
29341
41041
46657
52633
62745
63973
75361
101101
115921
126217
162401
172081
188461
252601
278545
294409
314821
334153
340561
399001
410041
449065
488881
512461
530881
552721
656601
658801
670033
748657
825265
838201
852841
997633
1024651
1033669
1050985
1082809
1152271
1193221
1461241
1569457
1615681
1773289
1857241
1909001
2100901
2113921
2433601
2455921
2508013
2531845
2628073
2704801
3057601
3146221
3224065
3581761
3664585
3828001
4335241
4463641
4767841
4903921
4909177
5031181
5049001
5148001
5310721
5444489
5481451
5632705
5968873
6049681
6054985
6189121
6313681
6733693
6840001
6868261
7207201
7519441
7995169
8134561
8341201
8355841
8719309
8719921
8830801
8927101
9439201
9494101
9582145
9585541
9613297
9890881
10024561
10267951
10402561
10606681
10837321
10877581
11119105
11205601
11921001
11972017
12261061
12262321
12490201
12945745
13187665
13696033
13992265
14469841
14676481
14913991
15247621
15403285
15829633
15888313
16046641
16778881
17098369
17236801
17316001
17586361
17812081
18162001
18307381
18900973

```

The code runs on a single core, without parallelism. No luck making it work on IDEONE though; time is checked on my machine running the process with `timeout 5s` UNIX command.
> 0 votes

# Answer
# Rust

Naive brute force. It can be optimized more.

Run it on Rust Playground!

```
extern crate rand;

fn gcd(a: i128, b: i128) -> i128 {
    if b == 0 {
        a
    } else {
        gcd(b, a % b)
    }
}

fn power(mut a: i128, mut b: i128, m: i128) -> i128 {
    let mut result = 1;
    a %= m;
    while b > 0 {
        if b % 2 != 0 {
            result = (result * a) % m;
        }
        b >>= 1;
        a = (a * a) % m;
    }
    result
}

fn is_prime_miller_rabin(n: i128) -> bool {
    if n < 2 {
        return false;
    }
    if n != 2 && n % 2 == 0 {
        return false;
    }
    let s = ((n - 1) as u128).trailing_zeros() as i128;
    let d = (n - 1) / 2i128.pow(s as u32);

    let witnesses: [i128; 12] = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37];
    'next_witness: for &a in witnesses.iter().take_while(|&&x| x <= n - 1) {
        let mut x = power(a, d, n);
        if x == 1 || x == n - 1 {
            continue;
        }
        for _r in 1..s {
            x = power(x, 2, n);
            if x == n - 1 {
                continue 'next_witness;
            }
        }
        return false;
    }
    true
}

fn fermat_test(n: i128) -> bool {
    for a in 2..n {
        if gcd(a, n) == 1 {
            if power(a, n - 1, n) != 1 {
                return false;
            }
        }
    }
    true
}

fn is_carmichael(n: i128) -> bool {
    !is_prime_miller_rabin(n) && fermat_test(n)
}

fn main() {
    let start = 561;
    let end = 100000;  // Adjust end limit based on time constraint
    for n in start..end {
        if is_carmichael(n) {
            println!("{}", n);
        }
    }
}

```
> 0 votes

# Answer
# Python 3, 1184 numbers

Somewhat simplified implementation of Richard Pinch's methods (Carmichael numbers up to 10^21, Carmichael numbers up to 10^16).

The idea is to enumerate all the Carmichael numbers less than \\$N\\$. We can prove that none of them will have prime factors \\$\>\sqrt{N}\\$, so we split into two cases based on the size of the largest prime factor (specifically if it is larger or smaller than \\$N^{0.3}\\$).

We'll let \\$n\\$ be the number, and \\$q\\$ it's largest factor. If there is a factor greater than \\$N^{0.3}\\$, then since \\$\frac{n}{q}\equiv1\mod{q-1}\\$, enumerating all possible values of \\$\frac{n}{q}\\$ given \\$q\\$ is \\$\mathcal{O}\left(N^{0.4}\right)\\$. So enumerate \\$q\\$ from \\$N^{0.3}\\$ to \\$N^{0.5}\\$, and test with a Fermat test to base 2 and Korselt's criterion.

In the other case, we use a depth first search on products of prime factors, where we can limit the branching factor to \\$N^{0.3}\\$.

Try it on ideone!

```
import math
from collections import deque
import itertools
# from gmpy2 import powmod as pow

def main():
    x = find_carmichael(5_000_000_000)
    for idx, i in enumerate(sorted(x), 1):# + len(s)):
        print(idx, i)

def find_carmichael(X):
    if X = q: return False
        x = i
        f.append(i)
    if len(f) >= 2:
        n1 = P*q-1
        return all(n1%(p-1)==0 for p in f)

def dfs(X, primes, n=1, t=(), i=0, lcm=1):
    if math.gcd(n, lcm) > 1: return
    if lcm * n * len(primes) > X:
        for Q in range(pow(n, -1, lcm), X//n + 1, lcm):
            if math.gcd(n, Q) > 1: continue
            N = n*Q
            if pow(2, N, Q) == 2:
                if is_carmichael0(n, Q, t):
                    yield N
        return
    for i1 in range(i, len(primes)):
        p = primes[i1]
        n1 = n*p
        if n1 > X: break
        if len(t) >= 2 and (n-1)%(p-1) == 0 and n1%lcm == 1:
            yield n1
        yield from dfs(X, primes, n1, t+(p,), i1+1, math.lcm(lcm, p-1))

def is_carmichael0(n, Q, t):
    x = 0
    f = list(t)
    for i in factor(Q):
        if i == x: return False
        if i in t: return False
        x = i
        f.append(i)
    if len(f) >= 3:
        n1 = n*Q-1
        return all(n1%(p-1)==0 for p in f)

class LazyList:
    def __init__(self, I):
        self.I = I
        self.l = []
    def __getitem__(self, i):
        try:
            return self.l[i]
        except IndexError:
            self.l.extend(itertools.islice(self.I, i+1 - len(self.l)))
            return self.l[i]
    def range(self, a, b=None):
        if b is None:
            a, b = 0, a
        for i in self:
            if i = b: break
            yield i

def primes_():
    yield 2
    l = [3,5,7,11,13,17,19,23]
    yield from l
    p2s = deque([(p, p*p) for p in l[2:]])
    n = 23
    stuff = [(0, 3), (4, 5)]
    while True:
        p2, next_p2 = p2s.popleft()
        R = range(n+4, next_p2, 2)
        z = [[] for _ in R]
        for i, p in stuff:
            z[i].append(p)
        stuff.clear()
        stuff.append((p2-1, p2))
        for i, (l, n) in enumerate(zip(z, R)):
            # print(i, l, n)
            if l:
                for p in l:
                    i1 = i + p
                    if i1 >= len(R):
                        stuff.append((i1 - len(R) - 1, p))
                    else:
                        z[i1].append(p)
            else:
                yield n
                p2s.append((n, n*n))

primes = LazyList(primes_())

def factor(n):
    for p in primes:
        while n%p == 0:
            yield p
            n //= p
        if p*p > n:
            if n > 1: yield n
            return

if __name__ == '__main__': main()

```
> 0 votes

---
Tags: sequence, fastest-code
---