The silo has circular base radius \(R\). Grain is poured from a point whose vertical projection onto the base is \(P=(x,0)\), where \(x\in[0,R]\). Because the heap settles at angle of repose \(\alpha\), the empty region between the flat roof level and the conical surface has volume \(V(x)\). We must find every offset \(x\) for which \(V(x)\) is a perfect square, and then sum those offsets.
Let the base disk be
$$D=\{(u,v)\in\mathbb{R}^2:u^2+v^2\le R^2\}.$$
For a base point \(Q=(u,v)\), the horizontal distance to the pouring point is
$$r=\sqrt{(u-x)^2+v^2}.$$
The grain surface rises linearly with slope \(\tan(\alpha)\), so the local empty-space height above \(Q\) is \(r\tan(\alpha)\). Therefore the wasted volume is the surface integral
$$V(x)=\tan(\alpha)\iint_D \sqrt{(u-x)^2+v^2}\,du\,dv.$$
This is the quantity evaluated by the implementations.
Use polar coordinates centered at \(P\):
$$u=x+\rho\cos\varphi,\qquad v=\rho\sin\varphi.$$
Along a fixed direction \(\varphi\), the ray leaves the disk when
$$\left(x+\rho\cos\varphi\right)^2+\left(\rho\sin\varphi\right)^2=R^2.$$
Solving this quadratic for the nonnegative root gives the radial intersection length
$$s(\varphi)=-x\cos\varphi+\sqrt{R^2-x^2\sin^2\varphi}.$$
The Jacobian contributes a factor \(\rho\), and the local height contributes another factor \(\rho\tan(\alpha)\). Hence
$$V(x)=\tan(\alpha)\int_0^{2\pi}\int_0^{s(\varphi)} \rho^2\,d\rho\,d\varphi =\frac{\tan(\alpha)}{3}\int_0^{2\pi}s(\varphi)^3\,d\varphi.$$
This is exactly the formula used in the C++, Python, and Java implementations.
At the center, \(x=0\), every ray has the same length \(s(\varphi)=R\). So
$$V(0)=\frac{\tan(\alpha)}{3}\int_0^{2\pi}R^3\,d\varphi=\frac{2\pi R^3}{3}\tan(\alpha).$$
At the wall, \(x=R\), the ray length becomes
$$s(\varphi)=-R\cos\varphi+R\lvert\cos\varphi\rvert,$$
which is zero on the outward half-plane and equals \(-2R\cos\varphi\) on the inward half-plane. Therefore
$$V(R)=\frac{\tan(\alpha)}{3}\int_{\pi/2}^{3\pi/2}\left(-2R\cos\varphi\right)^3\,d\varphi =\frac{32R^3}{9}\tan(\alpha).$$
For the actual parameters \(R=6\) and \(\alpha=40^\circ\), this gives
$$V(0)\approx 379.599730118848,\qquad V(R)\approx 644.428516744151.$$
So the only possible square volumes are
$$20^2,\,21^2,\,22^2,\,23^2,\,24^2,\,25^2.$$
This is the key simplification: instead of searching over a continuum of offsets, we only need to solve six scalar equations \(V(x)=k^2\).
The sample checkpoint used by the implementations is \(R=3\) and \(\alpha=30^\circ\). Then
$$V(0)=\frac{2\pi\cdot 3^3}{3}\tan(30^\circ)=\frac{18\pi}{\sqrt{3}}\approx 32.648388556,$$
$$V(R)=\frac{32\cdot 3^3}{9}\tan(30^\circ)=\frac{96}{\sqrt{3}}\approx 55.425625842.$$
The only square targets in this interval are \(36\) and \(49\). Solving \(V(x)=36\) and \(V(x)=49\) numerically yields
$$x\approx 1.114785284,\qquad x\approx 2.511167869,$$
which matches the checkpoint values verified by the C++ implementation.
No elementary antiderivative is used for the angular integral, so the implementation evaluates
$$\int_0^{2\pi}s(\varphi)^3\,d\varphi$$
with Simpson's rule using an even number \(N_q=4096\) of subintervals. If \(h=2\pi/N_q\) and \(f_i=f(ih)\), then
$$\int_0^{2\pi}f(\varphi)\,d\varphi\approx \frac{h}{3}\left(f_0+f_{N_q}+4\sum_{j=1}^{N_q/2}f_{2j-1}+2\sum_{j=1}^{N_q/2-1}f_{2j}\right).$$
After the endpoint values are known, the implementations enumerate every integer \(k\) with \(k^2\in[V(0),V(R)]\). For each target square they apply bisection on the interval \([0,R]\). The method relies on the smooth increase of \(V(x)\) from the center toward the wall, so each admissible square contributes one root. Using 120 bisection steps is far more than enough for nine correct decimal places.
The C++, Python, and Java implementations all follow the same pipeline. First they evaluate \(V(0)\) and \(V(R)\). Next they list all integer squares inside that volume range. For each target square they repeatedly bisect the offset interval and reevaluate the Simpson integral at the midpoint until the root is isolated. Finally they sum the resulting offsets and format the answer to nine decimal places. The C++ version also checks the sample values before solving the main instance.
Let \(K\) be the number of square targets, \(B\) the number of bisection iterations, and \(N_q\) the Simpson subinterval count. One evaluation of \(V(x)\) costs \(O(N_q)\), so the full computation costs \(O(KBN_q)\) time and \(O(1)\) extra memory. For the actual problem, \(K=6\), \(B=120\), and \(N_q=4096\), so the runtime is dominated by a modest number of accurate integral evaluations.
Die Basis des Silos ist eine Kreisscheibe mit Radius \(R\). Das Korn fällt aus einem Punkt, dessen lotrechte Projektion auf die Grundfläche \(P=(x,0)\) ist, wobei \(x\in[0,R]\). Wegen des Böschungswinkels \(\alpha\) entsteht zwischen der ebenen Dachhöhe und der kegelförmigen Oberfläche ein Leerraum mit Volumen \(V(x)\). Gesucht sind alle Offsets \(x\), für die \(V(x)\) ein perfektes Quadrat ist; anschließend wird die Summe dieser Offsets gebildet.
Die Grundfläche sei
$$D=\{(u,v)\in\mathbb{R}^2:u^2+v^2\le R^2\}.$$
Für einen Punkt \(Q=(u,v)\) der Basis ist der horizontale Abstand zum Einfüllpunkt
$$r=\sqrt{(u-x)^2+v^2}.$$
Die Körneroberfläche steigt mit Steigung \(\tan(\alpha)\), daher ist die lokale Leerhöhe gleich \(r\tan(\alpha)\). Also gilt für das verschwendete Volumen
$$V(x)=\tan(\alpha)\iint_D \sqrt{(u-x)^2+v^2}\,du\,dv.$$
Genau diese Größe wird in den Implementierungen ausgewertet.
Wir verwenden Polarkoordinaten mit Zentrum im Einfüllpunkt \(P\):
$$u=x+\rho\cos\varphi,\qquad v=\rho\sin\varphi.$$
Entlang einer festen Richtung \(\varphi\) verlässt der Strahl die Scheibe, sobald
$$\left(x+\rho\cos\varphi\right)^2+\left(\rho\sin\varphi\right)^2=R^2$$
erfüllt ist. Die nichtnegative Lösung lautet
$$s(\varphi)=-x\cos\varphi+\sqrt{R^2-x^2\sin^2\varphi}.$$
Aus dem Jacobi-Faktor \(\rho\) und der Höhe \(\rho\tan(\alpha)\) folgt dann
$$V(x)=\tan(\alpha)\int_0^{2\pi}\int_0^{s(\varphi)} \rho^2\,d\rho\,d\varphi =\frac{\tan(\alpha)}{3}\int_0^{2\pi}s(\varphi)^3\,d\varphi.$$
Das ist exakt die Formel, die in C++, Python und Java umgesetzt wird.
Im Zentrum \(x=0\) ist jede Strahlenlänge gleich \(R\). Daher
$$V(0)=\frac{\tan(\alpha)}{3}\int_0^{2\pi}R^3\,d\varphi=\frac{2\pi R^3}{3}\tan(\alpha).$$
An der Wand \(x=R\) gilt
$$s(\varphi)=-R\cos\varphi+R\lvert\cos\varphi\rvert,$$
also auf der äußeren Halbebene \(0\) und auf der inneren Halbebene \(-2R\cos\varphi\). Damit erhält man
$$V(R)=\frac{\tan(\alpha)}{3}\int_{\pi/2}^{3\pi/2}\left(-2R\cos\varphi\right)^3\,d\varphi =\frac{32R^3}{9}\tan(\alpha).$$
Für \(R=6\) und \(\alpha=40^\circ\) ergibt das
$$V(0)\approx 379.599730118848,\qquad V(R)\approx 644.428516744151.$$
Damit kommen nur die Quadrate
$$20^2,\,21^2,\,22^2,\,23^2,\,24^2,\,25^2$$
als Zielwerte in Frage. Das Problem reduziert sich also auf sechs Gleichungen \(V(x)=k^2\).
Die Implementierungen verwenden als Kontrollfall \(R=3\) und \(\alpha=30^\circ\). Dann gilt
$$V(0)=\frac{2\pi\cdot 3^3}{3}\tan(30^\circ)=\frac{18\pi}{\sqrt{3}}\approx 32.648388556,$$
$$V(R)=\frac{32\cdot 3^3}{9}\tan(30^\circ)=\frac{96}{\sqrt{3}}\approx 55.425625842.$$
In diesem Intervall liegen nur die Quadrate \(36\) und \(49\). Numerisch erhält man
$$x\approx 1.114785284,\qquad x\approx 2.511167869,$$
genau die Werte, die von der C++-Implementierung überprüft werden.
Für das Winkelintegral wird keine elementare Stammfunktion benutzt. Stattdessen wird
$$\int_0^{2\pi}s(\varphi)^3\,d\varphi$$
mit der Simpson-Regel und einer geraden Zahl \(N_q=4096\) von Teilintervallen angenähert. Für \(h=2\pi/N_q\) und \(f_i=f(ih)\) gilt
$$\int_0^{2\pi}f(\varphi)\,d\varphi\approx \frac{h}{3}\left(f_0+f_{N_q}+4\sum_{j=1}^{N_q/2}f_{2j-1}+2\sum_{j=1}^{N_q/2-1}f_{2j}\right).$$
Danach werden alle ganzen Zahlen \(k\) mit \(k^2\in[V(0),V(R)]\) betrachtet. Für jedes Quadrat wird auf \([0,R]\) eine Bisektion durchgeführt. Die Implementierungen nutzen dabei aus, dass \(V(x)\) vom Zentrum zur Wand glatt ansteigt; daher liefert jedes zulässige Quadrat genau eine Lösung. Mit 120 Bisektionsschritten ist die gewünschte Genauigkeit von neun Dezimalstellen problemlos erreichbar.
Die C++-, Python- und Java-Implementierungen folgen derselben Struktur. Zuerst werden \(V(0)\) und \(V(R)\) bestimmt. Danach werden alle Quadratwerte innerhalb dieses Bereichs aufgelistet. Für jedes Zielquadrat wird das zugehörige \(x\) per Bisektion gesucht, wobei an jedem Mittelpunkt das Simpson-Integral neu ausgewertet wird. Am Ende werden alle gefundenen Offsets addiert und mit neun Nachkommastellen formatiert. Die C++-Variante führt vor der Hauptrechnung zusätzlich die Kontrollwerte aus.
Seien \(K\) die Anzahl der Quadratziele, \(B\) die Zahl der Bisektionsschritte und \(N_q\) die Zahl der Simpson-Teilintervalle. Eine Auswertung von \(V(x)\) kostet \(O(N_q)\), somit liegt die Gesamtlaufzeit bei \(O(KBN_q)\), während der Zusatzspeicher \(O(1)\) bleibt. Im echten Problem gilt \(K=6\), \(B=120\) und \(N_q=4096\).
Silonun tabanı yarıçapı \(R\) olan bir dairedir. Tahıl, tabana izdüşümü \(P=(x,0)\) olan bir noktadan dökülür ve burada \(x\in[0,R]\) olur. Yığının eğimi dinlenme açısı \(\alpha\) ile belirlendiği için, düz tavan seviyesi ile konik yüzey arasında \(V(x)\) hacminde bir boşluk oluşur. İstenen şey, \(V(x)\) bir tam kare olduğunda ortaya çıkan tüm \(x\) değerlerini bulup bu değerleri toplamaktır.
Taban diski
$$D=\{(u,v)\in\mathbb{R}^2:u^2+v^2\le R^2\}$$
olsun. Taban üzerindeki \(Q=(u,v)\) noktası için döküm noktasına yatay uzaklık
$$r=\sqrt{(u-x)^2+v^2}$$
olur. Yüzey, \(\tan(\alpha)\) eğimiyle yükseldiğinden bu noktadaki boşluk yüksekliği \(r\tan(\alpha)\) olur. Dolayısıyla boşa giden hacim
$$V(x)=\tan(\alpha)\iint_D \sqrt{(u-x)^2+v^2}\,du\,dv$$
şeklindedir. Uygulamalar tam olarak bu integrali hesaplar.
Merkezi \(P\) olan kutupsal koordinatlara geçelim:
$$u=x+\rho\cos\varphi,\qquad v=\rho\sin\varphi.$$
Sabit bir \(\varphi\) yönünde ışın,
$$\left(x+\rho\cos\varphi\right)^2+\left(\rho\sin\varphi\right)^2=R^2$$
olduğunda diskten çıkar. Negatif olmayan kök
$$s(\varphi)=-x\cos\varphi+\sqrt{R^2-x^2\sin^2\varphi}$$
olur. Alan elemanındaki \(\rho\) çarpanı ve yerel yükseklikteki \(\rho\tan(\alpha)\) birlikte
$$V(x)=\tan(\alpha)\int_0^{2\pi}\int_0^{s(\varphi)} \rho^2\,d\rho\,d\varphi =\frac{\tan(\alpha)}{3}\int_0^{2\pi}s(\varphi)^3\,d\varphi$$
sonucunu verir. C++, Python ve Java çözümlerinin kullandığı temel formül budur.
Merkezde \(x=0\) iken bütün ışınların uzunluğu \(R\) olur. Bu yüzden
$$V(0)=\frac{\tan(\alpha)}{3}\int_0^{2\pi}R^3\,d\varphi=\frac{2\pi R^3}{3}\tan(\alpha).$$
Duvarda \(x=R\) iken
$$s(\varphi)=-R\cos\varphi+R\lvert\cos\varphi\rvert$$
elde edilir; dış yarı düzlemde sıfır, iç yarı düzlemde ise \(-2R\cos\varphi\) olur. Dolayısıyla
$$V(R)=\frac{\tan(\alpha)}{3}\int_{\pi/2}^{3\pi/2}\left(-2R\cos\varphi\right)^3\,d\varphi =\frac{32R^3}{9}\tan(\alpha).$$
Gerçek parametreler \(R=6\) ve \(\alpha=40^\circ\) için
$$V(0)\approx 379.599730118848,\qquad V(R)\approx 644.428516744151$$
bulunur. Bu aralıkta yalnızca
$$20^2,\,21^2,\,22^2,\,23^2,\,24^2,\,25^2$$
kareleri vardır. Yani sürekli bir \(x\) aralığını taramak yerine sadece altı adet \(V(x)=k^2\) denklemi çözülür.
Uygulamalardaki örnek kontrol durumu \(R=3\) ve \(\alpha=30^\circ\) değerleridir. Bu durumda
$$V(0)=\frac{2\pi\cdot 3^3}{3}\tan(30^\circ)=\frac{18\pi}{\sqrt{3}}\approx 32.648388556,$$
$$V(R)=\frac{32\cdot 3^3}{9}\tan(30^\circ)=\frac{96}{\sqrt{3}}\approx 55.425625842.$$
Aralık içindeki tek kareler \(36\) ve \(49\) olduğundan, sayısal çözüm
$$x\approx 1.114785284,\qquad x\approx 2.511167869$$
değerlerini verir. Bunlar C++ uygulamasının doğruladığı kontrol sonuçlarıyla aynıdır.
Açısal integral için kapalı bir antitürev kullanılmaz. Bunun yerine
$$\int_0^{2\pi}s(\varphi)^3\,d\varphi$$
ifadesi Simpson kuralı ile, çift sayıda \(N_q=4096\) alt aralık kullanılarak hesaplanır. \(h=2\pi/N_q\) ve \(f_i=f(ih)\) için
$$\int_0^{2\pi}f(\varphi)\,d\varphi\approx \frac{h}{3}\left(f_0+f_{N_q}+4\sum_{j=1}^{N_q/2}f_{2j-1}+2\sum_{j=1}^{N_q/2-1}f_{2j}\right).$$
Daha sonra \(k^2\in[V(0),V(R)]\) koşulunu sağlayan her tamsayı \(k\) için \([0,R]\) aralığında ikili arama uygulanır. Çözümler, \(V(x)\)'in merkezden duvara doğru düzgün biçimde arttığı gözlemine dayanır; bu da her uygun kare için tek bir kök verilmesini sağlar. 120 ikili arama adımı dokuz ondalık basamak için fazlasıyla yeterlidir.
C++, Python ve Java uygulamaları aynı işlem hattını izler. Önce \(V(0)\) ve \(V(R)\) hesaplanır. Sonra bu aralıkta kalan tüm kare hacimler listelenir. Her hedef kare için orta noktada Simpson integrali yeniden hesaplanarak ikili arama yapılır ve uygun \(x\) değeri elde edilir. Son aşamada bu \(x\) değerleri toplanır ve sonuç dokuz ondalık basamakla yazdırılır. C++ sürümü ana problemi çözmeden önce örnek kontrol değerlerini de sınar.
\(K\) kare hedef sayısı, \(B\) ikili arama iterasyonu, \(N_q\) ise Simpson alt aralık sayısı olsun. Tek bir \(V(x)\) hesabı \(O(N_q)\) zaman alır. Bu nedenle toplam karmaşıklık \(O(KBN_q)\), ek bellek ise \(O(1)\) olur. Gerçek durumda \(K=6\), \(B=120\) ve \(N_q=4096\) kullanılır.
La base del silo es un disco de radio \(R\). El grano se vierte desde un punto cuya proyección vertical sobre la base es \(P=(x,0)\), con \(x\in[0,R]\). Debido al ángulo de reposo \(\alpha\), entre el techo plano y la superficie cónica del montón queda un volumen vacío \(V(x)\). Hay que encontrar todos los desplazamientos \(x\) para los cuales \(V(x)\) es un cuadrado perfecto y sumar esos valores.
Sea
$$D=\{(u,v)\in\mathbb{R}^2:u^2+v^2\le R^2\}$$
el disco de la base. Para un punto \(Q=(u,v)\), la distancia horizontal al punto de descarga es
$$r=\sqrt{(u-x)^2+v^2}.$$
La superficie del montón sube con pendiente \(\tan(\alpha)\), así que la altura vacía local es \(r\tan(\alpha)\). Por tanto, el volumen buscado es
$$V(x)=\tan(\alpha)\iint_D \sqrt{(u-x)^2+v^2}\,du\,dv.$$
Esa es exactamente la magnitud que evalúan las implementaciones.
Usamos coordenadas polares centradas en \(P\):
$$u=x+\rho\cos\varphi,\qquad v=\rho\sin\varphi.$$
En una dirección fija \(\varphi\), el rayo sale del disco cuando
$$\left(x+\rho\cos\varphi\right)^2+\left(\rho\sin\varphi\right)^2=R^2.$$
La raíz no negativa es
$$s(\varphi)=-x\cos\varphi+\sqrt{R^2-x^2\sin^2\varphi}.$$
El jacobiano aporta un factor \(\rho\), y la altura aporta otro \(\rho\tan(\alpha)\). Entonces
$$V(x)=\tan(\alpha)\int_0^{2\pi}\int_0^{s(\varphi)} \rho^2\,d\rho\,d\varphi =\frac{\tan(\alpha)}{3}\int_0^{2\pi}s(\varphi)^3\,d\varphi.$$
Ésta es la fórmula implementada en C++, Python y Java.
En el centro, \(x=0\), toda dirección tiene longitud \(R\). Por eso
$$V(0)=\frac{\tan(\alpha)}{3}\int_0^{2\pi}R^3\,d\varphi=\frac{2\pi R^3}{3}\tan(\alpha).$$
En la pared, \(x=R\), se obtiene
$$s(\varphi)=-R\cos\varphi+R\lvert\cos\varphi\rvert,$$
que vale \(0\) hacia fuera y \(-2R\cos\varphi\) hacia dentro. Así,
$$V(R)=\frac{\tan(\alpha)}{3}\int_{\pi/2}^{3\pi/2}\left(-2R\cos\varphi\right)^3\,d\varphi =\frac{32R^3}{9}\tan(\alpha).$$
Para \(R=6\) y \(\alpha=40^\circ\), esto da
$$V(0)\approx 379.599730118848,\qquad V(R)\approx 644.428516744151.$$
Por lo tanto, los únicos cuadrados posibles son
$$20^2,\,21^2,\,22^2,\,23^2,\,24^2,\,25^2.$$
La búsqueda continua en \(x\) se transforma así en sólo seis ecuaciones del tipo \(V(x)=k^2\).
La comprobación numérica usada por las implementaciones toma \(R=3\) y \(\alpha=30^\circ\). Entonces
$$V(0)=\frac{2\pi\cdot 3^3}{3}\tan(30^\circ)=\frac{18\pi}{\sqrt{3}}\approx 32.648388556,$$
$$V(R)=\frac{32\cdot 3^3}{9}\tan(30^\circ)=\frac{96}{\sqrt{3}}\approx 55.425625842.$$
Los únicos cuadrados dentro de ese intervalo son \(36\) y \(49\). Al resolver numéricamente \(V(x)=36\) y \(V(x)=49\), se obtiene
$$x\approx 1.114785284,\qquad x\approx 2.511167869,$$
que coincide con los valores de control verificados por la implementación en C++.
No se usa una primitiva cerrada para la integral angular. En su lugar,
$$\int_0^{2\pi}s(\varphi)^3\,d\varphi$$
se aproxima con la regla de Simpson usando un número par \(N_q=4096\) de subintervalos. Si \(h=2\pi/N_q\) y \(f_i=f(ih)\), entonces
$$\int_0^{2\pi}f(\varphi)\,d\varphi\approx \frac{h}{3}\left(f_0+f_{N_q}+4\sum_{j=1}^{N_q/2}f_{2j-1}+2\sum_{j=1}^{N_q/2-1}f_{2j}\right).$$
Una vez conocidos los extremos, las implementaciones enumeran todos los enteros \(k\) con \(k^2\in[V(0),V(R)]\). Para cada cuadrado, aplican bisección en \([0,R]\). El método se apoya en que \(V(x)\) crece suavemente desde el centro hasta la pared, por lo que cada cuadrado admisible aporta una única solución. Con 120 iteraciones de bisección sobra precisión para nueve decimales.
Las implementaciones en C++, Python y Java siguen la misma estrategia. Primero calculan \(V(0)\) y \(V(R)\). Después enumeran los cuadrados enteros dentro de ese intervalo. Para cada objetivo vuelven a evaluar la integral de Simpson en los puntos medios de la bisección hasta aislar el desplazamiento correcto. Finalmente suman todos los valores de \(x\) encontrados y formatean el resultado con nueve decimales. La versión en C++ además comprueba el caso de ejemplo antes de resolver la instancia principal.
Sea \(K\) el número de cuadrados candidatos, \(B\) el número de iteraciones de bisección y \(N_q\) el número de subintervalos de Simpson. Una evaluación de \(V(x)\) cuesta \(O(N_q)\), así que el coste total es \(O(KBN_q)\) en tiempo y \(O(1)\) en memoria adicional. En la instancia real, \(K=6\), \(B=120\) y \(N_q=4096\).
粮仓的底面是半径为 \(R\) 的圆盘。投料点在底面上的垂直投影记为 \(P=(x,0)\),其中 \(x\in[0,R]\)。由于谷物会以安息角 \(\alpha\) 形成圆锥形堆面,平顶高度与该堆面之间会留下一个空体积 \(V(x)\)。题目要求找出所有使 \(V(x)\) 成为完全平方数的偏移量 \(x\),然后求这些 \(x\) 的总和。
设底面圆盘为
$$D=\{(u,v)\in\mathbb{R}^2:u^2+v^2\le R^2\}.$$
对底面上一点 \(Q=(u,v)\),它到投料点的水平距离为
$$r=\sqrt{(u-x)^2+v^2}.$$
堆面按斜率 \(\tan(\alpha)\) 上升,所以该点处的空隙高度就是 \(r\tan(\alpha)\)。因此浪费体积满足
$$V(x)=\tan(\alpha)\iint_D \sqrt{(u-x)^2+v^2}\,du\,dv.$$
这正是实现所计算的核心量。
以 \(P\) 为原点建立极坐标:
$$u=x+\rho\cos\varphi,\qquad v=\rho\sin\varphi.$$
在固定方向 \(\varphi\) 上,射线离开圆盘边界时满足
$$\left(x+\rho\cos\varphi\right)^2+\left(\rho\sin\varphi\right)^2=R^2.$$
取非负根,可得该方向上的截距长度
$$s(\varphi)=-x\cos\varphi+\sqrt{R^2-x^2\sin^2\varphi}.$$
面积元带来一个 \(\rho\),高度又带来一个 \(\rho\tan(\alpha)\),于是
$$V(x)=\tan(\alpha)\int_0^{2\pi}\int_0^{s(\varphi)} \rho^2\,d\rho\,d\varphi =\frac{\tan(\alpha)}{3}\int_0^{2\pi}s(\varphi)^3\,d\varphi.$$
这与 C++、Python 和 Java 实现中的公式完全一致。
当 \(x=0\) 时,所有方向上的长度都等于 \(R\),所以
$$V(0)=\frac{\tan(\alpha)}{3}\int_0^{2\pi}R^3\,d\varphi=\frac{2\pi R^3}{3}\tan(\alpha).$$
当 \(x=R\) 时,
$$s(\varphi)=-R\cos\varphi+R\lvert\cos\varphi\rvert,$$
朝外半平面时为 \(0\),朝内半平面时为 \(-2R\cos\varphi\)。因此
$$V(R)=\frac{\tan(\alpha)}{3}\int_{\pi/2}^{3\pi/2}\left(-2R\cos\varphi\right)^3\,d\varphi =\frac{32R^3}{9}\tan(\alpha).$$
对实际参数 \(R=6\)、\(\alpha=40^\circ\),有
$$V(0)\approx 379.599730118848,\qquad V(R)\approx 644.428516744151.$$
所以可能的平方体积只有
$$20^2,\,21^2,\,22^2,\,23^2,\,24^2,\,25^2.$$
这一步非常重要,因为它把连续区间上的搜索压缩成了 6 个标量方程 \(V(x)=k^2\)。
实现中使用的校验样例是 \(R=3\)、\(\alpha=30^\circ\)。这时
$$V(0)=\frac{2\pi\cdot 3^3}{3}\tan(30^\circ)=\frac{18\pi}{\sqrt{3}}\approx 32.648388556,$$
$$V(R)=\frac{32\cdot 3^3}{9}\tan(30^\circ)=\frac{96}{\sqrt{3}}\approx 55.425625842.$$
该区间内只有两个平方数 \(36\) 和 \(49\)。数值求解 \(V(x)=36\) 与 \(V(x)=49\) 可得到
$$x\approx 1.114785284,\qquad x\approx 2.511167869,$$
与 C++ 实现内置的校验值一致。
角度积分没有使用初等闭式原函数,而是直接对
$$\int_0^{2\pi}s(\varphi)^3\,d\varphi$$
应用 Simpson 求积,并取偶数个子区间 \(N_q=4096\)。若 \(h=2\pi/N_q\),且 \(f_i=f(ih)\),则
$$\int_0^{2\pi}f(\varphi)\,d\varphi\approx \frac{h}{3}\left(f_0+f_{N_q}+4\sum_{j=1}^{N_q/2}f_{2j-1}+2\sum_{j=1}^{N_q/2-1}f_{2j}\right).$$
得到端点体积后,实现会枚举所有满足 \(k^2\in[V(0),V(R)]\) 的整数 \(k\),并在区间 \([0,R]\) 上做二分。算法利用了 \(V(x)\) 从中心向边界平滑增大的性质,因此每个允许的平方值对应一个根。使用 120 次二分迭代,足以稳定得到 9 位小数。
C++、Python 和 Java 实现遵循同一流程。先计算 \(V(0)\) 与 \(V(R)\),再枚举其中所有整数平方体积。随后对每个目标平方,在二分中不断重算中点处的 Simpson 积分,直到对应的偏移量被定位。最后把所有求得的 \(x\) 相加,并按 9 位小数输出。C++ 版本还会先验证样例检查点。
设平方目标个数为 \(K\),二分迭代次数为 \(B\),Simpson 子区间数为 \(N_q\)。单次 \(V(x)\) 计算代价是 \(O(N_q)\),因此总时间复杂度为 \(O(KBN_q)\),额外空间复杂度为 \(O(1)\)。在实际参数下,\(K=6\)、\(B=120\)、\(N_q=4096\)。
Основание силоса представляет собой круг радиуса \(R\). Зерно засыпается из точки, чья вертикальная проекция на основание равна \(P=(x,0)\), где \(x\in[0,R]\). Из-за угла естественного откоса \(\alpha\) между плоским верхним уровнем и конической поверхностью насыпи остаётся пустой объём \(V(x)\). Нужно найти все смещения \(x\), при которых \(V(x)\) является точным квадратом, а затем просуммировать эти значения.
Обозначим диск основания через
$$D=\{(u,v)\in\mathbb{R}^2:u^2+v^2\le R^2\}.$$
Для точки основания \(Q=(u,v)\) горизонтальное расстояние до точки засыпки равно
$$r=\sqrt{(u-x)^2+v^2}.$$
Поверхность насыпи поднимается с уклоном \(\tan(\alpha)\), поэтому локальная высота пустоты равна \(r\tan(\alpha)\). Следовательно, искомый объём задаётся интегралом
$$V(x)=\tan(\alpha)\iint_D \sqrt{(u-x)^2+v^2}\,du\,dv.$$
Именно эту величину вычисляют реализации.
Перейдём к полярным координатам с центром в точке \(P\):
$$u=x+\rho\cos\varphi,\qquad v=\rho\sin\varphi.$$
При фиксированном направлении \(\varphi\) луч выходит из диска, когда выполняется
$$\left(x+\rho\cos\varphi\right)^2+\left(\rho\sin\varphi\right)^2=R^2.$$
Ненулевая неотрицательная корневая ветвь даёт длину луча внутри диска:
$$s(\varphi)=-x\cos\varphi+\sqrt{R^2-x^2\sin^2\varphi}.$$
Якобиан даёт множитель \(\rho\), а высота добавляет ещё \(\rho\tan(\alpha)\). Поэтому
$$V(x)=\tan(\alpha)\int_0^{2\pi}\int_0^{s(\varphi)} \rho^2\,d\rho\,d\varphi =\frac{\tan(\alpha)}{3}\int_0^{2\pi}s(\varphi)^3\,d\varphi.$$
Это в точности формула, используемая в версиях на C++, Python и Java.
В центре, при \(x=0\), длина каждого луча равна \(R\). Поэтому
$$V(0)=\frac{\tan(\alpha)}{3}\int_0^{2\pi}R^3\,d\varphi=\frac{2\pi R^3}{3}\tan(\alpha).$$
У стенки, при \(x=R\), получаем
$$s(\varphi)=-R\cos\varphi+R\lvert\cos\varphi\rvert,$$
то есть \(0\) на внешней полуплоскости и \(-2R\cos\varphi\) на внутренней. Следовательно,
$$V(R)=\frac{\tan(\alpha)}{3}\int_{\pi/2}^{3\pi/2}\left(-2R\cos\varphi\right)^3\,d\varphi =\frac{32R^3}{9}\tan(\alpha).$$
Для фактических параметров \(R=6\) и \(\alpha=40^\circ\) это даёт
$$V(0)\approx 379.599730118848,\qquad V(R)\approx 644.428516744151.$$
Значит, возможными квадратами являются только
$$20^2,\,21^2,\,22^2,\,23^2,\,24^2,\,25^2.$$
Тем самым непрерывный поиск по \(x\) заменяется решением всего шести уравнений \(V(x)=k^2\).
В качестве контрольного случая реализации используют \(R=3\) и \(\alpha=30^\circ\). Тогда
$$V(0)=\frac{2\pi\cdot 3^3}{3}\tan(30^\circ)=\frac{18\pi}{\sqrt{3}}\approx 32.648388556,$$
$$V(R)=\frac{32\cdot 3^3}{9}\tan(30^\circ)=\frac{96}{\sqrt{3}}\approx 55.425625842.$$
В этом интервале лежат только квадраты \(36\) и \(49\). Численное решение уравнений \(V(x)=36\) и \(V(x)=49\) даёт
$$x\approx 1.114785284,\qquad x\approx 2.511167869,$$
что совпадает с контрольными значениями, проверяемыми в реализации на C++.
Для углового интеграла не используется элементарная первообразная. Вместо этого интеграл
$$\int_0^{2\pi}s(\varphi)^3\,d\varphi$$
аппроксимируется формулой Симпсона с чётным числом подотрезков \(N_q=4096\). При \(h=2\pi/N_q\) и \(f_i=f(ih)\) имеем
$$\int_0^{2\pi}f(\varphi)\,d\varphi\approx \frac{h}{3}\left(f_0+f_{N_q}+4\sum_{j=1}^{N_q/2}f_{2j-1}+2\sum_{j=1}^{N_q/2-1}f_{2j}\right).$$
После вычисления крайних объёмов реализации перебирают все целые \(k\), для которых \(k^2\in[V(0),V(R)]\), и для каждого квадрата применяют бисекцию на \([0,R]\). Алгоритм опирается на гладкий рост \(V(x)\) при движении от центра к стенке, поэтому каждому допустимому квадрату соответствует единственный корень. 120 шагов бисекции более чем достаточно для точности в 9 знаков после запятой.
Реализации на C++, Python и Java используют одну и ту же схему. Сначала вычисляются \(V(0)\) и \(V(R)\). Затем перечисляются все квадратные объёмы внутри этого диапазона. Для каждого целевого квадрата выполняется бисекция, а в середине каждого интервала заново вычисляется интеграл по формуле Симпсона. После этого найденные значения \(x\) суммируются, и результат выводится с девятью знаками после запятой. Версия на C++ дополнительно проверяет контрольный пример перед основной задачей.
Пусть \(K\) — число квадратных целей, \(B\) — количество шагов бисекции, а \(N_q\) — число подотрезков в правиле Симпсона. Одна оценка \(V(x)\) стоит \(O(N_q)\), поэтому полная трудоёмкость равна \(O(KBN_q)\), а дополнительная память — \(O(1)\). В реальной задаче используются \(K=6\), \(B=120\) и \(N_q=4096\).
قاعدة الصومعة قرص دائري نصف قطره \(R\). تُسكب الحبوب من نقطة إسقاطها العمودي على القاعدة هو \(P=(x,0)\)، حيث \(x\in[0,R]\). بسبب زاوية الاستقرار \(\alpha\)، يتكوّن بين المستوى العلوي المستوي والسطح المخروطي للكومة حجم فارغ مقداره \(V(x)\). المطلوب إيجاد كل قيم \(x\) التي تجعل \(V(x)\) مربعاً كاملاً ثم جمع هذه القيم.
لنرمز لقرص القاعدة بـ
$$D=\{(u,v)\in\mathbb{R}^2:u^2+v^2\le R^2\}.$$
ولأي نقطة \(Q=(u,v)\) على القاعدة، تكون المسافة الأفقية إلى نقطة السكب
$$r=\sqrt{(u-x)^2+v^2}.$$
يرتفع سطح الحبوب بميل مقداره \(\tan(\alpha)\)، لذا فإن ارتفاع الفراغ فوق هذه النقطة يساوي \(r\tan(\alpha)\). ومن ثم فإن الحجم المهدور يساوي
$$V(x)=\tan(\alpha)\iint_D \sqrt{(u-x)^2+v^2}\,du\,dv.$$
وهذه هي الكمية التي تحسبها جميع التطبيقات مباشرة.
نستخدم إحداثيات قطبية مركزها عند \(P\):
$$u=x+\rho\cos\varphi,\qquad v=\rho\sin\varphi.$$
وفي اتجاه ثابت \(\varphi\)، يخرج الشعاع من القرص عندما يتحقق
$$\left(x+\rho\cos\varphi\right)^2+\left(\rho\sin\varphi\right)^2=R^2.$$
الجذر غير السالب يعطي طول التقاطع داخل القرص:
$$s(\varphi)=-x\cos\varphi+\sqrt{R^2-x^2\sin^2\varphi}.$$
يعطي يعقوبي التحويل عاملاً \(\rho\)، ويضيف الارتفاع عاملاً آخر \(\rho\tan(\alpha)\). لذلك نحصل على
$$V(x)=\tan(\alpha)\int_0^{2\pi}\int_0^{s(\varphi)} \rho^2\,d\rho\,d\varphi =\frac{\tan(\alpha)}{3}\int_0^{2\pi}s(\varphi)^3\,d\varphi.$$
وهذه هي الصيغة نفسها المستخدمة في تطبيقات C++ وPython وJava.
عند المركز \(x=0\)، يكون طول كل شعاع مساوياً لـ \(R\)، ومن ثم
$$V(0)=\frac{\tan(\alpha)}{3}\int_0^{2\pi}R^3\,d\varphi=\frac{2\pi R^3}{3}\tan(\alpha).$$
وعند الجدار \(x=R\)، يصبح
$$s(\varphi)=-R\cos\varphi+R\lvert\cos\varphi\rvert,$$
أي إنه يساوي \(0\) في النصف الخارجي و\(-2R\cos\varphi\) في النصف الداخلي. لذلك
$$V(R)=\frac{\tan(\alpha)}{3}\int_{\pi/2}^{3\pi/2}\left(-2R\cos\varphi\right)^3\,d\varphi =\frac{32R^3}{9}\tan(\alpha).$$
وللقيم الفعلية \(R=6\) و\(\alpha=40^\circ\) نحصل على
$$V(0)\approx 379.599730118848,\qquad V(R)\approx 644.428516744151.$$
إذن المربعات الوحيدة داخل هذا المجال هي
$$20^2,\,21^2,\,22^2,\,23^2,\,24^2,\,25^2.$$
وبذلك لا نبحث في مجال مستمر من قيم \(x\)، بل نحل ست معادلات فقط من الشكل \(V(x)=k^2\).
حالة التحقق التي تستخدمها التطبيقات هي \(R=3\) و\(\alpha=30^\circ\). عندها
$$V(0)=\frac{2\pi\cdot 3^3}{3}\tan(30^\circ)=\frac{18\pi}{\sqrt{3}}\approx 32.648388556,$$
$$V(R)=\frac{32\cdot 3^3}{9}\tan(30^\circ)=\frac{96}{\sqrt{3}}\approx 55.425625842.$$
المربعان الوحيدان داخل هذا المجال هما \(36\) و\(49\). وبحل المعادلتين \(V(x)=36\) و\(V(x)=49\) عددياً نحصل على
$$x\approx 1.114785284,\qquad x\approx 2.511167869,$$
وهما القيمتان اللتان تتحقق منهما نسخة C++.
لا تُستخدم صيغة ابتدائية مغلقة للتكامل الزاوي. بدلاً من ذلك يُقرَّب التكامل
$$\int_0^{2\pi}s(\varphi)^3\,d\varphi$$
بقاعدة سمبسون مع عدد زوجي من المقاطع \(N_q=4096\). إذا كان \(h=2\pi/N_q\) و\(f_i=f(ih)\)، فإن
$$\int_0^{2\pi}f(\varphi)\,d\varphi\approx \frac{h}{3}\left(f_0+f_{N_q}+4\sum_{j=1}^{N_q/2}f_{2j-1}+2\sum_{j=1}^{N_q/2-1}f_{2j}\right).$$
بعد حساب القيمتين الطرفيتين، تعدّد التطبيقات كل عدد صحيح \(k\) يحقق \(k^2\in[V(0),V(R)]\)، ثم تطبق البحث الثنائي على الفترة \([0,R]\) لكل مربع هدف. وتعتمد الطريقة على أن \(V(x)\) يزداد بسلاسة عند الانتقال من المركز إلى الجدار، لذلك يعطي كل مربع مسموح به جذراً واحداً. واستخدام 120 خطوة من البحث الثنائي يكفي بسهولة للحصول على 9 منازل عشرية صحيحة.
تتبع تطبيقات C++ وPython وJava المسار نفسه. أولاً تُحسب \(V(0)\) و\(V(R)\). ثم تُسرد جميع المربعات الصحيحة الواقعة داخل هذا المجال. بعد ذلك يُعاد حساب تكامل سمبسون عند نقاط المنتصف أثناء البحث الثنائي حتى يُعزل كل إزاحة مطلوبة. وفي النهاية تُجمع قيم \(x\) الناتجة ويُنسّق الجواب حتى 9 منازل عشرية. كما أن نسخة C++ تتحقق من حالة المثال قبل حل المسألة الأساسية.
إذا كان \(K\) عدد المربعات الهدف، و\(B\) عدد خطوات البحث الثنائي، و\(N_q\) عدد مقاطع سمبسون، فإن حساباً واحداً لـ \(V(x)\) يكلف \(O(N_q)\). لذا تكون الكلفة الكلية \(O(KBN_q)\) والزمن الإضافي في الذاكرة \(O(1)\). وفي الحالة الفعلية لدينا \(K=6\)، و\(B=120\)، و\(N_q=4096\).