univariate Karatsuba in Maxima



Hi!

I am still trying to implement Karatsuba
for univariate polynomials in Maxima
(I am not using Common Lisp because I am not 
familiar with Lisp)
Sofar what I did is way too much slower
than ratexpand for any size of polynomials
I have tried.

I wonder whether ratexpand uses something
better than "school" multiplication.

I had previously implemented Karatsuba in C++
and it worked quite well:

template<class coeRing>
univariateDensePolynomial<coeRing>
  _Karatsuba(const univariateDensePolynomial<coeRing> &lhs,
             const univariateDensePolynomial<coeRing> &rhs,
             const long degree)
{
if (lhs.degree() <= THRESHOLD ||  rhs.degree() <= THRESHOLD)
  {
   return lhs * rhs;
  };
univariateDensePolynomial<coeRing> lhsL;
...
univariateDensePolynomial<coeRing> prodH;
univariateDensePolynomial<coeRing> res;

lhsL = lhs.truncatePower(degree/2);
lhsH = lhs.extractHighPower(degree/2);
rhsL = rhs.truncatePower(degree/2);
rhsH = rhs.extractHighPower(degree/2);
prodL = _Karatsuba(lhsL,rhsL,degree/2);
prodH = _Karatsuba(lhsH,rhsH,degree/2);
prodM = _Karatsuba((lhsL+lhsH),(rhsL+rhsH),degree/2) - prodL - prodH;
res = prodL+prodM.shiftPower(degree/2)+prodH.shiftPower(degree);
return res;
};

Note: the C++ version beats "school multiplication"
with a trade-off THRESHOLD starting from about 100.

My first attempt in Maxima is:

karatsuba(lhs,rhs,var,deg) :=
 block([lhsL,lhsH,rhsL,rhsH,prodL,prodM,prodH,res],
   lhsDeg : degree(lhs,var),
   rhsDeg : degree(rhs,var),
   if lhsDeg <= KARATSUBA_THRESHOLD or rhsDeg <= KARATSUBA_THRESHOLD then
     return(ratexpand(lhs*rhs)),    
   lhsL : truncatePol(lhs,var,deg/2-1,0),
   lhsH : truncatePol(lhs,var,lhsDeg,deg/2)/x^(deg/2),
   rhsL : truncatePol(rhs,var,deg/2-1,0),
   rhsH : truncatePol(rhs,var,rhsDeg,deg/2)/x^(deg/2),

   prodL : karatsuba(lhsL,rhsL,var,deg/2),
   prodH : karatsuba(lhsH,rhsH,var,deg/2),
   return(ratexpand(prodL+(karatsuba((lhsL+lhsH),
     (rhsL+rhsH),var,deg/2) - prodL - prodH)*var^(deg/2)+prodH*var^deg))
);

The algorithm is quite simple and I don't understand what I should
do to make it faster than ratexpand.

Should I use arrays to implement polynomials?

Is someone planning to have it in the next Maxima version?

  Fabrizio