4.1.1 Mathematica and Rubi grading function
(* Original version thanks to Albert Rich emailed on 03/21/2017 *)
(* ::Package:: *)
(* Nasser: April 7,2022. add second output which gives reason for the grade *)
(* Small rewrite of logic in main function to make it*)
(* match Maple's logic. No change in functionality otherwise*)
(* ::Subsection:: *)
(*GradeAntiderivative[result,optimal]*)
(* ::Text:: *)
(*If result and optimal are mathematical expressions, *)
(* GradeAntiderivative[result,optimal] returns*)
(* "F" if the result fails to integrate an expression that*)
(* is integrable*)
(* "C" if result involves higher level functions than necessary*)
(* "B" if result is more than twice the size of the optimal*)
(* antiderivative*)
(* "A" if result can be considered optimal*)
GradeAntiderivative[result_,optimal_] := Module[{expnResult,expnOptimal,leafCountResult,leafCountOptimal,finalresult},
expnResult = ExpnType[result];
expnOptimal = ExpnType[optimal];
leafCountResult = LeafCount[result];
leafCountOptimal = LeafCount[optimal];
(*Print["expnResult=",expnResult," expnOptimal=",expnOptimal];*)
If[expnResult<=expnOptimal,
If[Not[FreeQ[result,Complex]], (*result contains complex*)
If[Not[FreeQ[optimal,Complex]], (*optimal contains complex*)
If[leafCountResult<=2*leafCountOptimal,
finalresult={"A","none"}
,(*ELSE*)
finalresult={"B","Both result and optimal contain complex but leaf count is larger than twice the leaf count of optimal. $"<>ToString@leafCountResult<>"$ vs. $2("<>ToString@leafCountOptimal<>")="<>ToString[2*leafCountOptimal]<>"$."}
]
,(*ELSE*)
finalresult={"C","Result contains complex when optimal does not."}
]
,(*ELSE*)(*result does not contains complex*)
If[leafCountResult<=2*leafCountOptimal,
finalresult={"A","none"}
,(*ELSE*)
finalresult={"B","Leaf count is larger than twice the leaf count of optimal. $"<>ToString@leafCountResult<>"$ vs. $2("<>ToString@leafCountOptimal<>")="<>ToString[2*leafCountOptimal]<>"$."}
]
]
,(*ELSE*) (*expnResult>expnOptimal*)
If[FreeQ[result,Integrate] && FreeQ[result,Int],
finalresult={"C","Result contains higher order function than in optimal. Order "<>ToString[expnResult]<>" vs. order "<>ToString[expnOptimal]<>" in optimal."}
,
finalresult={"F","Contains unresolved integral."}
]
];
finalresult
]
(* ::Text:: *)
(*The following summarizes the type number assigned an *)
(*expression based on the functions it involves*)
(*1 = rational function*)
(*2 = algebraic function*)
(*3 = elementary function*)
(*4 = special function*)
(*5 = hyperpergeometric function*)
(*6 = appell function*)
(*7 = rootsum function*)
(*8 = integrate function*)
(*9 = unknown function*)
ExpnType[expn_] :=
If[AtomQ[expn],
1,
If[ListQ[expn],
Max[Map[ExpnType,expn]],
If[Head[expn]===Power,
If[IntegerQ[expn[[2]]],
ExpnType[expn[[1]]],
If[Head[expn[[2]]]===Rational,
If[IntegerQ[expn[[1]]] || Head[expn[[1]]]===Rational,
1,
Max[ExpnType[expn[[1]]],2]],
Max[ExpnType[expn[[1]]],ExpnType[expn[[2]]],3]]],
If[Head[expn]===Plus || Head[expn]===Times,
Max[ExpnType[First[expn]],ExpnType[Rest[expn]]],
If[ElementaryFunctionQ[Head[expn]],
Max[3,ExpnType[expn[[1]]]],
If[SpecialFunctionQ[Head[expn]],
Apply[Max,Append[Map[ExpnType,Apply[List,expn]],4]],
If[HypergeometricFunctionQ[Head[expn]],
Apply[Max,Append[Map[ExpnType,Apply[List,expn]],5]],
If[AppellFunctionQ[Head[expn]],
Apply[Max,Append[Map[ExpnType,Apply[List,expn]],6]],
If[Head[expn]===RootSum,
Apply[Max,Append[Map[ExpnType,Apply[List,expn]],7]],
If[Head[expn]===Integrate || Head[expn]===Int,
Apply[Max,Append[Map[ExpnType,Apply[List,expn]],8]],
9]]]]]]]]]]
ElementaryFunctionQ[func_] :=
MemberQ[{
Exp,Log,
Sin,Cos,Tan,Cot,Sec,Csc,
ArcSin,ArcCos,ArcTan,ArcCot,ArcSec,ArcCsc,
Sinh,Cosh,Tanh,Coth,Sech,Csch,
ArcSinh,ArcCosh,ArcTanh,ArcCoth,ArcSech,ArcCsch
},func]
SpecialFunctionQ[func_] :=
MemberQ[{
Erf, Erfc, Erfi,
FresnelS, FresnelC,
ExpIntegralE, ExpIntegralEi, LogIntegral,
SinIntegral, CosIntegral, SinhIntegral, CoshIntegral,
Gamma, LogGamma, PolyGamma,
Zeta, PolyLog, ProductLog,
EllipticF, EllipticE, EllipticPi
},func]
HypergeometricFunctionQ[func_] :=
MemberQ[{Hypergeometric1F1,Hypergeometric2F1,HypergeometricPFQ},func]
AppellFunctionQ[func_] :=
MemberQ[{AppellF1},func]