József Kürschák proved in 1908 that the function

is never an integer for 0 < m < n. In particular, the harmonic numbers

are never integers for n > 1.

The function f(m, n) can get arbitrarily close to any integer value by taking m and n large enough, but it can never exactly equal an integer.

For this post, I’d like to look at how close  f(m, n) comes to an integer value when 0 < m < n ≤ N for some large N, say N = 100,000.

Computation strategy

The most naive way to approach this would be to compute f(m, n) for all m and n and keep track of which values were closest to integers. This would be wasteful since it would recompute the same terms over and over. Instead, we could take advantage of the fact that

Instead of working with f(m, n), it will be convenient to work with just its fractional part

because it won’t hurt to throw away the integer parts as we go. The values of m and n minimizing g(m, n) will be the values for which f(m, n) comes closest to an integer from above. The values of m and n maximizing g(m, n) will be the values for which f(m, n) comes closest to an integer from below.

We could calculate a matrix with all values of g(m, n), but this would take O(N²) memory. Instead, for each n we will calculate g(m, n), save the maximum and minimum values, then overwrite that memory with g(m, n + 1). This approach will only take O(N) memory.

Floating point error

When we compute f(m, n) for large values of n, can we rely on floating point arithmetic?

If N = 100,000, f(m, n) < 16 = 2⁴. A floating point fraction has 53 bits, so we’d expect each addition to be correct to within an error of 2⁻⁴⁹ and so we’d expect our total error to be less than 2⁻⁴⁹ N.

Python code

The following code computes the values of g(m, n) closest to 0 and 1.

import numpy as np

N = 100_000
f_m = np.zeros(N+1) # working memory

# best values of m for each n
min_fm = np.zeros(N+1)
max_fm = np.zeros(N+1)

n = 2
f_m[1] = 1.5

for n in range(3, N+1):
    f_m[n-1] = 1/(n-1)
    f_m[1:n] += 1/n
    f_m[1:n] -= np.floor(f_m[1:n])
    min_fm[n] = np.min(f_m[1:n])
    max_fm[n] = np.max(f_m[1:n])    

print(min(min_fm[3:]))
print(max(max_fm))

This reports a minimum value of 5.2841953035454026e-11 and a maximum value of 0.9999999996613634. The minimum value is closer to 0 than our (pessimistic) error estimate, though the maximum value is further from 1 than our error estimate.

Modifying the code a bit shows that the minimum occurs at (27134, 73756), and this the input to g that is within our error estimate. So we can be confident that it is the minimum, though we can’t be confident of its value. So next we turn to Mathematica to find the exact value of g(27133, 73756) as a rational number, a fraction with 32024 digits in the numerator and denominator, and convert it to a floating point number. The result agrees with our estimate in magnitude and to four significant figures.

So in summary

with an error on the order of 10⁻¹¹, and this is the closest value of f(m, n) to an integer, for 0 < m < n ≤ 100,000.

Related posts

• Accurately computing harmonic numbers
• Extending harmonic numbers
• Numerators of harmonic numbers

The “Freshman’s dream” is the statement

(x + y)^p = x^p + y^p

It’s not true in general, but it is true mod p if p is a prime. It’s a cute result, but it’s also useful in applications, such as finite field computations in cryptography.

Here’s a demonstration of the Freshman’s dream in Python.

>>> p = 5
>>> [((x + y)**p - x**p - y**p) % p for x in range(p) for y in range(p)]
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]

Here’s an example using a prime too large for verify the results by looking at the output.

>>> import numpy as np
>>> p = 103
>>> v = [((x + y)**p - x**p - y**p) % p for x in range(p) for y in range(p)]
>>> np.all( np.array(v) == 0 )
True

You can use the same code to show that the Freshman’s dream is not true in general if p is not a prime, and it’s not true in general if p is a prime but the exponent is less than p.

I recently saw someone post [1] that 987654321/123456789 is very nearly 8, specifically 8.0000000729.

I wondered whether there’s anything distinct about base 10 in this. For example, would the ratio of 54321_six and 12345_six be close to an integer? The ratio is 4.00268, which is pretty close to 4.

What about a larger base? Let’s try base 16. The expression

0xFEDCBA987654321 / 0x123456789ABCDEF

in Python returns 14. The exact ratio is not 14, but it’s as close to 14 as a standard floating point number can be.

For a base b, let denom(b) to be the number formed by concatenating all the digits in ascending order and let num(b) be the number formed by concatenating all the digits in descending order.

Then for b > 2 we have

The following Python code demonstrates [2] that this is true for b up to 1000.

num = lambda b: sum([k*b**(k-1) for k in range(1, b)])
denom = lambda b: sum([(b-k)*b**(k-1) for k in range(1, b)])

for b in range(3, 1001):
    n, d = num(b), denom(b)
    assert(n // d == b-2)
    assert(n % d == b-1)

So for any base the ratio is nearly an integer, namely b − 2, and the fractional part is roughly 1/b^b−2.

When b = 16, as in the example above, the result is approximately

14 + 16⁻¹⁴ = 8 + 4 + 2 + 2⁻⁵⁶

which would take 60 bits to represent exactly, but a floating point fraction only has 53 bits. That’s why our calculation returned exactly 14 with no fractional part.

[1] I saw @ColinTheMathmo post it on Mastodon. He said he saw it on Fermat’s Library somewhere. I assume it’s a very old observation and that the analysis I did above has been done many times before.

[2] Why include a script rather than a proof? One reason is that the proof is straight-forward but tedious and the script is compact.

A more general reason that I give computational demonstrations of theorems is that programs are complementary to proofs. Programs and proofs are both subject to bugs, but they’re not likely to have the same bugs. And because programs made details explicit by necessity, a program might fill in gaps that aren’t sufficiently spelled out in a proof.

In 2013, John Conway and Alex Ryba published a brief note [1] on how to convert identities involving sine and cosine into identities involving Fibonacci and Lucas numbers.

Fibonacci and Lucas

The Fibonacci numbers F_n are defined by F₀ = 0, F₁ = 1, and

F_n+2 = F_n + F_n+1

for n > 1. Similarly, Lucas numbers L_n are defined by the same recurrence relation

L_n+2 = L_n + L_n+1

but starting with L₀ = 2 and L₁ = 1.

The recipe

The recipe for turning trigonometric identities into Fibonacci and Lucas number identities is as follws.

1. Arguments become subscripts.
2. Replace sines by ½ iⁿ F_n.
3. Replace cosines by ½ iⁿ L_n.
4. Insert a factor of (−5)^k in front of every term that contains 2k or 2k + 1 sine terms.

The first step is admittedly unclear. It’s unclear in [1] as well. If sine or cosine takes an argument of

pα + qβ

then this becomes the subscript

pa + qb

where the Latin letters are integers and the Greek letters are angles.

There’s no proof in [1] why the recipe should work, but an earlier paper [2] gives more details.

Examples

Double angle

For a simple example, let’s start with the double angle identity

sin 2θ = 2 sin θ cos θ.

This becomes

½ i²ⁿ F_2n = 2 (½ iⁿ F_n)(½ iⁿ L_n)

which simplifies to

F_2n = F_n L_n.

Sum angle

For a little more involved example, let’s look at the sum identity

cos(θ + φ) = cos θ cos φ − sin θ sin φ.

This becomes

½ i^{a + b} L_{a + b} = ½ i^a L_a ½ i^b L_b −(−5) ½ i^a F_a ½ i^b F_b

which simplifies to

2L_{a + b} = L_a L_b  + 5 F_a F_b.

Triple angle

In the previous examples the powers of i in the recipe above cancelled out. This example will show that they are necessary.

We start with the triple angle identity

sin 3θ = 3 sin θ − 4 sin³ θ.

This becomes

½ i³ⁿ F_3n = 3 (½) iⁿ F_n − 4(−5)(½ iⁿ F_n)³

which simplifies to

F_3n = 3 (−1)ⁿ F_n + 5 (F_n)³.

Demonstration

The following Python code demonstrates that the Fibonacci / Lucas identities above are correct.

def F(n):
    if n == 0: return 0
    if n == 1: return 1
    return F(n-1) + F(n-2)

def L(n):
    if n == 0: return 2
    if n == 1: return 1
    return L(n-1) + L(n-2)

for n in range(10):
    assert(F(2*n) == F(n) * L(n))

for n in range(10):
    for m in range(10):
        assert(2*L(m + n) == L(m)*L(n) + 5*F(m)*F(n))

for n in range(10):
    assert(F(3*n) == 3*(-1)**n * F(n) + 5*F(n)**3)

Related posts

• Convert trig identities to hyperbolic identities
• Inverse Fibonacci numbers
• Lucas numbers and Lucas pseudoprimes
• Trig mnemonic diagrams

[1] John Conway and Alex Ryba. Fibonometry. The Mathematical Gazette, November 2013, pp. 494–495

[2] Barry Lewis, Trigonometry and Fibonacci Numbers, The Mathematical Gazette, July 2007, pp. 216–226

This morning I took an old blog post and turned it into an X thread. I think the thread is easier to read. More expository and less rigorous.

The post and thread look at generalizations of the fact that every integer and its fifth power end in the same digit. The conclusion is that n and n^k end in the same digit base b if b is square-free and k = φ(b) + 1. So, for example, 10 is square-free (i.e. not divisible by a non-trivial square) and φ(10) + 1 = 5.

Benjamin Clark replied suggesting looking at λ(b) + 1 in addition to φ(n) + 1.

Euler and Carmichael

To back up a bit, φ(n) is Euler’s totient function, the number of positive integers less than and relatively prime to n. Robert Carmichael’s totient function λ(n) is a closely related function, one that replaced Euler’s function in implementations of RSA encryption—more specifically in RSA decryption—something I wrote about earlier this week.

Euler’s generalization of Fermat’s little theorem says if a is relatively prime to m, then

a^φ(n) = 1 (mod n).

Carmichael’s totient λ(n) is the smallest number such that

a^λ(n) = 1 (mod n)

when a is relatively prime to n.

Sometimes φ(n) = λ(n), namely for n = 1, 2, 4, and every odd prime power. Otherwise λ(n) is a proper divisor of φ(n).

Carmichael’s conjecture

Carmichael conjectured that for every n there is at least one other integer m ≠ n such that φ(m) = φ(n). He proved that this is true at least for n less than 10³⁷. The conjecture is now known to be true for n less than 10¹⁰¹⁰.

Fermat numbers

The nth Fermat number is defined by

Pierre Fermat conjectured that the F(n) were prime for all n, and they are for n = 0, 1, 2, 3, and 4, but Leonard Euler factored F(5), showing that it is not prime.

Tangent numbers

The nth tangent number is defined by the Taylor series for tangent:

Another way to put it is that the exponential generating function for T(n) is tan(z).

Fermat primes and tangent numbers

Here’s a remarkable connection between Fermat numbers and tangent numbers, discovered by Richard McIntosh as an undergraduate [1]:

F(n) is prime if and only if F(n) does not divide T(F(n) − 2).

That is, the nth Fermat number is prime if and only if it does not divide the (F(n) − 2)th tangent number.

We could duplicate Euler’s assessment that F(5) is not prime by showing that 4294967297 does not divide the 4294967295th tangent number. That doesn’t sound very practical, but it’s interesting.

Update: To see just how impractical the result in this post would be for testing whether a Fermat number is prime, I found an asymptotic estimate of tangent numbers on OEIS,  and estimated that the 4294967295th tangent number has about 80 billion digits.

[1] Richard McIntosh. A Necessary and Sufficient Condition for the Primality of Fermat Numbers. The American Mathematical Monthly, Vol. 90, No. 2 (Feb., 1983), pp. 98–99

A rational triangle is a triangle whose sides have rational length and whose area is rational. Can any two rational numbers be sizes of a rational triangle? Surprisingly no. You can always find a third side of rational length, but it might not be possible to do so while keeping the area rational.

The following theorem comes from [1].

Let q be a prime number not congruent to 1 or −1 mod 8. And let ω be an odd integer such that no prime factor of q^2ω − 1 is congruent to 1 mod 4. Then a = q^ω + 1 and b = q^ω − 1 are not the sides of any rational triangle.

To make the theorem explicitly clear, here’s a Python implementation.

from sympy import isprime, primefactors

def test(q, omega):
    if not isprime(q):
        return False
    if q % 8 == 1 or q % 8 == 7:
        return False
    if omega % 2 == 0:
        return False
    factors = primefactors(q**(2*omega) - 1)
    for f in factors:
        if f % 4 == 1:
            return False
    return True

We can see that q = 2 and ω = 1 passes the test, so there is no rational triangle with sides a = 2¹ + 1 = 3 and b = 2¹ − 1 = 1.

You could use the code above to search for (a, b) pairs systematically.

from sympy import primerange

for q in primerange(20):
    for omega in range(1, 20, 2):
        if test(q, omega):
            a, b = q**omega + 1, q**omega - 1
            print(a, b)

This outputs the following:

3 1
9 7
33 31
129 127
8193 8191
32769 32767
131073 131071
524289 524287
4 2
6 4
126 124
14 12

Note that one of the outputs is (4, 2). Since multiplying the sides of a rational triangle by a rational number gives another rational triangle, there cannot be a rational triangle in which one side is twice as long as another.

As the author in [1] points out, “There are undoubtedly many pairs not covered by [the theorem cited above]. It would be of interest to characterize all pairs.” The paper was published in 1976. Maybe there has been additional work since then.

Related posts

• Rational right triangles
• The congruent number problem
• Do perimeter and area determine a triangle?
• Approximation with Pythagorean triangles

[1] N. J. Fine. On Rational Triangles. The American Mathematical Monthly. Vol. 83. No. 7, pp. 517–521.