Погрешность операций с плавающей точкой

В общем понадобилось мне тут получить интерполяционный комплексный полином по большому количеству точек. Написал программу, и всё бы хорошо, но чем больше точек, тем больше погрешность вычисления (сравнил значение полинома в данных точках с заданным). На сколько я понимаю, существуют методы компенсации погрешности, но в гугле я нашёл только это и это. Попробовал применить эти методы в моём случае, но получилась джигурда. Что можно сделать в моём случае? (длинную арифметику не предлагать, она очень медленная).

Собственно код:

Код (C):

1. #define FloatType float//double
2.  
3. struct COMPLEX
4. {
5.     FloatType Re;
6.     FloatType Im;
7. };
8.  
9. struct COMPLEX_POINT
10. {
11.     COMPLEX x;
12.     COMPLEX y;
13. };
14.  
15. void AddComplex(COMPLEX* Dst, COMPLEX* a, COMPLEX* b)
16. {
17.     Dst->Re = a->Re + b->Re;
18.     Dst->Im = a->Im + b->Im;
19. }
20.  
21. void SubComplex(COMPLEX* Dst, COMPLEX* Minuent, COMPLEX* Subtrahend)
22. {
23.     Dst->Re = Minuent->Re - Subtrahend->Re;
24.     Dst->Im = Minuent->Im - Subtrahend->Im;
25. }
26.  
27. void MultComplex(COMPLEX* Dst, COMPLEX* a, COMPLEX* b)
28. {
29.     FloatType Re;
30.  
31.     Re = a->Re*b->Re - a->Im*b->Im;
32.     Dst->Im = a->Re*b->Im + a->Im*b->Re;
33.     Dst->Re = Re;
34. }
35.  
36. void DivComplex(COMPLEX* Dst, COMPLEX* Num, COMPLEX* Denom)
37. {
38.     COMPLEX s;
39.     s.Re = Denom->Re;
40.     s.Im = -Denom->Im;
41.  
42.     FloatType d = Denom->Re*Denom->Re + Denom->Im*Denom->Im;
43.  
44.     MultComplex(Dst, Num, &s);
45.     Dst->Re /= d;
46.     Dst->Im /= d;
47. }
48.  
49. void IntPowComplex(COMPLEX* Dst, COMPLEX* Base, int Exponent)
50. {
51.     Dst->Re = 1;
52.     Dst->Im = 0;
53.     if (Exponent != 0)
54.     {
55.         for (unsigned int i = 0; i < abs(Exponent); i++)
56.         {
57.             MultComplex(Dst, Dst, Base);
58.         }
59.  
60.         if (Exponent < 0)
61.         {
62.             COMPLEX tmp;
63.             tmp.Re = 1;
64.             tmp.Im = 0;
65.  
66.             DivComplex(Dst, &tmp, Dst);
67.         }
68.     }
69. }
70.  
71. void GetM(COMPLEX_POINT Points[], unsigned int StartIndex, unsigned int MiddleIndex, unsigned int EndIndex, COMPLEX* M)
72. {
73.     COMPLEX tmp;
74.  
75.     M->Re = 1;
76.     M->Im = 0;
77.  
78.     if (MiddleIndex > StartIndex)
79.     {
80.         for (unsigned int i = StartIndex; i < MiddleIndex; i++)
81.         {
82.             SubComplex(&tmp, &Points[i].x, &Points[MiddleIndex].x);
83.             MultComplex(M, M, &tmp);
84.         }
85.     }
86.  
87.     if (MiddleIndex < EndIndex)
88.     {
89.         for (unsigned int i = MiddleIndex + 1; i <= EndIndex; i++)
90.         {
91.             SubComplex(&tmp, &Points[MiddleIndex].x, &Points[i].x);
92.             MultComplex(M, M, &tmp);
93.         }
94.     }
95. }
96.  
97. void GetCoef(COMPLEX_POINT Points[], unsigned int NumOfPoints, COMPLEX LocArr[], unsigned int index, COMPLEX* Result)
98. {
99.     Result->Re = 0;
100.     Result->Im = 0;
101.  
102.     for (unsigned int i = 0; i <= index; i++)
103.     {
104.         COMPLEX M, tmp;
105.         GetM(Points, 0, i, index, &M);
106.  
107.         if (i % 2 == 1)
108.         {
109.             tmp.Re = -LocArr[i].Re;
110.             tmp.Im = -LocArr[i].Im;
111.         }
112.         else
113.         {
114.             tmp.Re = LocArr[i].Re;
115.             tmp.Im = LocArr[i].Im;
116.         }
117.         DivComplex(&tmp, &tmp, &M);
118.         AddComplex(Result, Result, &tmp);
119.     }
120. }
121.  
122. void Interpoly(COMPLEX_POINT Points[], unsigned int NumOfPoints, COMPLEX Poly[])
123. {
124.     if (NumOfPoints > 0)
125.     {
126.         COMPLEX* LocArr = (COMPLEX*)malloc(sizeof(COMPLEX)*NumOfPoints);
127.      
128.         for (unsigned int i = 0; i < NumOfPoints; i++)
129.         {
130.             LocArr[i] = Points[i].y;
131.         }
132.  
133.         for (unsigned int i = NumOfPoints; i>0; i--)
134.         {
135.             GetCoef(Points, NumOfPoints, LocArr, i - 1, &Poly[i - 1]);
136.  
137.             for (unsigned int j = 0; j < i; j++)
138.             {
139.                 COMPLEX tmp;
140.  
141.                 IntPowComplex(&tmp, &Points[j].x, i - 1);
142.                 MultComplex(&tmp, &tmp, &Poly[i - 1]);
143.                 SubComplex(&LocArr[j], &LocArr[j], &tmp);
144.             }
145.         }
146.  
147.         free(LocArr);
148.     }
149. }

Многочлен строит функция Interpoly, остальные - вспомогательные подфункции.

 

_qwe8013 https://duckduckgo.com/?q=interpolation+complex+polynom&atb=v12&ia=web пробегись по ссылям -- думаю, всё найдёшь вплоть до уже готовых решений. по идее, в matlab/maple/mathematica -- все эти фичи ужо забиты, тч материала для чтения/тестов/использования хватает.

 

Спасибо.

 

Сначала проверяйте SSE поддержку, интел маны в кучу

Код (Text):

1. int isCPUIDsupported (void)
2. {
3.  
4.   __asm
5.   {
6.        push ecx ; save ecx
7.      pushfd ; push original EFLAGS
8.        pop eax ; get original EFLAGS
9.        mov ecx, eax ; save original EFLAGS
10.      xor eax, 200000h ; flip bit 21 in EFLAGS
11.         push eax ; save new EFLAGS value on stack
12.        popfd ; replace current EFLAGS value
13.         pushfd ; get new EFLAGS
14.      pop eax ; store new EFLAGS in EAX
15.        xor eax, ecx ; Bit 21 of flags at 200000h will be 1 if CPUID exists
16.      shr eax, 21  ; Shift bit 21 bit 0 and return it
17.      push ecx
18.         popfd ; restore bit 21 in EFLAGS first
19.       pop ecx ; restore ecx
20.    }
21. }
22. __asm
23.         {
24.             _emit 0x0F
25.                 _emit 0x1F
26.                 _emit 0x00
27.                 _emit 0x0F
28.                 _emit 0x1F
29.                 _emit 0x00
30.                 mov eax, 0x1
31.                 _emit 0x0F
32.                 _emit 0x1F
33.                 _emit 0x00
34.                 cpuid
35.                 _emit 0x0F
36.                 _emit 0x1F
37.                 _emit 0x00
38.                 test edx, 1<<25
39.                 _emit 0x0F
40.                 _emit 0x1F
41.                 _emit 0x00
42.                 jnz _SSE
43.                 _emit 0x0F
44.                 _emit 0x1F
45.                 _emit 0x00
46.                 int 3
47.                 call ExitProcess
48.                 _emit 0x0F
49.                 _emit 0x1F
50.                 _emit 0x00
51. _SSE:
52.             _emit 0x0F
53.                 _emit 0x1F
54.                 _emit 0x00
55.                 _emit 0x0F
56.                 _emit 0x1F
57.                 _emit 0x00
58.                 _emit 0x0F
59.                 _emit 0x1F
60.                 _emit 0x00
61.                 _emit 0x0F
62.                 _emit 0x1F
63.                 _emit 0x00
64.         }

Сорри я в кучу написал - кому надо - поймет, это выдержка из моего криптора

 

Не сильно вдавался в алгоритм (не сформулированы исходные данные и способ решения), но при беглом осмотре кода мне тут кажется, сам алгоритм страдает.
Зачем возводить комплексные числа в степень, если можно поделить исходный полином на сопряжённую пару?
Кстати говоря, рекомендую интерполяционный полином представлять не в виде
a₀⋅xⁿ + a₁⋅x^n-1 + ... + a_N⋅x⁰,
а в виде произведения многочленов первой (вещественные корни) либо второй (комплексные корни) степени:
(a₀⋅x² + a₁⋅x + a₂)⋅(b₀⋅x² + b₁⋅x + b₂)⋅(c₀⋅x + c₁) + ...

Какая отсюда выгода:
1. каждый многочлен в произведении можно нормировать, вынеся a₀, b₀, c₀ за скобки → имеем прирост в точности при вычислении каждого многочлена
2. при вычислении квадратного многочлена нет большой потери точности
3. квадратные многочлены вычисляются быстро, при правильной организации цикла можно даже не вычислять каждый раз значения x и x², а, вычислив один раз, просто подставлять коэффициенты следующего многочлена на следующей итерации.

 

SadKo, любая методика тихо-мирно упирается в битность числа.. во многих случаях без длинной арифы совсем никак. а методы представления полинома не столько про "контроль погрешности", сколько про "скорость алгоса".

 

UbIvItS,
+1
тут надо SEE юзать, если есть поддержка
Причем в компиляторе указывать типа /fp:precise

 

Я с вами не согласен, потому что совсем недавно (где-то полгода назад) приходилось решать подобную задачу. Нужно было проектировать цифровые БИХ-фильтры 32-го порядка и рисовать их АЧХ. Не поверите, но задачу решил спокойно на самых обычных float'ах.

 

SadKo,
ну это от опций компилятора тоже зависит

 

дискретность делает своё ЧЁРНОЕ

 

Куда-то ещё стандартные complex float/double не завезли?!

 

Не факт. Важно правильный алгоритм выбрать. Так, оптимизируя каскады biquad-фильтров и применяя различные техники оптимизации вплоть до SSE- и AVX-ассемблирования, я получил где-то в 10 раз более быструю реализацию фильтров по сравнению с начальной.
Аналогичное можно сказать про FFT: самая тупая реализация свёртки с использованием самой тупой реализации алгоритма FFT вносит меньше ошибок в результат, нежели расчёт свёртки "в лоб".

Кстати, интересная вещь: в DSP фильтры конструируют из biquad-каскадов (полином 2 степени в числителе и полином 2 степени в знаменателе). Это объясняется тем, что при большем порядке полинома, чем 2, фильтр резко теряет устойчивость. Все мои попытки сделать фильтр на полиномах четвёртой степени потерпели фиаско: при определённых условиях фильтр просто терял устойчивость. Зато цепочка biquad-фильтров работает на ура.

Стандартные complex float/double плохо оптимизируются и векторизуются. Проще вообще отказаться от комплексных чисел и вывести формулы для real и imaginary частей.

 

https://godbolt.org/g/5C16Xm да вроде векторизуются

 

Дурацкий вопрос: а как мне получить эти корни?

 

Ну а сейчас как вы интерполяционный полином получаете? Рассказали бы хоть алгоритм, чтобы долго код не ковырять и не гадать, что там происходит.
Если составляете систему уравнений, то тут нужно подумать, как минимизировать возведение в степень N комплексного числа. К тому же, лучше при возведении, наверное, пользоваться полярной формой:
(A ⋅ e^jP)ⁿ = Aⁿ ⋅ e^jPn
Это, по идее, может позволить вынести мусор за скобки и при возможности сократить с нулём или оставить для конечных вычислений.

 

_qwe8013,

Abramowitz and Stegun: Handbook of Mathematical Functions

Там много коэффициентов для подобных фунок.

 

Составляю систему уравнений с матрицей Вандермонда.

Не хотел использовать лишний раз иррациональные функции, хотя возможно вы правы, и лучше в полярной форме.

 

Так фишка полярной формы, как раз, в том, что у вас целая степень всегда, поэтому exp(j *P) у вас единичный вектор и будет просто вращаться N раз на свой угол поворота против часовой стрелки. Вам даже иррациональные функции (кроме sqrt) применять не придётся. То есть, exp(j*P) можно по-прежнему хранить как пару чисел re + im, но нормированную до единицы по модулю, а сам модуль вынести за скобки.

 

Так, а синус-косинус придётся, там же сложение.

 

Мне видится что-то вроде такого:

Код (C++):

1.  
2. typedef struct complex_t {
3.     float a; // Аргумент комплексного числа
4.     float re, im; // Единичный нормированный вектор вращения
5. } complex_t;
6.  
7. void complex_init(complex_t *c, float re, float im)
8. {
9.    c->a = sqrtf(re*re + im*im);
10.    if (c->a > 0.0f)
11.    {
12.        c->re = re / c->a;
13.        c->im = im / c->a;
14.    }
15.    else
16.    {
17.       c->re = 0.0f;
18.       c->im = 0.0f;
19.    }
20. }
21.  
22. void complex_mul(complex_t *x, complex_t *a, complex_t *b)
23. {
24.     x->a = a->a * b->a;
25.     float re = a->re*b->re - a->im*b->im;
26.     float im = a->re*b->im + a->im*b->re;
27.     x->re = re;
28.     x->im = im;
29. }
30.  
31. void complex_add(complex_t *x, complex_t *a, complex_t *b)
32. {
33.    if (a->a <= 0.0f)
34.    {
35.       *x = *b;
36.       return;
37.    }
38.    else if (b->a <= 0.0f)
39.    {
40.       *x = *a;
41.       return;
42.    }
43.  
44.    // Здесь, на самом деле, хорошо бы заранее вынести общий множитель (скорее всего это будет sqrtf(a->a * b->a))  и преобразовать формулу, чтобы погрешность стала меньше
45.    float re = a->re * a->a + b->re * b->a;
46.    float im = a->im * a->a + b->im * b->a;
47.  
48.    float a = sqrtf (re*re + im*im);
49.    if (a <= 0.0f)
50.    {
51.       x->a = 0.0f;
52.       x->re = 0.0f;
53.       x->im = 0.0f;
54.       return;
55.    }
56.  
57.    x->a = a;
58.    x->re = re / a;
59.    x->im = re / a;
60. }
61.  
62.  

Соответственно, стараемся, чтобы поле a также было невелико, а держалось близко к 1.