solve とともに（nnsolve その3）
solve は非常に強力な関数で
(%i1) solve(x=y,[x,y]); (%o1) [[x = %r1, y = %r1]]
と解いてくれます（笑）．
これを抑えるには
(%i2) solve(x=y,[x,y]),linsolve_params:false; (%o2) [[x = y]]
とすればよいのですが，次数が上がると，聞く耳を持たないご様子です．
(%i3) solve(x^2=y^2,[x,y]),linsolve_params:false; (%o3) [[x = %r2, y = - %r2], [x = %r3, y = %r3]]
一方
(%i4) nnsolve(x^2=y^2); (%o4) (y - x = 0) %or (y + x = 0)
なのですが，solve による解のパラメトリック表示を，存在量化の形で残すスイッチ nnsqesolve も付けています．これを false にすると，qepmax と nns による最後の整形を行いません．
(%i5) nnsolve(x^2=y^2),nnsqesolve:false;
(%o5) (qe([[E, s1]], (x = - s1) %and (y = s1)))
%or (qe([[E, s1]], (x = s1) %and (y = s1)))
(%i6) %,eval;
(%o6) (y - x = 0) %or (y + x = 0)
(%i7) nnsolve(x^3+y^3+z^3=3*x*y*z);
(%o7) ((y - x = 0) %and (z - x = 0)) %or (z + y + x = 0)
(%i8) nnsolve(x^3+y^3+z^3=3*x*y*z),nnsqesolve:false;
(%o8) (qe([[E, s1]], (x = s1) %and (y = s1) %and (z = s1)))
%or (qe([[E, s1], [E, s2]], (x = (- s2) - s1) %and (y = s2) %and (z = s1)))
(%i9) %,eval;
(%o9) ((y - x = 0) %and (z - x = 0)) %or (z + y + x = 0)
dnf と cnf（nnsolve その2）
選言標準形，連言標準形は代数で言えば，それぞれ展開，因数分解に当たる基本的なものであり，前者は所謂（排他的とは限らない）場合分けなので，方程式系の解の構成には欠かせません．
その実装は簡単で，例えば，A ∧ (B ∨ C) を (A ∧ B) ∨ (A ∧ C) に可能な限り置き換えれば良いだけです．
が，．．．それは連言，選言が 2 項オペレーターの場合の話．%and，%or は nary なので，そのままではパターンにマッチしません．
というわけで，作ったのが以下の dnf と cnf です．
infix(Or,60,40)$
infix(And,60,40)$
/* n項 %and,%or を2項 And,Or の合成に */
orand2orand(f):=block([g],g:f,
if atom(g) then return(g),
if op(g)="%or" then g:xreduce("Or",args(g)),
if op(g)="%and" then g:xreduce("And",args(g)),g)$
/* 分配 */
matchdeclare([aa,bb,cc],true)$
defrule(ru05r, (aa) Or ((bb) And (cc)), ((aa) Or (bb)) And ((aa) Or (cc)))$
defrule(ru08r, ((bb) And (cc)) Or (aa), ((aa) Or (bb)) And ((aa) Or (cc)))$
defrule(ru05a, (aa) And ((bb) Or (cc)), ((aa) And (bb)) Or ((aa) And (cc)))$
defrule(ru08a, ((bb) Or (cc)) And (aa), ((aa) And (bb)) Or ((aa) And (cc)))$
/* subst の並列でも可 */
defrule(ru00r, (aa) Or (bb), (aa) %or (bb))$
defrule(ru00a, (aa) And (bb), (aa) %and (bb))$
cnf(f):=applyb2(applyb2(scanmap(orand2orand,f),ru05r,ru08r),ru00r,ru00a)$
dnf(f):=applyb2(applyb2(scanmap(orand2orand,f),ru05a,ru08a),ru00r,ru00a)$
実行例
Maxima 5.36.1 http://maxima.sourceforge.net
using Lisp CMU Common Lisp 20d (20D Unicode)
Distributed under the GNU Public License. See the file COPYING.
Dedicated to the memory of William Schelter.
The function bug_report() provides bug reporting information.
(%i1) (x=0)%or((x-1=0)%and(y=0)%and(z=0));
(%o1) (x = 0) %or ((x - 1 = 0) %and (y = 0) %and (z = 0))
(%i2) cnf(%);
(%o2) ((x - 1 = 0) %or (x = 0)) %and ((x = 0) %or (y = 0))
%and ((x = 0) %or (z = 0))
(%i3) dnf(%);
(%o3) (x = 0) %or ((x - 1 = 0) %and (x = 0))
%or ((x - 1 = 0) %and (x = 0) %and (y = 0))
%or ((x - 1 = 0) %and (x = 0) %and (z = 0))
%or ((x - 1 = 0) %and (y = 0) %and (z = 0)) %or ((x = 0) %and (y = 0))
%or ((x = 0) %and (y = 0) %and (z = 0)) %or ((x = 0) %and (z = 0))
(%i4) cnf(%);
(%o4) ((x - 1 = 0) %or (x = 0)) %and ((x - 1 = 0) %or (x = 0) %or (y = 0))
%and ((x - 1 = 0) %or (x = 0) %or (y = 0) %or (z = 0))
%and ((x - 1 = 0) %or (x = 0) %or (z = 0)) %and ((x = 0) %or (y = 0))
%and ((x = 0) %or (y = 0) %or (z = 0)) %and ((x = 0) %or (z = 0))
(%i5) dnf(%);
(%o5) (x = 0) %or ((x - 1 = 0) %and (x = 0))
%or ((x - 1 = 0) %and (x = 0) %and (y = 0))
%or ((x - 1 = 0) %and (x = 0) %and (y = 0) %and (z = 0))
%or ((x - 1 = 0) %and (x = 0) %and (z = 0))
%or ((x - 1 = 0) %and (y = 0) %and (z = 0)) %or ((x = 0) %and (y = 0))
%or ((x = 0) %and (y = 0) %and (z = 0)) %or ((x = 0) %and (z = 0))
(%i6) cnf(%);
(%o6) ((x - 1 = 0) %or (x = 0)) %and ((x - 1 = 0) %or (x = 0) %or (y = 0))
%and ((x - 1 = 0) %or (x = 0) %or (y = 0) %or (z = 0))
%and ((x - 1 = 0) %or (x = 0) %or (z = 0)) %and ((x = 0) %or (y = 0))
%and ((x = 0) %or (y = 0) %or (z = 0)) %and ((x = 0) %or (z = 0))
(%i7) nns(qe([],%));
(%o7) (x = 0) %or ((x - 1 = 0) %and (y = 0) %and (z = 0))
nnsolve その１
原理は単純で，入力 f を選言標準形 ∨_k f_k に直し，f_k を方程式系部分 p_k と，不等式・非等式系部分 q_k との連言 p_k ∧ q_k に分解，p_k を maxima の solve で処理した結果 r_k を用いて，∨_k (r_k ∧ q_k) を出力とするものです．
実装において問題となったのは，やはり to_poly_solve の %and，%or です．これらは単なるシンボルではなく，独自の判断で妄動，いや行動するオペレーターなので，．．．
能動的な簡約
QEPCAD B，従って，qepmax では
(%i1) qe([[A,x]],(x=1)%eq(x^2+a*x+b=0));
2
(%o1) (a - 4 b = 0) %and (b + a + 1 = 0)
(%i2) qe([[A,x]],(x^2+x=1)%implies(x^2+a*x=b^2));
4 2 2 2
(%o2) (a - 1 = 0) %and (b + a b - 3 b - a + a + 1 = 0)
のように簡約が充分でない（と人間には見える）出力になる場合があります．
ですが，我々は maxima 上なので，solve を用いれば
(%i3) solve(args(%o1)); (%o3) [[b = 1, a = - 2]] (%i4) solve(args(%o2)); (%o4) [[b = - 1, a = 1], [b = 1, a = 1]]
とできます．
では，不等式との連言をとってみましょう．
(%i5) qe([[A,x]],(x^2+x=1)%implies(x^2+a*x=b^2))%and(b<0);
4 2 2 2
(%o5) (a - 1 = 0) %and (b < 0) %and (b + a b - 3 b - a + a + 1 = 0)
(%i6) solve(args(%));
solve: more equations than unknowns.
Unknowns given :
[b, a]
Equations given:
4 2 2 2
[a - 1 = 0, b < 0, b + a b - 3 b - a + a + 1 = 0]
-- an error. To debug this try: debugmode(true);
もちろんダメです．
ということで，quantifier-free の半代数系に含まれる方程式の次数を下げる（or見易く整理する）関数 nnsolve を作りました．
(%i7) nnsolve(%o1); (%o7) (a + 2 = 0) %and (b - 1 = 0) (%i8) nnsolve(%o2); (%o8) (a - 1 = 0) %and ((b - 1 = 0) %or (b + 1 = 0)) (%i9) nnsolve(%o5); (%o9) (a - 1 = 0) %and (b + 1 = 0) (%i10) qe([],((x^2-1)*(x^2-4)>0)%and(y^3=y)%and(x<y)); (%o10) (x - 1 < 0) %and (x + 1 # 0) %and (x + 2 # 0) %and (((x < 0) %and (x + 1 > 0) %and (y = 0)) %or ((x + 1 > 0) %and (y - 1 = 0)) %or ((x + 2 < 0) %and (y - 1 = 0)) %or ((x + 2 < 0) %and (y = 0)) %or ((x + 2 < 0) %and (y + 1 = 0))) %and (y - 1 <= 0) %and (y + 1 >= 0) (%i11) nnsolvexx(%); (%o11) ((x - 1 < 0) %and (x + 1 > 0) %and (y - 1 = 0)) %or ((x < 0) %and (x + 1 > 0) %and (y = 0)) %or ((x + 2 < 0) %and (y - 1 = 0)) %or ((x + 2 < 0) %and (y = 0)) %or ((x + 2 < 0) %and (y + 1 = 0)) (%i12) qe([[A,x],[A,y]],(%o10)%eq(%o11)); (%o12) true
コードは少し長いので段階的に公開したいと思います．
nnscan
nns は最も外にある %and または %or に作用するのみなので
(%i1) nns(((x>=0)%and(x=0))%or((x>=1)%and(x=1))); (%o1) ((x = 0) %and (x >= 0)) %or ((x = 1) %and (x >= 1)) (%i2) nns((((x>=0)%and(x=0))%or((x>=1)%and(x=1)))%and((y>=0)%and(y=0))); (%o2) (((x = 0) %and (x >= 0)) %or ((x = 1) %and (x >= 1))) %and (y = 0)
となってしまいます．
ということで，scanmap 版を作りました．
nns(f):=
block([cl,L,fk],cl(f,g):=length(f)<=length(g),
if atom(f) then return(f)
elseif op(f)="%and" then
(L:sort(rest(full_listify(powerset(setify(args(f))))),cl),
for k:1 thru length(L) do
(fk:substpart("%and",part(L,k),0),
if qe([],(fk)%implies(f))=true then return(fk)))
elseif op(f)="%or" then
(L:sort(rest(full_listify(powerset(setify(args(f))))),cl),
for k:1 thru length(L) do
(fk:substpart("%or",part(L,k),0),
if qe([],(f)%implies(fk))=true then return(fk)))
else f)$
nnscan(f):=scanmap(nns,f)$
(%i3) nnscan(%o1); (%o3) (x = 0) %or (x = 1) (%i4) nnscan(%o2); (%o4) ((x = 0) %or (x = 1)) %and (y = 0)
scan ではなく，入力をまず連言or選言標準形に書き換えて nns するのが一般的でしょうが，果たして．．．
高速化
アルゴリズムを「簡約」し，コストを半減させました．
nns(f):=
block([cl,L,fk],cl(f,g):=length(f)<=length(g),
L:sort(rest(full_listify(powerset(setify(args(f))))),cl),
if op(f)="%and" then
for k:1 thru length(L) do
(fk:substpart("%and",part(L,k),0),
if qe([],(fk)%implies(f))=true then return(fk))
elseif op(f)="%or" then
for k:1 thru length(L) do
(fk:substpart("%or",part(L,k),0),
if qe([],(f)%implies(fk))=true then return(fk))
else f)$
