*
*   File contains all relevant routines for Mellin transform calculations.
*   First there are some declarations.
*
#include nndecl.h
CF  Li,Lix,logs,logx;
S   Mel,[d/dx],y1,y2,null;
#call tables(`SIZE')
*
*   Rule one: all powers of x are collected in a function pow(x,n) = x^n.
*   The program knows the following functions:
*       den:
*           den(1-x)    = 1/(1-x)    special meaning in makesum
*           den(1+x)    = 1/(1+x)    special meaning in makesum
*           den(anything) = 1/anything
*       logs:
*           logs(1,x,n) = ln_(x)^n
*           logs(2,x,n) = ln_(1-x)^n
*           logs(3,x,n) = ln_(1+x)^n
*       Li: these are polylogarithms
*           The following arguments are recognized:
*           1: x
*           2: -x
*           3: (1+x)/2
*           4: 1-x
*           5: 1/(1+x)
*           6: 2*x/(1+x)
*           7: (1-x)/(1+x)
*           8: -(1-x)/(1+x)
*
*           Li(n,narg,x,k) = (Li_n(narg(x)))^k
*               or the n-th polylogarithm of the narg-th argument of x
*               and that all to the power k.
*           ( Li_n(x) = sum_(j,1,inf,x^j/j^n) )
*           ( Li_1(x) = -ln_(1-x) )
*
*       We need: Values of the arguments in zero and in one.
*                These are in the tables arg0 and arg1.
*       Also for derivatives:
*       [d/dx]*Li(n?,arg?,x,1) = Li(n-1,arg,x,1)*dlnarg(arg)
*       with dlnarg(arg) = ([d/dx]*arg(x))/arg(x)
*               also these are collected in a table.
*       The fourth table is basically the expression for
*       Li(1,arg,x,1) which can be worked out because it is just
*           -ln_(1-arg(x))
*
*   Note that the tables can be extended rather easily
*
CTable dlnarg(1:8,x); * 1/arg(x)*([d/dx]*arg(x))
Fill dlnarg(1) = pow(x,-1);
Fill dlnarg(2) = pow(x,-1);
Fill dlnarg(3) = +den(1+x);
Fill dlnarg(4) = -den(1-x);
Fill dlnarg(5) = -den(1+x);
Fill dlnarg(6) = pow(x,-1)-den(1+x);
Fill dlnarg(7) = -den(1-x)-den(1+x);
Fill dlnarg(8) = -den(1-x)-den(1+x);
*
CTable li1(1:8,x); * -ln_(1-arg(x))
Fill li1(1) = -logs(2,x,1);
Fill li1(2) = -logs(3,x,1);
Fill li1(3) = -logs(2,x,1)+ln_(2);
Fill li1(4) = -logs(1,x,1);
Fill li1(5) = -logs(1,x,1)+logs(3,x,1);
Fill li1(6) = -logs(2,x,1)+logs(3,x,1);
Fill li1(7) = -logs(1,x,1)-ln_(2)+logs(3,x,1);
Fill li1(8) = -ln_(2)+logs(3,x,1);
*
CTable arg0(1:8);
Fill arg0(1) = 0;
Fill arg0(2) = 0;
Fill arg0(3) = 1/2;
Fill arg0(4) = 1;
Fill arg0(5) = 1;
Fill arg0(6) = 0;
Fill arg0(7) = 1;
Fill arg0(8) = -1;
*
CTable arg1(1:8);
Fill arg1(1) = 1;
Fill arg1(2) = -1;
Fill arg1(3) = 1;
Fill arg1(4) = 0;
Fill arg1(5) = 1/2;
Fill arg1(6) = 1;
Fill arg1(7) = 0;
Fill arg1(8) = 0;
*
#procedure derive
*
*   Make the derivative of a product of functions logs and Li
*
if ( count([d/dx],1) );
*
*   The derivatives of the logs should come first because when the
*   Li functions get derived the Li(1,...) can give more logs.
*   Alternatively they are substituted after the derivation proper.
*
    id  [d/dx]*logs(1,x,n?) = logs(1,x,n)*[d/dx]+n*logs(1,x,n-1)*pow(x,-1);
    id  [d/dx]*logs(2,x,n?) = logs(2,x,n)*[d/dx]-n*logs(2,x,n-1)*den(1-x);
    id  [d/dx]*logs(3,x,n?) = logs(3,x,n)*[d/dx]+n*logs(3,x,n-1)*den(1+x);
*
*   In the next loop we use an auxiliary function Lix to make sure
*   that everybody gets served only once.
*
    repeat;
        id  [d/dx]*Li(n?,iarg?,x,n1?) = [d/dx]*Lix(n,iarg,x,n1)
                    +n1*Lix(n,iarg,x,n1-1)*Lix(n-1,iarg,x,1)*dlnarg(iarg,x);
    endrepeat;
*
*   Finished with the derivation.
*   Lix can be put back. Special values can be substituted now
*   and powers of identical functions can be collected.
*
    id  [d/dx] = 0;
    Multiply replace_(Lix,Li);
    id  Li(?a,0) = 1;
    id  logs(?a,0) = 1;
    id  Li(1,n?,x,1) = li1(n,x);
    repeat id Li(?a,n1?)*Li(?a,n2?) = Li(?a,n1+n2);
    repeat id logs(?a,n1?)*logs(?a,n2?) = logs(?a,n1+n2);
    repeat id pow(x,n1?)*pow(x,n2?) = pow(x,n1+n2);
endif;
#endprocedure
*
#procedure fixvalue
*
*   Deals with the specific values of the arguments of the functions.
*   Note how everything is expressed in terms of the S functions,
*   also the divergencies.
*
id  pow(0,n?) = 0;
id  pow(1,n?) = 1;
*id  pow(x,0) = 1;
id  logs(1,1,n?) = 0;
id  logs(2,0,n?) = 0;
id  logs(3,0,n?) = 0;
id  logs(1,0,n?) = (-S(R(1),inf))^n;
id  logs(2,1,n?) = (-S(R(1),inf))^n;
id  logs(3,1,n?) = (-S(R(-1),inf))^n;
id  Li(n?,iarg?,0,n1?) = Li(n,arg0(iarg))^n1;
id  Li(n?,iarg?,1,n1?) = Li(n,arg1(iarg))^n1;
id  ln_(2)    = -S(R(-1),inf);
id  Li(n?,0) = 0;
id  Li(n?,1)  = S(R(n),inf);
id  Li(n?,-1) = S(R(-n),inf);
*
*   And S(R(-1,1,1,...,1),inf) = -Li_n(1/2) if there are n-1 indices 1
*       and one index -1 in front.
*       we use a trick to generate enough indices 1.
*
id  Li(n?,1/2) = -R(-1,n-1);
repeat id R(-1,?a,n?pos_) = R(-1,?a,1,n-1);
id  R(?a,0) = S(R(?a),inf);
#endprocedure
*
#procedure makesum
*
*   Replaces powers of 1/(1-x) or 1/(1+x) by their series expansion
*   in such a way that the sums nest properly.
*   The power of xsum tells which will be the next sum in the term.
*   All relevant quantities exist in sets and hence can be set up 
*   that way without further problems.
*   The replacement of null by zero is for a certain class of integrals.
*   This is explained in the file mel4.frm
*
while ( match(den(1-x)) );
    id,once,pow(x,n?)*den(1-x)*xsum^m? =
        xsum^m*xsum*sumset[m](jset[m],n,inf)*pow(x,jset[m]);
    multiply replace_(null,0);
endwhile;
while ( match(den(1+x)) );
    id,once,pow(x,n?)*den(1+x)*xsum^m? =
        xsum^m*xsum*sumset[m](jset[m],n,inf)*pow(x,jset[m])
                *sign(jset[m])*sign(n);
    multiply replace_(null,0);
endwhile;
*
#endprocedure
*
#procedure dosigns
*
*   Simplified version of the dosign routine.
*   Takes the sign function apart.
*
Splitarg,sign;
repeat id sign(n1?,n2?,?a) = sign(n1)*sign(n2,?a);
id  sign(x?number_) = sign_(x);
*
*   The next statement assumes that all arguments are integers.
*
id  sign(x?)*sign(x?) = 1;
#endprocedure
*
#procedure melsum(j)
*
*   Procedure does one summation step in the Mellin transforms
*   In this algorithm we first build up all the sums and then start
*   summing from the top down.
*   The sums here are of a rather simple nature. Hence there is no
*   need to call the general routines. Those are much slower because
*   they have to worry about general synchronization.
*
*   First make sure the signs are OK.
*
#call dosigns
*
*   Now separate the part of the argument that carries the j`j'.
*   This allows us to select only the proper functions den and S.
*   Next we mark these functions and we split fractions if possible.
*   (This may not be necessary)
*   After that we collect all powers of the relevant denominators
*
SplitArg,(j`j'),den,S;
id  den(x?,j`j') = den`j'(x,j`j');
repeat id den(x1?!{x2?},j`j')*den`j'(x2?!{x1?},j`j') =
            den(x2-x1)*(den(x1,j`j')-den`j'(x2,j`j'));
id  den`j'(x?,j`j') = den`j'(x,j`j',1);
repeat id den`j'(x?,j`j',n1?)*den`j'(x?,j`j',n2?) = den`j'(x,j`j',n1+n2);
id  S(R(?a),j`j') = S`j'(R(?a),0,j`j');
id  S(R(?a),x?,j`j') = S`j'(R(?a),x,j`j');
id  sign(j`j') = sign`j'(j`j');
*
*   Next is the only synchronization that can occur.
*
id  S`j'(R(n1?,?a),0,j`j')*den`j'(1,j`j',n?) =
        S`j'(R(n1,?a),1,j`j')*den`j'(1,j`j',n)
                -S`j'(R(?a),1,j`j')*den`j'(1,j`j',abs_(n1)+n)
                        *(-sign`j'(j`j'))^theta_(-n1);
id  S`j'(R,n?) = 1;
*
*   Make sure that this does not cause doubles. Pull them to the basis.
*
if ( count(S`j',1) > 1 );
    id  S`j'(x?,1,j`j') = S`j'(x,j`j'+1);
    #call basis(S`j')
    id  S`j'(x?,j`j'+1) = S`j'(x,1,j`j');
endif;
id  sign`j'(j`j')*sign`j'(j`j') = 1;
*
*   Now the summing
*   The if statements make sure that we only take what is needed
*   (technically they are not needed, but it makes things a bit safer)
*   If there is 'junk' it will be left in the output presumably
*   and it should be noticed.
*
if ( count(sum`j',1,den`j',1,S`j',1,sign`j',1) == 4 );
    id  sum`j'(j`j',n2?,n?)*den`j'(x?,j`j',n1?)*S`j'(R(?a),x?,j`j')*
        sign`j'(j`j') = (S(R(-n1,?a),n+x)-S(R(-n1,?a),n2+x-1))*sign(x);
elseif ( count(sum`j',1,den`j',1,S`j',1,sign`j',1) == 3 );
    id  sum`j'(j`j',n2?,n?)*den`j'(x?,j`j',n1?)*S`j'(R(?a),x?,j`j') =
            S(R(n1,?a),n+x)-S(R(n1,?a),n2+x-1);
    id  sum`j'(j`j',n2?,n?)*den`j'(x?,j`j',n1?)*sign`j'(j`j') =
            (S(R(-n1),n+x)-S(R(-n1),n2+x-1))*sign(x);
elseif ( count(sum`j',1,den`j',1,S`j',1,sign`j',1) == 2 );
    id  sum`j'(j`j',n2?,n?)*den`j'(x?,j`j',n1?) =
            S(R(n1),n+x)-S(R(n1),n2+x-1);
endif;
id  S(R,n?) = 1;
id  S(R(?a),n?neg0_) = 0;
id  S(x?,1+inf) = S(x,inf);
*
*   Collect the signs again
*
#call dosigns
*
*   Combine the leftover S functions to the basis
*
#call basis(S)
#endprocedure
*
#procedure melbreak
*
*   Procedure breaks down the functions by successive partial integration
*   of the powers of x in the function pow. Hence
*   int(pow(x,n)*f(x)) = (pow(x,n+1)*f(x)/(n+1))(0,1)
*                      - int(pow(x,n+1)*([d/dx]*f(x)))/(n+1)
*   The values at 0 and 1 can be substituted with the routine fixvalue.
*   Factors den(1-x) or den(1+x) can be expanded and give new sums.
*   When a Li becomes a Li(1,...) it can be simplified into
*   logarithms. The logarithms will give factors den(1-x) or den(1+x) also.
*   The 1/x from logs(1,x,n) go directly into the pow(x,k) function.
*   This will be repeated till nothing but a pow function is left.
*   The object Mel acts as the dx of the integral
*   xsum is used to keep track of how many sums there are in the term.
*   Its power tells what will be the next sum.
*
Multiply xsum;
#call makesum
*
*   $mcount tells in the output how far we are
*
#$mcount = 1;
*
*   The construction with the loop and the redefine of ibrk later
*   sets up a repeat that includes .sort instructions.
*   The program will remain in the loop as long as there are still
*   functions Li or logs and integration (Mel)
*
#do ibrk = 1,1
#ifdef `NOPOLE'
*
*   The quick way, but alas not always correct
*
    if ( match(logs(2,x,n?)) && ( count(logs,1,Li,1) == 1 ) );
        id  Mel*pow(x,n?)*logs(2,x,n1?) =
                sign_(n1)*fac_(n1)*den(n+1)*R(n1)*S(n+1);
        repeat id   R(n1?pos_,?a) = R(n1-1,1,?a);
        id  R(0,?a)*S(n?) = S(R(?a),n);
    endif;
    id  Mel*pow(x,n?) =
        den(n+1)*replace_(x,1)-Mel*pow(x,n+1)*den(n+1)*[d/dx];
#else
*
*   The proper treatment of powers of ln_(1-x). If the function has n1 of
*   these powers we use:
*       1: if that is all: a formula for Mel*pow(x,n)*ln_(1-x)^n1
*       2: if we have Mel*pow(x,n)*ln_(1-x)^n1*f(x) ->
*            Mel*ln_(1-x)^n1*pow(x,n)*(f(x)-f(1))+Mel*pow(x,n)*ln_(1-x)^n1*f(1)
*            The last term we know and the first one is finite.
*            This means that if during their calculation divergencies develop
*            they are identical divergencies and should cancel cleanly.
*
*       Note that due to the pow(x,n+1) the partial integration never
*       causes problems at x=0.
*
    if ( match(logs(2,x,n?)) );
        if ( count(logs,1,Li,1) == 1 );
            id  Mel*pow(x,n?)*logs(2,x,n1?) =
                    sign_(n1)*fac_(n1)*den(n+1)*R(n1)*S(n+1);
        else;
            id  Mel*pow(x,n?)*logs(2,x,n1?) =
                    -Mel*den(n+1)*pow(n+1)*[d/dx]*logs(n1)*(1-replace_(x,1))
                    +sign_(n1)*fac_(n1)*den(n+1)*R(n1)*S(n+1)*replace_(x,1);
            id  logs(n?) = logs(2,x,n);
            id  pow(n?) = pow(x,n);
        endif;
        repeat id   R(n1?pos_,?a) = R(n1-1,1,?a);
        id  R(0,?a)*S(n?) = S(R(?a),n);
    else;
        id  Mel*pow(x,n?) =
                den(n+1)*replace_(x,1)-Mel*pow(x,n+1)*den(n+1)*[d/dx];
    endif;
#endif
*
*   Substitute the values at 0 and 1.
*
    #call fixvalue
*
*   Do the derivation
*
    #call derive
    id  pow(x,0) = 1;
    .sort: derive-`$mcount';
*
*   set the loop parameter down if we are not yet finished
*
    if ( count(logs,1,Li,1) && count(Mel,1) ) redefine ibrk "0";
*
*   Expand the den(1-x) and the den(1+x)
*
    #call makesum
*
*   and clean up the signs
*
    #call dosigns
    .sort: makesum-`$mcount';
*
*   Raise the counter once for the next pass in the loop.
*
#$mcount = $mcount + 1;
#enddo
*
*   If there is only a power of x left we can do the transform easily:
*
id  Mel*pow(x,n?) = den(n+1);
id  pow(x,0) = 1;
*
*   Lower xsum once. Now it indicates the number of sums in each term.
*
Multiply 1/xsum;
#endprocedure
#procedure melsum2
*
*   Other algorithm to work out a mellin-type sum/integral
*   This method is specifically suitable for integrals.
*   Each time we have two sums we exchange the order and we can
*   do one of the sums. This way there will never be a build-up of
*   nested sums.
*   It is much faster for the integrals and slower for the Mellin transforms
*   (when it is addapted for those).
*   Hence it should only be used for the integrals.
*
if ( match(pow(x,n?)) == 0 ) Multiply pow(x,0);
#call makesum
#$mcount = 1;
*
*   Set up an indefinite loop
*
#do i = 1,1
*
*   Do the partial integration. See ln_(1-x) as a special object
*
    if ( match(logs(2,x,n?)) );
        if ( count(logs,1,Li,1) == 1 );
            id  Mel*pow(x,n?)*logs(2,x,n1?) =
                    sign_(n1)*fac_(n1)*den(n+1)*R(n1)*S(n+1);
        else;
            id  Mel*pow(x,n?)*logs(2,x,n1?) =
                    -Mel*den(n+1)*pow(n+1)*[d/dx]*logs(n1)*(1-replace_(x,1))
                    +sign_(n1)*fac_(n1)*den(n+1)*R(n1)*S(n+1)*replace_(x,1);
            id  logs(n?) = logs(2,x,n);
            id  pow(n?) = pow(x,n);
        endif;
        repeat id   R(n1?pos_,?a) = R(n1-1,1,?a);
        id  R(0,?a)*S(n?) = S(R(?a),n);
    else;
        id  Mel*pow(x,n?) =
                den(n+1)*replace_(x,1)-Mel*pow(x,n+1)*den(n+1)*[d/dx];
    endif;
    id  den(n?number_) = 1/n;
    .sort: partial-1-`$mcount';
    #call fixvalue
    #call derive
    .sort: partial-2-`$mcount';
*
*   Expand den(1-x) and den(1+x)
*
    #call makesum
    #call dosign
    id  sum1(j1,1,inf) = sum1(j1,0,inf) - replace_(j1,0);
    id  den(n?number_) = 1/n;
    id  sign(n?number_) = sign_(n);
    .sort: partial-3-`$mcount';
*
*   If there are two sums, exchange them.
*
if ( count(sum2,1) );
    id  sum1(j1,0,inf)*sum2(j2,j1+1,inf) = sum2(j2,1,inf)*sum1(j1,0,j2-1);
*
*   Now do the old outermost sum.
*
    #call melsum(1)
    Multiply 1/xsum;
    id  S(R(n1?,?a),j2) = S(R(n1,?a),j2+1)
                    -S(R(?a),j2+1)*den(j2+1)^abs_(n1)*(-sign(j2))^theta_(-n1);
    id  S(R,n?) = 1;
*
*   Rename the 2 variable to the 1 variable
*
    Multiply replace_(sum2,sum1,j2,j1);
    id  sum1(j1,1,inf) = sum1(j1,0,inf)-replace_(j1,0)/xsum;
    id  den(1) = 1;
    id  sign(0) = 1;
    #call subesses(S)
endif;
*
*   See whether we should continue the loop
*
if ( count(logs,1,Li,1) || match(den(1+x)) || match(den(1-x)) ) redefine i "0";
.sort: sum-`$mcount';
*
*   If we are going to put in values, we may as well do it now
*
#ifdef `SUBNUM'
  id  S(R(?a),inf) = SS(R(?a),inf);
  id  SS(R(?a),inf) = SS(?a);
  repeat id SS(?a,n?{2,...,`MAXWEIGHT'},?b) = SS(?a,0,n-1,?b);
  repeat id SS(?a,n?{<-2>,...,<-`MAXWEIGHT'>},?b) = SS(?a,0,n+1,?b);
  #do k = 1,`SIZE'
    id  SS(n1?,...,n`k'?) = Stab`k'(n1,...,n`k');
  #enddo
  .sort:subnum-`$mcount';
#endif
#$mcount = $mcount +1;
*
*   End of loop
*
#enddo
id  pow(x,n?)*Mel = den(n+1);
.sort:Final integral;
*
*   Now do the left over sum.
*
#call melsum(1)
id  xsum = 1;
#call dosign
id  S(R(?a),inf+1) = S(R(?a),inf);
.sort:Last sum;
*
*   Cleaning up if needed
*
#ifdef `SUBNUM'
  id  S(R(?a),inf) = SS(R(?a),inf);
  id  SS(R(?a),inf) = SS(?a);
  repeat id SS(?a,n?{2,...,`MAXWEIGHT'},?b) = SS(?a,0,n-1,?b);
  repeat id SS(?a,n?{<-2>,...,<-`MAXWEIGHT'>},?b) = SS(?a,0,n+1,?b);
  #do k = 1,`SIZE'
    id  SS(n1?,...,n`k'?) = Stab`k'(n1,...,n`k');
  #enddo
#else
  #call basis(S)
#endif
*
*   And substitute values for fixed arguments
*
#call subesses(S)
.sort: sub S values;
#endprocedure