С.И.Хашин
Для нахождения рационального приближения действительного числа применяют "непрерывные", или "цепные" дроби, см.,например, в здесь.
Рассмотрим, для примера, рациональные приближения для e и π:
π ≈ |
e ≈ |
Выпишем получающиеся дроби и их погрешности:
n | value | deviation | value | deviation |
---|---|---|---|---|
π | e | |||
1 | 3 | 0.14 | 2 | 0.71 |
2 | 22/7 | 0.0012 | 3 | 0.28 |
3 | 333/106 | 0.000083 | 8/3 | 0.051 |
4 | 355/113 | 0.26e-6 | 11/4 | 0.031 |
5 | 103993/33102 | 0.57e-9 | 19/7 | 0.004 |
6 | 104348/33215 | 0.33e-9 | 87/32 | 0.0004 |
7 | 208341/66317 | 0.12e-9 | 106/39 | 0.00033 |
8 | 312689/99532 | 0.29e-10 | 193/71 | 0.000028 |
9 | 833719/265381 | 0.87e-11 | 1264/465 | 0.22e-5 |
10 | 1146408/364913 | 0.16e-11 | 1457/536 | 0.17e-5 |
11 | 4272943/1360120 | 0.40e-12 | 2721/1001 | 0.11e-6 |
12 | 5419351/1725033 | 0.22e-13 | 23225/8544 | 0.67e-8 |
13 | 80143857/25510582 | 0.58e-15 | 25946/9545 | 0.55e-8 |
14 | 49171/18089 | 0.28e-9 | ||
15 | 517656/190435 | 0.13e-10 | ||
16 | 566827/208524 | 0.11e-10 | ||
17 | 1084483/398959 | 0.48e-12 | ||
18 | 13580623/4996032 | 0.20e-13 | ||
19 | 14665106/5394991 | 0.17e-13 | ||
20 | 28245729/10391023 | 0.61e-15 |
Реализация этого алгоритма на 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); // less0 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, π).Применим к нему функцию rationalize со значениями tol равными 0.1, 0.01, ....
tol | w[0] | w[1] | deviation: abs(π - w[1]/w[0]) |
---|---|---|---|
1e-1 | 1 | 3 | 0.14 |
1e-2 | 7 | 22 | 0.0012 |
1e-3 | 7 | 22 | 0.0012 |
1e-4 | 113 | 355 | 0.26e-6 |
1e-5 | 113 | 355 | 0.26e-6 |
1e-6 | 113 | 355 | 0.26e-6 |
1e-7 | 113 | 355 | 0.26e-6 |
1e-8 | 33215 | 104348 | 0.33e-9 |
1e-9 | 33215 | 104348 | 0.33e-9 |
1e-10 | 33215 | 104348 | 0.33e-9 |
1e-11 | 99532 | 312689 | 0.3e-10 |
1e-12 | 364913 | 1146408 | 0.2e-11 |
1e-13 | 1725033 | 5419351 | 0.2e-13 |
1e-14 | 1725033 | 5419351 | 0.2e-13 |
1e-15 | 25510582 | 80143857 | 0.6e-15 |
То есть, получающиеся результаты согласуются с аппроксимацией числа π, описанной выше.
Пусть вектор v имеет вид:
v = (1, π, e) ≈ (1, 3.1415926535897932385, 2.7182818284590452354).Применим к нему функцию rationalize со значениями tol равными 0.1, 0.01, ....
tol | w[0] | w[1] | w[2] |
---|---|---|---|
1e-1 | 1 | 3 | 3 |
1e-2 | 7 | 22 | 19 |
1e-3 | 7 | 22 | 19 |
1e-4 | 685 | 2152 | 1862 |
1e-5 | 685 | 2152 | 1862 |
1e-6 | 20446 | 64233 | 55578 |
1e-7 | 123368 | 387572 | 335349 |
1e-8 | 1131457 | 3554578 | 3075619 |
1e-9 | 1131457 | 3554578 | 3075619 |
1e-10 | 13968034 | 43881874 | 37969053 |
1e-11 | 150377371 | 472424444 | 408768075 |
1e-12 | 3738677499 | 11745401765 | 10162779108 |
1e-13 | 3738677499 | 11745401765 | 10162779108 |
1e-14 | 37084888791 | 116505614185 | 100807179311 |
1e-15 | 37084888791 | 116505614185 | 100807179311 |
Рассмотрим ещё один пример. Пусть у нас есть уравнение кривой второго порядка:
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);получим R=57.697 и, следовательно f(x,y) пропорционально
3*x^2 - 4*x*y + 5*y^2 + 7*x - 3*y - 11.
tol=
Vector v:
: