So I've been developing a way to compute the geometric mean of two numbers (defined by ##\sqrt(xy)##) and my original algorithm had a few snazzy properties:
-It is an add-shift algorithm
-No multiplication needed and so no risk of overflow
-No explicit square root operation performed

I played with the math a bit, and if I take it out of the realm of add-shift routines, and allow division and multiplication, I came up with an algorithm that is relatively easy to do by hand, converges quarticly (quadruples the digits of accuracy at each step), and still has no overflow risk. In practice, on my TI-84+, I get full 14-digit precision in 2 iterations. The next iteration should theoretically get 56 digits, and each iteration requires 2 divisions, 1 inverse, and 1 multiplication (plus a shift and a few adds). Here is an example to show how it works:

Say we wish to find the geometric mean of two numbers, x=209, y=654321. First, for notation, I am just going to refer to the geometric mean as {x,y}, so {209,654321}. The first step is to find the first 'n' so that 4ⁿ*x is bigger than y, so in this case, n=6.
Then {x2ⁿ,y2^-n}={x,y}, but now the numbers are much "closer." We get {13376,10223.765625}. Now factor out the largest power of 2 that fits in both numbers, which is 8192. We have 8192{1.6328125,1.248018265}.

All of that preemptive stuff has set up a new {x,y} where we have y<x<2y, so 1<x/y<2, or 1>y/x>.5. This will allow for some fine-tuned magic in a bit. What we should realize is that if we can get x=y somehow, then we have {x,y}={r,r}=##\sqrt(r^{2})=r##. So notice that if we choose some number ##a##, we have ##{xa,y/a}=\sqrt{xay/a}=\sqrt{xy}={x,y}##. So what we need to do is find some ##a## such that xa=y/a, so xa²=y, ##a=\sqrt{y/x}##. However, we want to avoid computing square roots in our computation of square roots, so what we can do is make an estimate for a. As luck would have it, based on our restrictions above, we can make a pretty decent estimate of the square root just by averaging the quotient with 1. In this case, that is ##\frac{\frac{1.248018265}{1.6328125}+1}{2}\approx .8821682725##. So then we get 8192{1.440415382,1.414716788}. Iterating:
8192{1.427566085,1.427450431}
8192{1.427508258,1.427508256}
8192{1.427508257,1.427508257}

And for safety, we could do one more iteration and average the two numbers in case of extra hidden precision, but we get the final result 11694.147638883, which is pretty close to 11694.147638883.

But that required 5 iterations requiring 2 divisions and 1 multiplication each. This is not the promised 2 iterations! So here is where we do some magic. At the division step for y/x, we want to find the square root, right? Why not try a couple of iterations at this step? We would get {1,y/x} as the square root, our estimate is still ##\frac{\frac{y}{x}+1}{2}##, so if we apply this, we get ##{\frac{\frac{y}{x}+1}{2},y/x\frac{2}{\frac{y}{x}+1}}\approx \frac{\frac{y}{x}+1}{4}+\frac{y/x}{\frac{y}{x}+1}##. So in pseudo-code, we could do:
1+y/x→a
1-1/a+a>>2→a
{xa,y/a}→{x,y}
This is the same as performing two iterations at the cost of 3 divisions, 1 multiplication, compared to 4 divisions, 2 multiplications. So we save 1 division, 1 multiplication.



So, for some working TI-BASIC code:
{max(Ans),min(Ans→L1
int(log(L1(1)/L1(2))/log(4
L1{2^-Ans,2^Ans→L1
2^int(log(Ans(1))/log(2→P
L1/Ans→L1
While E-10<abs(max(deltaList(L1    ;engineering E
1+L1(2)/L1(1
1-Ans^-1+.25Ans
L1{Ans,Ans^-1→L1
End
Pmean(L1


Basically, a geometric mean routine would be used in the super-cool "Arithmetic-Geometric Mean" class of algorithms. What I was doing was making an efficient alternative to multiplying, then taking the square root. It is not more efficient in TI-BASIC, but it can be much more efficient in other applications. Plus, there is no overflow risk.

So as an example, you could do this on the homescreen:
{976,587→L1
1+L1(2)/L1(1:1-Ans^-1+.25Ans:L1{Ans,Ans^-1→L1
Just hit enter one or two more times and voila, you have the square root of 976*587.

Actually, yes it could, as long as you had two (not necessarily integer) factors. That is what makes it cool. If you wanted to find the geometric mean of 9E99 and 9.99E99, the multiplication would overflow the calculator, but using the BASIC code given, it return 9.482088378E99 properly.

Necro-update.

I was going to go to sleep a couple hours ago, but as soon as I turned out the lights, for whatever reason, I thought of a way to tweak this algorithm to be even faster on the Z80. I've spent the last couple of hours refining my algorithm. I'll refrain from a full mathematical proof, because 1] I'm tired, and 2] I haven't made a formal proof yet

The principal is the same: We want some algorithm for making x_n and y_n converge to some "r". We also need to be able to relate x_n and x_n to the original x and y.

One further constraint is to use very, very simple operations (shift-and-add).

If x_n<y_n, and I define
x_n+1=x_n(1+2^-k_n)²
y_n+1=y_n(1-2^-k_n)²
then we have that x_n+1y_n+1=x_ny_n(1-4^-k_n)²

So if we keep track of the factors (1-4^-k_n)², we can just divide them out of the final result. So basically, start with a constant factor of 1. Each iteration, multiply that by your sqrt((1-4^-k_n)²)=(1-4^-k_n). After all of the iterations are done, we have x_n=y_n and xy(c²)=x_ny_n=x_n². So then:
xy(c²)=x_n²
xy=x_n²/c²
sqrt(xy)=x_n/c

The algorithm below appears to require n iterations to compute 4n bits of accuracy.
;think of all values here as 2.m fixed point numbers.
INPUTS: X,Y are on [1,4)
    c=1
    iterate n times
        a=y-x
        if underflow
            a=-a
            (x,y) = (y,x)       //swap x and y, so x<y
        n=-2+intlg(a)           //intlg(a) basically finds the negative of the index of the first '1' bit in the binary expansion of a. Basically, the "clz" instruction for some processors.
        c-=c>>2n
        x+=x>>n
        x+=x>>n
        y-=y>>n
        y-=y>>n
    return x/c


It almost certainly won't run faster (I haven't tried in a while). It's only useful in situations where excessively large numbers are used (the BASIC version).

However, if you are working in a 32-bit environment, there is no risk of overflowing into a 64-bit int via multiplication, so all operations can take place with 32 bits.

The algorithm is best suited for a machine that can perform arbitrary shifts (the Z80 can only shift one place to the left or right) and even better if it has a "clz" instruction. Then it averages roughly 5 shift-and-adds, 1 clz, and 1 subtraction for each 4 bits of accuracy.

The algorithm is best suited for fixed point math in the form of 2.n (so 2 integer bits, n fractional bits). This can be easily achieved for floats where the exponent can be set aside for later, and the mantissa is almost always already in a 1.n form. For integers, you just need to shift left an even number of times until the upper 2 bits are non-zero (so 01, 10,11).