% Written by Min-Hua Huang
% June 9, 2011
% Test arithmetic rules of "h" in "Indefinite Integration of the Gamma Integral".
clear; clc; format long; s=2.4; a=1.2; b=-0.3;
fprintf('Testing the formulas of addition, subtraction, mutiplication, \nand division for the h function.\n');
% ************************************************************************
% addition
answer1=h(s,a+b,0); check1=0;
for i=0:1:20
    check1=check1+b^i/factorial(i)*difference(s,a,i);
end    
fprintf('\nchecking the addition rule: h(a+b)=%8.5f, formula=%8.5f,\n',answer1,check1);
% ************************************************************************
% subtraction
answer2=h(s,a-b,0); check2=0;
for i=0:1:20
    check2=check2+(-b)^i/factorial(i)*difference(s,a,i);
end    
fprintf('checking the subtraction rule: h(a-b)=%8.5f, formula=%8.5f,\n',answer2,check2);
% ************************************************************************
% multiplication
answer3=h(s,a*b,0); check3=h(s,a,0); f1=1;
for i=0:1:20
    f1=f1/(s+1+i);
    temp=(-a*b)^i*f1;
    check3=check3+a*(1-b)*temp*h(s+1+i,a,0);
end    
fprintf('checking the multiplication rule: h(a*b)=%8.5f, formula=%8.5f,\n',answer3,check3);
% ************************************************************************
% division
answer4=h(s,a/b,0); check4=h(s,a,0); f2=1;
for i=0:1:20
    f2=f2/(s+1+i);
    temp=(-a*b^-1)^i*f2;
    check4=check4+a*(1-b^-1)*temp*h(s+1+i,a,0);
end    
fprintf('checking the division rule: h(a/b)=%8.5f, formula=%8.5f,\n\n',answer4,check4);



