reduce-historical/r38/packages/odesolve/odetop.red

module odetop$  % Top level ODESolve routines, exact ODEs, general
                % nonlinear ODE simplifications and utilities

% F.J.Wright@maths.qmw.ac.uk, Time-stamp: <11 August 2001>

% TO DO:
%    allow for non-trivial denominator in exact ODEs

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Code to support hooks for extending the functionality of ODESolve
% without needing to edit the main source code.  (Based on the hooks
% supported by Emacs.)  [Note that in CSL, each hook must be declared
% global (or fluid), even for interactive testing, otherwise boundp
% does not work as expected!]

% To do: Run hooks within an errorset for extra security?

symbolic procedure ODESolve!-Run!-Hook(hook, args);
   %% HOOK is the *name* of a hook; ARGS is a *list* of arguments.
   %% If HOOK is a function or is bound to a function then apply it to
   %% ARGS; if HOOK is bound to a list of functions then apply them in
   %% turn to ARGS until one of them returns non-nil and return that
   %% value.  Otherwise, return nil.
   if getd hook then apply(hook, args)
   else if boundp hook then <<
      hook := eval hook;
      if atom hook then
         getd hook and apply(hook, args)
      else
      begin scalar result;
         while hook and null(
            getd car hook and
            (result:=apply(car hook, args))) do
               hook := cdr hook;
         if hook then return result
      end
   >>$

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Code to break ODE simplifier loops
% ==================================

fluid '(odesolve!-interchange!-list!* !*odesolve!-norecurse)$

global '(ODESolve!-Standard!-x ODESolve!-Standard!-y)$

ODESolve!-Standard!-x := gensym()$
ODESolve!-Standard!-y := gensym()$

symbolic procedure ODESolve!-Standardize(ode, y, x);
   %% Return the numerator of ode in true prefix form and with
   %% standardized variable names.  (What about sign, etc.?)
   subst(ODESolve!-Standard!-y, y,
      subst(ODESolve!-Standard!-x, x, prepf numr simp!* ode))$

symbolic procedure ODESimp!-Interrupt(ode, y, x);
   begin scalar std_ode;
      ode := num !*eqn2a ode;           % Returns ode as expression.
      if member(std_ode:=ODESolve!-Standardize(ode, y, x),
         odesolve!-interchange!-list!*) then <<
            traceode "ODE simplifier loop interrupted! ";
            return t
         >>;
      odesolve!-interchange!-list!* := std_ode .
         odesolve!-interchange!-list!*
   end$


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Top-level classification of an ODE, primarily as linear or
% nonlinear.

global '(ODESolve_Before_Hook ODESolve_After_Hook)$
global '(ODESolve_Before_Non_Hook ODESolve_After_Non_Hook)$

algebraic procedure ODESolve!*0(ode, y, x);
   %% Top-level general ODE solver.  If no derivatives call solve?
   %% ***** DO NOT CALL RECURSIVELY *****
   symbolic begin scalar !*precise, solution, !*odesolve!-norecurse,
         odesolve!-interchange!-list!*, !*odesolve!-solvable!-xy;
      %% (odesolve!-interchange!-list!* and !*odesolve!-solvable!-xy
      %% are used to prevent infinite loops.)
      ode := num !*eqn2a ode;           % returns ode as expression
      if (solution := or(
         ODESolve!-Run!-Hook('ODESolve_Before_Hook, {ode,y,x}),
         ODESolve!*1(ode, y, x),
         %% Call ODESolve!-Diff once only, not in recursive loop?
         %% SHOULD apply only to nonlinear ODEs?
         not !*odesolve_fast and ODESolve!-Diff(ode, y, x),
         ODESolve!-Run!-Hook('ODESolve_After_Hook, {ode,y,x})))
      then return solution;
      traceode "ODESolve cannot solve this ODE!"
   end$

algebraic procedure ODESolve!*1(ode, y, x);
   %% Top-level discrimination between linear and nonlinear ODEs.
   %% May be called recursively.
   %% (NB: A product of linear factors is NONLINEAR!)
   symbolic if !*odesolve!-norecurse then
      traceode "ODESolve terminated: no recursion mode!"
   else if ODESimp!-Interrupt(ode, y, x) then nil else
   <<
      !*odesolve!-norecurse := !*odesolve_norecurse;
      traceode1 "Entering top-level general recursive solver ...";
      if ODE!-Linearp(ode, y) then      % linear
         ODESolve!-linear(ode, y, x)
      else  % nonlinear -- turn off basis solution
      algebraic begin scalar !*odesolve_basis, ode_factors, solns;
         %% Split into algebraic factors (which may lose exactness).
         %% For each algebraic factor, check its linearity and call
         %% appropriate (linear or nonlinear) main solver.
         %% Merge solution sets.
         traceode1 "Trying to factorize nonlinear ODE algebraically ...";
         ode_factors := factorize ode;
         %% { {factor, multiplicity}, ... }
         if length ode_factors = 1 and second first ode_factors = 1 then
            %% Guaranteed algebraically-irreducible nonlinear ODE ...
            return ODESolve!-nonlinear(ode, y, x);
         traceode "This is a nonlinear ODE that factorizes algebraically ",
            "and each distinct factor ODE will be solved separately ...";
         solns := {};
         while ode_factors neq {} do
         begin scalar fac;
            %% Discard repeated factors:
            if smember(y, fac := first first ode_factors) then
               %% Guaranteed algebraically-irreducible -- may be
               %% either algebraic or linear or nonlinear ODE ...
               if (fac := ODESolve!*2!*(fac, y, x)) then <<
                  solns := append(solns, fac);
                  ode_factors := rest ode_factors
               >>  else solns := ode_factors := {}
            else <<
               if depends(fac, x) or depends(fac, y) then
                  symbolic MsgPri("ODE factor", fac, "ignored", nil, nil);
               ode_factors := rest ode_factors
            >>;
         end;
         %% Finally check whether the UNFACTORIZED ode was exact:
         return
            if solns = {} then Odesolve!-Exact!*(ode, y, x)
            else solns
      end
   >>$

algebraic procedure ODESolve!-FirstOrder(ode, y, x);
   %% Solve an ARBITRARY first-order ODE.
   %% (Called from various other modules.)
   symbolic <<
      ode := num !*eqn2a ode;
      traceode ode = 0;
%%       if ODE!-Linearp(ode, y)        % nil <> 0 !!!
%%       then ODENon!-Linear1(ode, y, x)
%%       else ODESolve!-NonLinear1(ode, y, x)
      %% A nonlinear first-order ODE may need the full solver ...
      %% but could later arrange to pass the order rather than
      %% recompute it.
      ODESolve!*1(ode, y, x)
   >>$

algebraic procedure ODESolve!*2!*(ode, y, x);
   %% Internal discrimination between algebraic or differential factor.
   if smember(df, ode) then             % ODE
      ODESolve!*2(ode, y, x)
   else if ode = y then                 % Common special algebraic case,
      {y = 0}                           % e.g. solving autonomous ODEs.
   else solve(ode, y)$                  % General algebraic case.

algebraic procedure ODESolve!*2(ode, y, x);
   %% Internal discrimination between linear and nonlinear ODEs.  Like
   %% ODESolve!*1 but does not attempt any algebraic factorization.
   symbolic <<
      traceode1 "Entering top-level recursive solver ",
         "without algebraic factorization ...";
      traceode ode=0;
      if ODE!-Linearp(ode, y) then      % linear
         ODESolve!-linear(ode, y, x)
      else                              % nonlinear
         ODESolve!-nonlinear(ode, y, x)
   >>$


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% The entry point to the non-trivially nonlinear ODE solver
% =========================================================

algebraic procedure ODESolve!-nonlinear(ode, y, x);
   %% Attempt to solve an algebraically-irreducible nonlinear ODE.
   symbolic %% if ODESimp!-Interrupt(ode, y, x) then nil else
      begin scalar ode_order;
         ode_order := ODE!-Order(ode, y);
         traceode "This is a nonlinear ODE of order ", ode_order, ".";
         return or(
            ODESolve!-Run!-Hook(
               'ODESolve_Before_Non_Hook, {ode,y,x,ode_order}),
            (if ode_order = 1 then
               ODESolve!-nonlinear1(ode, y, x)
            else
               %% ODESolve!-Diff(ode, y, x) or % TEMPORARY
               ODESolve!-nonlinearn(ode, y, x)),
            ODESolve!-Exact(ode, y, x, ode_order),
            not !*odesolve_fast and ODESolve!-Alg!-Solve(ode, y, x),
            not !*odesolve_fast and ODESolve!-Interchange(ode, y, x),
            ODESolve!-Run!-Hook(
               'ODESolve_After_Non_Hook, {ode,y,x,ode_order}))
      end$

symbolic procedure ODESolve!-Interchange(ode, y, x);
   %% Interchange x <--> y and try to solve.
   %% PROBABLY NOT DESIRABLE FOR LINEAR ODES!
   if !*odesolve_noswap then
      traceode "ODESolve terminated: no variable swap mode!"
   else
   ( begin scalar !*precise;            % Can cause trouble here
      traceode
         "Interchanging dependent and independent variables ...";
      %% Should fully canonicalize ode before comparison!!!
      %% Temporarily, just use reval to at least ensure the same format.
      %% Cannot use aeval form because simplified flag gets reset.
%%       odesolve!-interchange!-list!* :=
%%          %% reval ode . odesolve!-interchange!-list!*;
%%          ODESolve!-Standardize(ode, y, x) . odesolve!-interchange!-list!*;
      depend1(x, y, t);
      algebraic begin scalar rules;
         rules := {odesolve!-df(y,x,~n) => 1/odesolve!-df(x,y)*
            odesolve!-df(odesolve!-df(y,x,n-1),y) when n > 1,
            odesolve!-df(y,x) => 1/odesolve!-df(x,y),
            odesolve!-df(y,x,1) => 1/odesolve!-df(x,y)};
         ode := sub(df = odesolve!-df, ode);
         ode := (ode where rules);
         ode := num sub(odesolve!-df = df, ode)
      end;
      depend1(y, x, nil);   % Necessary to avoid dependence loops
      %% Now ode is an ode for x as a function of y
      traceode ode;
%%       %% if member(reval ode, odesolve!-interchange!-list!*) then
%%       if member(ODESolve!-Standardize(ode, x, y),
%%          odesolve!-interchange!-list!*) then
%%          %% Give up -- we have already interchanged variables in this
%%          %% ode once!
%%          << !*odesolve_failed := t;  return algebraic {ode=0} >>;
      ode := ODESolve!*1(ode, x, y);    % Try again ..
      if ode then return
         makelist for each soln in cdr ode join
            if smember(y, soln) then {soln} else {}
   end ) where depl!* = depl!*$


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Exact equations
% ===============

% Solve an ODE if it is an exact first or second order ODE.  Exactness
% might be lost by factorizing an ode, so all exact ode routines are
% gathered together here under one master routine that can be called
% independently of any ODE simplification.

% Replace by one general routine for any order nonlinear ODE?

% The first-order code is based on code by Malcolm MacCallum.

algebraic procedure ODESolve!-Exact!*(ode, y, x);
   %% Solve an exact first or second order nonlinear ODE of unknown
   %% order.
   ODESolve!-Exact(ode, y, x, ODE!-Order(ode, y))$

algebraic procedure ODESolve!-Exact(ode, y, x, ode_order);
   %% Solve an exact first or second order nonlinear ODE of known
   %% order.
   begin scalar c, den_ode, result;
      traceode1 "Checking for an exact ode ...";
      c := coeff(num ode, df(y,x,ode_order));
      den_ode := den ode;
      if not depends(den_ode, x) then den_ode := 0;
      %% ... meaning den ode has no effect on exactness.
      %% if length c neq 2 or depends(c, df(y,x,n)) then return;
      %% NB: depends recurses indefinitely if x depends on y, i.e. after
      %% interchange at present.  But smember nearly suffices anyway!
      if length c neq 2 or smember(df(y,x,ode_order), c) then return;
      return if ode_order = 1 then
         symbolic ODESolve!-Exact!-1(c, den_ode, y, x)
      else if ode_order = 2 then
         symbolic ODESolve!-Exact!-2(c, den_ode, y, x)
   end$

symbolic procedure ODESolve!-Exact!-1(c, den_ode, y, x);
   %% Solves the ode if it is an exact (nonlinear) first order ode of
   %% the form = N dy/dx + M.
   ( algebraic begin scalar M, N;
      M := first c;  N := second c;
      symbolic depend1(y, x, nil);      % all derivatives partial
      if df(M,y) - df(N,x) and
         (not den_ode or df(M:=M/den_ode,y) - df(N:=N/den_ode,x))
      then return;
      %% traceode "This is an exact first-order ODE.";
      traceode "It is exact and is solved by quadrature.";
      return {exact1_pde(M, N, y, x) = 0}
   end ) where depl!* = depl!*$

algebraic procedure exact1_pde(M, N, y, x);
   %% Return phi(x,y) such that df(phi,x) = M(x,y), df(phi,y) =
   %% N(x,y), required to integrate first and second order exact odes.
   begin scalar int_M;  int_M := int(M, x);
      %% phi = int_M + f(y)
      %% => df(phi,y) = df(int_M,y) + df(f,y) = N
      %% => f = int(N - df(int_M,y), y)
      return num(int_M + int(N - df(int_M,y), y) + newarbconst())
   end$

symbolic procedure ODESolve!-Exact!-2(c, den_ode, y, x);
   %% Computes a first integral of ODE if it is an exact (nonlinear)
   %% second order ODE of the form f(x,y,y') y'' + g(x,y,y') = 0.
   %% *** EXTEND THIS GENERAL CODE TO HIGHER ORDER ??? ***
   ( algebraic begin scalar p, f, g, h, first_int, h_x, h_y;
      p := gensym();
      f := sub(df(y,x) = p, second c);
      g := sub(df(y,x) = p, first c);
      symbolic depend1(y, x, nil);      % all derivatives partial
      if ODESolve!-Exact!-2!-test(f, g, p, y, x)
         and (not den_ode or
            ODESolve!-Exact!-2!-test(f:=f/den_ode, g:=g/den_ode, p, y, x))
      then return;
      %% ODE is exact
      %% traceode "This is an exact second-order ODE for which ",
      %%   "a first integral can be constructed:";
      traceode "It is exact and a first integral can be constructed ...";
      h := gensym();
      symbolic depend1(h, x, t);  symbolic depend1(h, y, t);
      first_int := int(f, p) + h;
      c := df(first_int,x) + df(first_int,y)*p - g; % = 0
      %% Should be linear in p by construction -- equate coeffs:
      c := coeff(num c, p);
      if length c neq 2 or depends(c, p) then return
         traceode "but ODESolve cannot determine the arbitrary function!";
      %% MUST be linear in h_x and h_y by construction, so ...
      h_x := coeff(first c, df(h,x));  h_x := -first h_x / second h_x;
      h_y := coeff(second c, df(h,y));  h_y := -first h_y / second h_y;
      h_x := exact1_pde(h_x, h_y, y, x);
      symbolic depend1(y, x, t);
      first_int := sub(h = h_x, p = df(y,x), first_int);
      %% traceode first_int = 0;
      first_int := ODESolve!-FirstOrder(first_int, y, x);
      return
         if first_int then ODESolve!-Simp!-ArbConsts(first_int, y, x)
         else traceode "But ODESolve cannot solve it!"
   end ) where depl!* = depl!*$

algebraic procedure ODESolve!-Exact!-2!-test(f, g, p, y, x);
   if ( (df(f,x,2) + 2p*df(f,x,y) + p^2*df(f,y,2)) -
      (df(g,x,p) + p*df(g,y,p) - df(g,y)) or
         (df(f,x,p) + p*df(f,y,p) + 2df(f,y)) - df(g,p,2) ) then 1$


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

switch odesolve_diff$                   % TEMPORARY?
symbolic(!*odesolve_diff := t)$         % TEMPORARY?

fluid '(!*arbvars)$

algebraic procedure ODESolve!-Diff(ode, y, x);
   %% If the derivative of ode factorizes then try to solve each
   %% factor and return the solutions, otherwise return nil.
   %% This is the inverse of detecting an exact ode!
   if symbolic !*odesolve_diff then     % TEMPORARY?
   begin scalar ode_factors, solns;
      load_package solve;  % to allow overriding !*arbvars := t
      traceode1
         "Trying to factorize derivative of ODE algebraically ...";
      ode_factors := factorize num df(ode, x);
      %% { {factor, multiplicity}, ... }
      if length ode_factors = 1 and second first ode_factors = 1 then return;
      traceode "The derivative of the ODE factorizes algebraically ",
         "and each distinct factor ODE will be solved separately ...";
      solns := {};
      while ode_factors neq {} do
      begin scalar fac, deriv_orders, first!!arbconst, arbconsts,
         !*arbvars;
         fac := first first ode_factors;  ode_factors := rest ode_factors;
         deriv_orders := get_deriv_orders(fac, y);
         %% Check for purely algebraic factor:
         if deriv_orders = {} then return;  % no y -- ignore
         if deriv_orders = {0} then return
            for each s in solve(fac, y) do
               if sub(s, ode) = 0 then solns := (s = 0) . solns;
         first!!arbconst := !!arbconst + 1;
         fac := ODESolve!*2(fac, y, x); % to avoid nasty loops
         if not fac then return solns := ode_factors := {};
         arbconsts :=
            for i := first!!arbconst : !!arbconst collect arbconst i;
         %% ***** THIS WILL WORK ONLY FOR EXPLICIT SOLUTIONS *****
         for each soln in fac do
            for each s in solve(sub(soln, ode), arbconsts) do
               solns := sub(s, soln) . solns
      end;
      if solns neq {} then return solns;
      traceode "... but cannot solve all factor ODEs.";
   end$

algebraic procedure ODESolve!-Alg!-Solve(ode, y, x);
   %% Try to solve algebraically for a single derivative and then
   %% solve each solution ode directly.
   begin scalar deriv, L, R, d, root_odes, solns;
      scalar !*fullroots, !*trigform, !*precise;
      %% symbolic(!*fullroots := t);       % Can be VERY slow!
      traceode1
         "Trying to solve algebraically for a single derivative ...";
      deriv := delete(0, get_deriv_orders(ode, y));
      if length deriv neq 1 then return; % not a single deriv
      %% Now ode is an expression in df(y,x,ord) involving no other
      %% derivatives.  Try to solve it algebraically for the
      %% derivative.
      deriv := df(y, x, first deriv);
      if not( smember(deriv, L:=lcof(ode,deriv)) or
              smember(deriv, R:=reduct(ode,deriv)) ) then
         if (d:=deg(ode,deriv)) = 1 then
            return                      % linear in single deriv
         else
            root_odes :=                % single integer power
               { num(deriv - (-R/L)^(1/d)*newroot_of_unity(d)) }
               %% Expand roots of unity later.
      else <<
         root_odes := solve(ode, deriv);
         if not(length root_odes > 1 or
            first root_multiplicities > 1) then return;
         %% Eventually, replace above 3 lines with this:
         %% root_odes := SolvePM(ode, deriv); % `use `plus_or_minus'
         root_odes := for each ode in root_odes collect
            num if symbolic eqcar(caddr ode, 'root_of) then
               sub(part(rhs ode, 2)=deriv, part(rhs ode, 1))
            else lhs ode - rhs ode
      >>;
      traceode "It can be (partially) solved algebraically ",
         "for the single-order derivative ",
         "and each `root ODE' will be solved separately ...";
      solns := {};
      while root_odes neq {} do
      begin scalar soln;
         if (soln := ODESolve!*2(first root_odes, y, x)) then <<
            solns := append(solns, soln);
            root_odes := rest root_odes
         >> else solns := root_odes := {}
      end;
      if solns neq {} then return solns
   end$


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Utility procedures
% ==================

% Linearity and order tests, which are best kept separate!

%% NB: smember in ODE!-Linearp should probably be depends!

symbolic operator ODE!-Linearp$
symbolic procedure ODE!-Linearp(ode, y);
   %% ODE is assumed to be an expression (not an equation).
   %% Returns t if ODE is linear in y, nil otherwise.
   %% Assumes on exp, mcd.
   ODE!-Lin!-Form!-p(numr simp!* ode, !*a2k y)$

symbolic procedure ODE!-Lin!-Form!-p(sf, y);
   %% A standard (polynomial) form `sf' is linear if each of its terms
   %% is linear:
   domainp sf or
      (ODE!-Lin!-Term!-p(lt sf, y) and ODE!-Lin!-Form!-p(red sf, y))$

symbolic procedure ODE!-Lin!-Term!-p(st, y);
   %% A standard (polynomial) term `st' is linear if either (a) its
   %% leading power is linear and its coefficient is independent of y,
   %% or (b) its leading power is independent of y and its coefficient
   %% is linear:
   begin scalar knl;  knl := tvar st;
      return if knl eq y or (eqcar(knl, 'df) and cadr knl eq y) then
         %% Kernel knl is either y or a derivative of y (df y ...)
         tdeg st eq 1 and not depends(tc st, y)
      else if not depends(knl, y) then ODE!-Lin!-Form!-p(tc st, y)
   end$


symbolic operator ODE!-Order$
symbolic procedure ODE!-Order(u, y);
   %% u is initially an ODE, assumed to be an expression (not an
   %% equation).  Returns its order wrt. y.
   if atom u then 0
   else if car u eq 'df and cadr u eq y then
      %% u = (df y x n) or (df y x)
      (if cdddr u then cadddr u else 1)
   else
      max(ODE!-Order(car u, y), ODE!-Order(cdr u, y))$


symbolic operator get_deriv_orders$
symbolic procedure get_deriv_orders(ode, y);
%    %% Return range of orders of derivatives df(y,x,n) in ode as the
%    %% algebraic list {min_ord, min_d_ord, max_ord} where min_ord
%    %% includes 0, and min_d_ord excludes 0.
   %% Return the SET of all orders of derivatives df(y,x,n) in ode as
   %% an unsorted algebraic list.  Empty if ode freeof y.
   begin scalar result;
      ode := kernels numr simp!* ode;
      if null ode then return makelist nil;
      result := get_deriv_ords_knl(car ode, y);
      for each knl in cdr ode do
         result := union(get_deriv_ords_knl(knl, y), result);
%       return {'list, min_ord,
%          if zerop min_ord then min!* delete(0, result) else min_ord,
%             max!* result} where min_ord = min!* result
      return makelist result
   end$

symbolic procedure get_deriv_ords_knl(knl, y);
   %% Return a list of all orders of derivatives df(y,x,n) in kernel
   %% knl, treating y as df(y,x,0).
   if atom knl then (if knl eq y then (0 . nil))
   else if car knl eq 'df then
      (if cadr knl eq y then
         (if cdddr knl then cadddr knl else 1) . nil)
   else
      ( if in_car then union(in_car, in_cdr) else in_cdr )
         where in_car = get_deriv_ords_knl(car knl, y),
            in_cdr = get_deriv_ords_knl(cdr knl, y)$


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Support for an n'th root of unity operator
% ==========================================

algebraic operator root_of_unity, plus_or_minus$

fluid '(!*intflag!*)$                   % true when in the integrator

% Simplify powers of these operators, but only when not in the
% integrator, which seems to be upset by this:
algebraic let (plus_or_minus(~tag))^2 => 1 when symbolic not !*intflag!*,
   (root_of_unity(~n, ~tag))^n => 1 when symbolic not !*intflag!*$
% Should really be more general, e.g.
%    (root_of_unity(~n, ~tag))^nn => 1 when fixp(nn/n)

algebraic procedure newroot_of_unity(n);
   if n = 0 then RedErr "zeroth roots of unity undefined"
   else if numberp n and (n:=abs num n) = 1 then 1
   else if n = 2 then
      plus_or_minus(newroot_of_unity_tag())
   else
      root_of_unity(n, newroot_of_unity_tag())$

algebraic procedure newplus_or_minus;
   %% Like this for immediate evaluation, especially in symbolic mode:
   symbolic {'plus_or_minus, newroot_of_unity_tag()}$

%% fluid '(!!root_of_unity)$  !!root_of_unity := 0$

algebraic procedure newroot_of_unity_tag;
   %% symbolic mkid('tag_, !!root_of_unity := add1 !!root_of_unity)$
   symbolic mkrootsoftag()$             % defined in module solve/solve1

define expand_plus_or_minus = expand_roots_of_unity$
define expand_root_of_unity = expand_roots_of_unity$

symbolic operator expand_roots_of_unity$
flag('(expand_roots_of_unity), 'noval)$

symbolic procedure expand_roots_of_unity u;
   begin scalar !*NoInt;  !*NoInt := t;  u := aeval u;
      return makelist union(            % discard repeats
         for each uu in (if rlistp u then cdr u else {u}) join
            cdr expand_roots_of_unity1 makelist {uu}, nil)
   end$

symbolic procedure expand_roots_of_unity1 u; % u is an rlist
   ( if r then expand_roots_of_unity1 makelist append(
      (if car r eq 'plus_or_minus then
         cdr subeval{{'equal, r, -1}, u}
      else
      begin scalar n, n!-1;
         if not fixp(n := numr simp!* cadr r) then
            TypErr(n, "root of unity");
         n!-1 := sub1 n;
         return for m := 1 : n!-1 join
            cdr algebraic sub(r = exp(i*2*pi*m/n), u)
      end), cdr subeval{{'equal, r, 1}, u} )
   else u ) where r = find_root_of_unity cdr u$

symbolic procedure find_root_of_unity u; % u is a list
   if atom u then nil
   else if car u eq 'plus_or_minus then u
   else if car u eq 'root_of_unity and evalnumberp cadr u then u
   else find_root_of_unity car u or find_root_of_unity cdr u$

endmodule$

end$