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.
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.
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\).
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}.$$
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.
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.
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.
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.
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.
Für \(0<q<1\) ist die gesuchte Größe die Lambert-Reihe
$$\mathcal{S}_k(q)=\sum_{d\ge 1}\frac{d^k q^d}{1-q^d}.$$
Liegt \(q\) extrem nahe bei \(1\), dann bleibt der Nenner \(1-q^d\) für viele Werte von \(d\) sehr klein, und die direkte Summation konvergiert quälend langsam. Deshalb verwenden die Implementierungen eine direkte Summe nur für einen kleinen Kontrollfall und berechnen den eigentlichen Zielwert über eine Transformation von Eisenstein-Reihen.
Der entscheidende Punkt ist, dass Lambert-Reihen mit ungeradem Exponenten als Fourier-Entwicklungen von Eisenstein-Reihen gerader Ordnung auftreten. Dadurch lässt sich eine langsam konvergierende Reihe nahe \(q=1\) in einen modularen Ausdruck umschreiben, der sofort auswertbar ist.
Wir entwickeln zunächst den geometrischen Nenner:
$$\frac{q^d}{1-q^d}=\sum_{r\ge 1} q^{dr}.$$
Damit folgt
$$\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,$$
wobei
$$\sigma_k(n)=\sum_{d\mid n} d^k$$
die Potenz-Teilerschumme ist. Das Problem besteht also darin, eine \(q\)-Reihe mit Koeffizienten \(\sigma_k(n)\) effizient auszuwerten. Der direkte Kontrollfall in den Implementierungen ist \(k=1\) mit \(q=1-2^{-4}=15/16\).
Für den modularen Zweig verwenden die Implementierungen ungerade \(k\) der Form
$$k=2m-1,$$
wobei die unterstützten Gewichte \(2m\in\{4,8,16\}\) sind. Die normierte Eisenstein-Reihe erfüllt
$$E_{2m}(\tau)=1-\frac{4m}{B_{2m}}\sum_{n\ge 1}\sigma_{2m-1}(n)e^{2\pi i n\tau},$$
wobei \(B_{2m}\) die Bernoulli-Zahl zum Gewicht \(2m\) ist. Mit
$$q=e^{2\pi i\tau}$$
erhält man direkt
$$\mathcal{S}_{2m-1}(q)=\frac{B_{2m}}{4m}\bigl(1-E_{2m}(\tau)\bigr).$$
Genau diese Identität trägt den Schnellweg. Verwendet werden die Konstanten
$$B_4=B_8=-\frac{1}{30},\qquad B_{16}=-\frac{3617}{510}.$$
Die Zielparameter setzen
$$q=1-2^{-p}.$$
Dann definieren wir
$$t=-\log q,\qquad y=\frac{t}{2\pi}.$$
Daraus folgt
$$q=e^{-t}=e^{-2\pi y}=e^{2\pi i(iy)},$$
also wird die Eisenstein-Reihe am rein imaginären Punkt
$$\tau=iy$$
ausgewertet. Für großes \(p\) ist \(2^{-p}\) winzig, also sind auch \(t\) und \(y\) winzig. Genau dann ist die ursprüngliche Lambert-Reihe am schlechtesten direkt zu summieren.
Eisenstein-Reihen erfüllen
$$E_{2m}\!\left(-\frac{1}{\tau}\right)=\tau^{2m}E_{2m}(\tau).$$
Mit \(\tau=iy\) wird daraus
$$E_{2m}(iy)=(iy)^{-2m}E_{2m}(i/y).$$
In den implementierten Fällen gilt \(2m\in\{4,8,16\}\), also ist \(m\) gerade und damit
$$ (iy)^{-2m}=y^{-2m}. $$
Da \(y\) sehr klein ist, besitzt \(i/y\) einen sehr großen Imaginärteil. Die Fourier-Entwicklung von \(E_{2m}(i/y)\) ist daher exponentiell nah an \(1\):
$$E_{2m}(i/y)=1+O\!\left(e^{-2\pi/y}\right).$$
Somit dominiert
$$E_{2m}(iy)\approx y^{-2m},$$
und daraus folgt die im Code verwendete Näherung
$$\mathcal{S}_{2m-1}(q)\approx \frac{B_{2m}}{4m}\left(1-y^{-2m}\right).$$
Hier ist
$$m=\frac{k+1}{2}=2,\qquad q=1-2^{-8}=\frac{255}{256}.$$
Dann erhält man
$$y=\frac{-\log(255/256)}{2\pi}\approx 6.229164237228602\times 10^{-4}.$$
Mit \(B_4=-1/30\) wird die Näherung zu
$$\mathcal{S}_3(q)\approx \frac{-1/30}{8}\left(1-y^{-4}\right)=\frac{y^{-4}-1}{240}.$$
Numerisch ergibt sich
$$\mathcal{S}_3(q)\approx 2.767385314772\times 10^{10},$$
genau der Kontrollwert der Implementierung. Dasselbe Verfahren mit Gewicht \(16\) liefert den Zielwert; dort ist \(y\) so klein, dass \(y^{-16}\) die Größe praktisch vollständig bestimmt.
Die C++-, Python- und Java-Implementierungen arbeiten alle mit hochpräziser Dezimalarithmetik, damit die Rundung in wissenschaftlicher Schreibweise stabil bleibt. Zuerst berechnen sie \(m=(k+1)/2\), daraus das Gewicht \(2m\), dann \(q=1-2^{-p}\), anschließend \(t=-\log q\) und schließlich \(y=t/(2\pi)\).
Danach wählt die Implementierung die Bernoulli-Konstante für Gewicht \(4\), \(8\) oder \(16\), wertet
$$\frac{B_{2m}}{4m}\left(1-y^{-2m}\right)$$
aus und formatiert das Resultat in normierter wissenschaftlicher Schreibweise mit 12 Nachkommastellen in der Mantisse. Die C++-Implementierung enthält zusätzlich Kontrollwerte: eine direkte Summation für \(\mathcal{S}_1(15/16)\approx 3.872155809243\times 10^2\) sowie modulare Tests bei \((k,p)=(3,8)\) und \((7,15)\). Die Java-Implementierung berechnet \(-\log(1-2^{-p})\) über seine konvergente Potenzreihe, während C++ und Python einen hochpräzisen Logarithmus direkt verwenden; mathematisch evaluieren alle drei denselben Ausdruck.
Der direkte Kontrollfall ist linear in der Zahl der aufsummierten Terme, denn es wird summiert, bis der nächste Beitrag vernachlässigbar klein ist. Der modulare Zweig benötigt nur eine konstante Zahl hochpräziser Operationen: einen Logarithmus, eine Division durch \(2\pi\), eine große negative Potenz und einige Multiplikationen und Divisionen. Der Speicherbedarf ist \(O(1)\). Für die Zielparameter wird damit eine praktisch unbrauchbar langsame Reihenberechnung in eine effektiv konstante Auswertung umgewandelt.
\(0<q<1\) için hesaplanan büyüklük şu Lambert tipi seridir:
$$\mathcal{S}_k(q)=\sum_{d\ge 1}\frac{d^k q^d}{1-q^d}.$$
\(q\) değeri \(1\)'e çok yaklaştığında \(1-q^d\) paydası birçok \(d\) için çok küçük kalır ve doğrudan toplama son derece yavaş yakınsar. Bu yüzden implementasyon küçük bir kontrol örneğinde doğrudan toplam alır, asıl büyük hedefte ise Eisenstein serileri üzerinden gelen modüler dönüşümü kullanır.
Temel fikir şudur: tek kuvvetli Lambert serileri, çift ağırlıklı Eisenstein serilerinin Fourier açılımı olarak yazılabilir. Böylece \(q\to 1\) bölgesindeki yavaş seri, modüler tarafta çok daha hızlı hesaplanan bir ifadeye dönüşür.
Önce geometrik payda açılır:
$$\frac{q^d}{1-q^d}=\sum_{r\ge 1} q^{dr}.$$
Bunu yerine koyunca
$$\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$$
elde edilir. Burada
$$\sigma_k(n)=\sum_{d\mid n} d^k$$
\(n\)'in bölenlerinin \(k\). kuvvetlerinin toplamıdır. Yani sorun aslında \(\sigma_k(n)\) katsayılarına sahip bir \(q\)-serisini hesaplamaktır. Implementasyonun doğrudan kontrol ettiği küçük örnek \(k=1\) ve \(q=1-2^{-4}=15/16\) durumudur.
Modüler kolda kullanılan üsler
$$k=2m-1$$
şeklindedir ve desteklenen ağırlıklar \(2m\in\{4,8,16\}\) olur. Normalize Eisenstein serisi
$$E_{2m}(\tau)=1-\frac{4m}{B_{2m}}\sum_{n\ge 1}\sigma_{2m-1}(n)e^{2\pi i n\tau}$$
eşitliğini sağlar. Eğer
$$q=e^{2\pi i\tau}$$
yazarsak
$$\mathcal{S}_{2m-1}(q)=\frac{B_{2m}}{4m}\bigl(1-E_{2m}(\tau)\bigr)$$
sonucuna ulaşırız. Kısayolun temeli budur. İmplementasyonda kullanılan Bernoulli sayıları da
$$B_4=B_8=-\frac{1}{30},\qquad B_{16}=-\frac{3617}{510}$$
şeklindedir.
Hedef girdilerde
$$q=1-2^{-p}$$
alınır. Şimdi
$$t=-\log q,\qquad y=\frac{t}{2\pi}$$
tanımını yapalım. Böylece
$$q=e^{-t}=e^{-2\pi y}=e^{2\pi i(iy)}$$
olur; yani Eisenstein serisini
$$\tau=iy$$
noktasında değerlendiriyoruz. \(p\) büyüdükçe \(2^{-p}\) çok küçülür; dolayısıyla hem \(t\) hem de \(y\) çok küçülür. İşte bu bölge, özgün Lambert serisinin terim terim toplanmasının en zor olduğu bölgedir.
Eisenstein serileri şu modüler bağıntıyı sağlar:
$$E_{2m}\!\left(-\frac{1}{\tau}\right)=\tau^{2m}E_{2m}(\tau).$$
\(\tau=iy\) koyarsak
$$E_{2m}(iy)=(iy)^{-2m}E_{2m}(i/y)$$
elde edilir. Kullanılan ağırlıklar \(4,8,16\) olduğundan \(m\) çifttir ve bu yüzden
$$ (iy)^{-2m}=y^{-2m}. $$
Ayrıca \(y\) çok küçük olduğundan \(i/y\) büyük sanal kısma sahiptir. Bu da Fourier açılımında
$$E_{2m}(i/y)=1+O\!\left(e^{-2\pi/y}\right)$$
yaklaşımını verir. Sonuçta baskın terim
$$E_{2m}(iy)\approx y^{-2m}$$
olur ve Lambert serisi
$$\mathcal{S}_{2m-1}(q)\approx \frac{B_{2m}}{4m}\left(1-y^{-2m}\right)$$
biçimine iner. Kodun modüler kolunda hesaplanan ifade tam olarak budur.
Bu örnekte
$$m=\frac{k+1}{2}=2,\qquad q=1-2^{-8}=\frac{255}{256}.$$
Buna göre
$$y=\frac{-\log(255/256)}{2\pi}\approx 6.229164237228602\times 10^{-4}.$$
\(B_4=-1/30\) olduğundan
$$\mathcal{S}_3(q)\approx \frac{-1/30}{8}\left(1-y^{-4}\right)=\frac{y^{-4}-1}{240}$$
yazılır. Sayısal değer
$$\mathcal{S}_3(q)\approx 2.767385314772\times 10^{10}$$
çıkar ve bu, implementasyondaki kontrol değeriyle birebir örtüşür. Aynı fikir ağırlık \(16\) için de kullanılır; orada \(y\) o kadar küçüktür ki \(y^{-16}\) toplam değerin baskın kısmını belirler.
C++, Python ve Java implementasyonlarının ortak noktası yüksek hassasiyetli ondalık aritmetik kullanmalarıdır; böylece bilimsel gösterimde son yuvarlama kararlı kalır. Önce \(m=(k+1)/2\) ve buna karşılık gelen ağırlık \(2m\) hesaplanır. Ardından \(q=1-2^{-p}\), \(t=-\log q\) ve \(y=t/(2\pi)\) bulunur.
Daha sonra implementasyon ağırlık \(4\), \(8\) veya \(16\) için uygun Bernoulli sabitini seçer,
$$\frac{B_{2m}}{4m}\left(1-y^{-2m}\right)$$
ifadesini hesaplar ve sonucu mantisada ondalık noktadan sonra 12 basamak olacak şekilde normalize edilmiş bilimsel gösterime dönüştürür. C++ implementasyonu ayrıca iki tür kontrol içerir: \(\mathcal{S}_1(15/16)\approx 3.872155809243\times 10^2\) için doğrudan toplam ve \((k,p)=(3,8)\) ile \((7,15)\) için modüler doğrulamalar. Java tarafı \(-\log(1-2^{-p})\) değerini yakınsak log serisinden üretir; C++ ve Python ise yüksek hassasiyetli logaritmayı doğrudan çağırır. Yine de üçü de aynı matematiksel formülü uygular.
Doğrudan kontrol örneği, toplanan terim sayısı kadar doğrusal zamana sahiptir; çünkü bir sonraki terim ihmal edilebilir küçüklüğe düşene kadar toplam sürdürülür. Modüler kol ise sabit sayıda yüksek hassasiyetli işlem yapar: bir logaritma, \(2\pi\)'ye bölüm, büyük bir negatif üs ve birkaç çarpma-bölme. Bellek kullanımı \(O(1)\)'dir. Hedef parametrelerde bu yaklaşım, pratikte kullanılamayacak kadar yavaş bir seri toplamını etkin olarak sabit boyutlu bir hesaba dönüştürür.
Para \(0<q<1\), la cantidad que se desea calcular es la serie de tipo Lambert
$$\mathcal{S}_k(q)=\sum_{d\ge 1}\frac{d^k q^d}{1-q^d}.$$
Cuando \(q\) está extremadamente cerca de \(1\), el denominador \(1-q^d\) sigue siendo muy pequeño para muchos valores de \(d\), así que la suma directa converge con mucha lentitud. Por eso las implementaciones usan suma directa solo en una comprobación pequeña y reservan una transformación con series de Eisenstein para el caso grande.
La observación central es que las series de Lambert con potencia impar aparecen como desarrollos de Fourier de series de Eisenstein de peso par. Eso permite sustituir una suma lenta cerca de \(q=1\) por una expresión modular mucho más manejable.
Primero se expande el denominador geométrico:
$$\frac{q^d}{1-q^d}=\sum_{r\ge 1} q^{dr}.$$
Entonces
$$\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,$$
donde
$$\sigma_k(n)=\sum_{d\mid n} d^k$$
es la suma de potencias de los divisores. Así, el problema se convierte en evaluar una \(q\)-serie cuyos coeficientes son \(\sigma_k(n)\). La comprobación directa incluida en las implementaciones corresponde a \(k=1\) y \(q=1-2^{-4}=15/16\).
En la rama modular, las implementaciones usan exponentes impares de la forma
$$k=2m-1,$$
con pesos admitidos \(2m\in\{4,8,16\}\). La serie de Eisenstein normalizada satisface
$$E_{2m}(\tau)=1-\frac{4m}{B_{2m}}\sum_{n\ge 1}\sigma_{2m-1}(n)e^{2\pi i n\tau},$$
donde \(B_{2m}\) es el número de Bernoulli de peso \(2m\). Si escribimos
$$q=e^{2\pi i\tau},$$
entonces
$$\mathcal{S}_{2m-1}(q)=\frac{B_{2m}}{4m}\bigl(1-E_{2m}(\tau)\bigr).$$
Esa es la identidad exacta que usa el atajo. Los valores de Bernoulli incorporados son
$$B_4=B_8=-\frac{1}{30},\qquad B_{16}=-\frac{3617}{510}.$$
Los parámetros objetivo toman
$$q=1-2^{-p}.$$
Definimos
$$t=-\log q,\qquad y=\frac{t}{2\pi}.$$
Con esto
$$q=e^{-t}=e^{-2\pi y}=e^{2\pi i(iy)},$$
de modo que la serie de Eisenstein se evalúa en
$$\tau=iy.$$
Si \(p\) es grande, entonces \(2^{-p}\), \(t\) y \(y\) son muy pequeños. Precisamente en ese régimen la suma original converge peor.
Las series de Eisenstein cumplen
$$E_{2m}\!\left(-\frac{1}{\tau}\right)=\tau^{2m}E_{2m}(\tau).$$
Al sustituir \(\tau=iy\) se obtiene
$$E_{2m}(iy)=(iy)^{-2m}E_{2m}(i/y).$$
En los pesos implementados \(4,8,16\), el valor de \(m\) es par, por lo que
$$ (iy)^{-2m}=y^{-2m}. $$
Además, \(y\) es muy pequeño, así que \(i/y\) tiene parte imaginaria enorme y la expansión de Fourier queda exponencialmente cerca de \(1\):
$$E_{2m}(i/y)=1+O\!\left(e^{-2\pi/y}\right).$$
Por tanto, domina el término
$$E_{2m}(iy)\approx y^{-2m},$$
y la serie de Lambert queda aproximada por
$$\mathcal{S}_{2m-1}(q)\approx \frac{B_{2m}}{4m}\left(1-y^{-2m}\right).$$
Esa es exactamente la fórmula que evalúa la implementación.
Aquí
$$m=\frac{k+1}{2}=2,\qquad q=1-2^{-8}=\frac{255}{256}.$$
Entonces
$$y=\frac{-\log(255/256)}{2\pi}\approx 6.229164237228602\times 10^{-4}.$$
Como \(B_4=-1/30\), la aproximación se convierte en
$$\mathcal{S}_3(q)\approx \frac{-1/30}{8}\left(1-y^{-4}\right)=\frac{y^{-4}-1}{240}.$$
Numéricamente,
$$\mathcal{S}_3(q)\approx 2.767385314772\times 10^{10},$$
que coincide con la comprobación incorporada en la implementación. El mismo procedimiento con peso \(16\) resuelve el caso final, donde \(y\) es tan pequeño que \(y^{-16}\) domina por completo.
Las implementaciones en C++, Python y Java usan aritmética decimal de alta precisión para que el redondeo final en notación científica sea estable. Primero calculan \(m=(k+1)/2\), luego el peso \(2m\), después forman \(q=1-2^{-p}\), obtienen \(t=-\log q\) y finalmente \(y=t/(2\pi)\).
Después, la implementación selecciona la constante de Bernoulli correspondiente al peso \(4\), \(8\) o \(16\), evalúa
$$\frac{B_{2m}}{4m}\left(1-y^{-2m}\right),$$
y normaliza la salida en notación científica con 12 cifras después del punto decimal en la mantisa. La implementación en C++ además contiene verificaciones explícitas: una suma directa para \(\mathcal{S}_1(15/16)\approx 3.872155809243\times 10^2\), y comprobaciones modulares para \((k,p)=(3,8)\) y \((7,15)\). La implementación en Java obtiene \(-\log(1-2^{-p})\) mediante su serie de potencias convergente, mientras que C++ y Python llaman un logaritmo de alta precisión; aun así, las tres siguen la misma fórmula matemática.
La comprobación directa es lineal en el número de términos retenidos, porque sigue sumando hasta que el siguiente término sea despreciable. La rama modular realiza solo una cantidad constante de operaciones de alta precisión: un logaritmo, una división por \(2\pi\), una potencia negativa grande y unas pocas multiplicaciones y divisiones. El uso de memoria es \(O(1)\). Para los parámetros del objetivo, esto reemplaza una suma impracticablemente lenta por un cálculo de tamaño esencialmente constante.
当 \(0<q<1\) 时，本题需要计算的对象可以写成一个 Lambert 型级数：
$$\mathcal{S}_k(q)=\sum_{d\ge 1}\frac{d^k q^d}{1-q^d}.$$
如果 \(q\) 非常接近 \(1\)，那么很多 \(d\) 对应的分母 \(1-q^d\) 都会非常小，逐项直接求和的收敛速度就会非常差。实现因此只把直接求和留给一个小型校验点，而把真正的大参数计算改写成 Eisenstein 级数的模变换问题。
核心事实是：奇次幂的 Lambert 级数恰好可以嵌入偶权 Eisenstein 级数的 Fourier 展开。这样一来，原本在 \(q\to 1\) 时很难直接求值的级数，就能转化成一个在模形式一侧几乎立刻可算的表达式。
先把几何分母展开：
$$\frac{q^d}{1-q^d}=\sum_{r\ge 1} q^{dr}.$$
代回原式后得到
$$\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,$$
其中
$$\sigma_k(n)=\sum_{d\mid n} d^k$$
表示 \(n\) 的所有正约数的 \(k\) 次幂之和。也就是说，程序真正面对的是一个以 \(\sigma_k(n)\) 为系数的 \(q\)-级数。实现里保留的直接校验例子是 \(k=1\) 且 \(q=1-2^{-4}=15/16\)。
在模方法分支里，代码处理的指数满足
$$k=2m-1,$$
并且支持的权重只有 \(2m\in\{4,8,16\}\)。规范化的 Eisenstein 级数满足
$$E_{2m}(\tau)=1-\frac{4m}{B_{2m}}\sum_{n\ge 1}\sigma_{2m-1}(n)e^{2\pi i n\tau},$$
其中 \(B_{2m}\) 是对应权重的 Bernoulli 数。如果令
$$q=e^{2\pi i\tau},$$
那么马上得到
$$\mathcal{S}_{2m-1}(q)=\frac{B_{2m}}{4m}\bigl(1-E_{2m}(\tau)\bigr).$$
这就是实现所依赖的精确恒等式。代码中写死的 Bernoulli 常数为
$$B_4=B_8=-\frac{1}{30},\qquad B_{16}=-\frac{3617}{510}.$$
目标输入采用
$$q=1-2^{-p}.$$
接着定义
$$t=-\log q,\qquad y=\frac{t}{2\pi}.$$
于是
$$q=e^{-t}=e^{-2\pi y}=e^{2\pi i(iy)},$$
因此 Eisenstein 级数需要在纯虚点
$$\tau=iy$$
处求值。当 \(p\) 很大时，\(2^{-p}\) 极小，所以 \(t\) 和 \(y\) 也极小；这正是原始 Lambert 级数最难直接求和的区域。
Eisenstein 级数满足模变换公式
$$E_{2m}\!\left(-\frac{1}{\tau}\right)=\tau^{2m}E_{2m}(\tau).$$
将 \(\tau=iy\) 代入，可得
$$E_{2m}(iy)=(iy)^{-2m}E_{2m}(i/y).$$
对本实现真正使用的权重 \(4,8,16\) 来说，\(m\) 都是偶数，因此
$$ (iy)^{-2m}=y^{-2m}. $$
另一方面，\(y\) 很小时，\(i/y\) 的虚部就非常大，于是 Fourier 展开里的指数项几乎全部被压扁，得到
$$E_{2m}(i/y)=1+O\!\left(e^{-2\pi/y}\right).$$
所以主导行为就是
$$E_{2m}(iy)\approx y^{-2m},$$
从而 Lambert 级数可近似为
$$\mathcal{S}_{2m-1}(q)\approx \frac{B_{2m}}{4m}\left(1-y^{-2m}\right).$$
这正是三份实现共同计算的公式。
此时
$$m=\frac{k+1}{2}=2,\qquad q=1-2^{-8}=\frac{255}{256}.$$
于是
$$y=\frac{-\log(255/256)}{2\pi}\approx 6.229164237228602\times 10^{-4}.$$
因为 \(B_4=-1/30\)，所以近似式化为
$$\mathcal{S}_3(q)\approx \frac{-1/30}{8}\left(1-y^{-4}\right)=\frac{y^{-4}-1}{240}.$$
数值上得到
$$\mathcal{S}_3(q)\approx 2.767385314772\times 10^{10},$$
与实现中的校验值完全一致。最终目标对应的权重是 \(16\)，那时 \(y\) 更小，\(y^{-16}\) 会成为绝对主导项，这也解释了结果数量级为何会非常大。
C++、Python 和 Java 实现都使用高精度十进制算术，以保证最后输出科学计数法时的舍入稳定。它们先计算 \(m=(k+1)/2\) 以及对应的权重 \(2m\)，然后构造 \(q=1-2^{-p}\)，求出 \(t=-\log q\)，再得到 \(y=t/(2\pi)\)。
接下来，程序根据权重 \(4\)、\(8\) 或 \(16\) 选取相应的 Bernoulli 常数，计算
$$\frac{B_{2m}}{4m}\left(1-y^{-2m}\right),$$
最后把结果整理成尾数保留 12 位小数的标准科学计数法。C++ 实现还额外放入了若干断言型校验：一是对 \(\mathcal{S}_1(15/16)\approx 3.872155809243\times 10^2\) 进行直接求和，二是检查模公式在 \((k,p)=(3,8)\) 和 \((7,15)\) 上的输出。Java 实现因为缺少内建高精度自然对数，所以通过其收敛的幂级数来计算 \(-\log(1-2^{-p})\)；而 C++ 与 Python 直接调用高精度对数。尽管实现细节不同，它们求值的数学表达式完全一致。
直接校验分支的时间复杂度与保留的项数成线性关系，因为它会一直累加到下一项足够小为止。模方法分支只做常数个高精度操作：一次对数、一次除以 \(2\pi\)、一次高次负幂，以及少量乘除法，空间复杂度为 \(O(1)\)。对目标参数而言，这相当于把一个几乎无法直接求和的慢收敛级数，改写成一个规模几乎固定的高精度算式。
При \(0<q<1\) требуется вычислить ламбертов ряд
$$\mathcal{S}_k(q)=\sum_{d\ge 1}\frac{d^k q^d}{1-q^d}.$$
Когда \(q\) очень близко к \(1\), знаменатель \(1-q^d\) остается малым для многих \(d\), поэтому прямое суммирование сходится крайне медленно. Поэтому реализации используют прямой подсчет только для небольшого контрольного примера, а основной случай сводят к преобразованию через ряды Эйзенштейна.
Главная идея состоит в том, что ламбертовы ряды с нечетной степенью входят в Fourier-разложения рядов Эйзенштейна четного веса. Благодаря этому медленно сходящийся ряд в режиме \(q\to 1\) можно заменить модульным выражением, которое вычисляется очень быстро.
Сначала раскрываем геометрический знаменатель:
$$\frac{q^d}{1-q^d}=\sum_{r\ge 1} q^{dr}.$$
Тогда
$$\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,$$
где
$$\sigma_k(n)=\sum_{d\mid n} d^k$$
есть сумма \(k\)-х степеней всех положительных делителей числа \(n\). Значит, задача на самом деле состоит в вычислении \(q\)-ряда с коэффициентами \(\sigma_k(n)\). Прямой контрольный пример в реализации соответствует \(k=1\) и \(q=1-2^{-4}=15/16\).
В модульной ветви реализации используются нечетные \(k\) вида
$$k=2m-1,$$
причем поддерживаемые веса равны \(2m\in\{4,8,16\}\). Нормированный ряд Эйзенштейна удовлетворяет формуле
$$E_{2m}(\tau)=1-\frac{4m}{B_{2m}}\sum_{n\ge 1}\sigma_{2m-1}(n)e^{2\pi i n\tau},$$
где \(B_{2m}\) обозначает число Бернулли соответствующего веса. Если положить
$$q=e^{2\pi i\tau},$$
то сразу получаем
$$\mathcal{S}_{2m-1}(q)=\frac{B_{2m}}{4m}\bigl(1-E_{2m}(\tau)\bigr).$$
Именно это точное тождество и используется как основа ускорения. В коде зашиты значения
$$B_4=B_8=-\frac{1}{30},\qquad B_{16}=-\frac{3617}{510}.$$
В целевых входах берется
$$q=1-2^{-p}.$$
Далее вводятся величины
$$t=-\log q,\qquad y=\frac{t}{2\pi}.$$
Тогда
$$q=e^{-t}=e^{-2\pi y}=e^{2\pi i(iy)},$$
то есть ряд Эйзенштейна нужно рассматривать в точке
$$\tau=iy.$$
Когда \(p\) велико, число \(2^{-p}\) мало, а вместе с ним малы и \(t\), и \(y\). Именно в этой зоне исходный ламбертов ряд особенно плохо поддается прямому суммированию.
Для рядов Эйзенштейна верно
$$E_{2m}\!\left(-\frac{1}{\tau}\right)=\tau^{2m}E_{2m}(\tau).$$
После подстановки \(\tau=iy\) получаем
$$E_{2m}(iy)=(iy)^{-2m}E_{2m}(i/y).$$
В используемых здесь весах \(4,8,16\) параметр \(m\) четный, поэтому
$$ (iy)^{-2m}=y^{-2m}. $$
Так как \(y\) очень мало, число \(i/y\) имеет огромную мнимую часть, а Fourier-разложение становится экспоненциально близким к \(1\):
$$E_{2m}(i/y)=1+O\!\left(e^{-2\pi/y}\right).$$
Следовательно, главный вклад равен
$$E_{2m}(iy)\approx y^{-2m},$$
и ламбертов ряд приближается выражением
$$\mathcal{S}_{2m-1}(q)\approx \frac{B_{2m}}{4m}\left(1-y^{-2m}\right).$$
Это именно та формула, которую вычисляет программа.
В этом случае
$$m=\frac{k+1}{2}=2,\qquad q=1-2^{-8}=\frac{255}{256}.$$
Тогда
$$y=\frac{-\log(255/256)}{2\pi}\approx 6.229164237228602\times 10^{-4}.$$
Поскольку \(B_4=-1/30\), имеем
$$\mathcal{S}_3(q)\approx \frac{-1/30}{8}\left(1-y^{-4}\right)=\frac{y^{-4}-1}{240}.$$
Численно это дает
$$\mathcal{S}_3(q)\approx 2.767385314772\times 10^{10},$$
что точно совпадает с контрольным значением в реализации. Тот же прием при весе \(16\) решает финальный случай; там \(y\) настолько мало, что величина \(y^{-16}\) полностью доминирует.
Реализации на C++, Python и Java используют высокоточную десятичную арифметику, чтобы итоговое округление в научной записи было устойчивым. Сначала вычисляются \(m=(k+1)/2\) и вес \(2m\), затем формируется \(q=1-2^{-p}\), находится \(t=-\log q\), после чего получается \(y=t/(2\pi)\).
Затем реализация выбирает константу Бернулли для веса \(4\), \(8\) или \(16\), вычисляет
$$\frac{B_{2m}}{4m}\left(1-y^{-2m}\right),$$
и приводит результат к нормализованной научной записи с 12 знаками после запятой в мантиссе. C++-реализация также содержит явные проверки: прямое суммирование для \(\mathcal{S}_1(15/16)\approx 3.872155809243\times 10^2\) и модульные проверки для \((k,p)=(3,8)\) и \((7,15)\). В Java величина \(-\log(1-2^{-p})\) вычисляется по сходящемуся степенному ряду, тогда как C++ и Python используют высокоточный логарифм напрямую; несмотря на это различие, математическая формула у всех трех одна и та же.
Прямой контрольный подсчет имеет линейную сложность по числу сохраненных членов, поскольку суммирование продолжается, пока следующий вклад не станет пренебрежимо малым. Модульная ветвь использует лишь постоянное число высокоточных операций: один логарифм, одно деление на \(2\pi\), одну большую отрицательную степень и несколько умножений и делений. Память равна \(O(1)\). Для целевых параметров это заменяет практически бесполезный медленно сходящийся ряд вычислением почти постоянного размера.
عندما يكون \(0<q<1\)، فإن الكمية المطلوب حسابها هي سلسلة من نوع لامبرت:
$$\mathcal{S}_k(q)=\sum_{d\ge 1}\frac{d^k q^d}{1-q^d}.$$
إذا كان \(q\) قريبًا جدًا من \(1\)، فإن المقام \(1-q^d\) يبقى صغيرًا لكثير من القيم \(d\)، ولهذا يصبح الجمع المباشر بطيئ التقارب جدًا. لذلك تستعمل التنفيذات الجمع المباشر فقط في نقطة تحقق صغيرة، ثم تنتقل في الحالة الكبيرة إلى تحويل يعتمد على سلاسل آيزنشتاين.
الفكرة الحاسمة هي أن سلاسل لامبرت ذات الأس الفردي تظهر داخل توسعات Fourier لسلاسل آيزنشتاين ذات الوزن الزوجي. وبذلك يمكن تحويل سلسلة بطيئة التقارب قرب \(q=1\) إلى صيغة معيارية يمكن تقييمها بسرعة كبيرة.
نبدأ بتوسيع المقام الهندسي:
$$\frac{q^d}{1-q^d}=\sum_{r\ge 1} q^{dr}.$$
وعند التعويض في الصيغة الأصلية نحصل على
$$\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,$$
حيث
$$\sigma_k(n)=\sum_{d\mid n} d^k$$
هي دالة مجموع القوى على قواسم \(n\). إذًا المسألة تتحول إلى تقييم متسلسلة \(q\) معاملاتُها هي \(\sigma_k(n)\). ونقطة التحقق المباشرة الموجودة في التنفيذ هي الحالة \(k=1\) مع \(q=1-2^{-4}=15/16\).
في الفرع المعياري تستعمل التنفيذات قيَمًا فردية من الشكل
$$k=2m-1,$$
مع الأوزان المدعومة \(2m\in\{4,8,16\}\). وتحقق سلسلة آيزنشتاين المعيارية المطَبَّعة العلاقة
$$E_{2m}(\tau)=1-\frac{4m}{B_{2m}}\sum_{n\ge 1}\sigma_{2m-1}(n)e^{2\pi i n\tau},$$
حيث \(B_{2m}\) هو عدد برنولي الموافق للوزن \(2m\). وإذا كتبنا
$$q=e^{2\pi i\tau},$$
فإننا نحصل مباشرة على
$$\mathcal{S}_{2m-1}(q)=\frac{B_{2m}}{4m}\bigl(1-E_{2m}(\tau)\bigr).$$
هذه هي المطابقة الدقيقة التي يعتمد عليها الاختصار. أما الثوابت المستعملة فعليًا فهي
$$B_4=B_8=-\frac{1}{30},\qquad B_{16}=-\frac{3617}{510}.$$
في القيم الهدفية نأخذ
$$q=1-2^{-p}.$$
ثم نعرّف
$$t=-\log q,\qquad y=\frac{t}{2\pi}.$$
وعندها يصبح
$$q=e^{-t}=e^{-2\pi y}=e^{2\pi i(iy)},$$
أي إن تقييم سلسلة آيزنشتاين يتم عند النقطة التخيلية الصرفة
$$\tau=iy.$$
عندما يكون \(p\) كبيرًا، تكون \(2^{-p}\) صغيرة جدًا، وبالتالي يكون كل من \(t\) و\(y\) صغيرين جدًا أيضًا. وهذه هي بالضبط المنطقة التي تصبح فيها السلسلة الأصلية صعبة للغاية إذا جُمعت حدًا حدًا.
سلاسل آيزنشتاين تحقق العلاقة المعيارية
$$E_{2m}\!\left(-\frac{1}{\tau}\right)=\tau^{2m}E_{2m}(\tau).$$
وبالتعويض عن \(\tau=iy\) نحصل على
$$E_{2m}(iy)=(iy)^{-2m}E_{2m}(i/y).$$
وفي الأوزان المستخدمة هنا \(4,8,16\)، يكون \(m\) زوجيًا، ولذلك
$$ (iy)^{-2m}=y^{-2m}. $$
ولأن \(y\) صغير جدًا، فإن \(i/y\) يملك جزءًا تخيليًا هائلًا، فتغدو توسعة Fourier قريبة من \(1\) اقترابًا أسيًا:
$$E_{2m}(i/y)=1+O\!\left(e^{-2\pi/y}\right).$$
ومن ثم يكون الحد الغالب
$$E_{2m}(iy)\approx y^{-2m},$$
فتأخذ سلسلة لامبرت الصورة التقريبية
$$\mathcal{S}_{2m-1}(q)\approx \frac{B_{2m}}{4m}\left(1-y^{-2m}\right).$$
وهذه هي الصيغة التي تحسبها التنفيذات فعلًا.
في هذا المثال
$$m=\frac{k+1}{2}=2,\qquad q=1-2^{-8}=\frac{255}{256}.$$
ومن ثم
$$y=\frac{-\log(255/256)}{2\pi}\approx 6.229164237228602\times 10^{-4}.$$
وبما أن \(B_4=-1/30\)، فإن الصيغة تصبح
$$\mathcal{S}_3(q)\approx \frac{-1/30}{8}\left(1-y^{-4}\right)=\frac{y^{-4}-1}{240}.$$
وقيمتها العددية
$$\mathcal{S}_3(q)\approx 2.767385314772\times 10^{10},$$
وهو نفس مقدار التحقق الموجود داخل التنفيذ. والفكرة نفسها مع الوزن \(16\) تعالج الحالة النهائية، حيث يصبح \(y^{-16}\) هو الحد المسيطر بالكامل.
تنفيذات C++ وPython وJava كلها تستعمل حسابًا عشريًا عالي الدقة حتى يبقى التقريب النهائي في الصيغة العلمية مستقرًا. أولًا يُحسب \(m=(k+1)/2\) ثم الوزن الموافق \(2m\)، وبعد ذلك تُكوَّن القيمة \(q=1-2^{-p}\)، ثم يُحسب \(t=-\log q\)، وأخيرًا \(y=t/(2\pi)\).
بعد هذا تختار الشيفرة ثابت برنولي الموافق للوزن \(4\) أو \(8\) أو \(16\)، ثم تقيّم
$$\frac{B_{2m}}{4m}\left(1-y^{-2m}\right),$$
ثم تنسّق الناتج في صيغة علمية معيارية مع 12 رقمًا بعد الفاصلة العشرية في المانتيسا. كما أن تنفيذ C++ يضيف اختبارات تحقق واضحة: جمعًا مباشرًا لـ \(\mathcal{S}_1(15/16)\approx 3.872155809243\times 10^2\)، إضافة إلى اختبارات معيارية عند \((k,p)=(3,8)\) و\((7,15)\). أما تنفيذ Java فيحسب \(-\log(1-2^{-p})\) من خلال متسلسلته القوية المتقاربة، بينما يستعمل C++ وPython لوغاريتمًا عالي الدقة مباشرة. ومع ذلك فإن التنفيذات الثلاثة تطبق الصيغة الرياضية نفسها تمامًا.
فرع التحقق المباشر خطي في عدد الحدود المجمعة، لأنه يستمر حتى يصبح الحد التالي مهمَلًا عمليًا. أما الفرع المعياري فيستخدم عددًا ثابتًا فقط من العمليات عالية الدقة: لوغاريتمًا واحدًا، وقسمة على \(2\pi\)، وقوة سالبة كبيرة، وبعض الضرب والقسمة. استهلاك الذاكرة هو \(O(1)\). بالنسبة إلى القيم الهدفية، فإن هذا يحول سلسلة بطيئة التقارب جدًا إلى حساب ذي حجم ثابت تقريبًا.