Author Topic: My way to compute Geometric Means (and square roots)  (Read 7536 times)

0 Members and 1 Guest are viewing this topic.

Offline Xeda112358

  • they/them
  • Moderator
  • LV12 Extreme Poster (Next: 5000)
  • ************
  • Posts: 4704
  • Rating: +719/-6
  • Calc-u-lator, do doo doo do do do.
    • View Profile
My way to compute Geometric Means (and square roots)
« on: March 31, 2014, 04:43:34 pm »
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 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##. 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 xa2=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:
Code: [Select]
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:
Code: [Select]
{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

Offline aeTIos

  • Nonbinary computing specialist
  • LV12 Extreme Poster (Next: 5000)
  • ************
  • Posts: 3915
  • Rating: +184/-32
    • View Profile
    • wank.party
Re: My way to compute Geometric Means (and square roots)
« Reply #1 on: March 31, 2014, 05:06:15 pm »
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
I'm not a nerd but I pretend:

Offline Xeda112358

  • they/them
  • Moderator
  • LV12 Extreme Poster (Next: 5000)
  • ************
  • Posts: 4704
  • Rating: +719/-6
  • Calc-u-lator, do doo doo do do do.
    • View Profile
Re: My way to compute Geometric Means (and square roots)
« Reply #2 on: March 31, 2014, 05:25:26 pm »
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: [Select]
{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.

Offline aeTIos

  • Nonbinary computing specialist
  • LV12 Extreme Poster (Next: 5000)
  • ************
  • Posts: 3915
  • Rating: +184/-32
    • View Profile
    • wank.party
Re: My way to compute Geometric Means (and square roots)
« Reply #3 on: March 31, 2014, 05:30:26 pm »
So could this be used to calc the square root of numbers that would normally overflow the calc? Ie 1e200 etc
I'm not a nerd but I pretend:

Offline Xeda112358

  • they/them
  • Moderator
  • LV12 Extreme Poster (Next: 5000)
  • ************
  • Posts: 4704
  • Rating: +719/-6
  • Calc-u-lator, do doo doo do do do.
    • View Profile
Re: My way to compute Geometric Means (and square roots)
« Reply #4 on: March 31, 2014, 10:50:33 pm »
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.

Offline aeTIos

  • Nonbinary computing specialist
  • LV12 Extreme Poster (Next: 5000)
  • ************
  • Posts: 3915
  • Rating: +184/-32
    • View Profile
    • wank.party
Re: My way to compute Geometric Means (and square roots)
« Reply #5 on: April 01, 2014, 02:39:42 am »
Aww I was planning to calculate the root of a super huge prime number :P
I'm not a nerd but I pretend:

Offline Xeda112358

  • they/them
  • Moderator
  • LV12 Extreme Poster (Next: 5000)
  • ************
  • Posts: 4704
  • Rating: +719/-6
  • Calc-u-lator, do doo doo do do do.
    • View Profile
Re: My way to compute Geometric Means (and square roots)
« Reply #6 on: June 30, 2015, 03:13:52 am »
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 :P

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

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

If xn<yn, and I define
xn+1=xn(1+2-kn)2
yn+1=yn(1-2-kn)2
then we have that xn+1yn+1=xnyn(1-4-kn)2

So if we keep track of the factors (1-4-kn)2, 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-kn)2)=(1-4-kn). After all of the iterations are done, we have xn=yn and xy(c2)=xnyn=xn2. So then:
xy(c2)=xn2
xy=xn2/c2
sqrt(xy)=xn/c

The algorithm below appears to require n iterations to compute 4n bits of accuracy.
Code: [Select]
;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

Offline c4ooo

  • LV5 Advanced (Next: 300)
  • *****
  • Posts: 252
  • Rating: +10/-1
  • The impossible chemical compound.
    • View Profile
Re: My way to compute Geometric Means (and square roots)
« Reply #7 on: June 30, 2015, 12:43:08 pm »
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
^That  *.*
And btw just interested, is that ti-basic code faster then just √(x*y)?
-German Kuznetsov
The impossible chemical compound.

Offline Xeda112358

  • they/them
  • Moderator
  • LV12 Extreme Poster (Next: 5000)
  • ************
  • Posts: 4704
  • Rating: +719/-6
  • Calc-u-lator, do doo doo do do do.
    • View Profile
Re: My way to compute Geometric Means (and square roots)
« Reply #8 on: July 01, 2015, 11:34:31 am »
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).