## Happy bases

A place to discuss the implementation and style of computer programs.

Moderators: phlip, Moderators General, Prelates

LucasBrown
Posts: 299
Joined: Thu Apr 15, 2010 2:57 am UTC
Location: Poway, CA

### Happy bases

A happy number is defined as follows:
• Given a number n, repeatedly replace n with the sum of the squares of its digits.
• This will eventually reach a cycle of numbers of at most 3 digits: the cycle will either be 1 --> 1 --> 1 --> etc, or consist of number(s) greater than 1.
• If the cycle is 1, the number is happy.
This is clearly a base-dependent property. For example, all numbers in bases 2 and 4 are happy, while in base 10, 2 is unhappy. A base in which all numbers are happy numbers is called a happy base.

2 and 4 are the only known happy bases. Wikipedia claims that searches for happy bases have gone as high as 500,000,000; the citation for this claim leads to OEIS A161872, where the claim is made with no further citation.

To determine whether a base b is a happy base, we first note that applying the happy transformation to any number of 4 or more digits yields a lesser number. It therefore suffices to check all numbers up to 3 digits. Furthermore, the maximum result of the happy transformation of a number of 3 or fewer digits is 3 * (b-1)2. It therefore suffices to check all numbers less than or equal to 3 * (b-1)2.

We can therefore write the following Python code to find the smallest unhappy number in base b, or return 0 if b is a happy base:

Code: Select all

`def happybase(b):    if b < 5: return (0, 2, 0, 2, 0)[b]    for n in xrange(2, 3 * (b-1)**2 + 1):        f, g = n, n        while g != 1:            f, x = divmod(f, b)            f, y = divmod(f, b)            f = f*f + y*y + x*x            g, x = divmod(g, b)            g, y = divmod(g, b)            g = g*g + y*y + x*x            g, x = divmod(g, b)            g, y = divmod(g, b)            g = g*g + y*y + x*x            if f == g != 1: return n    return 0`

We use the cycle-finding trick employed in Pollard's rho algorithm to determine when a happy sequence has entered a cycle. We can also do this in C for a significant speedup:

Code: Select all

`long unsigned happybase(long unsigned b) {     long unsigned n, f, g, x, y, z;    if (b == 1) return 2;    z = 3 * (b-1) * (b-1);    for (n = 2; n <= z; n++) {        f = g = n;        while (g != 1) {            x = f % b; f /= b;            y = f % b; f /= b;            f = f*f + y*y + x*x;            x = g % b; g /= b;            y = g % b; g /= b;            g = g*g + y*y + x*x;            x = g % b; g /= b;            y = g % b; g /= b;            g = g*g + y*y + x*x;            if (f == g && g != 1) return n;        }    }    return 0;}`

The C code has the obvious limitation that the variables will overflow when b is sufficiently large, and the 500,000,000 in that OEIS comment is getting close to that point. We also find that as b gets large, evaluating this function gets increasingly slower (in an average sense). So I'd like to do three things:

• Find a more efficient algorithm to determine whether a number is happy in a given base.
• Rewrite the C code to use arbitrary-precision integers. I have no idea how to do this.
• Find some theorems that can tell us whether a number n is happy in base b without having to evaluate n's happy sequence in base b. For example, if log2 ( logn b ) is an integer, then n is happy in base b.
Any help would be greatly appreciated. https://pypi.python.org/pypi/labmath
3239×21412529+1 is prime!

Tub
Posts: 475
Joined: Wed Jul 27, 2011 3:13 pm UTC

### Re: Happy bases

LucasBrown wrote:Find a more efficient algorithm to determine whether a number is happy in a given base.

One thing I've noticed is that you immediately forget all previous calculations. If you've verified all numbers from 1 to n-1 to be happy, you only need to calculate the sequence for n until reaching any number smaller than n. You don't need to reach 1.

Also, when evaluating the sequence of n, you've likely reached a few numbers greater than n. If n ends up being happy, you could mark all the larger numbers as happy - at least as long as there's enough memory.

Implementing both shortcuts means that, for a happy base, every step of the happy algorithm will mark exactly one number as happy. You could thus evaluate the base b calculating no more than 3*(b-1)^2 steps of the sequence. O(n^2) is not perfect, but being able to check each number in (amortized) constant time is probably as good as it'll get.

LucasBrown wrote:Rewrite the C code to use arbitrary-precision integers. I have no idea how to do this.

For arbitrary precision, I've used https://gmplib.org/ before.

You're probably OK with 128bit integers for quite a while, and those can be implemented more efficiently than arbitrary precision. Some compilers support them directly, there's support in boost, or one can trivially find a homebrew implementation, e.g. https://github.com/calccrypto/uint128_t (the latter even has a combined divmod operation you will find useful!)

Qaanol
The Cheshirest Catamount
Posts: 3069
Joined: Sat May 09, 2009 11:55 pm UTC

### Re: Happy bases

LucasBrown wrote:Furthermore, the maximum result of the happy transformation of a number of 3 or fewer digits is 3 * (b-1)2. It therefore suffices to check all numbers less than or equal to 3 * (b-1)2.

Moreover, when n ≥ b2 then happy(n) < n, so you only have to check numbers below b2.
wee free kings

Xanthir
My HERO!!!
Posts: 5423
Joined: Tue Feb 20, 2007 12:49 am UTC
Location: The Googleplex
Contact:

### Re: Happy bases

Proof:

The largest happy(n) value for an n less than b2 is happy(b2-1) = 2*(b-1)2. This is larger than b2, but less than 2*b2.

The next smallest n value that produces a larger happy value is happy(2*b2-1) = 2(b-1)2+1. The n value has more than doubled, but the happy result has only increased by 1, so happy(n) is now less than n. The next local maxima is 3*b2-1, whose happy result is only +4 over happy(b2-1), while n has increased by 50%. This continues until you reach b3-1, whose happy result is only 3*(b-1)2, or 50% larger than happy(b2-1), while n has increased by a factor of b. So all happy values between b2 and b3-1 (all 3-digit values) are less than n.

By observation, all 4-digit values give a happy result with less digits.

Thus, only 1 and 2 digit values can have happy(n) > n, and so you only need to test 1- and 2-digit values.
(defun fibs (n &optional (a 1) (b 1)) (take n (unfold '+ a b)))

jaap
Posts: 2094
Joined: Fri Jul 06, 2007 7:06 am UTC
Contact:

### Re: Happy bases

Tub wrote:
LucasBrown wrote:Find a more efficient algorithm to determine whether a number is happy in a given base.

One thing I've noticed is that you immediately forget all previous calculations. If you've verified all numbers from 1 to n-1 to be happy, you only need to calculate the sequence for n until reaching any number smaller than n. You don't need to reach 1.

Also, when evaluating the sequence of n, you've likely reached a few numbers greater than n. If n ends up being happy, you could mark all the larger numbers as happy - at least as long as there's enough memory.

If you are searching for happy bases (i.e. moving on the the next base as soon as you have found an unhappy number in that base) then I don't think these optimisations are worth it. As seen in OEIS A161872, the number 2 is unhappy in nearly every base. These optimisations only help in cases where 2 is happy, which are so few that it does not seem worth the effort.

notzeb
Without Warning
Posts: 629
Joined: Thu Mar 08, 2007 5:44 am UTC
Location: a series of tubes

### Re: Happy bases

This might help: there is a number n bigger than 1 with happy(n) = n in base b if and only if b^2+1 is composite.

Proof: write n = bx+y with 0 <= x,y < b, then we are trying to solve the equation bx+y = x^2 + y^2. Multiplying by 4 and rearranging, this is equivalent to b^2+1 = (b-2x)^2 + (2y-1)^2, and b^2+1 can be written as the sum of a square and an odd square in a nontrivial way if and only if it is composite.
Zµ«V­jÕ«ZµjÖ­Zµ«VµjÕ­ZµkV­ZÕ«VµjÖ­Zµ«V­jÕ«ZµjÖ­ZÕ«VµjÕ­ZµkV­ZÕ«VµjÖ­Zµ«V­jÕ«ZµjÖ­ZÕ«VµjÕ­ZµkV­ZÕ«ZµjÖ­Zµ«V­jÕ«ZµjÖ­ZÕ«VµjÕ­Z

Qaanol
The Cheshirest Catamount
Posts: 3069
Joined: Sat May 09, 2009 11:55 pm UTC

### Re: Happy bases

Xanthir wrote:Proof:

That is not a complete proof.
wee free kings

Xanthir
My HERO!!!
Posts: 5423
Joined: Tue Feb 20, 2007 12:49 am UTC
Location: The Googleplex
Contact:

### Re: Happy bases

Qaanol wrote:
Xanthir wrote:Proof:

That is not a complete proof.

I left the remainder as an exercise for the reader. Feel free to provide it; the proof is trivial.
(defun fibs (n &optional (a 1) (b 1)) (take n (unfold '+ a b)))

Qaanol
The Cheshirest Catamount
Posts: 3069
Joined: Sat May 09, 2009 11:55 pm UTC

### Re: Happy bases

You are talking to the person who left the entire proof as an exercise.

Regardless, you did not address the b2 ≤ n ≤ 2(b−1)2 case sufficiently.
wee free kings

Xanthir
My HERO!!!
Posts: 5423
Joined: Tue Feb 20, 2007 12:49 am UTC
Location: The Googleplex
Contact:

### Re: Happy bases

Ah, you're right, what I wrote wasn't quite enough for that range. Easy fix, tho.

Going from XX to 1XX increases the number by b², but the happy value by 1. The diff between number and happy value is high at b-1, decreases steadily, then, if the base is big enough, recovers and gets highest again at b²-1. At both b-1 and b²-1, the diff is less than b², so moving from XX to 1XX guarantees that the number will not be smaller than the happy value.

The rest of the proof then applies to 2XX, etc.
(defun fibs (n &optional (a 1) (b 1)) (take n (unfold '+ a b)))

LucasBrown
Posts: 299
Joined: Thu Apr 15, 2010 2:57 am UTC
Location: Poway, CA

### Re: Happy bases

notzeb wrote:This might help: ... b^2+1 can be written as the sum of a square and an odd square in a nontrivial way if and only if it is composite.

Could you explain that a bit more? https://pypi.python.org/pypi/labmath
3239×21412529+1 is prime!

notzeb
Without Warning
Posts: 629
Joined: Thu Mar 08, 2007 5:44 am UTC
Location: a series of tubes

### Re: Happy bases

LucasBrown wrote:
notzeb wrote:This might help: ... b^2+1 can be written as the sum of a square and an odd square in a nontrivial way if and only if it is composite.

Could you explain that a bit more?
I just came down with a major fever and can't really think straight, but here goes. The key fact is that the Gaussian integers - aka Z[i] - are a unique factorization domain. In Z[i], b^2+1 factors as (b+i)(b-i), and every way of factoring b^2+1 into a pair of conjugate factors in Z[i] corresponds to a way of writing b^2+1 as a sum of two squares. Since b and 1 are relatively prime, no factor of b^2+1 can be congruent to 3 (mod 4), so each prime factor is either 2 or a prime congruent to 1 (mod 4). Since 2 ramifies in Z[i] you have to be a little careful if b^2+1 is even, but in this case it's easy to see that b^2+1 can be written as the sum of a square and an odd square in at least two ways directly. If b^2+1 is odd, then every prime factor is 1 (mod 4), so now we use a fact (due to Fermat, I think?) that says that every prime congruent to 1 (mod 4) can be written in the form (x+yi)(x-yi) with x+yi, x-yi primes in Z[i], in a unique way up to multiplication by units, with x+yi and x-yi relatively prime in Z[i] - finally we get rid of the unit ambiguity to assume that x,y are positive and that y is odd. So if, for instance, (x+yi) divides b+i, say with b+i = (x+yi)(m+ni), then we get a new factorization of b^2+1 as ((x-yi)(m+ni))*((x+yi)(m-ni)), and this gives us a new way of writing b^2+1 as a sum of two squares.
Zµ«V­jÕ«ZµjÖ­Zµ«VµjÕ­ZµkV­ZÕ«VµjÖ­Zµ«V­jÕ«ZµjÖ­ZÕ«VµjÕ­ZµkV­ZÕ«VµjÖ­Zµ«V­jÕ«ZµjÖ­ZÕ«VµjÕ­ZµkV­ZÕ«ZµjÖ­Zµ«V­jÕ«ZµjÖ­ZÕ«VµjÕ­Z

LucasBrown
Posts: 299
Joined: Thu Apr 15, 2010 2:57 am UTC
Location: Poway, CA

### Re: Happy bases

That appears to be correct, but it's been several years since since I (briefly) studied quadratic fields, so I wouldn't go so far as to give it my seal of approval. https://pypi.python.org/pypi/labmath
3239×21412529+1 is prime!

Return to “Coding”

### Who is online

Users browsing this forum: No registered users and 11 guests