Algorithm for division by 2

[eudoxie]

I don't know if similar algorithms have been published before, but I've developed a linear time complexity division-by-two algorithm for balanced ternary numbers.

It's based on the recursive relation that

1/2 = 1/3 + (1/3)*(1/2)

(and division by 3 is achievable by shifting)

So, to divide some number ABCD, you do the following

Code: Select all

A B C D / 2
-----------
  A B C . D 
    A B . C D
 +    A . B C D
-----------------
[0, A, A+B, A+B+C].[B + C + D, C+D, D]

Where the right half is the residual (add the carry of that to the result to fix rounding errors)

Code: Select all

void bw_bisect(word target) {
	register int i;
	register int tmp = 0, sum = 0;
	register int residual = 0;
	register int power = 1;

	for(i = 0; i < WORD_DIGITS-1; i++, power*=3) {
		tmp += target[i];
		residual += power * tmp;
	}

	for(i = WORD_DIGITS-1; i >= 0; i--) {
		tmp = target[i];

		target[i] = sum;
		sum += tmp;
	}

	if(residual > 0) {
		target[0] += (residual + (NBWORD_MAX/3)/2) / (NBWORD_MAX/3);
	} else {
		target[0] += (residual - (NBWORD_MAX/3)/2) / (NBWORD_MAX/3);
	}


	bw_ripple(target);
}


/* Helper function. "smears" out a number so that
 * individual digits hold numbers that are in the
 * range -1,...,1 */
int bw_ripple(word target) {
	register int i;
	register int carry = 0;

	for(i = 0; i < WORD_DIGITS; i++) {
		target[i] += carry;

		if(target[i] > 0) 
			carry = (target[i] + 1) / 3;
		else 
			carry = (target[i] - 1) / 3;

		target[i] -= 3 * carry;
	}

	return carry;
}

[Shaos]

Great! In one Russian topic here we discussed ternary algorithms of multiplication and division by 10. See automatic translation of that topic:

http://translate.google.com/translate?j ... ry_state0=

[eudoxie]

Interesting...

I think it is probably possible to re-write my algorithm in a more efficient form. Looking over my code, it feels a bit awkward.

[Mac Buster]

Incredible! It is very interesting thing indeed

[eudoxie]

It seems possible to express any integer division as repeated shifting, addition and multiplication.

The method for finding how to do it is taking

1/n = 1/3 + (3 - n)/3 * (1/n)

So, you get
1/2 = 1/3 + 1/9 + 1/27 + 1/81 ...
1/4 = 1/3 - 1/9 + 1/27 - 1/81 ...
1/5 = 2/3 - 2/9 + 2/27 - 2/81 ...
1/7 = 4/3 - 4/9 + 4/27 - 4/81 ...
etc.

[Shaos]

What about x/10?

[eudoxie]

1/10 = 7/3 - 7/9 + 7/27 - 7/81 ...

You could express multiplication by 7 as multiplication by 9 and subtracting twice I guess.

But I think it's smarter to express division by 10 as repeated division by 9:

1/10 = 1/9 - 1/9 * 1/10 =>
1/10 = 1/9 - 1/81 + 1/729 - 1/6561 ...

You get much fewer terms, and don't need to worry about multiplication by 7.

If consideration is taken to the carry of the sum of the decimal part, this will produce exact results.

--------------

In it's most general form, division by some number n in terms of division by some other number b can be expressed like this:



The relationship is derived by evaluating 1/n - 1/b and finding that
1/n - 1/b = (b-n)/nb = (b-n)(1/b)(1/n)
Then shifting over 1/b to the right side
1/n = 1/b + (b-n)(1/b)(1/n)
Then you simply replace the (1/n) in the right hand part of the equation with the right hand equation (over and over again)
1/n = 1/b + (b-n)(1/b)(1/b + (b-n)(1/b)(1/b + (b-n)(1/b)(1/b + (b-n)(1/b)([...]))))
If you clean that mess up, you get the pretty relationship above.

[Shaos]

Finally I got first part. Because of
where |q|<1
and if we have q=1/3 then 1/(1-q)=1/(1-1/3)=1/(2/3)=3/2
It means we should dived by 3 both sides or
1/3*sum(1/3^n)=1/2 (sum from n=0)
or
sum(1/3^n)=1/2 (sum from n=1)

Code: Select all

#!/usr/local/bin/hopeless -f

dec round : num -> num;
--- round(x) <= if x-f < 0.5 then f
                else f+1 where f==floor(x);

dec tri10b : num -> num;
--- tri10b(x) <= round((x-round((x-round(x/9))/9))/9);

dec tri2b : num -> num;
--- tri2b(xx) <= let x==xx+1 in round((x+round((x+round((x+round((x+round((x+round(x/3))/3))/3))/3))/3))/3);

dec try : num # list(num # num) -> list(num # num);
--- try(x,l) <= let fx==tri10b(x) in 
                if round(x/10) = fx then try(x+1,ll) 
                else ll where ll==(x,fx)::l;

dec try2 : num # list(num # num) -> list(num # num);
--- try2(x,l) <= let fx==tri2b(x) in 
                if round(x/2) = fx then try2(x+1,ll) 
                else ll where ll==(x,fx)::l;
try2(0,nil);

algorithm tri2b is correct for values up to 729 (tri10b - up to 734)

P.S. correction +1 for argument x required in tri2b for "round-like" results (1/2=1, 3/2=2 etc), without it result will be "floor-like" (1/2=0, 3/2=1 etc)

[eudoxie]

That's pretty neat

I've cleaned up the derivation of a more general case, I hope it's easier to follow.

[Shaos]

It's brilliant
Thanks!

[Shaos]

again, no images

P.S. fixed