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