Triangular Fibonacci numbers

Saturday 17 June 2017This is more than seven years old. Be careful.

Yesterday in my post about 55, I repeated Wikipedia’s claim that 55 is the largest number that is both triangular and in the Fibonacci sequence. Chris Emerson commented to ask for a proof. After a moment’s thought, I realized I had no idea how to prove it.

The proof is in On Triangular Fibonacci Numbers, a dense 10-page excursion into number theory I don’t understand.

While I couldn’t follow the proof, I can partially test the claim empirically, which leads to fun with Python and itertools, something which is much more in my wheelhouse.

I started by defining generators for triangular numbers and Fibonacci numbers:

def tri():
    """Generate an infinite sequence of triangular numbers."""
    n = 0
    for i in itertools.count(start=1):
        n += i
        yield n
        
print(list(itertools.islice(tri(), 50)))
[1, 3, 6, 10, 15, 21, 28, 36, 45, 55, 66, 78, 91, 105, 120, 136, 153, 171,
 190, 210, 231, 253, 276, 300, 325, 351, 378, 406, 435, 465, 496, 528, 561,
 595, 630, 666, 703, 741, 780, 820, 861, 903, 946, 990, 1035, 1081, 1128,
 1176, 1225, 1275]
def fib():
    """Generate an infinite sequence of Fibonacci numbers."""
    a, b = 1, 1
    while True:
        yield a
        b, a = a, a+b
        
print(list(itertools.islice(fib(), 50)))
[1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377, 610, 987, 1597, 2584,
 4181, 6765, 10946, 17711, 28657, 46368, 75025, 121393, 196418, 317811,
 514229, 832040, 1346269, 2178309, 3524578, 5702887, 9227465, 14930352,
 24157817, 39088169, 63245986, 102334155, 165580141, 267914296, 433494437,
 701408733, 1134903170, 1836311903, 2971215073, 4807526976, 7778742049,
 12586269025, 20365011074]

The Fibonacci sequence grows much faster!

My first impulse was to make two sets of the numbers in the sequences, and intersect them, but building a very large set took too long. So instead I wrote a function that took advantage of the ever-increasing nature of the sequences to look for equal elements in two monotonic sequences:

def find_same(s1, s2):
    """Find equal elements in two monotonic sequences."""
    try:
        i1, i2 = iter(s1), iter(s2)
        n1, n2 = next(i1), next(i2)
        while True:
            while n1 < n2:
                n1 = next(i1)
            if n1 == n2:
                yield n1
                n1 = next(i1)
            while n2 < n1:
                n2 = next(i2)
            if n1 == n2:
                yield n1
                n1 = next(i1)
    except StopIteration:
        return

And a function to cut off an infinite sequence once a certain ceiling had been reached:

def upto(s, n):
    """Stop a monotonic sequence once it gets to n."""
    for i in s:
        if i > n:
            return
        yield i

Now I could ask, what values less than quadrillion are in both the triangular numbers and the Fibonacci sequence?:

>>> list(find_same(upto(fib(), 1e15), upto(tri(), 1e15)))
[1, 3, 21, 55]

This doesn’t prove the claim for all numbers, but it shows that 55 is the largest number under a quadrillion that is in both sequences.

Another way to do this is to take advantage of the asymmetry in growth rate. The Fibonacci sequence up a quadrillion is only 72 elements. Make that a set, then examine the triangular numbers up to quadrillion and keep the ones that are in the Fibonacci set. And I’m certain there are other techniques too.

I can’t explain why, but composable generators please me. This was fun. :)

Comments

[gravatar]
Terry Jan Reedy 8:29 PM on 17 Jun 2017
For anyone interested, the following may help get at least an idea of the linked proof.

* Easy algebra shows 'n = m*(m+1)//2 (triangular)' if and only if '(8n+1) = (2m+1)**2 (odd square)', as in 55 = 10*11//2 and 441 = 21*21. Squareness is easier to work with. The paper actually shows that fib(10) = 55 is the largest fib such that 'times 8 plus 1' is square.

* The Fib series is produced by the simplest non-trivial 2-term linear recurrence. The Lucas series produced by the next simplest that is not a multiple or translation of the Fib series. Both have many properties, and they also have some joint properties, as in (2) to (4).

* For fixed count n, mod(m, n) = m % n maps int m to one of n equivalence classes, represented by 0, 1, ... (n-1). Either the definition of some f(m) or a proof of some property thereof may depend on m % n. Case analysis is *extensively* used in section 3.

* Mathematicians are conservative about dropping old names and notation developed before the modern understanding of functions. A 'Jacobi Symbol' is not a symbol in the modern sense but a function mapping an int and a positive odd int to (-1, 0, 1). Various web pages discuss the properties, calculation, and usefulness of the function.

* The Jacobi Criterion theorem combines Fibonacci numbers, Lucas numbers, the mod function, and the Jacobi function. I believe "(a, vn) = 1" means that a and v sub n are relatively prime, gcd(a, vn) = 1. The Wikipedia page did not mention 'proper' versus 'improper' Jacobi symbols, so I don't know what 'proper' means here.
[gravatar]
Thanks for that Ned and Terry!
[gravatar]
I also like composable generators, especially when they're in the standard library.
import heapq
import itertools

def tri():
    return itertools.accumulate(itertools.count(1))

def fib():
    a, b = 1, 1
    while True:
        yield a
        b, a = a, a+b

def find_same(xs, ys):
    return (v for v, gp in itertools.groupby(heapq.merge(xs, ys))
            if len(list(gp)) > 1)

def upto(xs, n):
    return itertools.takewhile(lambda x: x
        
[gravatar]
@Thomas, thanks for other implementations. I'm not sure I would use the find_same, but the upto is nice... :)
[gravatar]
Not sure what happened to upto(), but I hope you get the idea ;-)

I would use my version of find_same() because (to me) it's obviously correct: merge the steams then look for repeated elements in the result. What I find surprising is that merge lives in heapq rather than itertools.
[gravatar]
Jurjen N. E. Bos 1:10 AM on 9 Jan 2023

Just a small optimization: You can define tri() simply as accumulate(count(1)).

Add a comment:

Ignore this:
Leave this empty:
Name is required. Either email or web are required. Email won't be displayed and I won't spam you. Your web site won't be indexed by search engines.
Don't put anything here:
Leave this empty:
Comment text is Markdown.