Chapter 5
Appendix

5.1 Listing of grading functions

5.1 Listing of grading functions

The following are the current version of the grading functions used for grading the quality of the antiderivative with reference to the optimal antiderivative included in the test suite.

There is a version for Maple and for Mathematica/Rubi. There is a version for grading Sympy and version for use with Sagemath.

The following are links to the current source code.

  1. Mathematica and Rubi grading function GradeAntiderivative.m
  2. Maple grading function GradeAntiderivative.mpl
  3. Sympy grading function grade_sympy.py
  4. Sagemath grading function grade_sagemath.py

The following are the listings of source code of the grading functions.

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",""} 
                    ,(*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",""} 
            ,(*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]

Maple grading function

# File: GradeAntiderivative.mpl 
# Original version thanks to Albert Rich emailed on 03/21/2017 
 
#Nasser 03/22/2017  Use Maple leaf count instead since buildin 
#Nasser 03/23/2017  missing 'ln' for ElementaryFunctionQ added 
#Nasser 03/24/2017  corrected the check for complex result 
#Nasser 10/27/2017  check for leafsize and do not call ExpnType() 
#                   if leaf size is "too large". Set at 500,000 
#Nasser 12/22/2019  Added debug flag, added 'dilog' to special functions 
#                   see problem 156, file Apostol_Problems 
#Nasser 4/07/2022   add second output which gives reason for the grade 
 
GradeAntiderivative := proc(result,optimal) 
local leaf_count_result, 
        leaf_count_optimal, 
        ExpnType_result, 
        ExpnType_optimal, 
        debug:=false; 
 
        leaf_count_result:=leafcount(result); 
        #do NOT call ExpnType() if leaf size is too large. Recursion problem 
        if leaf_count_result > 500000 then 
            return "B","result has leaf size over 500,000. Avoiding possible recursion issues."; 
        fi; 
 
        leaf_count_optimal := leafcount(optimal); 
        ExpnType_result  := ExpnType(result); 
        ExpnType_optimal := ExpnType(optimal); 
 
        if debug then 
            print("ExpnType_result",ExpnType_result," ExpnType_optimal=",ExpnType_optimal); 
        fi; 
 
# 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 
 
    #This check below actually is not needed, since I only 
    #call this grading only for passed integrals. i.e. I check 
    #for "F" before calling this. But no harm of keeping it here. 
    #just in case. 
 
 
    if not type(result,freeof('int')) then 
        return "F","Result contains unresolved integral"; 
    fi; 
 
 
    if ExpnType_result<=ExpnType_optimal then 
        if debug then 
             print("ExpnType_result<=ExpnType_optimal"); 
        fi; 
        if is_contains_complex(result) then 
            if is_contains_complex(optimal) then 
                if debug then 
                        print("both result and optimal complex"); 
                fi; 
                if leaf_count_result<=2*leaf_count_optimal then 
                    return "A"," "; 
                else 
                   return "B",cat("Both result and optimal contain complex but leaf count of result is larger than twice the leaf count of optimal. ", 
                                   convert(leaf_count_result,string)," vs. $2 (", 
                                   convert(leaf_count_optimal,string)," ) = ",convert(2*leaf_count_optimal,string),"$."); 
                end if 
            else #result contains complex but optimal is not 
                if debug then 
                        print("result contains complex but optimal is not"); 
                fi; 
                return "C","Result contains complex when optimal does not."; 
            fi; 
        else # result do not contain complex 
             # this assumes optimal do not as well. No check is needed here. 
            if debug then 
                   print("result do not contain complex, this assumes optimal do not as well"); 
            fi; 
            if leaf_count_result<=2*leaf_count_optimal then 
                if debug then 
                    print("leaf_count_result<=2*leaf_count_optimal"); 
                fi; 
                return "A"," "; 
            else 
                if debug then 
                    print("leaf_count_result>2*leaf_count_optimal"); 
                fi; 
                return "B",cat("Leaf count of result is larger than twice the leaf count of optimal. $", 
                                   convert(leaf_count_result,string),"$ vs. $2(", 
                                   convert(leaf_count_optimal,string),")=",convert(2*leaf_count_optimal,string),"$."); 
            fi; 
        fi; 
    else #ExpnType(result) > ExpnType(optimal) 
        if debug then 
            print("ExpnType(result) > ExpnType(optimal)"); 
        fi; 
        return "C",cat("Result contains higher order function than in optimal. Order ", 
                       convert(ExpnType_result,string)," vs. order ", 
                       convert(ExpnType_optimal,string),"."); 
    fi; 
 
end proc: 
 
# 
# is_contains_complex(result) 
# takes expressions and returns true if it contains "I" else false 
# 
#Nasser 032417 
is_contains_complex:= proc(expression) 
  return (has(expression,I)); 
end proc: 
 
# 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 := proc(expn) 
  if type(expn,'atomic') then 
    1 
  elif type(expn,'list') then 
    apply(max,map(ExpnType,expn)) 
  elif type(expn,'sqrt') then 
    if type(op(1,expn),'rational') then 
       1 
    else 
       max(2,ExpnType(op(1,expn))) 
    end if 
  elif type(expn,'`^`') then 
    if type(op(2,expn),'integer') then 
      ExpnType(op(1,expn)) 
    elif type(op(2,expn),'rational') then 
      if type(op(1,expn),'rational') then 
         1 
      else 
         max(2,ExpnType(op(1,expn))) 
      end if 
    else 
         max(3,ExpnType(op(1,expn)),ExpnType(op(2,expn))) 
    end if 
  elif type(expn,'`+`') or type(expn,'`*`') then 
    max(ExpnType(op(1,expn)),max(ExpnType(rest(expn)))) 
  elif ElementaryFunctionQ(op(0,expn)) then 
    max(3,ExpnType(op(1,expn))) 
  elif SpecialFunctionQ(op(0,expn)) then 
    max(4,apply(max,map(ExpnType,[op(expn)]))) 
  elif HypergeometricFunctionQ(op(0,expn)) then 
    max(5,apply(max,map(ExpnType,[op(expn)]))) 
  elif AppellFunctionQ(op(0,expn)) then 
    max(6,apply(max,map(ExpnType,[op(expn)]))) 
  elif op(0,expn)='int' then 
    max(8,apply(max,map(ExpnType,[op(expn)]))) else 
  9 
  end if 
end proc: 
 
 
ElementaryFunctionQ := proc(func) 
  member(func,[ 
        exp,log,ln, 
        sin,cos,tan,cot,sec,csc, 
        arcsin,arccos,arctan,arccot,arcsec,arccsc, 
        sinh,cosh,tanh,coth,sech,csch, 
        arcsinh,arccosh,arctanh,arccoth,arcsech,arccsch]) 
end proc: 
 
SpecialFunctionQ := proc(func) 
  member(func,[ 
        erf,erfc,erfi, 
        FresnelS,FresnelC, 
        Ei,Ei,Li,Si,Ci,Shi,Chi, 
        GAMMA,lnGAMMA,Psi,Zeta,polylog,dilog,LambertW, 
        EllipticF,EllipticE,EllipticPi]) 
end proc: 
 
HypergeometricFunctionQ := proc(func) 
  member(func,[Hypergeometric1F1,hypergeom,HypergeometricPFQ]) 
end proc: 
 
AppellFunctionQ := proc(func) 
  member(func,[AppellF1]) 
end proc: 
 
# u is a sum or product.  rest(u) returns all but the 
# first term or factor of u. 
rest := proc(u) local v; 
  if nops(u)=2 then 
     op(2,u) 
  else 
     apply(op(0,u),op(2..nops(u),u)) 
  end if 
end proc: 
 
#leafcount(u) returns the number of nodes in u. 
#Nasser 3/23/17  Replaced by build-in leafCount from package in Maple 
leafcount := proc(u) 
   MmaTranslator[Mma][LeafCount](u); 
end proc:

Sympy grading function

#Dec 24, 2019. Nasser M. Abbasi: 
#              Port of original Maple grading function by 
#              Albert Rich to use with Sympy/Python 
#Dec 27, 2019  Nasser. Added `RootSum`. See problem 177, Timofeev file 
#              added 'exp_polar' 
from sympy import * 
 
def leaf_count(expr): 
    #sympy do not have leaf count function. This is approximation 
    return round(1.7*count_ops(expr)) 
 
def is_sqrt(expr): 
    if isinstance(expr,Pow): 
        if expr.args[1] == Rational(1,2): 
            return True 
        else: 
            return False 
    else: 
        return False 
 
def is_elementary_function(func): 
    return func in [exp,log,ln,sin,cos,tan,cot,sec,csc, 
            asin,acos,atan,acot,asec,acsc,sinh,cosh,tanh,coth,sech,csch, 
            asinh,acosh,atanh,acoth,asech,acsch 
        ] 
 
def is_special_function(func): 
    return func in [ erf,erfc,erfi, 
             fresnels,fresnelc,Ei,Ei,Li,Si,Ci,Shi,Chi, 
             gamma,loggamma,digamma,zeta,polylog,LambertW, 
             elliptic_f,elliptic_e,elliptic_pi,exp_polar 
         ] 
 
def is_hypergeometric_function(func): 
    return func in [hyper] 
 
def is_appell_function(func): 
    return func in [appellf1] 
 
def is_atom(expn): 
    try: 
        if expn.isAtom or isinstance(expn,int) or isinstance(expn,float): 
           return True 
        else: 
           return False 
 
    except AttributeError as error: 
        return False 
 
def expnType(expn): 
    debug=False 
    if debug: 
        print("expn=",expn,"type(expn)=",type(expn)) 
 
    if is_atom(expn): 
       return 1 
    elif isinstance(expn,list): 
        return max(map(expnType, expn))   #apply(max,map(ExpnType,expn)) 
    elif  is_sqrt(expn): 
        if isinstance(expn.args[0],Rational): #type(op(1,expn),'rational') 
            return 1 
        else: 
            return max(2,expnType(expn.args[0]))  #max(2,ExpnType(op(1,expn))) 
    elif isinstance(expn,Pow):   #type(expn,'`^`') 
        if isinstance(expn.args[1],Integer):  #type(op(2,expn),'integer') 
            return expnType(expn.args[0])   #ExpnType(op(1,expn)) 
        elif isinstance(expn.args[1],Rational):  #type(op(2,expn),'rational') 
            if isinstance(expn.args[0],Rational): #type(op(1,expn),'rational') 
                return 1 
            else: 
                return max(2,expnType(expn.args[0]))  #max(2,ExpnType(op(1,expn))) 
        else: 
            return max(3,expnType(expn.args[0]),expnType(expn.args[1])) #max(3,ExpnType(op(1,expn)),ExpnType(op(2,expn))) 
    elif isinstance(expn,Add) or isinstance(expn,Mul): #type(expn,'`+`') or type(expn,'`*`') 
        m1 = expnType(expn.args[0]) 
        m2 = expnType(list(expn.args[1:])) 
        return max(m1,m2)  #max(ExpnType(op(1,expn)),max(ExpnType(rest(expn)))) 
    elif is_elementary_function(expn.func):  #ElementaryFunctionQ(op(0,expn)) 
        return max(3,expnType(expn.args[0]))  #max(3,ExpnType(op(1,expn))) 
    elif is_special_function(expn.func): #SpecialFunctionQ(op(0,expn)) 
        m1 = max(map(expnType, list(expn.args))) 
        return max(4,m1)   #max(4,apply(max,map(ExpnType,[op(expn)]))) 
    elif is_hypergeometric_function(expn.func): #HypergeometricFunctionQ(op(0,expn)) 
        m1 = max(map(expnType, list(expn.args))) 
        return max(5,m1)   #max(5,apply(max,map(ExpnType,[op(expn)]))) 
    elif is_appell_function(expn.func): 
        m1 = max(map(expnType, list(expn.args))) 
        return max(6,m1)   #max(5,apply(max,map(ExpnType,[op(expn)]))) 
    elif isinstance(expn,RootSum): 
        m1 = max(map(expnType, list(expn.args))) #Apply[Max,Append[Map[ExpnType,Apply[List,expn]],7]], 
        return max(7,m1) 
    elif str(expn).find("Integral") != -1: 
        m1 = max(map(expnType, list(expn.args))) 
        return max(8,m1)   #max(5,apply(max,map(ExpnType,[op(expn)]))) 
    else: 
        return 9 
 
#main function 
def grade_antiderivative(result,optimal): 
 
    #print ("Enter grade_antiderivative for sagemath") 
    #print("Enter grade_antiderivative, result=",result," optimal=",optimal) 
 
    leaf_count_result  = leaf_count(result) 
    leaf_count_optimal = leaf_count(optimal) 
 
    #print("leaf_count_result=",leaf_count_result) 
    #print("leaf_count_optimal=",leaf_count_optimal) 
 
    expnType_result  = expnType(result) 
    expnType_optimal = expnType(optimal) 
 
    if str(result).find("Integral") != -1: 
        grade = "F" 
        grade_annotation ="" 
    else: 
        if expnType_result <= expnType_optimal: 
            if result.has(I): 
                if optimal.has(I): #both result and optimal complex 
                    if leaf_count_result <= 2*leaf_count_optimal: 
                        grade = "A" 
                        grade_annotation ="" 
                    else: 
                        grade = "B" 
                        grade_annotation ="Both result and optimal contain complex but leaf count of result is larger than twice the leaf count of optimal. "+str(leaf_count_result)+" vs. $2 ("+str(leaf_count_optimal)+") = "+ str(2*leaf_count_optimal)+"$." 
                else: #result contains complex but optimal is not 
                    grade = "C" 
                    grade_annotation ="Result contains complex when optimal does not." 
            else: # result do not contain complex, this assumes optimal do not as well 
                if leaf_count_result <= 2*leaf_count_optimal: 
                    grade = "A" 
                    grade_annotation ="" 
                else: 
                    grade = "B" 
                    grade_annotation ="Leaf count of result is larger than twice the leaf count of optimal. "+str(leaf_count_result)+" vs. $2 ("+str(leaf_count_optimal)+") = "+ str(2*leaf_count_optimal)+"$." 
        else: 
            grade = "C" 
            grade_annotation ="Result contains higher order function than in optimal. Order "+str(ExpnType_result)+" vs. order "+str(ExpnType_optimal)+"." 
 
 
    #print("Before returning. grade=",grade, " grade_annotation=",grade_annotation) 
 
    return grade, grade_annotation

SageMath grading function

#Dec 24, 2019. Nasser: Ported original Maple grading function by 
#              Albert Rich to use with Sagemath. This is used to 
#              grade Fricas, Giac and Maxima results. 
#Dec 24, 2019. Nasser: Added 'exp_integral_e' and 'sng', 'sin_integral' 
#              'arctan2','floor','abs','log_integral' 
#June 4, 2022  Made default grade_annotation "none" instead of "" due 
#              issue later when reading the file. 
#July 14, 2022. Added ellipticF. This is until they fix sagemath, then remove it. 
 
from sage.all import * 
from sage.symbolic.operators import add_vararg, mul_vararg 
 
debug=False; 
 
def tree_size(expr): 
    r""" 
    Return the tree size of this expression. 
    """ 
    #print("Enter tree_size, expr is ",expr) 
 
    if expr not in SR: 
        # deal with lists, tuples, vectors 
        return 1 + sum(tree_size(a) for a in expr) 
    expr = SR(expr) 
    x, aa = expr.operator(), expr.operands() 
    if x is None: 
        return 1 
    else: 
        return 1 + sum(tree_size(a) for a in aa) 
 
def is_sqrt(expr): 
    if expr.operator() == operator.pow:   #isinstance(expr,Pow): 
        if expr.operands()[1]==1/2: #expr.args[1] == Rational(1,2): 
            if debug: print ("expr is sqrt") 
            return True 
        else: 
            return False 
    else: 
        return False 
 
def is_elementary_function(func): 
    #debug=False 
    m = func.name() in ['exp','log','ln', 
            'sin','cos','tan','cot','sec','csc', 
            'arcsin','arccos','arctan','arccot','arcsec','arccsc', 
            'sinh','cosh','tanh','coth','sech','csch', 
            'arcsinh','arccosh','arctanh','arccoth','arcsech','arccsch','sgn', 
        'arctan2','floor','abs' 
        ] 
    if debug: 
        if m: 
            print ("func ", func , " is elementary_function") 
        else: 
            print ("func ", func , " is NOT elementary_function") 
 
 
    return m 
 
def is_special_function(func): 
    #debug=False 
    if debug: 
        print ("type(func)=", type(func)) 
 
    m= func.name() in ['erf','erfc','erfi','fresnel_sin','fresnel_cos','Ei', 
           'Ei','Li','Si','sin_integral','Ci','cos_integral','Shi','sinh_integral' 
           'Chi','cosh_integral','gamma','log_gamma','psi','zeta', 
           'polylog','lambert_w','elliptic_f','elliptic_e','ellipticF', 
           'elliptic_pi','exp_integral_e','log_integral', 
           'weierstrassPInverse','weierstrass','weierstrassP','weierstrassZeta', 
           'weierstrassPPrime','weierstrassSigma'] 
 
    if debug: 
        print ("m=",m) 
        if m: 
            print ("func ", func ," is special_function") 
        else: 
            print ("func ", func ," is NOT special_function") 
 
 
    return m 
 
 
def is_hypergeometric_function(func): 
    return func.name() in ['hypergeometric','hypergeometric_M','hypergeometric_U'] 
 
def is_appell_function(func): 
    return func.name() in ['hypergeometric']   #[appellf1] can't find this in sagemath 
 
def is_atom(expn): 
 
    #debug=False 
    if debug: 
         print ("Enter is_atom, expn=",expn) 
 
    if  not hasattr(expn, 'parent'): 
        return False 
 
 
    #thanks to answer at https://ask.sagemath.org/question/49179/what-is-sagemath-equivalent-to-atomic-type-in-maple/ 
    try: 
        if expn.parent() is SR: 
            return expn.operator() is None 
        if expn.parent() in (ZZ, QQ, AA, QQbar): 
            return expn in expn.parent() # Should always return True 
        if hasattr(expn.parent(),"base_ring") and hasattr(expn.parent(),"gens"): 
            return expn in expn.parent().base_ring() or expn in expn.parent().gens() 
 
        return False 
 
    except AttributeError as error: 
        print("Exception,AttributeError in  is_atom") 
        print ("cought exception" , type(error).__name__ ) 
        return False 
 
 
def expnType(expn): 
 
    if debug: 
        print (">>>>>Enter expnType, expn=", expn) 
        print (">>>>>is_atom(expn)=", is_atom(expn)) 
 
    if is_atom(expn): 
        return 1 
    elif type(expn)==list:   #isinstance(expn,list): 
        return max(map(expnType, expn))   #apply(max,map(ExpnType,expn)) 
    elif  is_sqrt(expn): 
        if  type(expn.operands()[0])==Rational: #type(isinstance(expn.args[0],Rational): 
            return 1 
        else: 
            return max(2,expnType(expn.operands()[0]))  #max(2,expnType(expn.args[0])) 
    elif expn.operator() == operator.pow:   #isinstance(expn,Pow) 
        if type(expn.operands()[1])==Integer:  #isinstance(expn.args[1],Integer) 
            return expnType(expn.operands()[0])   #expnType(expn.args[0]) 
        elif type(expn.operands()[1])==Rational:  #isinstance(expn.args[1],Rational) 
            if type(expn.operands()[0])==Rational: #isinstance(expn.args[0],Rational) 
                return 1 
            else: 
                return max(2,expnType(expn.operands()[0]))  #max(2,expnType(expn.args[0])) 
        else: 
            return max(3,expnType(expn.operands()[0]),expnType(expn.operands()[1])) #max(3,expnType(expn.operands()[0]),expnType(expn.operands()[1])) 
    elif expn.operator() == add_vararg or expn.operator() == mul_vararg: #isinstance(expn,Add) or isinstance(expn,Mul) 
        m1 = expnType(expn.operands()[0]) #expnType(expn.args[0]) 
        m2 = expnType(expn.operands()[1:]) #expnType(list(expn.args[1:])) 
        return max(m1,m2)  #max(ExpnType(op(1,expn)),max(ExpnType(rest(expn)))) 
    elif is_elementary_function(expn.operator()):  #is_elementary_function(expn.func) 
        return max(3,expnType(expn.operands()[0])) 
    elif is_special_function(expn.operator()): #is_special_function(expn.func) 
        m1 = max(map(expnType, expn.operands()))      #max(map(expnType, list(expn.args))) 
        return max(4,m1)   #max(4,m1) 
    elif is_hypergeometric_function(expn.operator()): #is_hypergeometric_function(expn.func) 
        m1 = max(map(expnType, expn.operands()))       #max(map(expnType, list(expn.args))) 
        return max(5,m1)   #max(5,m1) 
    elif is_appell_function(expn.operator()): 
        m1 = max(map(expnType, expn.operands()))       #max(map(expnType, list(expn.args))) 
        return max(6,m1)   #max(6,m1) 
    elif str(expn).find("Integral") != -1: #this will never happen, since it 
                #is checked before calling the grading function that is passed. 
                #but kept it here. 
        m1 = max(map(expnType, expn.operands()))       #max(map(expnType, list(expn.args))) 
        return max(8,m1)   #max(5,apply(max,map(ExpnType,[op(expn)]))) 
    else: 
        return 9 
 
#main function 
def grade_antiderivative(result,optimal): 
 
 
    if debug: 
        print ("Enter grade_antiderivative for sagemath") 
        print("Enter grade_antiderivative, result=",result) 
        print("Enter grade_antiderivative, optimal=",optimal) 
        print("type(anti)=",type(result)) 
        print("type(optimal)=",type(optimal)) 
 
    leaf_count_result  = tree_size(result) #leaf_count(result) 
    leaf_count_optimal = tree_size(optimal) #leaf_count(optimal) 
 
    #if debug: print ("leaf_count_result=", leaf_count_result, "leaf_count_optimal=",leaf_count_optimal) 
 
 
    expnType_result  = expnType(result) 
    expnType_optimal = expnType(optimal) 
 
    if debug: print ("expnType_result=", expnType_result, "expnType_optimal=",expnType_optimal) 
 
    if expnType_result <= expnType_optimal: 
        if result.has(I): 
            if optimal.has(I): #both result and optimal complex 
                if leaf_count_result <= 2*leaf_count_optimal: 
                    grade = "A" 
                    grade_annotation ="none" 
                else: 
                    grade = "B" 
                    grade_annotation ="Both result and optimal contain complex but leaf count of result is larger than twice the leaf count of optimal. "+str(leaf_count_result)+" vs. $2 ("+str(leaf_count_optimal)+") = "+ str(2*leaf_count_optimal)+"$." 
            else: #result contains complex but optimal is not 
                grade = "C" 
                grade_annotation ="Result contains complex when optimal does not." 
        else: # result do not contain complex, this assumes optimal do not as well 
            if leaf_count_result <= 2*leaf_count_optimal: 
                grade = "A" 
                grade_annotation ="none" 
            else: 
                grade = "B" 
                grade_annotation ="Leaf count of result is larger than twice the leaf count of optimal. "+str(leaf_count_result)+" vs. $2 ("+str(leaf_count_optimal)+") = "+ str(2*leaf_count_optimal)+"$." 
    else: 
        grade = "C" 
        grade_annotation ="Result contains higher order function than in optimal. Order "+str(expnType_result)+" vs. order "+str(expnType_optimal)+"." 
 
 
    print("Before returning. grade=",grade, " grade_annotation=",grade_annotation) 
 
    return grade, grade_annotation