Detailed mathematical approach

Problem Summary

For \(0<q<1\), the quantity of interest is the Lambert-type series

$\mathcal{S}_k(q)=\sum_{d\ge 1}\frac{d^k q^d}{1-q^d}.$

When \(q\) is extremely close to \(1\), the denominator \(1-q^d\) stays tiny for many values of \(d\), so direct summation converges very slowly. The implementation therefore uses a direct sum only for a small checkpoint and evaluates the large target case through an Eisenstein-series transformation.

Mathematical Approach

The key observation is that odd-power Lambert series are Fourier expansions of even-weight Eisenstein series, so a slowly converging sum near \(q=1\) can be converted into a rapidly evaluable modular expression.

Step 1: Rewrite the Lambert Series with Divisor Sums

Expand the geometric denominator:

$\frac{q^d}{1-q^d}=\sum_{r\ge 1} q^{dr}.$

Substituting this into the original series gives

$\mathcal{S}_k(q)=\sum_{d\ge 1}\sum_{r\ge 1} d^k q^{dr}=\sum_{n\ge 1}\sigma_k(n)q^n,$

where

$\sigma_k(n)=\sum_{d\mid n} d^k$

is the divisor-power sum. So the problem is really about evaluating a \(q\)-series generated by \(\sigma_k(n)\). The direct checkpoint in the implementations is the case \(k=1\) with \(q=1-2^{-4}=15/16\).

Step 2: Match Odd Powers with an Eisenstein Series

For the modular branch, the implementations use odd \(k\) of the form

$k=2m-1,$

with the supported weights \(2m\in\{4,8,16\}\). The normalized Eisenstein series satisfies

$E_{2m}(\tau)=1-\frac{4m}{B_{2m}}\sum_{n\ge 1}\sigma_{2m-1}(n)e^{2\pi i n\tau},$

where \(B_{2m}\) is the Bernoulli number of weight \(2m\). If we write

$q=e^{2\pi i\tau},$

then the Lambert series becomes

$\mathcal{S}_{2m-1}(q)=\frac{B_{2m}}{4m}\bigl(1-E_{2m}(\tau)\bigr).$

This is the exact identity behind the shortcut. The implemented Bernoulli constants are

$B_4=B_8=-\frac{1}{30},\qquad B_{16}=-\frac{3617}{510}.$

Step 3: Convert \(q=1-2^{-p}\) into a Small Positive Height

The target inputs use

$q=1-2^{-p}.$

Set

$t=-\log q,\qquad y=\frac{t}{2\pi}.$

Then

$q=e^{-t}=e^{-2\pi y}=e^{2\pi i(iy)},$

so the Eisenstein series is evaluated at the purely imaginary point

$\tau=iy.$

For large \(p\), the quantity \(2^{-p}\) is tiny, hence \(t\) and \(y\) are tiny as well. That is exactly the regime where the original Lambert series is hardest to sum term by term.

Step 4: Use the Modular Transformation

Eisenstein series obey the modular relation

$E_{2m}\!\left(-\frac{1}{\tau}\right)=\tau^{2m}E_{2m}(\tau).$

Substituting \(\tau=iy\) gives

$E_{2m}(iy)=(iy)^{-2m}E_{2m}(i/y).$

In the implemented cases, \(2m\in\{4,8,16\}\), so \(m\) is even and therefore

$ (iy)^{-2m}=y^{-2m}. $

Now \(y\) is very small, so \(i/y\) has very large imaginary part. That forces the Fourier expansion of \(E_{2m}(i/y)\) to be exponentially close to \(1\):

$E_{2m}(i/y)=1+O\!\left(e^{-2\pi/y}\right).$

Hence the leading term is

$E_{2m}(iy)\approx y^{-2m},$

and the Lambert series is approximated by

$\mathcal{S}_{2m-1}(q)\approx \frac{B_{2m}}{4m}\left(1-y^{-2m}\right).$

This is exactly the expression evaluated in the modular branch of the code.

Step 5: Worked Example \((k,p)=(3,8)\)

Here

$m=\frac{k+1}{2}=2,\qquad q=1-2^{-8}=\frac{255}{256}.$

Then

$y=\frac{-\log(255/256)}{2\pi}\approx 6.229164237228602\times 10^{-4}.$

Because \(B_4=-1/30\), the approximation becomes

$\mathcal{S}_3(q)\approx \frac{-1/30}{8}\left(1-y^{-4}\right)=\frac{y^{-4}-1}{240}.$

Numerically this gives

$\mathcal{S}_3(q)\approx 2.767385314772\times 10^{10},$

which matches the checkpoint embedded in the implementation. The same method with weight \(16\) handles the final target, where \(y\) is so small that \(y^{-16}\) dominates the entire value.

How the Code Works

The C++, Python, and Java implementations all use high-precision decimal arithmetic so that the final scientific-notation rounding is stable. They compute \(m=(k+1)/2\), derive the corresponding weight \(2m\), form \(q=1-2^{-p}\), compute \(t=-\log q\), and then obtain \(y=t/(2\pi)\).

After that, the implementation selects the Bernoulli constant for weight \(4\), \(8\), or \(16\), evaluates

$\frac{B_{2m}}{4m}\left(1-y^{-2m}\right),$

and formats the result in normalized scientific notation with 12 digits after the decimal point. The C++ implementation also contains explicit checkpoint assertions: one direct summation for \(\mathcal{S}_1(15/16)\approx 3.872155809243\times 10^2\), and modular checks at \((k,p)=(3,8)\) and \((7,15)\). The Java implementation computes \(-\log(1-2^{-p})\) from its convergent power series, while the C++ and Python implementations call high-precision logarithms directly; despite that implementation detail, all three evaluate the same mathematical formula.

Complexity Analysis

The direct checkpoint is linear in the number of retained terms, because it keeps summing until the next term is negligible. The modular branch performs only a constant number of high-precision operations: one logarithm, one division by \(2\pi\), one large negative power, and a few multiplications and divisions. Memory usage is \(O(1)\). For the target parameters, this turns an impractically slow series evaluation into an effectively constant-size computation.

Footnotes and References

1. Problem page: https://projecteuler.net/problem=722
2. Lambert series: Wikipedia - Lambert series
3. Divisor function: Wikipedia - Divisor function
4. Bernoulli number: Wikipedia - Bernoulli number
5. Eisenstein series: Wikipedia - Eisenstein series