(* ::Package:: *)

(************************************************************************)
(* This file was generated automatically by the Mathematica front end.  *)
(* It contains Initialization cells from a Notebook file, which         *)
(* typically will have the same name as this file except ending in      *)
(* ".nb" instead of ".m".                                               *)
(*                                                                      *)
(* This file is intended to be loaded into the Mathematica kernel using *)
(* the package loading commands Get or Needs.  Doing so is equivalent   *)
(* to using the Evaluate Initialization Cells menu command in the front *)
(* end.                                                                 *)
(*                                                                      *)
(* DO NOT EDIT THIS FILE.  This entire file is regenerated              *)
(* automatically each time the parent Notebook file is saved in the     *)
(* Mathematica front end.  Any changes you make to this file will be    *)
(* overwritten.                                                         *)
(************************************************************************)



BeginPackage["Tornheim`"]


TornheimT::usage="TornheimT[m,k,n] is the Tornheim sum.
N[T[m,k,n], p] computes it with p digits of precision.";
TornheimY::usage="TornheimY[m,n] is the integral Y[m,n]*.
N[Y[m,k,n], p] computes it with p digits of precision.";


TornheimTriples::usage="TornheimTriples[w] gives a list of all admissible triples of (integer) weight w >= 3.";


StepReduce::usage="StepReduce[expr] applies to any Tornheim sum TornheimT[m,k,n] appearing in expr the basis identity T[m,k,n]=T[k,m-1,n+1]+T[m,k-1,n+1].";


HuardReduce::usage="HuardReduce[expr] expresses every Tornheim sum TornheimT[m,k,n] appearing in expr in terms of sums of the type TornheimT[i,0,j].";


IrreducibleTornheimBasis::usage="IrreducibleTornheimBasis[w] gives the irreducible basis of Tornheim sums corresponding to weight w.";


ComputeIrreducibleTornheimBasis::usage="IrreducibleTornheimBasis[w] computes the irreducible basis of Tornheim sums corresponding to weight w.";


ToIrreducibleTornheimBasis::usage="ToIrreducibleTornheimBasis[w] solves for all the sums T[m,0,n] of weight w in terms of zeta values and the corresponding irreducible basis.";


CompleteReduceToBasis::usage="CompleteReduceToBasis[expr] expresses every Tornheim sum TornheimT[m,k,n] appearing in expr in terms of the irreducible basis of sums of the type TornheimT[i,0,j].";


CompleteReduceToY::usage="CompleteReduceToY[expr] expresses every Tornheim sum TornheimT[m,k,n] appearing in expr in terms of the Y[i,j] corresponding to the irreducible basis.";


ClausenCl::usage="ClausenCl[n,q] is the Clausen function of order n.";


NSumTornheimT::usage="NSumTornheimT[m,k,n, NSum options] computes the numerical value of the Tornheim sum TornheimT[m,k,n] directly using NSum, with the given options.";


zeta::usage="zeta[n] symbolizes the unevaluated value Zeta[n], for n > 0."


Begin["`Private`"]


(* ::Text:: *)
(*Requires Combinatorica package*)


<<Combinatorica`;


(* ::Text:: *)
(*To generate all admissible triples for a given weight*)


(* ::Input:: *)
SetAttributes[TornheimTriples,ReadProtected];

TornheimTriples[w_?IntegerQ/;w>=3]:=
Select[Compositions[w,3],(#[[1]]>=#[[2]])&&(#[[1]]+#[[3]]>1)&&(#[[2]]+#[[3]]>1)&&(#[[1]]+#[[2]]+#[[3]]>2)&];


(* ::Input:: *)
SetAttributes[StepReduce,ReadProtected];

StepReduce[expr_]:=expr/.TornheimT[m_,k_,n_]/;m>=1&&k>=1->TornheimT[k,m-1,n+1]+TornheimT[m,k-1,n+1];


(* ::Text:: *)
(*To apply Huard' s reduction to an arbitrary T (m, k, n)*)


(* ::DisplayFormula:: *)
(*Usage: HuardReduce[expr_]*)


(* ::Text:: *)
(*where expr_ is any expression containing Tornheim sums T (m, k, n)*)


(* ::Input:: *)
SetAttributes[HuardReduce,ReadProtected];

HuardReduce[expr_]:=expr/.{TornheimT[m_Integer,k_Integer,n_]/;k>=1->With[{w=m+k+n},Sum[Binomial[m+k-i-1,m-i]TornheimT[i,0,w-i],{i,1,m}]+Sum[Binomial[m+k-i-1,k-i]TornheimT[i,0,w-i],{i,1,k}]]}


(* ::Text:: *)
(*Some known evaluations used to solve in terms of the irreducible T' s*)


myTeZ[m_,k_,0]:=Zeta[m]Zeta[k];
myTeZ[1,0,n_/;n>=2]:=1/2 (n Zeta[n+1]-Sum[Zeta[i]Zeta[n+1-i],{i,2,n-1}]);
simTZ[m_,n_]:=Zeta[m]Zeta[n]-Zeta[m+n];


myTe[m_Integer,0,n_Integer]/;OddQ[m+n]&&m>=1&&n>=2:=With[{w=m+n},(-1)^m Sum[Binomial[w-2j-1,m-1]zeta[2j]zeta[w-2j],{j,0,Floor[(n-1)/2]}]+(-1)^m Sum[Binomial[w-2j-1,n-1]zeta[2j]zeta[w-2j],{j,0,Floor[m/2]}]-1/2 zeta[w]]


(* ::Text:: *)
(*Function that gives the list of the T (m, 0, n) in the irreducible basis*)


(* ::Input:: *)
SetAttributes[IrreducibleTornheimBasis,ReadProtected];

IrreducibleTornheimBasis[w_?EvenQ]/;w>=3:=Table[TornheimT[w-2k,0,2k],{k,1,Floor[(w-2)/6]}];
IrreducibleTornheimBasis[w_?OddQ]/;w>=3:={};


(* ::Input:: *)
SetAttributes[ComputeIrreducibleTornheimBasis,ReadProtected];

ComputeIrreducibleTornheimBasis[w_?EvenQ]/;w>=4:=With[{allTs=ToIrreducibleTornheimBasis[w][[All,2]]},DeleteCases[Map[{#,FreeQ[allTs,#]}&,Table[TornheimT[m,0,w-m],{m,2,w-2}]],{_,True}]][[All,1]];

ComputeIrreducibleTornheimBasis[w_?OddQ]/;w>=3:={};


(* ::Text:: *)
(*Function to solve for all the T (m, 0, n) in terms of the irreducible basis*)


(* ::Input:: *)
SetAttributes[ToIrreducibleTornheimBasis,{NHoldAll,ReadProtected}];

Clear[ToIrreducibleTornheimBasis];

ToIrreducibleTornheimBasis[w_?EvenQ]/;w>=4:=
ToIrreducibleTornheimBasis[w]=Block[{pairs=Table[{w-k,k},{k,2,w/2}],reps={}},
Off[Solve::svars];
reps=Solve[Join[Map[HuardReduce[TornheimT[#[[1]],#[[2]],0]]==myTeZ[#[[1]],#[[2]],0]&,pairs],Map[TornheimT[#[[1]],0,#[[2]]]+TornheimT[#[[2]],0,#[[1]]]==simTZ[#[[1]],#[[2]]]&,pairs]
]/.{TornheimT[1,0,w-1]->myTeZ[1,0,w-1]},Table[TornheimT[k,0,w-k],{k,2,w-2}]
];
On[Solve::svars];
Prepend[reps[[1]],TornheimT[1,0,w-1]->myTeZ[1,0,w-1]]//Expand
]/.Tozeta;

ToIrreducibleTornheimBasis[w_?OddQ]/;w>=3:=ToIrreducibleTornheimBasis[w]=Table[TornheimT[j,0,w-j]->myTe[j,0,w-j],{j,1,w-2}];


(* ::Text:: *)
(*Functions to convert Zeta' s and \[Pi]' s to symbolic "zeta"*)


Tozeta={\[Pi]^n_?EvenQ->((-1)^(n/2+1) n!)/(2^(n-1) BernoulliB[n]) zeta[n],Zeta[n_]->zeta[n]};
zetazetaTozeta:={zeta[m_?EvenQ]zeta[n_?EvenQ]->(Zeta[m]Zeta[n])/Zeta[m+n] zeta[m+n]};
zeta[0]:=Zeta[0];


(* ::Input:: *)
SetAttributes[CompleteReduceToBasis,{NHoldAll,ReadProtected}];

CompleteReduceToBasis[expr_]:=(expr/.{TornheimT[0,0,n_]/;n>=3->zeta[n-1]-zeta[n],
TornheimT[m_,k_,n_]/;m+k>=1:>(HuardReduce[TornheimT[m,k,n]]/.ToIrreducibleTornheimBasis[m+k+n])})//Expand


(* ::Input:: *)
SetAttributes[CompleteReduceToY,{NHoldAll,ReadProtected}];

CompleteReduceToY[expr_]:=(expr/.TornheimT[m_,k_,n_]:> CompleteReduceToBasis[TornheimT[m,k,n]]/.TtoY/.zetazetaTozeta)//Expand


(* ::Text:: *)
(*Replacement rule to evaluate Tornheim sums T (m, 0, n) of even weight in terms of the Y*-integrals*)


TtoY:={TornheimT[m_?EvenQ/;m>=2,0,n_?EvenQ/;n>=2]->myT1[m,n]+myT2[n,m],TornheimT[m_?OddQ/;m>=2,0,n_?OddQ/;n>=2]->myT1[m,n]+myT2[m,n]}


myT1[m_,n_]:=zeta[m]zeta[n]-1/2 zeta[m+n];
myT2[m_,n_]:=-Sum[Binomial[m+n-2k-2,m-1]zeta[m+n-2k-1]zeta[2k+1],{k,1,(m+n)/2-2}]-(-1)^((m+n)/2) TornheimY[m,n];


(* ::Text:: *)
(*Explicit expression for Y in terms of the X - integrals*)


Ye[m_,n_]:=With[{w=m+n},(2(2\[Pi])^(w-2))/(m!(n-2)!) Sum[(-1)^k Binomial[m,k]X[k,w-2-k],{k,0,m}]+(-1)^(w/2-1) Binomial[w-2,m-1]Log[2\[Pi]]zeta[w-1]];


Xe[m_,n_]:=With[{n1=n+1},(-1)^Floor[n/2] n!/(2\[Pi])^n Defer[Integrate[LogGamma[q]BernoulliB[m,q]ClausenCl[n1,2\[Pi] q],{q,0,1}]]];


(* ::Section:: *)
(*Numerical Evaluations*)


(* ::Input:: *)
SetAttributes[NSumTornheimT,ReadProtected];

NSumTornheimT[m_,k_,n_,opts:OptionsPattern[]]:=NSum[1/(i^m j^k (i+j)^n),{i,1,\[Infinity]},{j,1,\[Infinity]},Evaluate[FilterRules[{opts}, Options[NSum]]]]


(* ::Input:: *)
SetAttributes[ClausenCl,ReadProtected];
ClausenCl[n_?EvenQ,q_]:=Re[1/2 \[ImaginaryI] (PolyLog[n,\[ExponentialE]^(-\[ImaginaryI] q)]-PolyLog[n,\[ExponentialE]^(\[ImaginaryI] q)])];
ClausenCl[n_?OddQ,q_]:=Re[1/2 (PolyLog[n,\[ExponentialE]^(-\[ImaginaryI] q)]+PolyLog[n,\[ExponentialE]^(\[ImaginaryI] q)])];


(* ::Input:: *)
SetAttributes[zeta,ReadProtected];

zeta/:N[zeta[n_]]:=N[Zeta[n]];
zeta/:N[zeta[n_],p_]:=N[Zeta[n],p];
zeta[0]=Zeta[0];


yy[m_,n_][q_]:=With[{w=m+n},2Sum[(-1)^Floor[(w-2+j)/2] (2\[Pi])^j/j! Binomial[w-2-j,n-2]BernoulliB[j,q]ClausenCl[w-1-j,2\[Pi] q],{j,1,m}]];

y2[m_,n_]:=With[{w=m+n},1/2 (-1)^((w-2)/2) Binomial[w-2,m]zeta[w]+(-1)^((w-2)/2) Binomial[w-2,m-1]zeta[w-1]Log[2]];


K[q_]:=-Log[Sin[\[Pi] q]];
yyK[m_,n_,p_:6]:=1/2 NIntegrate[yy[m,n][q]K[q],{q,0,1},WorkingPrecision->Max[p+8,$MachinePrecision],AccuracyGoal->p];


(* ::Input:: *)
SetAttributes[{TornheimT,TornheimY},{NHoldAll,ReadProtected}];
SetAttributes[{Tny,YnK,yyK,y2,TYval},NHoldAll];


(* ::Input:: *)
TornheimT/:N[TornheimT[m_,k_,n_],p_]:=Tny[m,k,n,p];
TornheimT/:N[TornheimT[m_,k_,n_]]:=Tny[m,k,n];


(* ::Input:: *)
TornheimY/:N[TornheimY[m_,n_],p_]:=YnK[m,n,p];
TornheimY/:N[TornheimY[m_,n_]]:=YnK[m,n];


(* ::Input:: *)
N[TornheimY[m_,n_],p_]:=YnK[m,n,p];
N[TornheimY[m_,n_]]:=YnK[m,n];


Tny[m_?IntegerQ,k_?IntegerQ,n_?IntegerQ,p_:6]:=Module[{TofY,tny},
TofY=CompleteReduceToY[TornheimT[m,k,n]]/.zeta->Zeta;
tny=TofY/.TornheimY[i_,j_]:>YnK[i,j,p+2];
N[tny,p]
]


YnK[m_,n_,p_:6]:=Module[{ty,pd},
pd=If[Head[p]===List,p[[1]],p];
(*Print["Function YnK called, with m=",m,",n=",n," and p=",pd];*)
If[ValueQ[TYval[m,n]]&&(pd<=TYval[m,n][[1]]),
(*Print["Old value used"];*)
ty=N[TYval[m,n][[2]],pd],
(*Print["New value computed"];*)
ty=N[y2[m,n],pd]+yyK[m,n,pd];TYval[m,n]={pd,ty}];
TYval[m,n][[2]]];


End[]


EndPackage[]