Добавил:
sergey123
Опубликованный материал нарушает ваши авторские права? Сообщите нам.
Вуз:
Предмет:
Файл:ООП / ООП / new_lab3_2 / Mathem
.cpp//---------------------------------------------------------------------------
#include "math.h"
#pragma hdrstop
#include "Mathem.h"
//---------------------------------------------------------------------------
#pragma package(smart_init)
typedef unsigned int nat;
//*********************** 13.05.2004 ************************//
//
// Кубический корень.
//
//*********************** 13.05.2004 ************************//
static double _cbrt ( double x )
{
double s = 1.;
while ( x < 1. )
{
x *= 8.;
s *= 0.5;
}
while ( x > 8. )
{
x *= 0.125;
s *= 2.;
}
double r = 1.5;
r -= 1./3. * ( r - x / ( r * r ) );
r -= 1./3. * ( r - x / ( r * r ) );
r -= 1./3. * ( r - x / ( r * r ) );
r -= 1./3. * ( r - x / ( r * r ) );
r -= 1./3. * ( r - x / ( r * r ) );
r -= 1./3. * ( r - x / ( r * r ) );
return r * s;
}
double cbrt ( double x )
{
if ( x > 0 ) return _cbrt ( x ); else
if ( x < 0 ) return-_cbrt (-x ); else
return 0.;
}
//*************************************************************//
//
// Решение алгебраических уравнений 2, 3 и 4 степени.
// Возвращаемое значение - количество найденных корней.
//
//*************************************************************//
// x^2 + a * x + b = 0
nat root2 ( double a, double b, double * x )
{
if ( b == 0 )
{
x[0] = -a;
x[1] = 0.;
}
else
{
a *= -0.5;
double d = a * a - b;
if ( d < 0 ) return 0;
d = sqrt ( d );
x[0] = a > 0 ? a + d : a - d;
x[1] = b / x[0];
}
return 2;
}
// x^3 + a * x + b = 0
static double cubic ( double a, double b )
{
double s = 1.;
while ( b + a > -1. )
{
a *= 4.;
b *= 8.;
s *= 0.5;
}
while ( b + a + a < -8. )
{
a *= 0.25;
b *= 0.125;
s *= 2.;
}
double x = 1.5;
x -= ( x * ( x * x + a ) + b ) / ( 3. * x * x + a );
x -= ( x * ( x * x + a ) + b ) / ( 3. * x * x + a );
x -= ( x * ( x * x + a ) + b ) / ( 3. * x * x + a );
x -= ( x * ( x * x + a ) + b ) / ( 3. * x * x + a );
x -= ( x * ( x * x + a ) + b ) / ( 3. * x * x + a );
x -= ( x * ( x * x + a ) + b ) / ( 3. * x * x + a );
x -= ( x * ( x * x + a ) + b ) / ( 3. * x * x + a );
x -= ( x * ( x * x + a ) + b ) / ( 3. * x * x + a );
x -= ( x * ( x * x + a ) + b ) / ( 3. * x * x + a );
return x * s;
}
// x^3 + p * x + q = 0
nat root3s ( double p, double q, double * x )
{
if ( q < 0 )
{
*x = cubic ( p, q );
}
else
if ( q > 0 )
{
*x = -cubic ( p, -q );
}
else
{
*x = 0;
}
return 1 + root2 ( *x, p + (*x)*(*x), x+1 );
}
// x^3 + a * x^2 + b * x + c = 0
nat root3 ( double a, double b, double c, double * x )
{
if ( c == 0 )
{
*x = 0;
}
else
{
const double a3 = a / 3.;
const double p = b - a3 * a;
const double q = c - ( a3 * a3 + p ) * a3;
if ( q < 0 )
{
*x = cubic ( p, q );
}
else
if ( q > 0 )
{
*x = -cubic ( p, -q );
}
else
{
*x = 0;
}
*x -= a3;
const double t = *x * ( *x * 3. + a * 2. ) + b;
if ( fabs ( t ) > 1e-3 )
{
*x -= ( *x * ( *x * ( *x + a ) + b ) + c ) / t;
}
a += *x;
b += *x * a;
}
return 1 + root2 ( a, b, x+1 );
}
// x^4 + p * x^2 + q * x + r = 0
nat root4s ( double p, double q, double r, double * x )
{
if ( r == 0 ) // это условие добавлено 21.10.2008
{
*x = 0;
return 1 + root3s ( p, q, x+1 );
}
if ( q == 0 )
{
double t[2];
nat n = root2 ( p, r, t );
if ( n == 0 ) return 0;
if ( t[0] >= 0 )
{
x[0] = sqrt ( t[0] );
x[1] = - x[0];
n = 2;
}
else
{
n = 0;
}
if ( t[1] >= 0 )
{
x[n] = sqrt ( t[1] );
x[n+1] = - x[n];
n += 2;
}
return n;
}
nat n = root3 ( p + p, p * p - 4. * r, -q * q, x );
double a = *x;
if ( n == 3 )
{
if ( a < x[1] ) a = x[1];
if ( a < x[2] ) a = x[2];
}
if ( a < 0 ) return 0;
p += a;
a = sqrt ( a );
const double b = q / a;
n = root2 ( a, 0.5 * ( p - b ), x );
n+= root2 ( -a, 0.5 * ( p + b ), x+n );
return n;
}
// x^4 + a * x^3 + b * x^2 + c * x + d = 0
nat root4 ( double a, double b, double c, double d, double * x )
{
if ( a == 0 )
{
return root4s ( b, c, d, x );
}
if ( d == 0 )
{
*x = 0;
return 1 + root3 ( a, b, c, x+1 );
}
const double e = a / 4;
const double h = e * e;
const double p = -6 * h + b;
const double q = 8 * h * e - 2 * b * e + c;
const double r = -3 * h * h + b * h - c * e + d;
const nat n = root4s ( p, q, r, x );
for ( nat i = 0; i < n; ++i )
{
x[i] -= e;
}
return n;
}
Соседние файлы в папке new_lab3_2