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

Тема в разделе "WASM.A&O", создана пользователем _qwe8013, 6 май 2017.

  1. _qwe8013

    _qwe8013 Active Member

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

    Собственно код:
    Код (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, остальные - вспомогательные подфункции.
     
    Последнее редактирование: 6 май 2017
  2. UbIvItS

    UbIvItS Well-Known Member

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

    _qwe8013 Active Member

    Публикаций:
    2
    Регистрация:
    30 ноя 2016
    Сообщения:
    123
    Спасибо.
     
  4. RET

    RET Well-Known Member

    Публикаций:
    17
    Регистрация:
    5 янв 2008
    Сообщения:
    789
    Адрес:
    Jabber: darksys@sj.ms
    Сначала проверяйте 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.         }
    Сорри я в кучу написал - кому надо - поймет, это выдержка из моего криптора
     
  5. SadKo

    SadKo Владимир Садовников

    Публикаций:
    8
    Регистрация:
    4 июн 2007
    Сообщения:
    1.610
    Адрес:
    г. Санкт-Петербург
    Не сильно вдавался в алгоритм (не сформулированы исходные данные и способ решения), но при беглом осмотре кода мне тут кажется, сам алгоритм страдает.
    Зачем возводить комплексные числа в степень, если можно поделить исходный полином на сопряжённую пару?
    Кстати говоря, рекомендую интерполяционный полином представлять не в виде
    a0xn + a1xn-1 + ... + aNx0,
    а в виде произведения многочленов первой (вещественные корни) либо второй (комплексные корни) степени:
    (a0x2 + a1x + a2)⋅(b0x2 + b1x + b2)⋅(c0x + c1) + ...

    Какая отсюда выгода:
    1. каждый многочлен в произведении можно нормировать, вынеся a0, b0, c0 за скобки → имеем прирост в точности при вычислении каждого многочлена
    2. при вычислении квадратного многочлена нет большой потери точности
    3. квадратные многочлены вычисляются быстро, при правильной организации цикла можно даже не вычислять каждый раз значения x и x2, а, вычислив один раз, просто подставлять коэффициенты следующего многочлена на следующей итерации.
     
    Mikl___ нравится это.
  6. UbIvItS

    UbIvItS Well-Known Member

    Публикаций:
    0
    Регистрация:
    5 янв 2007
    Сообщения:
    6.074
    SadKo, любая методика тихо-мирно упирается в битность числа.. во многих случаях без длинной арифы совсем никак. а методы представления полинома не столько про "контроль погрешности", сколько про "скорость алгоса".
     
  7. RET

    RET Well-Known Member

    Публикаций:
    17
    Регистрация:
    5 янв 2008
    Сообщения:
    789
    Адрес:
    Jabber: darksys@sj.ms
    UbIvItS,
    +1
    тут надо SEE юзать, если есть поддержка
    Причем в компиляторе указывать типа /fp:precise
     
  8. SadKo

    SadKo Владимир Садовников

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

    RET Well-Known Member

    Публикаций:
    17
    Регистрация:
    5 янв 2008
    Сообщения:
    789
    Адрес:
    Jabber: darksys@sj.ms
    SadKo,
    ну это от опций компилятора тоже зависит
     
  10. UbIvItS

    UbIvItS Well-Known Member

    Публикаций:
    0
    Регистрация:
    5 янв 2007
    Сообщения:
    6.074
    дискретность делает своё ЧЁРНОЕ :)
     
  11. CurryHowardIsomorphism

    CurryHowardIsomorphism Member

    Публикаций:
    0
    Регистрация:
    13 май 2017
    Сообщения:
    97
    Куда-то ещё стандартные complex float/double не завезли?!
     
  12. SadKo

    SadKo Владимир Садовников

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

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

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

    CurryHowardIsomorphism Member

    Публикаций:
    0
    Регистрация:
    13 май 2017
    Сообщения:
    97
    https://godbolt.org/g/5C16Xm да вроде векторизуются
     
  14. _qwe8013

    _qwe8013 Active Member

    Публикаций:
    2
    Регистрация:
    30 ноя 2016
    Сообщения:
    123
    Дурацкий вопрос: а как мне получить эти корни?
     
  15. SadKo

    SadKo Владимир Садовников

    Публикаций:
    8
    Регистрация:
    4 июн 2007
    Сообщения:
    1.610
    Адрес:
    г. Санкт-Петербург
    Ну а сейчас как вы интерполяционный полином получаете? Рассказали бы хоть алгоритм, чтобы долго код не ковырять и не гадать, что там происходит.
    Если составляете систему уравнений, то тут нужно подумать, как минимизировать возведение в степень N комплексного числа. К тому же, лучше при возведении, наверное, пользоваться полярной формой:
    (AejP)n = AnejPn
    Это, по идее, может позволить вынести мусор за скобки и при возможности сократить с нулём или оставить для конечных вычислений.
     
  16. Indy_

    Indy_ Well-Known Member

    Публикаций:
    4
    Регистрация:
    29 апр 2011
    Сообщения:
    4.775
  17. _qwe8013

    _qwe8013 Active Member

    Публикаций:
    2
    Регистрация:
    30 ноя 2016
    Сообщения:
    123
    Составляю систему уравнений с матрицей Вандермонда.
    Не хотел использовать лишний раз иррациональные функции, хотя возможно вы правы, и лучше в полярной форме.
     
  18. SadKo

    SadKo Владимир Садовников

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

    _qwe8013 Active Member

    Публикаций:
    2
    Регистрация:
    30 ноя 2016
    Сообщения:
    123
    Так, а синус-косинус придётся, там же сложение.
     
  20. SadKo

    SadKo Владимир Садовников

    Публикаций:
    8
    Регистрация:
    4 июн 2007
    Сообщения:
    1.610
    Адрес:
    г. Санкт-Петербург
    Мне видится что-то вроде такого:
    Код (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.