This is a work-in-progress, but the intent is for it to convert a token into it's individual chars. It first converts the token to ASCII (this is easy; the OS provides a call), then the hard part is converting the ASCII to individual tokens. Uppercase is easy, and I also added conversion for lowercase chars, space, and "(" so that should get most tokens. However, something like tan^-1( will not properly convert the ^-1.

I am trying to write some link port routines and I am having some issues. Every so often, data won't be correctly transmitted and I have deduced that this only occurs when the sent byte has bit 1,3, or 5 set, it sometimes causes the next bit up to be set (so if bit 3 is supposed to be set, it also sets bit 4). This is a chaotic occurence.

What I do know is when the odd bits are read, bit 0 of port 0 (so the tip) reads 1. As well, the sending calc had a low battery, so I am not sure if it is an issue of bringing the lines up. I am pretty sure this is not the problem, though, as I switched out batteries.

I also tried adding code that guaranteed that the lines were reading correctly, otherwise it would rewrite the value to port 0 until it did. This still didn't work.
I tried adding huge delays, it still was just as unreliable.
Interrupts and link assist are disabled.

Below is my code if you want to muddle through it. It is supposed to send a value in Ans to the other calc and the receiving calc receives in Ans (modulo 100).
#define bcall(x) rst 28h \ .dw x
.db $BB,$6D
.org $9D95
    bcall(4AD7h)
    bcall(4AEFh)
sendb:
    di
    ld c,a
    xor a
    out (8),a
    ld e,55h
sendloop:
    rrc c
    rla
    sla e
    ccf
    rla
    out (0),a
    ld b,6          ;this is a delay, and not claculated to be precise. We probably do not even need a delay
    djnz $
    jr nz,sendloop
    ld b,10         ;delay
    djnz $
    xor a
    out (0),a
    ret
#define bcall(x) rst 28h \ .dw x
.db $BB,$6D
.org $9D95
    call getb
    bcall(478Ch)
    bcall(4ABFh)
    ret
getb:
    di
    xor a
    out (8),a
    ld bc,$08AA
readloop:
    in a,(0)
    xor c
    rra
    jr c,readloop
    rra             ;bits cycled in are masked with 0x55 as a sideeffect. Need to invert anyways, so mask at the end with 0xAA
    rr c
    djnz readloop
    ld a,c
    and 1
    rrca
    xor c
    xor $AA
    ret


I don't want to put this in the Optimized Routines thread as I haven't had it long enough for it to be ultra-optimized.

This code is intended as a method of testing if HL is divisible by DE. It's average run time for random 16-bit inputs is 163cc, and where DE<=HL (more likely in practice), it is 229cc. For a code that is doing prime-checking, it will more likely take about 707.5cc (Since DE will be roughly the square root of HL or less).
isDivisible:
;;Inputs: DE,HL
;;Outputs: c flag set if HL is not divisible by DE, else c flag is set.
;;         HL is 0 if true.
;;1087cc worst case (HL=65535, DE=1).
;;Average for random inputs: 163cc
;;Average for DE<=HL: 229cc.
    ld a,d \ or e \ ccf \ ret z         ;remove this if DE is always guaranteed non-zero
;step 1
    ld a,e \ or l \ rra \ jr c,step2    ;\
    srl d \ rr e \ rr h \ rr l          ; |
    ld a,e \ or l \ rra \ jr nc,$-11    ; |Remove these if DE is always guaranteed odd at input.
step2:                                  ; |
    ld a,e \ rra \ ccf \ ret c          ;/
;steps 3, 4, and 5
    ld a,l
    or a
loop:
    sbc hl,de \ ret c \ ret z
    rr h \ rra \ bit 0,a \ jr z,$-5
    ld l,a
    jp loop

Motivation and Development
I often find myself in a situation where I need to find the factors of a number, but I have no technology around to aid me. This means I need to use... mental arithmetic! I've been doing this for 15 years, so I have refined my mental process quite a bit. It is still a trial division algorithm, but with a very obfuscated "division" technique. We don't need to do 1131/7 to see if it is divisible by 7, we just need to see if 7 divides 1131 and this is what my algorithm does. Interestingly, testing divisibility at the algorithmic level is a little faster than division. Not by much, but it is also non-negligible.
The Algorithm
The core algorithm is designed around checking that (A mod B == 0) is true or false. We also make the assumption that B is odd and by extension, non-zero. The case where B is non-zero and even will be discussed later.

Since B is odd, 2 does not divide B. This means that if A is even:
    (A mod B == 0) if and only if  (A/2 mod B)==0.
We also know by the definition of divisibility that
    (A mod B) == (A+c*B mod B)
where c is any integer. Combining all this, we have an algorithm:

• Remove all factors of 2 from A
• With A now odd, do A=A-B
  
  • If the result is zero, that means (A mod B == 0)
  • If the result underflow (becomes "negative", or on the Z80, sets the carry flag), it means that A was somewhere on [1,B-1], so it is not divisible by B.
• Continue back at 1.

Now suppose B is allowed to be non-zero and even. Then B is of the form d*2^k where d is odd. This just means there are some factors of 2 that can be removed from B until it is odd. The only way A is divisible by B, is if it has the same number or more of factors of 2 as B. If we factor out common factors of 2 and find B is still even, then A is not divisible by B. Otherwise we have an odd number and only need to check the new (A mod d) for which we can use the "odd algorithm" above.
So putting it all together:

• If B==0, return FALSE.
• Remove common factors of 2 from A and B.
• If B is even, return FALSE.
• Remove all factors of 2 from A.
• Subtract B from A (A=A-B).
  
  • If the result is zero, return TRUE.
    
  • If the result is "negative" (setting the carry flag on many processors), return FALSE.
• Repeat at 4.

The overhead steps are (1) to (3). The iterated steps are (4) and (5).
Because (5) always produces an even number, when it then performs step 4, it always divides by at least one factor of 2. This means the algorithm takes at most 1+ceil(log2(A))-floor(log2(B) iterations. For example, if A is a 37-bit number and B is a 13-bit number,this takes at most 38-13 = 25 iterations. However, in practice it is usually fewer iterations.
Example Time:
Say I wanted to test if 1337 is divisible by 17.
Since 17 is odd, we can proceed.
1337 is odd, so no factors of 2 to remove.
1337-17 == 1320.
1320/2 == 660
660/2 == 330
330/2 == 165
165-17 == 148
148/2 == 74
74/2 == 37
37-17 == 20
20/2 == 10
10/2 == 5
5-17 == -12
So 1337 is not divisible by 17.

Now test divisibility by 7:
1337 => 1330
=>665
=>658
=>329
=>322
=>161
=>154
=>77
=>70
=>35
=>28
=>14
=>7
=>0

So 1337 is divisible by 7.


The worst-case timing is 66*16+31 = 1087

A few years ago when Builderboy wrote his fire graphics turorial, I remember thinking it was super cool and I went on to make my own fire engines designed around it. However, one thing always peeved me-- the design was very clever and made for speed, but it seemed too artificial for my liking-- particularly in the way it masked out pixels.

The idea was to kill a pixel 1/8 of the time, so a PRNG was used to select from one of the following masks:
%01111111
%10111111
%11011111
%11101111
%11110111
%11111011
%11111101
%11111110
Of course an AND mask will kill one bit in a byte, so exactly 1/8 of the pixels in a byte. A few weeks ago, this was bugging me and preventing my sleep. I set out to devise an algorithm that would generate masks that didn't necessarily have a 1/8 kill rate, but over time the probability converges to 1/8. For example, one mask might be %11111111 (0% killed) and the next could be %10110111 (25% killed).

The Idea
Suppose you 3 random bits and you OR them together. The only way this can result in a 0 is if every input was a 0. The probability of that happening is .5^3 = .125 = 1/8. So for our fire mask, we can generate three random bytes, OR them together, and each bit has a 1/8 chance of being reset.
Speed
The disadvantage here is in speed. Builderboy's method requires one pseudo-random number, this method requires three. However, if we consider that we are almost certainly going to use a fast pseudo-random number generator, and that we will want a long (ish) period, we have room to take advantage of the PRNG and achieve almost the same speed. For example, suppose you generate a 24-bit pseudo-random number-- with this method, you can just OR the three bytes generated (12cc) versus using an LUT (Builder's method):
ld hl,LUT
and 7
add a,l
ld l,a
jr nc,$+3
inc hl
ld a,(hl)
;43cc
In the example code I will give, I use a 16-bit PRNG, so I generate three of these numbers (6 8-bit values) for 2 masks, making generation take 1.5 times longer than Builder's as opposed to 3 times as long.
Considerations
In order for this to work, you need an RNG or PRNG in which every bit has a 50% chance of being set or reset. I went with a Lehmer LCG that was fast to compute and had a period of 2^16 and had this property.
Example
The following code works at 6MHz or 15MHz and the LCD will provide almost no bottleneck. Sorry, no screenie:
smc = 1    ;1 for SMC, 0 for no SMC (use 1 if code is in RAM; it's faster)
plotSScreen = 9340h
saveSScreen = 86ECh


;==============================
#IF smc = 0
seed = saveSScreen
#ENDIF
;==============================



.db $BB,$6D
.org $9D95
fire:
;;first, set up. We will be writing bytes to the LCD left to right then down
    di
    ld a,7      ;LCD y-increment
    out (16),a
;;setup the keyboard port to read keys [Enter]...[Clear]
    ld a,%11111101
    out (1),a
;make the bottom row of pixel;s black to feed the flames
    ld hl,plotSScreen+756
    ld bc,$0CFF
    ld (hl),c \ inc hl \ djnz $-2
fireloopmain:
    ld ix,plotSScreen+12
    in a,(1)
    and %01000000
    ret z
    call LCG
    ld b,63
fireloop:
;wait for LCD delay
    in a,(16)
    rla
    jr c,$-3
    ld a,80h+63
    sub b
    out (16),a
    push bc
    call LCG
    ld a,20h
    out (16),a
    call fire_2bytes+3
    call fire_2bytes
    call fire_2bytes
    call fire_2bytes
    call fire_2bytes
    call fire_2bytes
    pop bc
    djnz fireloop
    jp fireloopmain   
fire_2bytes:
    call lcg
    push hl
    call lcg
    pop de
    ld a,e
    or d
    or l
    and (ix)
    out (17),a
    ld (ix-12),a
    inc ix
    push hl
    call lcg
    pop af
    or h
    or l
    and (ix)
    out (17),a
    ld (ix-12),a
    inc ix
    ret
lcg:
;240cc or 241cc, condition dependent (-6cc using SMC)
;;uses the Lehmer RNG used by the Sinclair ZX81
#IF SMC = 1
seed = $+1
    ld hl,1
#ELSE
    ld hl,(seed)
#ENDIF
;multiply by 75
    ld c,l
    ld b,h
    xor a
    ld d,a
    adc hl,hl
    jr nz,$+9
    ld hl,-74
    ld (seed),hl
    ret
    rla
    add hl,hl \ rla
    add hl,hl \ rla \ add hl,bc \ adc a,d
    add hl,hl \ rla
    add hl,hl \ rla \ add hl,bc \ adc a,d
    add hl,hl \ rla \ add hl,bc \ adc a,d
;mod by 2^16 + 1 (a prime)
;current form is A*2^16+HL
;need:
;  (A*2^16+HL) mod (2^16+1)
;add 0 as +1-1
;  (A*(2^16+1-1)+HL) mod (2^16+1)
;distribute
;  (A*(2^16+1)-A+HL) mod (2^16+1)
;A*(2^16+1) mod 2^16+1 = 0, so remove
;  (-A+HL) mod (2^16+1)
;Oh hey, that's easy! :P
;I use this trick everywhere, you should, too.
    ld e,a
    sbc hl,de       ;No need to reset the c flag since it is already
    jr nc,$+3
    inc hl
    ld (seed),hl
    ret
EDIT: .gif attatched. It looks slower in the screenshot than my actual calc, though.
EDIT2: On my actual calc, it is roughly 17.2FPS

I'm really peeved that I cannot seem to find an algorithm for computing arcsine that works like what I have. This algorithm is based on my work from years ago to compute sine and I have no idea why I never reversed it before now. Anyways, the algorithm:
ArcSine(x), ##x\in[-1,1]##
    x=2x
    s=z=sign(x)
iterate n times
    x=x*x-2
    if x<0
        s=-s
    z=2z+s
return pi*z*2^(-n-2)
This algorithm extracts one bit each iteration, so for 16 bit of accuracy, you would iterate 16 times. I think it is a pretty cute and compact algorithm, so where is it? @_@

As a note, at the endpoints {-1,1} it may not give the expected answer. In both cases, if allowed to run infinitely, it would return (in binary): ##\frac{\pi}{2}.1111111111111..._{2}## which of course is equivalent to ##\frac{\pi}{2}##.

Hi all, I was writing a binary search algorithm and since I cannot test it just yet, I was hoping people might be able to spot any problems or optimizations. My intent was to optimize this for the eZ80. Feel free to offer suggestions that include instructions from the eZ80 instruction set.

First, the input is a "key" to search for. It is zero terminated, so it may be something like .db "Hello",0
The table to look for it in is comprised of indices to the actual data. So the table might look like:
table_start:
.dw (table_end-table_start)/2-1
.dw expr0
.dw expr1
.dw expr2
...
.dw exprlast
.table_end:


expr0: .db "AAAaaa",0
;other data
expr1: .db "BBBbbb",0
;other data
...
The actual strings do not need to be in order, just the order in which they are indexed. The strings must be zero-terminated, too.

So, passing a pointer to the keyword in DE and the table start in HL:
BinarySearch:
;expected in Z80 mode
    ld (key_ptr),de
    ld e,(hl) \ inc hl
    ld d,(hl) \ inc hl
    ld b,h \ ld c,l
    ld (max),de
    or a \ sbc hl,hl \ ld (min),hl
.loop0:
;average max (de) and min (hl) with truncation
    add hl,de \ rr h \ rr l
    push hl
;2 bytes per index, and BC points to the base of the table
    add hl,hl \ add hl,bc
;load the pointer to the expression to compare to in HL
    ld hl,(hl) \ ld de,(key_ptr)
.loop1:
    ld a,(de) \ cp (hl) \ jr c,less
;if they are the same, make sure to test if the bytes are zero (termination)
    jr z,$+4
    jr nc,more
;if a != 0 then the string is not terminated, but we don't know which string is "bigger" so keep looping
    inc de \ inc hl \ or a \ jr nz,loop1
.match:
    pop bc \ ret
.less:
; "or a" is just to reset the c flag for future code
;if [index]<key, add 1 to the average of min and max and set that as the new min
    or a \ pop hl \ inc hl \ ld (min),hl \ ld de,(max) \ jr loop0_end
.more:
;[index]>key set max to the average of max and min.
    pop de \ ld (max),de \ ld hl,(min)
.loop0_end:
    sbc hl,de \ add hl,de \ jr nz,loop0
    or 1
    ret
The output should have z flag if a match was found, else nz. In the case of a match being found, BC is the index number, DE points to the byte after the key, HL points to the byte after the match. If no match was found, HL and DE point to the index in front of where the match should have been and BC points to the table data start. In any case, the c flag is reset.

EDIT: See the necro-update below-- it provides a much more efficient algorithm, even 1 multiply better than most suggested implementations of Newton-Raphson division

---

In the study of numerical algorithms, quadratic convergence refers to an iterative algorithm that approximates a function and each iteration doubles the digits of accuracy. In this post, I will expose an algorithm for division that can double, triple, quadruple, etc. the number of digits of accuracy at each iteration. I will then show that the optimal algorithm is that which offers cubic convergence (digits of accuracy tripled).

What does it mean to multiply the number of correct digits by a?
Suppose that we have a sequence ##x_{0},x_{1},x_{2},x_{3},...## and suppose that ##x_{n}\rightarrow c## Then by the definition of convergence, ##\lim\limits_{n\rightarrow\infty}{|x_{n}-c|}=0##. You should always think of |x-y| as a "distance" between two points, and in this case, that distance is the error of our approximations, ##x_{n}##.

If the number of correct digits is multiplied by a, then that means ##\log(|x_{n+1}-c|)=a\cdot\log|x_{n}-c|##, or rewritten using logarithm tricks, ##|x_{n+1}-c|=|x_{n}-c|^{a}##. This assumes that the error is initially less than 1.

The Algorithm
The algorithm is very simple, but to break it down we need to look at "obvious" bits of math. Firstly, ##\frac{N}{D}=\frac{N\cdot c}{D\cdot c}##, and secondly, all real numbers can be written in the form of ##x=y\cdot 2^{m}, y\in(.5,1]##. You may have said "duh" at the first statement, but if you need convincing of the second, it basically says that if you divide or multiply a number by 2 enough times, then it will eventually get between 1/2 and 1. Or in another way, shown by example: 2¹¹<3767<2¹², so dividing it all by 2¹², .5<3767/4096<1. All real numbers can be bounded by consecutive power of 2, therefore they can all be written in the above form.

The final piece of the puzzle is to recognize that ##\frac{N}{1}={N}##. What we will do is recursively multiply D by some constant c at each iteration in order to drive it closer to 1. If we use range reduction techniques to get D on (.5,1] (multiply both N and D by the power of 2 that gets D in that range). Then what we want is to choose c so that ##|1-D_{n}\cdot c|= |1-D_{n}|^{a}##. If ##0\leq D_{n}\leq 1##, then we have ##1-D_{n}\cdot c = (1-D_{n})^{a}##. When a is a natural number greater than 1 (if it is 1, then there is no convergence ), we have the following:
##1-D_{n}\cdot c = \sum\limits_{k=0}^{a}{{a \choose k}(-D_{n})^{k}}##
##D_{n}\cdot c = 1-\sum\limits_{k=0}^{a}{{a \choose k}(-D_{n})^{k}}##
##D_{n}\cdot c = 1-(1+\sum\limits_{k=1}^{a}{{a \choose k}(-D_{n})^{k}})##
##D_{n}\cdot c = -\sum\limits_{k=1}^{a}{{a \choose k}(-D_{n})^{k}}##
##D_{n}\cdot c = {a \choose 1}D_{n}-{a \choose 2}D_{n}^{2}+{a \choose 3}D_{n}^{3}+\cdots##
##c = {a \choose 1}-{a \choose 2}D_{n}+{a \choose 3}D_{n}^{2}+\cdots##
Using a technique similar to Horner's Method, we can obtain:
##c = {a \choose 1}-D_{n}({a \choose 2}-D_{n}({a \choose 3}-D_{n}({a \choose 4}-D_{n}\cdots##
Then the selection for c requires a-2 multiplications and a-1 subtractions. If we take ##N_{n+1}=N_{n}\cdot c## and ##D_{n+1}=D_{n}\cdot c##, then the number of digits of accuracy is multiplied by a using a total of a multiplications and a-1.

Example
Lets say we want quadratic convergence. That is, a=2. Then ##c_{n+1}=2-D_{n}## That's right-- all you have to do is multiply by 2 minus the denominator at each iteration and if you started with 1 digit of accuracy, you next get 2 digits, then 4, then 8, 16, 32,.... Let's perform 49/39. Divide numerator and denominator by 2⁶=64 and we get .765625/.609375.

First iteration:
c=2-D=1.390625
Nc=1.064697266
Dc=.8474121094

Second iteration:
c=1.152587891
Nc=1.227157176
Dc=.9767169356

next iterations in the form {c,N,D}:
{1.023283064,1.255729155,.9994578989}
{1.000542101,1.256409887,.9999997061}
{1.000000294,1.256410256,.9999999999}
{1.000000294,1.256410256,1.000000000}

And in fact, 49/39 is 1.256410256 according to my calc.

Can we find an optimal choice for a?
The only variable over which we have control is a-- how quickly the algorithm converges. Therefore we should find an optimal choice of a. In this case, "optimal" means getting the most accuracy for the least work (computation time).

First, ##1-D_{0}<.5=2^{-1}## So then after n iterations, we have at least aⁿ bits of accuracy and since ##D_{0}## may get arbitrarily close to .5, this is the best lower bound on the number of bits of accuracy achieved for n iterations of a randomly chosen D. As well, at n iterations, we need to use ##a\cdot n## multiplications and the subtractions are treated as trivial. Then in order to get m bits of accuracy, we need roughly ##a\log_{a}(m)## multiplications. We want to minimize this function of a for a is an integer greater than 1. To do that, we find when ##0=\frac{d}{da}(a\log_{a}(m))##:
##0=\frac{d}{da}(a\log_{a}(m))##
##0=\frac{d}{da}(a\frac{\log(m)}{\log(a)})##
##0=\frac{d}{da}(\frac{a}{\log(a)})##
##0=\frac{1}{\log(a)}-\frac{1}{\log^{2}(a)}##
Since a>1:
##0=1-\frac{1}{\log(a)}##
##1=\frac{1}{\log(a)}##
##\log(a)=1##
##a=e##

However, e is not an integer, but it is bounded by 2 and 3, so check which is smaller: 2/log(2) or 3/log(3) and we find that a=3 is the optimal value that achieves the most accuracy for the least work.

Example:
Our algorithm body is:
    c=3-D*(3-D)
    N=N*c
    D=D*c
Then 49/39=.765625/.609375 and {c,N,D}=:
{1.543212891,1.181522369,.9403953552}
{1.063157358,1.256144201,.9997882418}
{1.000211803,1.256410256,1.000000000}


Conclusion
In reality, their is a specific M so that for all m>M, a=3 always requires less computational power than a=2. Before that cutoff, there are values for which they require the same amount of work or even the cost is in favor of a=2.  For example, 32768 bits of accuracy can be achieved with the same amount of work. The cutoff appears to be m=2^24 which is pretty gosh darn huge (16777216 bits of accuracy is over 5 million digits).

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!