В начало

Рационализация действительного вектора

С.И.Хашин

Для нахождения рационального приближения действительного числа применяют "непрерывные", или "цепные" дроби, см.,например, в здесь.

Рассмотрим, для примера, рациональные приближения для e и π:

π ≈

e ≈

Выпишем получающиеся дроби и их погрешности:

Реализация этого алгоритма на C++ дана в функции chainFract, использующей вспомогательную функцию chainFractStep:

bool chainFractStep( double x, int &p, int &q, int maxQ) // one step of chainFract
{
    int maxN1;
    if( q<=0 ) { q=1; p=round(x); return true; }
    maxN1 = maxQ/q;
    double err=q*x-p;
    if( fabs(err)*maxQ <1 ) return false;
    int n1 = round(1/err);
    if( abs(n1)> maxN1 ) return false;
    p = p*n1+1;
    q *=n1;
    if( n1<0 ) {q=-q; p=-p; }
    return true;
}
double chainFract( double x, int &p, int &q, int maxQ) // rat.approx p/q of x, return tolerance
{
    double err=1;
    for( q=1; chainFractStep(x,p,q,maxQ); ) ;
    return fabs(x-(double)p/q);
}

Как же быть, если нам надо одновременно приблизить несколько рациональных чисел? Конечно, можно приблизить каждое из них по-отдельности, но знаменатели получаться, скорее всего различные. А значит при приведении к общему знаменателю, числа окажутся слишком велики. Например, приближая таким образом пару (π, e) с точностью до 1e-10 получим:

Это, конечно же, не слишком хорошее приближение.

Приведем функцию на С++, позволяющую решать эту задачу.
rationalize.h - заголовочный файл
rationalize.cpp - реализация.

В приведенном модуле имеется, фактически, одна-единственная функция

double rationalize(int n, double *v, double tol); // calculate R: v[i]*R - (almost) integer with deviation tol. Return R

double vdbl_len   (int N, double *a);                   // |a|
double vdbl_scpod (int N, double *a, double *b);        // (a,b)
void   vdbl_add   (int N, double *a, double *b,double R);// a += R*b
double vdbl_minAbs(int N, double *a, int &i);           // minimal abs.value of non-zero elemtnts = a[i]
void   vdbl_zero  (int N, double *a, double eps);       // less 0
double deviationI (double R);                           // deviation R from integer
double vdbl_devi  (int n, double *a, double co);        // deviation a*co from integer

Функция rationalize для данного вектора v длины n находит действительное число R такое, что вектор R*v имеет почти целые координаты. Обозначим через w вектор R*v с координатами, округлёнными до целых. Тогда вектор w/R будет отличаться от исходного вектора v не более, чем на tol.

Рассмотрим несколько примеров работы функции.

Пусть вектор v имеет вид:

    v = (1, π).

То есть, получающиеся результаты согласуются с аппроксимацией числа π, описанной выше.

Пусть вектор v имеет вид:

    v = (1, π, e) ≈ (1,  3.1415926535897932385, 2.7182818284590452354).

Рассмотрим ещё один пример. Пусть у нас есть уравнение кривой второго порядка:

  f(x,y) = a[0]*x^2 + a[1]*x*y + a[2]*y^2 + a[3]*x + a[4]*y + a[5],

    a[0] = 0.051995332184145352452
    a[1] =-0.069327109578860469936
    a[2] = 0.086658886973575587420
    a[3] = 0.12132244176300582239
    a[4] =-0.051995332184145352452
    a[5] =-0.19064955134186629232

В результате выполнения команд

    double v[] = {  0.051995332184145352452, -0.069327109578860469936,  0.086658886973575587420,  
                    0.12132244176300582239 , -0.051995332184145352452, -0.19064955134186629232   };
    int n = sizeof(v) / sizeof(double);
    double R = rationalize(n, v, 1e-10);

  3*x^2 - 4*x*y + 5*y^2 + 7*x - 3*y - 11.

tol=
Vector v:

Рационализация: