Updated 2016-01-28 20:18:17 by mpr ▲

mpr 2016-01-28: Karatsuba multiplication is an 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 this paper 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 { 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 " + "]
 }

# 8x^3 + 7x^2 + 8x + 32
set A {32 8 7 8}

# 4x^3 + 8x^2 + 12
set B {12 0 8 4}

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


# outputs:
# 32x^6 + 92x^5 + 88x^4 + 288x^3 + 340x^2 + 96x^1 + 384