Efficient implementation of the divided difference

Question

Question. I am looking for more efficients ways to compute Newton's divided difference, i.e. the operator $$\partial : f(x,y) \mapsto \frac{f(x,y)-f(y,x)}{x-y} \ , $$

for the case where $f$ is a polynomial in $x,y$.

Of course the above can be implemented as

 Par = (# - ReplaceAll[#,{x->y,y->x}]) /(x - y) & ;

but for $f$ a polynomial I suspect that there are faster ways to get the same result. Indeed, since the numerator is antisymmetrized in $x,y$ it is divisible by $x-y$, so if $f$ is a polynomial then so is $\partial f$.

It is really this quotient-free form that I am after, and one might expect the division to be relatively time-costly.

NB. The reason I would like to speed this up is that I really want to work with the case of more variables and various $\partial$s, where things start to become slower quite rapidly. See below for an example. For my question, though, the case with two variables suffices.

---

To get things started let me give one alternative implementation (although it actually appears to be slower than the above Par in examples, see below).

One can check that on monomials we have $$\partial(x^m y^n) = \text{sgn}(m-n) \, (x \, y)^{\min(m,n)} \, \sum_{i=1}^{|m-n|} x^{|m-n|-i} \, y^{i-1} \ , $$ where $\text{sgn}(0)=0$:

 Table[
  Par[x^m y^n] == Sign[m-n] (x y)^Min[m,n] Sum[x^(Abs[m-n] - i) y^(i-1),{i,Abs[m-n]}],
  {m, 0, 3}, {n, 0, 3}
 ] //Simplify //Flatten //Apply[And]
 (* Out: True *)

This leads me to another definition, as follows. First implement additivity:

par[expr_Plus] := par /@ expr

Next we reduce to monomials in $x,y$:

par[term_] /; FreeQ[term, x] && FreeQ[term,y] := 0 (* discard terms not involving x or y *)
par[a___ b_ c___] /; FreeQ[b, x] && FreeQ[b, y] := b par[a c] (* extract factors independent of x,y *)
par[term_Times] := par @ Expand @ term (* expand remaing products involving x,y *)

To extract factors of $x\,y$ we can use the comment of AccidentalFourierTransform to improve my previous code:

par[x^m_. y^n_.] := Sign[m-n] (x y)^Min[m,n] par[x^Abs[m-n]]

where I use $\partial(y^n) = -\partial(x^n)$. Finally,

par[y^n_.] := -par[x^n]
par[x^m_.] := Sum[x^(m - i) y^(i - 1), {i,m}]

completes the definition.

Bonus question. Is there a nicer way to implement the above, which in particular gets rid of the need to distinguish between $x$ and higher powers of $x$, and likewise for $y$?

---

The example that I actually have in mind relates to the Hecke algebra which acts on the space of polynomials in $N$ variables $z_i$ as follows. Writing $\partial_i$ for the divided difference with respect to $z_i$ and $z_{i+1}$, the generators act as $T_i = -(q \,z_i-q^{-1}\,z_{i+1})\,\partial_i + q\ $ for $1\leq i\leq N-1$: these obey the braid relations along with the 'Hecke condition' $(T_i - q)(T_i+q^{-1})=0$. In any case, we define

Par[i_] := (# - ReplaceAll[#,{z[i] -> z[i+1],z[i+1]->z[i]}]) /(z[i] - z[i+1]) & ;
T[i_] := -(q z[i] - q^-1 z[i + 1]) Cancel @ Par[i] @ # + q # &

where I note that Factor is slower; Simplify and FullSimplify are much slower still and don't seem to remove all denominators. The above alternative now becomes

par[i_][expr_Plus] := par[i] /@ expr (* reduce to terms; the Expand avoids getting stuck with e.g. (z_i + q)(2 z_i +q) *) 
par[i_][term_] /; FreeQ[term, z[i]] && FreeQ[term, z[i+1]] := 0 (* reduce to terms involving z_i or z_{i+1} *)
par[i_][a___ b_ c___] /; FreeQ[b, z[i]] && FreeQ[b, z[i+1]] := b par[i][a c](* reduce to monomials in z_i,z_{i+1} *)
par[i_][term_Times]:= par[i] @ Expand @ term (* expand products still involving z_i,z_{i+1} *)
par[i_][z[j_]^m_. z[k_]^n_.] /; j==i && k==i+1 := Sign[m-n] (z[i] z[i+1])^Min[m, n] par[i][z[i]^Abs[m - n]](* remove factors of (z_i z_{i+1}) *)
par[i_][z[j_]^m_.] := (-1)^(i-j) Sum[z[i]^(m-n) z[i+1]^(n-1), {n, m}]

and

t[i_] := -(q z[i] - q^-1 z[i + 1]) par[i] @ # + q # &

Now I want to compute things like $T_{j-1} \cdots T_2 \, T_1 \, z_1^r$. On my laptop (Win 10 x64, Mma v11.3) I get

j = 10; r=8;
z[1]^r // Composition @@ Table[ T[k], {k,j-1,1,-1}] // AbsoluteTiming ;
z[1]^r // Composition @@ Table[ t[k], {k,j-1,1,-1}] // AbsoluteTiming ;
{%[[1]], %%[[1]], %[[2]] == %%[[2]] // Factor}
Clear[j,r]
(* Out:
{8.59344, 13.2408, True}
*)

Answer 1

Okay, here goes a wild attempt at using linear algebra to deal with the problem. The idea is to work with the CoefficientList of the polynomial instead of the symbolic expression. In addition, I wish to turn the CoefficientList of the polynomial in $x,y$, which will be a 2D List, into a vector representing the polynomial:

demoPoly = Sum[c[i, j] x^i y^j, {i, 0, 8}, {j, 0, 8}];
polyVector = Flatten[CoefficientList[demoPoly, {x, y}], 1];

The position of $c_{i, j}$ in the vector is (8 + 1)*i + j + 1.

Now, we know that if $m > n$, then $\partial(x^m y^n) = \sum_{i=1}^{m-n} x^{m-i} y^{n+i-1}$. We also know that if $m < n$, we get the same, except for a minus sign. Thus, we can find the matrix of the operator $\partial$ in the space spanned by monomials $x^m y^n$ with maximum exponent $\alpha-1$:

dMatrix[maxExpo_] := Block[
  {alpha, dim, mat, j, j2, sign, min, diff},
  alpha = maxExpo + 1;
  dim = alpha^2;
  mat = ConstantArray[0, {dim, dim}];
  Do[
   j = alpha m + n + 1;
   Do[
    mat[[alpha (m - i) + n + i, j]] = 1
    , {i, 1, m - n}
    ];
   j2 = alpha n + m + 1;
   mat[[All, j2]] = -mat[[All, j]];
   , {n, 0, maxExpo - 1}, {m, n + 1, maxExpo}
   ];
  mat
  ]

The nice thing about this matrix is that it can be applied to lots of different polynomial vectors as long as the largest exponent doesn't exceed $\alpha - 1$.

We can make a helper function that applies this to a symbolic polynomial and returns a symbolic polynomial:

fromCoefficientList[cl_] := Expand[Fold[FromDigits[Reverse[#1], #2] &, cl, {x, y}]

divDiff[maxExpo_][poly_] := fromCoefficientList[
  Partition[dMatrix[maxExpo].Flatten[
     PadRight[
      CoefficientList[poly, {x, y}], {maxExpo + 1, maxExpo + 1}]
     , 1], maxExpo + 1]]

The AbsoluteTiming of Expand@Cancel[Par@demoPoly] is

0.088310

and the AbsoluteTiming of divDiff[8][demoPoly] is only

0.006377

Let me know if this helps you at all! I thought it was a fun problem anyways :)