(* ::Package:: *)

BeginPackage["pfq`"]

Print["*-*-*-*-*-* pfq 1.1 for Mathematica 6.0 *-*-*-*-*-*\n"]

Print["*-*-*-*-*-*Powered By V. Bytev, M. Kalmykov, B. Kniehl, \
01.03.2012*-*-*-*-*-*\n"]


explicitForm::usage=" put the answer of simplifyPFQ[] for \
Mathematica numerically calcualble way "
ToGroebnerBasis::usage="the only prameters are the  parameters of \
Hypergeometric function, returns the operator of differential
  basis, final parameters of Hypergeometric function and overall factor"

Begin["`Private`"]





simplifyPFQAiUpMOD2Matrix[vector_List,index_?NumberQ]:=Module[
{p=Length[vector[[2]]]},
{
{
Table[Sum[SymmetricPolynomial[p-j,Map[(#-1)&,vector[[2]]]] (-vector[[1,index]])^(j-phi),{j,phi,p}],{phi,0,p}],
Table[-SymmetricPolynomial[p-i,Delete[vector[[1]],index]],{i,0,p}]
},
ReplacePart[vector,vector[[1,index]]+1,{1,index}],
1/Apply[Times,Map[(#-vector[[1,index]]-1)&,vector[[2]]]],
1
}
];


simplifyPFQBiDownMOD2Matrix[vector_List,index_?NumberQ]:=Module[
{p=Length[vector[[2]]]},
{
{
Table[-Sum[SymmetricPolynomial[p-j,vector[[1]]](1-vector[[2,index]])^(j-phi),{j,phi,p}],{phi,0,p}  ],
Prepend[Table[ SymmetricPolynomial[p-i, If[ p==1,{1},Map[ (#-1)&,Delete[vector[[2]],index]]]],{i,1,p}],0]
},
ReplacePart[vector,vector[[2,index]]-1,{2,index}],
(vector[[2,index]]-1)/Apply[Times,Map[(1+#-vector[[2,index]])&,vector[[1]]]],
-1
}
];



simplifyPFQAiUpBiDownMatrix[vector_List,index_?NumberQ]:=Module[{},
Switch[index
,_?(#<0&),{   {    Table[0,{i,1,Length[vector[[1]]]}],PadRight[Insert[Table[SymmetricPolynomial[-j,Table[i,{i,index+1,0}]],{j,index+1,0}],0,1],Length[vector[[1]]],0]   },
                                   ReplacePart[vector+index,vector[[3]],3],Times@@Map[Pochhammer[#,-index]&,vector[[2]]+index]/Times@@Map[Pochhammer[#,-index]&,vector[[1]]+index]  , index
                              }
,_?(#<=Length[vector[[1]]]&),simplifyPFQAiUpMOD2Matrix[vector,index]
,_,simplifyPFQBiDownMOD2Matrix[vector,index-Length[vector[[1]]]]
]
];





permutationMatrix[list_List,variable_,power_?NumberQ]:=Module[
{
lengthList=Length[list[[1]]],
lengthtotal=Length[list]
},
Table[ If[(i-j+1)<=0||((i-j+1)>lengthList),0,list[[1,i-j+1]]]+Sum[       If[(i-k+1)<=0||((i-k+1)>lengthList)||lengthtotal==1,0,variable^power power^(j-k) Binomial[j-1,k-1] list[[2,i-k+1]]]   ,{k,1,j}  ],{j,1,lengthList},{i,1,2lengthList-1}]

];



DerivativeForPFQPowerPhi[list_List,phi_?NumberQ,l_?NumberQ]:=Module[
{p=Length[list[[2]]]},
(
listofBim1=Map[(#-1)&,list[[2]]];
list[[3]]  Sum[Binomial[phi,l-j]SymmetricPolynomial[p-j+1,list[[1]]],{j,If[(l-phi)>=0,l-phi,0],If[l<=p+1,l,p+1]}]-If[l>=phi+1,SymmetricPolynomial[p-l+phi+1,listofBim1],0]
)
];

DerivativeForPFQMatrix[list_List]:=Module[
{p=Length[list[[2]]]},
(
Join[Table[If[i==j,1,0],{j,1,p+1},{i,1,p+1}],Nest[Append[#,(phi=Length[#];   Table[DerivativeForPFQPowerPhi[list,phi,l                  ],{l,0,p    }]          +Apply[Plus,Table[#[[i]] DerivativeForPFQPowerPhi[list,phi,p+i],{i,1,phi}]]           )/(1-list[[3]])   ]&,{},p]]
)
];


(*--------------------------------------------------------------------------*)
(*  Проверка на приавильность заполнения вектора гипергеом. функции и вектора сдвижки  *)

listCheckingVector[vector_]:=Module[{},
{
simplifyPFQ::wrongParameters="wrong parameters in `1` vector";
If[(Length[vector[[1]]]!=Length[vector[[2]]]+1)||(Length[vector]!=3)||!FreeQ[vector,Null]||!FreeQ[vector[[1]],_List,1]||!FreeQ[vector[[2]],_List,1]||!FreeQ[vector[[3]],_List],{Message[simplifyPFQ::wrongParameters,vector];Abort[]}];
}];
listCheckingchangeList[vector_,changeList_]:=Module[{},
{
simplifyPFQ::wrongParameters="wrong parameters in `1` vector";
If[Length[changeList]!=2Length[vector[[1]]]-1||!ArrayQ[changeList,_,IntegerQ],{Message[simplifyPFQ::wrongParameters,changeList];Abort[]}];
}];

 listChecking[vector_,changeList_]:={listCheckingVector[vector];listCheckingchangeList[vector,changeList]};


(*--------------------------------------------------------------------------*)
(*Делает из вектора сдвижки набор, где на первом месте номер сдвигаемой переменной, на втором куда ее сдвигать (+ или - число). Нули отбрасываются
listOfseparatechanges[{2,1,-3}]------>{{1,2},{1,2},{2,1},{3,-3},{3,-3},{3,-3}}*)

listOfseparatechanges[vector_List]:=Module[{},
Partition[Flatten[Map[Table[{#[[1]],#[[2]]},{i,1,Abs[#[[2]]]}]&,DeleteCases[MapIndexed[If[#1!=0,{#2[[1]],#1}]&,vector],Null]]],2]
];


(*--------------------------------------------------------------------------*)
(* Выделяет из вектора сдвижки номера парамеров,которые сдвигаются вверх (Ai) или  вниз (Bi)
listOfApandBdchanges[{2,-3,4,6,-3},0]--------->{1,1,3,3,3,3,5,5,5} 
Отдельное число показвыает "общее" число сдвижек и меняет вектор сдвижек и добавляет в конец ответа нули, кот потом будут использоваься для 
вызова соотв. оператора:  listOfApandBdchanges[{-1,-1,-1,-1,-3},-1]-------->{5,5,-1}
  *)


listOfApandBdchanges[vector_List,number_]:=Module[{listOfchangesseparate= listOfseparatechanges[vector-number],total=(Length[vector]-1)/2},

Join[
Map[#[[1]]&,Select[listOfchangesseparate,(        #[[1]]<=(Length[vector]+1)/2&&#[[2]]>0
                                                                                          ||(#[[1]]>(Length[vector]+1)/2&&#[[2]]<0)
                                                                                          )&]]
,Table[-total,{i,1,IntegerPart[-number/total]}]
,If[Mod[-number,total] total!=0,{-Mod[-number,total]},{}]
]
];


(*--------------------------------------------------------------------------*)
(*Получене конечного вектора:   начальный вектор+ сдвижка    *)
vectorAfterChange[vector_List,changeVector_List]:=Module[{},
{
MapIndexed[(#1+changeVector[[      #2[[1]]                                                   ]])&,vector[[1]]],
MapIndexed[(#1+changeVector[[      #2[[1]]+Length[vector[[1]]]  ]])&,vector[[2]]],
vector[[3]]
}];

(*--------------------------------------------------------------------------*)
(* Получение списка сдвижек для Ai down и Biup операторов (простые). Это просто список с номерами сдвигаемых параметров
listOfAdownAndBupchanges[{2,-2,0,-2,2}]----------->{2,2,5,5} *)
listOfAdownAndBupchanges[changeList_List]:=Module[
{lengthAi=(Length[changeList]+1)/2},
Flatten[DeleteCases[MapIndexed[(If[(#1<0&&#2[[1]]<=lengthAi)||(#1>0&&#2[[1]]>lengthAi),Table[#2,{i,1,Abs[#1]}]])&,changeList],Null]]
];


(*--------------------------------------------------------------------------*)
(*Получение вектора после сдвижки и фактора  для  Ai down или Biup (простые операторы), номер сдвигаемой переменной идет вторым параметром: AdBupvector[vector,5]---->{b2,{{a1,a2,a3},{b1,1+b2},x}}
*)
AdBupvector[vector_List,number_]:=Module[{lengthAi=Length[vector[[1]]]},
If[number<=lengthAi,{vector[[1,number                 ]]-1,ReplacePart[vector,vector[[1,number                       ]]-1,{1,number                      }]},
                                              {vector[[2,number-lengthAi]],ReplacePart[vector,vector[[2,number-lengthAi]]+1,{2,number-lengthAi}]}
  ]
];


(*--------------------------------------------------------------------------*)
(* Выдает число, сообтветсвующее общей сдвижке всех параметров, второе число дает макс. длину вектора из theta
numberOfTotalChangings[{0,0,-2,-2,-2,-1,-1}]---->-1
*)
numberOfTotalChangings[changeList_List]:=Module[
{n=(Length[changeList]+1)/2,
orderedMinus=Sort[Select[changeList,#<0&],Greater],
numberOfchange
},
numberOfchange=((n-1)^2-1)/n-Apply[Plus,Map[(UnitStep[#])&,changeList]];
Switch[numberOfchange
,_?(#<0&),0
,_?(#>=0&&#<=Length[orderedMinus]-1&),orderedMinus[[numberOfchange+1]]
,_,spisok[[-1]]
]

];

(*--------------------------------------------------------------------------*)
(*Работа с одиноковыми переменными в Ai и Bi  
vector={{4,2,2,3,7},{2,2,5,2}};
changevector={1,2,3,4,5,6,7,8,9};
sameParRemove[vector,changevector]---->{{{4,3,7},{5,2}},{1,4,5,8,9}}
sameNumber[vector]----->{2,1}
 *)

sameNumber[vector_]:=Module[{},(Table[If[vector[[1,i]]==vector[[2,j]],Return[{i,j},Module]],{i,1,Length[vector[[1]]]},{j,Length[vector[[2]]]}];{})];

sameParRemove[vector_List,changeVector_List]:=Module[{tmp},FixedPoint[(tmp=sameNumber[#[[1]]];
If[tmp!={},Delete[#,{{1,1,tmp[[1]]},{1,2,tmp[[2]]},{2,tmp[[1]]},{2,Length[#[[1,1]]]+tmp[[2]]}}],#])&,{vector,changeVector}]];



(*--------------------------------------------------------------------------*)
(*Обработка исключительного случая неположительных параметров, а также случая только одной переменной и случая однаковых переменных*)

checkOfSpecialValues[vector_]:=Module[{},(
Or@@Map[(IntegerQ[#]&&#<=0)&,Flatten[Most[vector]]]
||Length[vector[[2]]]==0
||And@@Map[(IntegerQ[#])&,Flatten[Most[vector]]])];


simplifyPFQ[initialOper_,vector_?(checkOfSpecialValues[#]&),changeList_,initFactor_]:={explicitFormSimlifyPFQ[{initialOper,vector,1}],{},initFactor};

simplifyPFQ[initialOper_,vector_?(sameNumber[#]!={}&),changeList_,initFactor_]:=Module[{checkOfInitialData=listChecking[vector,changeList]},simplifyPFQ[initialOper,Sequence@@sameParRemove[vector,changeList],initFactor]];


(*--------------------------------------------------------------------------*)
(*Обработка общего случая неисключительных параметров. Случаи отриц параметров и одного обрабатываюся выше. Если при сдвижке получаются
исключительные праметры- то будет просто ошибка деления на ноль. За этим должен следить пользователь*)
simplifyPFQ[intialOper_,vector_,changeList_,initFactor_]:=Module[
{checkOfInitialData=listChecking[vector,changeList],
factor=initFactor,
initialVector={{{},vector,1,1,1}},
quantityOfTotalChangings,aiApBiDownChanges,listOfAdownAndBupch,numberAdBup
},
{
quantityOfTotalChangings=numberOfTotalChangings[changeList];
aiApBiDownChanges=listOfApandBdchanges[changeList,quantityOfTotalChangings ];
listOfAdownAndBupch=listOfAdownAndBupchanges[changeList-quantityOfTotalChangings];
numberAdBup=Ceiling[           Length[listOfAdownAndBupch]/Length[vector[[2]]]              ];


Apply[Dot,Prepend [Map[ (factor=factor*#[[3]];permutationMatrix[#[[1]],#[[2,3]],#[[4]]].DerivativeForPFQMatrix[    #[[2]] ]   )&   ,
Rest[
Nest[(Append[#,  simplifyPFQAiUpBiDownMatrix[#[[ -1,2 ]],aiApBiDownChanges[[Length[#]-numberAdBup]]]])&,


MapAt[ReplacePart[#,Join[#[[1,1]],Table[0,{i,1,Length[#[[2,1]]]-Length[#[[1,1]]]}]],{1,1}]&,
Map[ReplacePart[#,{Table[SymmetricPolynomial[-i,#[[1]]],{i,-Length[#[[1]]],0}]},1]&,
Nest[If[(numberOfAdBup=1+Sum[Length[#[[i,1]]],{i,1,Length[#]}]; 
           changevar=AdBupvector[#[[-1,2]],listOfAdownAndBupch[[numberOfAdBup]]];
           Mod[Length[#[[-1,1]]],Length[#[[-1,2,2]]]]==0
          ),Append[#, {{changevar[[1]]},changevar[[2]],1/changevar[[1]],0}],
ReplacePart[ReplacePart[Insert[#,changevar[[1]],{-1,1,1}],changevar[[2]], {-1,2}  ], #[[-1,3]]/ changevar[[1]],{-1,3}]]&  ,initialVector,Length[listOfAdownAndBupch]]]
,-1]


,Length[aiApBiDownChanges]]
]
]
,If[Length[intialOper]==Length[vector[[1]]], intialOper,            Table[If[i!=1,0,1],{i,1,Length[vector[[1]]]}]                         ]     ]
],
vectorAfterChange[vector,changeList],factor
}
];

(*--------------------------------------------------------------------------*)
explicitFormSimlifyPFQ[list_?(#[[2]]=={}&)]:=list[[1]] list[[3]];
explicitFormSimlifyPFQ[list_]:=Module[
{n=Length[list[[1]]],
funct=HypergeometricPFQ[list[[2,1]],list[[2,2]],list[[2,3]]]
},
(
explicitFormSimlifyPFQ::wrongParameter="For numerical checking the parameter `1`  in `2` should be a Symbol (it is used as a differentiating parameter) ";
If[ToString[Head[list[[2,3]]]]!="Symbol",{Message[explicitFormSimlifyPFQ::wrongParameter ,list[[2,3]],list[[2]]];Abort[]}];

 list[[3]] Apply[Plus,Join[{list[[1,1]]}funct,{Table[Sum[StirlingS2[i,j] list[[1,i+1]] list[[2,3]]^j,{i,i=j,n-1}],{j,1,n-1}].Table[D[funct,{list[[2,3]],i}],{i,1,n-1}]}]]
)];
(* явная форма для ToGroebnerBasis[]*)
explicitForm[List_]:=Plus@@Map[explicitFormSimlifyPFQ[#]&,List];
(*--------------------------------------------------------------------------*)
 
(* Считает коефициент при формуле, когда верхний параметр на целое чисто превышает нижний. Соотв. A(j1,j2,...jn)  в статье  *)
coeff[workingvector_,i_]:=Module[
{dualvector=Table[Unique[],{Length[workingvector[[1]]]-1}],
vector=Apply[Subtract,workingvector[[1]],1] ,
vectorOfJn ,tmp},(
vectorOfJn=Join[Nest[(Append[#,#[[-1]]+dualvector[[Length[#]]]])&,{i-Plus@@dualvector},Length[dualvector]],Table[i,{Evaluate[Length[workingvector[[2,2]] ]]  }]  ];

Sum[If[i-Plus@@dualvector>=0&&i-Plus@@dualvector<=vector[[1]], 
                                         Inner[Binomial,Sequence@@{vector,Prepend[dualvector,i-Plus@@dualvector]},Times]Inner[Pochhammer,Join[Rest[Extract[workingvector,Table[{1,k,1},{k,1,Length[workingvector[[1]]]}]]],workingvector[[2,1]]],vectorOfJn,Times]/Inner[Pochhammer,Join[Extract[workingvector,Table[{1,k,2},{k,1,Length[workingvector[[1]]]}]],workingvector[[2,2]]],vectorOfJn,Times]
   ,0],Evaluate[ If[dualvector=={},{tmp,0,0},Sequence@@MapIndexed[({dualvector[[#2[[1]]]],0,#1})&,Rest[vector]]]]
] 
)];

(*--------------------------------------------------------------------------*)
(*  Удаляет одинаковые параметры в векторе гипергеом функции   *)
sameParRemoveGroebnerBasis[vector_List]:=Module[{tmp},FixedPoint[(tmp=sameNumber[#];
If[tmp!={},Delete[#,{{1,tmp[[1]]},{2,tmp[[2]]}}],#])&,vector]];

(*--------------------------------------------------------------------------*)
(*Проверяет наличие случая, когда верхнийе параметры превышают нижние на целое число, выводят данные параметры и усеченный вектор гипергеом функции   *)
DiffIntPlusCheck[vector_List]:=Module[{},
(Table[If[IntegerQ[vector[[1,i]]-vector[[2,j]]]&&
                                  (vector[[1,i]]-vector[[2,j]]>0),Return[{i,j},Module]],{i,1,Length[vector[[1]]]},{j,Length[vector[[2]]]}];{})];

DiffIntRemove[vector_List]:=Module[{tmp},
FixedPoint[(tmp=DiffIntPlusCheck[#[[2]]];
If[tmp!={},(
                                            Delete[      Insert[#,{#[[ 2,1,tmp[[1]] ]], #[[2,2,tmp[[2]]   ]] },{1,1}]
,{{2,1,tmp[[1]]},          {2,2,tmp[[2]]}}]),#])&,{{},vector}]
];

(*--------------------------------------------------------------------------*)
(* ПРоверяет наличие  единички в верних параметрах, возвращает номер позиции единички checkForUnity[{{1,2+7 eps},{2+eps},x}]*)
checkForUnity[vector_?(#==={}&)]:=0;
checkForUnity[vector_]:=Module[{pos=Position[vector[[1]],1,1]},If[pos==={},0,pos[[1,1]]]];
(*Формула для сокращения верхней производной когда один из верхних параметров гипергеом функции равен 1 *)
UnityKilling[answer_?(checkForUnity[#[[2]]]==0&)]:={answer};
UnityKilling[answer_?(Length[#[[2,2]]]==1&)]:=Module[{unitPos=checkForUnity[answer[[2]]],ai,bi,per},{
 ai=Delete[answer[[2,1]],unitPos];bi=answer[[2,2]]-1; per=answer[[2,3]];
Print["случай 1   ",unitPos];
ReplacePart[Delete[answer,{{1,2}}],1->{answer[[1,1]]-answer[[1,2]] (-per ai[[1]]+bi[[1]])/(1-per)}],
{bi[[1]]/(1-per)answer[[1,2]],{},answer[[3]]}
}];

UnityKilling[answer_]:=Module[{unitPos=checkForUnity[answer[[2]]],ai,bi,per,p,i},{
 ai=Delete[answer[[2,1]],unitPos];
bi=answer[[2,2]]-1;
 per=answer[[2,3]];
p=Length[bi];

ReplacePart[answer,1->Delete[answer[[1]],{-1}]+Prepend[Table[SymmetricPolynomial[p-i,bi]-per SymmetricPolynomial[p-i,ai] ,{i,1,p-1}],Times@@bi-per Times@@ai]answer[[1,-1]]/(per-1)],
{Times@@bi answer[[1,-1]]/(1-per),{},answer[[3]]}

}];
(*--------------------------------------------------------------------------*)
(* Ищет целые числа в векторе гипергеом функции и выводит их в отдельный вектор. Надо для определния сдвижки к базису   *)
LookForInteger[x_,flag_]:=If[ToString[Head[x]]=="Plus",Plus@@Map[If[NumberQ[#],IntegerPart[#],0]&,List@@x], If[ToString[Head[x]]=="Rational"||ToString[Head[x]]=="Real"||ToString[Head[x]]=="Integer",If[flag==1,IntegerPart[x],Ceiling[x]],0]];

IntegerUnits[vector_]:=Join[Map[LookForInteger[#,1]&,Apart[vector[[1]]]],Map[LookForInteger[#,2]&,Apart[vector[[2]]]]];

ToGroebnerBasis[vector_?(checkOfSpecialValues[#]&)]:={{Print["special case"];explicitFormSimlifyPFQ[{{1},vector,1}],{},1}};

ToGroebnerBasis[vector_?(sameNumber[#]!={}&)]:=Module[{checkOfInitialData=listCheckingVector[vector]},ToGroebnerBasis[sameParRemoveGroebnerBasis[vector] ]  ];

(*--------------------------------------------------------------------------*)
ToGroebnerBasis[vector_]:=Module[
{checkvector=listCheckingVector[vector],
workingvector=DiffIntRemove[vector],
changeVector
},
(
(* определение базиса:*)
basis=Join[Table[0,{i,1,Length[ workingvector[[2,1]]    ]}],Table[1,{i,1,Length[workingvector[[2,1]]]-1}]  ];
If[Or@@Map[IntegerQ[#]&,workingvector[[2,1]]],basis=basis+1   ];
changeVector=-IntegerUnits[workingvector[[2]]]+basis;
Print["IntegerPart=",IntegerUnits[vector],"    changeVector=",changeVector,"    workingvector=",workingvector];

(* неисключительный случай:*)
(*If[workingvector[[1]]=={},simplifyPFQ[{1},workingvector[[2]],changeVector,1]]*)

(*Считает случай когда верхниеиндексы превышают нижний на целое число. В случае отсутствия таковых сводится к неисключительному случаю,
закомментированному выше*)
UnityKilling[
Join[{Plus@@Map[#[[1]]#[[3]]&,Join[Table[simplifyPFQ[PadRight[{workingvector[[2,3]]^i},Length[workingvector[[2,1]]],0],Insert[Most[workingvector[[2]]+i],workingvector[[2,3]],3],Table[-i,{Length[changeVector]}]+changeVector,coeff[workingvector,i]],{i,1,Subtract@@Plus@@workingvector[[1]]}],{simplifyPFQ[{1},workingvector[[2]],changeVector,1]}]]},If[workingvector[[2,2]]==={},{{},1},{workingvector[[2]]+Insert[Delete[Partition[Evaluate[Insert[changeVector,0,-1]],Length[workingvector[[2,1]]]],{-1,-1}],0,3],1}]]
                                      ]
)];
(*--------------------------------------------------------------------------*)


End[]
EndPackage[]

