Calculate an integer sequence derived from prime factors

Create a function, expression, or program which does the following:

1. Take the prime factors of any number and sum them. For example, the prime factors of 28 are 2 2 7, summed to 11.
2. Multiply the result by the number of prime factors for the given number. E.g, 28 has 3 prime factors which sum to 11. 11 * 3 is 33.
3. Repeat the process recursively, storing the resulting list (which starts with the original number), until you reach a number that is already included in the list. Halt without adding that final number, so that the list contains no duplicates. The progression for 28 is 28 33, because 33 results in 28 again.
4. Count the elements in the resulting list. In the case of 28, the answer is 2.

Here are the results for 0<n<=10, so you can check your algorithm.

2 1 1 10 1 11 1 9 5 10

(As balpha pointed out, the answer for higley(1) is 2, from the list 1 0. I originally had 1, due to a bug in my original algorithm written in J.)

Since I'm a conceited SOB and haven't found this on OEIS, let's call this the "Higley Sequence," at least for the duration of this round of code golf. As an added bonus, find the first two n having the lowest higley(n) where n is not prime and n>1. (I think there are only two, but I can't prove it.)

This is standard code golf, so as usual the fewest keystrokes win, but I ask that you please upvote clever answers in other languages, even if they are verbose.

• Why is highley(1) == 1? One has no prime factors, so the resulting list in 4) is [1, 0], so highley(1) == 2 as I see it.
– balpha
Apr 26, 2011 at 10:49
• Can we assume that the input number and intermediate values will be no greater than 2^31-1 (i.e. fits in a signed 32-bit integer)? Apr 26, 2011 at 11:43
• @Peter Taylor Sure. Apr 26, 2011 at 15:03
• In case anyone finds it helpful, the OEIS sequences which are vaguely related and may provide some inspiration are A001414, A001222, and A002217. Apr 26, 2011 at 16:12
• since you haven't commented I assume you haven't noticed: I proved that there are only the two non-prime fixpoints and added it as an appendix to my post. Apr 30, 2011 at 7:35

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

• Wow! A J solution that's shorter than a GolfScript solution. First one I've seen. (I'm a big fan of J.) Apr 26, 2011 at 23:09
• You can shorten this considerably by using a slightly different algorithm: #@((~.@,((+/*#)@:q:)@{:)^:_)2:@.(=&1), which is 38 characters. Apr 26, 2011 at 23:22
• Wow, I tried to figure out how to do it with q: but was trying to manhandle it into my 2 p: solution so I didn't get it. Obvious in retrospect. Apr 26, 2011 at 23:28
• The fact that you can look at that explosion of characters and say it "obvious in retrospect" simply blows my mind. One of these days I should check out Golfscript or J. Apr 27, 2011 at 1:10
• @Casey I felt the same way at one time, but the more J you learn and use, the more it sort of "leaps out at you," though I still see things I have to puzzle out. One helpful thing to know about J is that if you add a . or : after a symbol, it creates a new symbol, e.g, {, {., and {: all mean different things, but {- (for example) is definitely a sequence of two things, { and -. Apr 27, 2011 at 5:12

[.]({[.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. • Found those two same fixpoints using pencil-and-paper: 27 and 30. I would agree with OP it seems those are the only two. Apr 26, 2011 at 21:23 • Next interesting question might be. Are there infinitely many higley(x) = 2? How about is there a way to generate arbitrary higley(x), such as higley(x) = 100? Apr 26, 2011 at 21:29 • Very nice! I'm a J guy but just might have to learn GolfScript. Apr 26, 2011 at 21:53 • @mellamokb I think there are a number of interesting questions with this sequence. For instance, if we consider the sequence of numbers generated for each n before it is counted, are there any non-prime n after 49 for which said sequence does not end in 28? Apr 26, 2011 at 22:00 • Another interesting question to ask is whether there is a simple function of n which bounds higley(n) above. (That would allow simplifying the loop considerably - just iterate f(n) times and then discard duplicates). Apr 26, 2011 at 23:02 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]  Octave - 109 chars l=[input('')];while size_equal(unique(l),l);n=factor(l(1));l=[sum(n)*length(n) l];endwhile;disp(length(l)-1);  Husk, 9 bytes LU¡o§*LΣp  Explanation LU¡o§*LΣp input → n ¡o create an infinite list: p get the prime factors L length §* times Σ sum U longest unique prefix of the infinite list L length of the prefix  Jelly, 9 bytes ÆfS×LƊ$ƬL


Try it online!

Explanation

ÆfS×LƊ$ƬL Main monadic link Ƭ Repeat and collect results until a duplicate occurs$     (
Æf            Prime factorization
Ɗ        (
S             Sum
×            Multiplied by
L           Length
Ɗ        )
\$     )
L   Length


MATL, 19 bytes

ivt0)Yftnws*yy-]xn


Try it online!

Pyth, 10 bytes

l.u*sJPNlJ


Try it online!

Explanation:

 .u       Q | N = input, then N =
s       |  sum of
PN    |   list of prime factors of N
J      |   (call this list J)
*        |  times
lJ  |  length of J
.u         | until a previous result is repeated. return the list of results
l           | print the length of the list