## Version 7 of karatsuba polynomial multiplication

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