The candy core is the ellipsoid of revolution
$$\frac{x^2+y^2}{a^2}+\frac{z^2}{b^2}=1,$$
equivalently \(b^2x^2+b^2y^2+a^2z^2=a^2b^2\). A chocolate layer of thickness \(1\) mm is added along the outward surface normal, so the outer boundary is the parallel surface at distance \(1\). The required quantity \(C(a,b)\) is the volume of that outer body minus the volume of the original ellipsoid.
The important geometric point is that a normal offset is not obtained by simply replacing \((a,b)\) with \((a+1,b+1)\). We must derive the offset surface explicitly.
Because the body is rotationally symmetric around the \(z\)-axis, it is enough to work in the meridian plane \((\rho,z)\), where \(\rho=\sqrt{x^2+y^2}\). The inner ellipse is
$$\frac{\rho^2}{a^2}+\frac{z^2}{b^2}=1.$$
A convenient parametrization of its upper half is
$$\rho(\theta)=a\sin\theta,\qquad z(\theta)=b\cos\theta,\qquad 0\le\theta\le\frac{\pi}{2}.$$
Rotating this arc around the \(z\)-axis recovers the whole ellipsoid.
Write the ellipse implicitly as
$$F(\rho,z)=\frac{\rho^2}{a^2}+\frac{z^2}{b^2}-1=0.$$
An outward normal is proportional to \(\nabla F\), so on the parametrized curve we get the normal direction
$$\left(\frac{\rho}{a^2},\frac{z}{b^2}\right)=\left(\frac{\sin\theta}{a},\frac{\cos\theta}{b}\right).$$
Define
$$E(\theta)=\sqrt{\frac{\sin^2\theta}{a^2}+\frac{\cos^2\theta}{b^2}}.$$
Then the outward unit normal in the meridian plane is
$$\mathbf{n}(\theta)=\frac{1}{E(\theta)}\left(\frac{\sin\theta}{a},\frac{\cos\theta}{b}\right).$$
Moving one unit along the outward normal sends the meridian point \((\rho(\theta),z(\theta))\) to the outer profile \((R(\theta),Z(\theta))\):
$$R(\theta)=a\sin\theta+\frac{\sin\theta}{aE(\theta)}=\sin\theta\left(a+\frac{1}{aE(\theta)}\right),$$
$$Z(\theta)=b\cos\theta+\frac{\cos\theta}{bE(\theta)}=\cos\theta\left(b+\frac{1}{bE(\theta)}\right).$$
This is exactly the geometry evaluated by the implementations.
For the volume integral we need \(dZ\). First differentiate \(E(\theta)\):
$$E'(\theta)=\frac{\sin\theta\cos\theta}{E(\theta)}\left(\frac{1}{a^2}-\frac{1}{b^2}\right).$$
Now differentiate \(Z(\theta)=b\cos\theta+\frac{\cos\theta}{bE(\theta)}\). After simplification,
$$-Z'(\theta)=\sin\theta\left(b+\frac{1}{bE(\theta)}+\frac{\cos^2\theta}{b}\frac{\frac{1}{a^2}-\frac{1}{b^2}}{E(\theta)^3}\right).$$
The integrand used numerically is therefore \(R(\theta)^2(-Z'(\theta))\).
The outer surface is symmetric with respect to the plane \(z=0\), so we integrate the upper half and double. Using the disk method for a solid of revolution,
$$V_{\text{outer}}=2\pi\int_0^{\pi/2} R(\theta)^2\bigl(-Z'(\theta)\bigr)\,d\theta.$$
The inner ellipsoid volume is the standard formula
$$V_{\text{inner}}=\frac{4\pi a^2 b}{3}.$$
Hence the chocolate volume is
$$\boxed{C(a,b)=2\pi\int_0^{\pi/2} R(\theta)^2\bigl(-Z'(\theta)\bigr)\,d\theta-\frac{4\pi a^2 b}{3}.}$$
When \(a=b=1\), the core is the unit sphere. Then \(E(\theta)=1\), so
$$R(\theta)=2\sin\theta,\qquad Z(\theta)=2\cos\theta,$$
which means the outer body is a sphere of radius \(2\). Therefore
$$C(1,1)=\frac{4\pi}{3}(2^3-1^3)=\frac{28\pi}{3}.$$
For the second checkpoint, numerical evaluation gives
$$C(2,1)\approx 60.35475635.$$
Applying the same integral to the target parameters yields
$$C(3,1)\approx 103.37870096,$$
which is the value rounded to eight decimal places by the implementations.
The C++, Python, and Java implementations all follow the same pipeline. They evaluate the smooth meridian integrand derived above, apply Simpson's rule on \([0,\pi/2]\), and then refine the interval recursively with adaptive Simpson integration until the local error estimate is below the requested tolerance or the recursion depth limit is reached.
Once the integral is stable, the implementation multiplies by \(2\pi\), subtracts the ellipsoid volume \(\frac{4\pi a^2b}{3}\), and formats the result to eight digits after the decimal point. The two published checkpoints are used as numerical sanity checks for the geometry and the quadrature.
Adaptive Simpson integration is data-dependent. If the interval is split into \(m\) accepted subintervals, the total work is \(O(m)\) function evaluations up to a constant factor, and the extra memory is \(O(d)\), where \(d\) is the recursion depth. For this problem the integrand is smooth on \([0,\pi/2]\), so convergence is fast and memory usage is negligible.
Der Kern der Süßigkeit ist das Rotationsellipsoid
$$\frac{x^2+y^2}{a^2}+\frac{z^2}{b^2}=1,$$
äquivalent zu \(b^2x^2+b^2y^2+a^2z^2=a^2b^2\). Darauf wird eine Schokoladenschicht der Dicke \(1\) mm aufgetragen, gemessen entlang der äußeren Flächennormale. Gesucht ist \(C(a,b)\), also das Volumen des äußeren Körpers minus das Volumen des inneren Ellipsoids.
Wesentlich ist, dass ein Normalenabstand nicht einfach durch Ersetzen von \((a,b)\) mit \((a+1,b+1)\) entsteht. Die äußere Fläche muss explizit als Abstandsoffset hergeleitet werden.
Wegen der Rotationssymmetrie um die \(z\)-Achse genügt die Meridianebene \((\rho,z)\) mit \(\rho=\sqrt{x^2+y^2}\). Die innere Ellipse erfüllt
$$\frac{\rho^2}{a^2}+\frac{z^2}{b^2}=1.$$
Eine praktische Parametrisierung der oberen Hälfte ist
$$\rho(\theta)=a\sin\theta,\qquad z(\theta)=b\cos\theta,\qquad 0\le\theta\le\frac{\pi}{2}.$$
Durch Rotation dieses Bogens um die \(z\)-Achse erhält man das gesamte Ellipsoid.
Schreibe die Ellipse implizit als
$$F(\rho,z)=\frac{\rho^2}{a^2}+\frac{z^2}{b^2}-1=0.$$
Eine äußere Normale ist proportional zu \(\nabla F\). Auf der Parametrisierung ergibt sich daher die Richtungsnormale
$$\left(\frac{\rho}{a^2},\frac{z}{b^2}\right)=\left(\frac{\sin\theta}{a},\frac{\cos\theta}{b}\right).$$
Definiere
$$E(\theta)=\sqrt{\frac{\sin^2\theta}{a^2}+\frac{\cos^2\theta}{b^2}}.$$
Dann lautet die äußere Einheitsnormale in der Meridianebene
$$\mathbf{n}(\theta)=\frac{1}{E(\theta)}\left(\frac{\sin\theta}{a},\frac{\cos\theta}{b}\right).$$
Wenn man den Punkt \((\rho(\theta),z(\theta))\) um eine Einheit entlang der äußeren Normalen verschiebt, erhält man das äußere Profil \((R(\theta),Z(\theta))\):
$$R(\theta)=a\sin\theta+\frac{\sin\theta}{aE(\theta)}=\sin\theta\left(a+\frac{1}{aE(\theta)}\right),$$
$$Z(\theta)=b\cos\theta+\frac{\cos\theta}{bE(\theta)}=\cos\theta\left(b+\frac{1}{bE(\theta)}\right).$$
Genau diese Geometrie wird in den Implementierungen ausgewertet.
Für das Volumenintegral benötigen wir \(dZ\). Zunächst gilt
$$E'(\theta)=\frac{\sin\theta\cos\theta}{E(\theta)}\left(\frac{1}{a^2}-\frac{1}{b^2}\right).$$
Leitet man nun \(Z(\theta)=b\cos\theta+\frac{\cos\theta}{bE(\theta)}\) ab und vereinfacht, erhält man
$$-Z'(\theta)=\sin\theta\left(b+\frac{1}{bE(\theta)}+\frac{\cos^2\theta}{b}\frac{\frac{1}{a^2}-\frac{1}{b^2}}{E(\theta)^3}\right).$$
Der numerisch integrierte Ausdruck ist also \(R(\theta)^2(-Z'(\theta))\).
Die äußere Form ist zur Ebene \(z=0\) symmetrisch. Daher integriert man die obere Hälfte und verdoppelt. Mit der Scheibenmethode für Rotationskörper gilt
$$V_{\text{outer}}=2\pi\int_0^{\pi/2} R(\theta)^2\bigl(-Z'(\theta)\bigr)\,d\theta.$$
Das Volumen des inneren Ellipsoids ist bekanntlich
$$V_{\text{inner}}=\frac{4\pi a^2 b}{3}.$$
Damit folgt für das Schokoladenvolumen
$$\boxed{C(a,b)=2\pi\int_0^{\pi/2} R(\theta)^2\bigl(-Z'(\theta)\bigr)\,d\theta-\frac{4\pi a^2 b}{3}.}$$
Für \(a=b=1\) ist der Kern die Einheitskugel. Dann ist \(E(\theta)=1\) und damit
$$R(\theta)=2\sin\theta,\qquad Z(\theta)=2\cos\theta,$$
also ist der Außenkörper eine Kugel mit Radius \(2\). Daher
$$C(1,1)=\frac{4\pi}{3}(2^3-1^3)=\frac{28\pi}{3}.$$
Für den zweiten Prüfwert ergibt die numerische Auswertung
$$C(2,1)\approx 60.35475635.$$
Für die Zielparameter erhält man
$$C(3,1)\approx 103.37870096,$$
und genau auf acht Nachkommastellen gerundet geben die Implementierungen diesen Wert aus.
Die C++-, Python- und Java-Implementierungen verwenden denselben Ablauf. Zuerst wird der oben hergeleitete glatte Meridianintegrand ausgewertet. Danach wird auf dem Intervall \([0,\pi/2]\) die Simpson-Regel angewendet, und das Intervall wird per adaptiver Simpson-Quadratur rekursiv weiter unterteilt, bis die lokale Fehlerschätzung klein genug ist oder die maximale Rekursionstiefe erreicht ist.
Sobald das Integral stabil ist, multipliziert die Implementierung mit \(2\pi\), subtrahiert das Ellipsoidvolumen \(\frac{4\pi a^2b}{3}\) und formatiert das Ergebnis auf acht Stellen nach dem Dezimalpunkt. Die beiden bekannten Kontrollwerte dienen dabei als numerischer Konsistenztest für Geometrie und Quadratur.
Adaptive Simpson-Integration ist datenabhängig. Wenn das Intervall in \(m\) akzeptierte Teilintervalle zerlegt wird, beträgt der Aufwand \(O(m)\) Funktionsauswertungen bis auf einen konstanten Faktor. Der zusätzliche Speicher ist \(O(d)\), wobei \(d\) die Rekursionstiefe ist. Für dieses Problem ist der Integrand auf \([0,\pi/2]\) glatt, daher konvergiert das Verfahren schnell und der Speicherbedarf bleibt vernachlässigbar.
Şeker çekirdeği şu dönel elipsoiddir:
$$\frac{x^2+y^2}{a^2}+\frac{z^2}{b^2}=1,$$
eşdeğer olarak \(b^2x^2+b^2y^2+a^2z^2=a^2b^2\). Bunun üzerine, dış yüzey normali boyunca kalınlığı \(1\) mm olan bir çikolata tabakası eklenir. İstenen büyüklük \(C(a,b)\), dış gövdenin hacminden iç elipsoid hacminin çıkarılmasıdır.
Buradaki temel nokta şudur: normal doğrultudaki ofset, yarıçapları doğrudan \((a+1,b+1)\) yapmakla elde edilmez. Dış yüzeyi gerçekten türetmek gerekir.
Cisim \(z\)-ekseni etrafında dönel simetriye sahip olduğundan \((\rho,z)\) meridyen düzleminde çalışmak yeterlidir; burada \(\rho=\sqrt{x^2+y^2}\).
İç elips
$$\frac{\rho^2}{a^2}+\frac{z^2}{b^2}=1$$
eşitliğiyle verilir. Üst yarı için uygun parametreleme
$$\rho(\theta)=a\sin\theta,\qquad z(\theta)=b\cos\theta,\qquad 0\le\theta\le\frac{\pi}{2}.$$
Bu yayı \(z\)-ekseni etrafında döndürmek tüm elipsoidi verir.
Elipsi örtük biçimde
$$F(\rho,z)=\frac{\rho^2}{a^2}+\frac{z^2}{b^2}-1=0$$
olarak yazalım. Dış normal \(\nabla F\) ile aynı doğrultudadır. Parametreleme üzerinde bu yön
$$\left(\frac{\rho}{a^2},\frac{z}{b^2}\right)=\left(\frac{\sin\theta}{a},\frac{\cos\theta}{b}\right)$$
olur. Şimdi
$$E(\theta)=\sqrt{\frac{\sin^2\theta}{a^2}+\frac{\cos^2\theta}{b^2}}$$
tanımını yapalım. Böylece meridyen düzlemindeki dış birim normal
$$\mathbf{n}(\theta)=\frac{1}{E(\theta)}\left(\frac{\sin\theta}{a},\frac{\cos\theta}{b}\right)$$
şeklindedir.
\((\rho(\theta),z(\theta))\) noktasını dış normal boyunca bir birim taşırsak dış profil \((R(\theta),Z(\theta))\) elde edilir:
$$R(\theta)=a\sin\theta+\frac{\sin\theta}{aE(\theta)}=\sin\theta\left(a+\frac{1}{aE(\theta)}\right),$$
$$Z(\theta)=b\cos\theta+\frac{\cos\theta}{bE(\theta)}=\cos\theta\left(b+\frac{1}{bE(\theta)}\right).$$
İmplementasyonların sayısal olarak kullandığı geometri tam olarak budur.
Hacim integrali için \(dZ\) gerekir. Önce
$$E'(\theta)=\frac{\sin\theta\cos\theta}{E(\theta)}\left(\frac{1}{a^2}-\frac{1}{b^2}\right)$$
elde edilir. Sonra \(Z(\theta)=b\cos\theta+\frac{\cos\theta}{bE(\theta)}\) türetilirse sadeleşmiş biçim
$$-Z'(\theta)=\sin\theta\left(b+\frac{1}{bE(\theta)}+\frac{\cos^2\theta}{b}\frac{\frac{1}{a^2}-\frac{1}{b^2}}{E(\theta)^3}\right)$$
olur. Sayısal olarak entegre edilen ifade dolayısıyla \(R(\theta)^2(-Z'(\theta))\)'dir.
Dış şekil \(z=0\) düzlemine göre simetriktir. Bu yüzden üst yarı entegre edilip ikiyle çarpılır. Dönel cisimler için disk yöntemiyle
$$V_{\text{outer}}=2\pi\int_0^{\pi/2} R(\theta)^2\bigl(-Z'(\theta)\bigr)\,d\theta$$
yazılır. İç elipsoid hacmi standart olarak
$$V_{\text{inner}}=\frac{4\pi a^2 b}{3}$$
olduğundan çikolata hacmi
$$\boxed{C(a,b)=2\pi\int_0^{\pi/2} R(\theta)^2\bigl(-Z'(\theta)\bigr)\,d\theta-\frac{4\pi a^2 b}{3}.}$$
\(a=b=1\) durumunda çekirdek birim küredir. Bu halde \(E(\theta)=1\) olur ve
$$R(\theta)=2\sin\theta,\qquad Z(\theta)=2\cos\theta.$$
Yani dış cisim yarıçapı \(2\) olan bir küredir. Buradan
$$C(1,1)=\frac{4\pi}{3}(2^3-1^3)=\frac{28\pi}{3}$$
çıkar. İkinci kontrol değeri için sayısal integral
$$C(2,1)\approx 60.35475635$$
verir. Hedef parametreler için ise
$$C(3,1)\approx 103.37870096$$
elde edilir; implementasyonlar sekiz ondalık basamağa yuvarlayarak bu değeri üretir.
C++, Python ve Java implementasyonları aynı akışı izler. Önce yukarıda türetilen düzgün meridyen integrandını hesaplarlar. Ardından \([0,\pi/2]\) aralığında Simpson kuralını uygularlar ve adaptif Simpson ile aralığı özyinelemeli biçimde bölerek yerel hata tahmini yeterince küçülene veya derinlik sınırına ulaşılana kadar devam ederler.
İntegral kararlı hale geldikten sonra uygulama sonucu \(2\pi\) ile çarpar, \(\frac{4\pi a^2b}{3}\) iç hacmini çıkarır ve çıktıyı sekiz ondalıkla biçimlendirir. Verilen iki kontrol değeri hem geometriyi hem de sayısal integrasyonu doğrulamak için kullanılır.
Adaptif Simpson yöntemi veri bağımlıdır. Aralık \(m\) kabul edilmiş alt aralığa bölünürse toplam iş, sabit katsayılar dışında \(O(m)\) fonksiyon değerlendirmesidir. Ek bellek \(O(d)\)'dir; burada \(d\) özyineleme derinliğidir. Bu problemde integrand \([0,\pi/2]\) üzerinde düzgün olduğundan yakınsama hızlıdır ve bellek kullanımı ihmal edilebilir düzeydedir.
El núcleo del caramelo es el elipsoide de revolución
$$\frac{x^2+y^2}{a^2}+\frac{z^2}{b^2}=1,$$
equivalentemente \(b^2x^2+b^2y^2+a^2z^2=a^2b^2\). Encima se añade una capa de chocolate de grosor \(1\) mm medida a lo largo de la normal exterior. La cantidad pedida \(C(a,b)\) es el volumen del cuerpo exterior menos el volumen del elipsoide interior.
El detalle geométrico clave es que este recubrimiento normal no se obtiene sustituyendo simplemente \((a,b)\) por \((a+1,b+1)\). Hay que construir de forma explícita la superficie desplazada.
Como el sólido tiene simetría de revolución alrededor del eje \(z\), basta trabajar en el plano meridiano \((\rho,z)\), donde \(\rho=\sqrt{x^2+y^2}\). La elipse interior satisface
$$\frac{\rho^2}{a^2}+\frac{z^2}{b^2}=1.$$
Una parametrización conveniente de la mitad superior es
$$\rho(\theta)=a\sin\theta,\qquad z(\theta)=b\cos\theta,\qquad 0\le\theta\le\frac{\pi}{2}.$$
Al girar este arco alrededor del eje \(z\) se recupera el elipsoide completo.
Escribimos la elipse de forma implícita como
$$F(\rho,z)=\frac{\rho^2}{a^2}+\frac{z^2}{b^2}-1=0.$$
Una normal exterior es proporcional a \(\nabla F\), así que sobre la curva parametrizada obtenemos la dirección
$$\left(\frac{\rho}{a^2},\frac{z}{b^2}\right)=\left(\frac{\sin\theta}{a},\frac{\cos\theta}{b}\right).$$
Definimos
$$E(\theta)=\sqrt{\frac{\sin^2\theta}{a^2}+\frac{\cos^2\theta}{b^2}}.$$
Entonces la normal exterior unitaria en el plano meridiano es
$$\mathbf{n}(\theta)=\frac{1}{E(\theta)}\left(\frac{\sin\theta}{a},\frac{\cos\theta}{b}\right).$$
Al mover el punto \((\rho(\theta),z(\theta))\) una unidad en la dirección normal exterior, obtenemos el perfil exterior \((R(\theta),Z(\theta))\):
$$R(\theta)=a\sin\theta+\frac{\sin\theta}{aE(\theta)}=\sin\theta\left(a+\frac{1}{aE(\theta)}\right),$$
$$Z(\theta)=b\cos\theta+\frac{\cos\theta}{bE(\theta)}=\cos\theta\left(b+\frac{1}{bE(\theta)}\right).$$
Ésta es exactamente la geometría evaluada por las implementaciones.
Para la integral de volumen necesitamos \(dZ\). Primero,
$$E'(\theta)=\frac{\sin\theta\cos\theta}{E(\theta)}\left(\frac{1}{a^2}-\frac{1}{b^2}\right).$$
Ahora derivamos \(Z(\theta)=b\cos\theta+\frac{\cos\theta}{bE(\theta)}\). Tras simplificar, resulta
$$-Z'(\theta)=\sin\theta\left(b+\frac{1}{bE(\theta)}+\frac{\cos^2\theta}{b}\frac{\frac{1}{a^2}-\frac{1}{b^2}}{E(\theta)^3}\right).$$
Por tanto, el integrando evaluado numéricamente es \(R(\theta)^2(-Z'(\theta))\).
La figura exterior es simétrica respecto del plano \(z=0\), así que integramos la mitad superior y duplicamos. Mediante el método de discos para sólidos de revolución,
$$V_{\text{outer}}=2\pi\int_0^{\pi/2} R(\theta)^2\bigl(-Z'(\theta)\bigr)\,d\theta.$$
El volumen del elipsoide interior es la fórmula clásica
$$V_{\text{inner}}=\frac{4\pi a^2 b}{3}.$$
En consecuencia, el volumen del chocolate es
$$\boxed{C(a,b)=2\pi\int_0^{\pi/2} R(\theta)^2\bigl(-Z'(\theta)\bigr)\,d\theta-\frac{4\pi a^2 b}{3}.}$$
Cuando \(a=b=1\), el núcleo es la esfera unitaria. Entonces \(E(\theta)=1\), y
$$R(\theta)=2\sin\theta,\qquad Z(\theta)=2\cos\theta,$$
de modo que el cuerpo exterior es una esfera de radio \(2\). Por tanto,
$$C(1,1)=\frac{4\pi}{3}(2^3-1^3)=\frac{28\pi}{3}.$$
Para el segundo punto de control, la evaluación numérica da
$$C(2,1)\approx 60.35475635.$$
Aplicando la misma integral a los parámetros objetivo se obtiene
$$C(3,1)\approx 103.37870096,$$
que es el valor redondeado a ocho decimales por las implementaciones.
Las implementaciones en C++, Python y Java siguen exactamente la misma secuencia. Evalúan el integrando suave del perfil meridiano deducido arriba, aplican la regla de Simpson en \([0,\pi/2]\) y luego refinan el intervalo de forma recursiva con Simpson adaptativo hasta que la estimación local del error cae por debajo de la tolerancia pedida o se alcanza el límite de profundidad.
Una vez estabilizada la integral, la implementación multiplica por \(2\pi\), resta el volumen del elipsoide \(\frac{4\pi a^2b}{3}\) y formatea la salida con ocho cifras decimales. Los dos valores de control publicados sirven como verificación numérica de la geometría y de la cuadratura.
La integración adaptativa de Simpson depende de los datos. Si el intervalo termina dividido en \(m\) subintervalos aceptados, el trabajo total es \(O(m)\) evaluaciones de la función, salvo constantes, y la memoria adicional es \(O(d)\), donde \(d\) es la profundidad recursiva. En este problema el integrando es suave en \([0,\pi/2]\), así que la convergencia es rápida y el uso de memoria es despreciable.
糖果内核是旋转椭球
$$\frac{x^2+y^2}{a^2}+\frac{z^2}{b^2}=1,$$
也就是 \(b^2x^2+b^2y^2+a^2z^2=a^2b^2\)。在它外面再包一层厚度为 \(1\) mm 的巧克力，这个厚度是沿着外法线方向测量的。所求的 \(C(a,b)\) 就是外层实体体积减去原始椭球体积。
这里最容易误解的一点是：法线方向的平移并不等于把半轴直接改成 \((a+1,b+1)\)。必须显式写出偏移后的外表面。
由于该物体绕 \(z\) 轴旋转对称，只需在子午面 \((\rho,z)\) 中工作，其中 \(\rho=\sqrt{x^2+y^2}\)。内层椭圆满足
$$\frac{\rho^2}{a^2}+\frac{z^2}{b^2}=1.$$
其上半部分可以写成
$$\rho(\theta)=a\sin\theta,\qquad z(\theta)=b\cos\theta,\qquad 0\le\theta\le\frac{\pi}{2}.$$
把这条弧线绕 \(z\) 轴旋转，就恢复出整个椭球。
把椭圆写成隐式形式
$$F(\rho,z)=\frac{\rho^2}{a^2}+\frac{z^2}{b^2}-1=0.$$
外法线方向与 \(\nabla F\) 平行，所以在参数曲线上，其方向可以写成
$$\left(\frac{\rho}{a^2},\frac{z}{b^2}\right)=\left(\frac{\sin\theta}{a},\frac{\cos\theta}{b}\right).$$
定义
$$E(\theta)=\sqrt{\frac{\sin^2\theta}{a^2}+\frac{\cos^2\theta}{b^2}}.$$
于是子午面中的外侧单位法向量为
$$\mathbf{n}(\theta)=\frac{1}{E(\theta)}\left(\frac{\sin\theta}{a},\frac{\cos\theta}{b}\right).$$
把点 \((\rho(\theta),z(\theta))\) 沿外法线平移 \(1\) 个单位，就得到外层轮廓 \((R(\theta),Z(\theta))\)：
$$R(\theta)=a\sin\theta+\frac{\sin\theta}{aE(\theta)}=\sin\theta\left(a+\frac{1}{aE(\theta)}\right),$$
$$Z(\theta)=b\cos\theta+\frac{\cos\theta}{bE(\theta)}=\cos\theta\left(b+\frac{1}{bE(\theta)}\right).$$
实现里数值计算的几何对象正是这个偏移后的子午线。
体积积分需要 \(dZ\)。先对 \(E(\theta)\) 求导：
$$E'(\theta)=\frac{\sin\theta\cos\theta}{E(\theta)}\left(\frac{1}{a^2}-\frac{1}{b^2}\right).$$
再对 \(Z(\theta)=b\cos\theta+\frac{\cos\theta}{bE(\theta)}\) 求导并整理，得到
$$-Z'(\theta)=\sin\theta\left(b+\frac{1}{bE(\theta)}+\frac{\cos^2\theta}{b}\frac{\frac{1}{a^2}-\frac{1}{b^2}}{E(\theta)^3}\right).$$
因此数值积分的被积函数就是 \(R(\theta)^2(-Z'(\theta))\)。
外层形状关于平面 \(z=0\) 对称，所以只积上半部分再乘 \(2\)。用旋转体的圆盘法可得
$$V_{\text{outer}}=2\pi\int_0^{\pi/2} R(\theta)^2\bigl(-Z'(\theta)\bigr)\,d\theta.$$
内层椭球体积是标准公式
$$V_{\text{inner}}=\frac{4\pi a^2 b}{3}.$$
于是巧克力体积为
$$\boxed{C(a,b)=2\pi\int_0^{\pi/2} R(\theta)^2\bigl(-Z'(\theta)\bigr)\,d\theta-\frac{4\pi a^2 b}{3}.}$$
当 \(a=b=1\) 时，内核是单位球。此时 \(E(\theta)=1\)，于是
$$R(\theta)=2\sin\theta,\qquad Z(\theta)=2\cos\theta,$$
说明外层实体就是半径 \(2\) 的球。所以
$$C(1,1)=\frac{4\pi}{3}(2^3-1^3)=\frac{28\pi}{3}.$$
第二个校验值通过数值积分得到
$$C(2,1)\approx 60.35475635.$$
对目标参数应用同一公式可得
$$C(3,1)\approx 103.37870096,$$
这就是实现最终按八位小数输出的数值。
C++、Python 和 Java 实现遵循同一条计算链路。它们先计算上面推导出的光滑子午线被积函数，再在区间 \([0,\pi/2]\) 上应用 Simpson 公式，然后使用自适应 Simpson 递归细分区间，直到局部误差估计足够小或者达到递归深度上限。
积分稳定之后，程序将结果乘以 \(2\pi\)，减去椭球体积 \(\frac{4\pi a^2b}{3}\)，最后把答案格式化为八位小数。题目给出的两个检查点被用来验证几何推导和数值积分都没有偏差。
自适应 Simpson 的复杂度依赖于实际细分情况。若区间最终被接受为 \(m\) 个子区间，则总工作量在常数因子下为 \(O(m)\) 次函数求值，额外空间为 \(O(d)\)，其中 \(d\) 是递归深度。对本题而言，被积函数在 \([0,\pi/2]\) 上光滑，因此收敛很快，内存消耗也很小。
Ядро конфеты имеет форму эллипсоида вращения
$$\frac{x^2+y^2}{a^2}+\frac{z^2}{b^2}=1,$$
то есть \(b^2x^2+b^2y^2+a^2z^2=a^2b^2\). Снаружи наносится слой шоколада толщины \(1\) мм, измеренной вдоль внешней нормали к поверхности. Требуется найти \(C(a,b)\): объем внешнего тела минус объем исходного эллипсоида.
Главный геометрический нюанс состоит в том, что нормальное смещение нельзя получить простой заменой \((a,b)\) на \((a+1,b+1)\). Внешнюю поверхность нужно вывести как настоящую параллельную поверхность.
Так как тело симметрично относительно вращения вокруг оси \(z\), достаточно работать в меридианной плоскости \((\rho,z)\), где \(\rho=\sqrt{x^2+y^2}\). Внутренняя эллиптическая кривая задается уравнением
$$\frac{\rho^2}{a^2}+\frac{z^2}{b^2}=1.$$
Ее верхнюю половину удобно параметризовать так:
$$\rho(\theta)=a\sin\theta,\qquad z(\theta)=b\cos\theta,\qquad 0\le\theta\le\frac{\pi}{2}.$$
При вращении этой дуги вокруг оси \(z\) получается весь эллипсоид.
Запишем эллипс неявно:
$$F(\rho,z)=\frac{\rho^2}{a^2}+\frac{z^2}{b^2}-1=0.$$
Внешняя нормаль пропорциональна \(\nabla F\), поэтому на параметризованной кривой ее направление равно
$$\left(\frac{\rho}{a^2},\frac{z}{b^2}\right)=\left(\frac{\sin\theta}{a},\frac{\cos\theta}{b}\right).$$
Введем обозначение
$$E(\theta)=\sqrt{\frac{\sin^2\theta}{a^2}+\frac{\cos^2\theta}{b^2}}.$$
Тогда внешняя единичная нормаль в меридианной плоскости имеет вид
$$\mathbf{n}(\theta)=\frac{1}{E(\theta)}\left(\frac{\sin\theta}{a},\frac{\cos\theta}{b}\right).$$
Если сдвинуть точку \((\rho(\theta),z(\theta))\) на единицу вдоль внешней нормали, получаем внешний профиль \((R(\theta),Z(\theta))\):
$$R(\theta)=a\sin\theta+\frac{\sin\theta}{aE(\theta)}=\sin\theta\left(a+\frac{1}{aE(\theta)}\right),$$
$$Z(\theta)=b\cos\theta+\frac{\cos\theta}{bE(\theta)}=\cos\theta\left(b+\frac{1}{bE(\theta)}\right).$$
Именно эта геометрия численно используется во всех реализациях.
Для интеграла объема требуется \(dZ\). Сначала
$$E'(\theta)=\frac{\sin\theta\cos\theta}{E(\theta)}\left(\frac{1}{a^2}-\frac{1}{b^2}\right).$$
После дифференцирования \(Z(\theta)=b\cos\theta+\frac{\cos\theta}{bE(\theta)}\) и упрощения получаем
$$-Z'(\theta)=\sin\theta\left(b+\frac{1}{bE(\theta)}+\frac{\cos^2\theta}{b}\frac{\frac{1}{a^2}-\frac{1}{b^2}}{E(\theta)^3}\right).$$
Следовательно, численно интегрируемая функция равна \(R(\theta)^2(-Z'(\theta))\).
Внешняя форма симметрична относительно плоскости \(z=0\), поэтому достаточно проинтегрировать верхнюю половину и умножить результат на \(2\). По методу дисков для тел вращения
$$V_{\text{outer}}=2\pi\int_0^{\pi/2} R(\theta)^2\bigl(-Z'(\theta)\bigr)\,d\theta.$$
Объем внутреннего эллипсоида задается стандартной формулой
$$V_{\text{inner}}=\frac{4\pi a^2 b}{3}.$$
Поэтому объем шоколадного слоя равен
$$\boxed{C(a,b)=2\pi\int_0^{\pi/2} R(\theta)^2\bigl(-Z'(\theta)\bigr)\,d\theta-\frac{4\pi a^2 b}{3}.}$$
При \(a=b=1\) ядро является единичной сферой. Тогда \(E(\theta)=1\), и
$$R(\theta)=2\sin\theta,\qquad Z(\theta)=2\cos\theta,$$
то есть внешнее тело есть сфера радиуса \(2\). Следовательно,
$$C(1,1)=\frac{4\pi}{3}(2^3-1^3)=\frac{28\pi}{3}.$$
Для второго контрольного значения численное интегрирование дает
$$C(2,1)\approx 60.35475635.$$
Для целевых параметров получаем
$$C(3,1)\approx 103.37870096,$$
и именно это число реализации выводят после округления до восьми знаков после запятой.
Реализации на C++, Python и Java используют одну и ту же схему. Сначала вычисляется гладкий меридианный подынтегральный множитель, выведенный выше. Затем на отрезке \([0,\pi/2]\) применяется формула Симпсона, после чего отрезок рекурсивно дробится методом адаптивного Симпсона до тех пор, пока локальная оценка ошибки не станет достаточно малой либо не будет достигнут предел глубины рекурсии.
После стабилизации интеграла реализация умножает результат на \(2\pi\), вычитает объем эллипсоида \(\frac{4\pi a^2b}{3}\) и форматирует ответ с восемью знаками после запятой. Два опубликованных контрольных значения служат проверкой и геометрии, и квадратурной схемы.
Сложность адаптивного метода Симпсона зависит от того, сколько раз придется дробить интервал. Если в итоге принято \(m\) подынтервалов, то суммарная работа составляет \(O(m)\) вычислений функции с точностью до постоянного множителя, а дополнительная память равна \(O(d)\), где \(d\) - глубина рекурсии. В этой задаче подынтегральная функция гладкая на \([0,\pi/2]\), поэтому сходимость быстрая, а память практически не расходуется.
لبّ الحلوى هو إهليلج دوراني معادلته
$$\frac{x^2+y^2}{a^2}+\frac{z^2}{b^2}=1,$$
أي بصورة مكافئة \(b^2x^2+b^2y^2+a^2z^2=a^2b^2\). ثم تضاف طبقة شوكولاتة سماكتها \(1\) مم، وهذه السماكة مقاسة على امتداد العمود الخارجي على السطح. المطلوب هو \(C(a,b)\)، أي حجم الجسم الخارجي مطروحًا منه حجم الإهليلج الداخلي.
النقطة الهندسية الأساسية هنا هي أن الإزاحة العمودية لا تعني ببساطة استبدال \((a,b)\) بـ \((a+1,b+1)\). لا بد من اشتقاق السطح الخارجي المزاح اشتقاقًا صريحًا.
بما أن الجسم متماثل دورانيًا حول محور \(z\)، يكفي أن نعمل في المستوى الميريدياني \((\rho,z)\)، حيث \(\rho=\sqrt{x^2+y^2}\). عندئذ يحقق القطع الناقص الداخلي المعادلة
$$\frac{\rho^2}{a^2}+\frac{z^2}{b^2}=1.$$
ويمكن تمثيل نصفه العلوي بالصيغة
$$\rho(\theta)=a\sin\theta,\qquad z(\theta)=b\cos\theta,\qquad 0\le\theta\le\frac{\pi}{2}.$$
وعند تدوير هذا القوس حول محور \(z\) نحصل على الإهليلج كاملًا.
لنكتب القطع الناقص على الصورة الضمنية
$$F(\rho,z)=\frac{\rho^2}{a^2}+\frac{z^2}{b^2}-1=0.$$
اتجاه العمود الخارجي يتناسب مع \(\nabla F\)، ولذلك على المنحنى الممثل بالبارامتر نحصل على الاتجاه
$$\left(\frac{\rho}{a^2},\frac{z}{b^2}\right)=\left(\frac{\sin\theta}{a},\frac{\cos\theta}{b}\right).$$
نعرّف الآن
$$E(\theta)=\sqrt{\frac{\sin^2\theta}{a^2}+\frac{\cos^2\theta}{b^2}}.$$
وعندها يكون متجه العمود الخارجي الواحدي في المستوى الميريدياني هو
$$\mathbf{n}(\theta)=\frac{1}{E(\theta)}\left(\frac{\sin\theta}{a},\frac{\cos\theta}{b}\right).$$
إذا حرّكنا النقطة \((\rho(\theta),z(\theta))\) مسافة واحدة على امتداد العمود الخارجي، فإننا نحصل على المقطع الخارجي \((R(\theta),Z(\theta))\):
$$R(\theta)=a\sin\theta+\frac{\sin\theta}{aE(\theta)}=\sin\theta\left(a+\frac{1}{aE(\theta)}\right),$$
$$Z(\theta)=b\cos\theta+\frac{\cos\theta}{bE(\theta)}=\cos\theta\left(b+\frac{1}{bE(\theta)}\right).$$
وهذه هي الهندسة نفسها التي تعتمد عليها جميع التطبيقات عند الحساب العددي.
من أجل تكامل الحجم نحتاج إلى \(dZ\). أولًا نحصل على
$$E'(\theta)=\frac{\sin\theta\cos\theta}{E(\theta)}\left(\frac{1}{a^2}-\frac{1}{b^2}\right).$$
ثم نشتق \(Z(\theta)=b\cos\theta+\frac{\cos\theta}{bE(\theta)}\)، وبعد التبسيط نجد
$$-Z'(\theta)=\sin\theta\left(b+\frac{1}{bE(\theta)}+\frac{\cos^2\theta}{b}\frac{\frac{1}{a^2}-\frac{1}{b^2}}{E(\theta)^3}\right).$$
إذن الدالة التي يجري تكاملها عدديًا هي \(R(\theta)^2(-Z'(\theta))\).
الشكل الخارجي متماثل بالنسبة إلى المستوى \(z=0\)، لذا نكامل النصف العلوي ثم نضاعف النتيجة. وباستخدام طريقة الأقراص للأجسام الدورانية نحصل على
$$V_{\text{outer}}=2\pi\int_0^{\pi/2} R(\theta)^2\bigl(-Z'(\theta)\bigr)\,d\theta.$$
أما حجم الإهليلج الداخلي فهو الصيغة المعروفة
$$V_{\text{inner}}=\frac{4\pi a^2 b}{3}.$$
وعليه فإن حجم الشوكولاتة يساوي
$$\boxed{C(a,b)=2\pi\int_0^{\pi/2} R(\theta)^2\bigl(-Z'(\theta)\bigr)\,d\theta-\frac{4\pi a^2 b}{3}.}$$
عندما \(a=b=1\) يكون اللب كرة الوحدة. عندها \(E(\theta)=1\)، ومن ثم
$$R(\theta)=2\sin\theta,\qquad Z(\theta)=2\cos\theta,$$
أي أن الجسم الخارجي يصبح كرة نصف قطرها \(2\). لذلك
$$C(1,1)=\frac{4\pi}{3}(2^3-1^3)=\frac{28\pi}{3}.$$
أما قيمة الفحص الثانية فتعطيها الطريقة العددية على الصورة
$$C(2,1)\approx 60.35475635.$$
وبتطبيق الصيغة نفسها على المعاملات المطلوبة نحصل على
$$C(3,1)\approx 103.37870096,$$
وهو العدد الذي تطبعه التطبيقات بعد التقريب إلى ثماني خانات عشرية.
تتبع تطبيقات C++ وPython وJava المسار نفسه. فهي تحسب أولًا الدالة الملساء الناتجة عن مقطع الميريديان كما اشتققناها أعلاه، ثم تطبق قاعدة Simpson على المجال \([0,\pi/2]\)، وبعد ذلك تقسّم المجال تقسيمًا تكراريًا بطريقة Simpson التكيفية حتى يصبح تقدير الخطأ المحلي أصغر من العتبة المطلوبة أو حتى الوصول إلى حد عمق الاستدعاء.
بعد استقرار التكامل، تضرب النتيجة في \(2\pi\)، وتطرح حجم الإهليلج \(\frac{4\pi a^2b}{3}\)، ثم تنسّق الناتج إلى ثماني خانات بعد الفاصلة. كما تُستخدم قيمتا التحقق المعروفتان للتأكد من صحة كل من الهندسة والتكامل العددي.
تعقيد Simpson التكيفية يعتمد على عدد التقسيمات التي يحتاجها التكامل. فإذا انتهى المجال إلى \(m\) مقاطع فرعية مقبولة، كان العمل الكلي \(O(m)\) من حيث عدد تقييمات الدالة حتى عامل ثابت، وكانت الذاكرة الإضافية \(O(d)\) حيث \(d\) هو عمق الاستدعاء التكراري. وفي هذه المسألة تكون الدالة المدمجة ملساء على \([0,\pi/2]\)، لذلك يكون التقارب سريعًا واستهلاك الذاكرة ضئيلًا جدًا.