#
#  Integration file => dedicated code for program Dirac1
#  ----------------

#  Code for definite integration has been improved and generalised
#  for Maple 4.3 and even more so in Maple 4.4

# Set all integration counters to zero
counter :=0: counter2:=0: counter3:=0: counter4:=0:

coeffsterms := proc(p) local i,dum;
	i := map(proc(x) if has(x,t) then x fi end,indets(p));
	collect(p,i,distributed);
end;
pickout:=proc(f)
if type(f,`^`) then 
	op(2,f);
else 
	1;
fi;
end;
I1:=proc(uexp,ndeg,Eideg,mdeg,t)
local I1temp;
if mdeg=0 then
	if (ndeg=0) and (uexp<>0) then
		I1temp:=-((exp(t)^(-uexp))*Ei(Eideg*t)-
					Ei((Eideg-uexp)*t))/uexp;
	else
		I1temp:='I1(uexp,ndeg,Eideg,mdeg,t)';	
	fi;
else
	I1temp:='I1(uexp,ndeg,Eideg,mdeg,t)';
fi;
I1temp;
end:
`diff/I1` := proc(uexp,ndeg,Eideg,mdeg,t1,t2)
	if t1=t2 then
		(exp(t1)^(-uexp))*Ei(Eideg*t1)*(ln(t1)^mdeg)/(t1^ndeg);
	fi;
end:
prean:=proc(t,e,cof,u,mdeg,ndeg,Eideg)
local etemp,g,bad,i,o1,o2,Eitemp,mtemp,ntemp,utemp,h,r,r2,model;
case3:=0;
utemp:=0;
mtemp:=0;
ntemp:=0;
Eitemp:=0;
etemp:=e; 
bad:=1;
if type(etemp,`*`) then 
	g := 1; 
	for i in etemp do
	if has(i,t) then g := g*i else bad := bad*i fi
	od;
	if type(g,`*`) then
		h := 1;
		for o1 in g do
			if has(o1,exp(t)) then
				utemp:=-pickout(o1);
			else 
				h := h*o1;
			fi;
		od;
	else
		if has(g,exp) then 
			utemp:=-pickout(g);
		else 
			case3:=2;
		fi;
	fi;
	if (type(h,polynom)) or (type(1/h,polynom)) then ntemp:=degree(h,t) fi;
	if has(h,ln(t)) then
		if type(h,`*`) then
			r :=1;
			for o2 in h do
				if has(o2,ln(t)) then
					mtemp:=pickout(o2);
				else
					r := r*o2;
				fi;
			od;
		else
			mtemp:=pickout(h);
		fi;
	else
		r:=h;
	fi;
	if has(r,Ei) then
		if type(r,`*`) then
			r2 :=1;
			for o2 in r do
				if has(o2,Ei) then
					Eitemp:=op(1,o2)/t;
				else
					r2:=r2*o2;		
				fi;
			od;
		else
			Eitemp:=op(1,r)/t;
		fi;
	else
		r2:=r;
		Eitemp:=0;
	fi;
	if (type(r2,polynom)) or (type(1/r2,polynom)) then 
		ntemp:=degree(r2,t) 
	else
		case3:=2;
	fi;
else
	if has(etemp,exp) then
		utemp:=-pickout(etemp);
	else
		case3:=2;
	fi;	
fi;
model:=(exp(t)^(-utemp))*(ln(t)^mtemp)*(t^ntemp);
if Eitemp<>0 then
	model:=model*Ei(Eitemp*t);
fi;
if has(expand(etemp/model),t) then case3:=2 fi;
cof:=bad:
u:=utemp:
mdeg:=mtemp:
ndeg:=ntemp:
Eideg:=Eitemp:
end:
#
Parint:=proc(f,g,t)
local v,dv,res;
	dv:=expand(f/g);
	v:=IInt3(dv,t);
	res:=v*g - IInt3(v*diff(g,t),t);
res;
end:
IInt:=proc(f,t)
local res,temp,tempt,i;
	if type(f,`+`) then
		lprint();
		lprint(`Indefinite Integration`);
		lprint();
		temp := coeffsterms(f);
		if not type(temp,`+`) then RETURN(IInt(temp,t)) fi;
		res := 0;
		for i to nops(temp) do res := res + IInt(op(i,temp),t); od;
	else 
		tempt := map(proc(x,t) if has(x,t) then x else 1 fi end,f,t);
		res:=IInt2(tempt,t);
		counter4:=counter4+1;
		res := expand(f/tempt*res);
	fi;
	res;
end:
IInt2:=proc(f,t)
option remember;
counter3:=counter3+1;
IInt3(expand(f),t);
end:
#
IInt3:=proc(f,t)
local i,cof,uexp,mdeg,ndeg,Eideg,temp,temp2,temp3,res;
#
case4:=0;
temp:=expand(f):
temp:=map(normal,temp):
if type(temp,`+`) then
	res:=map(IInt3,temp,t);
elif has(temp,I1) then
	for i to nops(temp) do
		temp2:=op(i,temp);
		if has(temp2,I1) then 
			uexp:=op(1,temp2); ndeg:=op(2,temp2); 
			Eideg:=op(3,temp2); mdeg:=op(4,temp2);
		fi;
	od;
	res:=Parint(temp,I1(uexp,ndeg,Eideg,mdeg,t),t);
elif not(has(temp,Ei) or has(temp,ln(t))) then
	case4:=2;
else
	temp:=expand(temp);
	prean(t,temp,'cof','uexp','mdeg','ndeg','Eideg'):
	if (case3=0) then 
		if Eideg=0 then
			if (ndeg=0) and (mdeg=1) then
				res:=cof/uexp*(Ei(-uexp*t)*(exp(t)^uexp)-
						ln(t))*(exp(t)^(-uexp));
			elif uexp<>0 then
				if (ndeg=-1) and (mdeg=1) then
					res:=cof*(Ei(-uexp*t)*ln(t) - 
						I1(0,-ndeg,-uexp,0,t));
				elif (mdeg>0) then
					res:=Parint(temp,ln(t),t);
				else
					case4:=2;
				fi;
			else
				case4:=2;
			fi;
		else
			temp2:=expand(temp/Ei(Eideg*t));
			if has(temp2,t) then
				if (mdeg>0) then 
					if ndeg<>-1 then
						res:=Parint(temp,ln(t)^mdeg,t);
					else
						res:=cof*I1(uexp,-ndeg,Eideg,mdeg,t);
					fi;
				elif mdeg=0 then
					if ndeg>0 then
					res:=Parint(temp,Ei(Eideg*t),t);
					elif ndeg=-1 then
					res:=cof*I1(uexp,-ndeg,Eideg,mdeg,t);
					elif ndeg<0 then
					temp3:=expand(temp/(t^ndeg));
					res:=Parint(temp,temp3,t);
					else 
						if (Eideg-uexp)<>0 then
						res:=cof*I1(uexp,-ndeg,Eideg,mdeg,t);
					else
						res:=cof*(-1/uexp)*((exp(t)^(-uexp))*
							Ei(Eideg*t)-ln(t));
					fi;
					fi;
				else
					case4:=2;
				fi;
			else
				case4:=2;
			fi;
		fi;
	else
		case4:=2;
	fi;
fi;
if case4=2 then res:=int(temp,t); fi;
case4:=0;
expand(res);
end:
#
intfact:=proc(u,f,var) 
local ndeg,sm;
	ndeg:=degree(f,var);
	sm:=coeff(f,var,ndeg)*(ndeg!/(u^(ndeg+1)));
sm;
end:
ngritty:=proc(m,n)
local z;
if m>0 then
	subs(z=n+1,diff(GAMMA(z),z$m)); eval(");
else
	n!;
fi;
end:
intlog:=proc(f,x,m,n)
local o1,r,mdeg,ndeg;
case2:=1;
mdeg :=0;
ndeg :=0;
r :=1;
if type(f,`*`) then
	for o1 in f do
		if has(o1,ln(x)) then
			mdeg:=pickout(o1);
		else
			r := r*o1;
		fi;
	od;
	if type(r,polynom) then
		ndeg:=degree(r,x);
	else
		case2:=2;
	fi;
else
	mdeg:=pickout(f);
fi;
m:=mdeg;
n:=ndeg;
end:
#
IIntd:=proc(f,t)
local res,i,temp,tempt;
	if type(f,`+`) then
		lprint();
		lprint(`Definite Integration`);
		lprint();
		temp := coeffsterms(f);
		if not type(temp,`+`) then RETURN(IIntd(temp,t)) fi;
		res := 0;
		for i to nops(temp) do res := res+IIntd(op(i,temp),t); od;
	else 
		tempt := map(proc(x,t) if has(x,t) then x else 1 fi end,f,t);
		counter2:=counter2+1;
		res:=IIntd2(tempt,t);
		res:=expand(f/tempt*res);
	fi;
	res;
end:

IIntd2 := proc(e,t)
local etemp,h,i,u,m,n,o1,res;
option remember;
#
counter:=counter+1:
etemp:=e;
case:=1;
if type(etemp,`*`) then 
	h := 1;
	u := 0;
	for o1 in etemp  do
		if has(o1,exp(t)) then
			u := -pickout(o1);
		else 
			h := h*o1;
		fi;
	od;
	if u>0 then 
	if type(h,polynom) then
		res:=intfact(u,h,t);
		elif has(h,ln(t)) then
			intlog(h,t,'m','n');
			if case2=2 then
				res:=int(h*(exp(t)^(-u)),t=0..infinity);
			elif u=1 then
				res:=ngritty(m,n);
			else
			res:=0:
			for i from 0 to m do
				res:=res+(m!/(i!*(m-i)!))*((-ln(u))^(m-i))*
					ngritty(i,n):
			od:
				res:=res/(u^(n+1));
			fi;
		else
			case:=2;
		fi;
	else
		case:=2;
	fi;
else 
	if has(etemp,exp(t)) then
		u:=-pickout(etemp);
		res:=1/u;
	else
		case:=2;
	fi;
fi;
if case=2 then 
	res:=int(etemp,t=0..infinity) fi;
res := eval(res):
res := expand(res);
res;
end:
