Version 1 of karatsuba polynomial multiplication

Updated 2016-01-28 20:03:32 by mpr

mpr 2016-01-28: Karatsuba multiplication is a algorithm invented by Anatoly Karatsuba in 1960 to multiply n-digit numbers in less than O(n^2) time. The algorithm was invented for use on integers, but can be applied to polynomials represented as coefficient arrays. Here is an implementation based on the paper this http://www.ijcse.com/docs/INDJCSE12-03-01-082.pdf .

proc odd { n } { expr { $n % 2 } }

    proc zeros { n } { lrepeat $n 0 }

    proc swap { x y } {
        uplevel [subst { set tmp $$x; set $x $$y; set $y \$tmp }]
    }

    # Implements Karatsuba fast multiplication on polynomials represented
    # by coefficient lists.
    proc karatsuba_poly { A B } {
        # make A the longest list (highest degree polynomial).
        if { [llength $B] > [llength $A] } { swap A B }

        # put max length in n,
        # min length in m.
        set n [llength $A]
        set m [llength $B]

        set D [zeros [expr { ($n+$m)/2 }]]
        set S [zeros [expr { $n+$m-1 }]]
        set T [zeros [expr { $n+$m-1 }]]
        set C [zeros [expr { $n+$m-1 }]]

        # Compute products of same-indexed coefficients.
        for {set i 0} {$i < $m} {incr i} {
            lset D $i [expr { [lindex $A $i] * [lindex $B $i] }]
        }

        for {set i 1} {$i < $n+$m-2} {incr i} {
            for {set p 0} {$p < min($m,$i)} {incr p} {
                set q [expr { $i - $p }]
                if { $p >= $q } { break }

                set ap [lindex $A $p]
                set aq [lindex $A $q]
                set bp [lindex $B $p]
                set dp [lindex $D $p]

                if { $q < $m && $p < $m } {
                    set bq [lindex $B $q]
                    set dq [lindex $D $q]

                    lset S $i [expr { [lindex $S $i] + (($ap + $aq) * ($bp + $bq)) }]
                    lset T $i [expr { [lindex $T $i] + $dp + $dq }]
                } elseif { $q >= $m && $q < $n } {
                    lset S $i [expr { [lindex $S $i] + (($ap + $aq) * $bp) }]
                    lset T $i [expr { [lindex $T $i] + $dp }]
                }
            }
        }

        # Set coeffients of resulting polynomial.
        for {set i 0} {$i < $n+$m-1} {incr i} {
            if { $i == 0 } {
                lset C $i [lindex $D $i]
            } elseif { $i == $n+$m-2 } {
                lset C $i [expr { [lindex $A end] * [lindex $B end] }]
            } elseif { [odd $i] } {
                lset C $i [expr { [lindex $S $i] - [lindex $T $i] }]
            } else {
                lset C $i [expr { [lindex $S $i] - [lindex $T $i] + [lindex $D [expr { $i/2 }]] }]
            }
        }

        set C
    }
 
    proc polly { coeffs var } {
        for {set i 1} {$i < [llength $coeffs]} {incr i} {
            lset coeffs $i "[lindex $coeffs $i]$var^$i"
        }
        lreverse [join $coeffs " + "]
    }

    set A {32 8 7 8}
    set B {12 0 8 4}

    puts [polly [karatsuba_poly $A $B] x]