Problem Summary

For each positive integer \(n\), let \(T(n)\) be the number of positive integer triples \((a,b,c)\) with \(b \le n\) that satisfy the necklace condition from the problem statement. The direct geometric search is far too large for \(n=10^9\), so the solution used by the implementations turns the geometry into a divisor-counting problem and then evaluates that arithmetic description with cached summatory functions.

Mathematical Approach

The three implementations all start from the same arithmetic reduction of the necklace condition. The relevant chain lengths are \(k=3\), \(k=4\), and \(k=6\), leading to the Diophantine families

$$k=3:\ (a-3b)(c-3b)=12b^2,$$

$$k=4:\ (a-b)(c-b)=2b^2,$$

$$k=6:\ (3a-b)(3c-b)=4b^2.$$

So the geometric configuration is reduced to counting admissible factorizations of \(2b^2\), \(4b^2\), \(12b^2\), and \(36b^2\), together with the local congruence restrictions forced by the primes \(2\) and \(3\).

Step 1: Separate the Squarefree Part Away from \(2\) and \(3\)

For primes \(p \ge 5\), only the squarefree divisor pattern matters. Define

$$u(m)=\sum_{\substack{d \mid m\\ \mu^2(d)=1\\ (d,6)=1}} 1,$$

the number of squarefree divisors of \(m\) that are coprime to \(6\). Its summatory function is

$$B(x)=\sum_{m \le x} u(m)=\sum_{\substack{d \le x\\ \mu^2(d)=1\\ (d,6)=1}}\left\lfloor\frac{x}{d}\right\rfloor.$$

This function is the basic building block of the whole computation. Intuitively, it captures the contribution of all primes other than \(2\) and \(3\), while the local behavior at those two exceptional primes is handled separately.

Step 2: Count Squarefree Integers Coprime to \(6\)

To evaluate \(B(x)\), the implementations first count squarefree numbers by Möbius inversion:

$$Q(x)=\sum_{k \le \sqrt{x}} \mu(k)\left\lfloor\frac{x}{k^2}\right\rfloor.$$

Now let \(Q_6(x)\) denote the number of squarefree integers \(\le x\) that are coprime to \(6\). Every squarefree integer can be written uniquely as \(e m\) with \(e \in \{1,2,3,6\}\) and \((m,6)=1\), so

$$Q(x)=Q_6(x)+Q_6\left(\left\lfloor\frac{x}{2}\right\rfloor\right)+Q_6\left(\left\lfloor\frac{x}{3}\right\rfloor\right)+Q_6\left(\left\lfloor\frac{x}{6}\right\rfloor\right).$$

Rearranging gives the recursion used in the code:

$$Q_6(x)=Q(x)-Q_6\left(\left\lfloor\frac{x}{2}\right\rfloor\right)-Q_6\left(\left\lfloor\frac{x}{3}\right\rfloor\right)-Q_6\left(\left\lfloor\frac{x}{6}\right\rfloor\right).$$

Once \(Q_6\) is memoized, \(B(x)\) follows from the floor-sum identity above.

Step 3: Local Factors at the Primes \(2\) and \(3\)

After the squarefree part away from \(2\) and \(3\) has been isolated, the remaining arithmetic is encoded in four summatory functions. If \(B\) is the common base term from the previous step, the local combinations used by the implementations are

$$P_2(x)=2B(x)+2B\left(\left\lfloor\frac{x}{3}\right\rfloor\right),$$

$$P_4(x)=3B(x)+3B\left(\left\lfloor\frac{x}{3}\right\rfloor\right)-B\left(\left\lfloor\frac{x}{2}\right\rfloor\right)-B\left(\left\lfloor\frac{x}{6}\right\rfloor\right),$$

$$P_{12}(x)=6B(x)-2B\left(\left\lfloor\frac{x}{2}\right\rfloor\right),$$

$$P_{36}(x)=9B(x)-3B\left(\left\lfloor\frac{x}{2}\right\rfloor\right)-3B\left(\left\lfloor\frac{x}{3}\right\rfloor\right)+B\left(\left\lfloor\frac{x}{6}\right\rfloor\right).$$

These are the exact linear combinations implemented in C++, Python, and Java. They represent the different local weights coming from the three Diophantine families once the odd squarefree part has been factored out.

Step 4: Turn the Outer Sums into Harmonic Floor Sums

For each \(D \in \{2,4,12,36\}\), write

$$P_D(x)=\sum_{m \le x} p_D(m).$$

The corresponding outer contribution is then

$$M_D(n)=\sum_{m \le n} p_D(m)\left\lfloor\frac{n}{m}\right\rfloor.$$

The program does not evaluate this sum term by term. Instead, it groups intervals on which \(\left\lfloor n/m \right\rfloor\) is constant. If \(\left\lfloor n/m \right\rfloor=v\) for all \(m \in [L,R]\), then that whole block contributes

$$v\left(P_D(R)-P_D(L-1)\right).$$

Because the quotient \(\left\lfloor n/m \right\rfloor\) takes only \(O(\sqrt{n})\) distinct values, this reduces a linear scan to a near square-root computation.

Step 5: The Modulo-\(3\) Character Correction

The case \(3 \nmid b\) needs one extra correction term. Introduce the nontrivial Dirichlet character modulo \(3\):

$$\chi(m)=\begin{cases} 1,&m\equiv 1\pmod 3,\\ -1,&m\equiv 2\pmod 3,\\ 0,&3\mid m. \end{cases}$$

For squarefree numbers, the weighted prefix sum

$$R(x)=\sum_{\substack{m \le x\\ \mu^2(m)=1}}\chi(m)$$

is again computed by Möbius inversion:

$$R(x)=\sum_{\substack{k \le \sqrt{x}\\ 3 \nmid k}} \mu(k)\sum_{t \le x/k^2}\chi(t).$$

The inner character sum is especially simple, because

$$\sum_{t \le y}\chi(t)=\begin{cases} 1,&y\equiv 1\pmod 3,\\ 0,&y\equiv 0,2\pmod 3. \end{cases}$$

Now define the divisor-weighted function

$$c(m)=\sum_{d \mid m}\mu^2(d)\chi(d).$$

Since \(\chi(d)=0\) whenever \(3 \mid d\), we have \(c(3m)=c(m)\). Therefore the summatory function restricted to \(3 \nmid m\) is obtained by subtraction, and the final correction term has the form

$$H(n)=\sum_{\substack{m \le n\\ \left\lfloor n/m \right\rfloor \equiv 1 \pmod 3}} c(m)\,\mathbf{1}_{3 \nmid m}.$$

This is exactly the extra term that appears in the \(3 \nmid b\) branch of the implementation.

Step 6: Final Assembly by the \(3\)-Adic Valuation of \(b\)

The total is split according to the power of \(3\) dividing \(b\). For the case \(3 \nmid b\), the contribution is

$$C_0(n)=\frac{M_4(n)-M_{36}\left(\left\lfloor\frac{n}{3}\right\rfloor\right)-H(n)}{2}.$$

If \(v_3(b)=t \ge 1\), the local multiplicity becomes \(2t-1\), so the higher \(3\)-adic layers contribute

$$C_{\ge 1}(n)=\sum_{t \ge 1}(2t-1)\left(M_4\left(\left\lfloor\frac{n}{3^t}\right\rfloor\right)-M_{36}\left(\left\lfloor\frac{n}{3^{t+1}}\right\rfloor\right)\right).$$

Putting everything together gives the exact formula implemented by the program:

$$\boxed{T(n)=M_2(n)+M_{12}(n)+C_0(n)+C_{\ge 1}(n).}$$

How the Code Works

The C++, Python, and Java implementations follow the same pipeline. They precompute Möbius values up to \(\lfloor\sqrt{n}\rfloor\), memoize the count of squarefree integers coprime to \(6\), build the four local summatory functions above, evaluate each outer divisor sum through floor-division blocks, and then add the character correction together with the \(3\)-adic layers. As a sanity check, the arithmetic is verified on small checkpoints such as \(T(1)=9\), \(T(20)=732\), and \(T(3000)=438106\) before the full target is evaluated.

Complexity Analysis

The Möbius sieve up to \(\lfloor\sqrt{n}\rfloor\) costs \(O(\sqrt{n})\) time and memory. Every cached summatory function is queried only at values of the form \(\left\lfloor n/k \right\rfloor\), and there are only \(O(\sqrt{n})\) distinct quotients. In practice the running time is dominated by a small number of harmonic block sums plus memoized recursive lookups, which is fast enough for \(n=10^9\).

Footnotes and References

  1. Problem page: https://projecteuler.net/problem=428
  2. Möbius function: Wikipedia - Möbius function
  3. Squarefree integer: Wikipedia - Squarefree integer
  4. Dirichlet character: Wikipedia - Dirichlet character
  5. Dirichlet hyperbola method: Wikipedia - Dirichlet hyperbola method

Problemzusammenfassung

Für jede positive ganze Zahl \(n\) sei \(T(n)\) die Anzahl positiver Tripel \((a,b,c)\) mit \(b \le n\), die die Necklace-Bedingung aus der Aufgabenstellung erfüllen. Eine direkte geometrische Suche ist für \(n=10^9\) völlig unrealistisch. Daher übersetzen die Implementierungen die Geometrie zuerst in Divisor-Gleichungen und berechnen diese anschließend mit zwischengespeicherten Summenfunktionen.

Mathematischer Ansatz

Alle drei Implementierungen beginnen mit derselben arithmetischen Reduktion der Necklace-Bedingung. Die relevanten Kettenlängen sind \(k=3\), \(k=4\) und \(k=6\), woraus die diophantischen Familien

$$k=3:\ (a-3b)(c-3b)=12b^2,$$

$$k=4:\ (a-b)(c-b)=2b^2,$$

$$k=6:\ (3a-b)(3c-b)=4b^2$$

folgen. Damit wird das geometrische Problem zu einer Zählung zulässiger Faktorisierungen von \(2b^2\), \(4b^2\), \(12b^2\) und \(36b^2\), ergänzt um lokale Kongruenzbedingungen an den Primzahlen \(2\) und \(3\).

Schritt 1: Den quadratfreien Anteil außerhalb von \(2\) und \(3\) isolieren

Für Primzahlen \(p \ge 5\) ist nur das Muster der quadratfreien Teiler relevant. Definiere

$$u(m)=\sum_{\substack{d \mid m\\ \mu^2(d)=1\\ (d,6)=1}} 1,$$

also die Anzahl quadratfreier Teiler von \(m\), die teilerfremd zu \(6\) sind. Die zugehörige Summenfunktion lautet

$$B(x)=\sum_{m \le x} u(m)=\sum_{\substack{d \le x\\ \mu^2(d)=1\\ (d,6)=1}}\left\lfloor\frac{x}{d}\right\rfloor.$$

Diese Funktion ist der Grundbaustein der gesamten Rechnung. Anschaulich bündelt sie alle Beiträge der Primzahlen außer \(2\) und \(3\), während die Sonderbehandlung dieser beiden Primzahlen später durch lokale Korrekturen erfolgt.

Schritt 2: Quadratfreie Zahlen zählen, die zu \(6\) teilerfremd sind

Zur Berechnung von \(B(x)\) werden zunächst quadratfreie Zahlen per Möbius-Inversion gezählt:

$$Q(x)=\sum_{k \le \sqrt{x}} \mu(k)\left\lfloor\frac{x}{k^2}\right\rfloor.$$

Sei nun \(Q_6(x)\) die Anzahl quadratfreier Zahlen \(\le x\), die zu \(6\) teilerfremd sind. Jede quadratfreie Zahl lässt sich eindeutig als \(e m\) mit \(e \in \{1,2,3,6\}\) und \((m,6)=1\) schreiben. Daher gilt

$$Q(x)=Q_6(x)+Q_6\left(\left\lfloor\frac{x}{2}\right\rfloor\right)+Q_6\left(\left\lfloor\frac{x}{3}\right\rfloor\right)+Q_6\left(\left\lfloor\frac{x}{6}\right\rfloor\right).$$

Umgestellt ergibt das genau die Rekursion aus dem Programm:

$$Q_6(x)=Q(x)-Q_6\left(\left\lfloor\frac{x}{2}\right\rfloor\right)-Q_6\left(\left\lfloor\frac{x}{3}\right\rfloor\right)-Q_6\left(\left\lfloor\frac{x}{6}\right\rfloor\right).$$

Sobald \(Q_6\) memoisiert ist, folgt \(B(x)\) direkt aus der obigen Floor-Summen-Formel.

Schritt 3: Lokale Faktoren an den Primzahlen \(2\) und \(3\)

Nachdem der quadratfreie Anteil außerhalb von \(2\) und \(3\) abgetrennt wurde, bleibt die lokale Arithmetik in vier Summenfunktionen übrig. Mit \(B\) als Basisterm aus dem vorigen Schritt verwendet die Implementierung die Kombinationen

$$P_2(x)=2B(x)+2B\left(\left\lfloor\frac{x}{3}\right\rfloor\right),$$

$$P_4(x)=3B(x)+3B\left(\left\lfloor\frac{x}{3}\right\rfloor\right)-B\left(\left\lfloor\frac{x}{2}\right\rfloor\right)-B\left(\left\lfloor\frac{x}{6}\right\rfloor\right),$$

$$P_{12}(x)=6B(x)-2B\left(\left\lfloor\frac{x}{2}\right\rfloor\right),$$

$$P_{36}(x)=9B(x)-3B\left(\left\lfloor\frac{x}{2}\right\rfloor\right)-3B\left(\left\lfloor\frac{x}{3}\right\rfloor\right)+B\left(\left\lfloor\frac{x}{6}\right\rfloor\right).$$

Genau diese linearen Kombinationen stehen in den C++-, Python- und Java-Implementierungen. Sie kodieren die unterschiedlichen lokalen Gewichte der drei diophantischen Familien, nachdem der ungerade quadratfreie Anteil ausgelagert wurde.

Schritt 4: Äußere Summen als harmonische Floor-Summen schreiben

Für jedes \(D \in \{2,4,12,36\}\) sei

$$P_D(x)=\sum_{m \le x} p_D(m).$$

Der zugehörige äußere Beitrag ist dann

$$M_D(n)=\sum_{m \le n} p_D(m)\left\lfloor\frac{n}{m}\right\rfloor.$$

Das Programm summiert nicht termweise, sondern blockweise über Intervalle mit konstantem Quotienten \(\left\lfloor n/m \right\rfloor\). Gilt \(\left\lfloor n/m \right\rfloor=v\) für alle \(m \in [L,R]\), so trägt der gesamte Block

$$v\left(P_D(R)-P_D(L-1)\right)$$

bei. Da der Quotient \(\left\lfloor n/m \right\rfloor\) nur \(O(\sqrt{n})\) verschiedene Werte annimmt, schrumpft eine lineare Summation auf eine Rechnung von Größenordnung \(\sqrt{n}\).

Schritt 5: Korrektur mit dem Charakter modulo \(3\)

Der Fall \(3 \nmid b\) benötigt einen zusätzlichen Korrekturterm. Dazu wird der nichttriviale Dirichlet-Charakter modulo \(3\) eingeführt:

$$\chi(m)=\begin{cases} 1,&m\equiv 1\pmod 3,\\ -1,&m\equiv 2\pmod 3,\\ 0,&3\mid m. \end{cases}$$

Für quadratfreie Zahlen betrachtet man die gewichtete Präfixsumme

$$R(x)=\sum_{\substack{m \le x\\ \mu^2(m)=1}}\chi(m),$$

die erneut per Möbius-Inversion berechnet wird:

$$R(x)=\sum_{\substack{k \le \sqrt{x}\\ 3 \nmid k}} \mu(k)\sum_{t \le x/k^2}\chi(t).$$

Die innere Charaktersumme ist besonders einfach, denn

$$\sum_{t \le y}\chi(t)=\begin{cases} 1,&y\equiv 1\pmod 3,\\ 0,&y\equiv 0,2\pmod 3. \end{cases}$$

Außerdem definiert man

$$c(m)=\sum_{d \mid m}\mu^2(d)\chi(d).$$

Weil \(\chi(d)=0\) gilt, sobald \(3 \mid d\), folgt \(c(3m)=c(m)\). Damit erhält man die auf \(3 \nmid m\) eingeschränkte Summenfunktion durch Subtraktion, und der Korrekturterm lautet

$$H(n)=\sum_{\substack{m \le n\\ \left\lfloor n/m \right\rfloor \equiv 1 \pmod 3}} c(m)\,\mathbf{1}_{3 \nmid m}.$$

Genau dieser Zusatz erscheint im Programm im Zweig für \(3 \nmid b\).

Schritt 6: Endzusammensetzung nach der \(3\)-adischen Bewertung von \(b\)

Die Gesamtzahl wird nach der Potenz von \(3\) zerlegt, die \(b\) teilt. Für den Fall \(3 \nmid b\) ist der Beitrag

$$C_0(n)=\frac{M_4(n)-M_{36}\left(\left\lfloor\frac{n}{3}\right\rfloor\right)-H(n)}{2}.$$

Falls \(v_3(b)=t \ge 1\), lautet die lokale Vielfachheit \(2t-1\). Damit tragen die höheren \(3\)-adischen Schichten

$$C_{\ge 1}(n)=\sum_{t \ge 1}(2t-1)\left(M_4\left(\left\lfloor\frac{n}{3^t}\right\rfloor\right)-M_{36}\left(\left\lfloor\frac{n}{3^{t+1}}\right\rfloor\right)\right)$$

bei. Insgesamt ergibt sich damit exakt die im Programm benutzte Formel

$$\boxed{T(n)=M_2(n)+M_{12}(n)+C_0(n)+C_{\ge 1}(n).}$$

Wie der Code arbeitet

Die C++-, Python- und Java-Implementierungen verwenden dieselbe Pipeline: Zuerst werden die Möbius-Werte bis \(\lfloor\sqrt{n}\rfloor\) vorab berechnet. Dann werden die Anzahlen quadratfreier Zahlen teilerfremd zu \(6\) memoisiert, daraus die vier lokalen Summenfunktionen aufgebaut und alle äußeren Divisorsummen per Blockzerlegung nach konstantem Floor-Quotienten ausgewertet. Anschließend kommen die Charakterkorrektur und die \(3\)-adischen Schichten hinzu. Zur Kontrolle werden kleine Prüfpunkte wie \(T(1)=9\), \(T(20)=732\) und \(T(3000)=438106\) berechnet, bevor das eigentliche Ziel ausgewertet wird.

Komplexitätsanalyse

Das Möbius-Sieb bis \(\lfloor\sqrt{n}\rfloor\) benötigt \(O(\sqrt{n})\) Zeit und Speicher. Jede gecachte Summenfunktion wird nur an Werten der Form \(\left\lfloor n/k \right\rfloor\) abgefragt, und davon gibt es nur \(O(\sqrt{n})\) verschiedene. Praktisch wird die Laufzeit von wenigen harmonischen Blocksummen und memoisierten rekursiven Nachschlägen dominiert; damit ist \(n=10^9\) gut beherrschbar.

Fußnoten und Quellen

  1. Problemseite: https://projecteuler.net/problem=428
  2. Möbius-Funktion: Wikipedia - Möbiusfunktion
  3. Quadratfreie Zahl: Wikipedia - Quadratfreie Zahl
  4. Dirichlet-Charakter: Wikipedia - Dirichlet-Charakter
  5. Dirichlet-Hyperbelmethode: Wikipedia - Dirichlet hyperbola method

Problem Özeti

Her pozitif \(n\) için \(T(n)\), \(b \le n\) koşulunu sağlayan ve problemdeki necklace koşulunu taşıyan pozitif tam sayı üçlülerinin \((a,b,c)\) sayısıdır. \(n=10^9\) için doğrudan geometrik arama mümkün değildir. Bu yüzden çözüm, geometriyi önce bölen sayımına indirger, sonra da bu aritmetik yapıyı önbelleğe alınmış toplam fonksiyonlarıyla hesaplar.

Matematiksel Yaklaşım

Üç uygulamanın da çıkış noktası aynıdır: necklace koşulu aritmetik eşitliklere indirgenir. İlgili zincir uzunlukları \(k=3\), \(k=4\) ve \(k=6\) olup şu Diofant denklemlerini verir:

$$k=3:\ (a-3b)(c-3b)=12b^2,$$

$$k=4:\ (a-b)(c-b)=2b^2,$$

$$k=6:\ (3a-b)(3c-b)=4b^2.$$

Böylece geometrik yapı, \(2b^2\), \(4b^2\), \(12b^2\) ve \(36b^2\) sayılarının uygun çarpan çiftlerini saymaya dönüşür. Fakat bu sayımda özellikle \(2\) ve \(3\) asallarındaki yerel kongruans koşulları ayrıca ele alınmalıdır.

Adım 1: \(2\) ve \(3\) dışındaki karefree kısmı ayır

\(p \ge 5\) olan asallar için yalnızca karefree bölen deseni önemlidir. Şunu tanımlayalım:

$$u(m)=\sum_{\substack{d \mid m\\ \mu^2(d)=1\\ (d,6)=1}} 1,$$

yani \(m\)'in \(6\) ile aralarında asal olan karefree bölen sayısı. Bunun kümülatif toplamı

$$B(x)=\sum_{m \le x} u(m)=\sum_{\substack{d \le x\\ \mu^2(d)=1\\ (d,6)=1}}\left\lfloor\frac{x}{d}\right\rfloor$$

şeklindedir. Tüm hesabın temel taşı bu fonksiyondur. Sezgisel olarak, \(2\) ve \(3\) dışındaki tüm asal katkılar burada toplanır; \(2\) ve \(3\) için özel davranış ise ayrı yerel katsayılarla işlenir.

Adım 2: \(6\) ile aralarında asal karefree sayıları say

\(B(x)\)'i hesaplamak için önce Möbius terslemesi ile karefree sayılar sayılır:

$$Q(x)=\sum_{k \le \sqrt{x}} \mu(k)\left\lfloor\frac{x}{k^2}\right\rfloor.$$

Şimdi \(Q_6(x)\), \(x\)'e kadar olan ve \(6\) ile aralarında asal karefree sayıların sayısı olsun. Her karefree sayı tek biçimde \(e m\) olarak yazılır; burada \(e \in \{1,2,3,6\}\) ve \((m,6)=1\). Dolayısıyla

$$Q(x)=Q_6(x)+Q_6\left(\left\lfloor\frac{x}{2}\right\rfloor\right)+Q_6\left(\left\lfloor\frac{x}{3}\right\rfloor\right)+Q_6\left(\left\lfloor\frac{x}{6}\right\rfloor\right).$$

Buradan kodda kullanılan özyinelemeli bağıntı gelir:

$$Q_6(x)=Q(x)-Q_6\left(\left\lfloor\frac{x}{2}\right\rfloor\right)-Q_6\left(\left\lfloor\frac{x}{3}\right\rfloor\right)-Q_6\left(\left\lfloor\frac{x}{6}\right\rfloor\right).$$

\(Q_6\) önbelleğe alındığında, \(B(x)\) yukarıdaki floor-toplam özdeşliğinden hemen elde edilir.

Adım 3: \(2\) ve \(3\) asallarındaki yerel faktörler

\(2\) ve \(3\) dışındaki karefree kısım ayrıldıktan sonra geriye dört kümülatif yerel fonksiyon kalır. Önceki adımdaki ortak taban \(B\) ile uygulamaların kullandığı doğrusal birleşimler şunlardır:

$$P_2(x)=2B(x)+2B\left(\left\lfloor\frac{x}{3}\right\rfloor\right),$$

$$P_4(x)=3B(x)+3B\left(\left\lfloor\frac{x}{3}\right\rfloor\right)-B\left(\left\lfloor\frac{x}{2}\right\rfloor\right)-B\left(\left\lfloor\frac{x}{6}\right\rfloor\right),$$

$$P_{12}(x)=6B(x)-2B\left(\left\lfloor\frac{x}{2}\right\rfloor\right),$$

$$P_{36}(x)=9B(x)-3B\left(\left\lfloor\frac{x}{2}\right\rfloor\right)-3B\left(\left\lfloor\frac{x}{3}\right\rfloor\right)+B\left(\left\lfloor\frac{x}{6}\right\rfloor\right).$$

Bunlar genel amaçlı formüller değil, bu probleme ait C++, Python ve Java uygulamalarında tam olarak yer alan yerel ağırlıklardır. Üç diofantik ailenin \(2\) ve \(3\) üzerindeki davranışını kodlarlar.

Adım 4: Dış toplamları harmonik floor-toplamlarına çevir

Her \(D \in \{2,4,12,36\}\) için

$$P_D(x)=\sum_{m \le x} p_D(m)$$

olsun. Buna karşılık gelen dış katkı

$$M_D(n)=\sum_{m \le n} p_D(m)\left\lfloor\frac{n}{m}\right\rfloor$$

şeklindedir. Program bu toplamı tek tek terimler üzerinden değil, \(\left\lfloor n/m \right\rfloor\) değerinin sabit kaldığı bloklar üzerinden hesaplar. Eğer tüm \(m \in [L,R]\) için \(\left\lfloor n/m \right\rfloor=v\) ise, o bloğun katkısı

$$v\left(P_D(R)-P_D(L-1)\right)$$

olur. \(\left\lfloor n/m \right\rfloor\) yalnızca \(O(\sqrt{n})\) farklı değer aldığı için, doğrusal maliyet karekök mertebesine iner.

Adım 5: Modulo \(3\) karakter düzeltmesi

\(3 \nmid b\) durumu için ek bir düzeltme gerekir. Bunun için \(3\) modülündeki önemsiz olmayan Dirichlet karakterini tanımlayalım:

$$\chi(m)=\begin{cases} 1,&m\equiv 1\pmod 3,\\ -1,&m\equiv 2\pmod 3,\\ 0,&3\mid m. \end{cases}$$

Karefree sayılar üzerinde ağırlıklı önek toplamı

$$R(x)=\sum_{\substack{m \le x\\ \mu^2(m)=1}}\chi(m)$$

yine Möbius terslemesiyle hesaplanır:

$$R(x)=\sum_{\substack{k \le \sqrt{x}\\ 3 \nmid k}} \mu(k)\sum_{t \le x/k^2}\chi(t).$$

İç karakter toplamı çok basittir, çünkü

$$\sum_{t \le y}\chi(t)=\begin{cases} 1,&y\equiv 1\pmod 3,\\ 0,&y\equiv 0,2\pmod 3. \end{cases}$$

Buna ek olarak

$$c(m)=\sum_{d \mid m}\mu^2(d)\chi(d)$$

tanımı yapılır. \(\chi(d)=0\) olduğu her durumda \(3 \mid d\) geçerlidir; dolayısıyla \(c(3m)=c(m)\) elde edilir. Böylece \(3 \nmid m\) kısıtlı toplama fark alarak ulaşılır ve düzeltme terimi

$$H(n)=\sum_{\substack{m \le n\\ \left\lfloor n/m \right\rfloor \equiv 1 \pmod 3}} c(m)\,\mathbf{1}_{3 \nmid m}$$

biçimini alır. Kodun \(3 \nmid b\) dalındaki ek terim tam olarak budur.

Adım 6: Son birleşim, \(b\)'nin \(3\)-adik değeri üzerinden

Toplam sonuç, \(b\)'yi bölen \(3\) kuvvetine göre ayrılır. \(3 \nmid b\) için katkı

$$C_0(n)=\frac{M_4(n)-M_{36}\left(\left\lfloor\frac{n}{3}\right\rfloor\right)-H(n)}{2}$$

olur. Eğer \(v_3(b)=t \ge 1\) ise yerel katsayı \(2t-1\)'e dönüşür ve üst \(3\)-adik katmanların katkısı

$$C_{\ge 1}(n)=\sum_{t \ge 1}(2t-1)\left(M_4\left(\left\lfloor\frac{n}{3^t}\right\rfloor\right)-M_{36}\left(\left\lfloor\frac{n}{3^{t+1}}\right\rfloor\right)\right)$$

şeklinde yazılır. Sonuçta programın kullandığı tam formül

$$\boxed{T(n)=M_2(n)+M_{12}(n)+C_0(n)+C_{\ge 1}(n).}$$

Kodun Çalışma Mantığı

C++, Python ve Java uygulamaları aynı boru hattını izler: önce Möbius değerleri \(\lfloor\sqrt{n}\rfloor\)'ye kadar hesaplanır; sonra \(6\) ile aralarında asal karefree sayı sayımları önbelleğe alınır; ardından yukarıdaki dört yerel kümülatif fonksiyon kurulur; her dış bölen toplamı floor-bölüm bloklarıyla hesaplanır; en sonda da karakter düzeltmesi ile \(3\)-adik katmanlar eklenir. Tam hedefe geçmeden önce \(T(1)=9\), \(T(20)=732\) ve \(T(3000)=438106\) gibi küçük kontrol değerleri de hesaplanır.

Karmaşıklık Analizi

\(\lfloor\sqrt{n}\rfloor\)'ye kadar Möbius eleği \(O(\sqrt{n})\) zaman ve bellek gerektirir. Önbelleğe alınan her toplam fonksiyonu yalnızca \(\left\lfloor n/k \right\rfloor\) biçimindeki değerlere sorulur; bunların farklı sayısı \(O(\sqrt{n})\)'dir. Pratikte çalışma süresi birkaç harmonik blok toplamı ile memoize özyinelemeli çağrı tarafından belirlenir; bu da \(n=10^9\) için yeterince hızlıdır.

Dipnotlar ve Kaynakça

  1. Problem sayfası: https://projecteuler.net/problem=428
  2. Möbius fonksiyonu: Wikipedia - Möbius fonksiyonu
  3. Karefree sayı: Wikipedia - Squarefree integer
  4. Dirichlet karakteri: Wikipedia - Dirichlet character
  5. Dirichlet hiperbol yöntemi: Wikipedia - Dirichlet hyperbola method

Resumen del Problema

Para cada entero positivo \(n\), \(T(n)\) es el número de ternas positivas \((a,b,c)\) con \(b \le n\) que satisfacen la condición necklace del enunciado. Una búsqueda geométrica directa es imposible para \(n=10^9\). Por eso la solución usada en las implementaciones convierte primero la geometría en un problema de conteo de divisores y luego evalúa esa descripción aritmética mediante funciones sumatorias con memoria.

Enfoque Matemático

Las tres implementaciones parten de la misma reducción aritmética de la condición necklace. Las longitudes de cadena relevantes son \(k=3\), \(k=4\) y \(k=6\), lo que produce las familias diofánticas

$$k=3:\ (a-3b)(c-3b)=12b^2,$$

$$k=4:\ (a-b)(c-b)=2b^2,$$

$$k=6:\ (3a-b)(3c-b)=4b^2.$$

Así, la geometría se transforma en contar factorizaciones admisibles de \(2b^2\), \(4b^2\), \(12b^2\) y \(36b^2\), junto con las restricciones locales de congruencia impuestas por los primos \(2\) y \(3\).

Paso 1: Separar la parte libre de cuadrados fuera de \(2\) y \(3\)

Para los primos \(p \ge 5\), solo importa el patrón de divisores squarefree. Definimos

$$u(m)=\sum_{\substack{d \mid m\\ \mu^2(d)=1\\ (d,6)=1}} 1,$$

que cuenta los divisores squarefree de \(m\) coprimos con \(6\). Su función sumatoria es

$$B(x)=\sum_{m \le x} u(m)=\sum_{\substack{d \le x\\ \mu^2(d)=1\\ (d,6)=1}}\left\lfloor\frac{x}{d}\right\rfloor.$$

Esta es la pieza base de todo el cálculo. En esencia, recoge la contribución de todos los primos distintos de \(2\) y \(3\), mientras que el comportamiento local en esos dos primos se corrige aparte.

Paso 2: Contar enteros squarefree coprimos con \(6\)

Para evaluar \(B(x)\), primero se cuentan los números squarefree mediante inversión de Möbius:

$$Q(x)=\sum_{k \le \sqrt{x}} \mu(k)\left\lfloor\frac{x}{k^2}\right\rfloor.$$

Sea \(Q_6(x)\) el número de enteros squarefree \(\le x\) coprimos con \(6\). Todo entero squarefree se escribe de manera única como \(e m\), con \(e \in \{1,2,3,6\}\) y \((m,6)=1\). Por tanto,

$$Q(x)=Q_6(x)+Q_6\left(\left\lfloor\frac{x}{2}\right\rfloor\right)+Q_6\left(\left\lfloor\frac{x}{3}\right\rfloor\right)+Q_6\left(\left\lfloor\frac{x}{6}\right\rfloor\right).$$

Al reorganizar, aparece exactamente la recurrencia usada por la implementación:

$$Q_6(x)=Q(x)-Q_6\left(\left\lfloor\frac{x}{2}\right\rfloor\right)-Q_6\left(\left\lfloor\frac{x}{3}\right\rfloor\right)-Q_6\left(\left\lfloor\frac{x}{6}\right\rfloor\right).$$

Una vez memorizada \(Q_6\), la función \(B(x)\) sale directamente de la identidad de suma con pisos.

Paso 3: Factores locales en los primos \(2\) y \(3\)

Después de aislar la parte squarefree fuera de \(2\) y \(3\), la aritmética restante queda codificada en cuatro funciones sumatorias. Usando \(B\) como término base, las implementaciones emplean

$$P_2(x)=2B(x)+2B\left(\left\lfloor\frac{x}{3}\right\rfloor\right),$$

$$P_4(x)=3B(x)+3B\left(\left\lfloor\frac{x}{3}\right\rfloor\right)-B\left(\left\lfloor\frac{x}{2}\right\rfloor\right)-B\left(\left\lfloor\frac{x}{6}\right\rfloor\right),$$

$$P_{12}(x)=6B(x)-2B\left(\left\lfloor\frac{x}{2}\right\rfloor\right),$$

$$P_{36}(x)=9B(x)-3B\left(\left\lfloor\frac{x}{2}\right\rfloor\right)-3B\left(\left\lfloor\frac{x}{3}\right\rfloor\right)+B\left(\left\lfloor\frac{x}{6}\right\rfloor\right).$$

Estas combinaciones lineales son exactamente las que aparecen en C++, Python y Java. Representan los pesos locales distintos de las familias diofánticas una vez separado el componente squarefree impar.

Paso 4: Reescribir las sumas exteriores como sumas armónicas con pisos

Para cada \(D \in \{2,4,12,36\}\), escribimos

$$P_D(x)=\sum_{m \le x} p_D(m).$$

La contribución exterior correspondiente es

$$M_D(n)=\sum_{m \le n} p_D(m)\left\lfloor\frac{n}{m}\right\rfloor.$$

El programa no suma término por término. Agrupa intervalos donde \(\left\lfloor n/m \right\rfloor\) permanece constante. Si \(\left\lfloor n/m \right\rfloor=v\) para todo \(m \in [L,R]\), la contribución del bloque es

$$v\left(P_D(R)-P_D(L-1)\right).$$

Como el cociente \(\left\lfloor n/m \right\rfloor\) toma solo \(O(\sqrt{n})\) valores distintos, la suma pasa de lineal a casi \(\sqrt{n}\).

Paso 5: Corrección con el carácter módulo \(3\)

El caso \(3 \nmid b\) necesita una corrección adicional. Se introduce el carácter de Dirichlet no trivial módulo \(3\):

$$\chi(m)=\begin{cases} 1,&m\equiv 1\pmod 3,\\ -1,&m\equiv 2\pmod 3,\\ 0,&3\mid m. \end{cases}$$

Sobre los números squarefree se considera la suma prefijo ponderada

$$R(x)=\sum_{\substack{m \le x\\ \mu^2(m)=1}}\chi(m),$$

que vuelve a calcularse con inversión de Möbius:

$$R(x)=\sum_{\substack{k \le \sqrt{x}\\ 3 \nmid k}} \mu(k)\sum_{t \le x/k^2}\chi(t).$$

La suma interior del carácter es especialmente simple, porque

$$\sum_{t \le y}\chi(t)=\begin{cases} 1,&y\equiv 1\pmod 3,\\ 0,&y\equiv 0,2\pmod 3. \end{cases}$$

Luego se define

$$c(m)=\sum_{d \mid m}\mu^2(d)\chi(d).$$

Como \(\chi(d)=0\) siempre que \(3 \mid d\), se obtiene \(c(3m)=c(m)\). Así, la versión restringida a \(3 \nmid m\) sale por diferencia, y la corrección final es

$$H(n)=\sum_{\substack{m \le n\\ \left\lfloor n/m \right\rfloor \equiv 1 \pmod 3}} c(m)\,\mathbf{1}_{3 \nmid m}.$$

Ese es exactamente el término adicional que usa la implementación en la rama \(3 \nmid b\).

Paso 6: Ensamblaje final según la valuación \(3\)-ádica de \(b\)

El total se separa según la potencia de \(3\) que divide a \(b\). Cuando \(3 \nmid b\), la contribución es

$$C_0(n)=\frac{M_4(n)-M_{36}\left(\left\lfloor\frac{n}{3}\right\rfloor\right)-H(n)}{2}.$$

Si \(v_3(b)=t \ge 1\), la multiplicidad local pasa a ser \(2t-1\), de modo que las capas superiores aportan

$$C_{\ge 1}(n)=\sum_{t \ge 1}(2t-1)\left(M_4\left(\left\lfloor\frac{n}{3^t}\right\rfloor\right)-M_{36}\left(\left\lfloor\frac{n}{3^{t+1}}\right\rfloor\right)\right).$$

Por tanto, la fórmula exacta implementada por el programa es

$$\boxed{T(n)=M_2(n)+M_{12}(n)+C_0(n)+C_{\ge 1}(n).}$$

Cómo Funciona el Código

Las implementaciones en C++, Python y Java siguen la misma secuencia: primero precalculan los valores de Möbius hasta \(\lfloor\sqrt{n}\rfloor\); luego memorizan el conteo de enteros squarefree coprimos con \(6\); construyen las cuatro funciones sumatorias locales; evalúan cada suma exterior agrupando bloques con el mismo cociente entero; y finalmente añaden la corrección con el carácter y las capas \(3\)-ádicas. Como comprobación interna, calculan valores pequeños como \(T(1)=9\), \(T(20)=732\) y \(T(3000)=438106\) antes de atacar el objetivo final.

Complejidad

La criba de Möbius hasta \(\lfloor\sqrt{n}\rfloor\) cuesta \(O(\sqrt{n})\) tiempo y memoria. Cada función sumatoria memorizada solo se consulta en valores de la forma \(\left\lfloor n/k \right\rfloor\), y hay únicamente \(O(\sqrt{n})\) cocientes distintos. En la práctica, el coste está dominado por unas pocas sumas armónicas por bloques y búsquedas recursivas con caché, lo cual hace viable el caso \(n=10^9\).

Notas y Referencias

  1. Página del problema: https://projecteuler.net/problem=428
  2. Función de Möbius: Wikipedia - Función de Möbius
  3. Entero squarefree: Wikipedia - Squarefree integer
  4. Carácter de Dirichlet: Wikipedia - Carácter de Dirichlet
  5. Método de la hipérbola de Dirichlet: Wikipedia - Dirichlet hyperbola method

问题概述

对每个正整数 \(n\),记 \(T(n)\) 为所有满足 necklace 条件且 \(b \le n\) 的正整数三元组 \((a,b,c)\) 的个数。真正的目标是 \(n=10^9\),直接从几何配置暴力搜索完全不可行。所以实现先把几何闭链条件改写成算术分解问题,再用若干带缓存的求和函数来计算。

数学方法

三个实现都从同一个算术化步骤出发。与闭链条件对应的相关链长只有 \(k=3\)、\(k=4\)、\(k=6\),于是得到三族丢番图方程:

$$k=3:\ (a-3b)(c-3b)=12b^2,$$

$$k=4:\ (a-b)(c-b)=2b^2,$$

$$k=6:\ (3a-b)(3c-b)=4b^2.$$

这样一来,几何问题就变成了对 \(2b^2\)、\(4b^2\)、\(12b^2\)、\(36b^2\) 的可行因子分解进行计数,只是其中在素数 \(2\) 和 \(3\) 处还带有额外的局部同余限制。

步骤 1:把 \(2\) 和 \(3\) 之外的平方自由部分单独拿出来

对所有 \(p \ge 5\) 的素数,真正重要的是平方自由因子的选择方式。定义

$$u(m)=\sum_{\substack{d \mid m\\ \mu^2(d)=1\\ (d,6)=1}} 1,$$

也就是 \(m\) 的、与 \(6\) 互素的平方自由因子个数。它的前缀和写成

$$B(x)=\sum_{m \le x} u(m)=\sum_{\substack{d \le x\\ \mu^2(d)=1\\ (d,6)=1}}\left\lfloor\frac{x}{d}\right\rfloor.$$

这就是整个程序的基础模块。直观上,它统一处理了除 \(2\)、\(3\) 之外所有素数的贡献,而 \(2\)、\(3\) 的特殊局部行为则在后面单独补进去。

步骤 2:统计与 \(6\) 互素的平方自由整数

为了计算 \(B(x)\),实现先用 Möbius 反演计算平方自由整数个数:

$$Q(x)=\sum_{k \le \sqrt{x}} \mu(k)\left\lfloor\frac{x}{k^2}\right\rfloor.$$

再设 \(Q_6(x)\) 表示不超过 \(x\) 且与 \(6\) 互素的平方自由整数个数。每个平方自由整数都能唯一写成 \(e m\),其中 \(e \in \{1,2,3,6\}\),并且 \((m,6)=1\)。因此

$$Q(x)=Q_6(x)+Q_6\left(\left\lfloor\frac{x}{2}\right\rfloor\right)+Q_6\left(\left\lfloor\frac{x}{3}\right\rfloor\right)+Q_6\left(\left\lfloor\frac{x}{6}\right\rfloor\right).$$

移项后得到程序里真正使用的递推式:

$$Q_6(x)=Q(x)-Q_6\left(\left\lfloor\frac{x}{2}\right\rfloor\right)-Q_6\left(\left\lfloor\frac{x}{3}\right\rfloor\right)-Q_6\left(\left\lfloor\frac{x}{6}\right\rfloor\right).$$

一旦 \(Q_6\) 被缓存起来,\(B(x)\) 就能通过上面的整除分块恒等式快速得到。

步骤 3:素数 \(2\) 和 \(3\) 处的局部因子

把 \(2\)、\(3\) 之外的平方自由部分剥离之后,剩下的局部算术被编码成四个前缀函数。以前一步的 \(B\) 为共同基底,程序实际使用的是

$$P_2(x)=2B(x)+2B\left(\left\lfloor\frac{x}{3}\right\rfloor\right),$$

$$P_4(x)=3B(x)+3B\left(\left\lfloor\frac{x}{3}\right\rfloor\right)-B\left(\left\lfloor\frac{x}{2}\right\rfloor\right)-B\left(\left\lfloor\frac{x}{6}\right\rfloor\right),$$

$$P_{12}(x)=6B(x)-2B\left(\left\lfloor\frac{x}{2}\right\rfloor\right),$$

$$P_{36}(x)=9B(x)-3B\left(\left\lfloor\frac{x}{2}\right\rfloor\right)-3B\left(\left\lfloor\frac{x}{3}\right\rfloor\right)+B\left(\left\lfloor\frac{x}{6}\right\rfloor\right).$$

这些并不是一般性的数论公式,而是这个问题的 C++、Python、Java 实现中精确编码的局部权重组合。它们反映了三族丢番图方程在 \(2\) 和 \(3\) 处的不同局部行为。

步骤 4:把外层求和改写成调和分块的 floor-sum

对每个 \(D \in \{2,4,12,36\}\),记

$$P_D(x)=\sum_{m \le x} p_D(m).$$

对应的外层贡献就是

$$M_D(n)=\sum_{m \le n} p_D(m)\left\lfloor\frac{n}{m}\right\rfloor.$$

程序不会逐项累加,而是把所有使 \(\left\lfloor n/m \right\rfloor\) 相同的区间放到同一个块里。如果对区间 \(m \in [L,R]\) 都有 \(\left\lfloor n/m \right\rfloor=v\),那么这一整块的贡献为

$$v\left(P_D(R)-P_D(L-1)\right).$$

由于商 \(\left\lfloor n/m \right\rfloor\) 只会取到 \(O(\sqrt{n})\) 个不同值,因此时间复杂度从线性降到了接近 \(\sqrt{n}\) 的量级。

步骤 5:模 \(3\) 角色修正项

当 \(3 \nmid b\) 时,还需要额外的修正。这里引入模 \(3\) 的非平凡 Dirichlet 角色:

$$\chi(m)=\begin{cases} 1,&m\equiv 1\pmod 3,\\ -1,&m\equiv 2\pmod 3,\\ 0,&3\mid m. \end{cases}$$

在平方自由整数上考虑带权前缀和

$$R(x)=\sum_{\substack{m \le x\\ \mu^2(m)=1}}\chi(m),$$

它同样由 Möbius 反演给出:

$$R(x)=\sum_{\substack{k \le \sqrt{x}\\ 3 \nmid k}} \mu(k)\sum_{t \le x/k^2}\chi(t).$$

内层角色和特别简单,因为

$$\sum_{t \le y}\chi(t)=\begin{cases} 1,&y\equiv 1\pmod 3,\\ 0,&y\equiv 0,2\pmod 3. \end{cases}$$

再定义

$$c(m)=\sum_{d \mid m}\mu^2(d)\chi(d).$$

由于 \(3 \mid d\) 时 \(\chi(d)=0\),于是有 \(c(3m)=c(m)\)。所以限制在 \(3 \nmid m\) 的版本可以通过前缀差得到,最终修正项写成

$$H(n)=\sum_{\substack{m \le n\\ \left\lfloor n/m \right\rfloor \equiv 1 \pmod 3}} c(m)\,\mathbf{1}_{3 \nmid m}.$$

这正是实现里在 \(3 \nmid b\) 分支中出现的附加项。

步骤 6:按 \(b\) 的 \(3\)-进赋值做最终合并

总计数按 \(b\) 被 \(3\) 除的次数来拆分。当 \(3 \nmid b\) 时,贡献为

$$C_0(n)=\frac{M_4(n)-M_{36}\left(\left\lfloor\frac{n}{3}\right\rfloor\right)-H(n)}{2}.$$

如果 \(v_3(b)=t \ge 1\),局部重数变成 \(2t-1\),所以更高的 \(3\)-进层贡献为

$$C_{\ge 1}(n)=\sum_{t \ge 1}(2t-1)\left(M_4\left(\left\lfloor\frac{n}{3^t}\right\rfloor\right)-M_{36}\left(\left\lfloor\frac{n}{3^{t+1}}\right\rfloor\right)\right).$$

于是程序真正实现的总公式就是

$$\boxed{T(n)=M_2(n)+M_{12}(n)+C_0(n)+C_{\ge 1}(n).}$$

代码如何对应这个公式

C++、Python 和 Java 三个实现遵循完全相同的流程:先预处理到 \(\lfloor\sqrt{n}\rfloor\) 的 Möbius 值;再缓存与 \(6\) 互素的平方自由整数计数;接着构造上面的四个局部前缀函数;然后用整除分块计算所有外层求和;最后加上模 \(3\) 的角色修正和各个 \(3\)-进层。为了确认推导没有偏差,程序还会先验证一些小检查点,例如 \(T(1)=9\)、\(T(20)=732\)、\(T(3000)=438106\)。

复杂度分析

预处理到 \(\lfloor\sqrt{n}\rfloor\) 的 Möbius 筛需要 \(O(\sqrt{n})\) 时间和 \(O(\sqrt{n})\) 内存。所有带缓存的前缀函数只会在形如 \(\left\lfloor n/k \right\rfloor\) 的点上被查询,而这样的不同取值只有 \(O(\sqrt{n})\) 个。实际运行时间主要耗在少量调和分块求和和若干递归缓存查询上,因此足以处理 \(n=10^9\)。

参考资料

  1. 题目页面:https://projecteuler.net/problem=428
  2. Möbius 函数:Wikipedia - 莫比乌斯函数
  3. 平方自由整数:Wikipedia - Squarefree integer
  4. Dirichlet 角色:Wikipedia - Dirichlet character
  5. Dirichlet 双曲线方法:Wikipedia - Dirichlet hyperbola method

Краткое описание задачи

Для каждого положительного \(n\) обозначим через \(T(n)\) число положительных троек \((a,b,c)\) с условием \(b \le n\), удовлетворяющих necklace-условию из задачи. Прямой геометрический перебор для \(n=10^9\) невозможен. Поэтому реализация сначала переводит геометрию в задачу о подсчете делителей, а затем вычисляет эту арифметическую формулу с помощью кэшируемых сумматорных функций.

Математический подход

Во всех трех реализациях используется одно и то же арифметическое преобразование necklace-условия. Существенными оказываются длины цепочки \(k=3\), \(k=4\) и \(k=6\), что дает семейства диофантовых уравнений

$$k=3:\ (a-3b)(c-3b)=12b^2,$$

$$k=4:\ (a-b)(c-b)=2b^2,$$

$$k=6:\ (3a-b)(3c-b)=4b^2.$$

Тем самым геометрическая конфигурация сводится к подсчету допустимых факторизаций чисел \(2b^2\), \(4b^2\), \(12b^2\) и \(36b^2\) с дополнительными локальными условиями при простых \(2\) и \(3\).

Шаг 1: Выделить квадратсвободную часть вне \(2\) и \(3\)

Для простых \(p \ge 5\) важна только схема квадратсвободных делителей. Введем

$$u(m)=\sum_{\substack{d \mid m\\ \mu^2(d)=1\\ (d,6)=1}} 1,$$

то есть число квадратсвободных делителей числа \(m\), взаимно простых с \(6\). Его сумматорная функция равна

$$B(x)=\sum_{m \le x} u(m)=\sum_{\substack{d \le x\\ \mu^2(d)=1\\ (d,6)=1}}\left\lfloor\frac{x}{d}\right\rfloor.$$

Это основной строительный блок всей схемы. По сути, он собирает вклад всех простых, кроме \(2\) и \(3\), а локальное поведение в этих двух точках учитывается отдельно.

Шаг 2: Подсчитать квадратсвободные числа, взаимно простые с \(6\)

Чтобы вычислить \(B(x)\), сначала считается число квадратсвободных чисел с помощью обращения Мёбиуса:

$$Q(x)=\sum_{k \le \sqrt{x}} \mu(k)\left\lfloor\frac{x}{k^2}\right\rfloor.$$

Пусть \(Q_6(x)\) обозначает количество квадратсвободных чисел \(\le x\), взаимно простых с \(6\). Каждое квадратсвободное число единственным образом представляется как \(e m\), где \(e \in \{1,2,3,6\}\) и \((m,6)=1\). Поэтому

$$Q(x)=Q_6(x)+Q_6\left(\left\lfloor\frac{x}{2}\right\rfloor\right)+Q_6\left(\left\lfloor\frac{x}{3}\right\rfloor\right)+Q_6\left(\left\lfloor\frac{x}{6}\right\rfloor\right).$$

Перенос членов дает рекурсию, которая и используется в программе:

$$Q_6(x)=Q(x)-Q_6\left(\left\lfloor\frac{x}{2}\right\rfloor\right)-Q_6\left(\left\lfloor\frac{x}{3}\right\rfloor\right)-Q_6\left(\left\lfloor\frac{x}{6}\right\rfloor\right).$$

После мемоизации \(Q_6\) функция \(B(x)\) сразу получается из тождественной суммы с полом.

Шаг 3: Локальные множители при простых \(2\) и \(3\)

Когда квадратсвободная часть вне \(2\) и \(3\) отделена, оставшаяся локальная арифметика кодируется четырьмя сумматорными функциями. Если \(B\) обозначает базовый член из предыдущего шага, то реализации используют комбинации

$$P_2(x)=2B(x)+2B\left(\left\lfloor\frac{x}{3}\right\rfloor\right),$$

$$P_4(x)=3B(x)+3B\left(\left\lfloor\frac{x}{3}\right\rfloor\right)-B\left(\left\lfloor\frac{x}{2}\right\rfloor\right)-B\left(\left\lfloor\frac{x}{6}\right\rfloor\right),$$

$$P_{12}(x)=6B(x)-2B\left(\left\lfloor\frac{x}{2}\right\rfloor\right),$$

$$P_{36}(x)=9B(x)-3B\left(\left\lfloor\frac{x}{2}\right\rfloor\right)-3B\left(\left\lfloor\frac{x}{3}\right\rfloor\right)+B\left(\left\lfloor\frac{x}{6}\right\rfloor\right).$$

Это не абстрактные формулы из учебника, а точные линейные комбинации, зашитые в реализации на C++, Python и Java. Они отражают разные локальные веса, возникающие из трех диофантовых семейств после отделения нечетной квадратсвободной части.

Шаг 4: Преобразовать внешние суммы в гармонические floor-суммы

Для каждого \(D \in \{2,4,12,36\}\) положим

$$P_D(x)=\sum_{m \le x} p_D(m).$$

Тогда соответствующий внешний вклад равен

$$M_D(n)=\sum_{m \le n} p_D(m)\left\lfloor\frac{n}{m}\right\rfloor.$$

Программа не считает эту сумму по одному слагаемому. Она объединяет интервалы, на которых \(\left\lfloor n/m \right\rfloor\) постоянно. Если \(\left\lfloor n/m \right\rfloor=v\) для всех \(m \in [L,R]\), то вклад целого блока равен

$$v\left(P_D(R)-P_D(L-1)\right).$$

Поскольку частное \(\left\lfloor n/m \right\rfloor\) принимает лишь \(O(\sqrt{n})\) различных значений, линейный перебор заменяется вычислением порядка \(\sqrt{n}\).

Шаг 5: Поправка с характером по модулю \(3\)

Случай \(3 \nmid b\) требует дополнительной поправки. Вводится нетривиальный характер Дирихле по модулю \(3\):

$$\chi(m)=\begin{cases} 1,&m\equiv 1\pmod 3,\\ -1,&m\equiv 2\pmod 3,\\ 0,&3\mid m. \end{cases}$$

На квадратсвободных числах рассматривается взвешенная префиксная сумма

$$R(x)=\sum_{\substack{m \le x\\ \mu^2(m)=1}}\chi(m),$$

которая снова вычисляется через обращение Мёбиуса:

$$R(x)=\sum_{\substack{k \le \sqrt{x}\\ 3 \nmid k}} \mu(k)\sum_{t \le x/k^2}\chi(t).$$

Внутренняя сумма характера особенно проста, потому что

$$\sum_{t \le y}\chi(t)=\begin{cases} 1,&y\equiv 1\pmod 3,\\ 0,&y\equiv 0,2\pmod 3. \end{cases}$$

Далее определяется функция

$$c(m)=\sum_{d \mid m}\mu^2(d)\chi(d).$$

Так как \(\chi(d)=0\) при \(3 \mid d\), выполняется \(c(3m)=c(m)\). Поэтому ограничение на \(3 \nmid m\) получается вычитанием, а итоговая поправка имеет вид

$$H(n)=\sum_{\substack{m \le n\\ \left\lfloor n/m \right\rfloor \equiv 1 \pmod 3}} c(m)\,\mathbf{1}_{3 \nmid m}.$$

Именно этот дополнительный член появляется в ветви \(3 \nmid b\) у реализации.

Шаг 6: Итоговая сборка по \(3\)-адической норме числа \(b\)

Общий ответ разбивается по степени \(3\), делящей \(b\). В случае \(3 \nmid b\) вклад равен

$$C_0(n)=\frac{M_4(n)-M_{36}\left(\left\lfloor\frac{n}{3}\right\rfloor\right)-H(n)}{2}.$$

Если \(v_3(b)=t \ge 1\), локальная кратность становится равной \(2t-1\), и верхние \(3\)-адические слои дают

$$C_{\ge 1}(n)=\sum_{t \ge 1}(2t-1)\left(M_4\left(\left\lfloor\frac{n}{3^t}\right\rfloor\right)-M_{36}\left(\left\lfloor\frac{n}{3^{t+1}}\right\rfloor\right)\right).$$

Следовательно, точная формула, реализованная в программе, такова:

$$\boxed{T(n)=M_2(n)+M_{12}(n)+C_0(n)+C_{\ge 1}(n).}$$

Как это отражено в коде

Реализации на C++, Python и Java следуют одной и той же схеме: сначала предварительно вычисляются значения функции Мёбиуса до \(\lfloor\sqrt{n}\rfloor\); затем мемоизируется подсчет квадратсвободных чисел, взаимно простых с \(6\); после этого строятся четыре локальные сумматорные функции; все внешние суммы считаются блоками по одинаковому значению целой части; и в конце добавляются поправка с характером и \(3\)-адические слои. Для самопроверки перед основным запуском вычисляются небольшие контрольные значения \(T(1)=9\), \(T(20)=732\) и \(T(3000)=438106\).

Анализ сложности

Решето Мёбиуса до \(\lfloor\sqrt{n}\rfloor\) требует \(O(\sqrt{n})\) времени и памяти. Каждая кэшируемая сумматорная функция вызывается только на значениях вида \(\left\lfloor n/k \right\rfloor\), а таких различных значений лишь \(O(\sqrt{n})\). На практике время работы определяется несколькими гармоническими блочными суммами и мемоизированными рекурсивными вызовами, поэтому случай \(n=10^9\) оказывается вычислимым.

Ссылки и литература

  1. Страница задачи: https://projecteuler.net/problem=428
  2. Функция Мёбиуса: Wikipedia - Функция Мёбиуса
  3. Квадратсвободное число: Wikipedia - Квадратсвободное число
  4. Характер Дирихле: Wikipedia - Характер Дирихле
  5. Метод гиперболы Дирихле: Wikipedia - Dirichlet hyperbola method

ملخص المسألة

لكل عدد صحيح موجب \(n\)، نرمز بـ \(T(n)\) إلى عدد الثلاثيات الموجبة \((a,b,c)\) التي تحقق شرط necklace في نص المسألة مع القيد \(b \le n\). البحث الهندسي المباشر مستحيل عمليًا عندما يكون \(n=10^9\). لذلك تعتمد الحلول المطبقة على تحويل البنية الهندسية أولًا إلى مسألة عدّ قواسم، ثم حساب تلك الصياغة الحسابية بواسطة دوال جمع تراكمي مع تخزين مؤقت.

المنهج الرياضي

تبدأ التطبيقات الثلاثة من الاختزال الحسابي نفسه لشرط necklace. أطوال السلسلة المهمة هي \(k=3\) و\(k=4\) و\(k=6\)، وهذا يعطي عائلات المعادلات الديوفانتية

$$k=3:\ (a-3b)(c-3b)=12b^2,$$

$$k=4:\ (a-b)(c-b)=2b^2,$$

$$k=6:\ (3a-b)(3c-b)=4b^2.$$

وبذلك تتحول البنية الهندسية إلى عدّ التحليلات العاملية المسموح بها للأعداد \(2b^2\) و\(4b^2\) و\(12b^2\) و\(36b^2\)، مع شروط توافق محلية إضافية عند العددين الأوليين \(2\) و\(3\).

الخطوة 1: فصل الجزء الخالي من المربعات خارج \(2\) و\(3\)

بالنسبة إلى الأعداد الأولية \(p \ge 5\)، المهم هو نمط القواسم الخالية من المربعات فقط. نعرّف

$$u(m)=\sum_{\substack{d \mid m\\ \mu^2(d)=1\\ (d,6)=1}} 1,$$

أي عدد القواسم الخالية من المربعات للعدد \(m\) والتي هي أولية نسبيًا مع \(6\). ودالة المجموع التراكمي لها هي

$$B(x)=\sum_{m \le x} u(m)=\sum_{\substack{d \le x\\ \mu^2(d)=1\\ (d,6)=1}}\left\lfloor\frac{x}{d}\right\rfloor.$$

هذه هي اللبنة الأساسية في الحساب كله. فهي تجمع أثر جميع الأعداد الأولية غير \(2\) و\(3\)، بينما يعالج السلوك المحلي عند \(2\) و\(3\) في مرحلة مستقلة لاحقًا.

الخطوة 2: عدّ الأعداد الخالية من المربعات والأولية نسبيًا مع \(6\)

لحساب \(B(x)\)، تبدأ التطبيقات بعدّ الأعداد الخالية من المربعات باستخدام قلب Möbius:

$$Q(x)=\sum_{k \le \sqrt{x}} \mu(k)\left\lfloor\frac{x}{k^2}\right\rfloor.$$

ولتكن \(Q_6(x)\) عدد الأعداد الخالية من المربعات التي لا تتجاوز \(x\) وتكون أولية نسبيًا مع \(6\). كل عدد خالٍ من المربعات يمكن كتابته بشكل وحيد على الصورة \(e m\)، حيث \(e \in \{1,2,3,6\}\) و\((m,6)=1\). ومن ثم

$$Q(x)=Q_6(x)+Q_6\left(\left\lfloor\frac{x}{2}\right\rfloor\right)+Q_6\left(\left\lfloor\frac{x}{3}\right\rfloor\right)+Q_6\left(\left\lfloor\frac{x}{6}\right\rfloor\right).$$

وبإعادة الترتيب نحصل على العلاقة التكرارية المستخدمة في الكود:

$$Q_6(x)=Q(x)-Q_6\left(\left\lfloor\frac{x}{2}\right\rfloor\right)-Q_6\left(\left\lfloor\frac{x}{3}\right\rfloor\right)-Q_6\left(\left\lfloor\frac{x}{6}\right\rfloor\right).$$

وبعد حفظ قيم \(Q_6\) في الذاكرة، تُستخرج \(B(x)\) مباشرة من هوية مجموع القيم الصحيحة السابقة.

الخطوة 3: العوامل المحلية عند \(2\) و\(3\)

بعد فصل الجزء الخالي من المربعات خارج \(2\) و\(3\)، يبقى حساب محلي يُشفَّر في أربع دوال تراكميّة. باستخدام \(B\) بوصفها الحد الأساسي من الخطوة السابقة، فإن التطبيقات تستعمل التراكيب

$$P_2(x)=2B(x)+2B\left(\left\lfloor\frac{x}{3}\right\rfloor\right),$$

$$P_4(x)=3B(x)+3B\left(\left\lfloor\frac{x}{3}\right\rfloor\right)-B\left(\left\lfloor\frac{x}{2}\right\rfloor\right)-B\left(\left\lfloor\frac{x}{6}\right\rfloor\right),$$

$$P_{12}(x)=6B(x)-2B\left(\left\lfloor\frac{x}{2}\right\rfloor\right),$$

$$P_{36}(x)=9B(x)-3B\left(\left\lfloor\frac{x}{2}\right\rfloor\right)-3B\left(\left\lfloor\frac{x}{3}\right\rfloor\right)+B\left(\left\lfloor\frac{x}{6}\right\rfloor\right).$$

هذه ليست صيغًا نظرية عامة، بل هي بالضبط التراكيب الخطية المشفرة في تطبيقات C++ وPython وJava لهذه المسألة. وهي تمثل الأوزان المحلية المختلفة الناتجة من العائلات الديوفانتية الثلاث بعد إخراج الجزء الفردي الخالي من المربعات.

الخطوة 4: تحويل المجاميع الخارجية إلى مجاميع توافقية من نوع floor-sum

لكل \(D \in \{2,4,12,36\}\)، نكتب

$$P_D(x)=\sum_{m \le x} p_D(m).$$

وعندئذ يكون الإسهام الخارجي الموافق

$$M_D(n)=\sum_{m \le n} p_D(m)\left\lfloor\frac{n}{m}\right\rfloor.$$

لا يحسب البرنامج هذا المجموع حدًا حدًا، بل يجمع الفترات التي تكون فيها القيمة \(\left\lfloor n/m \right\rfloor\) ثابتة. فإذا كانت \(\left\lfloor n/m \right\rfloor=v\) لكل \(m \in [L,R]\)، فإن إسهام الكتلة كلها هو

$$v\left(P_D(R)-P_D(L-1)\right).$$

وبما أن خارج القسمة \(\left\lfloor n/m \right\rfloor\) لا يأخذ إلا \(O(\sqrt{n})\) من القيم المختلفة، فإن الكلفة تنخفض من رتبة خطية إلى رتبة قريبة من \(\sqrt{n}\).

الخطوة 5: حد التصحيح باستعمال المحرف modulo \(3\)

الحالة \(3 \nmid b\) تحتاج إلى حد تصحيح إضافي. نعرّف المحرف الديريشليه غير التافه modulo \(3\):

$$\chi(m)=\begin{cases} 1,&m\equiv 1\pmod 3,\\ -1,&m\equiv 2\pmod 3,\\ 0,&3\mid m. \end{cases}$$

وعلى الأعداد الخالية من المربعات نأخذ مجموعًا تراكميًا موزونًا هو

$$R(x)=\sum_{\substack{m \le x\\ \mu^2(m)=1}}\chi(m),$$

ويُحسب مرة أخرى بقلب Möbius:

$$R(x)=\sum_{\substack{k \le \sqrt{x}\\ 3 \nmid k}} \mu(k)\sum_{t \le x/k^2}\chi(t).$$

أما مجموع المحرف الداخلي فهو بسيط جدًا، لأن

$$\sum_{t \le y}\chi(t)=\begin{cases} 1,&y\equiv 1\pmod 3,\\ 0,&y\equiv 0,2\pmod 3. \end{cases}$$

ثم نعرّف

$$c(m)=\sum_{d \mid m}\mu^2(d)\chi(d).$$

وبما أن \(\chi(d)=0\) كلما كان \(3 \mid d\)، نحصل على \(c(3m)=c(m)\). لذلك يمكن استخراج النسخة المقيدة بشرط \(3 \nmid m\) بالطرح، ويصبح حد التصحيح النهائي

$$H(n)=\sum_{\substack{m \le n\\ \left\lfloor n/m \right\rfloor \equiv 1 \pmod 3}} c(m)\,\mathbf{1}_{3 \nmid m}.$$

وهذا هو الحد الإضافي نفسه الذي يظهر في فرع \(3 \nmid b\) داخل التنفيذ.

الخطوة 6: التجميع النهائي حسب التقييم \(3\)-الأدي للعدد \(b\)

يقسم المجموع الكلي بحسب قوة \(3\) التي تقسم \(b\). في حالة \(3 \nmid b\) يكون الإسهام

$$C_0(n)=\frac{M_4(n)-M_{36}\left(\left\lfloor\frac{n}{3}\right\rfloor\right)-H(n)}{2}.$$

أما إذا كان \(v_3(b)=t \ge 1\)، فإن المضاعف المحلي يصبح \(2t-1\)، ولذلك تسهم الطبقات \(3\)-الأدية الأعلى بالمقدار

$$C_{\ge 1}(n)=\sum_{t \ge 1}(2t-1)\left(M_4\left(\left\lfloor\frac{n}{3^t}\right\rfloor\right)-M_{36}\left(\left\lfloor\frac{n}{3^{t+1}}\right\rfloor\right)\right).$$

ومن ثم فإن الصيغة الدقيقة التي تطبقها البرامج هي

$$\boxed{T(n)=M_2(n)+M_{12}(n)+C_0(n)+C_{\ge 1}(n).}$$

كيف يظهر ذلك في الكود

تسير تطبيقات C++ وPython وJava على المسار نفسه تمامًا: تُحسب أولًا قيم Möbius حتى \(\lfloor\sqrt{n}\rfloor\)، ثم تُخزَّن أعداد الأعداد الخالية من المربعات والأولية نسبيًا مع \(6\)، وبعدها تُبنى الدوال التراكمية المحلية الأربع، وتُحسب جميع المجاميع الخارجية بتجميع فترات القسمة الصحيحة الثابتة، وفي النهاية يُضاف حد المحرف مع الطبقات \(3\)-الأدية. وكاختبار اتساق، تتحقق البرامج أولًا من قيم صغيرة مثل \(T(1)=9\) و\(T(20)=732\) و\(T(3000)=438106\) قبل حساب الهدف الكبير.

تحليل التعقيد

غربال Möbius حتى \(\lfloor\sqrt{n}\rfloor\) يكلف \(O(\sqrt{n})\) زمنًا وذاكرة. وكل دالة تراكميّة محفوظة في الذاكرة تُستدعى فقط عند قيم من الشكل \(\left\lfloor n/k \right\rfloor\)، وعدد هذه القيم المختلفة لا يتجاوز \(O(\sqrt{n})\). عمليًا، يهيمن على زمن التنفيذ عدد صغير من المجاميع التوافقية المجمعة مع استدعاءات تكرارية مخزنة مؤقتًا، وهذا يجعل \(n=10^9\) قابلًا للحساب.

مراجع

  1. صفحة المسألة: https://projecteuler.net/problem=428
  2. دالة Möbius: Wikipedia - دالة موبيوس
  3. عدد خالٍ من المربعات: Wikipedia - Squarefree integer
  4. محرف ديريشليه: Wikipedia - Dirichlet character
  5. طريقة قطع Dirichlet: Wikipedia - Dirichlet hyperbola method