# 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 Jan 25 '17 at 2:37
• @LegionMammal978 With the best runtime performance. It's a fastest-code challenge. – Noah L Jan 25 '17 at 2:40
• oeis.org/A037274 – HyperNeutrino Jan 25 '17 at 2:45
• 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. – JungHwan Min Jan 25 '17 at 2:59
• Is using a probabilistic primality test such as Miller-Rabin allowed? – miles Jan 25 '17 at 18:04

# 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... – JungHwan Min Jan 25 '17 at 3:19
• 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. – Noah L Jan 25 '17 at 3:57
• Since we're not golfing: there's CompositeQ for !PrimeQ (which also ensures that your answer doesn't loop forever on input 1). – Martin Ender Jan 25 '17 at 8:58
• 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. – Mario Jan 26 '17 at 14:34
• @NoahL Mathematica can be tested online. meta.codegolf.stackexchange.com/a/1445/34718 – mbomb007 Jan 26 '17 at 15:05

# PyPy 5.4.1 64bit (linux), HP(80) ~ 1.54s

The 32bit version will time slightly slower.

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

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

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

home-prime.py

import math
import my_math
import mpqs

max_trial = 1e10
max_pollard = 1e22

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

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

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

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

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

def ec_multiply(k, p, C, n):
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  # 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. – primo Jan 26 '17 at 17:50 • 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 – mbomb007 Jan 26 '17 at 21:23 • 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. – primo Jan 26 '17 at 22:17 • You should help him out. :) – mbomb007 Jan 26 '17 at 22:30 • 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 Jan 26 '17 at 22:35 # 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 Jan 25 '17 at 18:09 • @P.Ktinos I think so too... anyway it would be too big to include. – Mario Jan 25 '17 at 18:17 # 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. – primo Jan 28 '17 at 20:52
• Also, use feature "say";. – primo Jan 28 '17 at 21:08