Calculate an integer sequence derived from prime factors
J, 47 45
#@((~.@,[:(+/@{:*+/@:*/)2 p:{:)^:_)`2:@.(=&1)
It's possible this would be much shorter without using ^:_
, but my brain is sufficiently fried already.
Edit: (47->45) Double coupon day.
Usage:
higley =: #@((~.@,(+/@{:*+/@:*/)@(2&p:)@{:)^:_)`2:@.(=&1)
higley 1
2
higley"0 (1 + i. 10)
2 1 1 10 1 11 1 9 5 10
Golfscript, 68 67 62 61 chars
[.]({[.2@{1$1$%{)}{\1$/1$}if}*;;].,*0+{+}*.2$?@@.@+\@)!}do;,(
This is an expression: it takes n
on the stack and leaves the result on the stack. To turn it into a program which takes n
from stdin and prints the result to stdout, replace the leading [
with ~
The heart of it is [.2@{1$1$%{)}{\1$/1$}if}*;;]
(28 chars) which takes the top number on the stack and (by an incredibly inefficient algorithm) generates a list of its prime factors. C-style pseudocode equivalent:
ps = [], p = 2;
for (int i = 0; i < n; i++) {
if (n % p == 0) {
ps += p;
n /= p;
}
else p++;
}
The 0+
just before {+}*
is to handle the special case n==1
, because Golfscript doesn't like folding a binary operation over the empty list.
One of the non-prime fixpoints is 27; I found this without using the program by considering the mapping (pa ->
a2p), which is a fixpoint if a == p(a-1)/2, and trying small a
. (a==1
gives the fixpointedness of primes).
Searching with the program turns up a second fixpoint: 30 = (2+3+5)*3
Appendix: proof that there are only two non-prime fixpoints
Notation: sopfr(x)
is the sum of prime factors of x
, with repetition (A001414). Omega(x)
is the number of prime factors of x
(A001222). So the Higley successor function is h(x) = sopfr(x) Omega(x)
Suppose we have a fixpoint N = h(N)
which is a product of n=Omega(N)
primes.
N = p_0 ... p_{n-1} = h(N) = n (p_0 + ... + p_{n-1})
Basic number theory: n
divides into p_0 ... p_{n-1}
, so w=Omega(n)
of those primes are the prime factors of n
. Wlog we'll take them to be the last w
. So we can divide both sides by n
and get
p_0 ... p_{n-w-1} = p_0 + ... + p_{n-1}
or
p_0 ... p_{n-w-1} = p_0 + ... + p_{n-w-1} + sopfr(n)
Given that all of the primes p_0
to p_{n-w-1}
are greater than 1, increasing any of them increases the LHS more than the RHS. So for a given n
, we can enumerate all candidate solutions.
In particular, there can be no solutions if the LHS is greater than the RHS setting all the "free" primes to 2. I.e. there are no solutions if
2^{n-w} > 2 (n-w) + sopfr(n)
Since sopfr(n) <= n
(with equality only for n=4 or n prime), we can make the weaker statement that there are no fixpoints if
2^{n-w} > 3 n - 2 w
Holding w
fixed we can select different values of n
satisfying w=Omega(n)
. The smallest such n
is 2^w
. Note that if 2^{n-w}
is at least 3 (i.e. if n-w>1
, which is true if n>2
) then increasing n
while holding w
constant will increase the LHS more than the RHS. Note also that for w>2
and taking the smallest possible n
the inequality is satisfied, and there are no fixpoints.
That leaves us with three cases: w = 0
and n = 1
; w = 1
and n
is prime; or w = 2
and n
is semi-prime.
Case w = 0
. n = 1
, so N
is any prime.
Case w = 1
. If n = 2
then N = 2p
and we require p = p + 2
, which has no solutions. If n = 3
then we have pq = p + q + 3
and two solutions, (p=2, q=5)
and (p=3, q=3)
. If n = 5
then 2^4 > 3 * 5 - 2 * 1
, so there are no further solutions with w = 1
.
Case w = 2
. If n = 4
then N = 4pq
and we require pq = p + q + 4
. This has integer solution p=2, q=6
, but no prime solutions. If n = 6
then 2^4 > 3 * 6 - 2 * 2
, so there are no further solutions with w = 2
.
All cases are exhausted, so the only non-prime fixpoints are 27 and 30.
Ruby, 92 characters
f=->i{r=[i];(x=s=0;(2..i).map{|j|(s+=j;x+=1;i/=j)while i%j<1};r<<i=s*x)until r.uniq!;r.size}
This solution assumes higley(1) is actually 2, not 1 (see balpha's comment above):
(1..10).map &f
=> [2, 1, 1, 10, 1, 11, 1, 9, 5, 10]