Problem Summary

Problem 780 asks for the arithmetic counting function \(G(N)\) associated with toriangulations, evaluated at \(N=10^9\) modulo \(10^9+7\). The C++, Python, and Java implementations show that the geometric statement can be rewritten as an exact count of primitive directions in the triangular lattice, so the task becomes a problem about gcd filters, divisor sums, and the quadratic form \(u^2+uv+v^2\).

Instead of enumerating candidate configurations directly, the solution splits the total into an axial contribution, a strip-hyperbola contribution, and a hexagonal correction. Each of those pieces can be evaluated with arithmetic tools that scale roughly like \(O(\sqrt N\,\mathrm{polylog}\,N)\) rather than linearly in \(N\).

Mathematical Approach

The implementations reduce the final answer to an exact arithmetic decomposition. Define

$$m_0=\left\lfloor\frac N2\right\rfloor,\qquad \lambda=\left\lfloor\frac{m_0}{\sqrt3}\right\rfloor,\qquad x_0=\left\lfloor\frac N4\right\rfloor,$$

and let

$$D(t)=\sum_{k=1}^{t}\left\lfloor\frac{t}{k}\right\rfloor.$$

Then the count is organized as

$$G(N)=2D(m_0)+4\bigl(A(\lambda)-B(\lambda)\bigr)-4C(x_0).$$

Step 1: Separate the Three Arithmetic Pieces

The term \(2D(m_0)\) accounts for the two axial families that are already visible in one-dimensional divisor counting. The remaining work comes from genuinely two-dimensional directions in the triangular lattice.

The strip part is captured by two sums over positive integer pairs with \(uv\le \lambda\):

$$A(\lambda)=\sum_{uv\le \lambda}\left\lfloor\frac{N}{2\gcd(u,v)}\right\rfloor,$$

$$B(\lambda)=\sum_{uv\le \lambda}\left\lfloor\frac{\sqrt3\,uv}{\gcd(u,v)}\right\rfloor.$$

The first counts admissible widths after separating out the common divisor of \((u,v)\), while the second subtracts the exact \(\sqrt3\)-boundary coming from the same primitive data.

Step 2: Rewrite Coprimality with Möbius Inversion

To evaluate the strip sums efficiently, write

$$u=d\,a,\qquad v=d\,b,\qquad \gcd(a,b)=1.$$

For fixed \(v\) and fixed divisor \(d\mid v\), the reduced parameter \(b=v/d\) is fixed, and the admissible values of \(a\) lie in an interval determined by \(uv\le \lambda\). The coprimality condition is converted into an inclusion-exclusion sum

$$\mathbf{1}_{\gcd(a,b)=1}=\sum_{q\mid b}\mu(q)\,\mathbf{1}_{q\mid a},$$

where \(\mu\) is the Möbius function. This turns a scan over individual \(a\) values into a sum over the squarefree divisors of \(b\), which is exactly what the implementations precompute.

Step 3: Evaluate the \(\sqrt3\)-Weighted Boundary Exactly

After the Möbius step, the weighted part of the strip term becomes sums of the form

$$\sum_{k=1}^{n}\left\lfloor c\,k\sqrt3\right\rfloor,$$

for integer \(c\). These are Beatty-type floor sums. The implementations do not approximate them numerically; instead they apply an exact recursive floor-sum transformation for quadratic irrationals.

That matters because the final answer is sensitive to unit-size errors. The recursive reduction keeps everything in integer arithmetic, using only corrected integer square roots and algebraic transformations, so the strip contribution remains exact.

Step 4: Add the Hexagonal Norm Correction

The second geometric regime is controlled by the Eisenstein norm

$$Q(u,v)=u^2+uv+v^2.$$

The correction term can be written as

$$C(x_0)=D(x_0)+2\sum_{\substack{u>v\ge 1\\Q(u,v)\le x_0\\\gcd(u,v)=1\\u\not\equiv v\pmod 3}} D\!\left(\left\lfloor\frac{x_0}{Q(u,v)}\right\rfloor\right).$$

The condition \(u\not\equiv v\pmod 3\) removes the extra common factor coming from the Eisenstein prime above \(3\). In the implementations, that residue-class exclusion is handled arithmetically: first count all coprime values, then subtract the bad congruence class modulo \(3\).

Step 5: Compress the Hexagonal Sum by Quotient Plateaus

For fixed \(v\), the quotient

$$\left\lfloor\frac{x_0}{Q(u,v)}\right\rfloor$$

stays constant on intervals of \(u\). The implementations therefore move through the admissible \(u\)-range in blocks rather than one value at a time. Each block contributes a multiplicity times one cached divisor-summatory value.

This is the same philosophy as the Dirichlet hyperbola method: group equal floor quotients, evaluate each group once, and reuse the result.

Worked Example: \(N=10\)

For a small sanity check, take \(N=10\). Then

$$m_0=\left\lfloor\frac{10}{2}\right\rfloor=5,\qquad \lambda=\left\lfloor\frac{5}{\sqrt3}\right\rfloor=2,\qquad x_0=\left\lfloor\frac{10}{4}\right\rfloor=2.$$

Also

$$D(5)=\left\lfloor\frac51\right\rfloor+\left\lfloor\frac52\right\rfloor+\left\lfloor\frac53\right\rfloor+\left\lfloor\frac54\right\rfloor+\left\lfloor\frac55\right\rfloor=5+2+1+1+1=10.$$

The ordered pairs with \(uv\le 2\) are \((1,1)\), \((1,2)\), and \((2,1)\). Since all three pairs have gcd \(1\),

$$A(2)=3\left\lfloor\frac{10}{2}\right\rfloor=15.$$

For the weighted term,

$$B(2)=\lfloor\sqrt3\rfloor+\lfloor2\sqrt3\rfloor+\lfloor2\sqrt3\rfloor=1+3+3=7.$$

Finally, \(Q(2,1)=7>2\), so the inner hexagonal sum is empty and

$$C(2)=D(2)=\left\lfloor\frac21\right\rfloor+\left\lfloor\frac22\right\rfloor=3.$$

Therefore

$$G(10)=2\cdot 10+4(15-7)-4\cdot 3=20+32-12=40.$$

This example mirrors the exact decomposition used at the full input size.

How the Code Works

The C++, Python, and Java implementations first compute the three cutoffs \(m_0\), \(\lambda\), and \(x_0\). They then precompute smallest prime factors up to \(\max(\lfloor\sqrt{\lambda}\rfloor,\lfloor\sqrt{x_0}\rfloor)\), together with two divisor tables: all divisors of each integer, and squarefree divisors carrying Möbius signs.

For the strip contribution, the implementation loops over the second coordinate \(v\), splits by each divisor \(d\mid v\), and uses Möbius inversion to count reduced partners \(a\) with \(\gcd(a,v/d)=1\). The corresponding \(\sqrt3\)-weighted floor sums are evaluated by the exact Beatty recursion instead of direct iteration.

For the hexagonal correction, the implementation scans pairs ordered by \(v\) and advances \(u\) in maximal ranges where \(\left\lfloor x_0/Q(u,v)\right\rfloor\) is constant. Each distinct divisor-summatory query is cached, so repeated values are computed once and reused. The final arithmetic combination is then reduced modulo \(10^9+7\).

The three language versions follow the same mathematics. The C++ implementation additionally spreads the strip accumulation across multiple worker threads, while the Python and Java versions keep the same logic in serial form.

Complexity Analysis

Let \(P=\max(\lfloor\sqrt{\lambda}\rfloor,\lfloor\sqrt{x_0}\rfloor)\). Building the factor tables costs \(O(P\log\log P)\) time and \(O(P)\) memory. The strip and hexagonal phases both exploit divisor compression, quotient plateaus, and cached divisor-summatory values, so they avoid scanning all candidate pairs explicitly.

In practice the total running time behaves like \(O(\sqrt N\,\mathrm{polylog}\,N)\), with modest polylogarithmic factors coming from divisor enumeration and recursive floor-sum evaluation. Memory usage remains \(O(\sqrt N)\).

Footnotes and References

  1. Problem page: https://projecteuler.net/problem=780
  2. Möbius inversion formula: Wikipedia — Möbius inversion formula
  3. Divisor summatory function: Wikipedia — Divisor summatory function
  4. Beatty sequence: Wikipedia — Beatty sequence
  5. Eisenstein integer: Wikipedia — Eisenstein integer
  6. Hexagonal lattice: Wikipedia — Hexagonal lattice

Problemzusammenfassung

Problem 780 verlangt die arithmetische Zählfunktion \(G(N)\) zu den Toriangulations für \(N=10^9\) modulo \(10^9+7\). Die C++-, Python- und Java-Implementierungen zeigen, dass sich die geometrische Aufgabenstellung exakt in eine Zählung primitiver Richtungen im Dreiecksgitter umformen lässt. Dadurch wird das Problem zu einer Kombination aus ggT-Filtern, Teilersummen und der quadratischen Form \(u^2+uv+v^2\).

Statt Kandidaten direkt zu enumerieren, zerlegt die Lösung den Gesamtwert in einen axialen Anteil, einen Streifen-Hyperbel-Anteil und eine hexagonale Korrektur. Jede dieser Komponenten lässt sich mit arithmetischen Methoden berechnen, deren Aufwand ungefähr bei \(O(\sqrt N\,\mathrm{polylog}\,N)\) liegt.

Mathematischer Ansatz

Die Implementierungen reduzieren die gesuchte Anzahl auf eine exakte arithmetische Zerlegung. Definiere

$$m_0=\left\lfloor\frac N2\right\rfloor,\qquad \lambda=\left\lfloor\frac{m_0}{\sqrt3}\right\rfloor,\qquad x_0=\left\lfloor\frac N4\right\rfloor,$$

und

$$D(t)=\sum_{k=1}^{t}\left\lfloor\frac{t}{k}\right\rfloor.$$

Dann hat die Antwort die Form

$$G(N)=2D(m_0)+4\bigl(A(\lambda)-B(\lambda)\bigr)-4C(x_0).$$

Schritt 1: Die drei arithmetischen Bausteine trennen

Der Term \(2D(m_0)\) beschreibt die beiden axialen Familien, die bereits in einer eindimensionalen Teilerzählung sichtbar sind. Die übrige Arbeit stammt von echten zweidimensionalen Richtungen im Dreiecksgitter.

Der Streifenanteil wird durch zwei Summen über positive Paare mit \(uv\le \lambda\) beschrieben:

$$A(\lambda)=\sum_{uv\le \lambda}\left\lfloor\frac{N}{2\gcd(u,v)}\right\rfloor,$$

$$B(\lambda)=\sum_{uv\le \lambda}\left\lfloor\frac{\sqrt3\,uv}{\gcd(u,v)}\right\rfloor.$$

Die erste Summe zählt zulässige Breiten nach dem Herausziehen des gemeinsamen Teilers von \((u,v)\), die zweite zieht die exakte \(\sqrt3\)-Randkorrektur derselben primitiven Daten ab.

Schritt 2: Koprimalität mit Möbius-Inversion umschreiben

Zur effizienten Auswertung der Streifenterme schreiben wir

$$u=d\,a,\qquad v=d\,b,\qquad \gcd(a,b)=1.$$

Für festes \(v\) und einen festen Teiler \(d\mid v\) ist der reduzierte Parameter \(b=v/d\) konstant, und die zulässigen Werte von \(a\) liegen in einem Intervall, das aus \(uv\le \lambda\) folgt. Die Koprimalitätsbedingung wird dann per Inklusion-Exklusion zu

$$\mathbf{1}_{\gcd(a,b)=1}=\sum_{q\mid b}\mu(q)\,\mathbf{1}_{q\mid a},$$

wobei \(\mu\) die Möbius-Funktion ist. Damit ersetzt man eine Iteration über einzelne \(a\)-Werte durch eine Summe über die quadratfreien Teiler von \(b\), genau die Struktur, die die Implementierungen vorab berechnen.

Schritt 3: Die \(\sqrt3\)-gewichtete Grenze exakt summieren

Nach der Möbius-Umformung besteht der gewichtete Teil des Streifenbeitrags aus Summen der Form

$$\sum_{k=1}^{n}\left\lfloor c\,k\sqrt3\right\rfloor,$$

mit ganzzahligem \(c\). Das sind Beatty-artige Floorsummen. Die Implementierungen approximieren sie nicht numerisch, sondern verwenden eine exakte rekursive Floor-Summen-Transformation für quadratische Irrationalitäten.

Das ist entscheidend, weil der Endwert bereits auf Fehler der Größe \(1\) reagiert. Die rekursive Reduktion bleibt vollständig in ganzzahliger Arithmetik und benutzt nur korrigierte Ganzzahl-Quadratwurzeln sowie algebraische Transformationen.

Schritt 4: Die hexagonale Normkorrektur hinzufügen

Das zweite geometrische Regime wird durch die Eisenstein-Norm

$$Q(u,v)=u^2+uv+v^2$$

gesteuert. Der Korrekturterm lautet

$$C(x_0)=D(x_0)+2\sum_{\substack{u>v\ge 1\\Q(u,v)\le x_0\\\gcd(u,v)=1\\u\not\equiv v\pmod 3}} D\!\left(\left\lfloor\frac{x_0}{Q(u,v)}\right\rfloor\right).$$

Die Bedingung \(u\not\equiv v\pmod 3\) entfernt den zusätzlichen gemeinsamen Faktor, der vom Eisenstein-Primteiler über \(3\) herrührt. In den Implementierungen geschieht das rein arithmetisch: zuerst werden alle koprimen Werte gezählt, danach wird die verbotene Restklasse modulo \(3\) wieder abgezogen.

Schritt 5: Die hexagonale Summe über Quotienten-Plateaus komprimieren

Für festes \(v\) bleibt der Quotient

$$\left\lfloor\frac{x_0}{Q(u,v)}\right\rfloor$$

auf ganzen Intervallen von \(u\) konstant. Deshalb laufen die Implementierungen den zulässigen \(u\)-Bereich blockweise durch, statt jeden Wert einzeln zu behandeln. Jeder Block liefert eine Multiplizität mal genau einen zwischengespeicherten Wert der Teilersummenfunktion.

Das folgt derselben Idee wie die Dirichlet-Hyperbelmethode: gleiche Floor-Quotienten gruppieren, jede Gruppe einmal auswerten und das Ergebnis wiederverwenden.

Durchgerechnetes Beispiel: \(N=10\)

Für eine kleine Plausibilitätsprüfung setzen wir \(N=10\). Dann gilt

$$m_0=\left\lfloor\frac{10}{2}\right\rfloor=5,\qquad \lambda=\left\lfloor\frac{5}{\sqrt3}\right\rfloor=2,\qquad x_0=\left\lfloor\frac{10}{4}\right\rfloor=2.$$

Außerdem ist

$$D(5)=\left\lfloor\frac51\right\rfloor+\left\lfloor\frac52\right\rfloor+\left\lfloor\frac53\right\rfloor+\left\lfloor\frac54\right\rfloor+\left\lfloor\frac55\right\rfloor=5+2+1+1+1=10.$$

Die geordneten Paare mit \(uv\le 2\) sind \((1,1)\), \((1,2)\) und \((2,1)\). Da alle drei Paare ggT \(1\) haben, folgt

$$A(2)=3\left\lfloor\frac{10}{2}\right\rfloor=15.$$

Für den gewichteten Term erhält man

$$B(2)=\lfloor\sqrt3\rfloor+\lfloor2\sqrt3\rfloor+\lfloor2\sqrt3\rfloor=1+3+3=7.$$

Schließlich ist \(Q(2,1)=7>2\), also ist die innere hexagonale Summe leer und

$$C(2)=D(2)=\left\lfloor\frac21\right\rfloor+\left\lfloor\frac22\right\rfloor=3.$$

Damit ergibt sich

$$G(10)=2\cdot 10+4(15-7)-4\cdot 3=20+32-12=40.$$

Dieses Beispiel spiegelt exakt dieselbe Zerlegung wider, die auch für das volle Eingabewert verwendet wird.

Wie der Code arbeitet

Die C++-, Python- und Java-Implementierungen berechnen zuerst die drei Grenzwerte \(m_0\), \(\lambda\) und \(x_0\). Danach werden kleinste Primfaktoren bis \(\max(\lfloor\sqrt{\lambda}\rfloor,\lfloor\sqrt{x_0}\rfloor)\) vorab bestimmt, zusammen mit zwei Teilertabellen: allen Teilern jeder Zahl und den quadratfreien Teilern inklusive Möbius-Vorzeichen.

Für den Streifenanteil iteriert die Implementierung über die zweite Koordinate \(v\), zerlegt nach jedem Teiler \(d\mid v\) und zählt reduzierte Partner \(a\) mit \(\gcd(a,v/d)=1\) über Möbius-Inversion. Die zugehörigen \(\sqrt3\)-gewichteten Floorsummen werden mittels exakter Beatty-Rekursion statt durch lineares Durchlaufen berechnet.

Für die hexagonale Korrektur werden die Paare nach \(v\) geordnet durchlaufen, und \(u\) wird in maximalen Bereichen erhöht, in denen \(\left\lfloor x_0/Q(u,v)\right\rfloor\) konstant bleibt. Jeder verschiedene Wert der Teilersummenfunktion wird zwischengespeichert, sodass Wiederholungen nur einmal berechnet werden. Anschließend werden alle Beiträge modulo \(10^9+7\) zusammengesetzt.

Alle drei Sprachversionen folgen derselben Mathematik. Die C++-Variante verteilt den Streifenteil zusätzlich auf mehrere Threads, während Python und Java dieselbe Logik seriell ausführen.

Komplexitätsanalyse

Sei \(P=\max(\lfloor\sqrt{\lambda}\rfloor,\lfloor\sqrt{x_0}\rfloor)\). Der Aufbau der Faktortabellen kostet \(O(P\log\log P)\) Zeit und \(O(P)\) Speicher. Sowohl die Streifen- als auch die Hexagonphase nutzen Teilerkompression, Quotienten-Plateaus und zwischengespeicherte Teilersummenwerte, sodass kein vollständiger Scan aller Kandidatenpaare nötig ist.

Praktisch verhält sich die Gesamtlaufzeit wie \(O(\sqrt N\,\mathrm{polylog}\,N)\). Die zusätzlichen Faktoren stammen vor allem aus der Teilerenumeration und der rekursiven Floorsummen-Auswertung. Der Speicherbedarf bleibt bei \(O(\sqrt N)\).

Fußnoten und Referenzen

  1. Problemseite: https://projecteuler.net/problem=780
  2. Möbius-Inversionsformel: Wikipedia — Möbius inversion formula
  3. Divisor summatory function: Wikipedia — Divisor summatory function
  4. Beatty-Folge: Wikipedia — Beatty sequence
  5. Eisenstein-Ganzzahlen: Wikipedia — Eisenstein integer
  6. Hexagonales Gitter: Wikipedia — Hexagonal lattice

Problem Özeti

Problem 780, toriangulations ile ilişkili aritmetik sayım fonksiyonu \(G(N)\)'i \(N=10^9\) için \(10^9+7\) modunda istemektedir. C++, Python ve Java implementasyonları, geometrik tanımın üçgensel kafesteki ilkel yönlerin tam sayımına dönüştürülebildiğini gösteriyor. Böylece problem, EBOB filtreleri, bölen toplamları ve \(u^2+uv+v^2\) kuadratik formu üzerinden çözülen bir aritmetik probleme dönüşüyor.

Doğrudan aday üretmek yerine çözüm, toplamı eksenel katkı, şerit-hiperbol katkısı ve altıgensel düzeltme olarak ayırır. Bu parçaların her biri yaklaşık \(O(\sqrt N\,\mathrm{polylog}\,N)\) ölçüsünde çalışan aritmetik araçlarla hesaplanabilir.

Matematiksel Yaklaşım

İmplementasyonların kullandığı tam ayrışım şu şekildedir. Tanımlayalım:

$$m_0=\left\lfloor\frac N2\right\rfloor,\qquad \lambda=\left\lfloor\frac{m_0}{\sqrt3}\right\rfloor,\qquad x_0=\left\lfloor\frac N4\right\rfloor,$$

ve

$$D(t)=\sum_{k=1}^{t}\left\lfloor\frac{t}{k}\right\rfloor.$$

O zaman sonuç

$$G(N)=2D(m_0)+4\bigl(A(\lambda)-B(\lambda)\bigr)-4C(x_0)$$

şeklinde düzenlenir.

Adım 1: Üç aritmetik parçayı ayır

\(2D(m_0)\) terimi, bir boyutlu bölen sayımında zaten görülen iki eksenel aileyi temsil eder. Geriye kalan iş, üçgensel kafesteki gerçekten iki boyutlu yönlerden gelir.

Şerit kısmı, \(uv\le \lambda\) koşulunu sağlayan pozitif tamsayı çiftleri üzerinde iki toplamla ifade edilir:

$$A(\lambda)=\sum_{uv\le \lambda}\left\lfloor\frac{N}{2\gcd(u,v)}\right\rfloor,$$

$$B(\lambda)=\sum_{uv\le \lambda}\left\lfloor\frac{\sqrt3\,uv}{\gcd(u,v)}\right\rfloor.$$

İlk toplam, \((u,v)\) çiftinin ortak böleni ayrıldıktan sonra kalan izinli genişlikleri sayar. İkinci toplam ise aynı ilkel veriden gelen tam \(\sqrt3\) sınırını çıkarır.

Adım 2: Aralarında asallığı Möbius tersleme ile aç

Şerit toplamlarını verimli hesaplamak için

$$u=d\,a,\qquad v=d\,b,\qquad \gcd(a,b)=1$$

yazarız. Sabit bir \(v\) ve sabit bir \(d\mid v\) için indirgenmiş parametre \(b=v/d\) sabittir; uygun \(a\) değerleri ise \(uv\le \lambda\) koşulundan gelen bir aralık içinde kalır. Aralarında asallık koşulu şu dahil et-çıkar eşitliğiyle yazılır:

$$\mathbf{1}_{\gcd(a,b)=1}=\sum_{q\mid b}\mu(q)\,\mathbf{1}_{q\mid a},$$

burada \(\mu\) Möbius fonksiyonudur. Böylece tek tek \(a\) taramak yerine, \(b\)'nin kareden arınmış bölenleri üzerinden toplama geçilir. İmplementasyonların önceden hazırladığı yapı tam olarak budur.

Adım 3: \(\sqrt3\) ağırlıklı sınırı tam hesapla

Möbius adımından sonra şerit katkısının ağırlıklı kısmı

$$\sum_{k=1}^{n}\left\lfloor c\,k\sqrt3\right\rfloor$$

biçimindeki toplamların hesabına iner. Burada \(c\) bir tamsayıdır. Bunlar Beatty tipi floor toplamlarıdır. İmplementasyonlar bunları sayısal yaklaşıkla değil, kuadratik irrasyoneller için tam çalışan özyineli bir floor-toplam dönüşümüyle hesaplar.

Bu önemlidir; çünkü sonucun tek birimlik hata kaldırma payı yoktur. Özyineli indirgeme, düzeltilmiş tam sayı karekökleri ve cebirsel dönüşümler dışında her şeyi tam sayı aritmetiğinde tutar.

Adım 4: Altıgensel norm düzeltmesini ekle

İkinci geometrik rejim, Eisenstein normu olan

$$Q(u,v)=u^2+uv+v^2$$

ile kontrol edilir. Düzeltme terimi

$$C(x_0)=D(x_0)+2\sum_{\substack{u>v\ge 1\\Q(u,v)\le x_0\\\gcd(u,v)=1\\u\not\equiv v\pmod 3}} D\!\left(\left\lfloor\frac{x_0}{Q(u,v)}\right\rfloor\right)$$

şeklindedir. \(u\not\equiv v\pmod 3\) koşulu, \(3\) üzerindeki Eisenstein asalından gelen ek ortak çarpanı dışlar. İmplementasyonlar bunu aritmetik olarak uygular: önce tüm aralarında asal değerleri sayar, sonra mod \(3\) altında yasak kalıntı sınıfını çıkarır.

Adım 5: Altıgensel toplamı bölüm platolarıyla sıkıştır

Sabit \(v\) için

$$\left\lfloor\frac{x_0}{Q(u,v)}\right\rfloor$$

değeri, \(u\)'nun belirli aralıklarında sabit kalır. Bu yüzden implementasyonlar uygun \(u\) aralığını tek tek değil bloklar halinde ilerletir. Her blok, bir çarpan sayısı ile tek bir önbelleğe alınmış bölen-toplam değeri üretir.

Bu, Dirichlet hiperbol yöntemindeki temel fikrin aynısıdır: aynı floor bölümünü veren aralıkları grupla, her grubu bir kez hesapla ve sonucu yeniden kullan.

Çözümlü Örnek: \(N=10\)

Küçük bir doğrulama için \(N=10\) alalım. Bu durumda

$$m_0=\left\lfloor\frac{10}{2}\right\rfloor=5,\qquad \lambda=\left\lfloor\frac{5}{\sqrt3}\right\rfloor=2,\qquad x_0=\left\lfloor\frac{10}{4}\right\rfloor=2.$$

Ayrıca

$$D(5)=\left\lfloor\frac51\right\rfloor+\left\lfloor\frac52\right\rfloor+\left\lfloor\frac53\right\rfloor+\left\lfloor\frac54\right\rfloor+\left\lfloor\frac55\right\rfloor=5+2+1+1+1=10.$$

\(uv\le 2\) koşulunu sağlayan sıralı çiftler \((1,1)\), \((1,2)\) ve \((2,1)\)'dir. Üçünün de EBOB'u \(1\) olduğundan

$$A(2)=3\left\lfloor\frac{10}{2}\right\rfloor=15.$$

Ağırlıklı toplam için

$$B(2)=\lfloor\sqrt3\rfloor+\lfloor2\sqrt3\rfloor+\lfloor2\sqrt3\rfloor=1+3+3=7.$$

Son olarak \(Q(2,1)=7>2\) olduğu için iç altıgensel toplam boştur ve

$$C(2)=D(2)=\left\lfloor\frac21\right\rfloor+\left\lfloor\frac22\right\rfloor=3.$$

Dolayısıyla

$$G(10)=2\cdot 10+4(15-7)-4\cdot 3=20+32-12=40.$$

Bu küçük örnek, gerçek girdi boyutunda kullanılan ayrışımın aynısını gösterir.

Kod Nasıl Çalışıyor

C++, Python ve Java implementasyonları önce \(m_0\), \(\lambda\) ve \(x_0\) eşiklerini hesaplar. Ardından \(\max(\lfloor\sqrt{\lambda}\rfloor,\lfloor\sqrt{x_0}\rfloor)\) sınırına kadar en küçük asal bölen tablosu hazırlanır; buna ek olarak her sayı için tüm bölenler ve Möbius işaretli kareden arınmış bölenler önceden üretilir.

Şerit katkısı için implementasyon ikinci koordinat \(v\) üzerinde döner, her \(d\mid v\) bölenine göre parçalar ve \(\gcd(a,v/d)=1\) koşulunu Möbius tersleme ile sayar. Karşılık gelen \(\sqrt3\) ağırlıklı floor toplamları ise doğrusal tarama yerine tam Beatty özyinelemesiyle hesaplanır.

Altıgensel düzeltmede çiftler \(v\)'ye göre sıralı biçimde gezilir ve \(\left\lfloor x_0/Q(u,v)\right\rfloor\) ifadesinin sabit kaldığı en büyük \(u\) blokları üzerinden ilerlenir. Her farklı bölen-toplam sorgusu önbelleğe alınır; böylece tekrar eden değerler tek kez hesaplanır. Son adımda tüm katkılar \(10^9+7\) modunda birleştirilir.

Üç dildeki sürüm aynı matematiksel planı izler. C++ sürümü buna ek olarak şerit toplamını birden fazla iş parçacığına dağıtır; Python ve Java sürümleri aynı mantığı seri olarak uygular.

Karmaşıklık Analizi

\(P=\max(\lfloor\sqrt{\lambda}\rfloor,\lfloor\sqrt{x_0}\rfloor)\) olsun. Çarpan tablolarının kurulması \(O(P\log\log P)\) zaman ve \(O(P)\) bellek ister. Hem şerit hem de altıgensel faz, bölen sıkıştırması, bölüm platoları ve önbelleğe alınmış bölen-toplam değerleri kullandığı için tüm aday çiftleri tek tek taramaz.

Pratikte toplam süre \(O(\sqrt N\,\mathrm{polylog}\,N)\) davranışı gösterir. Ek çarpanlar daha çok bölenleri dolaşma ve özyineli floor-toplam değerlendirmesinden gelir. Bellek kullanımı \(O(\sqrt N)\) düzeyinde kalır.

Dipnotlar ve Kaynaklar

  1. Problem sayfası: https://projecteuler.net/problem=780
  2. Möbius tersleme formülü: Wikipedia — Möbius inversion formula
  3. Bölen toplam fonksiyonu: Wikipedia — Divisor summatory function
  4. Beatty dizisi: Wikipedia — Beatty sequence
  5. Eisenstein tamsayıları: Wikipedia — Eisenstein integer
  6. Altıgensel kafes: Wikipedia — Hexagonal lattice

Resumen del Problema

El problema 780 pide la función de conteo aritmético \(G(N)\) asociada a las toriangulations, evaluada en \(N=10^9\) módulo \(10^9+7\). Las implementaciones en C++, Python y Java muestran que la formulación geométrica puede reescribirse como un conteo exacto de direcciones primitivas en la red triangular. Por eso la tarea acaba expresándose en filtros de mcd, sumas de divisores y la forma cuadrática \(u^2+uv+v^2\).

En lugar de enumerar configuraciones directamente, la solución divide el total en una contribución axial, una contribución de franja-hipérbola y una corrección hexagonal. Cada bloque puede evaluarse con herramientas aritméticas cuyo costo es aproximadamente \(O(\sqrt N\,\mathrm{polylog}\,N)\).

Enfoque Matemático

Las implementaciones reducen la respuesta final a una descomposición aritmética exacta. Definimos

$$m_0=\left\lfloor\frac N2\right\rfloor,\qquad \lambda=\left\lfloor\frac{m_0}{\sqrt3}\right\rfloor,\qquad x_0=\left\lfloor\frac N4\right\rfloor,$$

y

$$D(t)=\sum_{k=1}^{t}\left\lfloor\frac{t}{k}\right\rfloor.$$

Entonces la cantidad buscada se organiza como

$$G(N)=2D(m_0)+4\bigl(A(\lambda)-B(\lambda)\bigr)-4C(x_0).$$

Paso 1: Separar las tres piezas aritméticas

El término \(2D(m_0)\) recoge las dos familias axiales que ya aparecen en el conteo unidimensional de divisores. El resto del trabajo proviene de direcciones verdaderamente bidimensionales en la red triangular.

La parte de franja se describe mediante dos sumas sobre pares positivos con \(uv\le \lambda\):

$$A(\lambda)=\sum_{uv\le \lambda}\left\lfloor\frac{N}{2\gcd(u,v)}\right\rfloor,$$

$$B(\lambda)=\sum_{uv\le \lambda}\left\lfloor\frac{\sqrt3\,uv}{\gcd(u,v)}\right\rfloor.$$

La primera cuenta anchuras admisibles después de separar el divisor común de \((u,v)\), y la segunda resta el borde exacto con factor \(\sqrt3\) asociado a los mismos datos primitivos.

Paso 2: Reescribir la coprimalidad con inversión de Möbius

Para evaluar con eficiencia las sumas de la franja, escribimos

$$u=d\,a,\qquad v=d\,b,\qquad \gcd(a,b)=1.$$

Si fijamos \(v\) y un divisor \(d\mid v\), el parámetro reducido \(b=v/d\) queda fijado, y los valores admisibles de \(a\) viven en un intervalo determinado por \(uv\le \lambda\). La condición de coprimalidad se transforma en

$$\mathbf{1}_{\gcd(a,b)=1}=\sum_{q\mid b}\mu(q)\,\mathbf{1}_{q\mid a},$$

donde \(\mu\) es la función de Möbius. Así se sustituye un barrido sobre valores individuales de \(a\) por una suma sobre divisores libres de cuadrados de \(b\), exactamente la estructura que las implementaciones precalculan.

Paso 3: Calcular exactamente el borde ponderado por \(\sqrt3\)

Después del paso de Möbius, la parte ponderada de la franja se convierte en sumas del tipo

$$\sum_{k=1}^{n}\left\lfloor c\,k\sqrt3\right\rfloor,$$

con \(c\) entero. Son sumas de tipo Beatty. Las implementaciones no las aproximan numéricamente; aplican una transformación recursiva exacta de sumas con piso para irracionales cuadráticos.

Eso importa porque el resultado final no tolera errores de tamaño \(1\). La reducción recursiva mantiene todo en aritmética entera, usando únicamente raíces cuadradas enteras corregidas y transformaciones algebraicas.

Paso 4: Añadir la corrección de norma hexagonal

El segundo régimen geométrico está gobernado por la norma de Eisenstein

$$Q(u,v)=u^2+uv+v^2.$$

La corrección puede escribirse como

$$C(x_0)=D(x_0)+2\sum_{\substack{u>v\ge 1\\Q(u,v)\le x_0\\\gcd(u,v)=1\\u\not\equiv v\pmod 3}} D\!\left(\left\lfloor\frac{x_0}{Q(u,v)}\right\rfloor\right).$$

La condición \(u\not\equiv v\pmod 3\) elimina el factor común extra asociado al primo de Eisenstein por encima de \(3\). En las implementaciones esto se hace aritméticamente: primero se cuentan todos los valores coprimos y después se resta la clase de residuos prohibida módulo \(3\).

Paso 5: Comprimir la suma hexagonal por mesetas de cociente

Para \(v\) fijo, el cociente

$$\left\lfloor\frac{x_0}{Q(u,v)}\right\rfloor$$

permanece constante en intervalos de \(u\). Por eso las implementaciones recorren el rango admisible de \(u\) por bloques en lugar de hacerlo valor por valor. Cada bloque aporta una multiplicidad por un único valor memorizado de la función sumatoria de divisores.

La idea es la misma que en el método de la hipérbola de Dirichlet: agrupar cocientes iguales, evaluar cada grupo una sola vez y reutilizar el resultado.

Ejemplo trabajado: \(N=10\)

Para una comprobación pequeña, tomemos \(N=10\). Entonces

$$m_0=\left\lfloor\frac{10}{2}\right\rfloor=5,\qquad \lambda=\left\lfloor\frac{5}{\sqrt3}\right\rfloor=2,\qquad x_0=\left\lfloor\frac{10}{4}\right\rfloor=2.$$

Además,

$$D(5)=\left\lfloor\frac51\right\rfloor+\left\lfloor\frac52\right\rfloor+\left\lfloor\frac53\right\rfloor+\left\lfloor\frac54\right\rfloor+\left\lfloor\frac55\right\rfloor=5+2+1+1+1=10.$$

Los pares ordenados con \(uv\le 2\) son \((1,1)\), \((1,2)\) y \((2,1)\). Como los tres tienen mcd igual a \(1\), se obtiene

$$A(2)=3\left\lfloor\frac{10}{2}\right\rfloor=15.$$

Para el término ponderado,

$$B(2)=\lfloor\sqrt3\rfloor+\lfloor2\sqrt3\rfloor+\lfloor2\sqrt3\rfloor=1+3+3=7.$$

Por último, \(Q(2,1)=7>2\), así que la suma hexagonal interior es vacía y

$$C(2)=D(2)=\left\lfloor\frac21\right\rfloor+\left\lfloor\frac22\right\rfloor=3.$$

Por tanto,

$$G(10)=2\cdot 10+4(15-7)-4\cdot 3=20+32-12=40.$$

Este ejemplo reproduce exactamente la misma descomposición que se usa para la entrada grande.

Cómo Funciona el Código

Las implementaciones en C++, Python y Java calculan primero los tres umbrales \(m_0\), \(\lambda\) y \(x_0\). Después precalculan los menores factores primos hasta \(\max(\lfloor\sqrt{\lambda}\rfloor,\lfloor\sqrt{x_0}\rfloor)\), junto con dos tablas de divisores: todos los divisores de cada entero y los divisores libres de cuadrados con signo de Möbius.

Para la contribución de franja, la implementación recorre la segunda coordenada \(v\), se descompone por cada divisor \(d\mid v\) y usa inversión de Möbius para contar compañeros reducidos \(a\) con \(\gcd(a,v/d)=1\). Las sumas con piso ponderadas por \(\sqrt3\) se evalúan mediante la recursión exacta de Beatty, no mediante iteración directa.

Para la corrección hexagonal, la implementación examina pares ordenados por \(v\) y avanza \(u\) por intervalos máximos donde \(\left\lfloor x_0/Q(u,v)\right\rfloor\) permanece constante. Cada consulta distinta a la función sumatoria de divisores se guarda en caché, de modo que los valores repetidos se calculan una sola vez. Al final, la combinación aritmética completa se reduce módulo \(10^9+7\).

Las tres versiones siguen la misma matemática. La implementación en C++ además paraleliza la acumulación de la parte de franja, mientras que Python y Java mantienen la misma lógica en forma secuencial.

Análisis de Complejidad

Sea \(P=\max(\lfloor\sqrt{\lambda}\rfloor,\lfloor\sqrt{x_0}\rfloor)\). Construir las tablas de factorización cuesta \(O(P\log\log P)\) tiempo y \(O(P)\) memoria. Tanto la fase de franja como la fase hexagonal explotan compresión por divisores, mesetas de cociente y valores en caché de la función sumatoria de divisores, evitando un barrido explícito de todos los pares candidatos.

En la práctica, el tiempo total se comporta como \(O(\sqrt N\,\mathrm{polylog}\,N)\). Los factores extra vienen sobre todo de enumerar divisores y de evaluar recursivamente las sumas con piso. El uso de memoria se mantiene en \(O(\sqrt N)\).

Notas y Referencias

  1. Página del problema: https://projecteuler.net/problem=780
  2. Fórmula de inversión de Möbius: Wikipedia — Möbius inversion formula
  3. Función sumatoria de divisores: Wikipedia — Divisor summatory function
  4. Sucesión de Beatty: Wikipedia — Beatty sequence
  5. Entero de Eisenstein: Wikipedia — Eisenstein integer
  6. Red hexagonal: Wikipedia — Hexagonal lattice

问题概述

第 780 题要求计算与 toriangulations 对应的算术计数函数 \(G(N)\),并在 \(N=10^9\) 时取模 \(10^9+7\)。从 C++、Python 和 Java 的实现可以看出,原始的几何描述能够被严格改写成三角晶格中“原始方向”的计数问题,因此核心变成了 gcd 过滤、除数求和以及二次型 \(u^2+uv+v^2\) 的处理。

解法并不直接枚举候选构型,而是把总数拆成三个部分:轴向贡献、条带-双曲线贡献,以及六角区域修正。这样每一部分都可以用数论工具精确计算,整体代价大约是 \(O(\sqrt N\,\mathrm{polylog}\,N)\),远小于按 \(N\) 线性扫描。

数学方法

实现最终使用的是一个完全精确的算术分解。定义

$$m_0=\left\lfloor\frac N2\right\rfloor,\qquad \lambda=\left\lfloor\frac{m_0}{\sqrt3}\right\rfloor,\qquad x_0=\left\lfloor\frac N4\right\rfloor,$$

再设

$$D(t)=\sum_{k=1}^{t}\left\lfloor\frac{t}{k}\right\rfloor.$$

那么答案整理成

$$G(N)=2D(m_0)+4\bigl(A(\lambda)-B(\lambda)\bigr)-4C(x_0).$$

步骤 1:先把三个算术部分分开

\(2D(m_0)\) 对应两族轴向情形,这一部分本质上已经退化为一维的除数计数。真正困难的部分来自三角晶格里的二维原始方向。

条带部分由满足 \(uv\le \lambda\) 的正整数对给出:

$$A(\lambda)=\sum_{uv\le \lambda}\left\lfloor\frac{N}{2\gcd(u,v)}\right\rfloor,$$

$$B(\lambda)=\sum_{uv\le \lambda}\left\lfloor\frac{\sqrt3\,uv}{\gcd(u,v)}\right\rfloor.$$

第一项表示把 \((u,v)\) 的公共因子提出之后还能容纳多少合法宽度,第二项则减去同一批原始数据对应的精确 \(\sqrt3\) 边界。

步骤 2:用 Möbius 反演处理互素条件

为了高效求条带和,把

$$u=d\,a,\qquad v=d\,b,\qquad \gcd(a,b)=1$$

写出来。固定 \(v\) 和一个除数 \(d\mid v\) 之后,约化参数 \(b=v/d\) 就固定了,而允许的 \(a\) 落在由 \(uv\le \lambda\) 决定的某个区间里。互素条件可以改写成

$$\mathbf{1}_{\gcd(a,b)=1}=\sum_{q\mid b}\mu(q)\,\mathbf{1}_{q\mid a},$$

其中 \(\mu\) 是 Möbius 函数。这样就不需要逐个检查 \(a\),而是转化为对 \(b\) 的平方自由除数做容斥,这正是实现里预处理出来的数据结构。

步骤 3:精确处理带有 \(\sqrt3\) 的边界和

在 Möbius 反演之后,条带部分里带权的那一项会变成

$$\sum_{k=1}^{n}\left\lfloor c\,k\sqrt3\right\rfloor$$

这样的和,其中 \(c\) 是整数。这类对象属于 Beatty 型 floor 和。实现并没有把它当成浮点近似问题,而是使用了针对二次无理数的精确递归 floor-sum 变换。

这一步必须精确,因为最终答案对单位级误差都敏感。递归过程中始终保留整数算术,只使用经过修正的整数平方根和代数等价变换,因此不会引入近似误差。

步骤 4:加入六角范数修正项

第二个几何区域由 Eisenstein 范数

$$Q(u,v)=u^2+uv+v^2$$

控制。对应的修正项是

$$C(x_0)=D(x_0)+2\sum_{\substack{u>v\ge 1\\Q(u,v)\le x_0\\\gcd(u,v)=1\\u\not\equiv v\pmod 3}} D\!\left(\left\lfloor\frac{x_0}{Q(u,v)}\right\rfloor\right).$$

条件 \(u\not\equiv v\pmod 3\) 用来去掉与 \(3\) 上方 Eisenstein 素元相关的额外公共因子。实现中的做法是先数出所有互素情形,再减去模 \(3\) 下那一类不合法的剩余类。

步骤 5:按商值平台压缩六角和

对固定的 \(v\),商

$$\left\lfloor\frac{x_0}{Q(u,v)}\right\rfloor$$

在一整段 \(u\) 上都是常数。因此实现并不是一个个地推进 \(u\),而是一次跳过一个最大区间。每个区间只贡献“该区间里的点数”乘上一个缓存好的除数求和值。

这与 Dirichlet 双曲线方法的思想一致:把相同的 floor 商分组,每组只算一次,再把结果复用。

算例:\(N=10\)

为了做一个小型核对,取 \(N=10\)。这时

$$m_0=\left\lfloor\frac{10}{2}\right\rfloor=5,\qquad \lambda=\left\lfloor\frac{5}{\sqrt3}\right\rfloor=2,\qquad x_0=\left\lfloor\frac{10}{4}\right\rfloor=2.$$

同时

$$D(5)=\left\lfloor\frac51\right\rfloor+\left\lfloor\frac52\right\rfloor+\left\lfloor\frac53\right\rfloor+\left\lfloor\frac54\right\rfloor+\left\lfloor\frac55\right\rfloor=5+2+1+1+1=10.$$

满足 \(uv\le 2\) 的有序整数对只有 \((1,1)\)、\((1,2)\)、\((2,1)\)。这三个对的 gcd 都等于 \(1\),所以

$$A(2)=3\left\lfloor\frac{10}{2}\right\rfloor=15.$$

加权项则为

$$B(2)=\lfloor\sqrt3\rfloor+\lfloor2\sqrt3\rfloor+\lfloor2\sqrt3\rfloor=1+3+3=7.$$

再看六角修正,因为 \(Q(2,1)=7>2\),所以内部求和为空,于是

$$C(2)=D(2)=\left\lfloor\frac21\right\rfloor+\left\lfloor\frac22\right\rfloor=3.$$

最终得到

$$G(10)=2\cdot 10+4(15-7)-4\cdot 3=20+32-12=40.$$

这个小例子和大规模输入时采用的公式分解完全一致。

代码如何工作

C++、Python 和 Java 实现首先计算三个阈值 \(m_0\)、\(\lambda\) 和 \(x_0\)。然后它们会预处理到 \(\max(\lfloor\sqrt{\lambda}\rfloor,\lfloor\sqrt{x_0}\rfloor)\) 为止的最小素因子,并同时构建两张表:每个整数的全部除数,以及带 Möbius 符号的平方自由除数。

在条带部分,实现按第二坐标 \(v\) 枚举,再按每个 \(d\mid v\) 分块,通过 Möbius 反演统计满足 \(\gcd(a,v/d)=1\) 的约化伙伴 \(a\)。与之对应的 \(\sqrt3\) 加权 floor 和,不是直接循环求和,而是交给精确的 Beatty 递归去完成。

在六角修正部分,实现按照 \(v\) 扫描,并让 \(u\) 在 \(\left\lfloor x_0/Q(u,v)\right\rfloor\) 保持常数的最大区间上跳跃。每一种不同的除数求和查询都会被缓存,因此重复值只计算一次。最后再把所有算术贡献组合起来,并对 \(10^9+7\) 取模。

三种语言的实现遵循完全相同的数学结构。C++ 版本另外把条带累加拆分到多个工作线程上,而 Python 和 Java 版本则保持串行但逻辑一致。

复杂度分析

记 \(P=\max(\lfloor\sqrt{\lambda}\rfloor,\lfloor\sqrt{x_0}\rfloor)\)。因子表预处理需要 \(O(P\log\log P)\) 时间和 \(O(P)\) 内存。条带阶段和六角阶段都利用了除数压缩、商值平台以及缓存的除数求和,因此避免了对所有候选点对进行显式扫描。

实际运行时间大致表现为 \(O(\sqrt N\,\mathrm{polylog}\,N)\)。额外的多对数因子主要来自除数枚举和递归 floor-sum 计算。内存占用保持在 \(O(\sqrt N)\) 量级。

脚注与参考资料

  1. 题目页面: https://projecteuler.net/problem=780
  2. Möbius 反演公式: Wikipedia — Möbius inversion formula
  3. 除数求和函数: Wikipedia — Divisor summatory function
  4. Beatty 序列: Wikipedia — Beatty sequence
  5. Eisenstein 整数: Wikipedia — Eisenstein integer
  6. 六角晶格: Wikipedia — Hexagonal lattice

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

В задаче 780 требуется вычислить арифметическую функцию подсчета \(G(N)\), связанную с toriangulations, при \(N=10^9\) по модулю \(10^9+7\). Из реализаций на C++, Python и Java видно, что исходную геометрическую постановку можно точно переписать как подсчет примитивных направлений в треугольной решетке. Поэтому вся трудность сводится к фильтрам по НОД, суммам по делителям и квадратичной форме \(u^2+uv+v^2\).

Вместо прямого перебора конфигураций решение раскладывает ответ на осевой вклад, вклад полосы-гиперболы и шестиугольную поправку. Каждая часть вычисляется чисто арифметическими методами с трудоемкостью порядка \(O(\sqrt N\,\mathrm{polylog}\,N)\).

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

Реализации сводят ответ к точному арифметическому разложению. Введем

$$m_0=\left\lfloor\frac N2\right\rfloor,\qquad \lambda=\left\lfloor\frac{m_0}{\sqrt3}\right\rfloor,\qquad x_0=\left\lfloor\frac N4\right\rfloor,$$

а также

$$D(t)=\sum_{k=1}^{t}\left\lfloor\frac{t}{k}\right\rfloor.$$

Тогда искомая величина записывается как

$$G(N)=2D(m_0)+4\bigl(A(\lambda)-B(\lambda)\bigr)-4C(x_0).$$

Шаг 1: Выделить три арифметических блока

Слагаемое \(2D(m_0)\) отвечает двум осевым семействам, которые уже видны в одномерном подсчете делителей. Оставшаяся часть связана с настоящими двумерными направлениями в треугольной решетке.

Полосовая часть задается двумя суммами по положительным парам с условием \(uv\le \lambda\):

$$A(\lambda)=\sum_{uv\le \lambda}\left\lfloor\frac{N}{2\gcd(u,v)}\right\rfloor,$$

$$B(\lambda)=\sum_{uv\le \lambda}\left\lfloor\frac{\sqrt3\,uv}{\gcd(u,v)}\right\rfloor.$$

Первая сумма считает допустимые ширины после выделения общего делителя пары \((u,v)\), а вторая вычитает точную \(\sqrt3\)-границу, зависящую от тех же примитивных данных.

Шаг 2: Заменить взаимную простоту инверсией Мёбиуса

Чтобы эффективно вычислять полосовые суммы, представим

$$u=d\,a,\qquad v=d\,b,\qquad \gcd(a,b)=1.$$

При фиксированных \(v\) и делителе \(d\mid v\) приведенный параметр \(b=v/d\) фиксирован, а допустимые значения \(a\) лежат в интервале, задаваемом неравенством \(uv\le \lambda\). Условие взаимной простоты переписывается как

$$\mathbf{1}_{\gcd(a,b)=1}=\sum_{q\mid b}\mu(q)\,\mathbf{1}_{q\mid a},$$

где \(\mu\) — функция Мёбиуса. Это заменяет перебор отдельных значений \(a\) суммированием по квадратсвободным делителям \(b\), именно такую структуру реализации и строят заранее.

Шаг 3: Точно посчитать \(\sqrt3\)-взвешенную границу

После применения инверсии Мёбиуса взвешенная часть полосового вклада превращается в суммы вида

$$\sum_{k=1}^{n}\left\lfloor c\,k\sqrt3\right\rfloor,$$

где \(c\) — целое число. Это floor-суммы типа Битти. Реализации не используют численные приближения, а применяют точное рекурсивное преобразование для квадратичных иррациональностей.

Это важно, потому что итоговая формула чувствительна даже к ошибке размера \(1\). Рекурсия сохраняет всю арифметику в целых числах и использует только корректированные целочисленные квадратные корни и алгебраические преобразования.

Шаг 4: Добавить поправку с шестиугольной нормой

Второй геометрический режим управляется нормой Эйзенштейна

$$Q(u,v)=u^2+uv+v^2.$$

Поправочный член имеет вид

$$C(x_0)=D(x_0)+2\sum_{\substack{u>v\ge 1\\Q(u,v)\le x_0\\\gcd(u,v)=1\\u\not\equiv v\pmod 3}} D\!\left(\left\lfloor\frac{x_0}{Q(u,v)}\right\rfloor\right).$$

Условие \(u\not\equiv v\pmod 3\) убирает лишний общий множитель, связанный с простым элементом Эйзенштейна над \(3\). В реализациях это делается арифметически: сначала считаются все взаимно простые пары, затем вычитается запрещенный класс вычетов по модулю \(3\).

Шаг 5: Сжать шестиугольную сумму по плато частного

При фиксированном \(v\) выражение

$$\left\lfloor\frac{x_0}{Q(u,v)}\right\rfloor$$

остается постоянным на целых интервалах по \(u\). Поэтому реализации двигаются не по одному значению \(u\), а сразу по максимальным блокам. Каждый блок дает кратность, умноженную на одно закэшированное значение сумматорной функции делителей.

Это та же идея, что и в методе гиперболы Дирихле: сгруппировать одинаковые floor-частные, вычислить каждую группу один раз и переиспользовать результат.

Разобранный пример: \(N=10\)

Для небольшой проверки возьмем \(N=10\). Тогда

$$m_0=\left\lfloor\frac{10}{2}\right\rfloor=5,\qquad \lambda=\left\lfloor\frac{5}{\sqrt3}\right\rfloor=2,\qquad x_0=\left\lfloor\frac{10}{4}\right\rfloor=2.$$

Кроме того,

$$D(5)=\left\lfloor\frac51\right\rfloor+\left\lfloor\frac52\right\rfloor+\left\lfloor\frac53\right\rfloor+\left\lfloor\frac54\right\rfloor+\left\lfloor\frac55\right\rfloor=5+2+1+1+1=10.$$

Упорядоченные пары с \(uv\le 2\) — это \((1,1)\), \((1,2)\) и \((2,1)\). У всех трех пар НОД равен \(1\), поэтому

$$A(2)=3\left\lfloor\frac{10}{2}\right\rfloor=15.$$

Для взвешенного члена получаем

$$B(2)=\lfloor\sqrt3\rfloor+\lfloor2\sqrt3\rfloor+\lfloor2\sqrt3\rfloor=1+3+3=7.$$

Наконец, \(Q(2,1)=7>2\), значит внутренняя шестиугольная сумма пуста и

$$C(2)=D(2)=\left\lfloor\frac21\right\rfloor+\left\lfloor\frac22\right\rfloor=3.$$

Следовательно,

$$G(10)=2\cdot 10+4(15-7)-4\cdot 3=20+32-12=40.$$

Этот пример повторяет ту же самую точную схему разложения, которая используется и для полного входного значения.

Как работает код

Реализации на C++, Python и Java сначала вычисляют три порога \(m_0\), \(\lambda\) и \(x_0\). Затем они предвычисляют минимальные простые делители до \(\max(\lfloor\sqrt{\lambda}\rfloor,\lfloor\sqrt{x_0}\rfloor)\), а также строят две таблицы: все делители каждого числа и квадратсвободные делители со знаками Мёбиуса.

Для полосового вклада реализация перебирает вторую координату \(v\), разбивает задачу по каждому делителю \(d\mid v\) и через инверсию Мёбиуса считает приведенных партнеров \(a\), удовлетворяющих \(\gcd(a,v/d)=1\). Соответствующие \(\sqrt3\)-взвешенные floor-суммы вычисляются точной рекурсией Битти, а не прямым проходом.

Для шестиугольной поправки реализация просматривает пары в порядке по \(v\) и увеличивает \(u\) сразу на максимальные отрезки, на которых \(\left\lfloor x_0/Q(u,v)\right\rfloor\) постоянно. Каждый новый запрос к сумматорной функции делителей кэшируется, поэтому повторяющиеся значения считаются только один раз. После этого все арифметические части объединяются по модулю \(10^9+7\).

Во всех трех языках используется одна и та же математика. C++-версия дополнительно распараллеливает накопление полосового вклада, тогда как версии на Python и Java сохраняют ту же логику в последовательном виде.

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

Пусть \(P=\max(\lfloor\sqrt{\lambda}\rfloor,\lfloor\sqrt{x_0}\rfloor)\). Построение факторизационных таблиц требует \(O(P\log\log P)\) времени и \(O(P)\) памяти. И полосовая, и шестиугольная фазы используют сжатие по делителям, плато равных частных и кэширование значений сумматорной функции делителей, поэтому не перебирают все пары явно.

На практике полное время работы ведет себя как \(O(\sqrt N\,\mathrm{polylog}\,N)\). Дополнительные множители появляются главным образом из перечисления делителей и рекурсивного вычисления floor-сумм. Память остается порядка \(O(\sqrt N)\).

Сноски и ссылки

  1. Страница задачи: https://projecteuler.net/problem=780
  2. Формула инверсии Мёбиуса: Wikipedia — Möbius inversion formula
  3. Сумматорная функция делителей: Wikipedia — Divisor summatory function
  4. Последовательность Битти: Wikipedia — Beatty sequence
  5. Целые числа Эйзенштейна: Wikipedia — Eisenstein integer
  6. Шестиугольная решетка: Wikipedia — Hexagonal lattice

ملخص المسألة

تطلب المسألة 780 حساب دالة العد الحسابية \(G(N)\) المرتبطة بـ toriangulations عند \(N=10^9\) بترديد \(10^9+7\). وتُظهر تطبيقات ++C وPython وJava أن الصياغة الهندسية يمكن تحويلها بدقة إلى عدّ للاتجاهات الأولية في الشبكة المثلثية. لذلك تصبح المسألة معتمدة على مرشحات القاسم المشترك الأكبر، ومجاميع القواسم، والصيغة التربيعية \(u^2+uv+v^2\).

بدلًا من تعداد جميع التراكيب مباشرة، يقسم الحل المجموع إلى مساهمة محورية، ومساهمة شريطية-زائدية، وتصحيح سداسي. ويمكن حساب كل جزء بأدوات عددية تعمل تقريبًا بكلفة \(O(\sqrt N\,\mathrm{polylog}\,N)\) بدلًا من كلفة خطية في \(N\).

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

تختزل التطبيقات الجواب إلى تفكيك حسابي دقيق. نعرّف

$$m_0=\left\lfloor\frac N2\right\rfloor,\qquad \lambda=\left\lfloor\frac{m_0}{\sqrt3}\right\rfloor,\qquad x_0=\left\lfloor\frac N4\right\rfloor,$$

وكذلك

$$D(t)=\sum_{k=1}^{t}\left\lfloor\frac{t}{k}\right\rfloor.$$

وعندها يُرتَّب الناتج على الصورة

$$G(N)=2D(m_0)+4\bigl(A(\lambda)-B(\lambda)\bigr)-4C(x_0).$$

الخطوة 1: فصل الأجزاء الحسابية الثلاثة

الحد \(2D(m_0)\) يمثّل العائلتين المحوريتين اللتين تظهران أصلًا في عدّ القواسم أحادي البعد. أما بقية العمل فتأتي من الاتجاهات الثنائية الحقيقية داخل الشبكة المثلثية.

جزء الشريط يُكتب على شكل مجموعين فوق الأزواج الموجبة التي تحقق \(uv\le \lambda\):

$$A(\lambda)=\sum_{uv\le \lambda}\left\lfloor\frac{N}{2\gcd(u,v)}\right\rfloor,$$

$$B(\lambda)=\sum_{uv\le \lambda}\left\lfloor\frac{\sqrt3\,uv}{\gcd(u,v)}\right\rfloor.$$

المجموع الأول يحسب السعات المسموح بها بعد فصل القاسم المشترك للزوج \((u,v)\)، أما الثاني فيطرح الحدّ الدقيق الموزون بـ \(\sqrt3\) الناتج عن البنية الأولية نفسها.

الخطوة 2: إعادة كتابة شرط التباين النسبي بعكس Möbius

لكي نحسب مجاميع الشريط بكفاءة نكتب

$$u=d\,a,\qquad v=d\,b,\qquad \gcd(a,b)=1.$$

عند تثبيت \(v\) ومقسوم \(d\mid v\)، يصبح المعامل المختزل \(b=v/d\) ثابتًا، وتقع القيم المسموح بها لـ \(a\) داخل مجال يحدده الشرط \(uv\le \lambda\). ويُحوَّل شرط التباين النسبي إلى

$$\mathbf{1}_{\gcd(a,b)=1}=\sum_{q\mid b}\mu(q)\,\mathbf{1}_{q\mid a},$$

حيث \(\mu\) هي دالة Möbius. وبهذا ننتقل من فحص كل قيمة من \(a\) على حدة إلى جمع فوق القواسم الخالية من المربعات لـ \(b\)، وهي بالضبط البنية التي تهيئها التطبيقات مسبقًا.

الخطوة 3: حساب الحدّ الموزون بـ \(\sqrt3\) حسابًا دقيقًا

بعد خطوة Möbius، يتحول الجزء الموزون من مساهمة الشريط إلى مجاميع من الشكل

$$\sum_{k=1}^{n}\left\lfloor c\,k\sqrt3\right\rfloor,$$

حيث \(c\) عدد صحيح. هذه مجاميع من نوع Beatty. لا تعتمد التطبيقات على تقريب عددي عائم، بل تستخدم تحويلًا تكراريًا دقيقًا لمجاميع الجزء الصحيح الخاصة باللاعداد الجذرية التربيعية.

وهذا ضروري لأن الجواب النهائي حساس حتى لخطأ مقداره \(1\). وتبقى الحسابات داخل الحساب الصحيح، مع استعمال جذور تربيعية صحيحة مصححة وتحويلات جبرية مكافئة فقط.

الخطوة 4: إضافة تصحيح المعيار السداسي

النظام الهندسي الثاني تحكمه قيمة معيار أيزنشتاين

$$Q(u,v)=u^2+uv+v^2.$$

ويُكتب حدّ التصحيح على الصورة

$$C(x_0)=D(x_0)+2\sum_{\substack{u>v\ge 1\\Q(u,v)\le x_0\\\gcd(u,v)=1\\u\not\equiv v\pmod 3}} D\!\left(\left\lfloor\frac{x_0}{Q(u,v)}\right\rfloor\right).$$

الشرط \(u\not\equiv v\pmod 3\) يزيل العامل المشترك الإضافي المرتبط بالأولية الأيزنشتاينية فوق \(3\). وفي التطبيقات يُنفَّذ هذا حسابيًا: أولًا يُحصى جميع الأزواج المتباينة نسبيًا، ثم تُطرح فئة البواقي الممنوعة modulo \(3\).

الخطوة 5: ضغط المجموع السداسي بواسطة هضاب خارج القسمة

عند تثبيت \(v\)، تبقى القيمة

$$\left\lfloor\frac{x_0}{Q(u,v)}\right\rfloor$$

ثابتة على فترات كاملة من \(u\). لذلك لا تتحرك التطبيقات قيمةً قيمةً، بل تقفز عبر أكبر كتل يبقى فيها خارج القسمة ثابتًا. كل كتلة تعطي مضاعفًا عدديًا مضروبًا في قيمة واحدة مخزنة لدالة مجموع القواسم.

هذه هي الفكرة نفسها في طريقة القطع الزائد لـ Dirichlet: نجمع القيم التي تعطي خارج قسمة واحدًا، ونحسب كل مجموعة مرة واحدة، ثم نعيد استعمال النتيجة.

مثال محلول: \(N=10\)

لإجراء فحص صغير، خذ \(N=10\). عندها

$$m_0=\left\lfloor\frac{10}{2}\right\rfloor=5,\qquad \lambda=\left\lfloor\frac{5}{\sqrt3}\right\rfloor=2,\qquad x_0=\left\lfloor\frac{10}{4}\right\rfloor=2.$$

كذلك

$$D(5)=\left\lfloor\frac51\right\rfloor+\left\lfloor\frac52\right\rfloor+\left\lfloor\frac53\right\rfloor+\left\lfloor\frac54\right\rfloor+\left\lfloor\frac55\right\rfloor=5+2+1+1+1=10.$$

الأزواج المرتبة التي تحقق \(uv\le 2\) هي \((1,1)\) و\((1,2)\) و\((2,1)\). وبما أن القاسم المشترك الأكبر في الحالات الثلاث يساوي \(1\)، فإن

$$A(2)=3\left\lfloor\frac{10}{2}\right\rfloor=15.$$

أما الحد الموزون فيساوي

$$B(2)=\lfloor\sqrt3\rfloor+\lfloor2\sqrt3\rfloor+\lfloor2\sqrt3\rfloor=1+3+3=7.$$

وأخيرًا، لأن \(Q(2,1)=7>2\)، فالمجموع السداسي الداخلي فارغ، وبالتالي

$$C(2)=D(2)=\left\lfloor\frac21\right\rfloor+\left\lfloor\frac22\right\rfloor=3.$$

ومن ثم

$$G(10)=2\cdot 10+4(15-7)-4\cdot 3=20+32-12=40.$$

هذا المثال الصغير يعكس البنية الدقيقة نفسها التي تستعملها الخوارزمية عند قيمة الإدخال الكاملة.

كيف تعمل الشيفرة

تبدأ تطبيقات ++C وPython وJava بحساب الحدود الثلاثة \(m_0\) و\(\lambda\) و\(x_0\). بعد ذلك تُجهَّز أصغر العوامل الأولية حتى \(\max(\lfloor\sqrt{\lambda}\rfloor,\lfloor\sqrt{x_0}\rfloor)\)، إلى جانب جدولين: جميع القواسم لكل عدد، والقواسم الخالية من المربعات مرفقة بإشارات Möbius.

في مساهمة الشريط، تدور الخوارزمية على الإحداثي الثاني \(v\)، ثم تقسم العمل بحسب كل \(d\mid v\)، وتستخدم عكس Möbius لعدّ الشركاء المختزلين \(a\) الذين يحققون \(\gcd(a,v/d)=1\). أما مجاميع الجزء الصحيح الموزونة بـ \(\sqrt3\) فتُحسب بواسطة تكرار Beatty الدقيق بدلًا من حلقة مباشرة.

في التصحيح السداسي، تُفحَص الأزواج بحسب \(v\)، ويُزاد \(u\) على شكل كتل قصوى تبقى فيها \(\left\lfloor x_0/Q(u,v)\right\rfloor\) ثابتة. كل استعلام مميز لدالة مجموع القواسم يُخزَّن في ذاكرة مؤقتة، ولذلك تُحسب القيم المكررة مرة واحدة فقط. ثم تُجمع كل المساهمات ويؤخذ الناتج modulo \(10^9+7\).

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

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

ليكن \(P=\max(\lfloor\sqrt{\lambda}\rfloor,\lfloor\sqrt{x_0}\rfloor)\). بناء جداول التحليل إلى عوامل يكلف \(O(P\log\log P)\) زمنًا و\(O(P)\) ذاكرة. كما أن مرحلتي الشريط والتصحيح السداسي تستفيدان من ضغط القواسم، وهضاب خارج القسمة، وقيم دالة مجموع القواسم المخزنة مؤقتًا، ولذلك لا تُفحَص جميع الأزواج المرشحة فحصًا صريحًا.

عمليًا يتصرف الزمن الكلي مثل \(O(\sqrt N\,\mathrm{polylog}\,N)\). وتأتي العوامل الإضافية أساسًا من تعداد القواسم ومن التقييم التكراري لمجاميع الجزء الصحيح. أما استهلاك الذاكرة فيبقى في حدود \(O(\sqrt N)\).

حواشٍ ومراجع

  1. صفحة المسألة: https://projecteuler.net/problem=780
  2. صيغة عكس Möbius: Wikipedia — Möbius inversion formula
  3. دالة مجموع القواسم: Wikipedia — Divisor summatory function
  4. متتالية Beatty: Wikipedia — Beatty sequence
  5. أعداد أيزنشتاين الصحيحة: Wikipedia — Eisenstein integer
  6. الشبكة السداسية: Wikipedia — Hexagonal lattice