Happy bases

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

Moderators: phlip, Moderators General, Prelates

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

Happy bases

Postby LucasBrown » Tue Jun 16, 2015 10:37 pm UTC

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.

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

Re: Happy bases

Postby Tub » Wed Jun 17, 2015 3:36 pm UTC

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!)

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

Re: Happy bases

Postby Qaanol » Thu Jun 18, 2015 3:21 am UTC

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

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

Re: Happy bases

Postby Xanthir » Thu Jun 18, 2015 5:47 am UTC

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)))

User avatar
jaap
Posts: 2073
Joined: Fri Jul 06, 2007 7:06 am UTC
Contact:

Re: Happy bases

Postby jaap » Thu Jun 18, 2015 5:53 am UTC

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.

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

Re: Happy bases

Postby notzeb » Thu Jun 18, 2015 8:00 am UTC

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

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

Re: Happy bases

Postby Qaanol » Fri Jun 19, 2015 1:11 am UTC

Xanthir wrote:Proof:

That is not a complete proof.
wee free kings

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

Re: Happy bases

Postby Xanthir » Sat Jun 20, 2015 3:31 am UTC

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)))

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

Re: Happy bases

Postby Qaanol » Sat Jun 20, 2015 4:15 am UTC

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

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

Re: Happy bases

Postby Xanthir » Sat Jun 20, 2015 5:39 pm UTC

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)))

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

Re: Happy bases

Postby LucasBrown » Sun Jun 21, 2015 1:33 am UTC

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?

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

Re: Happy bases

Postby notzeb » Sun Jun 21, 2015 10:53 pm UTC

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

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

Re: Happy bases

Postby LucasBrown » Mon Jun 22, 2015 2:41 am UTC

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.


Return to “Coding”

Who is online

Users browsing this forum: No registered users and 8 guests