В начало

Rationalization of vector of double

Sergei Khashin

For the rational approximation of real numbers one can use "continued", or "chain" fraction, see,for example, here.

Consider, for example, a rational approximation for e and π:

π ≈

e ≈

Let us write the resulting fractions and their errors (deviations):

Implementation of this algorithm is given in the C ++ function chainFract, which using an auxiliary function 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);
}

How can that be, if we must at the same time bring some rational numbers? Of course, you can zoom each of them separately, but the denominators to get, most likely different. So when casting to a common denominator, the numbers will be too great. For example, bringing thus the pair (π, e) with the accuracy up to 1e-10 will receive:

This, of course, not a very good approximation.

Here is a function in C ++ that allows to solve this problem.
rationalize.h     - header
rationalize.cpp - realization.

In the unit there is, in fact, the one and only function:

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

In fact, there are several auxiliary, very simple functions, implementing some basic operations with vectors of double numbers:

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

The rationalize for the given vector v of length n finds the actual number R such that the vector R*v has an almost integer coordinates. Denote by w vector R*v with coordinates rounded to a whole number. Then the vector w/R will be different from the original vector v not more than tol.

Consider a few examples.

Let vector v is a:

    v = (1, π).

Apply function rationalize with value of tol equals to 0.1, 0.01, ....

That is, the obtained results are consistent with an approximation of π, as described above.

Let vector v is:

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

Apply function rationalize with value of tol equals to 0.1, 0.01, ....

Consider another example. Suppose we have a second-order equation of the curve:

  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

As a result of execution:

    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:

Rationalization: