# CAS Best Practices

**URL:** <https://discourse.julialang.org/t/cas-best-practices/58092>\
**Category:** Machine Learning\
**Tags:** symbolic\
**Created:** [March 27, 2021, 8:51pm UTC](https://discourse.julialang.org/t/cas-best-practices/58092 "2021-03-27T20:51:11Z")\
**Posts on this page:** 10\
**Page:** 5

<div class="post-metadata">

**Author:** ![rfateman](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rfateman/32/22507_2.png) [@rfateman](https://discourse.julialang.org/u/rfateman)\
**Post date:** [April 6, 2021, 6:22pm UTC](https://discourse.julialang.org/t/cas-best-practices/58092/89 "2021-04-06T18:22:40Z")

</div>

One has to be rather careful, and I don’t know what the Julia philosophy dictates,.

As has been pointed out, the notion of principal square root of **a non-negative real number** is just fine. A very careful CAS might allow you to create and manipulate a “rootof” expression, e.g. RootOf(y^2=x,y) . This represents one of the square roots of x. You might even distinguish the two roots by an index, RootOf(y^2=x,1) and RootOf(y^2=x,2). Though which one is the “positive” or “principal” root is not identifiable. You do know that their sum is 0 and their product is x.

If sqrt() is a name for a program that operates only on non-neg reals, then identifying it as positive real principle square root is certainly ok. It has a certain relationship with the square-root concept in algebra, or in a CAS, but it seems to me that they are not the same.

So some of the comments are, in my view, right on track. Here’s some food for thought.

If z= x^2+y^2-2_x_y, what do you want for square\_root(z)?  
It could be x-y, y-x, abs(x-y), RootOf(w=z,w), or a set {x-y,y-x}. There are various arguments for each result. Maxima gives you something of a choice. Leaving it as RootOf, “unevaluated” might seem to be the most correct, but it makes further computation difficult and potentially ambiguous.

Oh, it is also possible sometimes to distinguish RootOfs by intervals. e.g. RootOf(y^2=9,y, inside [-inf,0]) which would presumably be unambiguously -3.

Have fun,.  
RJF

---

<div class="post-metadata">

**Author:** ![goretkin](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/goretkin/32/167_2.png) [@goretkin](https://discourse.julialang.org/u/goretkin)\
**Post date:** [April 7, 2021, 6:40pm UTC](https://discourse.julialang.org/t/cas-best-practices/58092/90 "2021-04-07T18:40:19Z")

</div>

> [@ChrisRackauckas](#):
>
> It’s worth a discussion on the SymbolicUtils.jl repo. Open an issue.

> <https://github.com/JuliaSymbolics/SymbolicUtils.jl/issues/269>
>
> (paraphrased and excerpted from \[a discussion on discourse\](https://discourse.ju…lialang.org/t/cas-best-practices/58092/78))
> 
> Over the domain \`x::AbstractFloat\`, the rule that rewrites \`x - x\` to \`zero(x)\` is wrong, since the domain includes \`NaN\`s and \`Inf\`s. However:
> 
> \`\`\`julia
> julia\> y = SymbolicUtils.Sym{AbstractFloat}(:y)
> y
> 
> julia\> y - y
> 0
> \`\`\`
> 
> Over the "real numbers" (or any \[field\](https://en.wikipedia.org/wiki/Field\_(mathematics))), the rule is correct, and that is where that rule is coming from.
> 
> It strikes me as un-Julian design to have that rule defined on \`Sym{Real}\` (or \`Sym{Number}\`), and then have a specialization that turns it off for \`AbstractFloat\`. Specializing behavior is essential in the Julia ecosystem for generic programming, but specializations are not, for lack of a better word, done "willy nilly". They are supposed to be consistent with a generic definition, but perhaps more efficient or producing a more compact (symbolic!) representation. For example:
> 
> \`\`\`julia
> julia\> typeof(1:4)
> UnitRange{Int64}
> 
> julia\> typeof(1:4) \<: AbstractArray
> true
> 
> julia\> (1:4) .+ 5
> 6:9
> 
> julia\> collect(1:4) .+ 5
> 4-element Vector{Int64}:
> 6
> 7
> 8
> 9
> 
> julia\> collect((1:4) .+ 5) == collect(1:4) .+ 5 # consistency
> true
> \`\`\`
> 
> If broadcasted addition on \`UnitRange{Int64}\` were defined in such a way that the “# consistency” property did not hold, that would be bad for generic programming!
> 
> Does that badness apply in this CAS setting? What consistency should exist across rules that are defined on types that have a subtype relationship among them? I am not sure. That’s what the question is.
> 
> My instinct is that this CAS should introduce a new type, e.g. \`Symbolics.Real\`, not \`Base.Real\` (though perhaps Symbolics.Real \<: \`Base.Real\`). Defining that simplification rule on \`Base.Real\` seems wrong since I can give you an \`x::Base.Real\` for which the rule is wrong. But I should not be able to give you an \`y::Symbolics.Real\` for which the rule is wrong.
> 
> It's already not the case that \`Sym{Real}\<:Real\`. What is the benefit of re-using this abstract type?

---

<div class="post-metadata">

**Author:** ![oscarbenjamin](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscarbenjamin/32/22868_2.png) [@oscarbenjamin](https://discourse.julialang.org/u/oscarbenjamin)\
**Post date:** [April 13, 2021, 10:48pm UTC](https://discourse.julialang.org/t/cas-best-practices/58092/91 "2021-04-13T22:48:36Z")

</div>

> [@Fredrik\_Johansson](#):
>
> The more general issue here is that computer algebra systems historically are very confused about what it means to solve something. Does solve() return some particular solution, the set of solutions, a set containing some of the solutions, a superset of the solutions (possibly containing spurious solutions), real solutions, complex solutions, or something random depending on context or the phase of the moon? I’ve seen all combinations.

Another fundamental question in the interpretation of a “solution” is how a solver handles equations parametrised by symbols e.g.:

```julia
julia> @variables a, b, x
(a, b, x)

julia> Symbolics.solve_for([a*x ~ b], [x])
1-element Array{Num,1}:
 b*(a^-1)

```

Most CAS will return something similar but what if `a` is zero?

More generally given a situation where solutions exist conditionally depending on the unknowable values of symbolic parameters how do you define formally what the interpretation should be on the returned result? Does it apply for “most” values of `a` or would the solver only return something if it was known to be correct for all values of `a`? What if an expression for the solution is correct for `a < 0` but no solution exists for `a > 0`?

Different users want different things in different situations so there isn’t a one size fits all approach. Many users would be extremely disappointed if `solve` was unable to get `b/a` from such a simple equation but in some contexts (especially for internal use) something more robust is needed. The next level up is to have a way to indicate any assumptions that were made in deriving the returned result so what you return is

```julia
>>> solve([a*x ~ b], [x])
b/a provided a != 0

```

At minimum the primitive solvers on which others are based need to have this capability so that robust higher-level solvers can be built on top. A friendly user API like `solve` might not do this but for power-users and internal use this is needed.

What most people as “users” want is to get the solution case that corresponds to the generic case for the parameters _except_ when it comes to determining that no solution exists e.g.:

```julia
In [5]: linsolve_cond([Eq(x + y, 1), Eq(x + y, a), Eq(x, 1)], [x, y])
Out[5]: 
⎧{(1, 0)} for a - 1 = 0
⎨                       
⎩ ∅ otherwise

```

Here the solution is generically (for almost all values of `a`) the empty set but many people would consider it a bug to return that: why can’t the solver figure out that `a=1`?

Another approach is to try to consider all possible cases but it doesn’t scale well. I made a sketch of this for SymPy here:

> <https://github.com/sympy/sympy/issues/16861>
>
> There are a number of open issues about proper solutions of under/over-determine…d systems of equations. I think that the fundamental problem is that we don't have a way to return conditional solution sets (from solve or solveset). I've written a linear equation solver that (ab)uses Piecewise to do this and I'll show the code below but examples first.
> 
> Consider the equation \`a\*x=b\`. We can solve this easily:
> \`\`\`julia
> In \[2\]: a, b, c, d, e = symbols('a, b, c, d, e')                                                                                  
> 
> In \[3\]: x, y, z = symbols('x, y, z')                                                                                              
> 
> In \[4\]: solve(a\*x-b, x)                                                                                                           
> Out\[4\]: 
> ⎡b⎤
> ⎢─⎥
> ⎣a⎦
> 
> In \[5\]: solveset(a\*x-b, x)                                                                                                        
> Out\[5\]: 
> ⎧b⎫
> ⎨─⎬
> ⎩a⎭
> \`\`\`
> The problem with this is that these solutions are not valid for \`a=0\`. The solver I have written gives
> \`\`\`julia
> In \[7\]: s = linsolve\_cond(\[a\*x-b\], \[x\])                                                                                           
> 
> In \[8\]: s                                                                                                                         
> Out\[8\]: 
> ⎧ ℂ for a = 0 ∧ b = 0
> ⎪                         
> ⎪ ∅ for a = 0 ∧ b ≠ 0
> ⎨                         
> ⎪⎧⎛b ⎞⎫                   
> ⎪⎨⎜─,⎟⎬ for a ≠ 0    
> ⎩⎩⎝a ⎠⎭         
> \`\`\`
> This correctly handles the three possible cases for the solution set. We can use subs with the result:
> \`\`\`julia
> In \[9\]: s.subs(a, 1)                                                                                                              
> Out\[9\]: {(b,)}
> 
> In \[10\]: s.subs(a, 0)                                                                                                             
> Out\[10\]: 
> ⎧ℂ for b = 0
> ⎨            
> ⎩∅ otherwise
> 
> In \[11\]: s.subs(a, 0).subs(b, 1)                                                                                                  
> Out\[11\]: ∅
> \`\`\`
> That latter example fails with the output of solve/solveset:
> \`\`\`julia
> In \[19\]: solveset(a\*x-b, x)                                                                                                       
> Out\[19\]: 
> ⎧b⎫
> ⎨─⎬
> ⎩a⎭
> 
> In \[20\]: solveset(a\*x-b, x).subs(a, 0)                                                                                            
> Out\[20\]: {zoo⋅b}
> 
> In \[21\]: solveset(a\*x-b, x).subs(a, 0).subs(b, 1)                                                                                 
> Out\[21\]: {zoo}
> 
> In \[22\]: solve(a\*x-b, x)\[0\].subs(a, 0).subs(b, 1)                                                                                 
> Out\[22\]: zoo
> \`\`\`
> 
> The solver I have written uses split-recursion. It performs row-reduction and splits recursively each time it encounters a pivot that isn't known to be zero. This means it \*should\* give a correct but possibly complicated answer for any system of linear equations although I haven't tested it extensively. Here's some examples:
> \`\`\`julia
> In \[23\]: linsolve\_cond(\[a\*x+y, x+y\], \[x, y\])                                                                                      
> Out\[23\]: 
> ⎧ {(0, 0)} for a ≠ 1
> ⎪                                
> ⎨⎧⎛-τ₀ ⎞ ⎫           
> ⎪⎨⎜────, τ₀⎟ | τ₀ ∊ ℂ⎬ otherwise
> ⎩⎩⎝ a ⎠ ⎭           
> 
> In \[24\]: linsolve\_cond(\[a\*x+y, x+y+a-1\], \[x, y\])                                                                                  
> Out\[24\]: 
> ⎧ {(1, 0)} for a = 0             
> ⎪                                                          
> ⎪⎧⎛-τ₀ ⎞ ⎫                                     
> ⎪⎨⎜────, τ₀⎟ | τ₀ ∊ ℂ⎬ for a = 1             
> ⎨⎩⎝ a ⎠ ⎭                                     
> ⎪                                                          
> ⎪⎧⎛-(1 - a) 1 - a⎞⎫                                      
> ⎪⎨⎜─────────, ─────⎟⎬ for (a \> 0 ∧ a \< 1) ∨ a \> 1 ∨ a \< 0
> ⎩⎩⎝a⋅(a - 1) a - 1⎠⎭                                      
> 
> In \[25\]: linsolve\_cond(\[x+y, x+y+a\], \[x, y\])                                                                                      
> Out\[25\]: 
> ⎧{(-τ₀, τ₀) | τ₀ ∊ ℂ} for a = 0
> ⎨                               
> ⎩ ∅ otherwise
> 
> In \[26\]: linsolve\_cond(\[x+y, x+y+a, x-y\], \[x, y\])                                                                                 
> Out\[26\]: 
> ⎧{(0, 0)} for a = 0
> ⎨                   
> ⎩ ∅ otherwise
> 
> In \[27\]: linsolve\_cond(\[a\*x+b\*y, c\*x+d\*y\], \[x, y\])                                                                                
> Out\[27\]: 
> ⎧ 2                                                                                                                   
> ⎪ ℂ for a = 0 ∧ b = 0 ∧ c = 0 ∧ d = 0                            
> ⎪                                                                                                                               
> ⎪ {(τ₀, 0) | τ₀ ∊ ℂ} for a = 0 ∧ c = 0 ∧ (b ≠ 0 ∨ d ≠ 0)                           
> ⎪                                                                                                                               
> ⎪⎧⎛-d⋅τ₀ ⎞ ⎫                                                                                                        
> ⎪⎨⎜──────, τ₀⎟ | τ₀ ∊ ℂ⎬ for a = 0 ∧ b = 0 ∧ c ≠ 0                                
> ⎨⎩⎝ c ⎠ ⎭                                                                                                        
> ⎪                                                                                                                               
> ⎪ {(0, 0)} for (a = 0 ∨ a⋅d - b⋅c ≠ 0) ∧ (a ≠ 0 ∨ b ≠ 0) ∧ (a ≠ 0 ∨ c ≠ 0) ∧ (b ≠ 0 ∨ a⋅d - b⋅c ≠ 0) ∧ (c ≠ 0 ∨ a⋅d - b⋅c ≠ 0)
> ⎪                                                                                                                               
> ⎪⎧⎛-b⋅τ₀ ⎞ ⎫                                                                                                        
> ⎪⎨⎜──────, τ₀⎟ | τ₀ ∊ ℂ⎬ for a⋅d - b⋅c = 0 ∧ a ≠ 0                                
> ⎪⎩⎝ a ⎠ ⎭                                                                                                        
> ⎩                                                                                                                               
> \`\`\` 
> 
> I think that SymPy really needs something like this although I'm not sure exactly how it should be integrated:
> 1. Should this be what linsolve (or solveset) does?
> 2. Should solveset generally return solutions of this form?
> 3. Piecewise is a subset of Expr whereas Set is a subset of Basic so I don't think Piecewise should be used here (but a similar class could be written)
> 4. The conditions in the output could be simpler in some cases. I'm not sure exactly how to simplify them though.
> 
> The code for the solver is here:
> \`\`\`python
> from sympy import \*
> 
> 
> def linsolve\_cond(eqs, unknowns, unique=False):
> if not eqs:
> return S.Complexes\*\*len(unknowns)
> 
> # Preprocessing
> A, b = linear\_eq\_to\_matrix(eqs, unknowns)
> Aaug = Matrix.hstack(A, b).tolist()
> 
> # Main workhorse:
> sols\_conds = \_linsolve\_cond(Aaug)
> 
> # sols\_conds is a list of 3-tuples:
> # \[(solset, pivot\_conds, consistency\_conds),...\]
> #
> # solset: solution set as a FiniteSet or ImageSet
> # pivot\_conds: list of conditions (e.g. a!=0) assumed in pivoting
> # consistency\_conds: list of conditions needed for existence of solutions
> 
> # Build all the separate cases into a Piecewise:
> sets\_conds = \[\]
> for solset, pivot\_conds, consistency\_conds in sols\_conds:
> pivot\_cond = And(\*pivot\_conds)
> consistency\_cond = And(\*consistency\_conds)
> if consistency\_cond is not S.false:
> sets\_conds.append((solset, pivot\_cond & consistency\_cond))
> if consistency\_cond is not S.true:
> sets\_conds.append((S.EmptySet, pivot\_cond & Not(consistency\_cond)))
> 
> sets\_conds\_d = {}
> for ss, conds in sets\_conds:
> if ss not in sets\_conds\_d:
> sets\_conds\_d\[ss\] = conds
> else:
> sets\_conds\_d\[ss\] = Or(sets\_conds\_d\[ss\], conds)
> 
> if unique:
> sets\_conds\_d = {s: c for s, c in sets\_conds\_d.items() if isinstance(s, FiniteSet)}
> 
> return Piecewise(\*sets\_conds\_d.items())
> 
> 
> def \_linsolve\_cond(Aaug, \_recurse=None):
> Nr, Nc = len(Aaug), len(Aaug\[0\])
> 
> Aorig = Matrix(Aaug)
> 
> if \_recurse is None:
> row, col, pivots, pivot\_conds = 0, 0, \[\], \[\]
> else:
> row, col, pivots, pivot\_conds = \_recurse
> 
> if pivots:
> row, col = pivots\[-1\]
> row += 1
> col += 1
> else:
> row, col = 0, 0
> 
> sols\_conds = \[\]
> 
> # Call self recursively for alternate pivots
> def recurse\_zero\_pivot(r, c):
> pivot = Aaug\[r\]\[c\]
> Aaugr = \[\[Arc.subs(pivot, 0) for Arc in Arow\] for Arow in Aaug\]
> pivot\_condsr = pivot\_conds\[:\] + \[Eq(pivot, 0)\]
> \_recurse = (r, c, pivots\[:\], pivot\_condsr)
> sols\_conds.extend(\_linsolve\_cond(Aaugr, \_recurse=\_recurse))
> 
> while row \< Nr and col \< Nc-1:
> # Find pivot row and swap into position
> for r in range(row, Nr):
> is\_zero = Aaug\[r\]\[col\].is\_zero
> if not is\_zero:
> if is\_zero is None:
> # Recurse for the case that the pivot is zero
> recurse\_zero\_pivot(r, col)
> pivot\_conds.append(Ne(Aaug\[r\]\[col\], 0))
> if r != row:
> Aaug\[r\], Aaug\[row\] = Aaug\[row\], Aaug\[r\]
> break
> else:
> # All zeros, next column
> col += 1
> continue
> 
> if pivots:
> assert pivots\[-1\]\[0\] != row
> 
> pivots.append((row, col))
> pivot\_row = Aaug\[row\]
> pivot\_div = Aaug\[row\]\[col\]
> for r in range(row+1, Nr):
> pivot\_mul = Aaug\[r\]\[col\]
> if pivot\_mul.is\_zero:
> continue
> Aaug\[r\]\[col\] = S.Zero
> for c in range(col+1, Nc):
> Aaug\[r\]\[c\] = Aaug\[r\]\[c\]\*pivot\_div - pivot\_row\[c\]\*pivot\_mul
> 
> # Next row/column...
> row += 1
> col += 1
> 
> # Back substitute and list of possibilities
> sol\_set, consistency\_conds = \_back\_substitute(Aaug, pivots)
> 
> sols\_conds.append((sol\_set, pivot\_conds, consistency\_conds))
> 
> return sols\_conds
> 
> 
> def \_back\_substitute(Aaug, pivots):
> Nc = len(Aaug\[0\])
> 
> # Check conditions for existence of solutions then find solutions by
> # back-substitution below
> consistency\_conds = \[\]
> for row in reversed(range(len(Aaug))):
> is\_zero = \[e.is\_zero for e in Aaug\[row\]\]
> if not all(x is True for x in is\_zero\[:-1\]):
> break
> elif is\_zero\[-1\] is False:
> consistency\_conds.append(S.false)
> elif is\_zero\[-1\] is None:
> consistency\_conds.append(Eq(Aaug\[row\]\[-1\], 0))
> 
> assert (row == 0 and not pivots) or row == pivots\[-1\]\[0\]
> 
> # Matrix of all zeros?
> if not pivots:
> solset = S.Complexes\*\*(Nc-1)
> return solset, consistency\_conds
> 
> gen = numbered\_symbols('tau')
> params = \[\]
> sol = \[None\] \* (Nc-1)
> 
> pivots\_cols = {c:r for r, c in pivots}
> for col in reversed(range(Nc-1)):
> if col in pivots\_cols:
> r = pivots\_cols\[col\]
> lhsterms = (Aaug\[r\]\[c\]\*sol\[c\] for c in range(col+1, Nc-1))
> sol\[col\] = (Aaug\[r\]\[-1\] - Add(\*lhsterms)) / Aaug\[r\]\[col\]
> else:
> # Non-pivot gets a free symbol
> sym = next(gen)
> params.append(sym)
> sol\[col\] = sym
> 
> if params:
> solset = ImageSet(Lambda(tuple(params), tuple(sol)), \*\[S.Complexes\]\*len(params))
> else:
> solset = FiniteSet(tuple(sol))
> 
> return solset, consistency\_conds
> 
> 
> x, y, z = symbols('x, y, z')
> a, b, c, d, e = symbols('a, b, c, d, e', finite=True)
> 
> unknowns = (x, y, z)
> eqs = \[sqrt(3)\*x+y, sqrt(2)\*z\]
> 
> sol = linsolve\_cond(eqs, unknowns)
> pprint(sol)
> 
> 
> M = Matrix(symbols('M:9')).reshape(3, 3)
> xs = Matrix(symbols('x:3'))
> b = Matrix(symbols('b:3'))
> sol3 = linsolve\_cond(list(M\*xs - b), list(xs), unique=True)
> \`\`\`

With that you get:

```julia
In [3]: linsolve_cond([Eq(a*x, b)], [x])
Out[3]: 
⎧⎧⎛b ⎞⎫           
⎪⎨⎜─,⎟⎬ for a ≠ 0
⎪⎩⎝a ⎠⎭           
⎪                 
⎨ ∅ for b ≠ 0
⎪                 
⎪ 1             
⎪ ℂ otherwise
⎩  

```

The problem is that the number of cases grows exponentially in the number of symbols:

```julia
In [2]: a, b, c, d, e, f, x, y = symbols('a:f, x:y')

In [3]: linsolve_cond([Eq(a*x + b*y, e), Eq(c*x + d*y, f)], [x, y])
Out[3]: 
⎧⎧⎛ b⋅(a⋅f - c⋅e) ⎞⎫                                                                                                                                                                                                                        
⎪⎪⎜- ───────────── + e ⎟⎪                                                                                                                                                                                                                        
⎪⎨⎜ a⋅d - b⋅c a⋅f - c⋅e⎟⎬                                                                                                                                                                                                                        
⎪⎪⎜───────────────────, ─────────⎟⎪ for a ≠ 0 ∧ a⋅d - b⋅c ≠ 0                                                                                               
⎪⎩⎝ a a⋅d - b⋅c⎠⎭                                                                                                                                                                                                                        
⎪                                                                                                                                                                                                                                                          
⎪ ∅ for (a⋅d - b⋅c = 0 ∧ a ≠ 0 ∧ a⋅f - c⋅e ≠ 0) ∨ (a = 0 ∧ b = 0 ∧ c ≠ 0 ∧ e ≠ 0) ∨ (a = 0 ∧ c = 0 ∧ b ≠ 0 ∧ b⋅f - d⋅e ≠ 0) ∨ (a = 0 ∧ b = 0 ∧ c = 0 ∧ d = 0 ∧ ¬(e = 0 ∧ f = 0)) ∨ (a = 0 ∧ b = 0 ∧ c = 0 ∧ d ≠ 0 ∧ e ≠ 0)
⎪                                                                                                                                                                                                                                                          
⎪ ⎧⎛-b⋅τ₀ + e ⎞ │ ⎫                                                                                                                                                                                                                            
⎪ ⎨⎜─────────, τ₀⎟ │ τ₀ ∊ ℂ⎬ for a⋅d - b⋅c = 0 ∧ a⋅f - c⋅e = 0 ∧ a ≠ 0                                                                                       
⎪ ⎩⎝ a ⎠ │ ⎭                                                                                                                                                                                                                            
⎪                                                                                                                                                                                                                                                          
⎪ ⎧⎛ e⎞ │ ⎫                                                                                                                                                                                                                                
⎪ ⎨⎜τ₀, ─⎟ │ τ₀ ∊ ℂ⎬ for a = 0 ∧ c = 0 ∧ b⋅f - d⋅e = 0 ∧ b ≠ 0                                                                                       
⎪ ⎩⎝ b⎠ │ ⎭                                                                                                                                                                                                                                
⎪                                                                                                                                                                                                                                                          
⎨ ⎧⎛ f⎞ │ ⎫                                                                                                                                                                                                                                
⎪ ⎨⎜τ₀, ─⎟ │ τ₀ ∊ ℂ⎬ for a = 0 ∧ b = 0 ∧ c = 0 ∧ e = 0 ∧ d ≠ 0                                                                                       
⎪ ⎩⎝ d⎠ │ ⎭                                                                                                                                                                                                                                
⎪                                                                                                                                                                                                                                                          
⎪ ⎧⎛-d⋅τ₀ + f ⎞ │ ⎫                                                                                                                                                                                                                            
⎪ ⎨⎜─────────, τ₀⎟ │ τ₀ ∊ ℂ⎬ for a = 0 ∧ b = 0 ∧ e = 0 ∧ c ≠ 0                                                                                           
⎪ ⎩⎝ c ⎠ │ ⎭                                                                                                                                                                                                                            
⎪                                                                                                                                                                                                                                                          
⎪ ⎧⎛ d⋅e ⎞⎫                                                                                                                                                                                                                                  
⎪ ⎪⎜f - ─── ⎟⎪                                                                                                                                                                                                                                  
⎪ ⎨⎜ b e⎟⎬                                                                                                                                                                                                                                  
⎪ ⎪⎜───────, ─⎟⎪ for a = 0 ∧ b ≠ 0 ∧ c ≠ 0                                                                                               
⎪ ⎩⎝ c b⎠⎭                                                                                                                                                                                                                                  
⎪                                                                                                                                                                                                                                                          
⎪ 2                                                                                                                                                                                                                                        
⎪ ℂ for a = 0 ∧ b = 0 ∧ c = 0 ∧ d = 0 ∧ e = 0 ∧ f = 0                                                                                   
⎩  

```

Although the number of cases explodes it is possible to implement a solver this way lazily so that if you know what kind of cases you want then they can be extracted efficiently without computing all of the others. Enumerating all cases is impractical but there needs to be some way to know what assumptions were made in deriving the solution.

Attaching validity conditions to what is returned by solvers needs to be an API and implementation consideration from the outset and it needs to start with the lowest level primitive routines. This is not really a feature that can be added to an implementation retrospectively.

Either way the most important thing is to define (and document!) in formal terms what the solver does and how the returned solutions should be interpreted. That needs to be clear from the outset or the implementation will end up choosing based on the phases of the moon as Fredrik says.

---

<div class="post-metadata">

**Author:** ![rfateman](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rfateman/32/22507_2.png) [@rfateman](https://discourse.julialang.org/u/rfateman)\
**Post date:** [April 15, 2021, 8:17pm UTC](https://discourse.julialang.org/t/cas-best-practices/58092/92 "2021-04-15T20:17:38Z")

</div>

The problem cited above about the proliferation of cases is, I think, only one aspect of the issue. True, you get more cases (we called them provisos… XXX provided YYY), and their numbers can grow exponentially as solutions bifurcate previous provisos. We used this to encode material from tables of integrals in a program Tilu.

The adjacent issue I see is to simplify and combine the provisos to a minimal set, and should you be so lucky as to conclude something like

y=1 provided a \>0  
and  
y=1 provided a\<=0

then computationally maybe conclude that y=1 [do we have to specify … provided a is real??]

I think it is useful to talk about formal computations in a domain which allow (say) cancelling of polynomials gcd in A/B; on the other hand, certain manipulations that are less well structured and are true “generically” but false “generally” are troublesome. Simple examples of these are sometimes posed as puzzles: Find the flaw in the proof … (where the proof concludes that 1=2, say)  
Protesting that a solution has (perhaps) an infinite set of places where it is false contrasts with the certainty that “but the program came up with this solution”. Mathematica does a fair amount of this; I don’t know about SymPy.

There are some tools in Mathematica that are more careful that Solve. I think Reduce and Eliminate are available. I expect that Mathematica’s further manipulation of the results of these programs falls short of the ideal, but I haven’t tried mucking about with their results.

The idea of leaving things “unevaluated” is pretty much inherent in Macsyma/Maxima as a kind of default -when you can’t do an integration, leave the result as “integrate(…)”.  
Leaving a result as a program, unevaluated, by default, seems a bit cruel to the user. e.g.

a+b becomes … if a is a number and b is a number then if a+b is a representable number, then add(a,b).  
a/b becomes … if b is not 0 or NaN …

and sqrt(x) checks that x is real, non-negative …

It works better if your business is not “computer algebra” but “generating code snippets”.

sorry if I’ve misrepresented what I’ve read but not studied in detail here.  
Have fun.

---

<div class="post-metadata">

**Author:** ![goretkin](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/goretkin/32/167_2.png) [@goretkin](https://discourse.julialang.org/u/goretkin)\
**Post date:** [April 15, 2021, 11:09pm UTC](https://discourse.julialang.org/t/cas-best-practices/58092/93 "2021-04-15T23:09:23Z")

</div>

> [@rfateman](#):
>
> I think it is useful to talk about formal computations in a domain which allow (say) cancelling of polynomials gcd in A/B; on the other hand, certain manipulations that are less well structured and are true “generically” but false “generally” are troublesome. Simple examples of these are sometimes posed as puzzles: Find the flaw in the proof … (where the proof concludes that 1=2, say)

Is there a way to get a handle on this distinction between “generally” and “generically”? I think what you’re referring to is e.g. the rule that  
`b * a / b == a`

with or without (respectively) checking that  
`b ≠ 0`.

---

<div class="post-metadata">

**Author:** ![rfateman](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rfateman/32/22507_2.png) [@rfateman](https://discourse.julialang.org/u/rfateman)\
**Post date:** [April 16, 2021, 12:01am UTC](https://discourse.julialang.org/t/cas-best-practices/58092/94 "2021-04-16T00:01:54Z")

</div>

See, for example, the wikipedia article about the field of rational functions

> **[Rational function](https://en.wikipedia.org/wiki/Rational_function)**
>
> In mathematics, a rational function is any function that can be defined by a rational fraction, which is an algebraic fraction such that both the numerator and the denominator are polynomials. The coefficients of the polynomials need not be rational numbers; they may be taken in any field K. In this case, one speaks of a rational function and a rational fraction over K. The values of the variables may be taken in any field L containing K. Then the domain of the function is the set of the values...

for a discussion which allows that one can consider the equivalence class of functions including  
R=(f(x)\*g(x)) / g(x) as well as f(x), by “cancelling” g(x). Even if g(x) might, for some values of x, evaluate to zero, R is computationally equivalent to f(x).

There’s a much broader scope of expressions that might come out of “solve” but are difficult to justify except if you are willing to avert your gaze and say …  
( [https://www.youtube.com/watch?v=vJXU7EVXs2A](https://www.youtube.com/watch?v=vJXU7EVXs2A) )

---

<div class="post-metadata">

**Author:** ![goretkin](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/goretkin/32/167_2.png) [@goretkin](https://discourse.julialang.org/u/goretkin)\
**Post date:** [April 16, 2021, 1:25am UTC](https://discourse.julialang.org/t/cas-best-practices/58092/95 "2021-04-16T01:25:30Z")

</div>

> [@rfateman](#):
>
> Even if g(x) might, for some values of x, evaluate to zero, R is computationally equivalent to f(x).

It certainly doesn’t feel right to say they are “computationally” equivalent. They indeed express different computations and are not [extensionally nor intensionally](https://en.wikipedia.org/wiki/Deductive_lambda_calculus#Intensional_versus_extensional_equality) equal. But I get they are equivalent possibly up to removable singularities.

The article warns about inadvertently losing removable singularities.

> The rational function f(x)={\tfrac{x}{x}} is equal to 1 for all _x_ except 0, where there is a [removable singularity](https://en.wikipedia.org/wiki/Removable_singularity). The sum, product, or quotient (excepting division by the zero polynomial) of two rational functions is itself a rational function. However, the process of reduction to standard form may inadvertently result in the removal of such singularities unless care is taken. Using the definition of rational functions as equivalence classes gets around this, since _x_ / _x_ is equivalent to 1/1.

My original question was about “general” vs “generic”, but I didn’t understand the dichotomy here. Perhaps it’s that “equivalence up to removable singularity” is more “general” and less “generic” than “equivalence, and strictly so”.

---

<div class="post-metadata">

**Author:** ![rfateman](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rfateman/32/22507_2.png) [@rfateman](https://discourse.julialang.org/u/rfateman)\
**Post date:** [April 16, 2021, 5:31pm UTC](https://discourse.julialang.org/t/cas-best-practices/58092/96 "2021-04-16T17:31:36Z")

</div>

My view is that this rational field division by maybe zero is “small potatoes”. There is probably a mathematical excuse buried in the literature for excusing the cancellation of a potentially-zero denominator. Maybe in “valuation” theory. Not something I am willing to spend time studying.

What is potentially far more annoying are situations like this. We can agree that limit(1/x, x-\>0) is something. Call it infinity. Maybe there’s +,-, complex infinity. But it is something. Call it q. Now for an ordinary symbol x that might possibly have a value q, we would like x-x to be 0. 2\*x to be different from x unless x==0, etc etc. This is not true for inf. So we have a problem. Maybe don’t allow limit? Maybe don’t allow operations on something that might later be assigned the value inf? Come up with an arithmetic that allows infinities and is consistent with the rest of the system? (There are arguments that interval arithmetic can do some of this.)  
Anyway, there are a number of different resolutions to this, and they aren’t obvious. E.g. you might say, of course we have to allow limits. Uh, no, that’s one resolution. Maybe not a favorite.  
Now for systems that are heavily oriented toward arithmetic, note that there are IEEE float representations for ±inf, NaN (many of them), signed zero, as well as traps, rounding modes, …  
and if you want to model them in a computer algebra “system”, it requires thought.  
You are aware, but don’t obsess over it, but you know you need to think about arithmetic even if it doesn’t overflow, underflow, etc… You can still make grevious mistakes by subtraction of two nearly-equal quantities, which elevates round-off and truncation noise to prominence.

I would not rely on the usual “mathematics” sources to come up with a unified theory of all of these, including objects infinities, undefineds, indeterminate, NaNs. Nor would I expect this to come out as a free consequence of computer language types, compiling, multiple-dispatch-object-oriented-functional-parallel-cloud-AI-whatever. Ultimately, I do expect that a computer system could mimic or improve on what human mathematicians try to do, just it is not a “free” consequence of something not specifically addressing these issues. For people who are concerned with “reliable computation” that may sort of be adjacent to the symbolic computation in pursuit of better numerics, there is the activity in interval arithmetic which may be helpful. A place to start may be this …  
[https://standards.ieee.org/standard/1788\_1-2017.html](https://standards.ieee.org/standard/1788_1-2017.html)

RJF

---

<div class="post-metadata">

**Author:** ![goretkin](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/goretkin/32/167_2.png) [@goretkin](https://discourse.julialang.org/u/goretkin)\
**Post date:** [April 16, 2021, 5:56pm UTC](https://discourse.julialang.org/t/cas-best-practices/58092/97 "2021-04-16T17:56:03Z")

</div>

> [@rfateman](#):
>
> Call it q. Now for an ordinary symbol x that might possibly have a value q, we would like x-x to be 0. 2\*x to be different from x unless x==0, etc etc. This is not true for inf. So we have a problem.

I am certainly underthinking this. But the problem seems to be an absence of types. Or perhaps “domains”. `x - x == 0` is a perfect rule for `Int` but [subtly] broken for `Float64`. The kind of `x` for which the rule holds, that kind should have a name. Maybe the name is _“`x - x == 0` holds”_, but I suspect there’s a better name.

And then `q` is not that kind of thing. It is a different thing, and as such, the rule `x - x == 0` does not apply to it.

The types here are geared toward safety / soundness. That’s in contrast to the types in Julia, which arguably exist _primarily_ to direct dispatch. I’m sure it’s not easy to figure out those types, but it does seem crucial to think in those terms.

---

<div class="post-metadata">

**Author:** ![rfateman](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rfateman/32/22507_2.png) [@rfateman](https://discourse.julialang.org/u/rfateman)\
**Post date:** [April 16, 2021, 6:41pm UTC](https://discourse.julialang.org/t/cas-best-practices/58092/98 "2021-04-16T18:41:59Z")

</div>

If you wish to constrain what the computer algebra system can compute by making everything type-safe, then it seems to me that you might declare that q is type “real” but then  
q:=limit(f(x), x-\>0) is type unsafe since f(x)=1/x returns inf, and that’s no good.

So q is perhap of type “algebraic expression”. Indeed, that covers Float64, Int, BigInt, arithmetic expressions like x+y, etc.  
I’m not saying that I have a recipe for fixing all this. I think experience has demonstrated that “types” at least in the conventional programming language situation, is not the solution. As you say, it’s not easy to figure out.

[Previous page](https://discourse.julialang.org/t/cas-best-practices/58092.md?page=4)