module definth;
fluid '(mellin!-transforms!* mellin!-coefficients!*);
algebraic <<
operator indefint2;
indefint2_rules :=
{ indefint2((~f1+~~f2)/~~f3,~x,~y) =>
indefint2(f1/f3,x,y)+indefint2(f2/f3,x,y) when not(f2=0),
indefint2(~n,~f1-~f2,~x,~y) =>
indefint2(n,f1,x,y)-indefint2(n,f2,x,y),
indefint2(~n,~f1+~f2,~x,~y) =>
indefint2(n,f1,x,y)+indefint2(n,f2,x,y),
indefint2(1/~x^(~~a),~f,~x,~y) => transf(defint_choose(f,x),-a,y,x),
indefint2(~x^(~~b)*sqrt(~x),~f,~x,~y) =>
transf(defint_choose(f,x),b+1/2,y,x),
indefint2(sqrt(~x),~f,~x,~y) =>
transf(defint_choose(f,x),1/2,y,x),
indefint2(~x^(~~a),~f,~x,~y) => transf(defint_choose(f,x),a,y,x),
indefint2(~b*~ff,~f,~x,~y) => b*indefint2(ff,f,x,y) when freeof (b,x),
indefint2(~b/(~~c*~ff),~f,~x,~y) => b/c*indefint2(1/ff,f,x,y)
when freeof (b,x) and freeof (c,x) and not(b=1 and c=1),
indefint2(~ff/~b,~f,~x,~y) =>
1/b*indefint2(ff,f,x,y) when freeof (b,x),
indefint2(~b*~ff,~f,~x,~y) => b*indefint2(ff,f,x,y) when freeof (b,x),
indefint2(~ff/~b,~f,~x,~y) =>
1/b*indefint2(ff,f,x,y) when freeof (b,x),
indefint2(~~b,~f,~x,~y) => b*indefint2(f,x,y)
when freeof (b,x) and not(b=1),
indefint2(~f,~x,~y) => transf(defint_choose(f,x),0,y,x)
};
let indefint2_rules;
symbolic procedure simpinteg(u);
begin scalar ff1,alpha,y,var,chosen_num,coef,!*uncached;
!*uncached := t;
ff1 := prepsq simp car u;
if ff1 = 'UNKNOWN then return simp 'UNKNOWN;
alpha := cadr u;
y := caddr u;
var := cadddr u;
chosen_num := cadr ff1;
if chosen_num = 0 then << alpha := caddr ff1;
return simp reval algebraic(alpha*y)>>
else
<< put('f1,'g,getv(mellin!-transforms!*,chosen_num));
coef := getv(mellin!-coefficients!*,chosen_num);
if coef then MELLINCOEF:= coef else MELLINCOEF :=1;
return simp list('new_mei,'f1 . cddr ff1,alpha,y,var)>>;
end$
put('new_mei,'simpfn,'new_meijer)$
symbolic procedure new_meijer(u);
begin scalar f,y,mellin,new_mellin,m,n,p,q,old_num,old_denom,temp,a1,
b1,a2,b2,alpha,num,denom,n1,temp1,temp2,coeff,v,var,new_var,new_y,
new_v,k;
f := prepsq simp car u;
y := caddr u;
mellin := bastab(car f,cddr f);
temp := car cddddr mellin;
var := cadr f;
if not idp VAR then RETURN error(99,'FAIL); % something is rotten, if not...
% better give up
temp := reval algebraic(sub(x=var,temp));
mellin := {car mellin,cadr mellin,caddr mellin,cadddr mellin,temp};
temp := reduce_var(cadr u,mellin,var);
alpha := simp!* car temp;
new_mellin := cdr temp;
if car cddddr new_mellin neq car cddddr mellin then
<< k := car cddddr mellin;
y := reval algebraic(sub(var=y,k));
new_y := simp y>>
else
<< new_var := car cddddr new_mellin;
new_y := simp reval algebraic(sub(x=y,new_var))>>;
n1 := addsq(alpha,'(1 . 1));
temp1 := {'expt,y,prepsq n1};
temp2 := cadddr new_mellin;
coeff := simp!* reval algebraic(temp1*temp2);
m := caar new_mellin;
n := cadar new_mellin;
p := caddar new_mellin;
q := car cdddar new_mellin;
old_num := cadr new_mellin;
old_denom := caddr new_mellin;
for i:=1 :n do
<< if old_num = nil then a1 := append(a1,{simp!* old_num })
else << a1 := append(a1,{simp!* car old_num});
old_num := cdr old_num>>;
>>;
for j:=1 :m do
<< if old_denom = nil then b1 := append(b1,{simp!* old_denom })
else << b1 := append(b1,{simp!* car old_denom});
old_denom := cdr old_denom>>;
>>;
a2 := listsq old_num;
b2 := listsq old_denom;
if a1 = nil and a2 = nil then
num := list({negsq alpha})
else if a2 = nil then num := list(append(a1,{negsq alpha}))
else
<< num := append(a1,{negsq alpha}); num := append({num},a2)>>;
if b1 = nil and b2 = nil then
denom := list({subtrsq(negsq alpha,'(1 . 1))})
else if b2 = nil then
denom := list(b1,subtrsq(negsq alpha,'(1 . 1)))
else
<< denom := list(b1,subtrsq(negsq alpha,'(1 . 1)));
denom := append(denom,b2)>>;
v := gfmsq(num,denom,new_y);
if v = 'fail then return simp 'fail
else v := prepsq v;
if eqcar(v,'meijerg) then new_v := v else new_v := simp v;
return multsq(new_v,coeff);
end$
symbolic procedure reduce_var(u,v,var1);
% Reduce Meijer G functions of powers of x to x
begin scalar var,m,n,coef,alpha,beta,alpha1,alpha2,expt_flag,k,temp1,
temp2,const,new_k;
var := car cddddr(v);
beta := 1;
% If the Meijer G-function is is a function of a variable which is not
% raised to a power then return initial function
if length var = 0 then return u . v
else
<< k := u; coef := cadddr v;
for each i in var do
<< if listp i then
<< if car i = 'expt then
<< alpha := caddr i; expt_flag := 't>>
else if car i = 'sqrt then
<< beta := 2; alpha := 1; expt_flag := 't>>
else if car i = 'times then
<< temp1 := cadr i; temp2 := caddr i;
if listp temp1 then
<< if car temp1 = 'sqrt then
<< beta := 2; alpha1 := 1; expt_flag := 't>>
else if car temp1 = 'expt and listp caddr temp1 then
<< beta := cadr cdaddr temp1;
alpha1 := car cdaddr temp1;
expt_flag := 't>>;
>>;
if listp temp2 and car temp2 = 'expt then
<< alpha2 := caddr temp2; expt_flag := 't>>;
if alpha1 neq 'nil then
alpha := reval algebraic(alpha1 + beta*alpha2)
else alpha := alpha2;
>>;
>>
else
<< if i = 'expt then << alpha := caddr var; expt_flag := 't>>;
>>;
>>;
% If the Meijer G-function is is a function of a variable which is not
% raised to a power then return initial function
if expt_flag = nil then return u . v
% Otherwise reduce the power by using the following formula :-
%
% y (c*y)^(m/n)
% / /
% | n |
% |t^k*F((c*t)^(m/n))dt = --------- |z^[((k + 1)*n - m)/m]*F(z)dz
% | m*c^(k+1) |
% / /
% 0 0
else
<< if listp alpha then << m := cadr alpha; n := caddr alpha;
n := reval algebraic(beta*n)>>
else << m := alpha; n := beta>>;
const := reval algebraic(sub(var1=1,var));
const := reval algebraic(1/(const^(n/m)));
new_k := reval algebraic(((k + 1)*n - m)/m);
coef := reval algebraic((n/m)*coef*(const)^(k+1));
var := reval algebraic(var^(n/m));
return {new_k,car v,cadr v, caddr v,coef,var}>>;
>>;
end$
>>;
endmodule;
end;