#=============================================================================
# Some useful Maple procedures.
#=============================================================================
# Marc A. Murison
# Astronomical Applications Dept.
# U.S. Naval Observatory
# 3450 Massachusetts Ave., NW
# Washignton, DC 20392
# murison@aa.usno.navy.mil
# http://aa.usno.navy.mil/murison/
#=============================================================================


utils := module()

	 description "a package of handy utilities";
	 option package;
	 export tmsg, debug_print, sniffmem, disremember, sniffmax, freemem,
			remove_subscripts, restore_subscripts, Subs, invsubs, invalgsubs,
			matsubs, location, Collect, `Collect/locsubs`, `Collect/algsubs`, simp,
			srsimp, exp_to_10, topsqrts, topsqrt, signsqrt, getsqrts, remove_sqrts, 
			restore_sqrts, rootfunc, patfunc, trigsimp, cosfix, sinfix, coshfix, sinhfix,
			CombineTrig, FunctionCalls, `CombineTrig/optimize`, CombTrigOpt, normalize,
			annular_average, is_real_parts, realsols, sortsols, combined_cost,
			sort_by_cost, rotaxis, rotvec, funcops, csquare, applyfunc, termfunc,
			pullout, mat2eqs, factormat, polytest, small_divisors, common_factor,
			expansion, `simplify/fractional_power`, eqcollect, combine_angle, utilsinit;
	local	is_fractional_power;
			
utilsinit := proc()
	 global `type/fractional_power`;
	`type/fractional_power` := eval(is_fractional_power);
end proc;


#=============================================================================
# Define a "fractional_power" type.  An expression is of this type if
# it is raised to a fractional power.
#=============================================================================
is_fractional_power := proc( expr::{algebraic,equation,list,set,indexed,array} )
    local p;
    if type(expr,`^`) then
        if type(op(2,expr),fraction) then
            return true;
        end if;
    end if;
    false;
end proc;


#=============================================================================
# Print messages of the form
#
# <proc_name> [<elapsed time>]: <msg>
#
# to the terminal screen, where elapsed time is the current
# time minus the time stored in the optional third argument.
#=============================================================================
tmsg := proc( proc_name, msg, t0 )
    local dt;
    if nargs=3 then
        dt := time() - t0;
    else
        dt := time();
    end if;
    if proc_name <> `` then
        if dt < 1 then
            printf(""||(convert(proc_name,string))||"[0]: %s\n",msg);
        else
            printf(""||(convert(proc_name,string))||"[%.0f]: %s\n",dt,msg);
        end if;
    else
        if dt < 1 then
            printf(`[0] %s\n`,msg);
        else
            printf(`[%.0f] %s\n`,dt,msg);
        end if;
    end if;
end proc;



#=============================================================================
# pname = procname in calling procedure
#=============================================================================
debug_print := proc( pname, msg, verbosity_level )

    global verbosity, time0, doprint;

    doprint := false;
    if verbosity_level=0 then
        doprint := true;
    elif assigned( printlevel ) then
        if verbosity_level <= printlevel then
            doprint := true;
        end if;
    end if;
    
    if doprint then
        if assigned( time0 ) then
            tmsg( pname, msg, time0 );
        else
            tmsg( pname, msg );
        end if;
        if nargs > 3 then
            print( eval(args[4..nargs]) );
        else
            print();                    #force a flush
        end if;
    end if;

end proc;



#=============================================================================
# A pair of memory utilities posted by Rober Israel (and modified slightly
# to get rid of lexical scoping problem with N)
#=============================================================================
sniffmem := proc()
    local L, N;
    if nargs > 0 then 
        N := args[1];
    else 
        N := 1000;
    end if;
    L := select( (t,n) -> type(t,procedure) and n < length(op(4,eval(t))), [anames()], N );
    L := map(t -> [t,length(op(4,eval(t)))], L);
    sort(L, (a, b) -> evalb(a[2] < b[2]));
end proc;

# this one will cause the 'option remember' to be removed from a procedure
disremember := proc(p::procedure) 
    local L;
    L := {op(3,eval(p))} minus {remember};
    subsop(3=op(L),eval(p));
end proc;

#=============================================================================
# max memory pig as returned by sniffmem()
#=============================================================================
sniffmax := proc()
    local L, Lmax, k;
    L := sniffmem(0);   #returns [[proc1,n1],[proc2,n2],...]
    Lmax := L[1];
    for k from 2 to nops(L) do
        if L[k][2] > Lmax[2] then
            Lmax := L[k];
        end if;
    od;
    Lmax;
end proc;
  
#=============================================================================
# free the memory associated with the remember tables of routines 
# returned from sniffmem()
#=============================================================================
freemem := proc()
    local L, N, tmp;
    if nargs > 0 then 
        N := args[1];
    else 
        N := 1000;
    end if;
    if not assigned(forget) then
        readlib(forget);
    end if;
    L := sniffmem(N);   #returns [[proc1,n1],[proc2,n2],...]
    while L<>[] do
        L := map((x)->x[1],L);
        tmp := traperror(forget(op(L)));
        if tmp=lasterror then
            debug_print(procname,"Error in forget():",1,lasterror);
            gc();
            RETURN;
        end if;
        L := sniffmem(N);   #returns [[proc1,n1],[proc2,n2],...]
    od;
    gc();
end proc;
  


#=============================================================================
# Remove subscripts, e.g. v[x] -> vx.  Use restore_subscripts() to restore.
#=============================================================================
remove_subscripts := proc( expr )

    local   p1, p2;
    global  _subscripts;
    
    if type(expr,{list,set,matrix,equation,function,`+`,`*`,`^`}) then
        return  map(procname,args);
    end if;

    if not assigned(_subscripts) then
        _subscripts := {};
    end if;

    if type(expr,diff) then
        if type( op(1,expr), function ) then
            p1 := op([1,0],expr);
        else
            p1 := op(1,expr);
        end if;
        p2 := op(2,expr);
        _subscripts := _subscripts union {``||p1||p2=expr};
        ``||p1||p2;
    elif type(expr,indexed) then
        p1 := op(0,expr);
        p2 := op(1,expr);
        _subscripts := _subscripts union {``||p1||p2=expr};
        ``||p1||p2;
    else
        expr;
    end if;

end proc;



#=============================================================================
# Restore subscripts, e.g. vx -> v[x], after using remove_subscripts().
#=============================================================================
restore_subscripts := proc( expr )
    global  _subscripts;
    subs( _subscripts, expr );
    _subscripts := '_subscripts';
    %%;
end proc;



#=============================================================================
# A better subs()
# Optional arguments:
# - a variable or list of variables or subexpressions can be protected from
#   the effects of subs() by specifying protect=[...]
# - a function to be applied to the subs() result. Optional additional 
#   function arguments must follow the function name
#=============================================================================
Subs := proc( expr::anything, subslist::{`=`,list(`=`),set(`=`),list(list)} )

    local Q, slist, p, protectlist, argslist, _P, k, listseq;

    #handle a list of lists for the substitutions list
    if type(subslist,list(list)) then
        slist := [];
        for p in subslist do
            slist := [op(slist),p];
        od;
    else
        slist := [copy(subslist)];
    end if;
    
    #look for a list of variables to protect
    protectlist := [];
    argslist    := [args[3..nargs]];
    for p in args[3..nargs] do
        if type(p,`=`) then
            if lhs(p)=`protect` then
                if type(rhs(p),list) then
                    listseq := op(rhs(p));
                else
                    listseq := rhs(p);
                end if;
                protectlist := [op(protectlist),listseq];
                argslist    := remove(has,argslist,p);
                break;
            end if;
        end if;
    od;
    
    #rename the protected variables
    Q := copy(expr);
    if nops(protectlist) > 0 then
        for k from 1 to nops(protectlist) do
            Q := subs( protectlist[k]=_P[k], Q );
        od;
    end if;

    #make the desired 'normal' substitution 
    Q := subs( op(slist), Q );

    #apply optional user-specified function to the substitution result
    if nops(argslist) > 0 then
        Q := argslist[1](Q,op(argslist[2..nops(argslist)]));
        if nops(protectlist) > 0 then
            Q := invsubs( protectlist, Q );
        end if;
    end if;

    #restore the protected variables
    if nops(protectlist) > 0 then
        for k from 1 to nops(protectlist) do
            Q := subs( _P[k]=protectlist[k], Q );
        od;
    end if;
    
    eval(Q);
end proc;



#=============================================================================
# Do an "inverse" subs().  That is, for each x=y in subslist, do a
# substitution in expr of y=x.
#=============================================================================
invsubs := proc( subslist::{`=`,list(`=`),set(`=`)}, expr::anything )

    local oldsubs, newsubs, p;

    oldsubs := copy(subslist);
    if type(oldsubs,`=`) then
        oldsubs := [oldsubs];
    end if;

    newsubs := [];
    for p in oldsubs do
        newsubs := [op(newsubs),rhs(p)=lhs(p)];
    od;

    if type(oldsubs,set) then
        newsubs := convert(newsubs,set);
    end if;

    subs( newsubs, expr );
end proc;



#=============================================================================
# Do an "inverse" algsubs().  That is, for each x=y in subslist, do an
# "algebraic" substitution in expr of y=x.  An extension of the
# algsubs() syntax is that for invalgsubs() you can specify more than
# one substitution, as in subs() and invsubs().
#=============================================================================
invalgsubs := proc( subslist::{`=`,list(`=`),set(`=`)}, expr::anything )

    local oldsubs, newsubs, p, newexpr;

    oldsubs := copy(subslist);
    if type(oldsubs,`=`) then
        oldsubs := [oldsubs];
    end if;

    newsubs := [];
    for p in oldsubs do
        newsubs := [op(newsubs),rhs(p)=lhs(p)];
    od;

    newexpr := copy(expr);
    for p in newsubs do
        newexpr := algsubs( p, newexpr );
    od;

    newexpr;
end proc;



#=============================================================================
# Make column vector substitutions.  Example:
#
# > conicoid := z-(x^2+y^2)/(2*f+sqrt(4*f^2 - epsilon*(x^2+y^2)));
#
#                                         2    2
#                                        x  + y
#     conicoid := z - ---------------------------------------------
#                               2            2            2 1/2
#                     2 f + (4 f  - epsilon x  - epsilon y )
#
# > mat(x,y,z) = mat(x0,y0,z0) + t*mat(vx,vy,vz);
#
#                               [x]   [x0]     [vx]
#                               [y] = [y0] + t [vy]
#                               [z]   [z0]     [vz]
#
# > rootfunc( matsubs( conicoid, % ), collect, epsilon );
#
#                                       2              2
#                            (x0 + t vx)  + (y0 + t vy)
#   z0 + t vz - --------------------------------------------------------
#                                   2              2               2 1/2
#               2 f + ((-(x0 + t vx)  - (y0 + t vy) ) epsilon + 4 f )
#=============================================================================
matsubs := proc( expr, subsvec::`=` )

    local i, k, tmp, s, n;

    s := evalm(subsvec);

    if not type(lhs(s),'array'(1..integer,1..1))
    or not type(rhs(s),'array'(1..integer,1..1)) then
        ERROR(`: 2nd arg must be of the form array(1..n,1..1)=array(1..n,1..1)`);
    end if;

    n := rows(lhs(s));
    for i from 1 to n do
        tmp[i] := lhs(s)[i,1] = rhs(s)[i,1];
    od;

    subs( seq(tmp[k],k=1..n), expr );
end proc;



#=============================================================================
# Find the op() location of term in expr.  location returns a Maple list
# suitable for use in op().  Hence, if it is possible to return
# term via the op() function, then
#    op( location(expr,term), expr ) = term
#=============================================================================
location := proc( expr::anything, term::anything )

    local L, M, k, parts, i;

    if has(expr,term) then

        if type(expr,table) then
            L := [3];
            parts := [op(op(3,expr))];
        else
            L := [];
            parts := [op(expr)];
        end if;

        #check each whole term of this level of op(expr)
        if member( term, parts, k ) then
            return [k];
        else
            #descend to each term of this level from left to right
            for i from 1 to nops(parts) do
                M := procname( parts[i], term );
                if nops(M) > 0 then
                    return [op(L),i,op(M)];
                end if;
            od;
        end if;
    end if;

    return [];
end proc;



#=============================================================================
# A better collect().  Optional argument is 'loc' or 'alg'.  'loc' forces
# use of location() method, while 'alg' forces use of algsubs() method.  The
# default is 'alg'.
#
# Side effects:  global _T_ cleared.
#=============================================================================
Collect := proc( expr::{algebraic,list,set,relation,array},
                         subexpr::{algebraic,list,set} )

    local   Z, s, k, j, p, L, form, fargs, skip;
    global  _T_;

    if not type( eval(expr), algebraic ) then
        return map(procname,args);
    end if;

    if not type(subexpr,{list,set}) then
        s := [subexpr];
    else
        s := subexpr;
    end if;

    #
    # extract collection method and adjust function args
    #
    form := 'alg';
    for p in args[3..nargs] do
        if p='loc' then
            form := 'loc';
        end if;
    od;
    fargs := op( remove( has, [args[3..nargs]], {'alg','loc'} ) );

    #
    # check to see if we can use Maple's collect()
    #
    select( hastype, s, {`+`,`*`} );
    if % = [] then
        skip := false;
        for p in s do
            if type(p,fractional_power) then
                skip := true;
                break;
            end if;
        od;
        if not skip then
            return collect(expr,subexpr,fargs);
        end if;
    end if;

    Z   := expr;
    _T_ := '_T_';

    #
    # dummy substitutions
    #
    for k from 1 to nops(s) do
        if form='loc' then
            Z := `Collect/locsubs`( Z, s, k );
        else
            Z := `Collect/algsubs`( Z, s, k );
        end if;
    od;

    #
    # collect on the dummies
    #
    Z := collect( Z, [seq(_T_[j],j=1..nops(s))], fargs );

    #
    # subs back
    #
    for k from 1 to nops(s) do
        Z := subs( _T_[k]=s[k], Z );
    od;

    Z;

end proc;

`Collect/locsubs` := proc( Z, s, k )
    local   L, tmp;
    global  _T_;
    #
    # repeat subsop until no more subexpressions in Z
    #
    tmp := Z;
    L   := location( tmp, s[k] );
    while nops(L) > 0 do
        tmp := subsop( L=_T_[k], tmp );
        L   := location( tmp, s[k] );
    od;
    tmp;
end proc;

`Collect/algsubs` := proc( Z, s, k )
    local   Q, p, H, M, S, j;
    global  _T_;

    Q := Z;

    #
    # if any of the subexpressions is the argument of a diff,
    # then we have to substitute a dummy or else a bug
    # in algsubs will delete the derivative from the expression!
    #
    if hastype(Q,diff) then
        M := MakeCoeffs(Q);
        S := [];
        j := 0;
        for p in M[derivs] do
            if isdiff(p) and has(p,s[k]) then
                j := j + 1;
                S := [op(S),p];     #save the substituted diffs
                Q := subs(p=H[j],Q);
            end if;
        od;
    end if;

    #
    # now do the algebraic substitutions
    #
    Q := traperror( algsubs( s[k]=_T_[k], Q, exact ) );

    #
    # check if algsubs failed
    #
    if member(`cannot compute degree of pattern in`,{Q})  #algsubs failed
    or member(`no variables in pattern`,{Q}) then
        debug_print(procname,
                    `subexpression too complicated for algsubs:`,
                    1, lasterror, s[k] );
        debug_print(procname,`trying location() approach...`,1);
        return `Collect/locsubs`( Z, s, k );
    end if;

    #
    # subs back for the diffs
    #
    if nops(S) <> 0 then
        for j from 1 to nops(S) do
            Q := subs( H[j]=S[j], Q );
        od;
    end if;

    Q;
end proc;



#------------------------------------------------------------------------------
# A simplification command that tries harder.
#
# optional arguments:
#	- expand    perform expand(expr) before starting the simplification process
#   - symbolic  perform simplify(expr,symbolic)
#	- assume=x  where 
#               x = none
#               x = values valid in simplify(expr,assume=x)
#
# Note: default is assume=real.
#------------------------------------------------------------------------------
simp := proc( expr )
    local k, sexpr, tmp, assumespec, do_expand;
    if type(expr,{equation,list,set,array}) then
    	return map(procname,args);
	end if;
    assumespec := assume=real;
    do_expand  := false;
    for k from 2 to nargs do
        if type(args[k],`=`) and lhs(args[k])=`assume` then
            if rhs(args[k])=`none` then
                assumespec := NULL;
            else
                assumespec := args[k];
            end if;
        elif args[k]=`symbolic` then
        	assumespec := symbolic;
        elif args[k]=`expand` then
            do_expand := true;
        end if;
    od;
    debug_print(procname,"initial cost:",4,combined_cost(expr),expr);
    #knock down sqrt expressions first
    debug_print(procname,"initial rootfunc...",4);
    sexpr := rootfunc(expr,simplify,assumespec);
    if do_expand then
        debug_print(procname,"expand...",4,sexpr);
    	sexpr := expand(sexpr);
    end if;
    debug_print(procname,"cost:",4,combined_cost(sexpr));
    debug_print(procname,"simplify...",4,sexpr);
    tmp := simplify(sexpr,assumespec);
    if has( tmp, signum ) then
        tmp := convert(tmp,abs);
    end if;
    if combined_cost(tmp) < combined_cost(sexpr) then
        debug_print(procname,"cost is smaller",4);
        sexpr := tmp;
    end if;
    debug_print(procname,"cost:",4,combined_cost(sexpr));
    debug_print(procname,"factor...",4,sexpr);
    tmp := factor(sexpr);
    if combined_cost(tmp) < combined_cost(sexpr) then
        debug_print(procname,"cost is smaller",4);
        sexpr := tmp;
    end if;
    debug_print(procname,"cost:",4,combined_cost(sexpr));
    debug_print(procname,"combine+factor...",4,sexpr);
    tmp := factor(combine(sexpr));
    if combined_cost(tmp) < combined_cost(sexpr) then
        debug_print(procname,"cost is smaller",4);
        sexpr := tmp;
        tmp := simplify(sexpr,assumespec);
        debug_print(procname,"simplify after combine+factor...",4,sexpr);
		if combined_cost(tmp) < combined_cost(sexpr) then
            debug_print(procname,"cost is smaller",4);
			sexpr := tmp;
		end if;
    end if;
    debug_print(procname,"cost:",4,combined_cost(sexpr));

    debug_print(procname,"collect...",4,sexpr);
    tmp := collect(sexpr,[op(indets(sexpr,function))],factor);
    if combined_cost(tmp) < combined_cost(sexpr) then
        debug_print(procname,"cost is smaller",4);
        sexpr := tmp;
    end if;
    debug_print(procname,"2nd rootfunc...",4,sexpr);
    tmp := rootfunc(sexpr,factor);
    if combined_cost(tmp) < combined_cost(sexpr) then
        debug_print(procname,"cost is smaller",4);
        sexpr := tmp;
    end if;
    if combined_cost(expr) < combined_cost(sexpr) then
        debug_print(procname,"ORIGINAL cost is smaller",4);
        sexpr := expr;
    end if;
    debug_print(procname,"final cost:",4,combined_cost(sexpr),sexpr);
    sexpr;
end proc;



#=============================================================================
# Side relations simplification with
#   - temporary replacement of subscripts
#   - collection on the side relations
#   - optional function and function args to apply during collection
#=============================================================================
srsimp := proc( expr::{`=`,algebraic,array,list,set},
                        siderelslist::{`=`,algebraic,list,set} )

    local srlist, sublist, A, k;

    if not type(expr,algebraic) then
        return map(procname,args);
    end if;

    if not type(siderelslist,{list,set}) then
        srlist := [siderelslist];
    elif type(siderelslist,set) then
        srlist := convert(siderelslist,list);
    else
        srlist := copy(siderelslist);
    end if;
    srlist := remove_subscripts(srlist);

    sublist := [seq(srlist[k]=A||k,k=1..nops(srlist))];
    remove_subscripts(expr);
    simplify(%,convert(sublist,set));
    collect( %, [seq(A||k,k=1..nops(srlist))], args[3..nargs] );
    invsubs( sublist, % );
    restore_subscripts(%);
end proc;



#=============================================================================
# convert powers of e to powers of 10 in an expression
#=============================================================================
exp_to_10 := proc( expr )
	indets(expr,function);
	select(patmatch,%,exp(a::anything));
	convert(%,list);
	map( (x)->utils:-simp(expand(x)), % );
	subs( zip( (x,y)->x=y, %%, % ), expr );
end proc;



#=============================================================================
# Extracts the top-level fractional powers and
# returns them in a list.  For example,
#
# > solve(a+b*t+c*t^2,t);
# 
#                         2                          2
#              -b + sqrt(b  - 4 a c)      -b - sqrt(b  - 4 a c)
#          1/2 ---------------------, 1/2 ---------------------
#                        c                          c
# 
# > (sqrt(%[2]*A+B)+C)^2;
# 
#           /                       2                       \2
#           |           (-b - sqrt(b  - 4 a c)) A           |
#           |1/2 sqrt(2 ------------------------- + 4 B) + C|
#           \                       c                       /
# 
# > topsqrts(%);
# 
#                          2
#              (-b - sqrt(b  - 4 a c)) A               2
#      [sqrt(2 ------------------------- + 4 B), sqrt(b  - 4 a c)]
#                          c
# 
#=============================================================================
topsqrts := proc( expr::{algebraic,equation,list,set} )
    local x, p, rootlist;
    if type(expr,equation) then
        rootlist := [op(procname(lhs(expr))),op(procname(rhs(expr)))];
    elif not hastype(expr,fractional_power) then
#       debug_print(procname,`No sqrt terms in expression!`,5);
        return [];
    else
        rootlist := [];
        for p in expr do
            if type(p,fractional_power) and not type(op(1,p),numeric) then
                rootlist := [op(rootlist),p];
            elif type(p,{list,set,algebraic,equation}) then 
                rootlist := [op(rootlist),op(procname(p))];
            end if;
        od;
    end if;
    convert(rootlist,set);
    convert(%,list);
end proc;



#=============================================================================
# Extract the first fractional power term encountered in an expression.
#=============================================================================
topsqrt := proc( expr::{algebraic,`=`} )
    local p;
    if not hastype(expr,fractional_power) then
        debug_print(procname,`No sqrt term in expression!`,4,expr);
        return FAIL;
    end if;
    if type( expr, fractional_power ) then
        expr;
    else
        # check top level terms
        for p in expr do
            if type( p, fractional_power ) then
                return p;
            end if;
        od;
        # descend subexpressions
        for p in expr do
            if hastype(p,fractional_power) then
                return procname(p);
            end if;
        od;
    end if;
    return FAIL;
end proc;



#=============================================================================
# Return the sign of a term involving a fractional power.
# Use getsqrts() to isolate fractional power terms.  For example,
#
# > f := solve( a + b*t + c*t^2, t );
# > g := (sqrt(f[2]*A+B)+C)^2;
# > h := topsqrts(g);
#
#                   2         1/2
#               A (b  - 4 c a)
#   h := [- 1/2 -----------------,
#                      c
#
#       /               2         1/2      \1/2
#       |   A b     A (b  - 4 c a)         |
#       |-2 --- - 2 ----------------- + 4 B|    C]
#       \    c             c               /
#
# > signsqrt(h[1]);
#
#                                   -1
#=============================================================================
signsqrt := proc( expr::{`*`,`^`} )
    local p, cursign;
    if not hastype(expr,fractional_power) then
        debug_print(procname,`No sqrt terms in expression!`,4);
        return FAIL;
    end if;
    cursign := 1;
    for p in expr do
        cursign := cursign*sign(p);
    od;
    return cursign;
end proc;



#=============================================================================
# Extracts all fractional power terms in expr and returns them in a
# list structure.  The list structure is of the form [[s1,loc1],[s2,loc2],...]
# where s1, etc. are the terms and loc1, etc. are the expression parse tree
# locations as defined by the op function.  Hence, s1=op(loc1,expr) and so
# on.  For example,
#
# > f := solve( a + b*t + c*t^2, t );
# > g := (sqrt(f[2]*A+B)+C)^2;
#
#          /    /          2         1/2         \1/2    \2
#          |    |  (-b - (b  - 4 c a)   ) A      |       |
#      g = |1/2 |2 ------------------------ + 4 B|    + C|
#          \    \            c                   /       /
#
# > getsqrts(g);
#=============================================================================
getsqrts := proc( expr::{algebraic,equation,list,set,indexed,array} )

    local k, p, oplist1, oplist2;

    if nargs > 1 then
        oplist1 := args[2];
        oplist2 := args[3];
    else
        oplist1 := [];
        oplist2 := [];
    end if;

    if type(expr,fractional_power) then
        if oplist1 <> [] then
            oplist2 := [ op(oplist2), oplist1 ];
        end if;
        if hastype( op(1,expr), fractional_power ) then
            oplist2 := [ op(oplist2), procname(op(1,expr),oplist2,[]) ];
        end if;
        return oplist2;
    elif hastype(expr,fractional_power) then
        for k from 1 to nops(expr) do
            p := op(k,expr);
            if type(p,fractional_power) then
                oplist2 := [ op(oplist2), [op(oplist1),k] ];
                if hastype( op(1,p), fractional_power ) then
                    oplist2 := [ op(oplist2), procname(op(1,p),oplist2,[op(oplist1),k]) ];
                end if;
            elif hastype(p,fractional_power) then
                if oplist2 <> [] then
                    oplist2 := [ op(oplist2), procname(p,oplist2,[op(oplist1),k]) ];
                else
                    oplist2 := procname(p,oplist2,[op(oplist1),k]);
                end if;
            end if;
        od;
    end if;

    return oplist2;
end proc;



#=============================================================================
# Remove all fractional power terms from an expression, replacing them
# with dummy variables A1, A2, etc.
#
# Global variables used: _subslist, _subsvar
#=============================================================================
remove_sqrts := proc( expr::{algebraic,equation,list,set}, A::name )
    local k, sqrtlist;
    global _subslist, _subsvars;
    sqrtlist := topsqrts(expr);
    if nops(sqrtlist)=0 then
        return expr;
    end if;
    _subslist := [];
    _subsvars := [];
    for k from 1 to nops(sqrtlist) do
        _subslist := [op(_subslist),sqrtlist[k]=A||k];
        _subsvars := [op(_subsvars),A||k];
    od;
    subs( _subslist, expr );
end proc; 



#=============================================================================
# Replace dummy variables with corresponding fractional power expressions,
# as stored in _subslist by remove_sqrts().
#
# Global variables used: _subslist, _subsvar
#=============================================================================
restore_sqrts := proc( expr::{algebraic,equation,list,set} )
    local rexpr;
    global _subslist, _subsvars;
    if not assigned(_subslist) or not assigned(_subsvars) then
        return expr;
    end if;
    rexpr := invsubs( _subslist, expr );
    _subslist := '_subslist';
    _subsvars := '_subsvars';
    rexpr;
end proc;



#=============================================================================
# Simplify fractional power expressions by collecting on the fractional
# powers and factoring the coefficients. 
#
# NOTE: Despite the name, Maple's simplify() procedure, which is
#       expensive and which generally uglifies things, is NOT used.
#
# Procedures used:  remove_sqrts(), restore_sqrts()
# Global variables used: _subslist, _subsvar
#=============================================================================
`simplify/fractional_power` := proc( expr::{algebraic,list,set,`=`} )
    local k, Z;
    if type(expr,{list,set,`=`}) then
        return map(procname,args);
    elif hastype(expr,fractional_power) then
        restore_sqrts(
            collect( remove_sqrts(expr,Z), _subsvars, factor )
        );
    else
        expr;
    end if;
end proc;



#=============================================================================
# Recursively apply a procedure, func(), to subexpressions that are raised
# to a fractional power.  This is handy for collecting terms under square
# roots.  For example,
#
#   rootfunc( expr, collect, [var1, var2], factor )
#
# Any arguments beyond the second argument, expr, are treated as additional
# arguments of func().  In the above example, subexpressions are passed
# to collect as collect( expr, [var1, var2], factor ).  rootfunc() is able
# to handle multiply-nested fractional power expressions.
#=============================================================================
rootfunc := proc( expr::{algebraic,equation,list,set,indexed,array,`=`},
                          func::procedure )
    local p, z, loc;

    if not type(expr,algebraic) then
        return map(procname,args);
    elif type(expr,fractional_power) then
        return applyop(func,1,expr,args[3..nargs]);
    else
        z := expr;
        for p in expr do
            if (nops(p)>1 or type(p,function)) and hastype(p,fractional_power) then
                loc := location(z,p);
                if type(p,fractional_power) then
                    p := applyop(func,1,p,args[3..nargs]);
                else
                    p := procname(p,func,args[3..nargs]);
                end if;
                z := subsop(loc=p,z);
            end if;
        od;
        z;
    
    end if;

end proc;



#=============================================================================
# Look for expanded sin/cos(a+b) forms and combine them back to sin/cos(a+b)
#=============================================================================
combine_angle := proc( expr::{algebraic,`=`,list,set,'array'(algebraic)} )
	option remember, system;
	local T;
	if not hastype(expr,trig) then
		expr;
	elif not type(expr,algebraic) then
		map(procname,args);
	elif (nops(expr)=1 and not type(expr,function)) or type(expr,numeric) then
		expr;
	elif type(expr,function) then
	 	op(0,expr)( seq( procname(op(k,expr)), k=1..nops(expr) ) );
	else
		define( T,
			T(c::algebraic
				+ f::algebraic*sin(a::algebraic)*cos(b::algebraic)
				+ f::algebraic*cos(a::algebraic)*sin(b::algebraic))
				= T(c) + f*sin(a+b),
			T(c::algebraic
				- f::algebraic*sin(a::algebraic)*cos(b::algebraic)
				+ f::algebraic*cos(a::algebraic)*sin(b::algebraic))
				= T(c) - f*sin(a-b),
			T(c::algebraic
				+ f::algebraic*sin(a::algebraic)*cos(b::algebraic)
				- f::algebraic*cos(a::algebraic)*sin(b::algebraic))
				= T(c) + f*sin(a-b),
			T(c::algebraic
				- f::algebraic*sin(a::algebraic)*cos(b::algebraic)
				- f::algebraic*cos(a::algebraic)*sin(b::algebraic))
				= T(c) - f*sin(a+b),
			T(c::algebraic
				+ f::algebraic*cos(a::algebraic)*cos(b::algebraic)
				+ f::algebraic*sin(a::algebraic)*sin(b::algebraic))
				= T(c) + f*cos(a-b),
			T(c::algebraic
				- f::algebraic*cos(a::algebraic)*cos(b::algebraic)
				+ f::algebraic*sin(a::algebraic)*sin(b::algebraic))
				= T(c) - f*cos(a+b),
			T(c::algebraic
				+ f::algebraic*cos(a::algebraic)*cos(b::algebraic)
				- f::algebraic*sin(a::algebraic)*sin(b::algebraic))
				= T(c) + f*cos(a+b),
			T(c::algebraic
				- f::algebraic*cos(a::algebraic)*cos(b::algebraic)
				- f::algebraic*sin(a::algebraic)*sin(b::algebraic))
				= T(c) - f*cos(a-b),

			T(c::algebraic
				+ g::numeric*f::algebraic*sin(a::algebraic)*cos(b::algebraic)
				+ h::numeric*f::algebraic*cos(a::algebraic)*sin(b::algebraic))
				= T(c) + f*g*T( sin(a)*cos(b)+h/g*cos(a)*sin(b) ),
			T(c::algebraic
				+ g::numeric*f::algebraic*cos(a::algebraic)*cos(b::algebraic)
				+ h::numeric*f::algebraic*sin(a::algebraic)*sin(b::algebraic))
				= T(c) + f*g*T( cos(a)*cos(b)+h/g*sin(a)*sin(b) ),

			T(a::nonunit(algebraic)+b::nonunit(algebraic)) = T(a)+T(b),
			T(a::nonunit(algebraic)*b::nonunit(algebraic)) = T(a)*T(b),
			T(a::nonunit(algebraic)^b::nonunit(algebraic)) = T(a)^T(b),
			T(a::anything)=a
		);
		return T(expr);
	end if;
end proc:



#=============================================================================
# Apply a function to a subexpression matching a pattern.
# Optional additional arguments are assumed to be additional arguments
# of func.  For example,
#   patfunc( expr, pattern, collect, [terms], factor );
#=============================================================================
patfunc := proc( expr::{algebraic,`=`,list,set,'array'(algebraic)},
                         pattern::algebraic, func::procedure )
    local p, tmp;
    if not type(expr,algebraic) then
        return map( procname, args );
    end if;
    if (nops(expr)=1 and not type(expr,function)) or type(expr,numeric) then
        expr;
    elif type(expr,function) then
    	op(0,expr)( seq( procname(op(k,expr),args[2..nargs]), k=1..nops(expr) ) );
    elif patmatch(expr,pattern) then
        func(expr,args[4..nargs]);
    else
        tmp := copy(expr);
        for p in expr do
            if patmatch(p,pattern) then
                tmp := subsop( location(tmp,p)=func(p,args[4..nargs]), tmp );
            else
                tmp := subsop( location(tmp,p)=procname(p,args[2..nargs]), tmp );
            end if;
        od;
        eval(tmp);
    end if;
end proc;



#=============================================================================
# Simplify: sin^2+cos^2 = 1
#=============================================================================
trigsimp := proc( expr )
    patfunc( expr, sin(x::algebraic)^2+cos(y::algebraic)^2, simplify );
    patfunc( expr, -sin(x::algebraic)^2-cos(y::algebraic)^2, simplify );
end proc;



#=============================================================================
# Force 1 - cos^2 ==> sin^2.
#=============================================================================
cosfix := proc( expr )
    local F1, F2, F3, F4;
    F1 := proc(y)               #y will be 1-cos(x)^2
		op(1,op(indets(y,trig)));
        sin(%)^2;
    end;
    F2 := proc(y)               #y will be -1+cos(x)^2
		op(1,op(indets(y,trig)));
        -sin(%)^2;
    end;
    F3 := proc(y)               #y will be a+1-cos(x)^2
   		select(patmatch,{op(y)},a::numeric*cos(b::algebraic)^2);
   		op(op(indets(op(1,%),trig)));	#x
        y-1+cos(%)^2+sin(%)^2;
    end;
    F4 := proc(y)               #y will be a-1+cos(x)^2
    	select(patmatch,{op(y)},a::numeric*cos(b::algebraic)^2);
   		op(op(indets(op(1,%),trig)));	#x
		y+1-cos(%)^2-sin(%)^2;
    end;
    patfunc( expr, 1-cos(x::algebraic)^2, F1 );
    patfunc( %, cos(x::algebraic)^2-1, F2 );
    patfunc( %, a::nonunit(algebraic)+1-cos(x::algebraic)^2, F3 );
    patfunc( %, a::nonunit(algebraic)-1+cos(x::algebraic)^2, F4 );
end proc;



#=============================================================================
# Force 1 - sin^2 ==> cos^2.
#=============================================================================
sinfix := proc( expr )
    local F1, F2, F3, F4;
    F1 := proc(y)               #y will be 1-sin(x)^2
		op(1,op(indets(y,trig)));
        cos(%)^2;
    end;
    F2 := proc(y)               #y will be -1+sin(x)^2
		op(1,op(indets(y,trig)));
        -cos(%)^2;
    end;
    F3 := proc(y)               #y will be a+1-sin(x)^2
    	select(patmatch,{op(y)},a::numeric*sin(b::algebraic)^2);
   		op(op(indets(op(1,%),trig)));	#x
        y-1+sin(%)^2+cos(%)^2;
    end;
    F4 := proc(y)               #y will be a-1+sin(x)^2
    	select(patmatch,{op(y)},a::numeric*sin(b::algebraic)^2);
   		op(op(indets(op(1,%),trig)));	#x
        y+1-sin(%)^2-cos(%)^2;
    end;
    patfunc( expr, 1-sin(x::algebraic)^2, F1 );
    patfunc( %, sin(x::algebraic)^2-1, F2 );
    patfunc( %, a::nonunit({name,numeric,`+`,`*`,`^`})+1-sin(x::algebraic)^2, F3 );
    patfunc( %, a::nonunit({name,numeric,`+`,`*`,`^`})-1+sin(x::algebraic)^2, F4 );
end proc;



#=============================================================================
# Force cosh^2 - 1 ==> sinh^2.
#=============================================================================
coshfix := proc( expr )
    local F;
    F := proc(y)                #y will be cosh(x::algebraic)^2-1
        select( has, y, cosh ); #cosh(x)^2
        op(1,%);                #cosh(x)
        op(%);                  #x
        sinh(%)^2;
    end;
    patfunc( expr, cosh(x::algebraic)^2-1, F );
end proc;



#=============================================================================
# Force 1 + sinh^2 ==> cosh^2.
#=============================================================================
sinhfix := proc( expr )
    local F;
    F := proc(y)                #y will be 1+sinh(x::algebraic)^2
        select( has, y, sinh ); #sinh(x)^2
        op(1,%);                #sinh(x)
        op(%);                  #x
        cosh(%)^2;
    end;
    patfunc( expr, 1+sinh(x::algebraic)^2, F );
end proc;



#=============================================================================
# Given an expression composed of various trigonometric terms involving
# at least two angles theta[1] and theta[2], combine the terms containing
# sines and cosines of theta[1] and theta[2] into single sine and cosine
# terms.  Other sines and cosines involving other angles are protected
# from the Maple combine() procedure.  The input expression may also be a
# one or two dimensional array.
#=============================================================================
CombineTrig := proc( expr::{algebraic,`=`,array}, var1::algebraic, var2::algebraic )

    local p, k, j, sublist, A, kmax, f, f1, f2, g, isfun;

    f := copy(expr);
    if type(f,{'array'(1),'array'(2),`=`}) then
        return map( procname, f, var1, var2 );
    end if;
    if not ( has(var1,indets(f)) and has(var2,indets(f)) ) then
        return eval(f);
    end if;

    # f is a function, such as arcsin(stuff)
    if type(f,function) then
        g := f;
        f := op(f);         #do the simplification on "stuff"
        isfun := YES;
    else
        isfun := NO;
    end if;

#   f := simplify( expand(f) );
    if not hasfun(f,sin) and not hasfun(f,cos) then
        if isfun=YES then
            return eval(g);
        else
            return eval(f);
        end if;
    end if;

    #
    # Substitute for sines and cosines of other indeterminates
    # in order to hide them from combine().  Must be careful, since
    # if var1 or var2 are themselves functions of other variables
    # [e.g., lambda(t)] then Maple will classify them as functions;
    # hence the "if hasfun..." condition.
    #
    sublist := array;
    k := 1;
    for p in indets(f,function) do
        if not p=sin(var1) and not p=cos(var1) and
           not p=sin(var2) and not p=cos(var2) and
           not p=var1 and not p=var2 then
            if hasfun(p,cos) or hasfun(p,sin) then
                sublist[k] := p;
                k := k + 1;
            end if;
        end if;
    od;
    kmax := k-1;
    f := subs( seq( sublist[k]=A[k], k=1..kmax ), f );

    #
    # Simplify the remaining trig terms with combine()
    #
    f := combine(f);

    #
    # Restore previously substituted sines and cosines
    #
    f := subs( seq(A[k]=sublist[k], k=1..kmax ), f );

    if isfun=YES then
        f := op(0,g)(f);
    end if;

    f;

end proc;



#=============================================================================
# Returns the number of function calls required by expr.
#=============================================================================
FunctionCalls := proc( expr )
    local t, p;
    t := cost(expr);
    if not has(t,`functions`) then
        return 0;
    end if;
    if type( t, `+` ) then
        for p in t do
            if has( p, `functions` ) then
                if is( op(1,p), integer ) then  #e.g. p = 5 `functions`
                    return op(1,p);
                else                            #p = `functions`
                    return 1;
                end if;
            end if;
        od;
    elif type( t, `*` ) then
        return op(1,t);
    elif t=`functions` then
        return 1;
    else
        return 0;
    end if;
end proc;



#=============================================================================
# aliased to CombTrigOpt
#
# Use CombineTrig() and FunctionCalls() to find the combination of
# angles, two at a time, that produces the minimum number of function
# calls and return that most-efficient form of expr.
#
# optional arguments:
# verbose       If arg2 or arg3 is ON, then progress is printed.
# [a1,a2...]    List of angles to perform optimization on.  If omitted,
#               then optimization is performed on all trig function arguments.
#=============================================================================
`CombineTrig/optimize` := proc( expr )

    local p, j, k, kmax, vars, L1, L2, c0, c1, c2, f, f1, f2,
          n, verbose;

    if nargs = 2 and args[2]=ON then
        verbose := ON;
    elif nargs > 2 and (args[2]=ON or args[3]=ON) then
        verbose := ON;
    else
        verbose := OFF;
    end if;

    f := copy(expr);
    if type(f,'array'(1)) then
        return map( procname, f, args[2..nargs] );
    end if;
    if type(f,'array'(2)) then
        return map( procname, f, args[2..nargs] );
    end if;
    if type(f,`=`) then
        return map( procname, f, args[2..nargs] );
    end if;
    if type(f,{`*`,`+`}) and has(f,`array`) then        #A*[array] or A+[array]
        return map( procname, f, args[2..nargs] );
    end if;

    #
    # fiducial cost
    #
    c0 := FunctionCalls(expand(f));
    if verbose=ON then
        printf(`fiducial cost = %d functions`, c0 );
        print();    #force the printf to the screen
    end if;
    if c0=1 then
        return eval(f);
    end if;

    #
    # pick out the angles (or at least the indeterminates)
    #
    if nargs=1 or (nargs=2 and not type(args[2],list)) or
       (not type(args[2],list) and not type(args[3],list)) then
        k := 1;
        vars := [];
        for p in indets(f) do
            if not hasfun(p,cos) and not hasfun(p,sin) and not isdiff(p) then
                vars := [op(vars),p];
                k := k + 1;
            else        # include already-existing multiple angles
                if type( op(1,p), `+` ) then
                    vars := [op(vars),op(1,p)];
                    k := k + 1;
                end if;
            end if;
        od;
    elif type(args[2],list) then
        vars := copy(args[2]);
    elif nargs=3 and type(args[3],list) then
        vars := copy(args[3]);
    else
        ERROR(` wrong arg type(s)`);
    end if;
    vars := convert(vars,set);  # weed out any repeats
    kmax := nops(vars);
    if kmax < 2 then
        return simplify(f);
    end if;
    if verbose=ON then
        printf(`Trying %d combinations...`, kmax*(kmax-1)/2 );
        print();
    end if;

    #
    # loop over angle combinations
    #
    L1 := CombineTrig( f, vars[1], vars[2] );
    c1 := FunctionCalls( L1 );
    if c1=0 then
        return eval(L1);
    end if;
    if verbose=ON then
        printf( `  1: cost[%a][%a] = %d functions`, vars[1], vars[2], c1 );
        lprint(` `);
    end if;
    n := 2;
    for j from 1 to kmax do
        for k from j+1 to kmax do
            if k <> j and not (j=1 and k=2) then
                L2 := CombineTrig( f, vars[j], vars[k] );
                c2 := FunctionCalls( L2 );
                if verbose=ON then
                    printf( `%3d: cost[%a][%a] = %d functions`,
                             n, vars[j], vars[k], c2 );
                    print();
                    n := n + 1;
                end if;
                if c2 < c1 then
                    L1 := L2;
                    c1 := c2;
                end if;
            end if;
        od;
    od;

    if c0 <= c1 then
        if verbose=ON then
            printf( `Best form is original: cost = %d functions\n`, c0 );
        end if;
        return eval(f);  #no more efficient form than originl
    else
        if verbose=ON then
            printf( `Best cost = %d functions\n`, c1 );
        end if;
        return eval(L1);
    end if;

end proc;

CombTrigOpt := proc() `CombineTrig/optimize`(args); end proc;



#==============================================================================
# Normalize an expression, which has the effect of removing all divisions
# but one, and without destroying subexpressions like normal() does.
#==============================================================================
normalize := proc( expr::{algebraic,list,set,`=`} )

    local p, denoms, cfactor, k, j, Z;

    if not type(expr,algebraic) then
        return map(procname,args);
    end if;
    if type(expr,`*`) or type(expr,`^`) then
        return map(procname,args);
    end if;
    if not type(expr,`+`) then
        return expr;
    end if;

    cfactor := denom(expr);
    if cfactor=1 then
        return expr;
    end if;

    Z := 0;
    for p in expr do
        Z := Z + procname(p)*cfactor;
    od;

    Z/cfactor;

end proc;



#==============================================================================
# collect and cancel common factors from both sides of an equation
#==============================================================================
eqcollect := proc( expr::equation )
	local L, R, S, x, p;
	L := lhs(expr);
	R := rhs(expr);
	S := indets(L,function) intersect indets(R,function);
	debug_print(procname,"S",2,S);
	x := factor(L-R);
	if type(x,`*`) then
		{op(x)} intersect S;
		for p in % do
			x := x / p;
		end do;
	end if;
	collect(x,convert(S,list));
	select(has,%,S) = -remove(has,%,S);
end proc;



#==============================================================================
# Calculate the average of expr over a circular annulus of
# inner radius a and outer radius b.  The coordinates used
# must be contained in the list coords, radius first then
# azimuth -- e.g., coords := [rho,phi].  The list paramlist
# is a list of ordering variables for collect.
#
#                        b  2*Pi
#                        /\ /\
#                1       |  |
# <expr> := ------------ |  | r * expr(r,phi) dphi dr
#           Pi*(b^2-a^2) |  |
#                       \/ \/
#                       r=a phi=0
#==============================================================================
annular_average := proc( expr, coords, a, b, paramlist )

    local I_phi, I_rho, I1, I2, expr_avg;

    debug_print( procname, `Calculating average over annulus...`, 1 );
    I_phi := int( expr, coords[2] );
    I_phi := subs( coords[2]=2*Pi, I_phi ) - subs( coords[2]=0, I_phi );
    I_rho := int( I_phi*coords[1], coords[1] );
    I1    := subs( coords[1]=b, I_rho );
    I2    := subs( coords[1]=a, I_rho );
    expr_avg := 1/Pi/(b^2-a^2)*(I1-I2);
    expr_avg := collect( expr_avg, paramlist, 'recursive', simplify );

    return expr_avg;
end proc;



#=============================================================================
# Determine if any of the subexpressions of expr are explicitly
# multiplied by I=sqrt(-1).
#
# Returns true if I does not explicitly appear anywhere in expr,
# false if it does.  WARNING: indexed quantities are treated as real,
# since Maple does not allow indexed names to be assumed.
#=============================================================================
is_real_parts := proc( expr )

    local p, result;
    option remember;

    for p in expr do
        if p = (-1)^(1/2) then
            return false;
        elif nops(p) > 1 then
            result := procname(p)
        else
            result := is(p,real)
        end if;
        if result = false then
            return false;
        elif result = FAIL then
            if is(p,indexed) then
                return true;
            end if;
            debug_print( procname,
                         `Unable to determine if `||(convert(p,string))||
                         ` is real or not!`, 0 );
            ERROR();
        fi
    od;

    return true;
end proc;



#=============================================================================
# Return only purely real solutions from solve()
#=============================================================================
realsols := proc( expr::algebraic, var::name )
    local sols, allsols, k;
    allsols := [solve(expr,var)];
    sols    := [];
    for k from 1 to nops(allsols) do
        if( is_real_parts(allsols[k]) ) then
            sols := [op(sols),var=allsols[k]];
        end if;
    od;
    sols;
#   remove( has, [solve(expr,var)], I );
end proc;



#=============================================================================
# Sort a solution set.
# 
# Example:
# > {a[6] = 3/2+1*sqrt(21)/2, a[5] = -3+sqrt(21), a[3] = -3-sqrt(21), 
#    a[2] = -3/2-1*sqrt(21)/2, a[1] = 2, a[4] = 1}:
# > sortsols(%,[seq(a[k],k = 1 .. 6)]);
# 
#   [ a[1] = 2, a[2] = -3/2 - 1/2 sqrt(21), a[3] = -3 - sqrt(21),
#     a[4] = 1, a[5] = -3 + sqrt(21), a[6] = 3/2 + 1/2 sqrt(21) ]
#=============================================================================
sortsols := proc( solset::{list(equation),set(equation)}, vars::list(name) )
    zip( (x,y)->x=y, vars, subs(solset,vars) );
end proc;




#=============================================================================
# "Combined Cost" of an Expression
#=============================================================================
combined_cost := proc( expr::{algebraic,equation,list,set,'array'(1)} )
    local c, cc, p, k;
    if type(expr,{equation,list,set}) then
        cc := 0;
        for p in expr do
            cc := cc + procname(p);
        od;
        return cc;
    end if;
    if type(expr,'array'(1)) then
        cc := 0;
        for k from 1 to size(expr) do
            cc := cc + procname(expr[k]);
        od;
        return cc;
    end if;
    c  := cost(expr);
    cc := 0;
    #single--e.g., `2*additions, or `divisions`
    if type(c,`*`) or nops(c)=1 then
        if nops(c)=1 then
            if c=`multiplications` or c=`divisions` or c=`additions` or c=`functions` then
                cc := 1;
            end if;
        else
            if op(2,c)=`multiplications` or op(2,c)=`divisions` or op(2,c)=`additions` or 
            op(2,c)=`functions` then
                cc := op(1,c);
            end if;
        end if;
    else
        for p in c do
            if nops(p)=1 then
                if p=`multiplications` or c=`divisions` or p=`additions` or p=`functions` then
                    cc := cc + 1;
                end if;
            else
                if op(2,p)=`multiplications` or op(2,p)=`divisions` or op(2,p)=`additions` or 
                op(2,p)=`functions` then
                    cc := cc + op(1,p);
                end if;
            end if;
        od;
    end if;
    cc;
end proc;


#=============================================================================
# Sort a List of Expressions in Ascending Order by Complexity
#=============================================================================
sort_by_cost := proc( exprs::{list,set} )
	local costsort, z;
	costsort := proc(x)
		sort( [op(x)], 
			  (a,b)->if utils[combined_cost](a) 
					  < utils[combined_cost](b)
					 then true else false fi 
			);
	end proc;
	z := exprs;
	if type(z,{list(list),list(set)}) then
		z := map( procname, exprs );
	end if;
	costsort(z);
end proc;




#=============================================================================
# Return a 3x3 matrix which will rotate a vector ccw around an axis.
# Positive rotation is ccw when viewing from the positive axis direction
# toward the negative (right-hand rule).
#=============================================================================
rotaxis := proc( theta::algebraic, axis::string )
    local R;
    if axis=`x` or axis=`X` then
        R := array( 1..3, 1..3, [
                    [ 1,          0,           0 ],
                    [ 0, cos(theta), -sin(theta) ],
                    [ 0, sin(theta),  cos(theta) ]
                   ] );
     elif axis=`y` or axis=`Y` then
        R := array( 1..3, 1..3, [
                    [  cos(theta), 0, sin(theta) ],
                    [           0, 1,          0 ],
                    [ -sin(theta), 0, cos(theta) ]
                   ] );
     elif axis=`z` or axis=`Z` then
        R := array( 1..3, 1..3, [
                    [ cos(theta), -sin(theta), 0 ],
                    [ sin(theta),  cos(theta), 0 ],
                    [          0,           0, 1 ]
                   ] );
     else
         ERROR(`:  axis must be "x", "y", or "z"`);
     end if;
     return op(R);
end proc;



#=============================================================================
# Rotate a 3D vector by angle phi around the axis defined by w.  Positive
# rotation is defined by the right-hand rule (clockwise when viewing
# along w).
#
# Neither w nor vec needs to be a unit vector.
#=============================================================================
rotvec := proc( v::{vector,matrix}, w::{vector,matrix}, phi::algebraic )

    local   newvec, w_hat, V, W;

    if type(v,'array'(1..3,1..1)) then
        V := convert(eval(v),vector);
    else
        V := v;
    end if;
    if type(w,'array'(1..3,1..1)) then
        W := convert(eval(w),vector);
    else
        W := w;
    end if;
    if nops(eval(V))<>3 or nops(eval(W))<>3 then
        ERROR(`: input vectors must contain 3 elements only`);
    end if;

    w_hat  := evalm(W/mag(W));
    newvec := eval(V)*cos(phi) + (1-cos(phi))*(dot(V,w_hat))*eval(w_hat)
              - cross(V,w_hat)*sin(phi);

    if type(v,matrix) then
        eval(mat(evalm(newvec)));
    else
        eval(newvec);
    end if;
end proc;



#=============================================================================
# Apply func to the specified operands (as a group) of the expression.
#=============================================================================
funcops := proc( expr::`+`, func::procedure, oplist::{integer,list} )
    local part1, part2, i, k;
    part1 := 0;
    for i from 1 to nops(oplist) do
        if type(oplist,integer) then
            k := oplist;
        else
            k := oplist[i];
        end if;
        part1 := part1 + op(k,expr);
    od;
    part2 := expr - part1;
    part1 := func( part1, args[4..nargs] );
    return part2 + part1;
end proc;



#=============================================================================
# Complete the square in the (optional) variable var in the (2nd degree in
# var) polynomial expr.
#
# var can be a name, a function, a `+`, `*`, or `^` expression, or a
# set or list of functions, variable names, and `+`, `*`, `^` expressions.
#
# An additional argument -- a Maple procedure name -- may be included and,
# if present, will be applied to the remainder terms.  For example,
#
#    csquare( expr, [x,y], factor );
#=============================================================================
csquare := proc( expr::{algebraic,array,indexed,equation,numeric} )

    local   A, B, C, newexpr, var, t, i, L, R,
            subsflag, savevar, srvars, F, func, _H_, _arg, _v, _srv, _F;

    if type(expr,'numeric') then
        return expr;
    end if;

    #
    # grab the indeterminate if no second arg; check for function to apply
    #
    func := NULL;
    if nargs=1 then
        if nops(indets(expr))=0 then
            return expr;
        elif nops(indets(expr))=1 then
            var := op(indets(expr));
        else
            ERROR( `: Can't determine the indeterminate!` );
        end if;
    elif nargs=2 then
        if type(args[2],procedure) then
            func := args[2];
        else
            var := args[2];
        end if;
    elif nargs=3 then
        if type(args[2],procedure) then
            func := args[2];
            var  := args[3];
        elif type(args[3],procedure) then
            func := args[3];
            var  := args[2];
        else
            ERROR(`: wrong argument type(s)`);
        end if;
    else
        ERROR(`: wrong number arguments!`);
    end if;
    if not type(var,{name,set,list,`+`,`*`,`^`,function}) then
        ERROR( `: var must be one of `
               ||`{name,set,list,``+``,``*``,``^``,function}.` );
    end if;

    newexpr := copy(expr);
    if type(var,{set,list}) then
        for t in var do
            newexpr := procname( newexpr, t, func );
        od;
        return eval(newexpr);
    end if;

    #
    # complete the square
    #
    _H_      := '_H_';
    subsflag := OFF;                    #used to flag algebraic var
    if type( newexpr, {array,indexed,equation} ) then
        if type(var,{`+`,`*`,`^`,function}) then
            subsflag := ON;
            savevar  := var;
            var      := _H_;
            if type(savevar,function) then
                srvars := [savevar];
            else
                srvars := [op(indets(savevar))];
            end if;
            #
            # Define a mapping function to apply to each element of newexpr.
            # F does side relation simplification so as to substitute the
            # dummy variable _H_ for whatever expression is contained in var.
            #
            _F := proc(_arg,_v,_srv, _H_)
                if has(_arg,_v) then
                    simplify( _arg, {_v=_H_}, _srv );
                else
                    _arg;
                end if;
            end proc;
            newexpr := map( _F, newexpr, savevar, srvars, _H_ );
        end if;
        newexpr := map( procname, newexpr, var, func );
        if subsflag=ON then
            newexpr := subs( _H_=savevar, eval(newexpr) );
        end if;
    elif type( newexpr, polynom(anything,var) ) then
        if type(var,{`+`,`*`,`^`,function}) then
            subsflag := ON;
            savevar  := var;
            var      := _H_;
            newexpr  := simplify( newexpr, {savevar=_H_}, [op(indets(savevar))] );
        end if;
        newexpr := collect( newexpr, var );
        A := coeff( newexpr, var, 2 );  #form: A*var^2 + B*var + C
        B := coeff( newexpr, var, 1 );
        if A <> 0 and B <> 0 then
            C := coeff( newexpr, var, 0 );
            if func <> NULL and type(func,procedure) then
                newexpr := A*(var + B/(2*A))^2 + func(C-(B^2/(4*A)) );
            else
                newexpr := A*(var + B/(2*A))^2 + (C - (B^2/(4*A)) );
            end if;
        end if;
        if subsflag=ON then
            newexpr := subs( _H_=savevar, newexpr );
        end if;
    else
        return eval(expr);
    end if;

    return eval(newexpr);
end proc;



#=============================================================================
# Apply a function to an expression n times.
#=============================================================================
applyfunc := proc( expr, func::{procedure,function}, n::nonnegint )
    if n=0 then
        expr;
    else
        if nargs > 3 then
            func( procname(expr,func,n-1,args[4..nargs]), args[4..nargs] );
        else
            func( procname(expr,func,n-1) );
        end if;
    end if;
end proc;



#=============================================================================
# Apply func to individual (top-level) terms in expression expr.
#=============================================================================
termfunc := proc( expr::{algebraic,array(1),`=`}, func::procedure )

    local   tlist, N, i, newexpr, tmp;

    newexpr := copy( expr );

    if type(newexpr,{'array'(1),`=`}) then
        return map( procname, newexpr, func, args[3..nargs] );
    end if;

    if type(newexpr,{`+`,`*`}) then
        tlist := map( func, [op(newexpr)], args[3..nargs] );
        if type(newexpr,`*`) then
            newexpr := convert( tlist, `*` );
        else
            newexpr := convert( tlist, `+` );
        end if;
    elif type(newexpr,`^`) then
        N := nops(newexpr);
        for i from 1 to N do
            newexpr := subsop( [i]=procname(op(i,newexpr),func,args[3..nargs]),
                               newexpr );
            if newexpr=1 then break end if;
        od;
    else
        newexpr := func( newexpr, args[3..nargs] );
    end if;

    gc();
    return newexpr;
end proc;



#=============================================================================
# Selectively factor out a subexpression subexpr from an expression expr.
#=============================================================================
pullout := proc( expr::{name,`+`,`^`,`*`,`=`,function},
                         subexpr::{name,function,`+`,`*`,`^`,list} )

    local newexpr1, newexpr2, p, q, flag, r, allpunts;

    if type(expr,name) then
        return expr;
    end if;

    #
    # special case: expr is an equation
    #
    if type(expr,`=`) then
        return map(procname,expr,subexpr);
    end if;

    #
    # special case: we have a list of subexpressions -- do each in order
    #
    if type(subexpr,list) then
        newexpr1 := expr;
        for p in subexpr do
            if not type(p,{name,function,`+`,`^`,`*`}) then
                ERROR(`: subexpression `||(convert(p,string))||` not of type name, function, +, ^, or *`);
            end if;
            newexpr1 := procname( newexpr1, p );
        od;
        return newexpr1;
    end if;

    #
    # special case: expr is type function
    #
    if type(expr,function) then
        return expr;
    end if;

    #
    # special case: expr is type `*`
    #
    if type(expr,`*`) then
        newexpr1 := expr;
        #
        # make sure expr contains the entire subexpr
        #
        flag := ON;
        if type(subexpr,{`+`,`^`}) then
            if not has(expr,subexpr) then
                flag := OFF;
            end if;
        else
            for q in subexpr do
                if not has(expr,q) then
                    flag := OFF;
                    break;
                end if;
            od;
        end if;
        if flag=ON then
            #
            # factor the first term of expr that contains subexpr
            #
            for p in expr do
                if type(p,{`+`,`^`}) then
                    newexpr1 := expr/p * procname(p,subexpr);
                    break;
                end if;
            od;
        end if;
        return newexpr1;
    end if;

    #
    # special case: expr is type `^`
    #
    if type(expr,`^`) then
        newexpr1 := expr;
        if type(op(1,expr),`+`) then
            newexpr1 := procname(op(1,expr),subexpr)^op(2,expr);
        end if;
        return newexpr1;
    end if;
    if not type(expr,`+`) then
        return expr;
    end if;

    #
    # general case: expr is type `+`
    #
    # first, some loop setup
    #
    newexpr1 := 0;
    newexpr2 := 0;
    allpunts := YES;
    #
    # loop over the parts of expr
    #
    for p in expr do
        #
        # first make sure each part contains the entire subexpr
        #
        flag := ON;
        if type(subexpr,{name,function,`+`,`^`}) then
            if not has(p,subexpr) then
                flag := OFF;
            end if;
        else    # is a `*`
            for q in subexpr do
                if not has(p,q) then
                    flag := OFF;
                    break;
                end if;
            od;
        end if;
        #
        # now pull out the subexpr from those terms of
        # expr that contain it
        #
        if flag=ON then
            if type(p,{`*`,`^`}) then
                r := procname(p,subexpr);
                if factor(r-p)=0 then
                    newexpr1 := newexpr1 + r/subexpr;
                else
                    newexpr1 := newexpr1 + r;
                end if;
                allpunts := NO;
            #
            # no dice -- punt on this p
            #
            else
                newexpr1 := newexpr1 + (p/subexpr);
            end if;
        #
        # punt on this p since it doesn't contain subexpr
        #
        else
            newexpr2 := newexpr2 + p;
        end if;
    od;

    if allpunts=YES then
        expr;
    else
        newexpr2 + subexpr*newexpr1;
    end if;

end proc;



#=============================================================================
# convert a matrix equation of the form mat(x,y,z)=mat(eqx,eqy,eqz) to a
# list of equations [x=eqx, y=eqy, z=eqz]
#=============================================================================
mat2eqs := proc( Meqn::equation )
    local k, eqlist, Mleft, Mright;
    Mleft  := evalm(lhs(Meqn));
    Mright := evalm(rhs(Meqn));
    eqlist := [];
    for k from 1 to rows(Mleft) do
        eqlist := [ op(eqlist), Mleft[k,1]=Mright[k,1] ];
    od;
    eval(eqlist);
end proc;



#=============================================================================
# factor common terms from a vector or matrix and pull them out front
#=============================================================================
factormat := proc( expr::anything )
	local M, f, cf, s;
	if type(expr,{`=`,`+`,`*`,`list`,`set`}) then
		map( procname, expr, args[2..nargs] );
	elif type(expr,`array`) then
		#extract the common factors
		M  := map(factor,expr);
		cf := `intersect`( seq( {op(M[k])}, k=1..size(M) ) ) minus {0};
		#the factors that will go out in front
		s  := convert(%,`*`);
		#set each common factor to one in the original expression
		cf := map( x->x=1, cf );
		cf := convert(cf,list);
		#we're done
		s*subs( cf, M );
	else
		expr;
	end if;
#
#old version
#
#    local S, N, Q, r, c, k;
#    if type(M,{`=`,`+`,`*`,`list`,`set`}) then
#        return map( procname, M, args[2..nargs] );
#    elif type(M,`array`) then
#        S := {};
#        if type(M,`vector`) then
#            for r from 1 to size(M) do
#                S := S union {M[r]};
#            od;
#        else
#            for r from 1 to rows(M) do
#                for c from 1 to cols(M) do
#                    S := S union {M[r,c]};
#                od;
#            od;
#        end if;
#        S := S minus {0};
#        S := map(factor,S);
#        Q := {op(S[1])};
#        for k from 2 to nops(S) do
#            Q := Q intersect {op(S[k])}
#        od;
#        Q := convert(Q minus {-1},`*`);
#        N := copy(M);
#        N := map( (x,y)->x/y, N, Q );
#        if nargs > 1 then
#            N := map( args[2], N, args[3..nargs] );
#            Q := map( args[2], Q, args[3..nargs] );
#        end if;
#        Q * eval(N);
#    else
#        M;
#    end if;
end proc;



#=============================================================================
# Test expr for the presence of a restricted form of polynomial in
# any of the variables in varlist.  Returns 'true' or 'false'.
#
# The restriction in the definition of "polynomial" is that the
# expression or any of the parts of the expression must contain 2
# or more terms separated by an addition operator ('+' or '-'), at
# least one of which contains a power of at least one of the
# variables in varlist.  For example, suppose varlist:=X.  Then
#
#      2                  3
#   a*X    and   (a + b) X
#
# are NOT "polynomials" as far as polytest is concerned, but
#
#      2                        3
#   a*X  + b   and   (a + b X) X
#
# are.  polytest() is designed to be used by subxfunc(), which
# makes use of this restriction.
#=============================================================================
polytest := proc( expr::algebraic, varlist::{list,string} )

    local p;

    if type( expr, `+` ) then
        for p in varlist do
            if has(expr,p) and degree(expr,p) <> 0 then
                return true;
            end if;
        od;
    elif type( expr, {`*`,`^`} ) then
        for p in expr do
            if polytest(p,varlist) = true then
                return true;
            end if;
        od;
    end if;

    false;

end proc;



#=============================================================================
# Check an expression or a list of expressions for small divisors in
# the variables given in varlist.
#
# For example, if varlist:=[x,b] and A:=(a+x)/(x*c+b^2), then
#
#   small_divisors( A, varlist ) = true
#
# But if varlist:=[x,c], then
#
#   small_divisors( A, varlist ) = false
#
# The argument varlist is a list (e.g., [theta,f,C]) or a vector.
#
# small_divisors() makes use of common_factor().
#=============================================================================
small_divisors := proc( expr::{algebraic,list(algebraic),'array'(1),`=`},
                                varlist::{list,name} )

    local den, num, count, _H_, p, equ;

    assume( _H_, real );

    #if expr is a list of equations, do each equation
    if not type(expr,algebraic) then
        for equ in expr do
            if procname(equ,varlist) then
                return true;
            end if;
        od;
        return false;
    end if;

    count := 0;
    for p in varlist do
        if not has( expr, p ) then
            count := count + 1;
        end if;
    od;
    if count = nops(varlist) then       #no instances of varlist variables
        return false;
    end if;

    equ := normal( expr );

    #multiply all varlist variables by a dummy variable _H_ to catch cross terms
    for p in varlist do
        equ := subs( p=_H_*p, equ );
    od;

    for p in equ do
        den := denom( p );
        num := numer( p );
        if den <> 1 then
            if common_factor( den, _H_ ) = true then
                return true;
            end if;
        end if;
        if denom(num) <> 1 then
            if procname( num, _H_ ) = true then
                return true;
            end if;
        end if;
    od;

    false;

end proc;



#=============================================================================
# Utility function for small_divisors().
#
# Check an expression for a common multiplicative factor in the variables
# contained in varlist.  It is assumed that negative powers in the list
# variables have been removed (e.g., by feeding common_factor the numerator
# or denominator of an equation in reduced normal form).
#
# For example, if varlist:=x, then
#
#   common_factor( a*(1+x^2+x^3)+b*x, varlist ) = false
#
# But if varlist:=[x,a], then
#
#   common_factor( a*(1+x^2+x^3)+b*x, varlist ) = true
#
# The argument varlist is a Maple list (e.g., [theta,f,C]).
#=============================================================================
common_factor := proc( expr::algebraic, varlist::{list,name} )

    local equ, _H_, p;
    assume( _H_, real );

    equ := copy( expr );

    #multiply all variables in varlist by dummy variable _H_ to catch
    #cross terms
    for p in varlist do
        equ := subs( p=_H_*p, equ );
    od;

    if ldegree( collect(equ,_H_), _H_ ) > 0 then
        return true;
    end if;

    return false;
end proc;



#=============================================================================
# Expand 'expression' into a series up to and including terms of order
# 'expansion_order' in the small parameters 'params'.
#
# If _order < 1, then no expansions are performed.
#
# params is a Maple list (e.g., [theta,f,C]).
#
# Multiple independent expansions can be done in one step.  For example,
#
#   expansion( expr, theta, 3, [epsilon,psi,alpha], 2, [phi,beta], 5 );
#
# NOTE: there is a Maple bug such that sometimes s is not characterized
# as a "series" type if there are noninteger powers of params
# involved in expression. I haven't been able to figure out a
# workaround for this.  I have tried convert(s,series,params),
# which fails.  When the bug appears, conversion of s to type
# polynom does not remove the O() term.
#=============================================================================
# Return Value
# ------------
# The series expansion of expression, converted to the
# Maple type polynom.
#=============================================================================
expansion := proc( expression::{list,set,'array'(1),'array'(2),matrix,`=`,algebraic},
                   params::{set,name,list,function}, expansion_order::integer )

    local   s, s1, ord1, _H_, p, plist, degree_s1, i;

    if type(nargs,even) then
        ERROR(`: wrong number of arguments`);
    end if;

    if not type( eval(expression), algebraic ) then
        return map( procname, args );
    end if;

    if not has(expression,params) then
        return expression;
    end if;

    #
    # Handle multiple independent expansions
    #
    if nargs > 3 then
        s := procname( expression, params, expansion_order );
        for i from 1 to (nargs-3)/2 do
            s := procname( s, args[2*i+2], args[2*i+3] );
        od;
        return eval(s);
    end if;


    if expansion_order < 0 then
        debug_print( procname, `WARNING!  Expansion order is less than 0`, 4 );
        debug_print( procname, `Expression NOT expanded:`, 5, expression );
        return expression;
    end if;

    if expansion_order = 0 then
        return eval( expression, map((x)->x=0,params) );
    end if;

    assume( _H_, real );
    s    := copy(expression);
    ord1 := expansion_order + 1;

    #
    # Expand in a series in the small parameters, taking into account
    # cross terms by expanding on a dummy variable _H_.
    #
    if not type(params,list) then
        plist := [params];
    else
        plist := params;
    end if;
    for p in plist do
        s := subs( p=_H_*p, s );
    od;

    convert( series(s,_H_,ord1), polynom );
    s1 := collect( %, [_H_,op(params)] );

    #
    # Version V release 4 of Maple does not expand the expression
    # in a series with the intended *output* order.  (It only
    # guarantees *internal* calculations to the correct order.)
    # Hence sometimes we must keep increasing the order until the
    # degree of the output polynomial is correct.
    #
    degree_s1 := degree(s1,_H_);
    if degree_s1 = FAIL then
        debug_print( procname, `WARNING!  Unable to determine degree `
                     ||`of expanded expression.`, 2 );
        debug_print( procname, `Problem expression:`, 2, s1 );
    else
        #
        # The FAIL test is so that expressions whose degree in {params} is
        # already less than _order won't trigger the while loop.
        #
        if degree(s,convert(params,set))=FAIL and nops(s1) > 1
        and degree_s1 < expansion_order then
            while degree_s1 < expansion_order do
                ord1 := ord1 + 1;
                debug_print( procname, `WARNING!  Increasing expansion order to `
                             ||(convert(ord1-1,string))||`...`, 4 );
                convert( series(s,_H_,ord1), polynom );
                s1 := collect( %, [_H_,op(params)] );
                degree_s1 := degree(s1,_H_);
                if degree_s1 = FAIL then
                    debug_print( procname, `WARNING!  Unable to determine degree `
                                 ||`of expanded expression.`, 2 );
                    debug_print( procname, `Problem expression:`, 2, s1 );
                    break;
                end if;
                if ord1 > expansion_order+7 then
                    debug_print(procname,`WARNING! Too many order increases!`,4,ord1);
                    debug_print( procname, `Problem expression:`, 5, s1 );
                    break;
                end if;
            od;
        end if;

        #
        # For some expressions (like stuff under square roots), Maple's
        # series() expands beyond the power asked for by 2 or more orders.
        # This seems to be a bug.  Redoing the expansion a second time
        # produces the correct expansion.
        #
        if degree_s1 > expansion_order then
            debug_print( procname, `WARNING!  Encountered series() bug...`, 4 );
            debug_print( procname, `Buggered expression:`, 5, s1 );
            s1 := convert( series(s,_H_,ord1), polynom );
        end if;
    end if;

    #
    # Maple gives a bogus O(1) for the series expansion
    # under certain obscure conditions.  A workaround is to force
    # it to find the leading term; then something about the
    # execution history that Maple remembers will trigger the
    # correct series expansion.  It's a time waster, but it's
    # the only way I can think of for getting around the problem.
    #
    if s1=0 then
        debug_print( procname, `WARNING!  Fixing O(1) bug in series()...`, 4 );
        series( leadterm(s), _H_, 999 );
        s1 := convert( series(s,_H_,ord1), polynom );
    end if;

    subs( _H_=1, s1 );
    s := collect( %, params );

    gc();
    return s;
end proc;


end module:

#==============================================================================
savelibname := LIB_PATH:
savelib(utils):
savelibname := 'savelibname':
#==============================================================================