unusual recursive factoring algorithm in C

[USENET news]

An extraordinary recursive factoring algorithm has been discovered (?) by
Risto Lankinen (ri...@yj.data.nokia.fi) based on mathematical properties he
described in sci.math article <6...@tuura.UUCP>. The `Lankinen algorithm',
which finds all natural numbers u, v such that u * v = n for an odd n, is
as follows:

1. select an odd n > 1.
2. u <- 1, v <- 1, b <- 1, n <- (n - 1).
3. n <- (n / 2).
4. if n = 0, (then) u * v = n; return.
5. if n < 0, (then) u * v <> n; return.
6. b <- (b * 2).
7. if n is odd, (then)
7a. restart at 3 with n <- (n - u), v <- (v + b),
7b. and again with n <- (n - v), u <- (u + b).
8. (otherwise) for even n,
8a. restart at 3,
8b. and again with n <- (n - u - v - b), u <- (u + b), v <- (v + b).

Execution time is comparable to the technique of division by all odd numbers
less than the square root of n.

(For an explanation see the cited article or
the information at the end of this article.)

- - -

Following is an implementation of Lankinen's algorithm, with embellishments,
in the C language:

DIVISOR.C

1: #include <stdio.h>
2:
3: void factor(long remain,
4: long left,
5: long right,
6: long bit)
7:
8: {
9: if (remain >>= 1)
10: {
11: if (remain > 0)
12: {
13: bit <<= 1;
14: if (remain & 1)
15: {
16: factor(remain - left, left, right | bit, bit);
17: if (left != right)
18: factor(remain - right, right, left | bit, bit);
19:
20: }
21: else
22: {
23: factor(remain, left, right, bit);
24: factor(remain - left - right - bit, left | bit, right | bit, bit);
25: }
26: }
27: }
28: else
29: printf("%li %li\n", left, right);
30:
31: }
32:
33: void main(void)
34:
35: {
36: long number;
37:
38: scanf("%li", &number);
39: while (!(number & 1))
40: number >>= 1;
41: factor(number, 1, 1, 1);
42: }

The correspondence is as follows:

step -> lines | variables
---- ----- | ---------
1 38-40 | n -> number, remain
2 41,9 | u left
3 9 | v right
4 9,28-29 | b bit
5 11,21-25 |
6 13 |
7 14,15-20 |
7a 16 |
7b 18* |
8 14,21 |
8a 23 |
8b 24 |

* Embellishment, see notes below.

The C language is ideal for implementing Lankinen's algorithm because both are
optimized for binary machine languages. Hence many arithmetical operations in
the algorithm can be replaced with C's bitwise operators for maximized
computational efficiency:

This function at [line] replaces at [step] (respectively).
---- -------- ------ -------- ------
`|' (binary OR) 16,18,24 addition 7a,7b,8b
`>>' (shift right) 9, 40 division 3 (1 implicit)
'<<' (shift left) 13 multiplication 6

Moreover, the parity test at steps 7 and 8 can be reduced to simply branching
on the state of the rightmost bit of the variable at lines 14 and 21. Recall
the existence of the implicit nonzero test in C, utilized in lines 14 and 39,
which can even be combined with assignment as at 9 (here some compilers may
warn of a possibly unintentional assignment instead of comparison).

This C routine embellishes steps 7 and 8 in Lankinen's algorithm at lines 16-18
and 23-24:

* The conditional statement at line 17 decreases total processing time roughly
by a factor of two by avoiding a recursive call that results in finding
reversed (duplicate) pairs. Hence at line 29 the pair (left,right) is
unique; for no (u,v) will the program also report (v,u). (If the trial
factors are equivalent to some point, the recursive branches originating
there would be symmetric with respect to each other.)

* Lankinen left unspecified the order of the two recursive calls in steps 7
and 8 (they were assigned definite order above for exposition), but when 7 is
encoded as above, in concert with the test at 17 and the implied interchange
of trial factors at 18 (compare with step 7b), an ordering is established
such that for all divisor pairs (left,right) found at line 29, left <= right.

- - -

Lankinen's basic algorithm can be modified slightly to determine the complete
prime decomposition of a number:

1: #include <stdio.h>
2:
3: #define factor(what) find(what, 1, 1, 1, what)
4:
5: long find(long remain,
6: long left,
7: long right,
8: long bit,
9: long now)
10:
11: {
12: if (right > now)
13: return now;
14:
15: if (remain >>= 1)
16: {
17: if (remain > 0)
18: {
19: bit <<= 1;
20: if (remain & 1)
21: {
22: if (left != right)
23: now = find(remain - right, right, left | bit, bit, now);
24: return find(remain - left, left, right | bit, bit, now);
25:
26: }
27: else
28: return find(remain,
29: left,
30: right,
31: bit,
32: find(remain - left - right - bit,
33: left | bit,
34: right | bit,
35: bit,
36: now));
37: }
38: else
39: return now;
40: }
41: else
42: if (now > 1)
43: if (left == 1)
44: printf("%li ", right);
45: else
46: {
47: factor(left);
48: factor(right);
49: now /= right;
50: }
51:
52: return now;
53: }
54:
55: void main(void)
56:
57: {
58: long number;
59:
60: scanf("%li", &number);
61: while (!(number & 1))
62: {
63: printf("2 ");
64: number >>= 1;
65: }
66: if (number > 1)
67: factor(number);
68: }

The function `factor' is replaced by a macro that calls on the new function
`find' that decomposes an odd number into its prime factors. This program
utilizes `intercommunication' between the branches of the recursive tree.
Whenever `factor' is started, the macro (line 3) places the original n in
the `now' variable (line 9). This variable represents the current unfactored
part of n (the original number divided by all actual factors found at some
point). It is passed from `parent' to `child' process as a parameter and as
the function result in the opposite direction.

Then if the tentative factor under consideration is ever greater than `now',
it is impossibly large. The comparison for this test occurs at line 12.
At that point both `left' and `right' are trial factors but again right >= left
by ordering so only `right' need be tested. The division occurs at line 49.

With adjustment to the unspecified ordering in the algorithm at lines 22-24
and 28-36, there is an implicit guarantee that the divisor that is the number
itself (left = 1, right = n) will be found last. Hence, if an unfactored part
remains (line 42), either it is composite and all factors will be found
recursively at lines 47-49, or it is prime and reported at line 44.

The `find' routine can be rewritten in several other forms. In this version
the `else' keywords at lines 38 and 41 are redundant and should be removed
if they improve the efficiency of the compilation (they were included here
for clarity). Also, the routine is readily modified to use a single exit point
at line 52 by eliminating the other three at lines 13, 24, and 28 with block
scoping and assignment to `now' (the version above attempts to maximize speed).

- - -

For a more rigorous mathematical treatment of Lankinen's approach to factoring,
see his preliminary findings in <6...@tuura.UUCP>. The description that follows
is intended to be more intuitively appealing than logically precise.

Lankinen's algorithm can be thought of as the `inverse' of binary
multiplication (otherwise known as the `Russian Peasant's Method')
in the special case where both the multiplier and multiplicand are odd
(and hence the product also).

The algorithm simply explores all possibilities in u and v that could match n
by focusing on the 2nd bit from the right and advancing left one bit at a time.
(The rightmost bit in n is matched immediately because both factors are odd.)
The fundamental rule is employed that if the next bit in n is 1, then it is
matched only when the next bits in u and v are different from each other,
and if the digit is 0 the bits must match. In either case, there are only two
possibilities, each of which is pursued recursively.

As the trial divisors accumulate, so does a carry that has to be accounted for
in the subsequent comparisons. This is accomplished exactly and conveniently
by subtracting the newly accumulated amount from n prior to entering each
recursive branch. If n < 0, the product of the trial factors has accrued a
carry that exceeds the original n to be factored, thereby eliminating them.
Otherwise, if n = 0, they are actual factors.

ld23...@longs.LANCE.ColoState.EDU