Возведение комплексных чисел в степень (Формула Муавра)

[BratSinot]

Доброго времени суток!

Имеем комплексное число Z. По Формуле Муавра возвожу в натуральную степень:

Код: Выделить всё

#define my_cpown(Z, n) (powl(cabsl(Z), n)*(cos(atanl(cimagl(Z)/creall(Z))*n)+sin(atanl(cimagl(Z)/creall(Z))*n)*I))

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

Все, в этой статье все нашел:
http://arbuz.uz/x_complex.html

[aumit]

Стыдно не знать алгебру. Читай учебник алгебры за 10-11 класс.

[drBatty]

жуть какая... тут про "образование" говорили, а его получается теперь и нет вовсе...

[BratSinot]

Стыдно не знать алгебру. Читай учебник алгебры за 10-11 класс.

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

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

[aumit]

Стыдно не знать алгебру. Читай учебник алгебры за 10-11 класс.

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

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

Однако ошибка кроется в определение угла, а это простая тригонометрия.
Область определения обратного тангенса (-pi/2, pi/2), а значит по твоей формуле не входит вся левая полуплоскость включая значения в +-pi/2.
Если вещественная часть комплексного числа 0, то вылетает ошибка деления на 0.

[BratSinot]

Однако ошибка кроется в определение угла, а это простая тригонометрия.
Область определения обратного тангенса (-pi/2, pi/2), а значит по твоей формуле не входит вся левая полуплоскость включая значения в +-pi/2.

Так, а почему мы именно Re проверяем на отрицательность? Почему Im не проверяем?

Если вещественная часть комплексного числа 0, то вылетает ошибка деления на 0.

Там никогда не будет ровно 0, поэтому проблем не будет.

[Reboot]

Область определения обратного тангенса (-pi/2, pi/2)

Стыдно не знать алгебру.

Занудство и провал.

BratSinot, залог успеха лежит в atan2, почти всегда нужен он, а не atan.
http://beej.us/guide/bgc/output/html/multipage/atan.html

ещё говорят, что можно man 3 cpow в консоле написать и получить бесплатное счастье.

[kyzic]

жуть какая... тут про "образование" говорили, а его получается теперь и нет вовсе...

как страшно жить!

[drBatty]

Да и изивините меня, я про этот случай не мог найти в справочниках по Высшей Математике, что вы хотите от учебника?

раньше ЕМНИП было.

Да и изивините меня, я про этот случай не мог найти в справочниках по Высшей Математике, что вы хотите от учебника?

странно...

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

простите, а вы возведение в степень проходили? понятия "экспонента", "логарифм" вам знакомы? Или это тоже вычеркнули? Я про действительные числа. Или вы болели когда это все изучали? Просто ваша задача сводится к возведению e в дробную степень. И надо оно на практике тогда, когда степень дробная. Если степень целая, и небольшая, то никакой необходимости считать экспоненту просто нет - достаточно просто перемножить числа, что естественно проще и быстрее. В случае комплексного числа всё точно также, умножать сложнее немного, но вычислять экспоненту тем более не сахар - как вы заметили, арктангенс функция многозначная, и вычисляется компьютером только частично, с точностью +-PI/2, ЕМНИП. Это тоже было в школьной программе. Теперь этого нет?

залог успеха лежит в atan2, почти всегда нужен он, а не atan.

не. Надо было просто вспомнить, что арктангенс - многозначная функция, и для каждого аргумента имеет бесконечно много результатов. А следовательно библиотечной функции тут недостаточно - надо ещё и выбрать нужное значение из бесконечного множества.

[Фантом]

простите, а вы возведение в степень проходили? понятия "экспонента", "логарифм" вам знакомы? Или это тоже вычеркнули? Я про действительные числа. Или вы болели когда это все изучали?

В рамках нынешней школьной программы логарифмы проходят во втором полугодии 11 класса. Т.е. если BratSinot школьник, то он сейчас это вполне может не знать. С экспонентой, соответственно, ситуация аналогична. Обратные тригонометрические функции в программе упоминаются, такие "мелочи", как многозначность, в программу также не входят. Так что не стоит особенно возмущаться знаниями топикстартера.

[BratSinot]

простите, а вы возведение в степень проходили? понятия "экспонента", "логарифм" вам знакомы? Или это тоже вычеркнули? Я про действительные числа. Или вы болели когда это все изучали? Просто ваша задача сводится к возведению e в дробную степень. И надо оно на практике тогда, когда степень дробная. Если степень целая, и небольшая, то никакой необходимости считать экспоненту просто нет - достаточно просто перемножить числа, что естественно проще и быстрее. В случае комплексного числа всё точно также, умножать сложнее немного, но вычислять экспоненту тем более не сахар - как вы заметили, арктангенс функция многозначная, и вычисляется компьютером только частично, с точностью +-PI/2, ЕМНИП. Это тоже было в школьной программе. Теперь этого нет?

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

как вы заметили, арктангенс функция многозначная, и вычисляется компьютером только частично, с точностью +-PI/2, ЕМНИП. Это тоже было в школьной программе. Теперь этого нет?

не. Надо было просто вспомнить, что арктангенс - многозначная функция, и для каждого аргумента имеет бесконечно много результатов. А следовательно библиотечной функции тут недостаточно - надо ещё и выбрать нужное значение из бесконечного множества.

Эм... Я может и ошибаюсь, но atan() вроде не многозначная и лежит в [-pi/2; pi/2].

В рамках нынешней школьной программы логарифмы проходят во втором полугодии 11 класса.

Кто как, мы логарифмы прошли в 9-ом, из моих знакомых в 9-ом производные проходили, но не важно.

Так что не стоит особенно возмущаться знаниями топикстартера.

Я бы еще попросил в случае чего не сильно пинать. Если я что-то и не знаю, то не из-за того, что я "оболтус". Я в случае чего по справочникам "шастаю". Просто в этом случае, я с комплексными числами привык работать в алгебраической форме, поэтому и не подумал о том, какие подводные камни могут встретится.
Просто я при генерации еще проверял формулы "попроще", где все нормально было.


А теперь по теме, есть ли еще какие способы возведения в степень? И может ли формула Муавра дать выигрыш в скорости (скажем, при больших степенях)? Можно попробовать Бином Ньютона использовать, но что-то мне подсказывает что он будет медленнее даже формулы Муавра.

[drBatty]

Фантом
эх... вот оно как...

Эм... Я может и ошибаюсь, но atan() вроде не многозначная и лежит в [-pi/2; pi/2].

1. Ну извиняюсь тогда... Вы не виноваты, в таком "чудном" решении... Наверное вместо математики ввели Over9000 разной ненужной НЕХ типа работы с МСО...
2. функция atan() однозначная по определению. К сожалению, эта функция НЕ является обратной к tg(). Тут та-же беда, что с квадратными уравнениями - операция возведения в квадрат вполне однозначна, а обратная к ней - не однозначна. Это конечно решается определением, вот только при операциях с корнями нужно всегда помнить о том, что определение sqrt() является искусственным.
3. насколько я понимаю, если вы возводите в целую степень, то многозначность уничтожается. (кстати, любимый мой случай, когда НЕХ*0), и вы приходите к формулам, которые могли-бы получить и так, ручками, просто перемножая многочлены.
т.е. вместо того, что-бы считать (x+i*y)^2 вы переводите в r^2*exp(i*ph*2), что в общем-то одно и тоже, но считать проще, если возводить в дробную степень. Проблема в том, что вы вычисляете сначала ph = tg(y/x), и у вас тут всё нормально получается, формула однозначная, а потом вычисляете (y/x)^2=atg(ph*2) и вот здесь-то и проявляется неоднозначность - если ph < pi/4, то всё хорошо, а вот если ph > pi/4, то при умножении на 2 получается > pi/2, мы перепрыгиваем через НЕХ в точки pi/2 и попадаем на соседнюю тангенсоиду. А когда считаем atn() она об этом уже забыла, и считает, что мы были на центральной, нулевой тангенсоиде. Всего-то делов, добавить к результату единицу, и получить правильный арктангенс. Фишка в том, что при вычислении тангенса мы выходим из диапазона, а потом нам этот выход необходимо учесть. Точно тоже самое и с вашими синусами/косинусами.

А теперь по теме, есть ли еще какие способы возведения в степень? И может ли формула Муавра дать выигрыш в скорости (скажем, при больших степенях)?

может. По сравнению с перемножением в цикле. Бином Ньютона очевидно быстрее (криво работает с дробями). Вот только не очень понятно, зачем возводить в большие степени комплексные числа?

[KiWi]

А теперь по теме, есть ли еще какие способы возведения в степень? И может ли формула Муавра дать выигрыш в скорости (скажем, при больших степенях)? Можно попробовать Бином Ньютона использовать, но что-то мне подсказывает что он будет медленнее даже формулы Муавра.

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

P.S.: некоторую погрешность можно увидеть на (1000 + i)^2 -- косинус равен 0.999998000002, а синус -- 0.0019999980000020004. И в точности правильному ответу результат равен не будет.

[Portnov]

P.S.: некоторую погрешность можно увидеть на (1000 + i)^2 -- косинус равен 0.999998000002, а синус -- 0.0019999980000020004

Не проверял, но приведённая погрешность очень похожа на обыкновенную погрешность при работе с float-ами (и не связанную ни с комплексными числами, ни с тригонометрией, ни с формулой Муавра).

[BratSinot]

может. По сравнению с перемножением в цикле. Бином Ньютона очевидно быстрее (криво работает с дробями).

Подождите, это как? При вычислении степени по Биному Ньютона надо же будет значительно больше операций сделать?

Вот только не очень понятно, зачем возводить в большие степени комплексные числа?

Просто было интересно. На самом деле больше 5-6 не использовал.

[KiWi]

Не проверял, но приведённая погрешность очень похожа на обыкновенную погрешность при работе с float-ами (и не связанную ни с комплексными числами, ни с тригонометрией, ни с формулой Муавра).

Да, это именно она и есть.
Просто в случае формулы Муавра -- мы обязаны работать с floatом, потому что нужно считать синус, косинус и арктангенс.

[BratSinot]

Не проверял, но приведённая погрешность очень похожа на обыкновенную погрешность при работе с float-ами (и не связанную ни с комплексными числами, ни с тригонометрией, ни с формулой Муавра).

Да, это именно она и есть.
Просто в случае формулы Муавра -- мы обязаны работать с floatом, потому что нужно считать синус, косинус и арктангенс.

Ну я работаю с двойной точностью (с каких-то пор AMD позволила с ними работать, ибо железно они у меня не поддерживаются) и не занимаюсь точными математическими вычислениями, поэтом точностью часто можно пренебречь. Это может сказаться только при высоком приближении.

[aumit]

Подождите, это как? При вычислении степени по Биному Ньютона надо же будет значительно больше операций сделать?

В случае натуральных чисел алгоритм быстрого возведения в степень сделает за логарифмическое число операций.
А биномиальные коэффициенты можно вычислить за квадратичное время используя методы динамического программирования.

[NickLion]

Подождите, это как? При вычислении степени по Биному Ньютона надо же будет значительно больше операций сделать?

В случае натуральных чисел алгоритм быстрого возведения в степень сделает за логарифмическое число операций.
А биномиальные коэффициенты можно вычислить за квадратичное время используя методы динамического программирования.

За квадрат? o_O
Я, наверное, комбинаторику забыл, но у меня за линейное время вполне нормально получилось:
(bash)

Код: Выделить всё

n=8;c=1;echo $c;for((i=0;i<n;i++));do ((c*=n-i,c/=i+1)); echo $c; done

Другое дело, что логарифм лучше даже линейного.

[drBatty]

Подождите, это как? При вычислении степени по Биному Ньютона надо же будет значительно больше операций сделать?

я посчитал, у меня получилось ~n*4 для простого цикла. а в биноме ньютона нужно одно деление и два умножения на n. Я не знаю, что будет быстрее на конкретном CPU.

А вот по формуле Муавра время пропорционально O(1) операций, вот только нам понадобятся очень точные операции + очень длинные числа. А вот насколько точные, и насколько длинные - мне не сообразить. Пока всё влезает в сопроцессор, будет быстрее. Только надо прямо на асме писать.

мы обязаны работать с floatом, потому что нужно считать синус, косинус и арктангенс.

а если мы ограничимся x86, разве тогда промежуточные результаты не смогут иметь размер 80 бит?

[NickLion]

я посчитал, у меня получилось ~n*4 для простого цикла. а в биноме ньютона нужно одно деление и два умножения на n. Я не знаю, что будет быстрее на конкретном CPU.

Возвести в степень можно так:

Код: Выделить всё

#include <stdio.h>

int main(void)
{
    double R = 0.707106781187, I = 0.707106781187;
    double bR, bI;
    double pR = 1, pI = 0;
    int n = 1000000;
    bR = R; bI = I; // pow=1
    for( ; n; n >>= 1 ) {
        double tR;
        if( n & 1 ) {
            tR = pR * bR - pI * bI;
            pI = pI * bR + pR * bI;
            pR = tR;
        }
        tR = bR * bR - bI * bI;
        bI = 2 * bR * bI;
        bR = tR;
    }
    printf( "%lg %lg\n", pR, pI );
    return 0;
}

Тут даже для 1000000 будет всего 20 итераций цикла.

[BratSinot]

Брррр, ребят мы про GPU говорим. Там даже простой цикл умножения комплексных чисел работает медленее.

[NickLion]

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

[drBatty]

NickLion
а что это за алгоритм?

[aumit]

За квадрат? o_O
Я, наверное, комбинаторику забыл, но у меня за линейное время вполне нормально получилось:
(bash)

Код: Выделить всё

n=8;c=1;echo $c;for((i=0;i<n;i++));do ((c*=n-i,c/=i+1)); echo $c; done

Другое дело, что логарифм лучше даже линейного.

Молодец. Однако в биноме при биномиальных коэффициентах стоят степени того же комплексного числа.

[aumit]

Так, а почему мы именно Re проверяем на отрицательность? Почему Im не проверяем?

Потому что есть знаки тангенса и формулы приведения такие как tg(a+pi)=tg(a).

Там никогда не будет ровно 0, поэтому проблем не будет.

Во первых нужна писать корректные функции/макросы.
Во вторых сложно сказать будет или нет, т. к. точка может приближаться к вещественной оси.

[NickLion]

NickLion
а что это за алгоритм?

Названия я плохо запоминаю, если у него оно есть. Работает на основе двоичного представления числа. Допустим, если степень 100, то двоичное представление: 1100100, степени:

Код: Выделить всё

6 3 1
4 2 6 8 4 2 1
1 1 0 0 1 0 0

Тогда x¹⁰⁰=x⁴*x³²*x⁶⁴. А сами степени получаем возведением в квадрат.

За квадрат? o_O
Я, наверное, комбинаторику забыл, но у меня за линейное время вполне нормально получилось:
(bash)

Код: Выделить всё

n=8;c=1;echo $c;for((i=0;i<n;i++));do ((c*=n-i,c/=i+1)); echo $c; done

Другое дело, что логарифм лучше даже линейного.

Молодец. Однако в биноме при биномиальных коэффициентах стоят степени того же комплексного числа.

Стоп. Мы исходно имеем:
F=(R+I*i)ⁿ, отсюда
F=C_n⁰*R⁰*Iⁿ*iⁿ+C_n¹*R¹*I^n-1*i^n-1+...+C_nⁿ*Rⁿ*I⁰*i⁰.
где i^k раскрывается в зависимости от остатка деления на 4:
0 - +действительное
1 - +мнимое
2 - −действительное
3 - −мнимое
И всё. Другое дело, что тут точность может убежать очень и очень далеко - много операций, с возможно маленькими и большими числами.
Степени R и I каждый раз заново естественно считать не надо - можно домножить/разделить предыдущее. Опять же при делении возможна потеря точности можно немного улучшить, если пожертвовать массивчик размера N на убывающую степень.
Так что алгоритм линейный.

[BratSinot]

я посчитал, у меня получилось ~n*4 для простого цикла. а в биноме ньютона нужно одно деление и два умножения на n. Я не знаю, что будет быстрее на конкретном CPU.

Возвести в степень можно так:

Код: Выделить всё

#include <stdio.h>

int main(void)
{
    double R = 0.707106781187, I = 0.707106781187;
    double bR, bI;
    double pR = 1, pI = 0;
    int n = 1000000;
    bR = R; bI = I; // pow=1
    for( ; n; n >>= 1 ) {
        double tR;
        if( n & 1 ) {
            tR = pR * bR - pI * bI;
            pI = pI * bR + pR * bI;
            pR = tR;
        }
        tR = bR * bR - bI * bI;
        bI = 2 * bR * bI;
        bR = tR;
    }
    printf( "%lg %lg\n", pR, pI );
    return 0;
}

Тут даже для 1000000 будет всего 20 итераций цикла.

Что это за работающее чудо и откуда оно? o_O

Простите, забыл про сообщение выше.

В моем варианте это короче смотрится :

Код:

inline float2 cpow(float2 Z, unsigned int n) { float2 p=(float2)(1, 0); do { if(n&1) p=cmul(p, Z); Z=csqr(Z); }while(n>>=1); return p; }

[aumit]

так поиграться гладкая покраска, для двух фракталов.