The blancmange function, also called the Takagi function, is the continuous nowhere-differentiable curve
$$B(x)=\sum_{n=0}^{\infty}\frac{\phi(2^n x)}{2^n},\qquad \phi(u)=\operatorname{dist}(u,\mathbb{Z}).$$
Problem 226 asks for the area common to the region under this curve and the circle centered at \(\left(\tfrac14,\tfrac12\right)\) with radius \(\tfrac14\). The implementations do not search for a symbolic antiderivative. Instead, they turn the geometry into a one-dimensional integral and make every sample value of the blancmange function exact by exploiting its behavior on dyadic rationals.
The key observation is that Simpson's rule samples points of the form \(x=k/2^m\), and on exactly those points the blancmange series collapses to a finite binary-digit formula. That is why the numerical integration is efficient and reliable even though the curve itself is fractal.
If \(x=k/2^m\), then \(2^n x\) is an integer for every \(n \ge m\), so \(\phi(2^n x)=0\) from that point onward. Therefore
$$B\left(\frac{k}{2^m}\right)=\sum_{n=0}^{m-1}\frac{\phi\!\left(2^n\frac{k}{2^m}\right)}{2^n},$$
and the infinite series becomes a finite sum of \(m\) terms. The implementations go one step further and avoid even that \(m\)-term evaluation by converting the value into a formula involving binary digit sums.
Let \(s_2(r)\) be the number of 1-bits in the binary expansion of \(r\), and define the prefix sum
$$P(k)=\sum_{r=0}^{k-1}s_2(r).$$
For dyadic points one has the exact identity
$$B\left(\frac{k}{2^m}\right)=\frac{mk-2P(k)}{2^m}.$$
This is the central formula used by all three implementations. It immediately yields the recurrence
$$B\left(\frac{k+1}{2^m}\right)-B\left(\frac{k}{2^m}\right)=\frac{m-2s_2(k)}{2^m},$$
because \(P(k+1)=P(k)+s_2(k)\). So once one dyadic value is known, the next one is obtained by a constant-time update controlled only by the bit count of the current index.
The circle is
$$\left(x-\frac14\right)^2+\left(y-\frac12\right)^2=\left(\frac14\right)^2.$$
It occupies only the horizontal range \(0 \le x \le \tfrac12\). For such an \(x\), the corresponding vertical slice of the circle is
$$y_{\pm}(x)=\frac12 \pm \sqrt{\frac1{16}-\left(x-\frac14\right)^2}.$$
The blancmange region is the set
$$R_B=\{(x,y): 0 \le y \le B(x)\}.$$
Hence the common vertical height at position \(x\) is
$$h(x)=\max\!\left(0,\ \min\!\bigl(B(x),y_+(x)\bigr)-\max\!\bigl(0,y_-(x)\bigr)\right).$$
The required area is therefore
$$A=\int_0^{1/2} h(x)\,dx.$$
Outside \([0,\tfrac12]\) the circle contributes nothing, so the integral really is one-dimensional and compactly supported.
Take \(m=3\) and \(k=3\), so \(x=3/8\). The prefix bit-count sum is
$$P(3)=s_2(0)+s_2(1)+s_2(2)=0+1+1=2,$$
and therefore
$$B\left(\frac38\right)=\frac{3 \cdot 3-2 \cdot 2}{8}=\frac58.$$
This agrees with the direct Takagi sum:
$$B\left(\frac38\right)=\phi\left(\frac38\right)+\frac12\phi\left(\frac34\right)+\frac14\phi\left(\frac32\right)=\frac38+\frac18+\frac18=\frac58.$$
For the circle slice, \(x-\tfrac14=\tfrac18\), so
$$y_-\left(\frac38\right)=\frac12-\frac{\sqrt3}{8},\qquad y_+\left(\frac38\right)=\frac12+\frac{\sqrt3}{8}.$$
Since \(\tfrac58 < y_+(\tfrac38)\), the overlap height is
$$h\left(\frac38\right)=\frac58-\left(\frac12-\frac{\sqrt3}{8}\right)=\frac{1+\sqrt3}{8}.$$
This is exactly what the integrand measures at each sample point: the portion of the circle's vertical segment that lies below the blancmange curve and above the \(x\)-axis.
At refinement level \(m\), the interval \([0,\tfrac12]\) is divided into \(N=2^{m-1}\) equal subintervals of width \(h=2^{-m}\), so every node is dyadic:
$$x_i=\frac{i}{2^m},\qquad 0 \le i \le 2^{m-1}.$$
That means each sampled value \(B(x_i)\) is computed exactly by the formula above. The only approximation comes from the quadrature itself. Simpson's rule gives
$$A_m=\frac{h}{3}\left(h(x_0)+h(x_N)+4\sum_{\substack{1 \le i \le N-1 \\ i\text{ odd}}} h(x_i)+2\sum_{\substack{2 \le i \le N-2 \\ i\text{ even}}} h(x_i)\right).$$
The implementations evaluate this on successively finer dyadic grids and stop when two consecutive levels differ by less than the requested tolerance.
To evaluate \(P(k)\) quickly, the implementations count 1-bits one position at a time. For bit position \(b\), the pattern repeats every \(2^{b+1}\) integers, and the number of ones among \(0,1,\dots,k-1\) is
$$\left\lfloor\frac{k}{2^{b+1}}\right\rfloor 2^b+\max\!\left(0,\ k \bmod 2^{b+1}-2^b\right).$$
Summing that over bit positions yields \(P(k)\) in \(O(\log k)\) arithmetic time rather than by iterating through all previous integers.
During one Simpson sweep at fixed level \(m\), the implementations maintain the numerator \(mk-2P(k)\) rather than recomputing it from scratch at every node. Moving from index \(k\) to \(k+1\) changes this numerator by \(m-2s_2(k)\), so the next blancmange value is obtained with one bit count and a few arithmetic operations. This invariant is the reason a very fine dyadic grid is still practical.
For each sample point, the code computes the circle's lower and upper \(y\)-values, clips that interval against \([0,B(x)]\), and adds the resulting overlap height to either the odd or the even Simpson accumulator. The interior sample range is split into independent chunks so partial sums can be formed separately and then combined. The C++, Python, and Java implementations all use the same mathematics; the compiled versions divide the work across threads, and the Python version can distribute chunks across worker processes when that is worthwhile. After one level is finished, the next level is computed and compared against it; refinement stops as soon as the tolerance test passes.
At level \(m\), Simpson's rule on \([0,\tfrac12]\) uses \(2^{m-1}+1\) sample points, so one refinement level costs \(O(2^m)\) time. Because the dyadic blancmange values are streamed by the recurrence above, each additional node is only a constant amount of arithmetic plus a bit count and a square root.
The extra memory is \(O(T)\) for \(T\) worker partial sums, or \(O(1)\) in a serial evaluation. Since the code checks a short sequence of increasing dyadic levels, the total running time up to convergence is still on the order of the finest level that gets evaluated.
Die Blancmange- oder Takagi-Funktion ist die stetige, aber nirgends differenzierbare Kurve
$$B(x)=\sum_{n=0}^{\infty}\frac{\phi(2^n x)}{2^n},\qquad \phi(u)=\operatorname{dist}(u,\mathbb{Z}).$$
In Problem 226 wird die gemeinsame Fläche der Region unter dieser Kurve mit dem Kreis gesucht, dessen Mittelpunkt \(\left(\tfrac14,\tfrac12\right)\) und dessen Radius \(\tfrac14\) ist. Die Implementierungen suchen keine geschlossene Stammfunktion. Stattdessen wird die Geometrie auf ein eindimensionales Integral reduziert, und die Blancmange-Werte an allen Stützstellen werden dank der Dyadenstruktur exakt berechnet.
Der entscheidende Punkt ist, dass Simpson-Quadratur nur Punkte der Form \(x=k/2^m\) benötigt. Genau an diesen dyadischen Stellen lässt sich die Blancmange-Funktion durch eine Formel mit Binärziffernsummen auswerten, ohne die ganze unendliche Reihe direkt zu summieren.
Für \(x=k/2^m\) ist \(2^n x\) für alle \(n \ge m\) ganzzahlig, also verschwindet dann \(\phi(2^n x)\). Daher gilt
$$B\left(\frac{k}{2^m}\right)=\sum_{n=0}^{m-1}\frac{\phi\!\left(2^n\frac{k}{2^m}\right)}{2^n},$$
und aus der unendlichen Reihe wird eine endliche Summe mit \(m\) Summanden. Die Implementierungen gehen noch weiter und drücken diesen Wert direkt durch eine Präfixsumme von Popcounts aus.
Sei \(s_2(r)\) die Anzahl der Einsen in der Binärdarstellung von \(r\), und
$$P(k)=\sum_{r=0}^{k-1}s_2(r).$$
Dann gilt an dyadischen Punkten die exakte Identität
$$B\left(\frac{k}{2^m}\right)=\frac{mk-2P(k)}{2^m}.$$
Das ist die zentrale Formel aller drei Implementierungen. Daraus folgt sofort die Rekurrenz
$$B\left(\frac{k+1}{2^m}\right)-B\left(\frac{k}{2^m}\right)=\frac{m-2s_2(k)}{2^m},$$
denn \(P(k+1)=P(k)+s_2(k)\). Ist also ein dyadischer Funktionswert bekannt, so erhält man den nächsten mit einem einzigen Bitcount und wenigen arithmetischen Operationen.
Der Kreis ist gegeben durch
$$\left(x-\frac14\right)^2+\left(y-\frac12\right)^2=\left(\frac14\right)^2.$$
Er liegt nur im horizontalen Bereich \(0 \le x \le \tfrac12\). Für festes \(x\) ist der vertikale Kreisschnitt
$$y_{\pm}(x)=\frac12 \pm \sqrt{\frac1{16}-\left(x-\frac14\right)^2}.$$
Die Blancmange-Region ist
$$R_B=\{(x,y): 0 \le y \le B(x)\}.$$
Die gemeinsame vertikale Höhe ist also
$$h(x)=\max\!\left(0,\ \min\!\bigl(B(x),y_+(x)\bigr)-\max\!\bigl(0,y_-(x)\bigr)\right).$$
Gesucht ist damit
$$A=\int_0^{1/2} h(x)\,dx.$$
Außerhalb von \([0,\tfrac12]\) gibt es keinen Beitrag, weil dort der Kreis gar nicht existiert.
Setze \(m=3\) und \(k=3\), also \(x=3/8\). Dann ist
$$P(3)=s_2(0)+s_2(1)+s_2(2)=0+1+1=2,$$
und somit
$$B\left(\frac38\right)=\frac{3 \cdot 3-2 \cdot 2}{8}=\frac58.$$
Das stimmt mit der direkten Takagi-Summe überein:
$$B\left(\frac38\right)=\phi\left(\frac38\right)+\frac12\phi\left(\frac34\right)+\frac14\phi\left(\frac32\right)=\frac38+\frac18+\frac18=\frac58.$$
Für den Kreisschnitt gilt \(x-\tfrac14=\tfrac18\), also
$$y_-\left(\frac38\right)=\frac12-\frac{\sqrt3}{8},\qquad y_+\left(\frac38\right)=\frac12+\frac{\sqrt3}{8}.$$
Da \(\tfrac58 < y_+(\tfrac38)\), ist die Überlappungshöhe
$$h\left(\frac38\right)=\frac58-\left(\frac12-\frac{\sqrt3}{8}\right)=\frac{1+\sqrt3}{8}.$$
Genau diese Größe misst der Integrand: den Anteil des vertikalen Kreissegments, der unter der Blancmange-Kurve und zugleich oberhalb der \(x\)-Achse liegt.
Auf Stufe \(m\) wird das Intervall \([0,\tfrac12]\) in \(N=2^{m-1}\) gleich große Teilintervalle der Breite \(h=2^{-m}\) zerlegt. Jeder Stützpunkt ist also dyadisch:
$$x_i=\frac{i}{2^m},\qquad 0 \le i \le 2^{m-1}.$$
Deshalb wird jeder Blancmange-Wert \(B(x_i)\) exakt über die obige Formel berechnet; numerisch approximiert wird nur noch das Integral. Die Simpson-Näherung lautet
$$A_m=\frac{h}{3}\left(h(x_0)+h(x_N)+4\sum_{\substack{1 \le i \le N-1 \\ i\text{ ungerade}}} h(x_i)+2\sum_{\substack{2 \le i \le N-2 \\ i\text{ gerade}}} h(x_i)\right).$$
Die Implementierungen verfeinern das dyadische Gitter so lange, bis sich zwei aufeinanderfolgende Stufen um weniger als die vorgegebene Toleranz unterscheiden.
Um \(P(k)\) schnell zu erhalten, zählen die Implementierungen die Einsen bitweise. Für eine feste Bitposition \(b\) wiederholt sich das Muster alle \(2^{b+1}\) Zahlen, und die Anzahl der Einsen unter \(0,1,\dots,k-1\) ist
$$\left\lfloor\frac{k}{2^{b+1}}\right\rfloor 2^b+\max\!\left(0,\ k \bmod 2^{b+1}-2^b\right).$$
Durch Summation über die Bitpositionen erhält man \(P(k)\) in \(O(\log k)\) statt durch eine Schleife über alle kleineren Zahlen.
Während einer festen Simpson-Stufe halten die Implementierungen den Zähler \(mk-2P(k)\) als Invariante fest und berechnen ihn nicht an jeder Stelle neu. Der Übergang von \(k\) nach \(k+1\) ändert diesen Zähler nur um \(m-2s_2(k)\). So lässt sich das gesamte dyadische Gitter in einem Durchlauf abarbeiten.
Für jeden Stützpunkt werden die beiden Kreisgrenzen \(y_-(x)\) und \(y_+(x)\) bestimmt, mit dem Intervall \([0,B(x)]\) geschnitten und die resultierende Höhe entweder zum ungeraden oder zum geraden Simpson-Anteil addiert. Der innere Bereich der Stützstellen wird in unabhängige Blöcke aufgeteilt, damit Teilsummen getrennt berechnet und am Ende zusammengeführt werden können. Die C++-, Python- und Java-Implementierungen benutzen dieselbe Mathematik; die kompilierten Versionen verteilen die Arbeit auf Threads, und die Python-Version kann dafür mehrere Prozesse einsetzen. Danach wird die nächste Verfeinerungsstufe berechnet und mit der vorherigen verglichen.
Auf Stufe \(m\) verwendet Simpson auf \([0,\tfrac12]\) genau \(2^{m-1}+1\) Stützpunkte, also kostet eine Verfeinerungsstufe \(O(2^m)\) Zeit. Durch die Rekurrenz für \(mk-2P(k)\) fällt pro zusätzlichem Punkt nur konstanter Rechenaufwand plus Bitcount und Quadratwurzel an.
Der Zusatzspeicher beträgt \(O(T)\) für \(T\) partielle Arbeitssummen, beziehungsweise \(O(1)\) in einer seriellen Variante. Da nur wenige aufeinanderfolgende dyadische Stufen geprüft werden, bleibt die Gesamtlaufzeit bis zur Konvergenz von derselben Größenordnung wie die feinste tatsächlich ausgewertete Stufe.
Blancmange ya da Takagi fonksiyonu,
$$B(x)=\sum_{n=0}^{\infty}\frac{\phi(2^n x)}{2^n},\qquad \phi(u)=\operatorname{dist}(u,\mathbb{Z})$$
şeklinde tanımlanan, sürekli ama hiçbir noktada türevlenebilir olmayan klasik eğridir. Problem 226, bu eğrinin altındaki bölge ile merkezi \(\left(\tfrac14,\tfrac12\right)\), yarıçapı \(\tfrac14\) olan çemberin ortak alanını ister. Çözümler sembolik bir kapalı form aramaz; geometriyi tek boyutlu bir integrale indirger ve blancmange fonksiyonunu örnekleme noktalarında tam olarak hesaplar.
Asıl püf noktası şudur: Simpson yöntemi yalnızca \(x=k/2^m\) biçimindeki dyadik noktaları örnekler. Blancmange fonksiyonu da tam bu noktalar üzerinde sonsuz seri olmaktan çıkıp ikili basamak toplamlarıyla yazılabilen sonlu ve çok hızlı bir ifadeye dönüşür.
\(x=k/2^m\) ise, \(n \ge m\) için \(2^n x\) bir tam sayı olur; dolayısıyla \(\phi(2^n x)=0\). Bu yüzden
$$B\left(\frac{k}{2^m}\right)=\sum_{n=0}^{m-1}\frac{\phi\!\left(2^n\frac{k}{2^m}\right)}{2^n}$$
olur ve sonsuz seri yalnızca \(m\) terimli bir toplam hâline gelir. Uygulamalar bununla da yetinmez; bu değeri doğrudan bit sayıları üzerinden hesaplar.
\(s_2(r)\), \(r\)'nin ikili yazımındaki 1 bitlerinin sayısı olsun. Ayrıca
$$P(k)=\sum_{r=0}^{k-1}s_2(r)$$
tanımını yapalım. O zaman dyadik noktalarda tam olarak
$$B\left(\frac{k}{2^m}\right)=\frac{mk-2P(k)}{2^m}$$
elde edilir. Üç çözümün de dayandığı ana eşitlik budur. Buradan hemen
$$B\left(\frac{k+1}{2^m}\right)-B\left(\frac{k}{2^m}\right)=\frac{m-2s_2(k)}{2^m}$$
rekürrensi gelir; çünkü \(P(k+1)=P(k)+s_2(k)\). Yani bir dyadik değer biliniyorsa, bir sonraki değere geçmek için sadece mevcut indeksin bit sayısını bilmek yeterlidir.
Çember denklemi
$$\left(x-\frac14\right)^2+\left(y-\frac12\right)^2=\left(\frac14\right)^2$$
şeklindedir. Bu çember yalnızca \(0 \le x \le \tfrac12\) aralığında görünür. Sabit bir \(x\) için düşey kesit sınırları
$$y_{\pm}(x)=\frac12 \pm \sqrt{\frac1{16}-\left(x-\frac14\right)^2}$$
olur. Blancmange bölgesi ise
$$R_B=\{(x,y): 0 \le y \le B(x)\}$$
olarak tanımlanır. Dolayısıyla ortak düşey yükseklik
$$h(x)=\max\!\left(0,\ \min\!\bigl(B(x),y_+(x)\bigr)-\max\!\bigl(0,y_-(x)\bigr)\right)$$
şeklindedir. Aranan alan da
$$A=\int_0^{1/2} h(x)\,dx$$
olur. \([0,\tfrac12]\) dışındaki \(x\) değerleri katkı yapmaz; çünkü o bölgede çember yoktur.
\(m=3\) ve \(k=3\) alalım; yani \(x=3/8\). Önce prefix bit toplamı
$$P(3)=s_2(0)+s_2(1)+s_2(2)=0+1+1=2$$
olur. Buradan
$$B\left(\frac38\right)=\frac{3 \cdot 3-2 \cdot 2}{8}=\frac58$$
elde edilir. Bu değer doğrudan seriyle de aynıdır:
$$B\left(\frac38\right)=\phi\left(\frac38\right)+\frac12\phi\left(\frac34\right)+\frac14\phi\left(\frac32\right)=\frac38+\frac18+\frac18=\frac58.$$
Çember tarafında \(x-\tfrac14=\tfrac18\) olduğundan
$$y_-\left(\frac38\right)=\frac12-\frac{\sqrt3}{8},\qquad y_+\left(\frac38\right)=\frac12+\frac{\sqrt3}{8}$$
bulunur. \(\tfrac58 < y_+(\tfrac38)\) olduğuna göre örtüşme yüksekliği
$$h\left(\frac38\right)=\frac58-\left(\frac12-\frac{\sqrt3}{8}\right)=\frac{1+\sqrt3}{8}$$
olur. İntegrandın her örnekleme noktasında ölçtüğü şey tam olarak budur.
Seviye \(m\)'de \([0,\tfrac12]\) aralığı genişliği \(h=2^{-m}\) olan \(N=2^{m-1}\) eşit alt aralığa bölünür. Böylece her düğüm dyadiktir:
$$x_i=\frac{i}{2^m},\qquad 0 \le i \le 2^{m-1}.$$
Bu yüzden örneklenen her \(B(x_i)\) değeri yukarıdaki formülle tam hesaplanır; yaklaşık olan tek şey integralin kendisidir. Simpson yaklaşımı
$$A_m=\frac{h}{3}\left(h(x_0)+h(x_N)+4\sum_{\substack{1 \le i \le N-1 \\ i\text{ tek}}} h(x_i)+2\sum_{\substack{2 \le i \le N-2 \\ i\text{ çift}}} h(x_i)\right)$$
şeklindedir. Uygulamalar giderek daha ince dyadik ızgaralara geçer ve ardışık iki seviye arasındaki fark toleransın altına indiğinde durur.
\(P(k)\) değerini hızlı almak için çözümler 1 bitlerini tek tek sayılar üzerinde değil, bit konumları üzerinden sayar. Sabit bir \(b\) biti için desen her \(2^{b+1}\) sayıda tekrar eder ve \(0,1,\dots,k-1\) aralığındaki 1 sayısı
$$\left\lfloor\frac{k}{2^{b+1}}\right\rfloor 2^b+\max\!\left(0,\ k \bmod 2^{b+1}-2^b\right)$$
olur. Bütün bit konumları toplanınca \(P(k)\), tüm önceki sayıları tek tek dolaşmadan \(O(\log k)\) zamanda elde edilir.
Tek bir Simpson taramasında uygulamalar her düğümde \(mk-2P(k)\) ifadesini baştan kurmaz. Bunun yerine bu payı bir değişmez olarak taşır; \(k\)'den \(k+1\)'e geçerken artış yalnızca \(m-2s_2(k)\) kadar olur. Böylece çok ince bir dyadik ızgara bile tek geçişte işlenebilir.
Her örnekleme noktasında çemberin alt ve üst \(y\) sınırları hesaplanır, bu aralık \([0,B(x)]\) ile kesiştirilir ve çıkan yükseklik Simpson toplamındaki tek ya da çift katsayılı haneye eklenir. İç düğümler bağımsız parçalara bölünebildiği için kısmi toplamlar ayrı ayrı hesaplanıp sonunda birleştirilebilir. C++, Python ve Java uygulamaları aynı matematiği izler; derlenen sürümler işi iş parçacıklarına bölerek hızlandırır, Python sürümü ise uygun olduğunda süreç tabanlı paralellik kullanabilir. Sonra bir sonraki seviye hesaplanır ve bir öncekiyle karşılaştırılır.
Seviye \(m\)'de Simpson yöntemi \([0,\tfrac12]\) aralığında \(2^{m-1}+1\) örnekleme noktası kullanır; dolayısıyla bir seviye \(O(2^m)\) zamanda hesaplanır. Dyadik blancmange değerleri rekürrens ile aktarıldığı için her yeni düğüm yalnızca sabit sayıda aritmetik işlem, bir bit sayımı ve bir karekök gerektirir.
Ek bellek, \(T\) adet işçi için \(O(T)\) kısmi toplam kadardır; seri çalışmada bu \(O(1)\)'e iner. İncelenebilecek seviyeler kısa bir artan dizi oluşturduğu için toplam süre de esasen en ince değerlendirilen seviye ile aynı mertebededir.
La función blancmange, también llamada función de Takagi, es la curva continua y no derivable
$$B(x)=\sum_{n=0}^{\infty}\frac{\phi(2^n x)}{2^n},\qquad \phi(u)=\operatorname{dist}(u,\mathbb{Z}).$$
En el problema 226 se pide el área común entre la región situada bajo esta curva y el círculo de centro \(\left(\tfrac14,\tfrac12\right)\) y radio \(\tfrac14\). Las implementaciones no intentan encontrar una primitiva cerrada. En lugar de eso, convierten la geometría en una integral unidimensional y evalúan exactamente la función blancmange en los puntos muestreados aprovechando la estructura diádica del problema.
La idea central es que la regla de Simpson utiliza puntos de la forma \(x=k/2^m\). Justamente en esos puntos la serie blancmange se reduce a una fórmula finita basada en sumas de dígitos binarios, lo que hace viable una integración numérica muy fina.
Si \(x=k/2^m\), entonces \(2^n x\) es entero para todo \(n \ge m\), y por tanto \(\phi(2^n x)=0\) a partir de ese índice. Así,
$$B\left(\frac{k}{2^m}\right)=\sum_{n=0}^{m-1}\frac{\phi\!\left(2^n\frac{k}{2^m}\right)}{2^n},$$
de modo que la serie infinita pasa a ser una suma finita de \(m\) términos. Las implementaciones todavía simplifican más el cálculo y lo expresan mediante una suma prefija de popcounts.
Sea \(s_2(r)\) el número de bits iguales a 1 en la expansión binaria de \(r\), y definamos
$$P(k)=\sum_{r=0}^{k-1}s_2(r).$$
Entonces, en puntos diádicos se cumple la identidad exacta
$$B\left(\frac{k}{2^m}\right)=\frac{mk-2P(k)}{2^m}.$$
Esa es la fórmula decisiva en las tres implementaciones. Además produce la recurrencia
$$B\left(\frac{k+1}{2^m}\right)-B\left(\frac{k}{2^m}\right)=\frac{m-2s_2(k)}{2^m},$$
porque \(P(k+1)=P(k)+s_2(k)\). Una vez conocido un valor diádico, el siguiente se obtiene con una sola cuenta de bits.
El círculo viene dado por
$$\left(x-\frac14\right)^2+\left(y-\frac12\right)^2=\left(\frac14\right)^2.$$
Su soporte horizontal es solo \(0 \le x \le \tfrac12\). Para un \(x\) fijo, el corte vertical del círculo es
$$y_{\pm}(x)=\frac12 \pm \sqrt{\frac1{16}-\left(x-\frac14\right)^2}.$$
La región blancmange es
$$R_B=\{(x,y): 0 \le y \le B(x)\}.$$
Por tanto, la altura común en la columna vertical de abscisa \(x\) es
$$h(x)=\max\!\left(0,\ \min\!\bigl(B(x),y_+(x)\bigr)-\max\!\bigl(0,y_-(x)\bigr)\right).$$
El área pedida queda reducida a
$$A=\int_0^{1/2} h(x)\,dx.$$
Fuera del intervalo \([0,\tfrac12]\) no hay contribución, porque allí el círculo no aparece.
Tomemos \(m=3\) y \(k=3\), es decir, \(x=3/8\). Entonces
$$P(3)=s_2(0)+s_2(1)+s_2(2)=0+1+1=2,$$
y de ahí
$$B\left(\frac38\right)=\frac{3 \cdot 3-2 \cdot 2}{8}=\frac58.$$
Esto coincide con la suma directa de Takagi:
$$B\left(\frac38\right)=\phi\left(\frac38\right)+\frac12\phi\left(\frac34\right)+\frac14\phi\left(\frac32\right)=\frac38+\frac18+\frac18=\frac58.$$
En el círculo, \(x-\tfrac14=\tfrac18\), luego
$$y_-\left(\frac38\right)=\frac12-\frac{\sqrt3}{8},\qquad y_+\left(\frac38\right)=\frac12+\frac{\sqrt3}{8}.$$
Como \(\tfrac58 < y_+(\tfrac38)\), la altura de solape vale
$$h\left(\frac38\right)=\frac58-\left(\frac12-\frac{\sqrt3}{8}\right)=\frac{1+\sqrt3}{8}.$$
Ese ejemplo muestra con claridad qué representa el integrando: la parte del segmento vertical del círculo que queda por debajo de la curva blancmange y por encima del eje \(x\).
En el nivel \(m\), el intervalo \([0,\tfrac12]\) se divide en \(N=2^{m-1}\) subintervalos iguales de anchura \(h=2^{-m}\), de modo que todos los nodos son diádicos:
$$x_i=\frac{i}{2^m},\qquad 0 \le i \le 2^{m-1}.$$
Por eso cada valor \(B(x_i)\) se calcula exactamente con la fórmula anterior; la única aproximación numérica está en la cuadratura. La regla de Simpson produce
$$A_m=\frac{h}{3}\left(h(x_0)+h(x_N)+4\sum_{\substack{1 \le i \le N-1 \\ i\text{ impar}}} h(x_i)+2\sum_{\substack{2 \le i \le N-2 \\ i\text{ par}}} h(x_i)\right).$$
Las implementaciones refinan la malla diádica y se detienen cuando dos niveles consecutivos difieren menos que la tolerancia requerida.
Para obtener \(P(k)\) con rapidez, las implementaciones cuentan los unos bit a bit. En una posición \(b\), el patrón se repite cada \(2^{b+1}\) enteros, y el número de unos entre \(0,1,\dots,k-1\) es
$$\left\lfloor\frac{k}{2^{b+1}}\right\rfloor 2^b+\max\!\left(0,\ k \bmod 2^{b+1}-2^b\right).$$
Al sumar sobre las posiciones de bit, \(P(k)\) se obtiene en \(O(\log k)\) tiempo aritmético.
Durante una pasada de Simpson a nivel fijo, la implementación mantiene como invariante el numerador \(mk-2P(k)\) y no lo reconstruye desde cero en cada nodo. Al pasar de \(k\) a \(k+1\), ese numerador solo cambia en \(m-2s_2(k)\). Gracias a eso, toda la malla fina puede recorrerse en una sola pasada lineal.
En cada nodo se calculan las cotas inferior y superior del círculo, se recorta ese intervalo con \([0,B(x)]\) y la altura resultante se añade a la suma impar o a la suma par de Simpson. El rango de nodos interiores se divide en bloques independientes, de modo que las sumas parciales se puedan combinar al final. Las implementaciones en C++, Python y Java siguen exactamente la misma matemática; las versiones compiladas reparten los bloques entre hilos, y la versión en Python puede usar varios procesos cuando compensa. Al terminar un nivel, se calcula el siguiente y se comparan ambos resultados.
En el nivel \(m\), la regla de Simpson sobre \([0,\tfrac12]\) usa \(2^{m-1}+1\) puntos de muestra, así que cada nivel cuesta \(O(2^m)\) tiempo. Como los valores diádicos de la función blancmange se van actualizando con la recurrencia anterior, cada nuevo nodo requiere solo una cantidad constante de aritmética, más un conteo de bits y una raíz cuadrada.
La memoria adicional es \(O(T)\) para \(T\) sumas parciales de trabajadores, u \(O(1)\) en la versión secuencial. Como solo se recorre una corta escalera de niveles cada vez más finos, el tiempo total hasta converger sigue siendo del orden del nivel más fino evaluado.
Blancmange 曲线，也就是 Takagi 函数，定义为
$$B(x)=\sum_{n=0}^{\infty}\frac{\phi(2^n x)}{2^n},\qquad \phi(u)=\operatorname{dist}(u,\mathbb{Z}).$$
第 226 题要求的是这条曲线下方区域与一个圆的公共面积；该圆的圆心是 \(\left(\tfrac14,\tfrac12\right)\)，半径是 \(\tfrac14\)。实现并没有去寻找闭式积分，而是把几何问题改写成一个一维积分，再利用 dyadic 有理点上的特殊结构，把每个采样点处的 blancmange 函数值精确算出来。
整个方法的核心在于：Simpson 求积会在 \(x=k/2^m\) 这样的点上取样，而 blancmange 函数恰好在这些点上可以从无穷级数坍缩成有限的二进制位和公式。因此，最后的数值误差只来自积分公式本身，而不来自对曲线取值的近似。
若 \(x=k/2^m\)，那么当 \(n \ge m\) 时，\(2^n x\) 已经是整数，所以 \(\phi(2^n x)=0\)。于是
$$B\left(\frac{k}{2^m}\right)=\sum_{n=0}^{m-1}\frac{\phi\!\left(2^n\frac{k}{2^m}\right)}{2^n},$$
原本的无穷级数就缩成了只有 \(m\) 项的有限和。实现还继续向前推进，用一个与二进制 1 的个数有关的公式直接计算这个值，从而避免在每个采样点都重新把这 \(m\) 项相加。
记 \(s_2(r)\) 为整数 \(r\) 的二进制表示中 1 的个数，再定义前缀和
$$P(k)=\sum_{r=0}^{k-1}s_2(r).$$
那么在 dyadic 点上有精确恒等式
$$B\left(\frac{k}{2^m}\right)=\frac{mk-2P(k)}{2^m}.$$
这就是三份实现共同依赖的主公式。它立刻给出递推关系
$$B\left(\frac{k+1}{2^m}\right)-B\left(\frac{k}{2^m}\right)=\frac{m-2s_2(k)}{2^m},$$
因为 \(P(k+1)=P(k)+s_2(k)\)。也就是说，只要知道某个 dyadic 点的值，沿着网格向右走一步时，只需再做一次 bit count 就能得到下一个点的函数值。
题目中的圆满足
$$\left(x-\frac14\right)^2+\left(y-\frac12\right)^2=\left(\frac14\right)^2.$$
它只出现在水平区间 \(0 \le x \le \tfrac12\) 上。对固定的 \(x\)，圆的纵向截线边界是
$$y_{\pm}(x)=\frac12 \pm \sqrt{\frac1{16}-\left(x-\frac14\right)^2}.$$
而 blancmange 区域就是
$$R_B=\{(x,y): 0 \le y \le B(x)\}.$$
因此，在这个 \(x\) 处真正贡献面积的纵向重叠高度为
$$h(x)=\max\!\left(0,\ \min\!\bigl(B(x),y_+(x)\bigr)-\max\!\bigl(0,y_-(x)\bigr)\right).$$
要求的面积便化为
$$A=\int_0^{1/2} h(x)\,dx.$$
区间 \([0,\tfrac12]\) 之外没有贡献，因为那里的圆截线根本不存在。
取 \(m=3\)、\(k=3\)，于是 \(x=3/8\)。先算前缀 bit 数之和：
$$P(3)=s_2(0)+s_2(1)+s_2(2)=0+1+1=2.$$
带入主公式得到
$$B\left(\frac38\right)=\frac{3 \cdot 3-2 \cdot 2}{8}=\frac58.$$
直接从 Takagi 级数去算也完全一致：
$$B\left(\frac38\right)=\phi\left(\frac38\right)+\frac12\phi\left(\frac34\right)+\frac14\phi\left(\frac32\right)=\frac38+\frac18+\frac18=\frac58.$$
对圆来说，此时 \(x-\tfrac14=\tfrac18\)，所以
$$y_-\left(\frac38\right)=\frac12-\frac{\sqrt3}{8},\qquad y_+\left(\frac38\right)=\frac12+\frac{\sqrt3}{8}.$$
由于 \(\tfrac58 < y_+(\tfrac38)\)，纵向重叠高度就是
$$h\left(\frac38\right)=\frac58-\left(\frac12-\frac{\sqrt3}{8}\right)=\frac{1+\sqrt3}{8}.$$
这个例子把积分函数的含义说得很直观：在每个采样点上，它测量的是“圆的这条竖直截线中，有多少长度既落在圆内，又落在 blancmange 曲线下方”。
在细化层级 \(m\) 上，区间 \([0,\tfrac12]\) 被等分成 \(N=2^{m-1}\) 个小区间，步长为 \(h=2^{-m}\)，因此所有节点都是 dyadic 点：
$$x_i=\frac{i}{2^m},\qquad 0 \le i \le 2^{m-1}.$$
这意味着每个被采样到的 \(B(x_i)\) 都能用上面的公式精确得到，数值近似只发生在积分本身。Simpson 公式写成
$$A_m=\frac{h}{3}\left(h(x_0)+h(x_N)+4\sum_{\substack{1 \le i \le N-1 \\ i\text{ 为奇数}}} h(x_i)+2\sum_{\substack{2 \le i \le N-2 \\ i\text{ 为偶数}}} h(x_i)\right).$$
实现会逐层细化这个 dyadic 网格，并在相邻两层结果之差小于给定容差时停止。
为了快速得到 \(P(k)\)，实现不是把 \(0\) 到 \(k-1\) 全部枚举一遍，而是按 bit 位置分别统计。对某个固定的 bit 位置 \(b\)，1 的模式每 \(2^{b+1}\) 个整数循环一次，因此在 \(0,1,\dots,k-1\) 中，这一位为 1 的次数是
$$\left\lfloor\frac{k}{2^{b+1}}\right\rfloor 2^b+\max\!\left(0,\ k \bmod 2^{b+1}-2^b\right).$$
把所有 bit 位置的贡献相加，就能在 \(O(\log k)\) 的算术时间内得到 \(P(k)\)。
在固定层级的一次 Simpson 扫描中，实现维护的是分子 \(mk-2P(k)\) 这个不变量，而不是在每个节点重新计算 \(B(k/2^m)\)。当索引从 \(k\) 走到 \(k+1\) 时，分子只会增加 \(m-2s_2(k)\)。因此，只要初始化好第一个节点，后面整条 dyadic 网格都可以线性流过。
对每个采样点，代码先求出圆的上下边界，再与区间 \([0,B(x)]\) 相交，把得到的重叠高度加到 Simpson 的奇数项或偶数项累加器中。内部采样点可以自然拆成互不依赖的多个区块，因此部分和能够独立计算后再合并。C++、Python 和 Java 三个实现遵循的是同一套数学结构；编译型实现把这些区块分配给多个线程，Python 实现在合适时也可以把区块交给多个工作进程。每完成一层，就继续计算下一层，并比较两者差值是否已经低于容差。
在层级 \(m\) 上，Simpson 规则对 \([0,\tfrac12]\) 使用 \(2^{m-1}+1\) 个采样点，因此单层时间复杂度为 \(O(2^m)\)。由于 blancmange 的 dyadic 取值是用递推流式更新的，所以每多一个节点只需要常数次算术、一次 bit count 和一次平方根计算。
额外空间复杂度是 \(O(T)\)，其中 \(T\) 是并行工作者保存的部分和数量；若串行执行则可视为 \(O(1)\)。因为实际只会检查一小串逐步细化的层级，总运行时间仍与最终评估到的最细层级同阶。
Функция Бланманже, или функция Такэги, задается формулой
$$B(x)=\sum_{n=0}^{\infty}\frac{\phi(2^n x)}{2^n},\qquad \phi(u)=\operatorname{dist}(u,\mathbb{Z}).$$
В задаче 226 требуется найти площадь общей части области под этой кривой и круга с центром \(\left(\tfrac14,\tfrac12\right)\) и радиусом \(\tfrac14\). Реализации не ищут замкнутую формулу для площади. Вместо этого геометрия сводится к одномерному интегралу, а значения функции Бланманже в точках выборки вычисляются точно благодаря структуре диадических рациональных чисел.
Главное наблюдение состоит в том, что правило Симпсона использует точки вида \(x=k/2^m\). Именно на таких точках бесконечный ряд для функции Бланманже превращается в конечную формулу через суммы двоичных единиц, и потому вычисление на очень тонкой сетке остается быстрым.
Если \(x=k/2^m\), то для всех \(n \ge m\) число \(2^n x\) уже является целым, а значит \(\phi(2^n x)=0\). Следовательно,
$$B\left(\frac{k}{2^m}\right)=\sum_{n=0}^{m-1}\frac{\phi\!\left(2^n\frac{k}{2^m}\right)}{2^n},$$
и бесконечный ряд сжимается до конечной суммы из \(m\) членов. В реализациях делается еще один шаг: это значение выражается прямо через префиксную сумму popcount.
Пусть \(s_2(r)\) обозначает количество единиц в двоичной записи числа \(r\), а
$$P(k)=\sum_{r=0}^{k-1}s_2(r).$$
Тогда в диадических точках выполняется точное равенство
$$B\left(\frac{k}{2^m}\right)=\frac{mk-2P(k)}{2^m}.$$
Именно это тождество лежит в основе всех трех реализаций. Из него сразу получается рекуррентное соотношение
$$B\left(\frac{k+1}{2^m}\right)-B\left(\frac{k}{2^m}\right)=\frac{m-2s_2(k)}{2^m},$$
поскольку \(P(k+1)=P(k)+s_2(k)\). Поэтому, зная одно значение на сетке, следующее можно получить за постоянное число операций.
Круг описывается уравнением
$$\left(x-\frac14\right)^2+\left(y-\frac12\right)^2=\left(\frac14\right)^2.$$
Его горизонтальная проекция ограничена интервалом \(0 \le x \le \tfrac12\). Для фиксированного \(x\) верхняя и нижняя границы вертикального сечения круга равны
$$y_{\pm}(x)=\frac12 \pm \sqrt{\frac1{16}-\left(x-\frac14\right)^2}.$$
Область под функцией Бланманже есть
$$R_B=\{(x,y): 0 \le y \le B(x)\}.$$
Значит, общая вертикальная высота в точке \(x\) равна
$$h(x)=\max\!\left(0,\ \min\!\bigl(B(x),y_+(x)\bigr)-\max\!\bigl(0,y_-(x)\bigr)\right).$$
Искомая площадь принимает вид
$$A=\int_0^{1/2} h(x)\,dx.$$
Вне интервала \([0,\tfrac12]\) вклад нулевой, потому что там круга просто нет.
Возьмем \(m=3\) и \(k=3\), то есть \(x=3/8\). Тогда
$$P(3)=s_2(0)+s_2(1)+s_2(2)=0+1+1=2,$$
и потому
$$B\left(\frac38\right)=\frac{3 \cdot 3-2 \cdot 2}{8}=\frac58.$$
Это полностью совпадает с прямым суммированием ряда Такэги:
$$B\left(\frac38\right)=\phi\left(\frac38\right)+\frac12\phi\left(\frac34\right)+\frac14\phi\left(\frac32\right)=\frac38+\frac18+\frac18=\frac58.$$
Для круга имеем \(x-\tfrac14=\tfrac18\), следовательно
$$y_-\left(\frac38\right)=\frac12-\frac{\sqrt3}{8},\qquad y_+\left(\frac38\right)=\frac12+\frac{\sqrt3}{8}.$$
Поскольку \(\tfrac58 < y_+(\tfrac38)\), получаем высоту перекрытия
$$h\left(\frac38\right)=\frac58-\left(\frac12-\frac{\sqrt3}{8}\right)=\frac{1+\sqrt3}{8}.$$
Этот пример наглядно показывает смысл интегранта: он измеряет длину той части вертикального сечения круга, которая лежит ниже кривой Бланманже и выше оси \(x\).
На уровне \(m\) интервал \([0,\tfrac12]\) разбивается на \(N=2^{m-1}\) равных частей ширины \(h=2^{-m}\), поэтому все узлы являются диадическими:
$$x_i=\frac{i}{2^m},\qquad 0 \le i \le 2^{m-1}.$$
Значит, каждое значение \(B(x_i)\) вычисляется точно по приведенной формуле, и численная погрешность возникает только из квадратурной аппроксимации. Формула Симпсона имеет вид
$$A_m=\frac{h}{3}\left(h(x_0)+h(x_N)+4\sum_{\substack{1 \le i \le N-1 \\ i\text{ нечетное}}} h(x_i)+2\sum_{\substack{2 \le i \le N-2 \\ i\text{ четное}}} h(x_i)\right).$$
Реализации последовательно уточняют диадическую сетку и останавливаются, когда результаты двух соседних уровней различаются меньше заданной точности.
Чтобы быстро вычислить \(P(k)\), реализации считают единицы не по всем числам подряд, а по разрядам. Для фиксированного разряда \(b\) шаблон повторяется каждые \(2^{b+1}\) чисел, а число единиц среди \(0,1,\dots,k-1\) равно
$$\left\lfloor\frac{k}{2^{b+1}}\right\rfloor 2^b+\max\!\left(0,\ k \bmod 2^{b+1}-2^b\right).$$
Суммирование по разрядам дает \(P(k)\) за \(O(\log k)\) арифметических операций.
Внутри одного прохода Симпсона на фиксированном уровне код поддерживает инвариант \(mk-2P(k)\) и не пересчитывает его заново для каждого узла. При переходе от \(k\) к \(k+1\) этот числитель меняется лишь на \(m-2s_2(k)\). Благодаря этому вся тонкая сетка обрабатывается линейным проходом.
В каждом узле вычисляются нижняя и верхняя границы круга, затем этот отрезок пересекается с \([0,B(x)]\), и полученная высота добавляется либо к нечетной, либо к четной сумме Симпсона. Внутренние узлы естественным образом разбиваются на независимые блоки, поэтому частичные суммы можно считать отдельно и объединять в конце. Реализации на C++, Python и Java используют одну и ту же математику; компилируемые версии распределяют блоки между потоками, а версия на Python может использовать рабочие процессы. После завершения очередного уровня вычисляется следующий, и оба результата сравниваются.
На уровне \(m\) правило Симпсона на \([0,\tfrac12]\) использует \(2^{m-1}+1\) узлов, так что стоимость одного уровня равна \(O(2^m)\). Поскольку значения функции Бланманже на сетке обновляются по рекуррентной формуле, каждый новый узел требует только постоянного числа арифметических операций, одного bit count и одного извлечения квадратного корня.
Дополнительная память равна \(O(T)\) для \(T\) частичных сумм работников, либо \(O(1)\) в последовательном режиме. Так как проверяется лишь короткая цепочка все более тонких диадических уровней, суммарное время до сходимости остается порядка последнего, самого тонкого уровня.
دالة Blancmange، أو دالة Takagi، تُعرَّف بالمتسلسلة
$$B(x)=\sum_{n=0}^{\infty}\frac{\phi(2^n x)}{2^n},\qquad \phi(u)=\operatorname{dist}(u,\mathbb{Z}).$$
في المسألة 226 نريد مساحة الجزء المشترك بين المنطقة الواقعة تحت هذا المنحنى وبين دائرة مركزها \(\left(\tfrac14,\tfrac12\right)\) ونصف قطرها \(\tfrac14\). لا تبحث الحلول عن تكامل رمزي مغلق، بل تحوِّل المسألة الهندسية إلى تكامل أحادي البعد، ثم تستفيد من البنية الثنائية للنقاط العيّنية لتحسب قيم منحنى Blancmange بدقة تامة عند كل نقطة تُستعمل في التكامل.
الفكرة الحاسمة هي أن قاعدة Simpson تأخذ عيناتها عند نقاط من الشكل \(x=k/2^m\). وعند هذه النقاط بالذات تنكمش متسلسلة Blancmange اللانهائية إلى صيغة منتهية مكتوبة بدلالة مجموع أعداد الواحدات في التمثيل الثنائي، ولهذا تصبح عملية التكامل العددي دقيقة وفعالة في الوقت نفسه.
إذا كان \(x=k/2^m\)، فإن \(2^n x\) يصبح عددًا صحيحًا لكل \(n \ge m\)، وبالتالي تكون \(\phi(2^n x)=0\) ابتداءً من ذلك الحد. لذلك نحصل على
$$B\left(\frac{k}{2^m}\right)=\sum_{n=0}^{m-1}\frac{\phi\!\left(2^n\frac{k}{2^m}\right)}{2^n},$$
فتتحول المتسلسلة اللانهائية إلى مجموع من \(m\) حدود فقط. ثم تذهب التطبيقات أبعد من ذلك، فتستبدل هذا المجموع نفسه بصيغة مباشرة تعتمد على مجموع popcount التراكمي.
لنرمز بـ \(s_2(r)\) إلى عدد الواحدات في التمثيل الثنائي للعدد \(r\)، ولنعرّف
$$P(k)=\sum_{r=0}^{k-1}s_2(r).$$
عندئذٍ تتحقق على النقاط الثنائية الصيغة الدقيقة
$$B\left(\frac{k}{2^m}\right)=\frac{mk-2P(k)}{2^m}.$$
وهذه هي العلاقة الأساسية التي تعتمد عليها جميع التطبيقات. ومنها نحصل مباشرة على العلاقة التراجعية
$$B\left(\frac{k+1}{2^m}\right)-B\left(\frac{k}{2^m}\right)=\frac{m-2s_2(k)}{2^m},$$
لأن \(P(k+1)=P(k)+s_2(k)\). وهذا يعني أن معرفة قيمة واحدة على الشبكة تكفي للانتقال إلى القيمة التالية بعملية ثابتة الكلفة تعتمد فقط على عدد البتات المساوية لـ 1 في الفهرس الحالي.
معادلة الدائرة هي
$$\left(x-\frac14\right)^2+\left(y-\frac12\right)^2=\left(\frac14\right)^2.$$
وهذه الدائرة لا تظهر أفقيًا إلا على المجال \(0 \le x \le \tfrac12\). وعند تثبيت قيمة \(x\)، يكون المقطع العمودي للدائرة محصورًا بين
$$y_{\pm}(x)=\frac12 \pm \sqrt{\frac1{16}-\left(x-\frac14\right)^2}.$$
أما منطقة Blancmange فهي
$$R_B=\{(x,y): 0 \le y \le B(x)\}.$$
إذن فإن الارتفاع العمودي المشترك عند الموضع \(x\) يساوي
$$h(x)=\max\!\left(0,\ \min\!\bigl(B(x),y_+(x)\bigr)-\max\!\bigl(0,y_-(x)\bigr)\right).$$
ومن ثم تصبح المساحة المطلوبة
$$A=\int_0^{1/2} h(x)\,dx.$$
ولا يوجد أي إسهام خارج \([0,\tfrac12]\)، لأن الدائرة نفسها لا تغطي ذلك الجزء من المحور الأفقي.
لنأخذ \(m=3\) و\(k=3\)، أي \(x=3/8\). عندئذٍ
$$P(3)=s_2(0)+s_2(1)+s_2(2)=0+1+1=2,$$
ومن ثم
$$B\left(\frac38\right)=\frac{3 \cdot 3-2 \cdot 2}{8}=\frac58.$$
وهذا يطابق مباشرة جمع متسلسلة Takagi:
$$B\left(\frac38\right)=\phi\left(\frac38\right)+\frac12\phi\left(\frac34\right)+\frac14\phi\left(\frac32\right)=\frac38+\frac18+\frac18=\frac58.$$
وبالنسبة إلى الدائرة فإن \(x-\tfrac14=\tfrac18\)، لذا
$$y_-\left(\frac38\right)=\frac12-\frac{\sqrt3}{8},\qquad y_+\left(\frac38\right)=\frac12+\frac{\sqrt3}{8}.$$
وبما أن \(\tfrac58 < y_+(\tfrac38)\)، فإن ارتفاع التداخل يصبح
$$h\left(\frac38\right)=\frac58-\left(\frac12-\frac{\sqrt3}{8}\right)=\frac{1+\sqrt3}{8}.$$
وهذا المثال يوضح بصورة ملموسة ماذا يقيس integrand: طول الجزء من المقطع العمودي للدائرة الذي يقع تحت منحنى Blancmange وفوق المحور \(x\).
عند مستوى التنعيم \(m\)، يُقسَّم المجال \([0,\tfrac12]\) إلى \(N=2^{m-1}\) فترات متساوية طول كل منها \(h=2^{-m}\)، ولذلك تكون جميع العقد ثنائية الشكل:
$$x_i=\frac{i}{2^m},\qquad 0 \le i \le 2^{m-1}.$$
وبما أن كل عقدة هي نقطة ثنائية، فإن قيم \(B(x_i)\) تُحسب بدقة كاملة من الصيغة السابقة، وتبقى التقريبية الوحيدة في خطوة التكامل نفسها. وتكون صيغة Simpson
$$A_m=\frac{h}{3}\left(h(x_0)+h(x_N)+4\sum_{\substack{1 \le i \le N-1 \\ i\text{ odd}}} h(x_i)+2\sum_{\substack{2 \le i \le N-2 \\ i\text{ even}}} h(x_i)\right).$$
تقوم التطبيقات بحساب هذه القيمة على شبكات ثنائية أدق فأدق، وتتوقف عندما يصبح الفرق بين مستويين متتاليين أصغر من حد الدقة المطلوب.
لحساب \(P(k)\) بسرعة، لا تمر التطبيقات على الأعداد من \(0\) إلى \(k-1\) واحدًا واحدًا، بل تحسب المساهمات بتقسيمها حسب موضع البت. فعند موضع بت ثابت \(b\)، يتكرر النمط كل \(2^{b+1}\) عددًا، وعدد المرات التي تكون فيها هذه البتة مساوية لـ 1 بين \(0,1,\dots,k-1\) هو
$$\left\lfloor\frac{k}{2^{b+1}}\right\rfloor 2^b+\max\!\left(0,\ k \bmod 2^{b+1}-2^b\right).$$
وبجمع هذه المساهمات على جميع المواضع نحصل على \(P(k)\) في زمن حسابي \(O(\log k)\).
أثناء مرور Simpson عند مستوى ثابت، تحتفظ التطبيقات بالمقدار \(mk-2P(k)\) بوصفه ثابتًا بنيويًا، بدلًا من إعادة حسابه من البداية عند كل عقدة. وعند الانتقال من \(k\) إلى \(k+1\)، يتغير هذا المقدار فقط بمقدار \(m-2s_2(k)\). ولهذا يمكن مسح شبكة ثنائية دقيقة جدًا في مرور خطي واحد.
في كل عقدة تُحسب الحافتان السفلية والعليا للدائرة، ثم يُقصّ هذا المجال مع \([0,B(x)]\)، ويُضاف ارتفاع التداخل الناتج إلى مجمّع الحدود الفردية أو الزوجية في صيغة Simpson. ويمكن تقسيم مجال العقد الداخلية إلى مقاطع مستقلة تمامًا، ولذلك تُحسب المجاميع الجزئية منفصلة ثم تُدمج في النهاية. تتبع تطبيقات C++ وPython وJava الصياغة الرياضية نفسها؛ فالنسخ المترجمة توزع المقاطع على خيوط متعددة، بينما يمكن لنسخة Python أن تستخدم عمليات عاملة متعددة عند الحاجة. وبعد إنهاء مستوى ما، يُحسب المستوى التالي وتُقارن القيمتان للحكم على التقارب.
في المستوى \(m\)، تستعمل قاعدة Simpson على المجال \([0,\tfrac12]\) عددًا مقداره \(2^{m-1}+1\) من نقاط العينة، ولذلك تكون كلفة المستوى الواحد \(O(2^m)\) زمنًا. وبسبب استعمال العلاقة التراجعية لتحديث قيم Blancmange على الشبكة الثنائية، فإن كل عقدة إضافية تحتاج فقط إلى عدد ثابت من العمليات الحسابية مع bit count واحد وجذر تربيعي واحد.
أما الذاكرة الإضافية فهي \(O(T)\) لعدد \(T\) من المجاميع الجزئية الخاصة بالعاملين، أو \(O(1)\) إذا نُفذ الحساب تسلسليًا. وبما أن التطبيق يفحص عددًا صغيرًا من مستويات التنعيم المتزايدة فقط، فإن الزمن الكلي حتى التقارب يبقى من رتبة المستوى الأدق الذي جرى تقييمه فعليًا.