share/to_poly_solve/to_poly_solve.mac

   1 /*  Author Barton Willis
   2     University of Nebraska at Kearney
   3     Copyright (C) 2008, 2009, 2010, 2013 Barton Willis
   4
   5     This program is free software: you can redistribute it or modify
   6     it under the terms of the GNU General Public License version two.
   7
   8     This program is distributed in the hope that it will be useful,
   9     but WITHOUT ANY WARRANTY; without even the implied warranty of
  10     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
  11     GNU General Public License for details.
  12
  13     You should have received a copy of the GNU General Public License
  14     along with this program. If not, see <http://www.gnu.org/licenses/>. */
  15
  16 eval_when(translate,declare_translated(to_poly_solve,to_poly_solve_h, safer_subst,
  17    merge_solutions, solution_subst,lambertify, catch_easy_cases, expunge_nonrelevant_cnds));
  18
  19 /* Conditionally load to_poly (version 2 or later), grobner, lrats, opsubst, unwind_protect,
  20 and basic (push macro). */
  21
  22 if get('to_poly,'version) = false then load("to_poly");
  23 if get('to_poly,'version) < 2 then error("You need to update 'to_poly' to version 2 or greater");
  24 if not(?functionp('poly_reduced_grobner)) then load("grobner");
  25 block([inflag : true], if not(member('fullratsubst, map('op, functions))) then load("lrats"));
  26 block([inflag : true], if not(member('lratsubst, map('op, functions))) then load("lrats"));
  27 if not(?functionp('opsubst)) then load("opsubst");
  28 block([inflag : true], if not(member('push, map('op, macros))) then load("basic"));
  29 if not(?functionp('simpfuncall)) then load("simplifying");
  30 load("to_poly_solve_extra");
  31
  32 /* Return a lambda form that is the composition of the functions in the list l.
  33    For example, compose_functions([f,g]) --> lambda([gxxxx], f(g(gxxxx))). When
  34    a list member isn't a function, generally funmake will signal an error. */
  35
  36 compose_functions(l) := block([z : new_variable('general), f],
  37   if listp(l) then (
  38     l : reverse(l),
  39     f : z,
  40     for lk in l do f : funmake(lk, [f]),
  41     buildq([z,f], lambda([z], f)))
  42   else error("The argument to 'compose_functions' must be a list."));
  43
  44 /* The function dfloat is similar to float, but dfloat applies rectform if dfloat fails
  45    to produce a complex double float result. Because float is both a function and an option
  46    variable (default value is false), it's easier to include dfloat than float in a
  47    simpfuncs list.*/
  48
  49 dfloat(e) := block([ef : float(e)], if complex_number_p(ef) then ef else float(rectform(e)));
  50
  51 /* Return true if the argument is either a symbol (but not a system constant such as %pi), a
  52    subscripted variable, or a conjugated variable.*/
  53
  54 variable_p(e) := block([inflag : true],
  55   (symbolp(e) and not(featurep(e,'constant))) or
  56   subvarp(e) or
  57   opequal_p(e,'conjugate) and variable_p(first(e)));
  58
  59 /* Apply %and to the list or set l. Signal an error when the argument isn't a list or set.*/
  60
  61 simp_and(l) := if listp(l) or setp(l) then xreduce("%and",l, true)
  62   else error("The argument to 'simp_and' must be a list");
  63
  64 /* Return true if the operator of e is a member of f; otherwise return false. */
  65
  66 opequal_p(e,[f]) := block([inflag : true], not(mapatom(e)) and member(op(e),f));
  67
  68 /* Return a set of nonconstant variables that appear explicitly in the expression e. The variable x
  69    and conjugate(x) are distinct variables.*/
  70
  71 set_of_variables(e) := block([inflag : true],
  72   if variable_p(e) then set(e)
  73   else if mapatom(e) then set()
  74   else xreduce('union, (map('set_of_variables, args(e))), set()));
  75
  76 singular_solutions(eq, v, p, f, maxdepth) := block([cnd, free_vars, sol, vs, ps, acc : %union()],
  77   cnd : last(last(triangularize(jacobian(eq,v)))),
  78   cnd : first(elim_allbut(cons(cnd,eq),p)),
  79   sol : to_poly_solve(append(cnd,eq),append(v,p), 'simpfuncs = f, 'maxdepth = maxdepth),
  80   for sk in sol do (
  81     vs : map("=", v, subst(sk, v)),
  82     ps : map("=", p, subst(sk, p)),
  83     ps : simp_and(ps),
  84     acc : %union(%if(ps, vs, %union()), acc)),
  85   sol : algsys(eq,v),
  86   cnd : simp_and(map(lambda([s], s # 0),cnd)),
  87   for sk in sol do (
  88     vs : map("=", v, subst(sk, v)),
  89     acc : %union(%if(cnd, vs, %union()), acc)),
  90   acc);
  91
  92 /* Return true if the operator of e is "="; otherwise return false. */
  93 equation_p(e) := block([inflag : true], not(mapatom(e)) and is(op(e) = "="));
  94
  95 /* Return true if the operator of e is %if; otherwise return false. */
  96 if_p(e) := block([inflag : true], not(mapatom(e)) and is(op(e) = '%if));
  97
  98 /* Return true if the operator of e is %union; otherwise return false. */
  99 %union_p(e) :=  block([inflag : true], not(mapatom(e)) and is(op(e) = '%union));
 100
 101 solution_subst(sol, v, f) := (
 102   if is(sol = %union()) then sol
 103   else if if_p(sol) then %if(first(sol), solution_subst(second(sol),v,f), solution_subst(third(sol),v,f))
 104   else v = apply(f, [subst(sol,v)]));
 105
 106 /* Substitute l into e. When l or e is a conditional, do the substitution in a way that
 107    doesn't unnecessarily signal an error. Also, for substitutions such as safer_subst([x=y,y=z],x)
 108    keep doing the substitution until the output stops changing. This can lead to an infinite loop,
 109    but the solver should never give safer_subst such an infinite loop... */
 110
 111 safer_subst(l, e) := block([algebraic : true, r, cnd, e1, e2],
 112  if listp(e) or equation_p(e) or %union_p(e) then map(lambda([s], safer_subst(l,s)), e)
 113  else if (e = %union()) or (e = true) or (e = false) or (l = %union()) then e
 114  else if if_p(l) then (
 115    cnd : simp_inequality(first(l)),
 116    e1 : safer_subst(second(l),e),
 117    e2 : safer_subst(third(l),e),
 118    %union(%if(expunge_nonrelevant_cnds(cnd, e1), e1, %union()),
 119           %if(expunge_nonrelevant_cnds(cnd, e2), e2, %union())))
 120   else if if_p(e) then (
 121    /* When e is a conditional, substitute and simplify the predicate before substituting into the rest of e.*/
 122    cnd : simp_inequality(safer_subst(l, first(e))),
 123    if is(cnd = false) then safer_subst(l, third(e))
 124    else if is(cnd = true) then safer_subst(l, second(e))
 125    else %if(cnd, safer_subst(l, second(e)), safer_subst(l, third(e))))
 126  else (
 127    r : errcatch(subst(l,e)),
 128    while (r # [ ] and e # first(r)) do ( /* iterate to a fixed point */
 129      e : first(r),
 130      r : errcatch(subst(l,e))),
 131    if emptyp(r) then error("Unable to substitute") else first(r)));
 132
 133 list_conjugate_vars(e) := block([ ],
 134   listify(setify(sublist(xreduce('append, gatherargs(e, 'conjugate)), lambda([s], mapatom(s))))));
 135
 136 /* Let e be the output of to_poly_solve that involves the solution for some internal variables
 137    and let v be the list of user variables. The function expunge_nonrelevant_solutions deletes
 138    the solutions for the internal variables. */
 139
 140 expunge_nonrelevant_solutions(e,v) :=
 141   if mapatom(e) then e
 142   else if op(e) = '%union or op(e) = '%if then (
 143     map(lambda([s], expunge_nonrelevant_solutions(s,v)), e))
 144   else if listp(e) then (
 145     sublist(e, lambda([s], member(lhs(s),v))))
 146   else e;
 147
 148 /* Solve equations that involve both x and conjugate(x). This is not a user-level function.
 149      e = list of equations to be solved,
 150      v = list of solution variables,
 151      vc = list of solution variables that appear in the equations as conjugated.
 152
 153   The input vc is unused. */
 154
 155 conjugate_solver(esave,v,vc,funs,maxdepth) := block([e, ee, sub, nv, rv, sol, cnd, acc, sx1, sx2, cnd1, cnd2],
 156   e : setify(append(esave,conjugate(esave))),
 157   vc : list_conjugate_vars(e),
 158   sub : map(lambda([s], conjugate(s) = new_variable('complex)), vc),
 159   nv : map('rhs, sub),
 160   e : listify(subst(sub, e)),
 161   nv : append(v,nv),
 162   sol : to_poly_solve(e, nv, 'maxdepth = maxdepth, 'simpfuncs = funs, 'field = 'complex),
 163   if not mapatom(sol) and op(sol) = '%union then (
 164     sub : map(lambda([s],  rhs(s)-lhs(s)=0),sub),
 165     acc : %union(),
 166     for sx in sol do (
 167       if if_p(sx) then (
 168         cnd : first(sx),
 169         sx1 : second(sx),
 170         sx2 : third(sx),
 171         cnd1 : safer_subst(sx1,sub),
 172         cnd2 : safer_subst(sx2,sub),
 173         cnd1 : simp_and(cnd1),
 174         cnd2 : simp_and(cnd2),
 175         acc : %union(acc, %if(cnd,
 176             %if(cnd1, expunge_nonrelevant_solutions(sx1,v),%union()),
 177             %if(cnd2, expunge_nonrelevant_solutions(sx2,v),%union()))))
 178       else (
 179         cnd : subst(sx,sub),
 180         cnd : simp_and(cnd),
 181         acc : %union(acc, (%if(cnd, expunge_nonrelevant_solutions(sx,v),%union()))))),
 182     acc)
 183   else simpfunmake('%solve,[esave,v, 'maxdepth = maxdepth, 'simpfuncs = funs]));
 184
 185
 186 /* Convert log(constant) --> rectform(log(constant)). */
 187 rectform_log_if_constant(e) :=
 188   subst(lambda([s], if constantp(s) then rectform(log(s)) else log(s)), 'log, e);
 189
 190 /* Convert realpart and imagpart to conjugates. The calls to nounify are required by
 191 noun / verb bugs, I think. */
 192
 193 real_imagpart_to_conjugate(e) :=
 194    subst([nounify('realpart) = lambda([s], (s + conjugate(s))/2),
 195      nounify('imagpart) = lambda([s], (s - conjugate(s))/(2 * %i))], e);
 196
 197 merge_solutions(sol1, sol2) := block([cnd1, cnd2, a1, a2, b1, b2],
 198   if is(sol1 = %union()) then sol2
 199   else if is(sol2 = %union()) then sol1
 200   else (
 201     [cnd1, a1, b1] : if listp(sol1) then [true, sol1, %union()] else args(sol1),
 202     [cnd2, a2, b2] : if listp(sol2) then [true, sol2, %union()] else args(sol2),
 203     %if(simp_and([cnd1,cnd2]), append(a1,a2), merge_solutions(b1,b2))));
 204
 205 non_symbol_subst (e, v) := block([vars, listconstvars : false],
 206   vars : makelist(lambda([s], if variable_p(s) then s = s else s = new_variable('general))(k), k, v),
 207   e : fullratsubst(vars, e),
 208   [e, map('reverse, vars)]);
 209
 210 (one_to_one_reduce[a,b] := false,
 211   one_to_one_reduce['sin,'sin] : lambda([a,b],
 212       (a - b - 2 * %pi * new_variable('integer)) * (a + b - 2 * %pi * new_variable('integer) + %pi)),
 213
 214   one_to_one_reduce['cos,'cos] : lambda([a,b],
 215       (a - b - 2 * %pi * new_variable('integer)) * (a + b - 2 * %pi * new_variable('integer))),
 216
 217   one_to_one_reduce['sin,'cos] : lambda([a,b],
 218     (a/2 - b/2 - %pi/4 - %pi * new_variable('integer))*(a/2 + b/2 - %pi/4 - %pi * new_variable('integer))),
 219
 220   one_to_one_reduce['cos,'sin] :  lambda([a,b],
 221     (b/2 - a/2 - %pi * new_variable('integer) - %pi/4)*(b/2 + a/2 - %pi * new_variable('integer) - %pi/4)),
 222
 223   one_to_one_reduce['tan,'tan] : lambda([a,b], (a-b - %pi*new_variable('integer))/(cos(a)*cos(b))),
 224
 225   one_to_one_reduce['%exp,'%exp] : lambda([a,b],
 226       a - b - 2 * %pi * %i * new_variable('integer)));
 227
 228 one_to_one_crunch(e) := block([inflag : true, f, nv : set(), %exp, q, opl, opr],
 229   e : subst("^" = lambda([a,b], if a = %e then funmake('%exp, [b]) else a^b),e),
 230   e : if opequal_p(e,"+") then first(e) = -rest(e) else e,
 231   if opequal_p(e, "=") then (
 232     e : block([inflag : false], if not mapatom(lhs(e)) and op(lhs(e)) = "-" then -e else e),
 233     opl : if mapatom(lhs(e)) then false else op(lhs(e)),
 234     opr : if mapatom(rhs(e)) then false else op(rhs(e)),
 235     f : one_to_one_reduce[opl, opr],
 236     if f # false then (
 237       nv : setify(listofvars(e)),
 238       q : one_to_one_crunch(apply(f, [first(lhs(e)),first(rhs(e))])),
 239       e : first(q),
 240       nv : union(second(q), setdifference(setify(listofvars(e)), nv)))
 241     else (
 242       e : lhs(e) - rhs(e))),
 243   e : subst(%exp = lambda([a], %e^a),e),
 244   [e, nv]);
 245
 246
 247 (function_inverse[otherwise]:= false,
 248   function_inverse['acos] : lambda([s], cos(s)),
 249   function_inverse['acosh] : lambda([s], cosh(s)),
 250   function_inverse['asin] : lambda([s], sin(s)),
 251   function_inverse['asinh] : lambda([s], sinh(s)),
 252   function_inverse['atan] : lambda([s], tan(s)),
 253   function_inverse['atanh] : lambda([s], tanh(s)),
 254   function_inverse['log] : lambda([s], exp(s)));
 255
 256 (function_zeros[otherwise] := lambda([], false),
 257   function_zeros['abs] : lambda([],0),
 258   function_zeros['conjugate] : lambda([],0),
 259   function_zeros['sin] : lambda([], %pi*new_variable('integer)),
 260   function_zeros['cos] : lambda([], %pi*(2*new_variable('integer)+1)/2),
 261   function_zeros['log] : lambda([],1),
 262   function_zeros['acos] : lambda([],1),
 263   function_zeros['acosh] : lambda([],1),
 264   function_zeros['asin] : lambda([],0),
 265   function_zeros['asinh] : lambda([],0),
 266   function_zeros['atan] : lambda([], 0),
 267   function_zeros['atanh] : lambda([],0),
 268   function_zeros['acsc] : lambda([], []),
 269   function_zeros['acsch] : lambda([], []),
 270   function_zeros['acoth] : lambda([],[]),
 271   function_zeros['acot] : lambda([],[]));
 272
 273 (function_domain[otherwise] := lambda([s],set()),
 274   function_domain['log] : lambda([s], set(s # 0)));
 275
 276 top_level_ops(e) := block([inflag : false],
 277   if mapatom(e) then [ ]
 278   else if opequal_p(e,"+","*") then (
 279     xreduce('append, map('top_level_ops, args(e))))
 280   else [op(e)]);
 281
 282 catch_easy_cases (e, vars) := block([inflag : true, ops, finv, fzero, p, q, a, b,pk,pk_save,
 283   ok : true, cnds : set(), nv : set(), listconstvars : false],
 284   if opequal_p(e,"*") then (
 285     e : map(lambda([s], catch_easy_cases(s,vars)), args(e)),
 286     [xreduce("*", map('first, e)), xreduce('union,map('second,e)), xreduce('union,map('third,e))])
 287   else (
 288     ops : top_level_ops(e),
 289     for f in ops do (
 290       finv : function_inverse[f],
 291       fzero : function_zeros[f],
 292       fzero : apply(fzero,[]),
 293       nv : union(nv, setify(listofvars(fzero))),
 294       if (finv # false) or (fzero # false) then (
 295         p : gatherargs(e,f),
 296         p : sublist(p, lambda([s], not lfreeof(vars,s))),
 297         for pk in p while ok do (
 298           pk_save : pk,
 299           pk : funmake(f, pk),
 300           a : ratcoef(e, pk),
 301           b : e - a * pk,
 302           if is(equal(b,0)) = true and lfreeof(vars,a) and fzero # false then (
 303             ops : [],
 304             ok : false,
 305             if fzero = [] then [1, set(a # 0)]
 306             else (
 307               q : catch_easy_cases(first(pk)-fzero, vars),
 308               e : first(q),
 309               nv : union(nv, third(q)),
 310               cnds : union(cnds, second(q), apply(function_domain[f],pk_save)),
 311               cnds : adjoin(a # 0, cnds)))
 312           else if lfreeof(vars, a) and lfreeof(vars, b) and is(finv # false) then (
 313             ops : [],
 314             ok : false,
 315             q : catch_easy_cases(first(args(pk))-apply(finv, [-b/a]),vars),
 316             nv : union(nv, setdifference(setify(listofvars(q)), setify(listofvars(pk)))),
 317             e : first(q),
 318             cnds : adjoin(a # 0, second(q)))))),
 319     [e, cnds, nv]));
 320
 321 standardize_inverse_trig(e) :=
 322  e : subst(['acot = lambda([s], atan(1/s)),
 323             'acsc = lambda([s], asin(1/s)),
 324             'asec = lambda([s], acos(1/s)),
 325             'acoth = lambda([s], atanh(1/s)),
 326             'acsch = lambda([s], asinh(1/s)),
 327             'asech = lambda([s], acosh(1/s))], e);
 328
 329 /* Lambert Solver -- replace terms of the form w * exp(w) by a %Cxx variable and
 330 push %Cxx = lambert_w(w) into subs. Return both e and the list of subs.*/
 331
 332 lambertify(e,vars) := block([ee, eee, expargs : set(), cexpargs, nexpargs,
 333   %g : new_variable('complex), logargs : set(),
 334   unmatched : true, subs : set()],
 335
 336   ee : expand(lhs(e)-rhs(e),0,0),
 337   expargs : setify(gatherargs(ee, "^")),
 338   cexpargs : subset(expargs, lambda([s], lfreeof(vars,first(s)))),
 339   /* nexpargs : subset(expargs, lambda([s], is(equal(first(s),second(s))))), */
 340   nexpargs : subset(expargs, lambda([s], is(first(s) = second(s)))),
 341   logargs : setify(gatherargs(ee,'log)),
 342
 343   /* look for matches x * c^x = %g, log(c) * x = lambert_w(%g), where c is constant. */
 344
 345   for qk in cexpargs while unmatched do (
 346     eee : ratsubst(%g / second(qk), first(qk)^second(qk), ee),
 347     eee : ratnumer(eee),
 348     if eee # %g and polynomialp(eee,[%g], lambda([s], lfreeof(vars,s))) then (
 349       unmatched : false,
 350       e : eee,
 351       subs : set(lambert_w(%g) = log(first(qk)) * second(qk)))),
 352
 353   /* look for matches x^x = %g, x = log(%g) / lambert_w(log(%g)) */
 354   for qk in nexpargs while unmatched do (
 355     eee : ratsubst(%g, first(qk)^second(qk), ee),
 356     if polynomialp(eee,[%g], lambda([s], lfreeof(vars,s))) then (
 357       unmatched : false,
 358       e : eee,
 359       subs : set(log(%g) / lambert_w(log(%g)) = first(qk)))),
 360
 361   /* look for matches x * c^(-x) = %g, log(c) * x = lambert_w(-%g), where c is a constant. */
 362   for qk in cexpargs while unmatched do (
 363     eee : ratsubst(second(qk) / %g, first(qk)^second(qk), ee),
 364     eee : ratnumer(eee / second(qk)),
 365     if eee # %g and polynomialp(eee,[%g], lambda([s], lfreeof(vars,s))) then (
 366       unmatched : false,
 367       e : eee,
 368       subs : set(-lambert_w(-%g) = log(first(qk)) * second(qk)))),
 369
 370   /* look for matches x + log(x) = %g, x = lambert_w(exp(%g)) */
 371   for qk in logargs while unmatched do (
 372     qk : first(qk),
 373     eee : ratsubst(%g, log(qk) + qk, ee),
 374     eee : ratnumer(eee),
 375     if eee # %g and polynomialp(eee,[%g], lambda([s], lfreeof(vars,s))) then (
 376       unmatched : false,
 377       e : eee,
 378       subs : set(lambert_w(exp(%g)) = qk))),
 379
 380   /* look for matches x * log(x) = %g, x = %g / lambert_w(%g) */
 381   for qk in logargs while unmatched do (
 382     qk : first(qk),
 383     eee : ratsubst(%g, qk * log(qk), ee),
 384     eee : ratnumer(eee),
 385     if polynomialp(eee,[%g], lambda([s], lfreeof(vars,s))) then (
 386       unmatched : false,
 387       e : eee,
 388       subs : set(%g / lambert_w(%g) = qk))),
 389
 390   [e, subs]);
 391
 392 alias(%solve, to_poly_solve);
 393
 394 simp_%solve(e,v,[extraargs]) := block([vc, all_vars, sol, solvedomain, parameter_list,
 395                 algexact : assoc('algexact,extraargs,true),
 396                 alias_solver : false, use_grobner, maxdepth, cntx : false, simpfuncs, ratprint : false],
 397   simpfuncs : assoc('simpfuncs, extraargs, []),
 398   if not(member('rectform_log_if_constant, simpfuncs)) then
 399     simpfuncs : append(simpfuncs,['real_imagpart_to_conjugate,'rectform_log_if_constant]),
 400   solvedomain : assoc('field, extraargs, 'complex),
 401   use_grobner : assoc('use_grobner, extraargs, false),
 402   parameter_list : assoc('parameters, extraargs, false),
 403   maxdepth : assoc('maxdepth, extraargs, 5),
 404   if is(maxdepth < 0) then (
 405     error("Unable to solve"))
 406   else maxdepth : maxdepth - 1,
 407
 408   /* Convert e and v into sets */
 409   [e,v] : map(lambda([s], if listp(s) then setify(s) else if setp(s) then s else set(s)), [e,v]),
 410
 411   /* Check that each member of v is a variable. */
 412   alias_solver : not every('variable_p, v),
 413
 414   /* Decide if we need to dispatch the conjugate_solver. Convert real/imag part to conjugate.*/
 415
 416   e : real_imagpart_to_conjugate(e),
 417   vc : subset(set_of_variables(e), lambda([s], opequal_p(s,'conjugate))),
 418   v : setdifference(v,vc),
 419
 420   /* Listify e, v, and vc.*/
 421   [e,v,vc] : map('listify,[e,v, vc]),
 422
 423   /* Special case dispatcher.*/
 424
 425   if some(lambda([s], if_p(s)), e) then (
 426     error("Not yet implemented!"))
 427
 428   else if some(lambda([s], not(mapatom(s)) and opequal_p(s, "<","<=", "#", ">", ">=")), e) then (
 429     sol : fourier_elim(e, v),
 430     if opequal_p(sol, "or") then subst("or" = '%union, sol)
 431     else if sol = 'emptyset then %union()
 432     else if sol = 'universalset then (
 433       %union(xreduce('append, makelist([s = new_variable('real)], s, v))))
 434     else %union(sol))
 435
 436   else if alias_solver then (
 437     block([ee, nonvar_sub, sol],
 438       ee : non_symbol_subst(e,v),
 439       [ee, nonvar_sub] : ee,
 440       sol : apply('%solve, append([ee,  map('lhs, nonvar_sub)], extraargs)),
 441       if opequal_p(sol,'%solve) then sol else
 442       safer_subst(reverse(nonvar_sub), %solve(ee, map('lhs, nonvar_sub), 'maxdepth = maxdepth))))
 443
 444   else if every(lambda([s], polynomialp(rhs(s)-lhs(s), v, lambda([w], lfreeof(v,w)))),e) then (
 445     block([freevars : setify(%rnum_list), simpfuncs_lambda : compose_functions(simpfuncs)],
 446       if use_grobner then (
 447         e : map(lambda([s], rhs(s)-lhs(s)), e),
 448         e : poly_reduced_grobner(e,listofvars(e))),
 449       if listp(parameter_list) then (
 450         singular_solutions(e,v,parameter_list, simpfuncs, maxdepth))
 451       else (
 452         sol : errcatch(algsys(e, v)),
 453         if emptyp(sol) then (
 454           simpfunmake('%solve, cons(e,cons(v,extraargs))))
 455         else (
 456           sol : first(sol),
 457           freevars : setdifference(setify(%rnum_list), freevars),
 458           for vk in freevars do (
 459             sol : subst(vk = new_variable(solvedomain), sol)),
 460           sol : map(lambda([s], map("=", map('lhs, s), map(simpfuncs_lambda, map('rhs, s)))), sol),
 461           xreduce('%union, sol, %union())))))
 462
 463   else if not(emptyp(vc)) then (
 464     conjugate_solver(e,v,vc,simpfuncs, maxdepth))
 465
 466   else (
 467     sol : errcatch(to_poly_solve_h(e,v, simpfuncs, maxdepth, use_grobner,algexact)),
 468     if emptyp(sol) then simpfunmake('%solve, cons(e,(cons(v,extraargs)))) else first(sol)));
 469
 470 simplifying('%solve, 'simp_%solve);
 471
 472 /* Try to solve the equation(s) e for the the variable(s) v. Both 'e' and 'v' can be
 473    a single Maxima expression, a list, or a set. Members of 'e' that aren't equations
 474    are (not of the form a = b) assumed to vanish. Each solution gets simplified by
 475    the functions in 'simpfuncs.' If simpfuncs contains non-identity transformations,
 476    well, it's not my fault.
 477
 478    Users should call to_poly_solve, not to_poly_solve_h.
 479 */
 480
 481 to_poly_solve_h(e,v, simpfuncs, maxdepth, grob,use_algexact) := block([eqs : set(), nonalg : set(), listconstvars : false,
 482   %gvars, evars, sol, vars, freevars, eqs_vars, non_free, true_sol, %ivar, esave, %sol_vars, %cnds,
 483   xsol, nsol, inflag : true, cnds : set(), nonfree, eq, lambert_subs, var_exception, simpfuncs_lambda,
 484   nsol_vars],
 485
 486   /* save a list of the variables in e */
 487   evars : set_of_variables(e),
 488
 489   /* Convert e and v into sets */
 490
 491   e : if listp(e) then setify(e) else if setp(e) then e else set(e),
 492   v : if listp(v) then setify(v) else if setp(v) then v else set(v),
 493   v : listify(v),
 494
 495   /* Convert each member of e to polynomial form. The function to_poly extracts numerators and demands
 496   that the denominators be nonvanishing.*/
 497
 498   block([acc : []],
 499     for ek in e do (
 500       [ek, var_exception] : one_to_one_crunch(ek),
 501       ek : standardize_inverse_trig(ek),
 502       ek : catch_easy_cases(lhs(ek) - rhs(ek),v),
 503       var_exception : union(third(ek), var_exception),
 504       cnds : union(cnds, second(ek)),
 505       ek : first(ek),
 506       ek : lhs(ek) - rhs(ek),
 507       ek :  trigexpand(logcontract(rationalize(ek))),
 508       ek :  subst(lambda([a,b], if numberp(a) and not(numberp(b)) then exp(log(a)*b) else a^b), "^", ek),
 509       [ek, lambert_subs] : lambertify(ek, v),
 510       nonalg : union(nonalg, lambert_subs),
 511       ek : errcatch(to_poly(ek,v)),
 512       if emptyp(ek) then error("unable to solve") else ek : first(ek),
 513       push(ek, acc)),
 514     e : acc),
 515
 516   /* Gather equations, conditions, and non-algebraic equations into lists eqs, cnds, and nonalg,
 517   respectively.*/
 518
 519   for ek in e do (
 520     eqs : union(eqs, setify(first(ek))),
 521     cnds : union(cnds, setify(second(ek))),
 522     nonalg : union(nonalg, setify(third(ek)))),
 523
 524    /* convert cnds into a boolean expression */
 525    cnds : simp_and(cnds),
 526
 527    /* nonfree is the union of v (the set of variables to solve for) and the
 528       set of variables introduced by to_poly. The set nonfree is the new set
 529       of variables that we need to solve for.
 530    */
 531
 532    nonfree : union(setify(v), setdifference(set_of_variables(eqs), evars)),
 533    nonfree : setdifference(nonfree, var_exception),
 534    nonfree : subset(nonfree, lambda([s], not featurep(s,'integer))),
 535    esave : union(set_of_variables(eqs), nonfree),
 536    sol : to_poly_solve(eqs, nonfree, 'simpfuncs = simpfuncs, 'maxdepth = maxdepth-1, 'use_grobner = grob, 'algexact=use_algexact),
 537    if opequal_p(sol, '%solve) then error("Unable to solve"),
 538
 539    /* The solve process introduces new variables. Let %sol_vars be the set
 540    of variables introduced by to_poly_solve */
 541
 542    %sol_vars : setdifference(set_of_variables(sol), union(esave, var_exception)),
 543
 544    /* Update the set of variables we need to solve for */
 545    nonfree : union(setify(v), %sol_vars),
 546
 547    /* Each member of nonalg looks like %g = exp(blob), %g = log(blob), %g = variable^blob, or
 548       inverse_lambert_w(blob) = atom. Convert
 549         %g = exp(blob) --> log(%g) + 2 * %i * %pi * %ivar = blob,
 550         %g = x^blob --> x = %g^(1/blob),
 551         %g = log(blob) --> exp(%g) = blob.
 552       where %ivar eventually is replaced by an integer gentemp. */
 553
 554      nonalg : map(lambda([s],
 555          if opequal_p(rhs(s),"^")  and first(rhs(s)) = '%e then (
 556            log(lhs(s)) + 2 * %i * %pi * '%ivar - second(rhs(s)))
 557         else if opequal_p(rhs(s),"^") and lfreeof(v, second(rhs(s))) then (
 558            first(rhs(s)) = lhs(s)^(1/second(rhs(s))) * exp(2 * '%pi * '%i * '%ivar / second(rhs(s))))
 559          else if opequal_p(rhs(s),'log) then (
 560            exp(lhs(s)) - first(args(rhs(s))))
 561          else s),
 562        nonalg),
 563
 564    simpfuncs_lambda : compose_functions(simpfuncs),
 565
 566    /* Paste each solution into nonalg. The substitution can yield log(0), so use errcatch. */
 567    true_sol : set(),
 568    for sk in sol do (
 569      sk : real_imagpart_to_conjugate(sk),
 570      /* bug when sk isn't a list! */
 571
 572       if emptyp(nonalg) then (
 573         eq : set())
 574       else (
 575         eq : errcatch(subst(sk, subst('%ivar = new_variable('integer),nonalg))),
 576         eq : if emptyp(eq) then [1=0] else first(eq)),
 577
 578       xsol : to_poly_solve(eq, intersection(nonfree, set_of_variables(eq)),
 579         'simpfuncs = simpfuncs, 'maxdepth = maxdepth-1, 'use_grobner = grob),
 580
 581       if opequal_p(xsol, '%solve) then error("Unable to solve"),
 582       for xk in xsol do (
 583         xk : real_imagpart_to_conjugate(xk),
 584         xk : merge_solutions(sk,xk),
 585         nsol : map(lambda([s], solution_subst(xk, s, simpfuncs_lambda)), v),
 586         if listp(nsol) and if_p(first(nsol)) then nsol : first(nsol),
 587         if if_p(nsol) and not(listp(second(nsol))) then nsol : %if(first(nsol),
 588           [second(nsol)], third(nsol)),
 589         nsol_vars : set_of_variables(nsol),
 590         %cnds : safer_subst(xk, cnds),
 591         %cnds : expunge_nonrelevant_cnds(%cnds, nsol),
 592         nsol : %if(%cnds, nsol, %union()),
 593         true_sol : adjoin(nsol, true_sol))),
 594     xreduce('%union, true_sol, %union()));
 595
 596 /* Maybe all this should be fixed higher up; for now ...*/
 597
 598 expunge_nonrelevant_cnds (cnds, sol) := block([listconstvars : false, sol_vars : listofvars(sol)],
 599   if cnds = true or cnds = false then cnds
 600   else if disjunction_p(cnds) or conjunction_p(cnds) then (
 601     map(lambda([s], expunge_nonrelevant_cnds(s, sol)), cnds))
 602   else if lfreeof(sol_vars, cnds) and not constantp(cnds) then true
 603   else cnds);
 604