Low-order polylogarithms

[mikldk]

Hello,

I am trying to implement Low-order polylogarithms.

I have tried these rules:

SumFunc(_k, 1, _n, _k * (z_IsFreeOf(k)^k), z*( 1 - (n+1)*z^n + n*z^(n+1))/((1 - z)^2));
SumFunc(_k, 1, _n, _k^2 * (z_IsFreeOf(k)^k), z*( 1 + z - (n+1)^2*z^n + (2*n^2 + 2*n - 1)*z^(n+1))/((1 - z)^3));
SumFunc(_k, 1, _n, _k^_m * (z_IsFreeOf(k)^k), (D(z) Sum(k, 1, n, k^(m-1) * (z^k))) * z, IsPositiveInteger(m));

But the last rule seem to not match the lower order values of m because recursion does not seem to occur for m-1 etc. See example below, where I had hoped that the last line (Deriv(x,Sum(k,1,n,k^(m-1)*z^k))*x) would be processed further, but it seems like m-1 is not resolved to 2 and therefore the explicit rule for this is not matched.

Cheers, Mikkel.

This example:

x1 := Sum(k, 1, n, k * (x^k));
x2 := Sum(k, 1, n, k^2 * (x^k));
x3 := Sum(k, 1, n, k^3 * (x^k));

Print(" a ");
Print( Eval(x1) );
Print( Eval(x2) );
Print( Eval(x3) );

Print(" b ");
SumFunc(_k, 1, _n, _k * (z_IsFreeOf(k)^k), z*( 1 - (n+1)*z^n + n*z^(n+1))/((1 - z)^2));
SumFunc(_k, 1, _n, _k^1 * (z_IsFreeOf(k)^k), z*( 1 - (n+1)*z^n + n*z^(n+1))/((1 - z)^2));

Print(" c ");
Print( Eval(x1) );
Print( Eval(x2) );
Print( Eval(x3) );

Print(" d ");
SumFunc(_k, 1, _n, _k^2 * (z_IsFreeOf(k)^k), z*( 1 + z - (n+1)^2*z^n + (2*n^2 + 2*n - 1)*z^(n+1))/((1 - z)^3));

Print(" e ");
Print( Eval(x1) );
Print( Eval(x2) );
Print( Eval(x3) );

Print(" f ");
SumFunc(_k, 1, _n, _k^_m * (z_IsFreeOf(k)^k), (D(z) Sum(k, 1, n, k^(m-1) * (z^k))) * z, IsPositiveInteger(m));

Print(" g ");
Print( Eval(x1) );
Print( Eval(x2) );
Print( Eval(x3) );

Produces this:

" a "
Sum(k,1,n,k*x^k)
Sum(k,1,n,k^2*x^k)
Sum(k,1,n,k^3*x^k)
" b "
" c "
(x*(1-(n+1)*x^n+n*x^(n+1)))/(1-x)^2
Sum(k,1,n,k^2*x^k)
Sum(k,1,n,k^3*x^k)
" d "
" e "
(x*(1-(n+1)*x^n+n*x^(n+1)))/(1-x)^2
(x*(x-(n+1)^2*x^n+(2*n^2+2*n-1)*x^(n+1)+1))/(1-x)^3
Sum(k,1,n,k^3*x^k)
" f "
" g "
(x*(1-(n+1)*x^n+n*x^(n+1)))/(1-x)^2
(x*(x-(n+1)^2*x^n+(2*n^2+2*n-1)*x^(n+1)+1))/(1-x)^3
Deriv(x,Sum(k,1,n,k^(m-1)*z^k))*x

[grzegorzmazur]

Hi,

The issue stems from the fact that Sum() puts on hold (or prevents from evaluation if you like) the expression to be summed. Usually this is what you'd want, as it defers the evaluation to the point of the expression being actually summed. Here we'd like to evaluate the expression at least partially before passing it to sum. I'll try to come up with a workaround.

Cheers, Grzesiek

[mikldk]

That makes sense. Thanks for a good explanation.

I also tried with a For loop in the Sum body. But I couldn't get it working with multiple statements (temp var, For loop etc.). Would it be possible solving this way? Maybe with a helper function.

[mikldk]

I managed to get something working:

// _k apparantly different from _k^1 ....
SumFunc(_k, 1, _n, _k * (z_IsFreeOf(k)^k), z*( 1 - (n+1)*z^n + n*z^(n+1))/((1 - z)^2));

Function("testfunc", {n, z, m})
[
  res := (1 - z^(n+1)) / (1-z);
  For (i := 1, i <= m, i++) [
    res := (D(z) res) * z;
  ];
  res;
];
HoldArg(testfunc, z);

SumFunc(_k, 1, _n, _k^_m * (z_IsFreeOf(k)^k), testfunc(n, z, m), IsPositiveInteger(m));

Is it the same reason I need a SumFunc for _k * (z_IsFreeOf(k)^k); _k is not the same as _k^1?

[grzegorzmazur]

Hi,

Starting from the last question, basically yes. k^1 will immediately evaluate to just k, unless put on hold. Eg.

In> k^1
Out> k
In> Hold(k^1)
Out> k^1

The workaround you've found seems to be the best possible, I've tried to come up with something not involving a helper function but failed miserably. I've got only a few nitpicks:

1. I think that

   SumFunc(_k, 1, _n, _k * (z_IsFreeOf(k)^k), z*( 1 - (n+1)*z^n + n*z^(n+1))/((1 - z)^2));

   should read

   SumFunc(_k, 1, _n, _k * (z_IsFreeOf(k)^_k), z*( 1 - (n+1)*z^n + n*z^(n+1))/((1 - z)^2));

   and analogously in the other rule.

2. In the helper function res should be declared Local() (unfortunately, by default all variables are global); the slightly improved implementation could be eg

   Function("testfunc", {n, z, m})
   [
     Local(res);
     res := (1 - z^(n+1)) / (1-z);
     ForEach(i, 1 .. m)
       res := (Deriv(z) res) * z;
     res;
   ];

   (Deriv() instead of D() better handles some corner cases; D() is a macro which was implemented with ease of use in typical end-user cases in mind, implementing any yacas code you're better off with Deriv())

[mikldk]

@grzegorzmazur, thanks heaps. Would you be willing to consider inclusion of this in yacas as well as some of those I mention in #252, comment 4 Dec 2018 (is that case, I would try to make a pull request)? Or do you think they belong in end user code?

[grzegorzmazur]

Absolutely yes, please. While the yacas scripts library definitely could use some more structure and quite possibly stripping down a bit, the rules you've proposed (both here and in #252) are definitely good to have.

Thanks, Grzesiek

[mikldk]

Cheers.

To avoid to much duplication, I have been trying to make a general rule such that sums starting from 0 is first term for k = 0 and then the sum starting from 1. But I have not had much luck.

Pseudo code would look like this:

SumFunc(_k,0,_n,_expr, Eval(_expr, k = 0) + Sum(k, 1, n, expr) );

But that is not how Eval() works. Do you have a suggestions of how to make such a rule? Or would you prefer separate rules for sums starting from k=0 and k=1?

Cheers, Mikkel.

[grzegorzmazur]

Hi,

In principle this should not be necessary. Every SumFunc() defines Sum() from the given starting point p and, additionally, another Sum() from any integer q larger than the starting point. The latter Sum() is defined as difference between the original sum and Sum() from p to q. Hence, if your sum makes sense to be defined from 0, just do it and the sum starting from 1 should be automatically generated. From what I've just checked it seems to work for a couple of simple cases. Please let me know if it fails for you.

Cheers, Grzesiek

[mikldk]

Okay. I read it as: if the sum is defined if starting from 0 instead of 1, then SumFunc should define it starting from 0?

So the current rule

SumFunc(_k,0,_n,_k, n*(n+1)/2 );

I would redefine to

SumFunc(_k,0,_n,_k, n*(n+1)/2 );

I guess that's okay? Because then both

Sum(k,1,n,k)

and

Sum(k,0,n,k)

can be evaluated (which is not possible now).

Also: there is this chunk in sums/code.ys

// Some type of canonical form is needed so that these match when
// given in a different order, like x^k/k! vs. (1/k!)*x^k
// works !
SumFunc(_k,1,_n,_c + _d,
  Eval(UnList({Sum,sumvar'arg,1,n,c})) +
  Eval(UnList({Sum,sumvar'arg,1,n,d}))
);
SumFunc(_k,1,_n,_c*_expr,Eval(c*UnList({Sum,sumvar'arg,1,n,expr})), IsFreeOf(k,c) );
SumFunc(_k,1,_n,_expr/_c,Eval(UnList({Sum,sumvar'arg,1,n,expr})/c), IsFreeOf(k,c) );

Shouldn't all these also be converted to start from 0 instead of 1 then?

[grzegorzmazur]

Apart from some accidental typos (you've got 0 instead of 1 in some places in your post) this should be right. From some preliminary tests I've run previously it seems to work. Unfortunately, I've just managed to make a big mess of my working copy of yacas and can't test it more thoroughly now. But I believe it should just work after converting to start from 0 instead of 1.

[grzegorzmazur]

OK, I've run a couple of tests, and apparently there's something wrong with the difference mechanism in case of sums with additional condition. Eg

SumFunc(_k,0,_n,c_IsFreeOf(_k),c*(n+1));

works as expected:

In> Sum(k, 1, n, c)
Out> c*n
In> Sum(k, 0, n, c)
Out> c*(n+1)

while the original formulation (but with starting value changed to 0)

SumFunc(_k,0,_n,_c,c*(n+1),IsFreeOf(k,c) );

yields

In> Sum(k, 0, n, c)
Out> c*(n+1)
In> Sum(k, 1, n, c)
CommandLine(1) : Max evaluation stack depth reached.
Please use MaxEvalDepth to increase the stack size as needed.
CommandLine(1) : Max evaluation stack depth reached.
Please use MaxEvalDepth to increase the stack size as needed.

So while it seems that there's a way forward at least in the cases where condition can be moved to the expression itself, there's something going wrong with the external condition. I'll try to investigate it, but don't have a clue at the moment.

[grzegorzmazur]

OK, it's been a nasty typo. Should work now, see bb70c07d82354835dc837de30a3a6619c556224e

[mikldk]

Thanks! Good catch. Now, you identifying and fixing a bug, at least all my questions have not been a total waste :). I'll try to proceed.

[mikldk]

I now continue to try to implement the general low-order polylogarithms.

Just before returning the result of the helper function, I perform the following:

res := Simplify(res * (1-z)^(m+1)) / (1-z)^(m+1);

For m = 1 I get

(z*(z*n*z^n+z^(n+1)+1-z^(n+1)-z^n-n*z^n))/(1-z)^2

whereas Wikipedia gives

(z*(1-(n+1)*z^n+n*z^(n+1)))/(1-z)^2

According to Wolfram Alpha, the two are the same)%2F(1-z)%5E2+%3D+(z(1-(n%2B1)z%5En%2Bn*z%5E(n%2B1)))%2F(1-z)%5E2).

For m = 2 I get

(z*(2*z^2*n^2*z^(n-1)+z^2*n*z^n+2*z^2*n*z^(n-1)+z^(n+2)-(z^3*n*z^(n-1)+z^3*n^2*z^(n-1))+z-z^(n+2)-(z*n*z^(n-1)+z*n^2*z^(n-1))+1-z^(n+1)-z^n-n*z^n))/(1-z)^3

whereas Wikipedia gives

(z*(z-(n+1)^2*z^n+(2*n^2+2*n-1)*z^(n+1)-n^2*z^(n+2)+1))/(1-z)^3

I cannot get these into Wolfram Alpha.

But the expression I get out contains terms like z^2*n*z^n. I cannot understand why I don't get z^(n+2)*n here? I am hoping that by getting that, I would be able to return a simpler result.

[mikldk]

I have tried a local simplication rule in testfunc like this:

res := res /:: { 
    _z*_a*_z^_q <- a*z^(q+1), 
    _z^_p*_a*_z^_q <- a*z^(p+q) 
};

But it does not work. If I do this in the yacas program directly, it works:

In> { z*n*z^n, z^2*n*z^n } /:: { _z*_n*_z^_q <- n*z^(q+1), _z^_p*_n*_z^_q <- n*z^(p+q) };
Out> {n*z^(n+1),n*z^(n+2)}

I would think it has something to do with that it's inside a function, z is held (HoldArg()), I need to specify type of z such that * is commutative or something like that, but I cannot figure out exactly what's going on. I have tried replacing _z with z in the simplification rule without any luck.

[mikldk]

It is a bit confusing with z and n in the rules; if I change them to e.g. x and y I get the same as above:

In> e := { z*n*z^n, z^2*n*z^n, 2*z^2*n*z^n };
Out> {z*n*z^n,z^2*n*z^n,2*z^2*n*z^n}
In> r1 := { _x*_y*_x^_q <- y*x^(q+1), _x^_p*_y*_x^_q <- y*x^(p+q) };
Out> {_x*_y*_x^_q<-y*x^(q+1),_x^_p*_y*_x^_q<-y*x^(p+q)}
In> e /:: r1;
Out> {n*z^(n+1),n*z^(n+2),2*n*z^(n+2)}

Even with more terms:

In> e := { z*n*z^n+a, b+z^2*n*z^n, c-2*z^2*n*z^n };
Out> {z*n*z^n+a,b+z^2*n*z^n,c-2*z^2*n*z^n}
In> e /:: r1;
Out> {n*z^(n+1)+a,b+n*z^(n+2),c-2*n*z^(n+2)}

But applying the same rule in the helper function does not perform the same simplification.

[grzegorzmazur]

Let me start with the simple stuff - I've verified that your result for m=2 is exactly the same as the reference result from wikipedia (pen and napkin are really powerful tools ;)

With respect to z^2*n*z^n vs z^(n+2)*n, this is problem with getting multiplicative expressions to the canonical form. It should involve sorting the elements so the terms involving z are brought together and subsequently simplified. Unfortunately, this seems to be not that easy to implement correctly; the main issue is efficiency for large (thousands of terms) expressions. On the other hand the expressions you deal with are always univariate rational functions, so perhaps some simplified approach could be good enough. I'll try to come up with something.

I'll have a closer look at the rest of your comments and try to comment ASAP, but I'm afraid that it may take a few days (I'll be leaving for a conference next week and need to prepare to that).

[mikldk]

I agree :-); I simply focused on simplifying it first because I think it is a general and interesting problem, and this seemed like a nice place to start learning about it because of the simplicity.

No worries, this is in no way under time pressure, I just think it is a funny little problem and a good way to learn about CAS and yacas. You are also welcome to just give a few tips: I am thinking of trying to do something like this:

1. Expand() to transform the polynomial to expanded form
2. Split terms into a list, L
3. For each term T in L: collect factors of same variables (as well as constants) and simplify by multiplication

But I don't know yacas well enough to have any idea of how to do this (or something smarter). So it would be nice and interesting with help from you when you have time. Have a nice conference.