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 easy^{citation 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 e^z-s=e^ze^-s=e^z/e^s. So then if z-a-b-c-d-...=0, then we have e^z-a-b-c-...=e^z/(e^ae^be^c...), so then 1=e^z/(e^ae^be^c...), so e^ae^be^c...=e^z.

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 e^a=1+d2^-n. Multiplication by 2^-n is a right shift by n.

If e^a=1+d2^-n, then log(e^a)=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=magic^{citation 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 log_e or ln() not log₁₀.
##\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.
But there are lots of repeated values, so we actually need tables for:

For our purposes, we will note that to choose s and t, we will use s=sign(iPart(Re(z)*2^m+1)), t=sign(iPart(Im(z)*2^m+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:
1.0→v
1.0→w
For 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 2^k+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:
For the optimized code to 2m bits of accuracy,
1.0→v
1.0→w
For 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→v
Summary:
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 .__.

Good job! Do you think this can have applications in z80 programming (is it faster, more memory efficient, etc.)?

This would not be the best way to compute tangent (though it is decent). For 21 bits of accuracy, there are probably much better ways It's just neat and if you've never seen Horner's method, it can get you thinking about math programming in the future.

I think that if you really want speed, you should use a LUT for the sine function, calculate the cosine as sine+90°, and calculate the tangent as sine/cosine. This does require a bit of memory and isn't really related to math though.

Related: I previously had no idea how the sine, cosine and tangent were actually calculated, and good job on simplifying it so much!

EDIT: while the z80 will probably not be able to calculate that at high speed for use in games, I guess it might be usefull for CAS-related programs.

Yeah, LUT's are the fastest method. There are fairly efficient ways to use small LUTs to compute the fun functions (trig and exponentials, among others). Glad to show you something new!

This is really cool! I found using ries that the constant is arctan(pi)/pi (if it turns out to be a constant, which seems to be the case).

EDIT: that is probably not correct. I just saw that first and it had a tangent in it. Looking at the others, 4/(pi^2) looks like it could be promising. But those are all approximations anyway.

Yeah, 4/pi² looks most promising, but it might not converge to anything nice. Then again, I have lots of identities that I've derived about the Bernoulli numbers, so who knows! For example, if you sum them up to get 'C', then do c/(c-1), you get euler's number, and if you sum the absolute value of the terms, you get a function involving cotangent!

Wow this topic is really interesting!
It's pretty funny that I'm reading this now, because it's Math-day at school today

Math day is best day.


I worked on this a little last night, and I will probably work on this some more. I might actually figure out the constant and I would be more inclined to believe that it has something to do with e.

Hi again! This time I was fooling around with tangent (who knows why? I don't .__.) Anyways, I was having some fun and I figured out that Horner's method can be abused to approximate tangent by a similar logic I used to approximate arctangent: all you have to do is approximate infinitely many terms of the Maclaurin series with one simple function!


So before the fun, I'll show what Horner's method is (it's a super simple idea):
If you've never seen the identity that ##sin(x)=x-\frac{x^{3}}{3!}+\frac{x^{5}}{5!}-\frac{x^{7}}{7!}+...##, that's fine, just trust me (you'll get to that in calculus and it's crazy cool stuff). Horner's method is like what a true 31337 programmer would come up with if you want to efficiently approximate that polynomial. First, factor out the first term in the series (##x##):
##\sin(x)=x(1-\frac{x^{2}}{3!}+\frac{x^{4}}{5!}-\frac{x^{6}}{7!}+...##
Now in that inside term, factor out the ##-\frac{x^{2}}{3!}##:

##\sin(x)=x(1-\frac{x^{2}}{3!}(1-\frac{x^{2}}{4\cdot 5}+\frac{x^{4}}{4\cdot 5\cdot 6\cdot 7\cdot }-\frac{x^{6}}{4\cdot 5\cdot 6\cdot 7\cdot 8\cdot 9}+...##
Rinse, repeat:

##\sin(x)=x(1-\frac{x^{2}}{3!}(1-\frac{x^{2}}{4\cdot 5}(1-\frac{x^{2}}{6\cdot 7}(1-\frac{x^{2}}{8\cdot 9}-...##
So the programmer sees, "aw, shweet, I only need to compute x² once and reuse it a bunch of times with some constant multiplications!" So yeah, Horner's method is snazztacular, and I urge you to use it when possible.


Now let's look at tan(x). The Maclaurin series expansion for tan(x) is uuuugly. Each term uses Bernoulli numbers which are recursively computed and are intricately linked with some really cool functions like the Riemann Zeta function. For comparison (thanks to Mathworld for the table):

##\cos(x)=\sum\limits_{k=0}^{\infty}{\frac{(-1)^{k}}{(2k)!}x^{2k}}## (simple, cute, compact)
##\sin(x)=\sum\limits_{k=0}^{\infty}{\frac{(-1)^{k}}{(2k+1)!}x^{2k+1}}## (simple, cute, compact)
##\tan(x)=\sum\limits_{k=0}^{\infty}{\frac{(-1)^{k}4^{k+1}\left(4^{k+1}-1\right)B_{2k+2}}{(2k+2)!}x^{2k+1}}## (wtf are math)
That's what it looks like to divide sine by cosine. Not pretty (well, okay, it is actually pretty gorgeous, but whatevs). So the first few terms of tan(x) are:
##\tan(x)=x+\frac{x^{3}}{3}+\frac{2x^{5}}{15}+\frac{17x^{7}}{315}+\frac{62x^{9}}{2835}+\frac{1382x^{11}}{155925}+...##
Nasty gorgeous, but applying Horner's Method, we end up with:
##\tan(x)=x(1+\frac{x^{2}}{3}(1+\frac{2x^{2}}{5}(1+\frac{17x^{2}}{42}(1+\frac{62x^{2}}{153}(1+...##

So seeing this, I realized, "wait a minute, those coefficients sure don't look like they are converging to zero or one, but something in between!" And that made me think about how ##\frac{1}{1-ax^{2}}=(1+ax^{2}(1+ax^{2}(1+ax^{2}(1+ax^{2}(1+...##. Now, I am a mathematician, but I am also a programmer so what I am about to say is blasphemous and I am sorry, but I decided to numerically evaluate and assume those constants converge. I ended up with approximately .405285. What I am saying is:

##\tan(x)\approx x(1+\frac{x^{2}}{3}(1+\frac{2x^{2}}{5}(1+\frac{17x^{2}}{42}(1+\frac{62x^{2}}{153}(1+.405285x^{2}(1+.405285x^{2}(1+.405285x^{2}(1+...##
So then I ended up with:

##\tan(x)\approx x(1+\frac{x^{2}}{3}(1+\frac{2x^{2}}{5}(1+\frac{17x^{2}}{42}(1+\frac{62x^{2}}{153}\frac{1}{1-.405285x^{2}}))))##
##\tan(x)\approx x(1+\frac{x^{2}}{3}(1+\frac{2x^{2}}{5}(1+\frac{17x^{2}}{42}(1+\frac{62x^{2}}{153-62x^{2}}))))##
##\tan(x)\approx x(1+\frac{x^{2}}{3}(1+\frac{2x^{2}}{5}(1+\frac{17x^{2}}{42}\frac{153}{153-62x^{2}})))##It's still a bit of work to compute, but I had fun! No to do the mathy thing and actually prove the constants converge and maybe even find it's exact value!

Wow, this is cool! Will it work for grayscale models?

Well, I've been wanting to get back into graph theory, so I guess this will make an excuse to look into some of my books.

At the store where I work, I can buy  pack of 8 tr1p1eas for either $1 or $1.50 (I forget). They are cheap store-brand ones, though.

Well I could do it all except the CSE stuff, but I doubt I'll ever again have time or drive