#include <iostream>
#include <stdio.h>
#include <math.h>
using namespace std;
double x[12000] = {0}, y[12000] = {0}, z[12000] = {0};
double f(double x, double y, double z)
{
return 6*x;//2*z-y+x-2;
}
double Y(double x)
{
double C1 = 1, C2 = 1;
return pow(x, 3) + C1*x + C2;//C1*exp(x)+C2*x*exp(x)+x;
}
double L(int n)
{
double s = 0;
for (int i = 0; i < n; i++)
{
s += sqrt(pow(x[i]-x[i+1], 2) + pow(y[i]-y[i+1], 2));
}
return s;
}
void d2(double x0, double y0, double z0, double b, int n)
{
double tau = ((b-x0)+.0)/(n+.0);
x[0] = x0; y[0] = y0; z[0] = z0;
double k1, k2, k3, k4;
double m1, m2, m3, m4;
//cout << "\tn = " << n << "\ntau = " << tau << endl;
for (int i = 0; i <= n; i++)
{
x[i] = x0 + i*tau;
//cout << "x_i = " << x0+i*tau << endl;
m1 = z[i];
k1 = f(x[i], y[i], z[i]);
m2 = z[i]+tau*k1/2.;
k2 = f(x[i] + tau/2., y[i] + tau*m1/2., z[i] + tau*k1/2.);
m3 = z[i]+tau*k2/2.;
k3 = f(x[i] + tau/2., y[i] + tau*m2/2., z[i] + tau*k2/2.);
m4 = z[i]+tau*k3;
k4 = f(x[i] + tau, y[i] + tau*m3, z[i] + tau*k3);
y[i+1] = y[i] + tau * (m1 + 2.*m2 + 2.*m3 + m4)/6.;
z[i+1] = z[i] + tau * (k1 + 2.*k2 + 2.*k3 + k4)/6.;
}
}
double f2(double x0, double y0, double beta, double b, int n)
{
d2(x0, y0, beta, b, n);
return L(n);
}
double gs(double x0, double y0, double b, int n)
{
double x11 = -100, x22 = 100, x1, x2;
double md = (1 + sqrt(5)) / 2, eps = 1e-6;
int i = 0;
//std::cout << "i a b |a - b| \n";
while (fabs(x11 - x22) >= eps)
{
x1 = x22 - (x22 - x11) / md;
x2 = x11 + (x22 - x11) / md;
if (f2(x0, y0, x1, b, n) >= f2(x0, y0, x2, b, n)) x11 = x1; else x22 = x2;
//std::cout << "\n" << ++i << " " << x11 << " " << x22 << " " << fabs(x11 - x22) << "\n";
}
return (x11+x22)/2;
}
int main() {
double x0, y0, z0, b, n, eps = 1e-8;
//cout << "n = ";
cin >> n;
//cout << "x0 = ";
cin >> x0;
//cout << "y0 = ";
cin >> y0;
//cout << "z0 = ";
cin >> z0;
//cout << "b = ";
cin >> b;
d2(x0, y0, z0, b, n);
cout << "x_i \t y_i \t \t Y_i \t \t |Y_i - y_i|\n";
for (int i = 0; i <= 10; i++)
{
cout << x[i*int(n/10)] << "\t\t " << y[i*int(n/10)]
<< " \t\t " << Y(x[i*int(n/10)])
<< " \t " << fabs(Y(x[i*int(n/10)]) - y[i*int(n/10)]) << endl;
}
double zL, zL_;
zL_ = gs(x0, y0, b, n/2);
zL = gs(x0, y0, b, n);
while (fabs(zL-zL_) > eps)
{
zL_ = zL;n*=2;
zL = gs(x0, y0, b, n);
d2(x0, y0, zL, b, n);
//cout << "\n=\nn = " << n << "\nmin zL = " << zL << "\nmin zL_ = " << zL_
//<< "\nmin L = " << L(n) << endl;
printf("\n=\nn=%d\nzL=%.20lf\nL=%.20lf\n", n, zL, L(n));
}
return 0;
}
I2luY2x1ZGUgPGlvc3RyZWFtPgojaW5jbHVkZSA8c3RkaW8uaD4KI2luY2x1ZGUgPG1hdGguaD4KdXNpbmcgbmFtZXNwYWNlIHN0ZDsKCmRvdWJsZSB4WzEyMDAwXSA9IHswfSwgeVsxMjAwMF0gPSB7MH0sIHpbMTIwMDBdID0gezB9OwoKZG91YmxlIGYoZG91YmxlIHgsIGRvdWJsZSB5LCBkb3VibGUgeikKewoJcmV0dXJuIDYqeDsvLzIqei15K3gtMjsKfQoKZG91YmxlIFkoZG91YmxlIHgpCnsKCWRvdWJsZSBDMSA9IDEsIEMyID0gMTsKCXJldHVybiBwb3coeCwgMykgKyBDMSp4ICsgQzI7Ly9DMSpleHAoeCkrQzIqeCpleHAoeCkreDsKfQoKZG91YmxlIEwoaW50IG4pCnsKCWRvdWJsZSBzID0gMDsKCWZvciAoaW50IGkgPSAwOyBpIDwgbjsgaSsrKQoJewoJCXMgKz0gc3FydChwb3coeFtpXS14W2krMV0sIDIpICsgcG93KHlbaV0teVtpKzFdLCAyKSk7Cgl9CglyZXR1cm4gczsKfQoKdm9pZCBkMihkb3VibGUgeDAsIGRvdWJsZSB5MCwgZG91YmxlIHowLCBkb3VibGUgYiwgaW50IG4pCnsKCWRvdWJsZSB0YXUgPSAoKGIteDApKy4wKS8obisuMCk7Cgl4WzBdID0geDA7IHlbMF0gPSB5MDsgelswXSA9IHowOwoJZG91YmxlIGsxLCBrMiwgazMsIGs0OwoJZG91YmxlIG0xLCBtMiwgbTMsIG00OwoJLy9jb3V0IDw8ICJcdG4gPSAiIDw8IG4gPDwgIlxudGF1ID0gIiA8PCB0YXUgPDwgZW5kbDsKCWZvciAoaW50IGkgPSAwOyBpIDw9IG47IGkrKykKCXsKCQl4W2ldID0geDAgKyBpKnRhdTsKCQkvL2NvdXQgPDwgInhfaSA9ICIgPDwgeDAraSp0YXUgPDwgZW5kbDsKCQltMSA9IHpbaV07CgkJazEgPSBmKHhbaV0sIHlbaV0sIHpbaV0pOwoJCW0yID0geltpXSt0YXUqazEvMi47CgkJazIgPSBmKHhbaV0gKyB0YXUvMi4sIHlbaV0gKyB0YXUqbTEvMi4sIHpbaV0gKyB0YXUqazEvMi4pOwoJCW0zID0geltpXSt0YXUqazIvMi47CgkJazMgPSBmKHhbaV0gKyB0YXUvMi4sIHlbaV0gKyB0YXUqbTIvMi4sIHpbaV0gKyB0YXUqazIvMi4pOwoJCW00ID0geltpXSt0YXUqazM7CgkJazQgPSBmKHhbaV0gKyB0YXUsIHlbaV0gKyB0YXUqbTMsIHpbaV0gKyB0YXUqazMpOwoJCXlbaSsxXSA9IHlbaV0gKyB0YXUgKiAobTEgKyAyLiptMiArIDIuKm0zICsgbTQpLzYuOwoJCXpbaSsxXSA9IHpbaV0gKyB0YXUgKiAoazEgKyAyLiprMiArIDIuKmszICsgazQpLzYuOwoJfQp9Cgpkb3VibGUgZjIoZG91YmxlIHgwLCBkb3VibGUgeTAsIGRvdWJsZSBiZXRhLCBkb3VibGUgYiwgaW50IG4pCnsKCWQyKHgwLCB5MCwgYmV0YSwgYiwgbik7CglyZXR1cm4gTChuKTsKfQoKZG91YmxlIGdzKGRvdWJsZSB4MCwgZG91YmxlIHkwLCBkb3VibGUgYiwgaW50IG4pCnsKCWRvdWJsZSB4MTEgPSAtMTAwLCB4MjIgPSAxMDAsIHgxLCB4MjsKCWRvdWJsZSBtZCA9ICgxICsgc3FydCg1KSkgLyAyLCBlcHMgPSAxZS02OwoJaW50IGkgPSAwOwoJLy9zdGQ6OmNvdXQgPDwgImkgICBhICAgYiAgIHxhIC0gYnwgICBcbiI7Cgl3aGlsZSAoZmFicyh4MTEgLSB4MjIpID49IGVwcykKCXsKCQl4MSA9IHgyMiAtICh4MjIgLSB4MTEpIC8gbWQ7CgkJeDIgPSB4MTEgKyAoeDIyIC0geDExKSAvIG1kOwoJCWlmIChmMih4MCwgeTAsIHgxLCBiLCBuKSA+PSBmMih4MCwgeTAsIHgyLCBiLCBuKSkgeDExID0geDE7IGVsc2UgeDIyID0geDI7CgkJLy9zdGQ6OmNvdXQgPDwgIlxuIiA8PCArK2kgPDwgIiAiIDw8IHgxMSA8PCAiICIgPDwgeDIyIDw8ICIgIiA8PCBmYWJzKHgxMSAtIHgyMikgPDwgIlxuIjsKCX0KCXJldHVybiAoeDExK3gyMikvMjsKfQoKaW50IG1haW4oKSB7Cglkb3VibGUgeDAsIHkwLCB6MCwgYiwgbiwgZXBzID0gMWUtODsKCS8vY291dCA8PCAibiA9ICI7IAoJY2luID4+IG47CgkvL2NvdXQgPDwgIngwID0gIjsKCWNpbiA+PiB4MDsKCS8vY291dCA8PCAieTAgPSAiOwoJY2luID4+IHkwOwoJLy9jb3V0IDw8ICJ6MCA9ICI7CgljaW4gPj4gejA7CgkvL2NvdXQgPDwgImIgPSAiOwoJY2luID4+IGI7CglkMih4MCwgeTAsIHowLCBiLCBuKTsKCWNvdXQgPDwgInhfaSBcdCB5X2kgXHQgXHQgWV9pIFx0IFx0IHxZX2kgLSB5X2l8XG4iOwoJZm9yIChpbnQgaSA9IDA7IGkgPD0gMTA7IGkrKykKCXsKCQljb3V0IDw8IHhbaSppbnQobi8xMCldIDw8ICJcdFx0ICIgPDwgeVtpKmludChuLzEwKV0gCgkJPDwgIiBcdFx0ICIgPDwgWSh4W2kqaW50KG4vMTApXSkgCgkJPDwgIiBcdCAiIDw8IGZhYnMoWSh4W2kqaW50KG4vMTApXSkgLSB5W2kqaW50KG4vMTApXSkgPDwgZW5kbDsKCX0KCWRvdWJsZSB6TCwgekxfOwoJekxfID0gZ3MoeDAsIHkwLCBiLCBuLzIpOwoJekwgPSBncyh4MCwgeTAsIGIsIG4pOwoJd2hpbGUgKGZhYnMoekwtekxfKSA+IGVwcykKCXsKCQkKCQl6TF8gPSB6TDtuKj0yOwoJCXpMID0gZ3MoeDAsIHkwLCBiLCBuKTsKCQlkMih4MCwgeTAsIHpMLCBiLCBuKTsKCQkvL2NvdXQgPDwgIlxuPVxubiA9ICIgPDwgbiA8PCAiXG5taW4gekwgPSAiIDw8IHpMIDw8ICJcbm1pbiB6TF8gPSAiIDw8IHpMXwoJCS8vPDwgIlxubWluIEwgPSAiIDw8IEwobikgPDwgZW5kbDsKCQlwcmludGYoIlxuPVxubj0lZFxuekw9JS4yMGxmXG5MPSUuMjBsZlxuIiwgbiwgekwsIEwobikpOwoJfQoJCglyZXR1cm4gMDsKfQ==