Exponentiation by squaring (AWK/dc)

From LiteratePrograms
Jump to: navigation, search

[edit] the theory

The straightforward way to exponentiate is via repeated multiplication, using the recurrence

\begin{matrix} a^0 & = & 1 \\ a^n & = & a a^{n-1} \end{matrix}

But if we view the exponent as a binary number (which we can deconstruct bit by bit) instead of as a natural number (which we must deconstruct by decrementing), it's possible to build up the result with a much shorter sequence of multiplications. The tradeoff is that instead of always multiplying by the base, we will be squaring repeatedly, and multiplying by the base only in the cases where we find "1" bits.

a^0 & = & 1 \\ 
a^{s0} & = & { ( a^s ) }^2 \\
a^{s1} & = & a { ( a^s ) }^2 

If one is folding constants during code generation this is an easy tradeoff to make: a few more branches, taken once during compilation, can yield much shorter straight-line code paths, which will be executed many times.

To simulate this process, we will generate some code to evaluate an exponential (minimizing the multiplications) and then interpret it in dc.

<<evaluate expression>>=
function exponentiate(a,b)      { print genexpr(a,b) | "dc" }

Not only does the two-process split simulate the difference between compile-time calculation and run-time calculation, but it offers the additional advantage that dc performs arbitrary length integer arithmetic, while awk on its own would silently shift over to floating point. For instance, we can now evaluate:


[edit] the practice

We want to produce output, so our expression will at least contain the p command to print the result. (this is done here rather than elsewhere to finesse AWK's lack of local variables)

The loop is rather straightforward. We shift out bits from the bottom of the exponent and build up our expression. It would be even simpler to shift while(b > 0), but the rationale for special casing the final bit is given below.

Exercise 1: recode with a loop which shifts the exponent all the way to 0, treating all bits identically Exercise 2: add code to Ex. 1 to patch the final expression, removing excess multiplies

<<generate expression>>=
function genexpr(a,b) {
    e = "p"
    for(;;) {
        if(b < 2) {
           most significant bit
            return e
        lower order bits
        b = int(b/2)

For most bits, we square (d*) the running result from the higher order bits. If the bit in question is a 1 bit, we also need a multiply by the base.

NB. we shift out from the bottom, so we process the exponent right-to-left. But the final calculation is carried out left-to-right. This is not a problem, as there is nothing constraining us to build the expression in the interpetation order, so we simply build it in the reverse.

Question 3: can reverse-order generation be used to simplify tail-call elimination? "

<<lower order bits>>=
        if(b % 2)      e = "d*" a "*" e
        else           e = "d*"       e

Finally, we come to the most significant bit. We could treat it like all the other bits, prepending a "1 " to the expression to get things started (dc is not sophisticated enough to figure out the unit for multiplication, and will complain bitterly about empty stacks), but taking advantage of a little constant propagation

 1 d*a* = a
 1 d*   = 1

we can trim a multiply or two from the final expression.

<<most significant bit>>=
            if(b % 2)  e = a " " e
            else       e = 1 " " e

Note that we treat the value of a as opaque throughout.

Exercise 4: rewrite to generate the sequence of assembly instructions for an arbitrary base

[edit] wrapping up

Now we just need a little glue to handle rudimentary input

<<handle input>>=
BEGIN                           { FS="^" }
$1 ~ /[0-9]+/ && $2 ~ /[0-9]+/  { exponentiate($1,$2); next }
BEGIN                           { print "enter expressions of the form i^j" }

and we can bundle up the script

#!/usr/bin/awk -f
handle input
generate expression
evaluate expression
Download code