Project Euler Bonus Problem

Among all non-square integers $n$ with absolute value not exceeding $10^3$, find the value of $n$ such that $\cos(\pi \sqrt n)$ is closest to an integer.

Searching for the solution

This seems easy enough. We just compute $\cos(\pi \sqrt n)$ for all non-square integers $|n| \le 1000$ and keep track of which ones are closest to an integer.

First thing to notice is that when $n > 0$ then $\cos(\pi \sqrt n)$ just oscillates between $-1$ and $1$ but maybe it'll get close to an integer. When $n < 0$ we actually have

which grows exponentially.

Seeing that $\cosh \left( \pi \sqrt{1000} \right) \approx 7 \times 10^{42}$ we will need at least 43 digits of precision and maybe a nice 32 digit buffer for the fractional part. 75 digits worth of precision requires $75 \log_2(10) \approx 249$ bits of precision which we'll round up to $250$. So Float64 will not be enough and neither will 128-bit floats (which are provided by Quadmath.jl in Julia). We'll use BigFloat and calculate the precision required based on the limit we want to search up to.

function required_precision_bits(limit; fractional_digits=32)
    val = cosh(big(π) * sqrt(big(limit)))
    integer_digits = ceil(Int, log10(val))
    total_digits = integer_digits + fractional_digits
    return ceil(Int, total_digits * log2(10))
end

We can then code up the search. Each value of $n$ is independent so this is embarrassingly parallel: the range from $-\text{limit}$ to $\text{limit}$ is divided into one chunk per thread, each thread computes $\cos(\pi \sqrt n)$ and its distance to the nearest integer for its own chunk, and the results are concatenated and sorted by distance (with a plain loop as the fallback when Julia is started with a single thread):

function distance_to_nearest_integer(x)
    return abs(x - round(x))
end

function process_chunk(chunk_start, chunk_end)
    local_results = Vector{Tuple{Int,BigFloat,BigFloat}}()
    for n in chunk_start:chunk_end
        n == 0 && continue
        isqrt(abs(n))^2 == abs(n) && continue

        if n > 0
            val = cos(big(π) * sqrt(big(n)))
        else
            val = cosh(big(π) * sqrt(big(-n)))
        end

        dist = distance_to_nearest_integer(val)
        push!(local_results, (n, val, dist))
    end
    return local_results
end

function find_closest_cos_to_integer(limit)
    precision_bits = required_precision_bits(limit)
    @info "Using $precision_bits bits of precision for limit=$limit"

    # On Julia 1.10, precision is global: set it once and restore it only after
    # every worker finishes, including when a worker throws an exception.
    results = setprecision(BigFloat, precision_bits) do
        num_chunks = Threads.nthreads()

        if num_chunks == 1
            return process_chunk(-limit, limit)
        end

        # Multi-threaded: chunk the range from -limit to limit
        total_range = 2 * limit + 1
        chunk_size = cld(total_range, num_chunks)

        tasks = @sync map(1:num_chunks) do i
            chunk_start = -limit + (i - 1) * chunk_size
            chunk_end = min(chunk_start + chunk_size - 1, limit)
            Threads.@spawn process_chunk(chunk_start, chunk_end)
        end

        return reduce(vcat, fetch.(tasks))
    end

    sort!(results, by=x -> x[3])
    log_top_results(results)

    return results[1][1]
end

The top 10 are logged so we can see how close the runners-up get:

using Printf

function log_top_results(results)
    @info "Top 10 values of n where cos(π√n) is closest to an integer:"
    for i in 1:min(10, length(results))
        n, val, dist = results[i]
        @info @sprintf("%d: n = %d, value ≈ %.10e, distance ≈ %.10e", i, n, val, dist)
    end
end

Going up to $n = 10^3$ requires 250 bits of precision and produces the correct answer in 801.681 μs using all the threads on the machine.

Going up to $n = 10^6$ requires 4641 bits of precision and actually does not find a closer value after searching for 7.648 s!

Just knowing the solution

There are nine Heegner numbers. I don't fully understand them and arriving at them seems quite nontrivial so I won't try to describe them, but they are related to almost integers.

The largest Heegner number, 163, leads to the almost integer

Apparently this number is sometimes called the Ramanujan constant even though it was probably not found by Ramanujan himself, although he did find others like $e^{\pi \sqrt{58}}$ which is also an almost integer.

The existence of this almost integer explains the connection between Heegner numbers and this bonus problem, and may explain why searching up to $n = 10^6$ did not produce any closer answers.