General Discussion > Math and Science

My way to compute Geometric Means (and square roots)

(1/2) > >>

Xeda112358:
So I've been developing a way to compute the geometric mean of two numbers (defined by ##\sqrt(xy)##MathJax.Hub.Queue(["Typeset", MathJax.Hub, document.getElementById("bbclatex6630bcce3665f")]);console.log("Queued!");) 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 4n*x is bigger than y, so in this case, n=6.
Then {x2n,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##MathJax.Hub.Queue(["Typeset", MathJax.Hub, document.getElementById("bbclatex6630bcce3668b")]);console.log("Queued!");. So notice that if we choose some number ##a##MathJax.Hub.Queue(["Typeset", MathJax.Hub, document.getElementById("bbclatex6630bcce36696")]);console.log("Queued!");, we have ##{xa,y/a}=\sqrt{xay/a}=\sqrt{xy}={x,y}##MathJax.Hub.Queue(["Typeset", MathJax.Hub, document.getElementById("bbclatex6630bcce366a0")]);console.log("Queued!");. So what we need to do is find some ##a##MathJax.Hub.Queue(["Typeset", MathJax.Hub, document.getElementById("bbclatex6630bcce366aa")]);console.log("Queued!"); such that xa=y/a, so xa2=y, ##a=\sqrt{y/x}##MathJax.Hub.Queue(["Typeset", MathJax.Hub, document.getElementById("bbclatex6630bcce366bc")]);console.log("Queued!");. 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##MathJax.Hub.Queue(["Typeset", MathJax.Hub, document.getElementById("bbclatex6630bcce366c7")]);console.log("Queued!");. 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}##MathJax.Hub.Queue(["Typeset", MathJax.Hub, document.getElementById("bbclatex6630bcce366d1")]);console.log("Queued!");, 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}##MathJax.Hub.Queue(["Typeset", MathJax.Hub, document.getElementById("bbclatex6630bcce366da")]);console.log("Queued!");. So in pseudo-code, we could do:

--- Code: ---1+y/x→a
1-1/a+a>>2→a
{xa,y/a}→{x,y}

--- End code ---
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:

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

--- End code ---

aeTIos:
You're just a math goddess. I don't understand 50% of what you're saying here but whatever, might come in handy one time :D

Xeda112358:
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:

--- Code: ---{976,587→L1
1+L1(2)/L1(1:1-Ans^-1+.25Ans:L1{Ans,Ans^-1→L1

--- End code ---
Just hit enter one or two more times and voila, you have the square root of 976*587.

aeTIos:
So could this be used to calc the square root of numbers that would normally overflow the calc? Ie 1e200 etc

Xeda112358:
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.

Navigation

[0] Message Index

[#] Next page