[Maxima-commits] CVS: maxima/share/contrib to_poly_solver.mac, NONE, 1.1

[Barton W. <wil...@us...>]

Update of /cvsroot/maxima/maxima/share/contrib
In directory 23jxhf1.ch3.sourceforge.com:/tmp/cvs-serv18758

Added Files:
	to_poly_solver.mac 
Log Message:
o new file

--- NEW FILE: to_poly_solver.mac ---
/*  Author Barton Willis
    University of Nebraska at Kearney
    Copyright (C) 2008 Barton Willis

    This program is free software: you can redistribute it or modify
    it under the terms of the GNU General Public License version two.

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
    GNU General Public License for details.

    You should have received a copy of the GNU General Public License
    along with this program. If not, see <http://www.gnu.org/licenses/>.
 
*/

eval_when(translate,declare_translated(to_poly_solve,to_poly_solve_h, safer_subst,
   merge_solutions, solution_subst,lambertify, catch_easy_cases, expunge_nonrelavant_cnds));

/* Conditionally load to_poly (version 2 or later), grobner, lrats, opsubst, unwind_protect,
and basic (push macro). */

if get('to_poly,'version) = false then load("topoly.lisp");
if get('to_poly,'version) < 2 then error("You need to update 'to_poly' to version 2 or greater");
if not(?functionp('poly_reduced_grobner)) then load("grobner");
block([inflag : true], if not(member('fullratsubst, map('op, functions))) then load("lrats"));
if not(?functionp('opsubst)) then load("opsubst");
if not(?functionp('unwind_protect)) then load("unwind_protect");
block([inflag : true], if not(member('push, map('op, macros))) then load("basic"));
load("to_poly_solve_extra");

lratsubst(l, e) := (
  if not(listp(l)) then l : [l],
  for lk in l do (
    if equation_p(lk) then e : ratsubst(rhs(lk), lhs(lk), e) else
    error("The first argument to 'lratsubst' must be an equation or a list of equations")),
  e);
 
/* Return a lambda form that is the composition of the functions in the list l.
   For example, compose_functions([f,g]) --> lambda([gxxxx], f(g(gxxxx))). When
   a list member isn't a function, generally funmake will signal an error. */

compose_functions(l) := block([z : new_variable('general), f],
  if listp(l) then (
    l : reverse(l),
    f : z,
    for lk in l do f : funmake(lk, [f]),
    buildq([z,f], lambda([z], f)))
  else error("The argument to 'compose_functions' must be a list."));

/* The function dfloat is similar to float, but dfloat applies rectform if dfloat fails
   to produce a complex double float result. Because float is both a function and an option
   variable (default value is false), it's easier to include dfloat than float in a
   simpfuncs list.*/

dfloat(e) := block([ef : float(e)], if complex_number_p(ef) then ef else float(rectform(e)));

/* Return true if the argument is either a symbol (but not a system constant such as %pi), a
   subscripted variable, or a conjugated variable.*/

variable_p(e) := block([inflag : true],
  (symbolp(e) and not(member('constant, apply('properties, [e])))) or
  subvarp(e) or
  opequal_p(e,'conjugate) and variable_p(first(e)));

/* Apply %and to the list or set l. Signal an error when the argument isn't a list or set.*/

simp_and(l) := if listp(l) or setp(l) then xreduce("%and",l, true)
  else error("The argument to 'simp_and' must be a list");

/* Return true if the operator of e is a member of f; otherwise return false. */

opequal_p(e,[f]) := block([inflag : true], not(mapatom(e)) and member(op(e),f));

/* Return a set of nonconstant variables that appear explicitly in the expression e. The variable x
   and conjugate(x) are distinct variables.*/

set_of_variables(e) := block([inflag : true],
  if variable_p(e) then set(e)
  else if mapatom(e) then set()
  else xreduce('union, (map('set_of_variables, args(e))), set()));

singular_solutions(eq, v, p, f, maxdepth) := block([cnd, free_vars, sol, vs, ps, acc : %union()],
  cnd : last(last(triangularize(jacobian(eq,v)))),
  cnd : first(elim_allbut(cons(cnd,eq),p)),
  sol : to_poly_solve(append(cnd,eq),append(v,p), 'simpfuncs = f, 'maxdepth = maxdepth),
  for sk in sol do (
    vs : map("=", v, subst(sk, v)),
    ps : map("=", p, subst(sk, p)),
    ps : simp_and(ps),
    acc : %union(%if(ps, vs, %union()), acc)),
  sol : algsys(eq,v),
  cnd : simp_and(map(lambda([s], s # 0),cnd)),
  for sk in sol do (
    vs : map("=", v, subst(sk, v)),
    acc : %union(%if(cnd, vs, %union()), acc)),
  acc);

/* Return true if the operator of e is "="; otherwise return false. */
equation_p(e) := block([inflag : true], not(mapatom(e)) and is(op(e) = "=")); 

/* Return true if the operator of e is %if; otherwise return false. */
if_p(e) := block([inflag : true], not(mapatom(e)) and is(op(e) = %if));

/* Return true if the operator of e is %union; otherwise return false. */
%union_p(e) :=  block([inflag : true], not(mapatom(e)) and is(op(e) = '%union));

recip(e) := if listp(e) or equation_p(e) then map('recip, e)
               else if is(equal(e,0))=true then 'infinity else 1/e;

solution_subst(sol, v, f) := (
  if is(sol = %union()) then sol
  else if if_p(sol) then %if(first(sol), solution_subst(second(sol),v,f), solution_subst(third(sol),v,f))
  else v = apply('f, [subst(sol,v)]));

safer_subst(l, e) := block([algebraic : true, r, cnd, error_save : errormsg],
  unwind_protect(
    (errormsg : false,
    if listp(e) or equation_p(e) or %union_p(e) then map(lambda([s], safer_subst(l,s)), e)
    else if if_p(e) then (
      cnd : simp_inequality(safer_subst(l, first(e))),
      if is(cnd = false) then safer_subst(l, third(e))
      else if is(cnd = true) then safer_subst(l, second(e))
      else map(lambda([s], safer_subst(l,s)), e))
    else (
      r : errcatch(subst(l,e)),
      if emptyp(r) then recip(safer_subst(l, ratsimp(1/e))) else first(r))),
   errormsg : error_save));

/* Solve equations that involve both x and conjugate(x). This is not a user-level function:

   e = list of equations to be solved,
   vsave = list of solution variables,
   vc = list of solution variables that appear in the equations as conjugated,
   simpfuncs = same as for to_poly_solve.
*/

conjugate_solver(esave,vsave,vc, funs, maxdepth) := block([v,vv,acc : %union(), enew : [],
           real_flag : false, ec, e, subs, cnd, nv, xr,xi,sol,%c, cv, algebraic : true],
  [e,v,vc] : map('setify,[esave,vsave,vc]),
  for ek in e do (
    ec : conjugate(ek),
    if linearly_dependent_p(rhs(ek)-lhs(ek),rhs(ec)-lhs(ec)) then real_flag : true
    else (push(ec, enew), push(ek, enew))),

  e : enew,
  if is(real_flag = true) then (
    vv : intersection(v, conjugate(vc)),
    subs : set(),
    cnd : set(),
    nv : setdifference(v,vv),
    for vk in vv do (
      xr : new_variable('real),
      nv : adjoin(xr, nv),
      xi : new_variable('real),
      nv : adjoin(xi, nv),
      subs : adjoin(vk = xr + %i * xi, subs)),
    e : setify(lratsubst(listify(subs), esave)),
    e : union(e,map('conjugate,e)),
    sol : to_poly_solve(e, nv,'simpfuncs = funs, 'maxdepth = maxdepth, 'field = 'real),
    acc : %union(),
    subs : listify(subs),
    for sk in sol do (
      if if_p(sk) then (
        cnd : first(sk),
        sk : second(sk))
      else (
        cnd : true),
      cnd : xreduce(lambda([a,b], a %and isreal_p(b)), subst(sk, nv), cnd),
      acc : %union(acc, %if(cnd, subst(sk, subs), %union()))),
    acc)
      
  
  else (
    vv : intersection(v, conjugate(vc)),
    subs : set(),
    cnd : set(),
    nv : set(),
    for vk in vv do (
      %c : new_variable('complex),
      nv : adjoin(%c, nv),
      subs : adjoin(conjugate(vk) = %c, subs)),
    nv : union(v,nv),
    e : lratsubst(listify(subs), e),
    sol : to_poly_solve(e, nv,'simpfuncs = funs, 'maxdepth = maxdepth, 'field = 'complex),
    for sk in sol do (
      vv : map(lambda([s], solution_subst(sk, s, 'identity)), vsave),
      cv : listify(map(lambda([s], solution_subst(sk, s, 'identity)), vc)),
      cnd : map('rootscontract, listify(subst(sk, map('conjugate, subs)))),
      cnd : xreduce("%and", cnd),
      sk : map(lambda([s], solution_subst(sk, s, 'identity)), vsave),
      sk : %if(cnd, sk, %union()),
      acc : %union(acc, sk)),
    acc));


/* Convert log(constant) --> rectform(log(constant)). */
rectform_log_if_constant(e) :=
  subst(lambda([s], if constantp(s) then rectform(log(s)) else log(s)), 'log, e);

/* Convert realpart and imagpart to conjugates. The calls to nounify are required by
noun / verb bugs, I think. */
  
real_imagpart_to_conjugate(e) :=
   subst([nounify('realpart) = lambda([s], (s + conjugate(s))/2),
     nounify('imagpart) = lambda([s], (s - conjugate(s))/(2 * %i))], e);
  
merge_solutions(sol1, sol2) := block([cnd1, cnd2, a1, a2, b1, b2],
  if is(sol1 = %union()) then sol2
  else if is(sol2 = %union()) then sol1
  else (
    [cnd1, a1, b1] : if listp(sol1) then [true, sol1, %union()] else args(sol1),
    [cnd2, a2, b2] : if listp(sol2) then [true, sol2, %union()] else args(sol2),
    %if(simp_and([cnd1,cnd2]), append(a1,a2), merge_solutions(b1,b2))));

non_symbol_subst (e, v) := block([vars, listconstvars : false],
  vars : makelist(lambda([s], if variable_p(s) then s = s else s = new_variable('general))(k), k, v),
  e : fullratsubst(vars, e),
  [e, map('reverse, vars)]);

one_to_one_crunch(e) := block([ee, inflag : true, nv],
  ee : if opequal_p(e,"+") then first(e) = -rest(e) else e,
  if opequal_p(ee, "=") and opequal_p(lhs(ee),"^") and opequal_p(rhs(ee),"^") and
     is(first(lhs(ee)) = first(rhs(ee))) then (
       nv :  new_variable('integer),
       [second(lhs(ee)) = second(rhs(ee)) + 2 * %pi * %i * nv / log(first(lhs(ee))), set(nv)])
     else [e, set()]);

catch_easy_cases(e) := block([inflag : true],
  if opequal_p(e, 'log, 'acos) then catch_easy_cases(first(e)-1)
  elseif opequal_p(e, 'asin, 'atan, 'abs) then catch_easy_cases(first(e))
  else e);

/* Lambert Solver -- replace terms of the form w * exp(w) by a %Cxx variable and
push %Cxx = lambert_w(w) into subs. Return both e and the list of subs.*/

lambertify(e, vars) := block([a,asubs, b, bsubs, s, eqs : set(), exp_args, %g, cf,
  subs : set(), inflag : false],
  if mapatom(e) then [e, set()]
  else if op(ratsimp(e)) = "/" then (
    [a, asubs] : lambertify(ratnumer(e),vars),
    [b, bsubs] : lambertify(ratdenom(e),vars),
    [a/b, union(asubs,bsubs)])
  else (
    e : ratexpand(e),
    exp_args : gatherargs(e, "^"),
    for ak in exp_args do (
      if is(first(ak) = '%e) then (
        ak : second(ak),
        %g : new_variable('complex),
        cf : ratcoef(e, exp(ak)) / ak,
        if constantp(cf) and is(cf # 0) and not(lfreeof(vars, ak)) then (
          subs : adjoin(lambert_w(%g) = ak, subs),
          e : ratsubst(%g, ak * exp(ak), e)))),
    [e, subs]));

to_poly_solve(e,v,[extraargs]) := block([vc, all_vars, sol, solvedomain, parameter_list,
                        alias_solver : false, use_grobner, maxdepth, cntx : false, simpfuncs],
 
  simpfuncs : assoc('simpfuncs, extraargs, []),
  if not(member('rectform_log_if_constant, simpfuncs)) then
    simpfuncs : append(simpfuncs,['real_imagpart_to_conjugate,'rectform_log_if_constant]),
  solvedomain : assoc('field, extraargs, 'complex),
  use_grobner : assoc('use_grobner, extraargs, false),
  parameter_list : assoc('parameters, extraargs, false),
  maxdepth : assoc('maxdepth, extraargs, 5),
  if is(maxdepth < 0) then (
    error("Maximum recursion depth exceeded; unable to solve", funmake('to_poly_solve, [e,v])))
  else maxdepth : maxdepth - 1,
 
  /* Convert e and v into sets */
  [e,v] : map(lambda([s], if listp(s) then setify(s) else if setp(s) then s else set(s)), [e,v]),
 
  /* Check that each member of v is a variable. */
  alias_solver : not every('variable_p, v),

  /* Check that no member of e is an inequality. */
  if some(lambda([s], not(mapatom(s)) and opequal_p(s, "<","<=", "#", ">", ">=")), e) then
  error("The first argument to to_poly_solve must not contain any inequalities"),

  /* Decide if we need to dispatch the conjugate_solver. Convert real/imag part to conjugate.*/

  e : real_imagpart_to_conjugate(e),
  vc : subset(set_of_variables(e), lambda([s], opequal_p(s,'conjugate))),
  v : setdifference(v,vc),

  /* Listify e, v, and vc.*/
  [e,v,vc] : map('listify,[e,v, vc]),
  
  /* Special case dispatcher.*/

  if some(lambda([s], if_p(s)), e) then (
    error("Not yet implemented!"))
  
  else if alias_solver then (
    block([ee, nonvar_sub, sol],
      ee : non_symbol_subst(e,v),
      [ee, nonvar_sub] : ee,
      safer_subst(reverse(nonvar_sub), to_poly_solve(ee, map('lhs, nonvar_sub), 'maxdepth = maxdepth))))
  
  else if every(lambda([s], polynomialp(rhs(s)-lhs(s), v, lambda([w], lfreeof(v,w)))),e) then (
    block([freevars : setify(%rnum_list), simpfuncs_lambda : compose_functions(simpfuncs)],
      if use_grobner then (
        e : map(lambda([s], rhs(s)-lhs(s)), e),
        e : poly_reduced_grobner(e,listofvars(e))),      
      if listp(parameter_list) then (
        singular_solutions(e,v,parameter_list, simpfuncs, maxdepth))
      else (
        sol : algsys(e, v),
        freevars : setdifference(setify(%rnum_list), freevars),
        for vk in freevars do (
          sol : subst(vk = new_variable(solvedomain), sol)),
        sol : map(lambda([s], map("=", map('lhs, s), map(simpfuncs_lambda, map('rhs, s)))), sol),
        xreduce('%union, sol, %union()))))
  
  else if not(emptyp(vc)) then (
    conjugate_solver(e,v,vc,simpfuncs, maxdepth))
  
  else (
    to_poly_solve_h(e,v, simpfuncs, maxdepth, use_grobner)));
           
/* Try to solve the equation(s) e for the the variable(s) v. Both 'e' and 'v' can be
   a single Maxima expression, a list, or a set. Members of 'e' that aren't equations
   are (not of the form a = b) assumed to vanish. Each solution gets simplified by
   the functions in 'simpfuncs.' If simpfuncs contains non-identity transformations,
   well, it's not my fault.

   Users should call to_poly_solve, not to_poly_solve_h.
*/

to_poly_solve_h(e,v, simpfuncs, maxdepth, grob) := block([eqs : set(), nonalg : set(), listconstvars : false,
  %gvars, evars, sol, vars, freevars, eqs_vars, non_free, true_sol, %ivar, esave, %sol_vars, %cnds,
  xsol, nsol, inflag : true, cnds : set(), nonfree, eq, lambert_subs, var_exception, simpfuncs_lambda,
  nsol_vars],

  /* save a list of the variables in e */
  evars : set_of_variables(e),

  /* Convert e and v into sets */
  
  e : if listp(e) then setify(e) else if setp(e) then e else set(e),
  v : if listp(v) then setify(v) else if setp(v) then v else set(v),
  v : listify(v),
 
  /* Convert each member of e to polynomial form. The function to_poly extracts numerators and demands
  that the denominators be nonvanishing.*/
  
  block([acc : []],
    for ek in e do (
      [ek, var_exception] : one_to_one_crunch(ek),
      ek : catch_easy_cases(lhs(ek) - rhs(ek)),
      ek : lhs(ek) - rhs(ek),
      ek :  trigexpand(logcontract(rationalize(ek))),
      ek :  subst(lambda([a,b], if numberp(a) and not(numberp(b)) then exp(log(a)*b) else a^b), "^", ek),
      [ek, lambert_subs] : lambertify(ek, v), 
      nonalg : union(nonalg, lambert_subs),
      ek : errcatch(to_poly(ek,v)),
      if emptyp(ek) then error("unable to solve") else ek : first(ek),
      push(ek, acc)),
    e : acc),

  /* Gather equations, conditions, and non-algebraic equations into lists eqs, cnds, and nonalg,
  respectively.*/
 
  for ek in e do (
    eqs : union(eqs, setify(first(ek))),
    cnds : union(cnds, setify(second(ek))),
    nonalg : union(nonalg, setify(third(ek)))),
   
   /* convert cnds into a boolean expression */
   cnds : simp_and(cnds),
  
   /* nonfree is the union of v (the set of variables to solve for) and the
      set of variables introduced by to_poly. The set nonfree is the new set
      of variables that we need to solve for.
   */
  
   nonfree : union(setify(v), setdifference(set_of_variables(eqs), evars)),
   nonfree : setdifference(nonfree, var_exception),
   esave : union(set_of_variables(eqs), nonfree),
   sol : to_poly_solve(eqs, nonfree, 'simpfuncs = simpfuncs, 'maxdepth = maxdepth, 'use_grobner = grob),

   /* The solve process introduces new variables. Let %sol_vars be the set 
   of variables introduced by to_poly_solve */
   
   %sol_vars : setdifference(set_of_variables(sol), union(esave, var_exception)),

   /* Update the set of variables we need to solve for */
   nonfree : union(setify(v), %sol_vars),
     
   /* Each member of nonalg looks like %g = exp(blob), %g = log(blob), or
      inverse_lambert_w(blob) = atom. Convert
        %g = exp(blob) --> log(%g) + 2 * %i * %pi * %ivar = blob,
        %g = log(blob) --> exp(%g) = blob.
      where %ivar eventually is replaced by an integer gentemp. */

     nonalg : map(lambda([s],
       if opequal_p(rhs(s),"^") then (
         log(lhs(s)) + 2 * %i * %pi * '%ivar - second(rhs(s)))
       elseif opequal_p(rhs(s),'log) then (
         exp(lhs(s)) - first(args(rhs(s))))
       else s),
     nonalg),

   simpfuncs_lambda : compose_functions(simpfuncs),
   
   /* Paste each solution into nonalg. The substitution can yield log(0), so use errcatch. */
   true_sol : set(),
   for sk in sol do (
     sk : real_imagpart_to_conjugate(sk),
     /* bug when sk isn't a list! */

      if emptyp(nonalg) then (
        eq : set())
      else (
        eq : errcatch(subst(sk, subst('%ivar = new_variable('integer),nonalg))),
        eq : if emptyp(eq) then [1=0] else first(eq)),

      xsol : to_poly_solve(eq, intersection(nonfree, set_of_variables(eq)),
        'simpfuncs = simpfuncs, 'maxdepth = maxdepth, 'use_grobner = grob),
      
      for xk in xsol do (
        xk : real_imagpart_to_conjugate(xk),
        xk : merge_solutions(sk,xk),
        nsol : map(lambda([s], solution_subst(xk, s, simpfuncs_lambda)), v),
        if listp(nsol) and if_p(first(nsol)) then nsol : first(nsol),
        if if_p(nsol) and not(listp(second(nsol))) then nsol : %if(first(nsol),
          [second(nsol)], third(nsol)),
        nsol_vars : set_of_variables(nsol),
        %cnds : safer_subst(sk,cnds),
        %cnds : expunge_nonrelavant_cnds(%cnds, nsol),
        nsol : %if(%cnds, nsol, %union()),
        true_sol : adjoin(nsol, true_sol))),
    xreduce('%union, true_sol, %union()));

/* Maybe all this should be fixed higher up; for now ...*/

expunge_nonrelavant_cnds (cnds, sol) := block([listconstvars : false, sol_vars : listofvars(sol)],
  if is(cnds = true) or is (cnds = false) then cnds
  else if opequal_p(cnds, "%and", "%or") then (
    map(lambda([x], if lfreeof(sol_vars,x) then true else x), cnds))
  else if lfreeof(sol_vars, cnds) then true else cnds);