Нахождение функции по точкам
От: Olegator  
Дата: 27.02.06 00:33
Оценка:
Здравствуйте!

Хочу предложить поток моих ночных программистских мыслей. Следующий код по заданным точкам находит функцию (подбирает коэффициенты многочлена), график которой проходит через эти точки. Зачем это нужно? Не знаю . Ну к примеру, нам нужна формула для нахождения суммы квадратов натуральных чисел от 1 до N. Для этого берём 4 первых значения и вручную считаем суммы. Забиваем всё это в код и получаем нужную формулу. Или ещё можно предсказывать количество страниц в следующей книге о Гарри Поттере (как здесь http://www.cs.rit.edu/~ark/hpbs.shtml).

Также имеется реализация рациональных чисел.

#include <iostream>

using namespace std;

/******** это можно пропустить ********/

template<class IntegerType>
class rational
{
public:
    rational() : _numerator(0), _denominator(1)
    {}

    rational(IntegerType n) : _numerator(n), _denominator(1)
    {}

    rational(IntegerType n, IntegerType m) : _numerator(n), _denominator(m)
    { _normalize(); }

    rational(const rational& rhs) :
        _numerator(rhs.numerator()), _denominator(rhs.denominator())
    {}

private:
    IntegerType _numerator;
    IntegerType _denominator;

public:
    IntegerType numerator() const
    {
        return _numerator;
    }

    IntegerType denominator() const
    {
        return _denominator;
    }

public: // операторы сравнения
    bool operator==(const rational& rhs)
    {
        return (_numerator == rhs.numerator() && _denominator == rhs.denominator());
    }

    bool operator!=(const rational& rhs)
    {
        return (_numerator != rhs.numerator() || _denominator != rhs.denominator());
    }

    bool operator<(const rational& rhs)
    {
        return (_numerator * rhs.denominator() < _denominator * rhs.numerator());
    }

    bool operator>(const rational& rhs)
    {
        return (_numerator * rhs.denominator() > _denominator * rhs.numerator());
    }

    bool operator<=(const rational& rhs)
    {
        return (_numerator * rhs.denominator() <= _denominator * rhs.numerator());
    }

    bool operator>=(const rational& rhs)
    {
        return (_numerator * rhs.denominator() >= _denominator * rhs.numerator());
    }

public: // операторы присваивания
    rational& operator=(const rational& rhs)
    {
        _numerator   = rhs.numerator();
        _denominator = rhs.denominator();

        return *this;
    }

    rational& operator+=(const rational& rhs)
    {
        _numerator   = _numerator * rhs.denominator() + _denominator * rhs.numerator();
        _denominator = _denominator * rhs.denominator();

        _normalize();
        return *this;
    }

    rational& operator-=(const rational& rhs)
    {
        _numerator   = _numerator * rhs.denominator() - _denominator * rhs.numerator();
        _denominator = _denominator * rhs.denominator();

        _normalize();
        return *this;
    }

    rational& operator*=(const rational& rhs)
    {
        _numerator   *= rhs.numerator();
        _denominator *= rhs.denominator();

        _normalize();
        return *this;
    }

    rational& operator/=(const rational& rhs)
    {
        _numerator   *= rhs.denominator();
        _denominator *= rhs.numerator();

        _normalize();
        return *this;
    }

public: // арифметические операции
    rational operator+(const rational& rhs)
    {
        rational res = *this;
        res += rhs;
        return res;
    }

    rational operator-(const rational& rhs)
    {
        rational res = *this;
        res -= rhs;
        return res;
    }

    rational operator*(const rational& rhs)
    {
        rational res = *this;
        res *= rhs;
        return res;
    }

    rational operator/(const rational& rhs)
    {
        rational res = *this;
        res /= rhs;
        return res;
    }

public: // приведение к вещественным типам
    //template<class RealType>
    //operator RealType()
    //{
    //    return (static_cast<RealType>(_numerator)/_denominator);
    //}

private:
    void _normalize()
    {
        IntegerType d = gcd(_numerator, _denominator);

        if(_denominator < 0)
        {
            _numerator   = -_numerator;
            _denominator = -_denominator;
        }

        _numerator   /= d;
        _denominator /= d;
    }

public:
    static IntegerType gcd(IntegerType a, IntegerType b)
    {
        if(a < 0) a = -a;
        if(b < 0) b = -b;

        while(a != 0 && b != 0)
        {
            if(a > b) a %= b; else b %= a;
        }

        return (a + b);
    }
};

template<class IntegerType>
ostream& operator<<(ostream& os, const rational<IntegerType>& r)
{
    os << r.numerator();
    if(r.denominator() != 1)
        os << '/' << r.denominator();
    return os;
}

template<class IntegerType>
IntegerType power(IntegerType a, IntegerType b)
{
    IntegerType res = 1;
    while(b > 0)
    {
        if(b % 2 == 0)
        {
            a *= a;
            b /= 2;
        }
        else
        {
            res *= a;
            b--;
        }
    }
    return res;
}

/******** тут начинаем смотреть ********/

const int max_n = 10;

typedef rational<long long> rational_t;

rational_t equation_set[max_n][max_n+1];
rational_t coefficient [max_n];

int main()
{
    int arg[max_n] = {1,2,3,4};   // аргументы (x)
    int val[max_n] = {1,5,14,30}; // значения (y)
    int n = 4;                    // количество точек

    // делаем первоначальную систему уравнений
    for(int i = 0; i < n; i++)
    {
        for(int j = 0; j < n; j++)
            equation_set[i][j] = power(arg[i], j);
        equation_set[i][n] = val[i];
    }

    // вычитаем
    for(int i = 0; i < n - 1; i++)
    {
        for(int j = i + 1; j < n; j++)
        {
            rational_t c = equation_set[j][i] / equation_set[i][i];

            for(int k = i; k < n + 1; k++)
            {
                equation_set[j][k] -= c * equation_set[i][k];
            }
        }
    }

    // находим коэффициенты
    for(int i = n - 1; i >= 0; i--)
    {
        for(int j = n - 1; j > i; j--)
            equation_set[i][n] -= coefficient[j] * equation_set[i][j];
        coefficient[i] = equation_set[i][n] / equation_set[i][i];
    }

    // проверяем
    for(int i = 0; i < n; i++)
    {
        rational_t v = 0;
        for(int j = 0; j < n; j++)
            v += coefficient[j] * power(arg[i], j);
        cout << v << '\n';
    }

    for(int i = n - 1; i > 0; i--)
        cout << '(' << coefficient[i] << ")*" << "x^" << i << " + ";
    cout << '(' << coefficient[0] << ")\n";

    return 0;
}


Или лучше ночью мне спать?
... << RSDN@Home 1.1.4 stable rev. 510>>
 
Подождите ...
Wait...
Пока на собственное сообщение не было ответов, его можно удалить.