### Author Topic: Efficient Computation of the Complex Exponential Function  (Read 7522 times)

0 Members and 1 Guest are viewing this topic.

#### Xeda112358

• they/them
• Moderator
• LV12 Extreme Poster (Next: 5000)
•            • • Posts: 4704
• Rating: +719/-6
• Calc-u-lator, do doo doo do do do. ##### Efficient Computation of the Complex Exponential Function
« on: December 03, 2014, 12:02:16 pm »
Overview
I was reading the BKM algorithm and I thought it was similar to my previous attempt at making an algorithm for the complex exponential function. However, I tried implementing it and failed miserably. So I altered my previous algorithm with an idea from the BKM algorithm, and I made my own algo that works as quickly, uses a smaller LUT, works on a wider range of values, and has a more straightforward implementation.

Why is the complex exponential function useful?
The complex exponential function can be evaluated as $e^{x+iy}=e^{x}\left(\cos(y)+i \sin(y)\right)$. So for people who actually need a complex exponential function, this is useful, but it can also compute the real exponential function(take y=0), or real sine and cosine, simultaneously (take x=0).

Goals of the algorithm
Although I only advertise my algorithm as working on the rectangle [-1,1]+[-1,1]i, it actually works on approximately (-1.754365792,1.754365792)+(-1.213427912,1.213427912)i. Feel free to use this to your advantage.
As well, since shifting and adding are easycitation needed, this is a shift and add algorithm. It takes approximately n iterations for n bits of accuracy, and 8 tables of n/2 floats (versus 18 tables suggested by the BKM algorithm).

The Algorithm
We will first give simple pseudo code, then proceed to optimize it.
As in my old algorithm and the BKM algorithm, we start by observing that ez-s=eze-s=ez/es. So then if z-a-b-c-d-...=0, then we have ez-a-b-c-...=ez/(eaebec...), so then 1=ez/(eaebec...), so eaebec...=ez.

The trick is finding values for a,b,c,... so that the multiplications are "easy" as in, comprised of a shift and add.

Since we are working with complex numbers, take the complex number d where the real and imaginary parts are in the set {-1,0,1}. Multiplication by those values is even more trivial than a shift and add. What we will do is take ea=1+d2-n. Multiplication by 2-n is a right shift by n.

If ea=1+d2-n, then log(ea)=log(1+d2-n), so a=log(1+d2-n). There are 9 possible values for d, so you might think, "oh no Zeda, we need 9 tables of complex numbers which is basically 18 tables of floats," but we all know that math+programming=magiccitation not needed. So let's look at how the complex logarithm is defined and call the real part of d, s, and the imaginary part t. So d=s+t*i. By the way, references to "log" imply the natural logarithm in the math world, so loge or ln() not log10.
$\log(x+iy) = \log(\sqrt(x^{2}+y^{2}))+i\tan^{-1}(y/x)= \frac{1}{2}\log(x^{2}+y^{2})+i\tan^{-1}(\frac{y}{x})$ (side note: atan has infinitely many valid outputs-- just add 2*k*pi for integers k, but we are only interested in k=0, the principle Arg).

So our tables would be with x=1+s2-n, y=t2-n.
Spoiler For tables for n:
s=-1, t=-1: $\frac{1}{2}\log((1-2^{-n})^{2}+2^{-2n})-i\tan^{-1}(\frac{2^{-n}}{1-2^{-n}})$
s=-1, t=0: $\log(1-2^{-n})$
s=-1, t=1: $\frac{1}{2}\log((1-2^{-n})^{2}+2^{-2n})+i\tan^{-1}(\frac{2^{-n}}{1-2^{-n}})$
s=0, t=-1: $\frac{1}{2}\log(1+2^{-2n})-i\tan^{-1}(2^{-n})$
s=0, t=0: 0
s=0, t=1: $\frac{1}{2}\log(1+2^{-2n})+i\tan^{-1}(2^{-n})$
s=1, t=-1: $\frac{1}{2}\log((1+2^{-n})^{2}+2^{-2n})-i\tan^{-1}(\frac{2^{-n}}{1+2^{-n}})$
s=1, t=0: $\log(1+2^{-n})$
s=1, t=1: $\frac{1}{2}\log((1+2^{-n})^{2}+2^{-2n})+i\tan^{-1}(\frac{2^{-n}}{1+2^{-n}})$
But there are lots of repeated values, so we actually need tables for:
Spoiler For Actual tables needed:
$\frac{1}{2}\log((1-2^{-n})^{2}+4^{-n})$
$\frac{1}{2}\log((1+2^{-n})^{2}+4^{-n})$
$\tan^{-1}(\frac{2^{-n}}{1-2^{-n}})=\tan^{-1}(\frac{1}{2^{n}-1})$
$\tan^{-1}(\frac{1}{2^{n}+1})$
$\log(1-2^{-n})$
$\log(1+2^{-n})$
$\frac{1}{2}\log(1+4^{-n})$
$\tan^{-1}(2^{-n})$

For our purposes, we will note that to choose s and t, we will use s=sign(iPart(Re(z)*2m+1)), t=sign(iPart(Im(z)*2m+1))  where iPart() returns the integer part of the input, and sign() returns 1 if positive, -1 if negative, and 0 if 0. Without further ado, the algorithm:
Code: [Select]
1.0→v1.0→wFor k,0,n-1    sign(x)→s    sign(y)→t    if |x|<2^(1-k)        0→s    if |y|<2^(1-k)        0→t    if s≠0 & t≠0        v→a        if t=0            x+log(1+s2^-k)→x            v+s*v>>k→v        if s=0            x+log(1+4^-k)→x            y+t*atan(2^-k)→y            v-tw>>k→v            w+ta>>k→w        if s=-1            x+log((1-2^-k)^2+4^-k)→x            y+t*atan(1/(2^k-1))→y            v-v>>k-tw>>k→w            w-w>>k+ta>>k→v        if s=1            x+log((1+2^-k)^2+4^-k)→x            y+t*atan(1/(2^k+1))→y            v+v>>k-tw>>k→w            w+w>>k+ta>>k→v
The result is v+iw to n bits of accuracy. It seems complicated, but it is basically just taking care of each case for s and t. As well, we can multiply x, y by 2 each iteration so instead of checking the nth bit, we can check if they are less than 1. Then we would need to multiply each table entry by 2k+1. For further optimization, we note that to 2m bits of accuracy, when k>m-1, log(1+d2-(k+1))=.5*log(1+d2-k) So we can basically reuse the last table value, divide by 2 each iteration, then multiplying by 2 to get our 2-(k+1). But wait, x/2*2=x, so we can directly use the last value in each table!
To clean things up a bit, lets name the tables and use an index to address each value:
Spoiler For tables:

for n from 0 to m-1:
table0=$2^{n}\log((1-2^{-n})^{2}+4^{-n})$
table1=$2^{n+1}\frac{1}{2}\log((1+2^{-n})^{2}+4^{-n})$
table2=$2^{n+1}\tan^{-1}(\frac{2^{-n}}{1-2^{-n}})=2^{n+1}\tan^{-1}(\frac{1}{2^{n}-1})$
table3=$2^{n+1}\tan^{-1}(\frac{1}{2^{n}+1})$
table4=$2^{n+1}\log(1-2^{-n})$
table5=$2^{n+1}\log(1+2^{-n})$
table6=$2^{n}\log(1+4^{-n})$
table7=$2^{n+1}\tan^{-1}(2^{-n})$
For the optimized code to 2m bits of accuracy,
Code: [Select]
1.0→v1.0→wFor k,0,2m-1    if k<m        k→n    x<<1→x    y<<1→y    sign(x)→s    sign(y)→t    if |x|<1        0→s    if |y|<1        0→t    v→a    if s=0        if t=-1            x+table6[n]→x            y-table7[n]→y            v+w>>k→v            w-a>>k→w        if t=1            x+table6[n]→x            y+table7[n]→y            v-w>>k→w            w+a>>k→v    if s=-1        if t=-1            x+table0[n]→x            y-table2[n]→y            v-v>>k+w>>k→w            w-w>>k-a>>k→v        if t=0            x+table4[n]→x            v-v>>k→w            w-w>>k→v        if t=1            x+table0[n]→x            y+table2[n]→y            v-v>>k-w>>k→w            w-w>>k+a>>k→v    if s=1        if t=-1            x+table1[n]→x            y-table3[n]→y            v+v>>k+w>>k→w            w+w>>k-a>>k→v        if t=0            x+table1[n]→x            v+v>>k→w            w+w>>k→v        if t=1            x+table1[n]→x            y+table3[n]→y            v+v>>k-w>>k→w            w+w>>k+a>>k→vSummary:
For floating point math with floats in base 2, arbitrary shifting is trivial (just change the exponent). So each iteration requires at most 6 trivial shifts, 6 trivial compares, 6 non-trivial additions/subtractions. So for 2m bits of precison, the total speed cost is 12mA+epsilon, where A is the speed of an addition or subtraction, and epsilon is a small, positive number. The total space cost is 8m Floats. Compared to the BKM algorithm, the speed is the same (with just a change in epsilon), but the BKM uses 18m Floats, so we save 10 Floats. We also work on a larger space than the BKM algorithm which works on a trapezoid region of the complex pain, entirely contained in the set (a proper subset) on which this algorithm works.

Future Plans
I plan to write this up more formally, including the proofs as to why we can use m elements for 2m bits of accuracy, why we can just look at the sign and the nth bits of the real and complex component of the input to choose our d. This last bit was a difficult task for the creator of the BKM algorithm and in their paper, they mentioned that the bounds were theoretical and only as tight as they could manage. In practice they found that the bounds were loose and they suspected that they could expand their bounds, but they couldn't prove it. I found an easy way to prove this for my algorithm which I suspect will apply to the BKM algorithm since at this point they are nearly identical.

Gosh I hope I don't have typos, formatting issues, or errors in the algorithms I posted. I have been typing and revising this for hours .__.
« Last Edit: December 03, 2014, 12:12:08 pm by Xeda112358 »

#### TheCoder1998

• LV6 Super Member (Next: 500)
•      • • Posts: 434
• Rating: +20/-2
• my art is written in code, not in graphite ##### Re: Efficient Computation of the Complex Exponential Function
« Reply #1 on: December 03, 2014, 12:15:44 pm »
This is brilliant!
I really enjoy reading your posts   my ticalc acc:

http://www.ticalc.org/archives/files/authors/113/11365.html

Spoiler For The Best Song Ever:

follow me on tumblr www.rickdepizza.tumblr.com

check out my anilist http://anilist.co/animelist/29701/Rickdepizza

#### Xeda112358

• they/them
• Moderator
• LV12 Extreme Poster (Next: 5000)
•            • • Posts: 4704
• Rating: +719/-6
• Calc-u-lator, do doo doo do do do. ##### Re: Efficient Computation of the Complex Exponential Function
« Reply #2 on: December 03, 2014, 12:30:19 pm »
Thanks! Next I plan to make a similar algorithm for the complex logarithm using the same LUTs. It's pretty similar.

EDIT: Complex log will offer a logarithm function (real and complex) as well as an arctangent function (real).

#### Xeda112358

• they/them
• Moderator
• LV12 Extreme Poster (Next: 5000)
•            • • Posts: 4704
• Rating: +719/-6
• Calc-u-lator, do doo doo do do do. ##### Re: Efficient Computation of the Complex Exponential Function
« Reply #3 on: December 08, 2014, 01:52:53 pm »
So for the actual, cleaned up code:

First, take the following tables:
$L1[k]:=2^{k}\log((1-2^{-k})^{2}+4^{-k})$
$L2[k]:=2^{k+1}\log(1-2^{-k})$
$L3[k]:=2^{k}\log(1+4^{-k})$
$L4[k]:=2^{k}\log((1+2^{-k})^{2}+4^{-k})$
$L5[k]:=2^{k+1}\log(1+2^{-k})$
$A1[k]:=2^{k+1}\arctan(\frac{1}{2^{k}-1})$
$A2[k]:=2^{k+1}\arctan(2^{-k})$
$A3[k]:=2^{k+1}\arctan(\frac{1}{2^{k}+1})$
Code: [Select]
1→v0→w    for k,0,2m-1:        if k<m:k→n        2x→x        2y→y        sign(x)→s        sign(y)→t        if |x|<1:0→s        if |y|<1:0→t        v→a        if s=-1            if t=-1                x-L1[n]→x                y+A1[n]→y                v-v>>k+w>>k→v                w-w>>k-a>>k→w            if t=0                x-L2[n]→x                v-v>>k→v                w-w>>k→w            if t=1                x-L1[n]→x                y-A1[n]→y                v-v>>k-w>>k→v                w-w>>k+a>>k→w        if s=0            if t=-1                x-L3[n]→x                y+A2[n]→y                v+w>>k→v                w-a>>k→w            if t=1                x-L3[n]→x                y-A2[n]→y                v-w>>k→v                w+a>>k→w        if s=1            if t=-1                x-L4[n]→x                y+A3[n]→y                v+v>>k+w>>k→v                w+w>>k-a>>k→w            if t=0                x-L5[n]→x                v+v>>k→v                w+w>>k→w            if t=1                x-L4[n]→x                y-A3[n]→y                v+v>>k-w>>k→v                w+w>>k+a>>k→wIt might also be convenient to move all the v,w calculations outside of the cases, since they are all essentially the same with components multiplied by s or t which are just multiplications by {1,0,-1}.
Code: [Select]
1→v0→w    for k,0,2m-1:        if k<m:k→n        2x→x        2y→y        sign(x)→s        sign(y)→t        if |x|<1:0→s        if |y|<1:0→t        v→a        v+sv>>k-tw>>k→v        w+sw>>k+ta>>k→w        if s=-1            if t=-1                x-L1[n]→x                y+A1[n]→y            if t=0                x-L2[n]→x            if t=1                x-L1[n]→x                y-A1[n]→y        if s=0            if t=-1                x-L3[n]→x                y+A2[n]→y            if t=1                x-L3[n]→x                y-A2[n]→y        if s=1            if t=-1                x-L4[n]→x                y+A3[n]→y            if t=0                x-L5[n]→x            if t=1                x-L4[n]→x                y-A3[n]→yAlso, in the last m iterations, can be even more optimized if you are using a binary floating point format for x and y. Convert the mantissas so that you have m*2^0 (so 1.1011*2^-5 would have a mantissa of 110110000.... but would be converted to 000011011000...). Now use these as your new, converted x and y value, treating them as ints. If the function "MSb(x)" returns the most significant bit in x (so for 64-bit numbers, bit 63):
Code: [Select]
sign(x)→bsign(y)→cconvert(|x|)→xconvert(|y|)→y    for k,m,2m-1:        MSb(x)→s        MSb(y)→t        x<<1→x        y<<1→y        v→a        v+s*b*v>>k-t*c*w>>k→v        w+s*b*w>>k+t*c*a>>k→wIn the last part, s,t are either 0 or 1, and b,c are either -1, 0, or 1. If any of them are 0, of course the algorithm can be optimized further (just plug in 0 to see how it reduces). For example, suppose x=0 upon entering:
Code: [Select]
    for k,m,2m-1:        MSb(y)→t        y<<1→y        v→a        v-t*c*w>>k→v        w+t*c*a>>k→w
« Last Edit: December 08, 2014, 02:02:26 pm by Xeda112358 »

#### ivelegacy ##### Re: Efficient Computation of the Complex Exponential Function
« Reply #4 on: September 02, 2016, 10:10:02 am »
hi Xeda112358
have you already implemented this algorithm in C/Matlab ?
here I have already implemented CORDIC.sin-cos in VHDL
and tested on RTL simulator and fpga

let me know « Last Edit: September 02, 2016, 10:16:26 am by ivelegacy »

#### Xeda112358

• they/them
• Moderator
• LV12 Extreme Poster (Next: 5000)
•            • • Posts: 4704
• Rating: +719/-6
• Calc-u-lator, do doo doo do do do. ##### Re: Efficient Computation of the Complex Exponential Function
« Reply #5 on: September 03, 2016, 02:50:39 pm »
hi Xeda112358
have you already implemented this algorithm in C/Matlab ?
here I have already implemented CORDIC.sin-cos in VHDL
and tested on RTL simulator and fpga

let me know I have not. Are you the person who emailed me? I was about to point you here since I'm too exhausted to work on this at the moment (worked 16+ hours in a 24 hour span of time).