ООП / ООП / 3-2 / Mathem

//---------------------------------------------------------------------------

#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;
}