text:
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.

Define as quickly as possible.

@LegionMammal978 With the best runtime performance. It's a fastest-code challenge.

http://oeis.org/A037274

How do we know which code is faster? Some people may be testing on a five-year-old laptop (cough like me cough), while others may be using some high-end desktop/server. Also, performance varies among interpreters of the same language.

@JungHwanMin As per standard fastest-code rules, the host must test each program and time it on their computer. I suggest you use the fastest interpreter, and tell me which one you use.

@JungHwanMin all code is tested on the same system, so that's irrelevant

@JungHwanMin the host as in the OP

@TrojanByAccident I stand corrected... I misread that.

Is using a probabilistic primality test such as Miller-Rabin allowed?

@miles Yes, but only to boost the speed. If the program spits out a wrong answer, the program is disqualified.

Please post the correct answers for the test cases.

Just curious, is it known that this will always terminate for a given input?

@NoahL Are you going to test the answers anytime soon? You said "Program speed will be tested on a Raspberry Pi 2"

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:

n < 230 → factorization by trial division
n < 270 → Pollard's Rho (Brent's variant)
else → Lenstra ECF (using Montgomery curves, Suyama's parameters), test ~100 curves
if ECF failed → Multiple Polynomial Quadratic Sieve (using Montgomery's polynomials)

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

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.

The 1 at the end shouldn't be there...

I don't have Mathematica on my computer, meaning I'll have to test this (and the rest of the programs) on my Raspberry Pi 2.

Since we're not golfing: there's CompositeQ for !PrimeQ (which also ensures that your answer doesn't loop forever on input 1).

How is it possible that Mathematica performs HP(80) in so short time without having the primes hardcoded somewhere? My i7 laptop is taking hours to perform a primality check, as well as finding the prime factors, for HP(80) when it arrives at 47109211289720051.

@NoahL Mathematica can be tested online. http://meta.codegolf.stackexchange.com/a/1445/34718

@Mario you're allowed to use a list of primes, just not a list of home primes.

String manipulation is about 10% faster NestWhile[ToExpression[StringJoin[StringRepeat[ToString[#[[1]]],#[[2]]]&/@FactorInteger[#]]]&,#,CompositeQ]

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

Making an array with pre-determined primes/values is not allowed, I believe, since its a standard loophole.

@P.Ktinos I think so too... anyway it would be too big to include.

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.

pip install primefac worked for me, although 65 and 80 don't seem to run on windows, due to forking in the background.

Looking at the source for primefac was pretty funny, as there are a lot of comments with TODO or find out why this is throwing errors

I just did too. The author actually uses my mpqs! ...slightly modified. Line 551 # This occasionally throws IndexErrors. Yeah, because he removed the check that there are more smooths than primes factors used.

You should help him out. :)

I'll probably contact him when I'm done with this challenge, I intend to work on optimizing mpqs a bit (gotta beat mathmatica, am I right?).

@primo If only it were code-golf. I could win. :D

I can't figure out how to contact the guy. :(

@primo lucasbrown.cit@gmail.com Found here

Didn't see that it was a email link, thanks ;)

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.

I found that any factor found by P-1 was also found - sooner - by ECM, so I dropped it (same holds for Williams P+1). Maybe I'll try adding SQUFOF. Brilliant library, btw.

Also, use feature "say";.

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");
)
}

source: Stack Exchange

license: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/3.0/

all_licenses: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/4.0/

sort: votes


text:
Find the shortest Golomb rulers
Golomb rulers are sets of non-negative integers such that no two pairs of integers in the set are the same distance apart.
For example, [0, 1, 4, 6] is a Golomb ruler because all distances between two integers in this set are unique:
0, 1 -> distance 1
0, 4 -> distance 4
0, 6 -> distance 6
1, 4 -> distance 3
1, 6 -> distance 5
4, 6 -> distance 2

For the sake of simplicity in this challenge (and since translation is trivial), we impose that a Golomb ruler always contains the number 0 (which the previous example does).
Since this set is of length 4, we say that this is a Golomb ruler of order 4. The biggest distance in this set (or element, since 0 is always in the set) is 6, therefore we say that this is a Golomb Ruler of length 6.
Your task
Find Golomb rulers of order 50 to 100 (inclusive) that have as small lengths as you can find. The rulers you find need not be optimal (see below).
Optimality
A Golomb ruler of order N, is said to be optimal if there is no other Golomb ruler of order N which has a smaller length.
Optimal Golomb rulers are known for orders less than 28, though finding and proving optimality is harder and harder as the order increases.
Therefore, it is not expected that you find the optimal Golomb ruler for any of the orders between 50 and 100 (and even less expected that you can prove they are optimal).
There are no time limits in the execution of your program.
Baseline
The list below is the list of lengths of Golomb rulers from 50 to 100 (in order) evaluated with a naïve search strategy (Thanks to @PeterTaylor for this list):
[4850 5122 5242 5297 5750 5997 6373 6800 6924 7459 7546 7788 8219 8502 8729 8941 9881 10199 10586 10897 11288 11613 11875 12033 12930 13393 14046 14533 14900 15165 15687 15971 16618 17354 17931 18844 19070 19630 19669 20721 21947 22525 23290 23563 23880 24595 24767 25630 26036 26254 27218]

The sum of all those lengths is 734078.
Scoring
Your score will be the sum of the lengths of all your Golomb rulers between 50 and 100, divided by the sum of the lengths of Golomb rulers between 50 and 100 in the baseline: 734078.
In case you did not find a Golomb ruler for a specific order, you shall compute your score the same way, using the double of the length in the baseline for that specific order.
The answer with the lowest score wins.
In case of a tie, the lengths of the biggest order where the two answers differ are compared, and the shortest one wins. In case both answers have the same lengths for all orders, then the answer that was posted first wins.

Related. (Same challenge in 2D.)

And OEIS entry.

When you say rulers between 50 and 100, do you mean the range [50, 100)? So not including the order 100 ruler? Because the baseline only contains 50 elements.

@orlp The baseline length for 50 was missing. I have added it and fixed the total baseline length. So this is in the range [50, 100]. Thanks.

Side note: the smallest possible length of a Golomb ruler of order n is n(n-1)/2, since that's how many positive differences there are. Therefore the smallest possible score in this challenge is 147050/734078 > 0.2003193.

@GregMartin Thanks, though this is not quite the "smallest possible score" but rather a lower bound on that smallest possible score!

C#, 259421/734078 ~= 0.3534
Methods
I finally found a more-or-less readable explanation of the projective field method (Singer's method) in Constructions of Generalised Sidon Sets, although I still think it can be improved slightly. It turns out to be more similar to the affine field method (Bose's method) than the other papers I read had communicated.
This assumes knowledge of finite fields. Consider \$q = p^a\$ is a prime power, and let \$F(q)\$ be our base field.
The affine field method works over \$F(q^2)\$. Take a generator \$g_2\$ of \$F(q^2)\$ and a non-zero element \$k\$ of \$F(q)\$ and consider the set $$\{a : g_2{}^a - k g_2 \in F_q\}$$ Those values form a modular Golomb ruler mod \$q^2 - 1\$. Further rulers can be obtained by multiplying modulo \$q^2 - 1\$ by any number which is coprime with the modulus.
The projective field method works over \$F(q^3)\$. Take a generator \$g_3\$ of \$F(q^3)\$ and a non-zero element \$k\$ of \$F(q)\$ and consider the set $$\{0\} \cup \{a : g_3{}^a - k g_3 \in F_q\}$$ Those values form a modular Golomb ruler mod \$q^2 + q + 1\$. Further rulers can be obtained by modular multiplication in the same way as for the affine field method.
Note that these methods between them give the best known values for every length greater than 16. Tomas Rokicki and Gil Dogon are offering a $250 reward for anyone who beats them for lengths 36 to 40000. Therefore anyone who beats this answer is in for a monetary prize.
Code
The C# isn't very idiomatic, but I need it to compile with an old version of Mono. Also, despite the argument checking, this is not production quality code. I'm not happy with the types, but I don't think there's a really good solution to that in C#. Maybe in F# or C++ with insane templating.
using System;
using System.Collections.Generic;
using System.Linq;

namespace Sandbox {
    class Program {
        static void Main(string[] args) {
            var winners = ComputeRulerRange(50, 100);
            int total = 0;
            for (int i = 50; i <= 100; i++) {
                Console.WriteLine("{0}:\t{1}", i, winners[i][i - 1]);
                total += winners[i][i - 1];
            }
            Console.WriteLine("\t{0}", total);
        }

        static IDictionary<int, int[]> ComputeRulerRange(int min, int max) {
            var best = new Dictionary<int, int[]>();

            var naive = Naive(max);
            for (int i = min; i <= max; i++) best[i] = naive.Take(i).ToArray();

            var finiteFields = FiniteFields(max * 11 / 10).OrderBy(x => x.Size).ToArray();

            // The projective plane method generates rulers of length p^a + 1 for prime powers p^a.
            // We can then look at subrulers for a reasonable range, say down to two prime powers below.
            for (int ppi = 0; ppi < finiteFields.Length; ppi++) {
                // Range under consideration
                var field = finiteFields[ppi];
                int q = field.Size;
                int subFrom = Math.Max(min, ppi >= 2 ? finiteFields[ppi - 2].Size : 1);
                int subTo = Math.Min(max, q + 1);
                if (subTo < subFrom) continue;

                int m = q * q + q + 1;
                foreach (var ruler in ProjectiveRulers(field)) {
                    for (int sub = subFrom; sub <= subTo; sub++) {
                        var subruler = BestSubruler(ruler, sub, m);
                        if (subruler[sub - 1] < best[sub][sub - 1]) best[sub] = subruler;
                    }
                }
            }

            // Similarly for the affine plane method, which generates rulers of length p^a for prime powers p^a
            for (int ppi = 0; ppi < finiteFields.Length; ppi++) {
                // Range under consideration
                var field = finiteFields[ppi];
                int q = field.Size;
                int subFrom = Math.Max(min, ppi >= 2 ? finiteFields[ppi - 2].Size : 1);
                int subTo = Math.Min(max, q);
                if (subTo < subFrom) continue;

                int m = q * q - 1;
                foreach (var ruler in AffineRulers(field)) {
                    for (int sub = subFrom; sub <= subTo; sub++) {
                        var subruler = BestSubruler(ruler, sub, m);
                        if (subruler[sub - 1] < best[sub][sub - 1]) best[sub] = subruler;
                    }
                }
            }

            return best;
        }

        static int[] BestSubruler(int[] ruler, int sub, int m) {
            int[] expand = new int[ruler.Length + sub - 1];
            for (int i = 0; i < ruler.Length; i++) expand[i] = ruler[i];
            for (int i = 0; i < sub - 1; i++) expand[ruler.Length + i] = ruler[i] + m;

            int best = m, bestIdx = -1;
            for (int i = 0; i < ruler.Length; i++) {
                if (expand[i + sub - 1] - expand[i] < best) {
                    best = expand[i + sub - 1] - expand[i];
                    bestIdx = i;
                }
            }

            return expand.Skip(bestIdx).Take(sub).Select(x => x - ruler[bestIdx]).ToArray();
        }

        static IEnumerable<int[]> ProjectiveRulers(FiniteField field) {
            var q = field.Size;
            var fq3 = PowerField.Create(field, 3);
            var m = q * q + q + 1;
            var g = fq3.Generators.First();

            // Define the set T<k> = {0} \union {a \in [q^3-1] : g^a - kg \in F(q)} for 0 != k \in F(q)
            // This could alternatively be T<k> = {0} \union {log_g(b - kg) : b in F(q)} for 0 != k \in F(q)
            // Then T<k> % (q^2 + q + 1) gives a Golomb ruler.
            // For a given generator we seem to get the same ruler for every k.
            var t_k = new HashSet<int>();
            t_k.Add(0);
            var ga = fq3.One;
            for (int a = 1; a < fq3.Size; a++) {
                ga = ga * g;
                if (fq3.Convert(ga + g) < q) t_k.Add(a % m);
            }

            // TODO: optimise by detecting duplicates
            for (int s = 1; s < m; s++) {
                if (Gcd(s, m) == 1) yield return t_k.Select(x => x * s % m).OrderBy(x => x).ToArray();
            }
        }

        static IEnumerable<int[]> AffineRulers(FiniteField field) {
            var q = field.Size;
            var fq2 = PowerField.Create(field, 2);
            var m = q * q - 1;
            var g = fq2.Generators.First();

            // Define the set T<k> = {0} \union {a \in [q^2-1] : g^a - kg \in F(q)} for 0 != k \in F(q)
            // Then T<k> % (q^2 - 1) gives a Golomb ruler.
            var t_k = new HashSet<int>();
            var ga = fq2.One;
            for (int a = 1; a < fq2.Size; a++) {
                ga = ga * g;
                if (fq2.Convert(ga + g) < q) t_k.Add(a % m);
            }

            // TODO: optimise by detecting duplicates
            for (int s = 1; s < m; s++) {
                if (Gcd(s, m) == 1) yield return t_k.Select(x => x * s % m).OrderBy(x => x).ToArray();
            }
        }

        static int Gcd(int a, int b) {
            while (a != 0) {
                var t = b % a;
                b = a;
                a = t;
            }

            return b;
        }

        static int[] Naive(int size) {
            if (size == 0) return new int[0];
            if (size == 1) return new int[] { 0 };

            int[] ruler = new int[size];
            var diffs = new HashSet<int>();
            int i = 1, c = 1;
            while (true) {
                bool valid = true;
                for (int j = 0; j < i; j++) {
                    if (diffs.Contains(c - ruler[j])) { valid = false; break; }
                }

                if (valid) {
                    for (int j = 0; j < i; j++) diffs.Add(c - ruler[j]);
                    ruler[i++] = c;
                    if (i == size) return ruler;
                }

                c++;
            }
        }

        static IEnumerable<FiniteField> FiniteFields(int max) {
            bool[] isComposite = new bool[max + 1];
            for (int p = 2; p < isComposite.Length; p++) {
                if (!isComposite[p]) {
                     FiniteField baseField = new PrimeField(p); yield return baseField;
                    for (int pp = p * p, pow = 2; pp < max; pp *= p, pow++) yield return PowerField.Create(baseField, pow);
                    for (int pq = p * p; pq <= max; pq += p) isComposite[pq] = true;
                }
            }
        }
    }

    public abstract class FiniteField {
        private Lazy<FiniteFieldElement> _Zero;
        private Lazy<FiniteFieldElement> _One;

        public FiniteFieldElement Zero { get { return _Zero.Value; } }
        public FiniteFieldElement One { get { return _One.Value; } }
        public IEnumerable<FiniteFieldElement> Generators {
            get {
                for (int _g = 1; _g < Size; _g++) {
                    int pow = 0;
                    FiniteFieldElement g = Convert(_g), gpow = One;
                    while (true) {
                        pow++;
                        gpow = gpow * g;
                        if (gpow == One) break;
                        if (pow > Size) {
                            throw new Exception("Is this really a field? " + this);
                        }
                    }
                    if (pow == Size - 1) yield return g;
                }
            }
        }

        public abstract int Size { get; }
        internal abstract FiniteFieldElement Convert(int i);
        internal abstract int Convert(FiniteFieldElement f);

        internal abstract bool Eq(FiniteFieldElement a, FiniteFieldElement b);
        internal abstract FiniteFieldElement Negate(FiniteFieldElement a);
        internal abstract FiniteFieldElement Add(FiniteFieldElement a, FiniteFieldElement b);
        internal abstract FiniteFieldElement Mul(FiniteFieldElement a, FiniteFieldElement b);

        protected FiniteField() {
            _Zero = new Lazy<FiniteFieldElement>(() => Convert(0));
            _One = new Lazy<FiniteFieldElement>(() => Convert(1));
        }
    }

    public abstract class FiniteFieldElement {
        internal abstract FiniteField Field { get; }

        public static FiniteFieldElement operator -(FiniteFieldElement a) {
            return a.Field.Negate(a);
        }

        public static FiniteFieldElement operator +(FiniteFieldElement a, FiniteFieldElement b) {
            if (a.Field != b.Field) throw new ArgumentOutOfRangeException("b");
            return a.Field.Add(a, b);
        }

        public static FiniteFieldElement operator *(FiniteFieldElement a, FiniteFieldElement b) {
            if (a.Field != b.Field) throw new ArgumentOutOfRangeException("b");
            return a.Field.Mul(a, b);
        }

        public static bool operator ==(FiniteFieldElement a, FiniteFieldElement b) {
            if (Equals(a, null)) return Equals(b, null);
            else if (Equals(b, null)) return false;

            if (a.Field != b.Field) throw new ArgumentOutOfRangeException("b");
            return a.Field.Eq(a, b);
        }

        public static bool operator !=(FiniteFieldElement a, FiniteFieldElement b) { return !(a == b); }

        public override bool Equals(object obj) {
            return (obj is FiniteFieldElement) && (obj as FiniteFieldElement).Field == Field && this == (obj as FiniteFieldElement);
        }

        public override int GetHashCode() { return Field.Convert(this).GetHashCode(); }

        public override string ToString() { return Field.Convert(this).ToString(); }
    }

    public class PrimeField : FiniteField {
        private readonly int _Prime;
        private readonly PrimeFieldElement[] _Featherweight;

        internal int Prime { get { return _Prime; } }
        public override int Size { get { return _Prime; } }

        public PrimeField(int prime) {
            if (prime < 2) throw new ArgumentOutOfRangeException("prime");

            // TODO A primality test would be nice...

            _Prime = prime;
            _Featherweight = new PrimeFieldElement[Math.Min(prime, 256)];
        }

        internal override FiniteFieldElement Convert(int i) {
            if (i < 0 || i >= _Prime) throw new ArgumentOutOfRangeException("i");
            if (i >= _Featherweight.Length) return new PrimeFieldElement(this, i);
            if (Equals(_Featherweight[i], null)) _Featherweight[i] = new PrimeFieldElement(this, i);
            return _Featherweight[i];
        }

        internal override int Convert(FiniteFieldElement f) {
            if (f == null) throw new ArgumentNullException("f");
            if (f.Field != this) throw new ArgumentOutOfRangeException("f");

            return (f as PrimeFieldElement).Value;
        }

        internal override bool Eq(FiniteFieldElement a, FiniteFieldElement b) {
            if (a.Field != this) throw new ArgumentOutOfRangeException("a");
            if (b.Field != this) throw new ArgumentOutOfRangeException("b");

            return (a as PrimeFieldElement).Value == (b as PrimeFieldElement).Value;
        }

        internal override FiniteFieldElement Negate(FiniteFieldElement a) {
            if (a.Field != this) throw new ArgumentOutOfRangeException("a");
            var fa = a as PrimeFieldElement;
            return fa.Value == 0 ? fa : Convert(_Prime - fa.Value);
        }

        internal override FiniteFieldElement Add(FiniteFieldElement a, FiniteFieldElement b) {
            if (a.Field != this) throw new ArgumentOutOfRangeException("a");
            if (b.Field != this) throw new ArgumentOutOfRangeException("b");

            return Convert(((a as PrimeFieldElement).Value + (b as PrimeFieldElement).Value) % _Prime);
        }

        internal override FiniteFieldElement Mul(FiniteFieldElement a, FiniteFieldElement b) {
            if (a.Field != this) throw new ArgumentOutOfRangeException("a");
            if (b.Field != this) throw new ArgumentOutOfRangeException("b");

            return Convert(((a as PrimeFieldElement).Value * (b as PrimeFieldElement).Value) % _Prime);
        }

        public override string ToString() { return string.Format("F({0})", _Prime); }
    }

    internal class PrimeFieldElement : FiniteFieldElement {
        private readonly PrimeField _Field;
        private readonly int _Value;

        internal override FiniteField Field { get { return _Field; } }
        internal int Value { get { return _Value; } }

        internal PrimeFieldElement(PrimeField field, int val) {
            if (field == null) throw new ArgumentNullException("field");
            if (val < 0 || val >= field.Prime) throw new ArgumentOutOfRangeException("val");

            _Field = field;
            _Value = val;
        }
    }

    public class PowerField : FiniteField {
        private readonly FiniteField _BaseField;
        private readonly FiniteFieldElement[] _Polynomial;

        internal FiniteField BaseField { get { return _BaseField; } }
        internal int Power { get { return _Polynomial.Length; } }
        public override int Size { get { return (int)Math.Pow(_BaseField.Size, Power); } }

        public PowerField(FiniteField baseField, FiniteFieldElement[] polynomial) {
            if (baseField == null) throw new ArgumentNullException("baseField");
            if (polynomial == null) throw new ArgumentNullException("polynomial");
            if (polynomial.Length < 2) throw new ArgumentOutOfRangeException("polynomial");
            for (int i = 0; i < polynomial.Length; i++) if (polynomial[i].Field != baseField) throw new ArgumentOutOfRangeException("polynomial[" + i + "]");

            // TODO Check that the polynomial is irreducible over the base field.

            _BaseField = baseField;
            _Polynomial = polynomial.ToArray();
        }

        internal override FiniteFieldElement Convert(int i) {
            if (i < 0 || i >= Size) throw new ArgumentOutOfRangeException("i");

            var vec = new FiniteFieldElement[Power];
            for (int j = 0; j < vec.Length; j++) {
                vec[j] = BaseField.Convert(i % BaseField.Size);
                i /= BaseField.Size;
            }

            return new PowerFieldElement(this, vec);
        }

        internal override int Convert(FiniteFieldElement f) {
            if (f == null) throw new ArgumentNullException("f");
            if (f.Field != this) throw new ArgumentOutOfRangeException("f");

            var pf = f as PowerFieldElement;
            int i = 0;
            for (int j = Power - 1; j >= 0; j--) i = i * BaseField.Size + BaseField.Convert(pf.Value[j]);
            return i;
        }

        internal override bool Eq(FiniteFieldElement a, FiniteFieldElement b) {
            if (a.Field != this) throw new ArgumentOutOfRangeException("a");
            if (b.Field != this) throw new ArgumentOutOfRangeException("b");

            var fa = a as PowerFieldElement;
            var fb = b as PowerFieldElement;
            for (int i = 0; i < Power; i++) if (fa.Value[i] != fb.Value[i]) return false;
            return true;
        }

        internal override FiniteFieldElement Negate(FiniteFieldElement a) {
            if (a.Field != this) throw new ArgumentOutOfRangeException("a");
            return new PowerFieldElement(this, (a as PowerFieldElement).Value.Select(x => -x).ToArray());
        }

        internal override FiniteFieldElement Add(FiniteFieldElement a, FiniteFieldElement b) {
            if (a.Field != this) throw new ArgumentOutOfRangeException("a");
            if (b.Field != this) throw new ArgumentOutOfRangeException("b");

            var fa = a as PowerFieldElement;
            var fb = b as PowerFieldElement;
            var vec = new FiniteFieldElement[Power];
            for (int i = 0; i < Power; i++) vec[i] = fa.Value[i] + fb.Value[i];
            return new PowerFieldElement(this, vec);
        }

        internal override FiniteFieldElement Mul(FiniteFieldElement a, FiniteFieldElement b) {
            if (a.Field != this) throw new ArgumentOutOfRangeException("a");
            if (b.Field != this) throw new ArgumentOutOfRangeException("b");

            var fa = a as PowerFieldElement;
            var fb = b as PowerFieldElement;

            // We consider fa and fb as polynomials of a variable x and multiply modulo (x^Power - _Polynomial).
            // But to keep things simple we want to manage the cascading modulo.
            var vec = Enumerable.Repeat(BaseField.Zero, Power).ToArray();
            var fa_xi = fa.Value.ToArray();
            for (int i = 0; i < Power; i++) {
                for (int j = 0; j < Power; j++) vec[j] += fb.Value[i] * fa_xi[j];
                if (i < Power - 1) ShiftLeft(fa_xi);
            }

            return new PowerFieldElement(this, vec);
        }

        private void ShiftLeft(FiniteFieldElement[] vec) {
            FiniteFieldElement head = vec[vec.Length - 1];
            for (int i = vec.Length - 1; i > 0; i--) vec[i] = vec[i - 1] + head * _Polynomial[i];
            vec[0] = head * _Polynomial[0];
        }

        public static FiniteField Create(FiniteField baseField, int power) {
            if (baseField == null) throw new ArgumentNullException("baseField");
            if (power < 2) throw new ArgumentOutOfRangeException("power");

            // Since the field is cyclic, there is only one finite field of a given prime power order (up to isomorphism).
            // For most practical purposes that means that we can pick any arbitrary monic irreducible polynomial.
            // We can abuse PowerField to do polynomial multiplication in the base field.
            var fakeField = new PowerField(baseField, Enumerable.Repeat(baseField.Zero, power).ToArray());
            var excluded = new HashSet<FiniteFieldElement>();
            for (int lpow = 1; lpow <= power / 2; lpow++) {
                int upow = power - lpow;
                // Consider all products of a monic polynomial of order lpow with a monic polynomial of order upow.
                int xl = (int)Math.Pow(baseField.Size, lpow);
                int xu = (int)Math.Pow(baseField.Size, upow);
                for (int i = xl; i < 2 * xl; i++) {
                    var pi = fakeField.Convert(i);
                    for (int j = xu; j < 2 * xu; j++) {
                        var pj = fakeField.Convert(j);
                        excluded.Add(-(pi * pj));
                    }
                }
            }

            for (int p = baseField.Size; true; p++) {
                var pp = fakeField.Convert(p) as PowerFieldElement;
                if (!excluded.Contains(pp)) return new PowerField(baseField, pp.Value.ToArray());
            }
        }

        public override string ToString() {
            var sb = new System.Text.StringBuilder();
            sb.AppendFormat("GF({0}) with primitive polynomial x^{1} ", Size, Power);
            for (int i = Power - 1; i >= 0; i--) sb.AppendFormat("+ {0}x^{1}", _Polynomial[i], i);
            sb.AppendFormat(" over base field ");
            sb.Append(_BaseField);
            return sb.ToString();
        }
    }

    internal class PowerFieldElement : FiniteFieldElement {
        private readonly PowerField _Field;
        private readonly FiniteFieldElement[] _Vector; // The version of Mono I have doesn't include IReadOnlyList<T>

        internal override FiniteField Field { get { return _Field; } }
        internal FiniteFieldElement[] Value { get { return _Vector; } }

        internal PowerFieldElement(PowerField field, params FiniteFieldElement[] vector) {
            if (field == null) throw new ArgumentNullException("field");
            if (vector == null) throw new ArgumentNullException("vector");
            if (vector.Length != field.Power) throw new ArgumentOutOfRangeException("vector");
            for (int i = 0; i < vector.Length; i++) if (vector[i].Field != field.BaseField) throw new ArgumentOutOfRangeException("vector[" + i + "]");

            _Field = field;
            _Vector = vector.ToArray();
        }
    }
}

Results
Unfortunately adding the rulers would take me about 15k characters past the post size limit, so they're on pastebin.

Would you be so kind to post your rulers for [50, 100] somewhere? I have a genetic algorithm I want to try out, feeding it some seed values.

@orlp, added a link.

As I suspected, the evolutionary algorithm can extract nothing of use out of these superior specimen. Even though initially it seemed that evolutionary algorithms might work (it pretty much instantly moves from invalid rulers to actual rulers), there is too much global structure required for the evolutionary algorithm to work.

Python 3, score 603001 / 734078 = 0.82144
Naive search combined with Erdős–Turan construction:
$$2pk + (k^2\bmod p),\, k \in [0, p-1]$$
For odd primes p this gives an asymptotically optimal golomb ruler.
def isprime(n):
    if n < 2: return False
    if n % 2 == 0: return n == 2
    k = 3
    while k*k <= n:
         if n % k == 0: return False
         k += 2
    return True

rulers = []
ruler = []
d = set()
n = 0
while len(ruler) <= 100:
    order = len(ruler) + 1
    if order > 2 and isprime(order):
        ruler = [2*order*k + k*k%order for k in range(order)]
        d = {a-b for a in ruler for b in ruler if a > b}
        n = max(ruler) + 1
        rulers.append(tuple(ruler))
        continue

    nd = set(n-e for e in ruler)
    if not d & nd:
        ruler.append(n)
        d |= nd
        rulers.append(tuple(ruler))
    n += 1

isuniq = lambda l: len(l) == len(set(l))
isruler = lambda l: isuniq([a-b for a in l for b in l if a > b])

assert all(isruler(r) for r in rulers)

rulers = list(sorted([r for r in rulers if 50 <= len(r) <= 100], key=len))
print(sum(max(r) for r in rulers))

I don't think this construction is asymptotically optimal: it yields a Golomb ruler of order p and length about 2p^2, whereas there exist Golomb rulers of order n and length about n^2 asymptotically.

@GregMartin Asymptotically there is no difference between 2p^2 and p^2.

Depends on your definition of "asymptotically", I guess, but to me, in this context they're very different.

Mathematica, score 276235/734078 < 0.376302
ruzsa[p_, i_] := Module[{g = PrimitiveRoot[p]},
  Table[ChineseRemainder[{t, i PowerMod[g, t, p]}, {p - 1, p}], {t, 1, p - 1}] ]

reducedResidues[m_] := Select[Range@m, CoprimeQ[m, #] &]

rotate[set_, m_] := Mod[set - #, m] & /@ set

scaledRuzsa[p_] := Union @@ Table[ Sort@Mod[a b, p (p - 1)],
  {a, reducedResidues[p (p - 1)]}, {b, rotate[ruzsa[p, 1], p (p - 1)]}]

manyRuzsaSets = Join @@ Table[scaledRuzsa[Prime[n]], {n, 32}];

tryGolomb[set_, k_] := If[Length[set] < k, Nothing, Take[set, k]]

Table[First@MinimalBy[tryGolomb[#, k] & /@ manyRuzsaSets, Max], {k, 50, 100}]

The function ruzsa implements a construction of a Golobm ruler (also called a Sidon set) found in Imre Z. Ruzsa. Solving a linear equation in a set of integers. I. Acta Arith., 65(3):259–282, 1993. Given any prime p, this construction yields a Golomb ruler with p-1 elements contained in the integers modulo p(p-1) (that's an even stronger condition than being a Golomb ruler in the integers themselves).
Another advantage of working in the integers modulo m is that any Golomb ruler can be rotated (the same constant added to all elements modulo m), and scaled (all elements multiplied by the same constant, as long as that constant is relatively prime to m), and the result is still a Golomb ruler; sometimes the largest integer is decreased significantly by doing so. So the function scaledRuzsa tries all of these scalings and records the results. manyRuzsaSets contains the results of doing this construction and scaling for all of the first 32 primes (chosen a bit arbitrarily, but the 32nd prime, 131, is well larger than 100); there are nearly 57,000 Golomb rulers in this set, which takes a good several minutes to compute.
Of course, the first k elements of a Golomb ruler themselves form a Golomb ruler. So the function tryGolomb looks at such a ruler made from any of the sets computed above. The last line Table... selects the best Golomb ruler it can, of every order from 50 to 100, from all the Golomb rulers found in this way.
The lengths found were:
{2241, 2325, 2399, 2578, 2640, 2762, 2833, 2961, 3071, 3151, 3194, 3480, 3533, 3612, 3775, 3917, 4038, 4150, 4237, 4368, 4481, 4563, 4729, 4974, 5111, 5155, 5297, 5504, 5583, 5707, 5839, 6077, 6229, 6480, 6611, 6672, 6913, 6946, 7025, 7694, 7757, 7812, 7969, 8139, 8346, 8407, 8678, 8693, 9028, 9215, 9336}

I was originally going to combine this with two other constructions, those of Singer and of Bose; but it seems that Peter Taylor's answer has already implemented this, so presumably I would simply recover those lengths.

I'm confused by your claim that working in the integers modulo m you can rotate / scale freely. Look at [0, 1, 4, 6] mod 7. If I add 1 we get [0, 1, 2, 5], which is not a Golomb ruler.

That's because you have to start with a mod-7 Golomb ruler for it to work. [0, 1, 4, 6] is not a mod-7 Golomb ruler because 1 – 0 equals 0 – 6 modulo 7, for instance.

While I was writing and debugging my finite field implementation in C# I wished I knew Mathematica better. Definitely one of the right languages for the job.

all_licenses: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/3.0/


text:
Knit me an ASCII-Sock
Introduction
Long story short: a few days ago I accidentally started knitting a sock and the pretty logical structure of the plain stitches lead me to the idea: Why don't we just ASCII-knit a sock?
Input
The input is an even integer N in the range [2,30].
Output
The output is a sock, obviously.
Structure
Because we will only use plain stitches you will only use v for a downward stitch and > for a stitch to the right. Let me explain the full structure with an example for N=8:
vvvvvvvv
vvvvvvvv
vvvvvvvv
vvvvvvvv
vvvvvvvv
vvvvvvvv
vvvvvvvv
vvvvvvvv
vvvvvvvv
vvvvvvvv
vvvvvvvv
vvvvvvvv
vvvvvvv>>>
vvvvvv>>>>>>
vvvvv>>>>>>>>>
vvvv>>>>>>>>>>>>
vvvv>>>>>>>>>>>>>>>>>>>>>>>>
vvvv>>>>>>>>>>>>>>>>>>>>>>>>>>
vvvv>>>>>>>>>>>>>>>>>>>>>>>>>>
vvvv>>>>>>>>>>>>>>>>>>>>>>>>

The first block is the upper part. It has the width N and the height N*1.5
Afterwards the heel begins. You subtract one v per row from the right side until you have N/2 as width. Then you add N/2 rows with the width N/2 beneath to finish the downward knits. Note that it gives a total of N/2 + 1 rows with the width N/2

first part:
vvvvvvvv    //begin with width N
vvvvvvvv
vvvvvvvv
vvvvvvvv
vvvvvvvv
vvvvvvvv
vvvvvvvv
vvvvvvvv
vvvvvvvv
vvvvvvvv
vvvvvvvv
vvvvvvvv    // N*1.5 times
vvvvvvv     // reduce
vvvvvv
vvvvv
vvvv        // till here
vvvv        // repeat N/2 times
vvvv
vvvv
vvvv

Now the side wards knits begin. You start with the first reduced row, fill > till you have the width N and add  two more. Continue that process including the first row with N/2 width
Fill the rows beneath to that level and add N*1.5 times > to that
Let the outer parts as they are and add on the row beneath the first and over the lase two more >'s until it is finished

second part:
vvvvvvvv
vvvvvvv>>>  //fill to N and add two >
vvvvvv>>>>>>     //repeat
vvvvv>>>>>>>>>
vvvv>>>>>>>>>>>>   //till here
vvvv>>>>>>>>>>>>>>>>>>>>>>>>       //fill to the row above and add N*1.5 times >
vvvv>>>>>>>>>>>>>>>>>>>>>>>>>>     //add two more '>' than the row above
vvvv>>>>>>>>>>>>>>>>>>>>>>>>>>     //add two more '>' than the row beneath
vvvv>>>>>>>>>>>>>>>>>>>>>>>>

Rules
Trailing new lines are ok.
You can choose between uppercase or lowercase V but it has to be consistent.
This is code-golf so the shortest code in bytes wins.
Test cases
N=2
vv
vv
vv
v>>>
v>>>>>>

N=4
vvvv
vvvv
vvvv
vvvv
vvvv
vvvv
vvv>>>
vv>>>>>>
vv>>>>>>>>>>>>
vv>>>>>>>>>>>>

N=6
vvvvvv
vvvvvv
vvvvvv
vvvvvv
vvvvvv
vvvvvv
vvvvvv
vvvvvv
vvvvvv
vvvvv>>>
vvvv>>>>>>
vvv>>>>>>>>>
vvv>>>>>>>>>>>>>>>>>>
vvv>>>>>>>>>>>>>>>>>>>>
vvv>>>>>>>>>>>>>>>>>>

N=12
vvvvvvvvvvvv
vvvvvvvvvvvv
vvvvvvvvvvvv
vvvvvvvvvvvv
vvvvvvvvvvvv
vvvvvvvvvvvv
vvvvvvvvvvvv
vvvvvvvvvvvv
vvvvvvvvvvvv
vvvvvvvvvvvv
vvvvvvvvvvvv
vvvvvvvvvvvv
vvvvvvvvvvvv
vvvvvvvvvvvv
vvvvvvvvvvvv
vvvvvvvvvvvv
vvvvvvvvvvvv
vvvvvvvvvvvv
vvvvvvvvvvv>>>
vvvvvvvvvv>>>>>>
vvvvvvvvv>>>>>>>>>
vvvvvvvv>>>>>>>>>>>>
vvvvvvv>>>>>>>>>>>>>>>
vvvvvv>>>>>>>>>>>>>>>>>>
vvvvvv>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
vvvvvv>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
vvvvvv>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
vvvvvv>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
vvvvvv>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
vvvvvv>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>

I hope a trailing newline is okay? (also, how does one accidentally start knitting a sock?)

@GregMartin Visited my grandmother and she asked me. It's hard to say no to her so here we go :) Finished the first one two days ago. I think it's pretty fine but I really underestimated the amount of work. So be happy if someone made socks for you.

Can I use v or V interchangeably? I'm going to assume for the moment that v is the only one allowed, but allowing V will effect the byte count

@GoldenRatio I changed the rules, now both is allowed.

@GurupadMamadapur As you can read in 'Input', only even numbers are allowed as input. Your program doesn't need to handle odd numbers

All the images seem Brainfuck code to me!

Pyth - 93 Bytes
J/Q2K"V"P*s*Q1.5+*QKbVJp*+Q_+1NK*+3*3N">";InQ2VJp*JK*+s*2aaNc-J 1 2c-J 1 2+3*Q3">";IqQ2pK*6">

Explanation:
J/Q2                             sets J to half of input
K"V"                             sets K to "V"
P                                prevents newline after main block of sock
*                                multiply the next 2 items
s*Q1.5                           rounds up input*1.5 (columns in main block)
+                                adds the next 2 item (this item gets multiplied by the previous item)
*QK                              multiplies "V" by input (resulting in each row)
b                                New line
VJ                               For all the numbers 0-((Q/2)-1) (as represented by N)
p                                print
*                                multiply the next 2 items
+Q_+1N                           input + -1*(1+N)
K                                "V"
*+3*3N">"                        3N+3 many ">"'s
;                                ends for statement
InQ2                             if input is not 2
VJ                               For all the numbers 0-((input/2)-1) (as represented by N)
p                                print
*JK                              input/2 many ">"'s
*+s*2aaNc-J 1 2c-J 1 2+3*Q3">"   int(abs(abs(N-((input/2)-1)/2)-((input/2)-1)/2)))+(input*3)+3 many ">"'s
;                                end of for statement
IqQ2                             if input is 2
pK*6">                           print "V" and 6 ">"'s

You can try it here!

Vyxal, 54 bytes
1.*(\v*,)½(n-\v*₴n:d+\>*,)½½ʁ(½\v*?T\>*Jnd\>*J:⅛,)¾(¼,

Boy oh boy, this is my new favorite answer now.
The satisfaction when I got it down from 79 bytes to 54 bytes was immense.
I bet you can't golf this anymore.
Try it Online!
Explanation
1.*(\v*,)½(n-\v*₴n:d+\>*,)½½ʁ(½\v*?T\>*Jnd\>*J:⅛,)¾(¼,

1.*(\v*,)     <Top part>
1.            Evaluates to 1.5
  *           Times input i
   (    )     Loop (with loop var n)
    \v*,      Repeat "v" n times and print it

½(n-\v*₴n:d+\>*,)     <Diagonal part>
½                     Halve input i
 (              )     Loop (with loop var n)
  n-\v*₴              Repeat "v" i minus n times and print it (without a newline)
        n:d           Duplicate n and double it
           +\>*,      Repeat ">" 2n plus n times

½½ʁ(½\v*?3*\>*Jnd\>*J:⅛,)     <Top half bumpy part>
½½ʁ                           Create range 0 to input i/4
   (                    )     Loop
    ½\v*                      Repeat "v" i/2 times
        ?T\>*J                Repeat ">" 3i times and join it
               nd\>*J:        Repeat ">" 2n times, merge and duplicate it
                      ⅛,      Push it to the global array and print it

¾(¼,                          <Bottom half bumpy part>
¾(                            Loop through the length of the global array
  ¼,                          Pop from the global array (to reverse it)

Mathematica, 104 bytes
""<>{t=Table;c=t[{"v"~t~#,">"~t~#2,"
"},#3]&;h=#/2;c[#,0,3h],c[#-i,3i,z={i,h}],c[h,6h+2Min[i-1,h-i],z]}&

Unnamed function taking a positive even integer as input and returning a string (with a trailing newline). Note that the newline between the two lines of code above is part of the code. The main work is done by the function c, whose ungolfed definition
c[#1, #2, #3] = Table[{Table["v", #1], Table[">", #2], "\n"}, #3]

creates a table of lists, each one consisting of several "v"s in a row followed by several ">"s followed by a newline. Fortunately, the range #3 for the outer Table can have a variable name, which can be referred to in the expressions #1 and #2; this allows the function c to be called with both constant and variable arguments. The rest is just computation, and the string-joining operator ""<> flattens the nested lists that arise for free.

Python, 3.5 183 177 bytes
n=int(input())
h=n//2
z=h%2
l=['v'*h+'>'*(2*n-h+3*h+i*2)for i in range(1*z+n//4)]
*map(print,['v'*n]*(3*h)+['v'*i+'>'*(3*(n-i))for i in range(h,n)][::-1]+[l,l[:-1]][z]+l[::-1]),

Explanation incoming.
Try it Online!

Just so that you know, TIO Nexus works better than TIO v2 at the moment as v2 is in "the early alpha stage"

To save one byte, you could assign 'v' top a variable.

@JackBates Duly noted.

@nedla2004 Do you mean this way : n,v=int(input()),'v' ? Then it saves nothing.

Groovy, 168 bytes
{n->h=n/2;z=(int)h%2;l=(0..<1*z+n/4).collect{'v'*h+'>'*(2*n-h+3*h+it*2)};(['v'*n]*(3*h)+(n-1..h).collect{'v'*it+'>'*(3*(n-it))}+(z==0?l:l[0..<-1])+l[-1..0]).join('\n')}

This is an unnamed closure. The approach is completely similar to my Python 3.5 answer.
Try it Online!

Batch, 279 bytes
@echo off
set i=for /l %%i in (2,2,%1)do
set v=
%i%call set v=v%%v%%v
%i%echo %v%&echo %v%&echo %v%
%i%call:a
%i%call set v=%%v%%]]]
%i:,2,=,4,%call:c&if %%i lss %1 call set v=%%v%%]]
%i:2,=4,%call:b
exit/b
:a
set v=%v:~1%]]]]]
:b
set v=%v:~0,-2%
:c
echo %v:]=^>%

Note: the second line ends in a space. Batch has problems with echoing >s in variables so I have to a) use a substitute b) call a subroutine to echo the corrected string. In fact I have three subroutines:

a is used for the heel. One v is removed each time and three >s are added. Two extra >s are added so that the subroutine can fall through into the b subroutine.
b is used for the bottom of the foot. Two >s are removed each time. The subroutine then falls through into the c subroutine.
c prints the current row of the sock, replacing the placeholder character so that the >s print properly.

The fourth line handles the width of the leg while the fifth line handles its length. The seventh line handles the length of the foot while the eighth line handles the top half of the foot height by adding two >> to each line after printing it. This works for sock sizes that are multiples of 4, and an adjustment is made for other sock sizes so that the longest line does not get repeated.

dc, 269 bytes
1k?d1.5*sadsdsb[[[v]nlb1-dsb0<f]dsfx[]pldsbla1-dsa0<c]dscxldd1-sisb2su[[[v]nlb1-dsbldli-<g]dsgx[[>]nlb1-dsblu+0<t]dstx[]plu2+suli1-dsilddsb2/!>m]dsmxld3*dsdsolidswst[[[v]nlw1-dsw0!>x]dsxx[2-]sa[2+]sy[[>]nlo1-dso0<g]dsgxldltli1+2/<yltli1+2/>adsdsolisw[]plt1-dst0!>f]dsfx

All right, well...this is massive. It wasn't easy, I'll tell you that much, what with all the register references, comparisons, macros, lack of string manipulation abilities, etc. Nonetheless, I had fun putting this together, and that's all that really matters, right? ;)
Try It Online!

C (gcc), 185 bytes
#define L for(i=0;i<h;i++,p(10,1))
p(c,n){n&&putchar(--n?p(c,n):c);}
f(n,i,k,h){h=n/2;for(k=3;k--;)L p(118,n);L{p(118,n-i-1);p(62,3*i+3);}L{p(118,h);p(62,3*n);p(62,2*(4*i/n?h-i-1:i));}}

Try it online!
Mostly boring, except I use a recursive function to print repeated characters.

161 bytes

PHP, 229 bytes
a quite lazy approach
for($m=$n=$p=$argv[1];$i++<$q=$n*1.5;)echo($d=str_pad)("
",$n+1,v);for(;$p-->$n/2;)echo$s="
".$d($d("",$p,v),$m+=2,">");$r=[$s.=$d("",$q,">")];for($k=$n>>2;--$k;)$r[]=$s.=">>";echo join($r),$n&2?"$s>>":"",join(array_reverse($r));

prints a leading newline. Run with -nr.
breakdown
// print first part + preparing some variables for later
for($m=$n=$p=$argv[1];$i++<$q=$n*1.5;)echo str_pad("\n",$n+1,v);

// print second part + prepare $s for third part
for(;$p-->$n/2;)echo$s="\n".str_pad(str_pad("",$p,v),$m+=2,">");

// third part: create an array of strings ...
$r=[$s.=str_pad("",$q,">")];
for($k=$n>>2;--$k;)$r[]=$s.=">>";
// ... then print it and its reverse
echo join($r),$n&2?"$s>>":"",join(array_reverse($r));

license: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/4.0/


text:
OOP: Overlapping Oriented Programming
One of the lesser known programming paradigms which seems rather fitting for code golfing is Overlapping Oriented Programming (OOP) *. When writing partly identical code, many bytes can be saved by simply overlapping the identical parts and remembering in some way where the two original code lines start. Your task is to write two overlapping programs or functions compress and decompress with the following specification:
*Do not use in production code, probably.
compress
compress takes two strings in any convenient format and overlaps them as much as possible. That is a string s with minimal length is returned such that both input strings are substrings of s. Additionally, some output which identifies the start and end indices of both string is returned.
Examples: (The exact IO-format is up to you)
compress("abcd", "deab") -> "deabcd" ((2,5),(0,3))
compress("abcd", "bc")   -> "abcd" ((0,3),(1,2))
compress("abc", "def")   -> "abcdef" ((0,2),(3,5)) or "defabc" ((3,5),(0,2))

decompress
decompress computes the inverse function of compress, that is given a string and two start and end indices (in the format in which they are returned by your compress), return the two original strings. You need only handle valid inputs.
The following equality should hold for all strings s1, s2:
(s1, s2) == decompress (compress (s1, s2))

Examples: (inverses of compress examples)
decompress "deabcd" ((2,5),(0,3)) -> "abcd" "deab"
decompress "abcd" ((0,3),(1,2))   -> "abcd" "bc"

decompress "abcdef" ((0,2),(3,5)) -> "abc" "def"
 or (whichever version your "compress" generates)
decompress "defabc" ((3,5),(0,2)) -> "abc" "def"

Scoring
Your score is the size of the string returned by calling compress("<code of compress>", "<code of decompress>"). As this is code-golf a lower score is better.
Example:
Assume the code for your function compress is c=abcd and the code for decompress is d=efghi. Then compress("c=abcd", "d=efghi") yields "c=abcd=efghi" (and the indices, but those don't affect scoring) so the score is length "c=abcd=efghi" = 12.
Additional Rules

In the spirit of this challenge, your compress and decompress must overlap in at least one character. You may achieve this trivially by adding a comment, but note that doing so will increase your score and there might be shorter solutions using inherently overlapping code.
compress and decompress should be able to handle strings containing any printable ASCII characters as well as all characters you used to define compress and decompress.
The indices can be zero- or one-indexed.
Your programs or functions don't have to actually be named compress and decompress.

Can you use different command line arguments to run you compress and decompress code?

Sure. You have to give two programs and the site policy allows command line arguments as long as they are counted, so you can give different command line arguments for each of your programs.

GNU Prolog, 105 points
s(U,L/M,C):-prefix(A,C),length(A,M),suffix(U,A),length(U,L).
o(A-B,C-X-Y):-length(C,_),s(A,X,C),s(B,Y,C).

(This requires GNU Prolog because prefix and suffix are not portable.)
Prolog has one interesting, major advantage for this challenge; you can write a function to handle multiple call patterns (i.e. not only can you give the function an input to get a corresponding output, you can give the function an output to get the corresponding input). As such, we can define a function that can handle both compression and decompression, leading to a 105-byte submission that defines a function o that works both ways round. (Incidentally, I mostly wrote it as a decompressor – as it's simpler – and got the compressor "for free". ) In general, we could expect a very short program in Prolog for this task, if not for the fact that it's so bad at string handling (both in terms of missing primitives, and in terms of the primitives in question having terribly long names).
The first argument to o is a tuple of strings, e.g. "abcd"-"deab". The second argument has a form like "deabcd"-4/6-4/4; this is a fairly standard nested tuple, and means that the string is "deabcd", the first string has length 4 and ends at the sixth character, the second string has length 4 and ends at the fourth character. (Note that a string in GNU Prolog is just a list of character codes, which makes debugging annoying because the implementation prefers the latter interpretation by default.) If you give o one argument, it'll output the other for you (thus working as a compressor if you give the first argument, and a decompressor if you give the second). If you give it both arguments, it'll verify that the compressed representation matches the strings given. If you give it zero arguments, it'll generate output like this:
| ?- o(X,Y).
X = []-[]
Y = []-0/0-0/0 ? ;

X = []-[]
Y = [_]-0/0-0/0 ? ;

X = []-[A]
Y = [A]-0/0-1/1 ? ;

…many lines later…

X = [A]-[B,A,C]
Y = [B,A,C]-1/2-3/3 ? ;

The above description of the I/O format is pretty much entirely just a description of the program, too; there's not very much in the program at all. The only subtlety is to do with evaluation order hints; we need to ensure that the program doesn't just use a search strategy that's guaranteed to terminate, but also produces the shortest possible output string.
When compressing, we start with length(C,_) ("C has a length"), which is a trick that I've used in many Prolog and Brachylog answers; if this is the first thing Prolog sees, it'll cause it to prioritise reducing the length of C over anything else. This ensures that we have a minimum-length C. The ordering of the constraints in s is carefully chosen so that the search will take a finite time for each possible candidate length of C; A is constrained by C (we don't know C, but we do know the target value we have for its length), M by A, U by A, and L by U, so none of the searches can take unlimited time.
When decompressing, we're given C directly by the user. This again ensures the program will run in finite time, due to the same sequence of constraints. (People who are aware of Prolog's evaluation order will note that the definition of s is very inefficient when decompressing; placing length(A,M) and length(U,L) first would be faster, but moving length(A,M) to the start could cause an infinite loop when compressing because neither A nor M is bound by anything at the time.)

Brachylog, 50 46 bytes
{Ċ∧Lċ₂l∧Lgj:?z{tT∧?h~cṪhlI∧ṪbhTl:I+-₁:I↔}ᵐ:L}

Try it online!
Decompress:
~{Ċ∧Lċ₂l∧Lgj:?z{tT∧?h~cṪhlI∧ṪbhTl:I+-₁:I↔}ᵐ:L}

Try it online!
Saved 5 bytes thanks to @ais523
Explanation
The good side of declarative languages is that we can reuse the same code for both compressing and decompressing. As such, the code for compress is exactly the same as for decompress, with an additional ~ at the beginning.
This ~ tells Brachylog to reverse the order of arguments (that is, use the Input as output and the Output as input). Since compress has no ~, it actually runs the predicate in the standard order. Since decompress has only one, it runs it with the Input as output and the Output as input.
That way, we can use the same code (modulo that extra ~) to both compress and decompress: compressing is providing the two strings as input and a variable as output, and decompressing is providing the indexes and compressed string as output and a variable as Input.
Obviously this also means that we have to be a bit explicit about our compressing code, so that the interpreter is capable of running it "backwards". This is why the compressor itself is a bit long.
Here is a breakdown of the code for Compress (and thus also of decompress):
{……………………………………………………………………}   Call that predicate the normal way (with swapped arguments
                                 for decompress)
   Ċ                           Input has two elements
   ∧Lċ₂l                       L is a string of any length (measuring its length forces it to
                                 take a specific length from 0 to +inf)
   ∧Lgj                        The list [L,L]
       :?z                     The list [[L, First elem of Input],[L,second elem of input]]
          {………………………………}ᵐ:L    Final output is the [M,L] where M is the result of mapping
                                 the predicate below on both elements of the zip
           tT                  The second element of the input is T
           ∧?h~cṪ              Anticoncatenate the first element of the input into [A,B,C]
           hlI                 I = length(A)
           ∧ṪbhTl:I+-₁         J = length(T) + I - 1
           :I↔                 Output = [I,J]

Jelly, 58 50 bytes
-1 byte thanks to ais523 (use ⁾ for a two-byte string)
This may well be quite golfable...
Compression takes two string arguments and returns a list:
[[[startA, lengthA], [startB, lengthB]], compressedString]
w³;w⁴$
0;J⁸ḣ;€
ç;ç@ÑẠ$ÐfLÐṂṪµ³,⁴L€ż@Ñ,

Decompression takes one argument (such a list) and returns two* strings:
,
⁾ṫḣżFv
Ḣç€Ṫ

Overlapped code:
w³;w⁴$
0;J⁸ḣ;€
ç;ç@ÑẠ$ÐfLÐṂṪµ³,⁴L€ż@Ñ,
⁾ṫḣżFv
Ḣç€Ṫ

One-indexed.
* this may not be obvious due to Jelly's implicit print formatting, so the code at TryItOnline linked to above has an extra byte (a Y at the end) to insert a line feed between the two in the printed output.
I started with a single program, which used the depth of the first (or only) argument to decide between compression and decompression, but having an unused link in the decompression code (the first line) and a single overlapping byte is shorter by seven bytes.
How?
ç;ç@ÑẠ$ÐfLÐṂṪµ³,⁴L€ż@Ñ, - Compression: stringA, stringB
ç                       - call the last link (2) as a dyad
  ç@                    - call the last link (2) as a dyad with reversed arguments
 ;                      - concatenate (gives all overlapping strings)
       Ðf               - filter keep:
      $                 -     last two links as a monad
    Ñ                   -         call the next link (1) as a monad
     Ạ                  -         All? (no zeros exist in that result)
          ÐṂ            - filter keep with minimal:
         L              -     length
            Ṫ           - tail (for when more than one exists)
             µ          - monadic chain separation (we now have the compressed string)
              ³,⁴       - [stringA, stringB]
                 L€     - length of €ach
                   ż@   - zip with reversed arguments with
                     Ñ  - next link (1) as a monad with the compressed string
                      , - paired with the compressed string

J0;⁸ḣ;€ - Link 2, possible overlaps: stringL, stringR
J       - range(length(stringL)) - [1,2,...,length(stringL)]
 0;     - zero concatenate       - [0,1,2,...,length(stringL)]
   ⁸    - stringL
    ḣ   - head (vectorises)      - [empty string, first char, first two, ..., stringL]
     ;€ - concatenate €ach with stringR

w³;w⁴$ - Link 1, substring indexes: stringX
w³     - first index of first program argument in stringX or 0 if not found
  ;    - concatenated with
     $ - last two links as a monad
   w⁴  -     first index of second program argument in stringX or 0 if not found

Ḣñ€Ṫ - Decompression: [[[startA, lengthA], [startB, lengthB]], compressedString], ?
Ḣ    - head - [[startA, lengthA], [startB, lengthB]]
   Ṫ - tail - compressedString
 ç€  - call the last link (2) as a dyad for €ach of the left list
     -- extra Y atom at TIO joins the resulting list of two strings with a line feed.

⁾ṫḣżFv - Link 2, extract a substring: [start, length], string
⁾ṫḣ    - string "ṫḣ"
   ż   - zip with [start, length] to yield [['ṫ', start],['ḣ', length]]
    F  - flatten, making a list of characters
     v - evaluate as Jelly code with the string as an argument
       - this evaluates as string.tail(start).head(length) yielding the substring

, - Link 1: only here to make an overlap with the compression program.

“ṫḣ” can be golfed by 1 byte by using Jelly's syntax for 2-character strings.

This is a question completely unrelated to the answer per se, but do you write the explanation of the code by hand or is there a tool generating it from code?

@tfrascaroli I write it by hand


text:
Hilbert's Grand Hotel
Introduction
Some of you may have heard of Hilbert's Grand Hotel. The manager there has lost his list of where the guests are staying but he still has the order in which they checked in. Each guest can not stay in a room with a room number less than their value and if a guest is added to a lower room, all guests in higher rooms with no empty space between them and the new guest are shifted up one room. Can you help him find where each of the guests are staying?
Requirements
Write a program which receives an ordered list of natural numbers as input and places them at their index. If there is already a value in that index, it is shifted up to the next entry in the list. This process repeats until the first empty (0 or undefined) space is found. Any undefined spaces between the current highest index and any new input will be filled by adding 0s. As this is Hilbert's Grand Hotel, rooms higher than the current highest occupied index do not exist.
Input and Output
Input will be an ordered list of natural numbers (allowed to be read through any accepted form of input)
Each number in the input is considered one guest arriving at the hotel and is in order of arrival
Output will be the final arrangement of guests (numbers)
Examples

Input: 1 3 1
Output: 1 1 3
Step by step:
  1
  Create room at index 1 and place 1 in it
  1 0 3
  Create rooms up to index 3 and place 3 in room 3
  1 1 3
  Shift the contents of room 1 up one room and place 1 in room 1
Input: 1 4 3 1 2 1
Output: 1 1 2 1 3 4
Step by step:
  1
  Create room at index 1 and place 1 in it
  1 0 0 4
  Create rooms up to index 4 and place 4 in room 4
  1 0 3 4
  Place 3 in room 3
  1 1 3 4
  Shift contents of room 1 up one room and place 1 in room 1
  1 2 1 3 4
  Shift contents of rooms 2 to 4 up one room and place 2 in room 2
  1 1 2 1 3 4
  Shift contents of rooms 1 to 5 up one room and place 1 in room 1
Input: 10
Output: 0 0 0 0 0 0 0 0 0 0 10
Step by step:
  0 0 0 0 0 0 0 0 0 10
  Create rooms up to room 10 and place 10 in room 10
Notes:
  Working with 0 indexed is fine and you may insert a 0 at the front of the output in that case

Standard loopholes are forbidden, shortest code in bytes wins

Jelly, 28 bytes
ṫœṡ0F
’⁹ḣ;⁸;ç@ð0ẋ;@
0;ç@/œr0

Try it online!

Explain please?

PHP 93 bytes
for(;$r=$g?$r+1:$g=$argv[++$i];[$g,$h[$r]]=[$h[$r],$g])$h=array_pad($h?:[],$g,0);print_r($h);

0 indexed. Uses a 2 in 1 loop that looks up the next guest after it gets a 0 (or a null form going beyond the current final room). Use like:
php -r "for(;$r=$g?$r+1:$g=$argv[++$i];[$g,$h[$r]]=[$h[$r],$g])$h=array_pad($h?:[],$g,0);print_r($h);" 1 4 3 1 2 1

Ungolfed:
for( $hotel = []; $room = $guest = $argv[ ++$i ]; )
{
    for( $hotel = array_pad( $hotel, $guest, 0 ); $guest; $room++ )
    {
        $last = $hotel[ $room ];
        $hotel[ $room ] = $guest;
        $guest = $last;
    }
}
print_r( $hotel );

Hahaha, that almost looks like it could be perl!

Haskell, 107 bytes
h=foldl(\h n->let{(b,a)=splitAt n h;(c,d)=span(>0)a;e=0:e;f=take n$b++e;g(0:m)=m;g m=m}in f++[n]++c++g d)[]

Try it online!

Jelly, 24 22 bytes
aŻƤSƇḢßoɓ@ooƑ?
Ṭ×)ç@ƒ0

Try it online!
I still have the strangest feeling that something closer to Dennis's answer should beat this, but I haven't had any luck with that so far. I have twisted the first line into something that sounds like some kind of exploding can (from my original recursive solution; I thought it saved a byte but I didn't look at the output too carefully so it actually ties but it does look funnier).
aŻƤSƇḢßoɓ@ooƑ?    Helper link: add a pre-positioned guest to the hotel

        ɓ         With reversed arguments (hotel left, guest right):
           oƑ?    Unless both arrays have a nonzero at the same location,
          o       return with the guest placed in its empty room;
        ɓ         if they do:
         @        (guest left, hotel right)
a                 Overlay the guest's zeroes on the hotel.
 ŻƤ               Prepend another zero to each prefix,
   SƇḢ            and take the first with a nonzero sum.
      ß           Call this link again with that as the new guest
       o          and with the hotel having the previous guest replace the new one.

Ṭ×)ç@ƒ0    Main link

  )        For each guest value,
Ṭ          obtain a list of all zeroes then a one at that index,
 ×         and multiply that by the same guest value.
     ƒ0    Starting with 0, reduce that by
   ç@      the helper with reversed arguments.

Essentially just boots occupants out of their rooms in sequence until there's no collision, considering a new guest as equivalent to someone evicted from the last room they can't start in.
Dennis's answer with some simple golfs for fair comparison, and an explanation for those still wondering:
Jelly, 26 bytes
’⁹ḣ;⁸;ṫ@œṡ0Fɗð0ẋ;@
ç@ƒ0œr0

Try it online!
’⁹ḣ;⁸;ṫ@œṡ0Fɗð0ẋ;@    Helper link: add a guest on the left to the hotel on the right

             ð        Let the right argument to the rest of the link be
              0ẋ      zero repeated by the value of the guest
                ;@    concatenated to the end of the hotel.
 ⁹ḣ                   Take elements from that right argument numbering
’                     one less than the guest,
   ;⁸                 concatenate the guest,
     ;      ɗ         and concatenate:
      ṫ@              remove elements from that right argument numbering the guest value,
        œṡ0           split those remaining around the first zero,
           F          and flatten them back together.

ç@ƒ0œr0    Main link

ç@ƒ0       With the hotel starting as 0, add each guest to the hotel in sequence.
    œr0    Trim trailing zeroes.

You outgolfed Dennis? You absolute madlad. Have an upvote.

JavaScript (ES6), 144 120 bytes
Saved 20B thanks to Arnauld and 11B thanks to Neil
a=>a.map((d,c)=>b[d]?(b.splice(d,0,d),b.splice([...b,0].findIndex((e,g)=>g&&!e),1)):b[d]=d,b=[])&&[...b].map(c=>c|0)

Usage
You can assign the function to the variable f and the list should be given as an array. Example:
f=a=>a.map((d,c)=>b[d]?(b.splice(d,0,d),b.splice([...b,0].findIndex((e,g)=>g&&!e),1)):b[d]=d,b=[])&&[...b].map(c=>c|0)
f([1,4,3,1,2,1])

Output
The output is in an array as well. Since Javascript works zero-indexed, there is an extra 0 up front.
[0, 1, 1, 2, 1, 3, 4]

I didn't really test that, but could (c+'').split\,`.map(Number)` do the job?

Clever trick, and yes, it does work. Thanks.

Oops, I forgot to test testcase 2, and it turns out my function doesn't support adding things to the front of the array...

I believe you could do c.map(n=>n|0) rather than (c+'').split`,`.map(Number).

@ETHproductions The problem is that map() doesn't iterate at all on undefined values in the array. (That said, I'm pretty sure there's a shorter way than the one I suggested.)

Took some effort, but the bug is now fixed. Please let me know if you have any improvements...

Do you really need a for? a=>a.map((d,c)=>b[d]?(b.splice(d,0,d),b.splice([...b,0].findIndex((e,g)=>g&&!e),1)):b[d]=d,b=[])&&(b+'').split\,`.map(Number)` should work.

You're right, added. Thanks!

[...c].map(n=>n|0) should work.

You're right! Learned another trick today...

JavaScript (ES6), 86 bytes
f=([n,...a],r=[],p=n,m=r[p],_=r[p]=n)=>n?m?f([m,...a],r,p+1):f(a,r):[...r].map(n=>n|0)

Leading zero in the result because JavaScript is 0-indexed.

Mathematica, 98 bytes
Fold[If[Length@#<#2,PadRight@##~Append~#2,Join[Take@##,{#2},Drop@##/.{a___,0,b__}->{a,b}]]&,{0},#]&

Unnamed function taking a list of positive integers and returning a 0-indexed list of integers. The whole If function takes a partially completed list and the next integer to insert as arguments. If the next integer exceeds the length of the partial list, PadRight@##~Append~#2 increases the partial list accordingly; otherwise, Join[Take@##,{#2},Drop@##/.{a___,0,b__}->{a,b}]] inserts the next integer into its position, then throws away the first 0 found after it. Fold[...,{0},#] applies this function repeatedly to the original list, starting with the empty hotel {0}, and outputs the final hotel list.

JavaScript (ES6), 81
Using 0 indexing
l=>l.map(v=>(s=(p,v,w=r[p])=>(w&&s(p+1,w),r[p]=v))(v,v),r=[])&&[...r].map(x=>~~x)

Less golfed
l => {
  s = ( // recursively insert a value, shifting present values up until it finds an empty spot
       p, // insert position
       v, // value to insert
       w = r[p] // value a current position
      ) => (
         w && s(p+1,w), // if there is a value, put it in the next position
         r[p] = v // now the current place is empty
      );
  r = []; // initialize the result array
  l.map( // execute the following for each value in input
     v => s(v,v) // insert value
  );
  return [...r].map(x=>~~x) // fill the missing places with 0
}

Test

F=
l=>l.map(v=>(s=(p,v,w=r[p])=>(w&&s(p+1,w),r[p]=v))(v,v),r=[])&&[...r].map(x=>~~x)

out=x=>O.textContent+=x+'\n'

out([1,3,1]+' -> '+F([1,3,1]))
out([10]+' -> '+F([10]))

function go() {
   var v = I.value.match(/\d+/g).map(x=>+x)
   out(v +' -> ' + F(v))
}
go()
<input id=I value='1 4 3 1 2 1'><button onclick='go()'>go</button>

R, 133 bytes
a=scan()
z=rep(0,max(a))
for(i in a){b=match(0,z[-1:-i])+i
if(z[i])z[i:b]=c(i,z[i:(b-1)])else(z[i]=i)
z=c(z,0)}
z[1:max(which(z!=0))]

To avoid problems with bad indexing, I pad with some zeros, and then strip them at the end. This is maybe not the best solution, but it works.

Python, 134 125 116 bytes
def a(b):
 r=[]
 i=0
 for n in b:
    d=n
    while n:
     if i<d:r.extend([0]*(d-i));i=d
     n,r[d-1]=r[d-1],n;d+=1
 return r

Valid for both Python 2.7.13 and 3.6.0. This code functions by means of a swapping the held value with the value contained at each index until the held value is 0. If it reaches an index not yet in the array, it adds zeros to the end of the array until the array contains that index. Thanks to Wheat Wizard and xnor for golfing off 9 bytes each

Instead of using spaces for all your indents you could use a space for level one indents, a tab for level two and a tab followed by a space for level three.

I was working with a stupid editor which converted tab to 4 spaces XP

The conditions for while and if don't need parens. You can put multiple statements on one line separated by ; like if(i<d):r.extend([0]*(d-i));i=d unless there's control flow in the later statements.