@Nanowar and tr1p1ea: Wow, those are amazing o.o

GCD is the Greatest Common Divisor. So GCD(15,27) is 3 since 3 is the largest number that divides both 27 and 15. It isn't used often, bt if somebody wanted to make a 16-bit rational library *cough*I should*cough* it is extremely useful.

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.

As it is only a 16-bit by 16-bit multiplication, yes, only up to 16-bit factors

And the upper 16 bits of the result are in DE, lower 16 in BC.

No, the eZ80 has two modes, ADL mode (24-bit) and Z80 (16-bit). In Z80 mode, certain instructions are one or two cycles faster, so it is beneficial for me to use that mode since I don't use any ADL-specific instructions. As well, add hl,hl sets the c flag on 16-bit overflow in Z80 mode, which I needed, but doesn't set it in ADL mode (it has to overflow 24 bits). Since there is no easy way to access the top bits of HL, it would have to be performed in Z80 mode, making it take even more clock cycles.


The full 32-bit result is returned.

So do you have the coordinates of the objects?

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).