We must evaluate
$$S(N)=\sum_{n=1}^{N} f(n)^2 \pmod{10^9+7},$$
where
$$f(1)=1,\qquad f(2m)=2f(m),\qquad f(2m+1)=2m+1+2f(m)+\frac{f(m)}{m}\quad (m\ge 1).$$
A direct recurrence evaluation up to \(N=10^{16}\) is impossible. The key is to identify a closed form for \(f(n)\), then sum that closed form with a binary prefix method instead of iterating over all integers.
Write \(w(n)=\operatorname{popcount}(n)\), the number of ones in the binary expansion of \(n\). The whole solution is driven by the identity \(f(n)=n\,w(n)\).
The binary weight satisfies
$$w(2m)=w(m),\qquad w(2m+1)=w(m)+1.$$
Assume inductively that \(f(m)=m\,w(m)\). Then for even arguments,
$$f(2m)=2f(m)=2m\,w(m)=2m\,w(2m).$$
For odd arguments, the division term becomes harmless because \(f(m)/m=w(m)\):
$$\begin{aligned} f(2m+1)&=2m+1+2f(m)+\frac{f(m)}{m}\\ &=2m+1+2m\,w(m)+w(m)\\ &=(2m+1)(w(m)+1)\\ &=(2m+1)w(2m+1). \end{aligned}$$
Since \(f(1)=1=1\cdot w(1)\), induction yields
$$f(n)=n\,w(n).$$
Therefore the required sum is simply
$$S(N)=\sum_{n=1}^{N} n^2 w(n)^2.$$
For each \(m\ge 0\), let \(X_m=\{0,1,\dots,2^m-1\}\). For \(a,b\in\{0,1,2\}\), define
$$T_m^{a,b}=\sum_{x\in X_m} x^a w(x)^b.$$
These nine moments are exactly what we need, because any expansion of a square in the value and a square in the popcount can only produce powers \(x^a w(x)^b\) with \(a,b\le 2\).
The target quantity over a complete \(m\)-bit block is the special case
$$T_m^{2,2}=\sum_{x=0}^{2^m-1} x^2 w(x)^2.$$
The base case is immediate:
$$T_0^{0,0}=1,\qquad T_0^{a,b}=0\ \text{for}\ (a,b)\ne(0,0),$$
because the only \(0\)-bit suffix is \(x=0\).
Every \((m+1)\)-bit number is either \(2x\) or \(2x+1\) with \(x\in X_m\). Their popcounts satisfy
$$w(2x)=w(x),\qquad w(2x+1)=w(x)+1.$$
So each new moment obeys
$$T_{m+1}^{a,b}=\sum_{x\in X_m}(2x)^a w(x)^b+\sum_{x\in X_m}(2x+1)^a (w(x)+1)^b.$$
Because \(a,b\le 2\), the right-hand side always reduces to a linear combination of the same nine moments. The smallest updates are
$$T_{m+1}^{0,0}=2T_m^{0,0},\qquad T_{m+1}^{0,1}=2T_m^{0,1}+T_m^{0,0},$$
while the target second-order moment expands to
$$\begin{aligned} T_{m+1}^{2,2}={}&8T_m^{2,2}+8T_m^{2,1}+4T_m^{1,2}+8T_m^{1,1}\\ &+4T_m^{2,0}+4T_m^{1,0}+T_m^{0,2}+2T_m^{0,1}+T_m^{0,0}. \end{aligned}$$
This is the reason the implementations maintain exactly nine aggregate tables and update them level by level.
During the final summation, the higher bits of a number are fixed first. Suppose the already chosen prefix contributes a value \(p\) and contains \(r\) ones. If \(m\) lower bits are still free, every number in that block has the form
$$n=p+x,\qquad 0\le x\lt 2^m,$$
and its popcount is
$$w(n)=r+w(x),$$
because the suffix occupies positions strictly below the fixed prefix.
Hence the whole block contributes
$$B(p,r,m)=\sum_{x=0}^{2^m-1}(p+x)^2(r+w(x))^2.$$
Expanding both squares gives a fixed linear combination of the precomputed moments:
$$\begin{aligned} B(p,r,m)={}&p^2r^2T_m^{0,0}+2p^2rT_m^{0,1}+p^2T_m^{0,2}\\ &+2pr^2T_m^{1,0}+4prT_m^{1,1}+2pT_m^{1,2}\\ &+r^2T_m^{2,0}+2rT_m^{2,1}+T_m^{2,2}. \end{aligned}$$
So once the nine moments are known, an entire interval of length \(2^m\) is summed in constant time.
Write \(N\) in binary and scan from the most significant bit down to the least significant bit. Whenever the current bit of \(N\) is \(1\), there is a complete block of smaller numbers obtained by keeping all earlier bits equal to the current prefix, setting the current bit to \(0\), and letting all lower bits vary freely. That full block is added with \(B(p,r,m)\).
After adding the block, the scan turns the current bit on in the prefix and increases the prefix popcount by one. When the loop ends, the prefix itself equals \(N\), so we add the single endpoint term
$$N^2 w(N)^2.$$
Every integer from \(0\) to \(N\) appears exactly once in this decomposition, and the term for \(0\) is automatically zero.
Since \(13=1101_2\), the scan splits the range into four parts:
$$[0,7],\qquad [8,11],\qquad [12,12],\qquad \{13\}.$$
These four pieces correspond to
$$B(0,0,3),\qquad B(8,1,2),\qquad B(12,2,0),\qquad 13^2\cdot 3^2.$$
Their values are
$$742,\qquad 1877,\qquad 576,\qquad 1521,$$
so
$$S(13)=742+1877+576+1521=4716.$$
As another small checkpoint, the same closed form gives
$$S(10)=\sum_{n=1}^{10} n^2 w(n)^2=1530.$$
The C++, Python, and Java implementations follow the same pipeline. First they precompute the nine tables \(T_m^{a,b}\) for every bit length up to 64. In parallel, they precompute the powers of two modulo \(10^9+7\), so a fixed prefix value can be updated immediately when the scan accepts a new set bit.
Next they scan the bits of \(N\) from high to low. At each set bit, the implementation adds the complete suffix block determined by the current prefix, using the expanded formula for \(B(p,r,m)\). Then it extends the prefix by turning that bit on and increasing the prefix popcount. After all bits have been processed, it adds the single endpoint term \(N^2 w(N)^2\).
All arithmetic is performed modulo \(10^9+7\), and the algorithm never iterates through the integers \(1,2,\dots,N\) one by one.
Let \(B=\lfloor\log_2 N\rfloor+1\). Building the nine moment tables up to bit length \(B\) takes \(O(B)\) time and \(O(B)\) memory, because each new level is obtained from the previous one by a constant amount of arithmetic. The prefix scan over \(N\) also costs \(O(B)\) time. For the actual input size here, the implementations use a fixed 64-bit bound, so both time and memory are effectively \(O(64)\).
Gesucht ist
$$S(N)=\sum_{n=1}^{N} f(n)^2 \pmod{10^9+7},$$
wobei
$$f(1)=1,\qquad f(2m)=2f(m),\qquad f(2m+1)=2m+1+2f(m)+\frac{f(m)}{m}\quad (m\ge 1).$$
Eine direkte Auswertung der Rekursion bis \(N=10^{16}\) ist unmöglich. Der entscheidende Schritt besteht darin, eine geschlossene Form für \(f(n)\) zu erkennen und diese dann mit einer binären Präfixzerlegung zu summieren.
Schreibe \(w(n)=\operatorname{popcount}(n)\) für die Anzahl der Einsen in der Binärdarstellung von \(n\). Die gesamte Lösung beruht auf der Identität \(f(n)=n\,w(n)\).
Für das Binärgewicht gilt
$$w(2m)=w(m),\qquad w(2m+1)=w(m)+1.$$
Angenommen, induktiv gelte \(f(m)=m\,w(m)\). Dann folgt für gerade Argumente
$$f(2m)=2f(m)=2m\,w(m)=2m\,w(2m).$$
Für ungerade Argumente wird der Divisionsterm harmlos, weil \(f(m)/m=w(m)\):
$$\begin{aligned} f(2m+1)&=2m+1+2f(m)+\frac{f(m)}{m}\\ &=2m+1+2m\,w(m)+w(m)\\ &=(2m+1)(w(m)+1)\\ &=(2m+1)w(2m+1). \end{aligned}$$
Da \(f(1)=1=1\cdot w(1)\), liefert die Induktion
$$f(n)=n\,w(n).$$
Damit wird die gesuchte Summe zu
$$S(N)=\sum_{n=1}^{N} n^2 w(n)^2.$$
Für jedes \(m\ge 0\) sei \(X_m=\{0,1,\dots,2^m-1\}\). Für \(a,b\in\{0,1,2\}\) definieren wir
$$T_m^{a,b}=\sum_{x\in X_m} x^a w(x)^b.$$
Diese neun Momente reichen aus, weil bei der Ausmultiplizierung eines Quadrats im Zahlenwert und eines Quadrats im Bitgewicht nur Terme \(x^a w(x)^b\) mit \(a,b\le 2\) auftreten.
Die Zielgröße über einen vollständigen \(m\)-Bit-Block ist der Spezialfall
$$T_m^{2,2}=\sum_{x=0}^{2^m-1} x^2 w(x)^2.$$
Der Anfang ist unmittelbar:
$$T_0^{0,0}=1,\qquad T_0^{a,b}=0\ \text{für}\ (a,b)\ne(0,0),$$
denn der einzige \(0\)-Bit-Suffix ist \(x=0\).
Jede \((m+1)\)-Bit-Zahl entsteht entweder als \(2x\) oder als \(2x+1\) mit \(x\in X_m\). Dabei gilt
$$w(2x)=w(x),\qquad w(2x+1)=w(x)+1.$$
Also erfüllt jedes neue Moment die Rekursion
$$T_{m+1}^{a,b}=\sum_{x\in X_m}(2x)^a w(x)^b+\sum_{x\in X_m}(2x+1)^a (w(x)+1)^b.$$
Weil \(a,b\le 2\), reduziert sich die rechte Seite stets wieder auf eine Linearkombination derselben neun Momente. Die kleinsten Übergänge sind
$$T_{m+1}^{0,0}=2T_m^{0,0},\qquad T_{m+1}^{0,1}=2T_m^{0,1}+T_m^{0,0},$$
während sich das Zielmoment zweiter Ordnung zu
$$\begin{aligned} T_{m+1}^{2,2}={}&8T_m^{2,2}+8T_m^{2,1}+4T_m^{1,2}+8T_m^{1,1}\\ &+4T_m^{2,0}+4T_m^{1,0}+T_m^{0,2}+2T_m^{0,1}+T_m^{0,0} \end{aligned}$$
entwickelt. Deshalb führen die Implementierungen genau neun Aggregattabellen und aktualisieren sie Ebenenweise.
In der Endsumme werden die höheren Bits zuerst festgelegt. Angenommen, das bereits gewählte Präfix hat den Wert \(p\) und enthält \(r\) Einsen. Wenn \(m\) niedrigere Bits noch frei sind, hat jede Zahl in diesem Block die Form
$$n=p+x,\qquad 0\le x\lt 2^m,$$
und ihr Bitgewicht ist
$$w(n)=r+w(x),$$
weil der freie Suffix nur Positionen unterhalb des festen Präfixes belegt.
Damit ist der gesamte Blockbeitrag
$$B(p,r,m)=\sum_{x=0}^{2^m-1}(p+x)^2(r+w(x))^2.$$
Nach dem Ausmultiplizieren beider Quadrate erhält man eine feste Linearkombination der vorberechneten Momente:
$$\begin{aligned} B(p,r,m)={}&p^2r^2T_m^{0,0}+2p^2rT_m^{0,1}+p^2T_m^{0,2}\\ &+2pr^2T_m^{1,0}+4prT_m^{1,1}+2pT_m^{1,2}\\ &+r^2T_m^{2,0}+2rT_m^{2,1}+T_m^{2,2}. \end{aligned}$$
Sobald die neun Momente bekannt sind, lässt sich also ein ganzes Intervall der Länge \(2^m\) in konstanter Zeit summieren.
Schreibe \(N\) binär und scanne von der höchstwertigen Eins bis zum niederwertigsten Bit. Immer wenn das aktuelle Bit von \(N\) gleich \(1\) ist, gibt es einen vollständigen Block kleinerer Zahlen: Man behält alle früheren Bits wie im aktuellen Präfix, setzt das aktuelle Bit auf \(0\) und lässt alle niedrigeren Bits frei variieren. Dieser gesamte Block wird mit \(B(p,r,m)\) addiert.
Danach wird das aktuelle Bit im Präfix eingeschaltet und die Präfix-Popcount um eins erhöht. Am Ende der Schleife ist das Präfix selbst gleich \(N\), also kommt noch der einzelne Endpunktterm
$$N^2 w(N)^2$$
hinzu. Jede Zahl von \(0\) bis \(N\) erscheint genau einmal in dieser Zerlegung, und der Beitrag von \(0\) ist automatisch null.
Da \(13=1101_2\) gilt, zerlegt der Scan den Bereich in vier Teile:
$$[0,7],\qquad [8,11],\qquad [12,12],\qquad \{13\}.$$
Diese vier Stücke entsprechen
$$B(0,0,3),\qquad B(8,1,2),\qquad B(12,2,0),\qquad 13^2\cdot 3^2.$$
Die Werte dazu sind
$$742,\qquad 1877,\qquad 576,\qquad 1521,$$
also
$$S(13)=742+1877+576+1521=4716.$$
Als weiterer kleiner Kontrollwert ergibt die geschlossene Form
$$S(10)=\sum_{n=1}^{10} n^2 w(n)^2=1530.$$
Die C++-, Python- und Java-Implementierungen folgen demselben Ablauf. Zuerst werden die neun Tabellen \(T_m^{a,b}\) für alle Bitlängen bis 64 vorab berechnet. Parallel dazu werden die Zweierpotenzen modulo \(10^9+7\) vorbereitet, damit ein festes Präfix sofort aktualisiert werden kann, sobald der Scan ein neues gesetztes Bit übernimmt.
Anschließend werden die Bits von \(N\) von oben nach unten durchlaufen. Bei jedem gesetzten Bit addiert die Implementierung den vollständigen Suffix-Block, der durch das aktuelle Präfix bestimmt ist, mithilfe der ausmultiplizierten Formel für \(B(p,r,m)\). Danach wird das Bit im Präfix gesetzt und die Popcount des Präfixes um eins erhöht. Nach der letzten Position kommt noch der einzelne Endpunktterm \(N^2 w(N)^2\) hinzu.
Alle Rechnungen erfolgen modulo \(10^9+7\), und zu keinem Zeitpunkt muss man die Zahlen \(1,2,\dots,N\) einzeln durchlaufen.
Sei \(B=\lfloor\log_2 N\rfloor+1\). Das Aufbauen der neun Momenttabellen bis zur Bitlänge \(B\) benötigt \(O(B)\) Zeit und \(O(B)\) Speicher, weil jede neue Ebene mit konstant viel Arithmetik aus der vorherigen entsteht. Das Präfix-Scannen von \(N\) kostet ebenfalls \(O(B)\) Zeit. Für die konkrete Eingabegrenze in dieser Aufgabe verwenden die Implementierungen eine feste 64-Bit-Schranke; damit sind Laufzeit und Speicher praktisch \(O(64)\).
Hesaplanması gereken ifade
$$S(N)=\sum_{n=1}^{N} f(n)^2 \pmod{10^9+7},$$
ve burada
$$f(1)=1,\qquad f(2m)=2f(m),\qquad f(2m+1)=2m+1+2f(m)+\frac{f(m)}{m}\quad (m\ge 1).$$
Bu özyinelemeyi \(N=10^{16}\)'ya kadar doğrudan yürütmek mümkün değildir. Esas fikir, önce \(f(n)\) için kapalı formu bulmak, sonra bu kapalı formu ikili önek bloklarına ayırarak toplamaktır.
\(w(n)=\operatorname{popcount}(n)\) olsun; yani \(n\)'in ikilik yazımındaki 1 bit sayısı. Tüm çözüm, \(f(n)=n\,w(n)\) özdeşliğine dayanır.
İkilik ağırlık için
$$w(2m)=w(m),\qquad w(2m+1)=w(m)+1$$
eşitlikleri vardır. İndüksiyon varsayımı olarak \(f(m)=m\,w(m)\) olsun. O zaman çift sayılar için
$$f(2m)=2f(m)=2m\,w(m)=2m\,w(2m)$$
olur. Tek sayılarda ise bölme terimi artık sorun çıkarmaz; çünkü \(f(m)/m=w(m)\):
$$\begin{aligned} f(2m+1)&=2m+1+2f(m)+\frac{f(m)}{m}\\ &=2m+1+2m\,w(m)+w(m)\\ &=(2m+1)(w(m)+1)\\ &=(2m+1)w(2m+1). \end{aligned}$$
\(f(1)=1=1\cdot w(1)\) olduğuna göre indüksiyonla
$$f(n)=n\,w(n)$$
elde edilir. Böylece aranan toplam
$$S(N)=\sum_{n=1}^{N} n^2 w(n)^2$$
şekline indirgenir.
Her \(m\ge 0\) için \(X_m=\{0,1,\dots,2^m-1\}\) olsun. \(a,b\in\{0,1,2\}\) için
$$T_m^{a,b}=\sum_{x\in X_m} x^a w(x)^b$$
tanımını yapalım. Sayının karesini ve popcount'un karesini açtığımızda ortaya çıkabilecek tüm terimler \(x^a w(x)^b\) biçimindedir ve burada \(a,b\le 2\) olur. Bu yüzden tam olarak dokuz moment yeterlidir.
Tam bir \(m\)-bit blok üzerindeki hedef toplam özel olarak
$$T_m^{2,2}=\sum_{x=0}^{2^m-1} x^2 w(x)^2$$
ifadesidir. Başlangıç değeri de doğrudan bellidir:
$$T_0^{0,0}=1,\qquad T_0^{a,b}=0\ \text{eğer}\ (a,b)\ne(0,0),$$
çünkü tek \(0\)-bit kuyruk \(x=0\)'dır.
Her \((m+1)\)-bit sayı, \(x\in X_m\) için ya \(2x\) ya da \(2x+1\) biçimindedir. Popcount için
$$w(2x)=w(x),\qquad w(2x+1)=w(x)+1$$
olduğundan her yeni moment
$$T_{m+1}^{a,b}=\sum_{x\in X_m}(2x)^a w(x)^b+\sum_{x\in X_m}(2x+1)^a (w(x)+1)^b$$
bağıntısını sağlar. \(a,b\le 2\) olduğu için sağ taraf yine aynı dokuz momentin lineer kombinasyonuna indirgenir. En basit geçişler
$$T_{m+1}^{0,0}=2T_m^{0,0},\qquad T_{m+1}^{0,1}=2T_m^{0,1}+T_m^{0,0}$$
şeklindedir. Asıl hedef olan ikinci mertebe momenti ise
$$\begin{aligned} T_{m+1}^{2,2}={}&8T_m^{2,2}+8T_m^{2,1}+4T_m^{1,2}+8T_m^{1,1}\\ &+4T_m^{2,0}+4T_m^{1,0}+T_m^{0,2}+2T_m^{0,1}+T_m^{0,0} \end{aligned}$$
şeklinde gelişir. Uygulamanın tam olarak dokuz tablo tutmasının nedeni budur.
Son toplam alınırken önce yüksek bitler sabitlenir. Seçilmiş önekin sayısal değeri \(p\), bu önekteki 1 sayısı \(r\) olsun. Eğer altta hâlâ \(m\) serbest bit varsa, bu bloktaki her sayı
$$n=p+x,\qquad 0\le x\lt 2^m$$
şeklindedir ve popcount'u da
$$w(n)=r+w(x)$$
olur; çünkü serbest kuyruk, sabit önekin altındaki bit konumlarını doldurur.
Bu yüzden tüm blok katkısı
$$B(p,r,m)=\sum_{x=0}^{2^m-1}(p+x)^2(r+w(x))^2$$
şeklindedir. İki kareyi açınca önceden hesaplanan momentlerin sabit bir lineer birleşimi elde edilir:
$$\begin{aligned} B(p,r,m)={}&p^2r^2T_m^{0,0}+2p^2rT_m^{0,1}+p^2T_m^{0,2}\\ &+2pr^2T_m^{1,0}+4prT_m^{1,1}+2pT_m^{1,2}\\ &+r^2T_m^{2,0}+2rT_m^{2,1}+T_m^{2,2}. \end{aligned}$$
Yani dokuz moment hazır olduğunda, uzunluğu \(2^m\) olan bir aralık sabit zamanda toplanabilir.
\(N\)'yi ikilik yaz ve en yüksek bitten en düşüğe doğru tara. O anda bakılan bit \(1\) ise, daha küçük sayılardan oluşan tam bir blok vardır: Önceki tüm bitler mevcut önekle aynı kalır, o bit \(0\) seçilir ve daha düşük bitler serbest bırakılır. İşte bu blok \(B(p,r,m)\) ile eklenir.
Daha sonra önekte ilgili bit \(1\) yapılır ve önekin popcount'u bir artırılır. Döngü bittiğinde önek doğrudan \(N\)'ye eşittir; bu yüzden tek kalan uç nokta terimi
$$N^2 w(N)^2$$
toplama eklenir. \(0\) ile \(N\) arasındaki her sayı bu ayrışımda tam bir kez yer alır; \(0\)'ın katkısı da zaten sıfırdır.
\(13=1101_2\) olduğundan tarama aralığı dört parçaya böler:
$$[0,7],\qquad [8,11],\qquad [12,12],\qquad \{13\}.$$
Bu parçalar sırasıyla
$$B(0,0,3),\qquad B(8,1,2),\qquad B(12,2,0),\qquad 13^2\cdot 3^2$$
ifadelerine karşılık gelir. Sayısal değerler
$$742,\qquad 1877,\qquad 576,\qquad 1521$$
olduğundan
$$S(13)=742+1877+576+1521=4716$$
elde edilir. Bir başka küçük kontrol noktası olarak kapalı form
$$S(10)=\sum_{n=1}^{10} n^2 w(n)^2=1530$$
sonucunu verir.
C++, Python ve Java uygulamaları aynı akışı izler. Önce 64 bite kadar tüm \(m\) seviyeleri için dokuz adet \(T_m^{a,b}\) tablosu hazırlanır. Buna paralel olarak \(10^9+7\) modunda ikinin kuvvetleri de önceden hesaplanır; böylece tarama sırasında sabit önek değeri yeni bir set bit geldiğinde anında güncellenebilir.
Daha sonra \(N\)'nin bitleri yukarıdan aşağı taranır. Her set bitte uygulama, mevcut önek tarafından belirlenen tam kuyruk bloğunu \(B(p,r,m)\) açılımıyla ekler. Ardından o bit öneğe dahil edilir ve önekin popcount'u bir artırılır. Bütün bitler işlendiğinde tek kalan uç nokta terimi \(N^2 w(N)^2\) da eklenir.
Tüm aritmetik \(10^9+7\) modunda yapılır ve algoritma hiçbir zaman \(1,2,\dots,N\) sayılarını tek tek dolaşmaz.
\(B=\lfloor\log_2 N\rfloor+1\) olsun. Dokuz moment tablosunu \(B\) bit uzunluğuna kadar oluşturmak \(O(B)\) zaman ve \(O(B)\) bellek gerektirir; çünkü her yeni seviye bir öncekinden sabit sayıda aritmetik işlemle elde edilir. \(N\) üzerindeki önek taraması da \(O(B)\) zamandadır. Bu problemde gerçek girdi için sabit 64-bit sınırı kullanıldığı için zaman ve bellek pratikte \(O(64)\) düzeyindedir.
Debemos calcular
$$S(N)=\sum_{n=1}^{N} f(n)^2 \pmod{10^9+7},$$
donde
$$f(1)=1,\qquad f(2m)=2f(m),\qquad f(2m+1)=2m+1+2f(m)+\frac{f(m)}{m}\quad (m\ge 1).$$
Evaluar esa recurrencia directamente hasta \(N=10^{16}\) es inviable. La clave es descubrir una forma cerrada para \(f(n)\) y luego sumar esa forma cerrada mediante una descomposición binaria por prefijos.
Escribamos \(w(n)=\operatorname{popcount}(n)\), es decir, el número de unos en la expansión binaria de \(n\). Toda la solución descansa sobre la identidad \(f(n)=n\,w(n)\).
El peso binario satisface
$$w(2m)=w(m),\qquad w(2m+1)=w(m)+1.$$
Supongamos inductivamente que \(f(m)=m\,w(m)\). Entonces, para argumentos pares,
$$f(2m)=2f(m)=2m\,w(m)=2m\,w(2m).$$
Para argumentos impares, el término con división deja de ser problemático porque \(f(m)/m=w(m)\):
$$\begin{aligned} f(2m+1)&=2m+1+2f(m)+\frac{f(m)}{m}\\ &=2m+1+2m\,w(m)+w(m)\\ &=(2m+1)(w(m)+1)\\ &=(2m+1)w(2m+1). \end{aligned}$$
Como \(f(1)=1=1\cdot w(1)\), la inducción da
$$f(n)=n\,w(n).$$
Por tanto, la suma buscada se convierte en
$$S(N)=\sum_{n=1}^{N} n^2 w(n)^2.$$
Para cada \(m\ge 0\), sea \(X_m=\{0,1,\dots,2^m-1\}\). Para \(a,b\in\{0,1,2\}\), definimos
$$T_m^{a,b}=\sum_{x\in X_m} x^a w(x)^b.$$
Estos nueve momentos bastan porque al expandir un cuadrado en el valor y otro cuadrado en el popcount solo pueden aparecer términos de la forma \(x^a w(x)^b\) con \(a,b\le 2\).
La cantidad objetivo sobre un bloque completo de \(m\) bits es el caso particular
$$T_m^{2,2}=\sum_{x=0}^{2^m-1} x^2 w(x)^2.$$
La base es inmediata:
$$T_0^{0,0}=1,\qquad T_0^{a,b}=0\ \text{si}\ (a,b)\ne(0,0),$$
porque el único sufijo de \(0\) bits es \(x=0\).
Cada número de \((m+1)\) bits procede de \(2x\) o de \(2x+1\), con \(x\in X_m\). Sus popcounts cumplen
$$w(2x)=w(x),\qquad w(2x+1)=w(x)+1.$$
Por ello, todo momento nuevo satisface
$$T_{m+1}^{a,b}=\sum_{x\in X_m}(2x)^a w(x)^b+\sum_{x\in X_m}(2x+1)^a (w(x)+1)^b.$$
Como \(a,b\le 2\), el lado derecho vuelve a expresarse como combinación lineal de los mismos nueve momentos. Las actualizaciones más pequeñas son
$$T_{m+1}^{0,0}=2T_m^{0,0},\qquad T_{m+1}^{0,1}=2T_m^{0,1}+T_m^{0,0},$$
mientras que el momento objetivo de segundo orden queda
$$\begin{aligned} T_{m+1}^{2,2}={}&8T_m^{2,2}+8T_m^{2,1}+4T_m^{1,2}+8T_m^{1,1}\\ &+4T_m^{2,0}+4T_m^{1,0}+T_m^{0,2}+2T_m^{0,1}+T_m^{0,0}. \end{aligned}$$
Por eso las implementaciones mantienen exactamente nueve tablas agregadas y las actualizan nivel por nivel.
En la suma final, primero se fijan los bits altos. Supongamos que el prefijo ya elegido aporta un valor \(p\) y contiene \(r\) unos. Si todavía quedan \(m\) bits bajos libres, cada número de ese bloque tiene la forma
$$n=p+x,\qquad 0\le x\lt 2^m,$$
y su popcount vale
$$w(n)=r+w(x),$$
porque el sufijo libre ocupa posiciones estrictamente inferiores a las del prefijo fijo.
Entonces la contribución total del bloque es
$$B(p,r,m)=\sum_{x=0}^{2^m-1}(p+x)^2(r+w(x))^2.$$
Al expandir ambos cuadrados se obtiene una combinación lineal fija de los momentos precomputados:
$$\begin{aligned} B(p,r,m)={}&p^2r^2T_m^{0,0}+2p^2rT_m^{0,1}+p^2T_m^{0,2}\\ &+2pr^2T_m^{1,0}+4prT_m^{1,1}+2pT_m^{1,2}\\ &+r^2T_m^{2,0}+2rT_m^{2,1}+T_m^{2,2}. \end{aligned}$$
Así, una vez conocidos los nueve momentos, un intervalo completo de longitud \(2^m\) se suma en tiempo constante.
Escribimos \(N\) en binario y recorremos sus bits desde el más significativo hasta el menos significativo. Siempre que el bit actual de \(N\) sea \(1\), existe un bloque completo de números menores: se mantienen los bits anteriores iguales al prefijo actual, se pone el bit actual a \(0\) y se deja variar libremente a todos los bits inferiores. Ese bloque se añade con \(B(p,r,m)\).
Después se enciende ese bit dentro del prefijo y se incrementa en uno el popcount del prefijo. Al terminar el recorrido, el propio prefijo ya es \(N\), así que solo falta añadir el término final
$$N^2 w(N)^2.$$
Cada entero entre \(0\) y \(N\) aparece exactamente una vez en esta descomposición, y el término correspondiente a \(0\) vale automáticamente cero.
Como \(13=1101_2\), el recorrido divide el rango en cuatro partes:
$$[0,7],\qquad [8,11],\qquad [12,12],\qquad \{13\}.$$
Estas cuatro piezas corresponden a
$$B(0,0,3),\qquad B(8,1,2),\qquad B(12,2,0),\qquad 13^2\cdot 3^2.$$
Sus valores son
$$742,\qquad 1877,\qquad 576,\qquad 1521,$$
de modo que
$$S(13)=742+1877+576+1521=4716.$$
Como verificación pequeña adicional, la misma forma cerrada da
$$S(10)=\sum_{n=1}^{10} n^2 w(n)^2=1530.$$
Las implementaciones en C++, Python y Java siguen el mismo esquema. Primero precomputan las nueve tablas \(T_m^{a,b}\) para todas las longitudes de bit hasta 64. Al mismo tiempo, preparan las potencias de dos módulo \(10^9+7\), para que el valor numérico del prefijo pueda actualizarse en cuanto el recorrido acepta un nuevo bit encendido.
Después recorren los bits de \(N\) de arriba hacia abajo. En cada bit igual a uno, la implementación suma el bloque completo de sufijo determinado por el prefijo actual usando la fórmula expandida de \(B(p,r,m)\). Luego incorpora ese bit al prefijo y aumenta en uno el popcount del prefijo. Cuando todos los bits ya fueron procesados, añade el término final \(N^2 w(N)^2\).
Toda la aritmética se realiza módulo \(10^9+7\), y el algoritmo nunca recorre uno por uno los enteros \(1,2,\dots,N\).
Sea \(B=\lfloor\log_2 N\rfloor+1\). Construir las nueve tablas de momentos hasta longitud \(B\) requiere \(O(B)\) tiempo y \(O(B)\) memoria, porque cada nivel nuevo se obtiene del anterior con una cantidad constante de operaciones aritméticas. El recorrido por prefijos sobre \(N\) también cuesta \(O(B)\) tiempo. Para el tamaño real de entrada de este problema, las implementaciones usan un límite fijo de 64 bits, así que tanto el tiempo como la memoria son efectivamente \(O(64)\).
我们要求的是
$$S(N)=\sum_{n=1}^{N} f(n)^2 \pmod{10^9+7},$$
其中
$$f(1)=1,\qquad f(2m)=2f(m),\qquad f(2m+1)=2m+1+2f(m)+\frac{f(m)}{m}\quad (m\ge 1).$$
如果按递推式一直算到 \(N=10^{16}\)，复杂度根本无法接受。真正有效的做法，是先把 \(f(n)\) 化成显式公式，再用二进制前缀分块的方法一次性求整段区间的贡献。
记 \(w(n)=\operatorname{popcount}(n)\)，也就是 \(n\) 的二进制表示中 1 的个数。整个解法的核心是证明
$$f(n)=n\,w(n).$$
二进制位权满足
$$w(2m)=w(m),\qquad w(2m+1)=w(m)+1.$$
假设归纳地有 \(f(m)=m\,w(m)\)。那么对偶数项，直接得到
$$f(2m)=2f(m)=2m\,w(m)=2m\,w(2m).$$
对奇数项，除法项之所以安全，是因为此时 \(f(m)/m=w(m)\)：
$$\begin{aligned} f(2m+1)&=2m+1+2f(m)+\frac{f(m)}{m}\\ &=2m+1+2m\,w(m)+w(m)\\ &=(2m+1)(w(m)+1)\\ &=(2m+1)w(2m+1). \end{aligned}$$
再结合初值 \(f(1)=1=1\cdot w(1)\)，归纳可得
$$f(n)=n\,w(n).$$
于是原题就等价于求
$$S(N)=\sum_{n=1}^{N} n^2 w(n)^2.$$
对每个 \(m\ge 0\)，定义集合 \(X_m=\{0,1,\dots,2^m-1\}\)。对 \(a,b\in\{0,1,2\}\)，定义
$$T_m^{a,b}=\sum_{x\in X_m} x^a w(x)^b.$$
之所以只需要这九个量，是因为我们最终只会展开“数值平方”和“popcount 平方”，所以出现的项一定都形如 \(x^a w(x)^b\)，并且 \(a,b\le 2\)。
在一个完整的 \(m\) 位区间上，我们真正想要的是
$$T_m^{2,2}=\sum_{x=0}^{2^m-1} x^2 w(x)^2.$$
初始层非常简单：
$$T_0^{0,0}=1,\qquad T_0^{a,b}=0\ \text{当}\ (a,b)\ne(0,0),$$
因为唯一的 \(0\) 位后缀就是 \(x=0\)。
每个 \((m+1)\) 位数，都可以由某个 \(x\in X_m\) 生成成 \(2x\) 或 \(2x+1\)。对应的 popcount 满足
$$w(2x)=w(x),\qquad w(2x+1)=w(x)+1.$$
因此任意一个新矩都满足
$$T_{m+1}^{a,b}=\sum_{x\in X_m}(2x)^a w(x)^b+\sum_{x\in X_m}(2x+1)^a (w(x)+1)^b.$$
由于 \(a,b\le 2\)，右边展开后仍然只会回到同样的九个矩上。最简单的两个更新是
$$T_{m+1}^{0,0}=2T_m^{0,0},\qquad T_{m+1}^{0,1}=2T_m^{0,1}+T_m^{0,0},$$
而目标二阶矩的展开式是
$$\begin{aligned} T_{m+1}^{2,2}={}&8T_m^{2,2}+8T_m^{2,1}+4T_m^{1,2}+8T_m^{1,1}\\ &+4T_m^{2,0}+4T_m^{1,0}+T_m^{0,2}+2T_m^{0,1}+T_m^{0,0}. \end{aligned}$$
这也解释了为什么实现中恰好维护九张聚合表，并逐层从短后缀推到长后缀。
在最终求和时，先固定高位。设当前已经固定的前缀数值为 \(p\)，其中 1 的个数为 \(r\)。如果还有 \(m\) 个低位没有决定，那么这一整块中的每个数都可以写成
$$n=p+x,\qquad 0\le x\lt 2^m,$$
并且它的 popcount 是
$$w(n)=r+w(x),$$
因为自由后缀只占据固定前缀之下的那些位，不会和高位重叠。
于是整个区间块的贡献为
$$B(p,r,m)=\sum_{x=0}^{2^m-1}(p+x)^2(r+w(x))^2.$$
把两个平方完全展开后，就得到九个预处理矩的固定线性组合：
$$\begin{aligned} B(p,r,m)={}&p^2r^2T_m^{0,0}+2p^2rT_m^{0,1}+p^2T_m^{0,2}\\ &+2pr^2T_m^{1,0}+4prT_m^{1,1}+2pT_m^{1,2}\\ &+r^2T_m^{2,0}+2rT_m^{2,1}+T_m^{2,2}. \end{aligned}$$
也就是说，只要九个矩准备好了，长度为 \(2^m\) 的整段区间就可以 \(O(1)\) 求出贡献。
把 \(N\) 写成二进制，并从最高位一路扫到最低位。每当当前位是 \(1\) 时，就会出现一个“完整的更小区间块”：更高位与当前前缀完全相同，把当前位改成 \(0\)，而所有更低位都允许自由变化。这个整块的贡献就用 \(B(p,r,m)\) 一次加入。
加入之后，再把这一个 1 位真正并入前缀，同时把前缀的 popcount 加一。扫描结束时，前缀本身就已经等于 \(N\)，因此还要补上最后的单点项
$$N^2 w(N)^2.$$
这样一来，区间 \(0\) 到 \(N\) 中的每个整数都会被恰好计入一次，而 \(0\) 的贡献本身就是零。
因为 \(13=1101_2\)，按位扫描会把区间拆成四部分：
$$[0,7],\qquad [8,11],\qquad [12,12],\qquad \{13\}.$$
它们分别对应
$$B(0,0,3),\qquad B(8,1,2),\qquad B(12,2,0),\qquad 13^2\cdot 3^2.$$
数值上这四项是
$$742,\qquad 1877,\qquad 576,\qquad 1521,$$
所以
$$S(13)=742+1877+576+1521=4716.$$
再给出一个更小的校验点，闭式同样给出
$$S(10)=\sum_{n=1}^{10} n^2 w(n)^2=1530.$$
C++、Python 和 Java 三个实现遵循完全相同的流程。第一步是为所有不超过 64 的位长预处理九张 \(T_m^{a,b}\) 表。同时还会预处理二的幂在 \(10^9+7\) 模下的值，这样在扫描过程中只要接受一个新的高位 1，就能立刻更新前缀的数值。
接下来，程序从高位到低位扫描 \(N\) 的二进制位。每遇到一个 1 位，就把“当前前缀决定的完整后缀块”通过 \(B(p,r,m)\) 的展开式加入答案。然后把该位并入前缀，并让前缀的 popcount 加一。全部位处理完以后，再把单独的端点项 \(N^2 w(N)^2\) 加上。
整个过程中所有运算都在模 \(10^9+7\) 下进行，而且算法从头到尾都不需要逐个枚举 \(1,2,\dots,N\)。
设 \(B=\lfloor\log_2 N\rfloor+1\)。把九张矩表预处理到位长 \(B\) 需要 \(O(B)\) 时间和 \(O(B)\) 空间，因为每一层都只是由前一层经过常数次算术运算得到。之后对 \(N\) 做一次前缀按位扫描，同样是 \(O(B)\) 时间。对这道题的实际输入规模来说，实现直接采用固定的 64 位上界，因此时间和空间都可以视为 \(O(64)\)。
Требуется вычислить
$$S(N)=\sum_{n=1}^{N} f(n)^2 \pmod{10^9+7},$$
где
$$f(1)=1,\qquad f(2m)=2f(m),\qquad f(2m+1)=2m+1+2f(m)+\frac{f(m)}{m}\quad (m\ge 1).$$
Прямое применение этой рекурсии до \(N=10^{16}\) нереально. Эффективное решение сначала выводит явную формулу для \(f(n)\), а затем суммирует ее с помощью двоичного префиксного разбиения.
Обозначим через \(w(n)=\operatorname{popcount}(n)\) число единиц в двоичной записи \(n\). Вся идея держится на тождестве \(f(n)=n\,w(n)\).
Для двоичного веса верно
$$w(2m)=w(m),\qquad w(2m+1)=w(m)+1.$$
Предположим по индукции, что \(f(m)=m\,w(m)\). Тогда для четных аргументов
$$f(2m)=2f(m)=2m\,w(m)=2m\,w(2m).$$
Для нечетных аргументов деление тоже упрощается, потому что \(f(m)/m=w(m)\):
$$\begin{aligned} f(2m+1)&=2m+1+2f(m)+\frac{f(m)}{m}\\ &=2m+1+2m\,w(m)+w(m)\\ &=(2m+1)(w(m)+1)\\ &=(2m+1)w(2m+1). \end{aligned}$$
Так как \(f(1)=1=1\cdot w(1)\), индукция дает точную формулу
$$f(n)=n\,w(n).$$
Следовательно, исходная сумма равна
$$S(N)=\sum_{n=1}^{N} n^2 w(n)^2.$$
Для каждого \(m\ge 0\) положим \(X_m=\{0,1,\dots,2^m-1\}\). Для \(a,b\in\{0,1,2\}\) определим
$$T_m^{a,b}=\sum_{x\in X_m} x^a w(x)^b.$$
Именно этих девяти моментов достаточно, потому что при раскрытии квадрата по значению числа и квадрата по popcount возникают только выражения вида \(x^a w(x)^b\) с \(a,b\le 2\).
Целевой вклад на полном \(m\)-битном блоке есть частный случай
$$T_m^{2,2}=\sum_{x=0}^{2^m-1} x^2 w(x)^2.$$
База очевидна:
$$T_0^{0,0}=1,\qquad T_0^{a,b}=0\ \text{при}\ (a,b)\ne(0,0),$$
поскольку единственный \(0\)-битный суффикс это \(x=0\).
Любое \((m+1)\)-битное число получается либо как \(2x\), либо как \(2x+1\), где \(x\in X_m\). Для popcount имеем
$$w(2x)=w(x),\qquad w(2x+1)=w(x)+1.$$
Отсюда каждый новый момент удовлетворяет формуле
$$T_{m+1}^{a,b}=\sum_{x\in X_m}(2x)^a w(x)^b+\sum_{x\in X_m}(2x+1)^a (w(x)+1)^b.$$
Так как \(a,b\le 2\), правая часть снова выражается через те же девять моментов. Самые простые переходы таковы:
$$T_{m+1}^{0,0}=2T_m^{0,0},\qquad T_{m+1}^{0,1}=2T_m^{0,1}+T_m^{0,0},$$
а для целевого момента второго порядка получается
$$\begin{aligned} T_{m+1}^{2,2}={}&8T_m^{2,2}+8T_m^{2,1}+4T_m^{1,2}+8T_m^{1,1}\\ &+4T_m^{2,0}+4T_m^{1,0}+T_m^{0,2}+2T_m^{0,1}+T_m^{0,0}. \end{aligned}$$
Именно поэтому реализации хранят ровно девять агрегированных таблиц и обновляют их по длине суффикса.
При окончательном суммировании сначала фиксируются старшие биты. Пусть уже выбранный префикс дает численное значение \(p\) и содержит \(r\) единиц. Если еще свободны \(m\) младших битов, то каждое число в таком блоке имеет вид
$$n=p+x,\qquad 0\le x\lt 2^m,$$
а его popcount равен
$$w(n)=r+w(x),$$
потому что свободный суффикс занимает позиции строго ниже фиксированного префикса.
Значит, вклад всего блока равен
$$B(p,r,m)=\sum_{x=0}^{2^m-1}(p+x)^2(r+w(x))^2.$$
После раскрытия квадратов получаем фиксированную линейную комбинацию предвычисленных моментов:
$$\begin{aligned} B(p,r,m)={}&p^2r^2T_m^{0,0}+2p^2rT_m^{0,1}+p^2T_m^{0,2}\\ &+2pr^2T_m^{1,0}+4prT_m^{1,1}+2pT_m^{1,2}\\ &+r^2T_m^{2,0}+2rT_m^{2,1}+T_m^{2,2}. \end{aligned}$$
Следовательно, если девять моментов уже известны, целый интервал длины \(2^m\) суммируется за постоянное время.
Запишем \(N\) в двоичном виде и будем идти от старшего бита к младшему. Каждый раз, когда текущий бит числа \(N\) равен \(1\), существует полный блок меньших чисел: все более старшие биты совпадают с текущим префиксом, текущий бит ставится в \(0\), а все младшие биты могут быть любыми. Такой блок добавляется через \(B(p,r,m)\).
После этого текущий бит включается в префикс, а popcount префикса увеличивается на единицу. Когда цикл заканчивается, сам префикс уже равен \(N\), поэтому остается прибавить одиночный конечный член
$$N^2 w(N)^2.$$
Каждое число от \(0\) до \(N\) попадает в это разбиение ровно один раз, а вклад числа \(0\) автоматически равен нулю.
Поскольку \(13=1101_2\), сканирование разбивает диапазон на четыре части:
$$[0,7],\qquad [8,11],\qquad [12,12],\qquad \{13\}.$$
Этим частям соответствуют выражения
$$B(0,0,3),\qquad B(8,1,2),\qquad B(12,2,0),\qquad 13^2\cdot 3^2.$$
Их численные значения равны
$$742,\qquad 1877,\qquad 576,\qquad 1521,$$
поэтому
$$S(13)=742+1877+576+1521=4716.$$
Еще одна маленькая контрольная точка:
$$S(10)=\sum_{n=1}^{10} n^2 w(n)^2=1530.$$
Реализации на C++, Python и Java действуют одинаково. Сначала они предвычисляют девять таблиц \(T_m^{a,b}\) для всех длин битов до 64. Одновременно вычисляются степени двойки по модулю \(10^9+7\), чтобы при принятии очередного установленного бита можно было сразу обновить численное значение префикса.
Затем биты числа \(N\) просматриваются сверху вниз. На каждом установленном бите реализация добавляет полный суффиксный блок, который задается текущим префиксом, используя развернутую формулу для \(B(p,r,m)\). После этого бит включается в префикс, а popcount префикса увеличивается на единицу. После обработки всех битов прибавляется одиночный конечный член \(N^2 w(N)^2\).
Вся арифметика ведется по модулю \(10^9+7\), и алгоритм никогда не перебирает числа \(1,2,\dots,N\) по одному.
Пусть \(B=\lfloor\log_2 N\rfloor+1\). Построение девяти таблиц моментов до длины \(B\) требует \(O(B)\) времени и \(O(B)\) памяти, поскольку каждый новый уровень получается из предыдущего постоянным числом арифметических операций. Само префиксное сканирование числа \(N\) также занимает \(O(B)\) времени. Для реального ограничения задачи используется фиксированный 64-битный предел, поэтому время и память фактически равны \(O(64)\).
نريد حساب
$$S(N)=\sum_{n=1}^{N} f(n)^2 \pmod{10^9+7},$$
حيث
$$f(1)=1,\qquad f(2m)=2f(m),\qquad f(2m+1)=2m+1+2f(m)+\frac{f(m)}{m}\quad (m\ge 1).$$
ولا يمكن تنفيذ هذه العلاقة التكرارية مباشرة حتى \(N=10^{16}\). لذلك يعتمد الحل على استخراج صيغة مغلقة لـ \(f(n)\)، ثم جمعها باستعمال تفكيك ثنائي قائم على البادئات بدل المرور على كل عدد على حدة.
لنكتب \(w(n)=\operatorname{popcount}(n)\)، أي عدد الواحدات في التمثيل الثنائي للعدد \(n\). جوهر الحل كله هو الهوية
$$f(n)=n\,w(n).$$
يعطي الوزن الثنائي العلاقتين
$$w(2m)=w(m),\qquad w(2m+1)=w(m)+1.$$
ولنفترض استقرائيًا أن \(f(m)=m\,w(m)\). عندئذٍ في الحالة الزوجية نحصل على
$$f(2m)=2f(m)=2m\,w(m)=2m\,w(2m).$$
أما في الحالة الفردية فإن حد القسمة يصبح بسيطًا لأن \(f(m)/m=w(m)\):
$$\begin{aligned} f(2m+1)&=2m+1+2f(m)+\frac{f(m)}{m}\\ &=2m+1+2m\,w(m)+w(m)\\ &=(2m+1)(w(m)+1)\\ &=(2m+1)w(2m+1). \end{aligned}$$
ومع القيمة الابتدائية \(f(1)=1=1\cdot w(1)\)، ينتج بالاستقراء أن
$$f(n)=n\,w(n).$$
ومن ثم تصبح المسألة هي حساب
$$S(N)=\sum_{n=1}^{N} n^2 w(n)^2.$$
لكل \(m\ge 0\)، عرّف المجموعة \(X_m=\{0,1,\dots,2^m-1\}\). ثم لكل \(a,b\in\{0,1,2\}\) نضع
$$T_m^{a,b}=\sum_{x\in X_m} x^a w(x)^b.$$
هذه العزوم التسعة تكفي تمامًا، لأننا لا نحتاج إلا إلى توسيع مربع في قيمة العدد ومربع في عدد الواحدات، وبالتالي لا تظهر إلا حدود من الشكل \(x^a w(x)^b\) مع \(a,b\le 2\).
والكمية الهدف على كتلة كاملة من أعداد \(m\) بت هي الحالة الخاصة
$$T_m^{2,2}=\sum_{x=0}^{2^m-1} x^2 w(x)^2.$$
أما الحالة الابتدائية فهي مباشرة:
$$T_0^{0,0}=1,\qquad T_0^{a,b}=0\ \text{إذا كان}\ (a,b)\ne(0,0),$$
لأن اللاحقة الوحيدة ذات \(0\) بت هي \(x=0\).
كل عدد بطول \((m+1)\) بت يأتي إما من \(2x\) أو من \(2x+1\) حيث \(x\in X_m\). ويحقق popcount عندئذٍ
$$w(2x)=w(x),\qquad w(2x+1)=w(x)+1.$$
لذلك فإن كل عزم جديد يحقق العلاقة
$$T_{m+1}^{a,b}=\sum_{x\in X_m}(2x)^a w(x)^b+\sum_{x\in X_m}(2x+1)^a (w(x)+1)^b.$$
وبما أن \(a,b\le 2\)، فإن الطرف الأيمن يعود دائمًا إلى تركيب خطي من العزوم التسعة نفسها. وأبسط انتقالين هما
$$T_{m+1}^{0,0}=2T_m^{0,0},\qquad T_{m+1}^{0,1}=2T_m^{0,1}+T_m^{0,0},$$
بينما يتمدد عزم الرتبة الثانية الهدف إلى
$$\begin{aligned} T_{m+1}^{2,2}={}&8T_m^{2,2}+8T_m^{2,1}+4T_m^{1,2}+8T_m^{1,1}\\ &+4T_m^{2,0}+4T_m^{1,0}+T_m^{0,2}+2T_m^{0,1}+T_m^{0,0}. \end{aligned}$$
ولهذا السبب بالضبط تحتفظ التطبيقات بتسع جداول تجميعية لا أكثر ولا أقل.
في مرحلة الجمع النهائية تُثبَّت البتات العليا أولًا. لنفترض أن البادئة الحالية تعطي قيمة عددية \(p\) وتحتوي على \(r\) واحدات. وإذا بقيت \(m\) بتات دنيا حرة، فإن كل عدد في تلك الكتلة يكتب على الصورة
$$n=p+x,\qquad 0\le x\lt 2^m,$$
ويكون عدد الواحدات فيه
$$w(n)=r+w(x),$$
لأن اللاحقة الحرة تشغل مواضع أدنى من مواضع البادئة المثبتة.
إذن مساهمة الكتلة كلها هي
$$B(p,r,m)=\sum_{x=0}^{2^m-1}(p+x)^2(r+w(x))^2.$$
وعند توسيع المربعين نحصل على تركيب خطي ثابت من العزوم المحسوبة مسبقًا:
$$\begin{aligned} B(p,r,m)={}&p^2r^2T_m^{0,0}+2p^2rT_m^{0,1}+p^2T_m^{0,2}\\ &+2pr^2T_m^{1,0}+4prT_m^{1,1}+2pT_m^{1,2}\\ &+r^2T_m^{2,0}+2rT_m^{2,1}+T_m^{2,2}. \end{aligned}$$
وبالتالي، ما إن تكون العزوم التسعة جاهزة، يصبح جمع فترة كاملة طولها \(2^m\) عملية ثابتة الزمن.
اكتب \(N\) بالثنائي وامسح البتات من الأعلى إلى الأسفل. وكلما كان البت الحالي في \(N\) مساويًا لـ \(1\)، وُجدت كتلة كاملة من الأعداد الأصغر: نبقي جميع البتات الأعلى مساوية للبادئة الحالية، ونضع البت الحالي صفرًا، ونترك كل البتات الأدنى حرة. هذه الكتلة تضاف دفعة واحدة بواسطة \(B(p,r,m)\).
بعد ذلك نفعّل هذا البت داخل البادئة ونزيد عدد الواحدات في البادئة بمقدار واحد. وعند نهاية المسح تصبح البادئة نفسها مساوية لـ \(N\)، فلا يبقى إلا إضافة الحد المفرد الأخير
$$N^2 w(N)^2.$$
وبذلك يظهر كل عدد من \(0\) إلى \(N\) مرة واحدة بالضبط في هذا التفكيك، أما مساهمة الصفر فهي صفر تلقائيًا.
بما أن \(13=1101_2\)، فإن المسح يقسم المجال إلى أربع قطع:
$$[0,7],\qquad [8,11],\qquad [12,12],\qquad \{13\}.$$
وهذه القطع تقابل
$$B(0,0,3),\qquad B(8,1,2),\qquad B(12,2,0),\qquad 13^2\cdot 3^2.$$
وقيمها العددية هي
$$742,\qquad 1877,\qquad 576,\qquad 1521,$$
ومن ثم
$$S(13)=742+1877+576+1521=4716.$$
وكفحص صغير إضافي، تعطي الصيغة المغلقة كذلك
$$S(10)=\sum_{n=1}^{10} n^2 w(n)^2=1530.$$
تتبع تطبيقات C++ وPython وJava المسار نفسه. أولًا تُحسب مسبقًا الجداول التسعة \(T_m^{a,b}\) لكل أطوال البت حتى 64. وفي الوقت نفسه تُجهَّز قوى العدد 2 بترديد \(10^9+7\)، حتى يمكن تحديث القيمة العددية للبادئة فور قبول أي بت جديد يساوي \(1\).
بعد ذلك تُمسح بتات \(N\) من الأعلى إلى الأسفل. وعند كل بت يساوي واحدًا، تضيف الخوارزمية الكتلة الكاملة للواحق التي يحددها الوضع الحالي للبادئة، وذلك باستخدام الصيغة الموسعة لـ \(B(p,r,m)\). ثم تُمدَّد البادئة بإضافة هذا البت، ويزداد عدد الواحدات في البادئة بمقدار واحد. وبعد الانتهاء من جميع البتات، يضاف الحد النهائي المفرد \(N^2 w(N)^2\).
كل الحسابات تجرى بترديد \(10^9+7\)، ولا تحتاج الخوارزمية في أي لحظة إلى المرور على الأعداد \(1,2,\dots,N\) واحدًا واحدًا.
لنضع \(B=\lfloor\log_2 N\rfloor+1\). إن بناء الجداول التسعة حتى طول \(B\) يحتاج إلى زمن \(O(B)\) وذاكرة \(O(B)\)، لأن كل مستوى جديد يُشتق من المستوى السابق بعدد ثابت من العمليات الحسابية. كما أن مسح بادئة العدد \(N\) يستهلك أيضًا \(O(B)\) من الزمن. وبالنسبة لحجم الإدخال الفعلي في هذه المسألة، تستخدم التطبيقات حدًا ثابتًا مقداره 64 بت، لذا يصبح الزمن والذاكرة فعليًا \(O(64)\).