sumpolyk(P,a_,b_) = {
    local(p,m,Q,y);
    p = poldegree(P,k);
    y = [];
    for(m=1,p+2,y=concat(y,sum(k=1,m,eval(P))));
    Q = polinterpolate(y);
    Q = Q - subst(Q,x,a_-1);
    Q = subst(Q,x,b_);
    return(factorisation(Q));
}