#3799 fourier_elim returns non-(in)equations

[Ricardo Buring]

The command fourier_elim should return a list of equations or inequations, or disjunctions of such. The following is fine (solving for both variables):

(%i1) load("fourier_elim");
(%o1) /home/rburing/src/SageMath-9.2/local/share/maxima/5.44.0/share/fourier_e\
lim/fourier_elim.lisp
(%i2) fourier_elim([x0 + 2*x1 + 2 <= 0, -3*x0 + 4*x1 - 4 <= 0], [x0,x1]);
              8         1           4 x1 - 4         1
(%o2) [x0 = - -, x1 = - -] or [x0 = --------, x1 < - -]
              5         5              3             5
                               1      4 x1   4
 or [x0 = (- 2 x1) - 2, x1 < - -] or [---- - - < x0, x0 < (- 2 x1) - 2, 
                               5       3     3
       1
x1 < - -]
       5

Unfortunately the following is not (solving only for the first variable):

(%i3) fourier_elim([x0 + 2*x1 + 2 <= 0, -3*x0 + 4*x1 - 4 <= 0], [x0]);
                                             4 x1 - 4
(%o3) [x0 = (- 2 x1) - 2, 5 x1 + 1] or [x0 = --------, - 6 (5 x1 + 1) > 0]
                                                3
 or [x0 = (- 2 x1) - 2, - 2 (5 x1 + 1) > 0]
     4 x1   4
 or [---- - - < x0, x0 < (- 2 x1) - 2, - 6 (5 x1 + 1) > 0]
      3     3

Here the first list contains an expression without an equality or inequality sign, namely 5 x1 + 1. By substitution into the original inequation, we can see what the inequality should have been:

(%i4) expand(subst(x0=-2*x1 - 2, -3*x0 + 4*x1 - 4 <= 0));
(%o4)                           10 x1 + 2 <= 0

So, it should have been 5 x1 + 1 <= 0. Somehow, the <= 0 gets lost along the way in fourier_elim.

Originally reported in Ask SageMath question #57512.

[Ricardo Buring]

Or rather, it should be = 0 because the < 0 case is listed separately.

[Stavros Macrakis]

The problem appears to be the treatment of variables which are not listed in the variable list. Here is a minimal example:

fourier_elim([a=0],[a]) => [a=0]
fourier_elim([a=0],[]) => [a]

It's not clear to me what fourier_elim is supposed to do with variables that aren't listed in the variable list in the first place. Consider

fourier_elim([a=0,a>0],[a]) => emptyset        <<< correct
fourier_elim([a=0,a>0],[]) => [a > 0, a]       <<< ???

It seems to just return the inputs. Here's an example with multiple variables:

fourier_elim([a=k,a>k],[a]) => emptyset
fourier_elim([a=k,a>k],[k]) => emptyset
fourier_elim([a=k,a>k],[]) => [a - k > 0, k - a]

In this case, some of the returned relations are simplified if the variable is mentioned, but not if it isn't:
~~~
fourier_elim([0<=b+x,b<=x],[b,x]) =>
[b = 0, x = 0] or
[b = x, 0 < x] or
[b = - x, 0 < x] or
[- x < b, b < x, 0 < x]

fourier_elim([0<=b+x,b<=x],[b]) =>
[b = x, x] or
[b = x, 2 x > 0] or
[b = - x, 2 x > 0] or
[- x < b, b < x, 2 x > 0]
~~~

[Barton Willis]

Hi,

I'm testing a fix.

First, for an empty variable list, I think it's OK for fourier_elim to return its first argument (possibly "simplified.")

(%i28)  fourier_elim([a=0],[a]);
(%o28)  [a=0]

(%i29)  fourier_elim([a=0,a>0],[]);
(%o29)  [a>0,a=0]

Second, I'm not sure what fourier_elim is should do with variables that aren't listed in the variable list either. But one thing I'm sure of is that shouldn't return something that is wrong. And it should return something that can be feed back into fourier_elim (after subst for a variable):

(%i30)  fourier_elim([0<=b+x,b<=x],[b]) ;
(%o30)  [b=x,x=0] or [b=x,2*x>0] or [b=-x,2*x>0] or [-x<b,b<x,2*x>0]

If we now decided that x = 107, we can substitute and re-process:

(%i31)  subst(x=107,%);
(%o31)  [b=107,107=0] or [b=107,214>0] or [b=-107,214>0] or [-107<b,b<107,214>0]
(%i32)  fourier_elim(%,[b]);
(%o32)  [-107<b,b<107] or [b=-107] or [b=107]

As for the original bug, with my7 patch, I get

(%i53)  fourier_elim([x0 + 2*x1 + 2 <= 0, -3*x0 + 4*x1 - 4 <= 0], [x0]);
(%o53)  [x0=-2*x1-2,5*x1+1=0] or [x0=(4*x1-4)/3,-6*(5*x1+1)>0] or [x0=-2*x1-2,-2*(5*x1+1)>0] or [(4*x1)/3-4/3<x0,x0<-2*x1-2,-
6*(5*x1+1)>0]

~~~~
Arguably,   ``[x0=-2*x1-2,5*x1+1=0]``  should simplify, but since ``x1``  isn't listed as a variable to solve for, I think it's OK.  At least a user can substitute a value for ``x1`` and process again using fourier_elim:

(%i54) subst(x1=-3,%);
(%o54) [x0=4,-14=0] or [x0=-16/3,84>0] or [x0=4,28>0] or [-16/3<x0,x0<4,84>0]

(%i55) fourier_elim(%,[x0]);
(%o55) [-16/3<x0,x0<4] or [x0=4] or [x0=-16/3]
~~~

Thank you for your interest in fourier_elim and for the bug report. Please let me know what you all think of my putative fix.

[Barton Willis]

My putative fix:

(defun linear-elimination (l v)
  (let (($linsolve_params nil) ($backsubst t) ($programmode t) 
          ($linsolvewarn nil) ($globalsolve nil) (subs) (unsolved) (vars))

    (setq l ($elim l v))
    (cond (($member 1 ($first l)) '$emptyset)
      (t
       (setq subs ($linsolve ($second l) v))
       (setq vars (mapcar '$lhs (margs subs)))
       (setq vars (push '(mlist) vars))
     (setq unsolved ($first l))
     (setq unsolved (cons '(mlist) (mapcar #'(lambda (q) (take '(mequal) q 0)) (cdr ($first l)))))
     (simplifya (list '(mlist) subs unsolved vars) t)))))

[Barton Willis]

Fixed by commit a0e80d . Regression tests added; closing bug