Комплексные числа (complex.h) и очевидное отставание точности - PullRequest
Новая
       6

Комплексные числа (complex.h) и очевидное отставание точности

1 голос
/ 07 октября 2019

Я решил немного поэкспериментировать с complex.h и столкнулся с тем, что считаю очень любопытной проблемой. 

int mandelbrot(long double complex c, int lim)
{
    long double complex z = c;
    for(int i = 0; i < lim; ++i, z = cpowl(z,2)+c)
    {
        if(creall(z)*creall(z)+cimagl(z)*cimagl(z) > 4.0)
            return 0;
    }
    return 1;
}

int mandelbrot2(long double cr, long double ci, int lim)
{
    long double zr = cr;
    long double zi = ci;
    for(int i = 0; i < lim; ++i, zr = zr*zr-zi*zi+cr, zi = 2*zr*zi+ci)
    {
        if(zr*zr+zi*zi > 4.0)
            return 0;
    }
    return 1;
}

Эти функции не работают одинаково. Если мы введем -2.0 + 0.0i и предел выше 17, последний вернет 1, что является правильным для любого предела, тогда как первый вернет 0, по крайней мере, в моей системе. GCC 9.1.0, Ryzen 2700x.

Я не могу понять, как это может случиться. Я имею в виду, хотя я не совсем понимаю, как complex.h работает за кулисами, для этого конкретного примера не имеет смысла, что результаты должны отклоняться, как это.

Во время написания я замечаю cpowl (z, 2) + c, и попытался изменить его на z * z + c, что помогло, однако после быстрого теста я обнаружил, что поведение все еще отличается. Ex. -1,3 + 0,1 * I, lim = 18.

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

--- edit ---

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

#include <stdio.h>
#include <complex.h>

int mandelbrot(long double complex c, int lim)
{
    long double complex z = c;
    for(int i = 0; i < lim; ++i, z = z*z+c)
    {
        if(creall(z)*creall(z)+cimagl(z)*cimagl(z) > 4.0)
            return 0;
    }
    return 1;
}

int mandelbrot2(long double cr, long double ci, int lim)
{
    long double zr = cr;
    long double zi = ci;
    long double tmp;
    for(int i = 0; i < lim; ++i)
    {
        if(zr*zr+zi*zi > 4.0) return 0;
        tmp = zi;
        zi = 2*zr*zi+ci;
        zr = zr*zr-tmp*tmp+cr;
    }
    return 1;
}

int main()
{
    long double complex c = -2.0+0.0*I;
    printf("%i\n",mandelbrot(c,100));
    printf("%i\n",mandelbrot2(-2.0,0.0,100));
    return 0;
}

cpowl () все еще портит ситуацию, но, если я захочу, я мог бы просто создать свою собственную реализацию.

Ответы [ 2 ]

1 голос
/ 07 октября 2019

Вторая неправильная функция, а не первая.

В выражении в третьем предложении for: 

zr = zr*zr-zi*zi+cr, zi = 2*zr*zi+ci

Расчет ziиспользует новое значение zr, а не текущее. Вам нужно будет сохранить результаты этих двух вычислений во временных переменных, а затем присвоить их обратно zr и zi: 

int mandelbrot2(long double cr, long double ci, int lim)
{
    long double zr = cr;
    long double zi = ci;
    for(int i = 0; i < lim; ++i)
    {   
        printf("i=%d, z=%Lf%+Lfi\n", i, zr, zi);
        if(zr*zr+zi*zi > 4.0)
            return 0;
        long double new_zr = zr*zr-zi*zi+cr;
        long double new_zi = 2*zr*zi+ci;
        zr = new_zr; 
        zi = new_zi;
    }
    return 1;
}

Кроме того, использование cpowl для простого возведения в квадрат приведет кнеточности, которых можно избежать, если в этом случае просто использовать z*z.

0 голосов
/ 07 октября 2019

Разница для входа −2 + 0 i

cpowl является неточной. Экспонирование является сложной функцией для реализации, и при ее вычислении могут возникать различные ошибки. В macOS 10.14.6 z в подпрограмме mandelbrot принимает следующие значения в последовательных итерациях: 

z = -2 + 0 i.
z = 2 + 4.33681e-19 i.
z = 2 + 1.73472e-18 i.
z = 2 + 6.93889e-18 i.
z = 2 + 2.77556e-17 i.
z = 2 + 1.11022e-16 i.
z = 2 + 4.44089e-16 i.
z = 2 + 1.77636e-15 i.
z = 2 + 7.10543e-15 i.
z = 2 + 2.84217e-14 i.
z = 2 + 1.13687e-13 i.
z = 2 + 4.54747e-13 i.
z = 2 + 1.81899e-12 i.
z = 2 + 7.27596e-12 i.
z = 2 + 2.91038e-11 i.
z = 2 + 1.16415e-10 i.
z = 2 + 4.65661e-10 i.

Таким образом, после первоначальной ошибки получается 2 + 4,33681 • 10 −19 i, z продолжает расти (правильно, в результате математики, а не просто ошибок с плавающей запятой), пока не станет достаточно большим, чтобы пройти тест, сравнивая квадрат его абсолютного значения с 4. (тест не сразу фиксирует избыток, потому что квадрат мнимой части настолько мал, что при добавлении к квадрату реальной части он теряется при округлении.)

В отличие от этого, если мы заменим z = cpowl(z,2)+c на z = z*z + c, z остается 2 (то есть 2 + 0i). Как правило, операции в z*z также испытывают некоторые ошибки округления, но не так сильно, как при cpowl.

Разница для входа -1,3 + 0,1 i

Для этого входаразница вызвана неправильным вычислением на этапе обновления цикла for: 

++i, zr = zr*zr-zi*zi+cr, zi = 2*zr*zi+ci

, который использует новое значение zr при вычислении zi. Это можно исправить, вставив long double t; и изменив шаг обновления на 

++i, t = zr*zr - zi*zi + cr, zi = 2*zr*zi + ci, zr = t
...