#include <iostream>
#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) ;
cout << "\n =\n n = " << n << "\n min zL = " << zL << "min zL_ = " << zL_
<< "\n min L = " << L( n) << endl;
while ( fabs ( zL- zL_) > eps)
{
zL_ = zL; n* = 2 ;
zL = gs( x0, y0, b, n) ;
d2( x0, y0, zL, b, n) ;
cout << "\n =\n n = " << n << "\n min zL = " << zL << "min zL_ = " << zL_
<< "\n min L = " << L( n) << endl;
}
return 0 ;
}
I2luY2x1ZGUgPGlvc3RyZWFtPgojaW5jbHVkZSA8bWF0aC5oPgp1c2luZyBuYW1lc3BhY2Ugc3RkOwoKZG91YmxlIHhbMTIwMDBdID0gezB9LCB5WzEyMDAwXSA9IHswfSwgelsxMjAwMF0gPSB7MH07Cgpkb3VibGUgZihkb3VibGUgeCwgZG91YmxlIHksIGRvdWJsZSB6KQp7CglyZXR1cm4gNip4Oy8vMip6LXkreC0yOwp9Cgpkb3VibGUgWShkb3VibGUgeCkKewoJZG91YmxlIEMxID0gMSwgQzIgPSAxOwoJcmV0dXJuIHBvdyh4LCAzKSArIEMxKnggKyBDMjsvL0MxKmV4cCh4KStDMip4KmV4cCh4KSt4Owp9Cgpkb3VibGUgTChpbnQgbikKewoJZG91YmxlIHMgPSAwOwoJZm9yIChpbnQgaSA9IDA7IGkgPCBuOyBpKyspCgl7CgkJcyArPSBzcXJ0KHBvdyh4W2ldLXhbaSsxXSwgMikgKyBwb3coeVtpXS15W2krMV0sIDIpKTsKCX0KCXJldHVybiBzOwp9Cgp2b2lkIGQyKGRvdWJsZSB4MCwgZG91YmxlIHkwLCBkb3VibGUgejAsIGRvdWJsZSBiLCBpbnQgbikKewoJZG91YmxlIHRhdSA9ICgoYi14MCkrLjApLyhuKy4wKTsKCXhbMF0gPSB4MDsgeVswXSA9IHkwOyB6WzBdID0gejA7Cglkb3VibGUgazEsIGsyLCBrMywgazQ7Cglkb3VibGUgbTEsIG0yLCBtMywgbTQ7CgkvL2NvdXQgPDwgIlx0biA9ICIgPDwgbiA8PCAiXG50YXUgPSAiIDw8IHRhdSA8PCBlbmRsOwoJZm9yIChpbnQgaSA9IDA7IGkgPD0gbjsgaSsrKQoJewoJCXhbaV0gPSB4MCArIGkqdGF1OwoJCS8vY291dCA8PCAieF9pID0gIiA8PCB4MCtpKnRhdSA8PCBlbmRsOwoJCW0xID0geltpXTsKCQlrMSA9IGYoeFtpXSwgeVtpXSwgeltpXSk7CgkJbTIgPSB6W2ldK3RhdSprMS8yLjsKCQlrMiA9IGYoeFtpXSArIHRhdS8yLiwgeVtpXSArIHRhdSptMS8yLiwgeltpXSArIHRhdSprMS8yLik7CgkJbTMgPSB6W2ldK3RhdSprMi8yLjsKCQlrMyA9IGYoeFtpXSArIHRhdS8yLiwgeVtpXSArIHRhdSptMi8yLiwgeltpXSArIHRhdSprMi8yLik7CgkJbTQgPSB6W2ldK3RhdSprMzsKCQlrNCA9IGYoeFtpXSArIHRhdSwgeVtpXSArIHRhdSptMywgeltpXSArIHRhdSprMyk7CgkJeVtpKzFdID0geVtpXSArIHRhdSAqIChtMSArIDIuKm0yICsgMi4qbTMgKyBtNCkvNi47CgkJeltpKzFdID0geltpXSArIHRhdSAqIChrMSArIDIuKmsyICsgMi4qazMgKyBrNCkvNi47Cgl9Cn0KCmRvdWJsZSBmMihkb3VibGUgeDAsIGRvdWJsZSB5MCwgZG91YmxlIGJldGEsIGRvdWJsZSBiLCBpbnQgbikKewoJZDIoeDAsIHkwLCBiZXRhLCBiLCBuKTsKCXJldHVybiBMKG4pOwp9Cgpkb3VibGUgZ3MoZG91YmxlIHgwLCBkb3VibGUgeTAsIGRvdWJsZSBiLCBpbnQgbikKewoJZG91YmxlIHgxMSA9IC0xMDAsIHgyMiA9IDEwMCwgeDEsIHgyOwoJZG91YmxlIG1kID0gKDEgKyBzcXJ0KDUpKSAvIDIsIGVwcyA9IDFlLTY7CglpbnQgaSA9IDA7CgkvL3N0ZDo6Y291dCA8PCAiaSAgIGEgICBiICAgfGEgLSBifCAgIFxuIjsKCXdoaWxlIChmYWJzKHgxMSAtIHgyMikgPj0gZXBzKQoJewoJCXgxID0geDIyIC0gKHgyMiAtIHgxMSkgLyBtZDsKCQl4MiA9IHgxMSArICh4MjIgLSB4MTEpIC8gbWQ7CgkJaWYgKGYyKHgwLCB5MCwgeDEsIGIsIG4pID49IGYyKHgwLCB5MCwgeDIsIGIsIG4pKSB4MTEgPSB4MTsgZWxzZSB4MjIgPSB4MjsKCQkvL3N0ZDo6Y291dCA8PCAiXG4iIDw8ICsraSA8PCAiICIgPDwgeDExIDw8ICIgIiA8PCB4MjIgPDwgIiAiIDw8IGZhYnMoeDExIC0geDIyKSA8PCAiXG4iOwoJfQoJcmV0dXJuICh4MTEreDIyKS8yOwp9CgppbnQgbWFpbigpIHsKCWRvdWJsZSB4MCwgeTAsIHowLCBiLCBuLCBlcHMgPSAxZS04OwoJLy9jb3V0IDw8ICJuID0gIjsgCgljaW4gPj4gbjsKCS8vY291dCA8PCAieDAgPSAiOwoJY2luID4+IHgwOwoJLy9jb3V0IDw8ICJ5MCA9ICI7CgljaW4gPj4geTA7CgkvL2NvdXQgPDwgInowID0gIjsKCWNpbiA+PiB6MDsKCS8vY291dCA8PCAiYiA9ICI7CgljaW4gPj4gYjsKCWQyKHgwLCB5MCwgejAsIGIsIG4pOwoJY291dCA8PCAieF9pIFx0IHlfaSBcdCBcdCBZX2kgXHQgXHQgfFlfaSAtIHlfaXxcbiI7Cglmb3IgKGludCBpID0gMDsgaSA8PSAxMDsgaSsrKQoJewoJCWNvdXQgPDwgeFtpKmludChuLzEwKV0gPDwgIlx0XHQgIiA8PCB5W2kqaW50KG4vMTApXSAKCQk8PCAiIFx0XHQgIiA8PCBZKHhbaSppbnQobi8xMCldKSAKCQk8PCAiIFx0ICIgPDwgZmFicyhZKHhbaSppbnQobi8xMCldKSAtIHlbaSppbnQobi8xMCldKSA8PCBlbmRsOwoJfQoJZG91YmxlIHpMLCB6TF87Cgl6TF8gPSBncyh4MCwgeTAsIGIsIG4vMik7Cgl6TCA9IGdzKHgwLCB5MCwgYiwgbik7Cgljb3V0IDw8ICJcbj1cbm4gPSAiIDw8IG4gPDwgIlxubWluIHpMID0gIiA8PCB6TCA8PCAibWluIHpMXyA9ICIgPDwgekxfCgkJPDwgIlxubWluIEwgPSAiIDw8IEwobikgPDwgZW5kbDsKCXdoaWxlIChmYWJzKHpMLXpMXykgPiBlcHMpCgl7CgkJCgkJekxfID0gekw7bio9MjsKCQl6TCA9IGdzKHgwLCB5MCwgYiwgbik7CgkJZDIoeDAsIHkwLCB6TCwgYiwgbik7CgkJY291dCA8PCAiXG49XG5uID0gIiA8PCBuIDw8ICJcbm1pbiB6TCA9ICIgPDwgekwgPDwgIm1pbiB6TF8gPSAiIDw8IHpMXwoJCTw8ICJcbm1pbiBMID0gIiA8PCBMKG4pIDw8IGVuZGw7CgkJCgl9CgkKCXJldHVybiAwOwp9