This doesn’t work if x is an array — zero(x) * x will throw if x is not a square matrix. In the context of rings it is fine, of course, though it still won’t be equivalent to x-x for square matrices that may contain non-finite values.
I guess you could do something like zeroscalar(x) * x where
zeroscalar(x) = zero(x)
zeroscalar(::Type{<:AbstractArray{T}}) where {T} = zeroscalar(T)
zeroscalar(x::AbstractArray) = zeroscalar(eltype(x))
julia> using Symbolics
julia> Term(sqrt,[2])
sqrt(2)
I don’t get why you’d work harder than that. You can build any expression to be lazily like that, or use @register on others.
I thought the reason for it was because the poster wanted “a zero in the type of x”:
In which case trying to get 0 * x as they posted is an X-Y problem: just writing it as zero(x) is more direct and means what they actually meant using Base Julia functionality.
I think you’re right, I was mistaken. Term(-, [x,x]) produces x - x which seems like what’s desired. Why doesn’t x - x produce the same result? For @syms x::Any, x-x errors, and for @syms x::Num, x-x returns 0::Float64.
0 * x is not the same as zero(x) if x is not finite (NaN or ±Inf). That’s why x - x can be simplified to 0*x (for an appropriate 0) but not to zero(x) in general, and simplifying x-x is the context here.
As the originator of this thread, it is my hope and request that we ( that means you Julians 
) return to being about “CAS Best Practices”. That there are different considerations held by CAS experts is helpful information. Let it be, so we move forward. Thank you all for participating.
If you want it exact, don’t ask for the Float64 result?
julia> ex = Term(sqrt,[big(2)])
sqrt(2)
julia> eval(toexpr(ex))
1.414213562373095048801688724209698078569671875376948073176679737990732478462102
It just follows Julia’s semantics.
I see. Yeah there’s a more interesting question to be answered here, which is what kind of number is the default symbol? The default is @variables x::Real, specifically not ::AbstractFloat. Since these rules are all type-based, I think simplifying x-x -> 0 on ::Real is the semantics you’d expect, while x-x -> 0 * x should be what happens on ::AbstractFloat. This keeps two algebras distinct: real numbers obey commutative and associative laws, they are a field, etc. while floating point numbers have a different set of rules.
The reason for keeping this difference is that the set of transformations that keep exact floating point semantics is already in LLVM, so I don’t think there’s too much to gain there. But it would be a good opt-in: “make simplifications only do what’s allowed on floating point numbers” instead of “make simplifications only do what’s allowed on real numbers”. The latter seems to be what people normally want on models, or at least out of a CAS, though there are times when exact floating point semantics need to be followed (in which case we have @register for someone to register a function and keep some calculation intact). What I want to do is add a call to Herbie, so that in many cases we aren’t just applying rules that change floating point behavior but actually improve it, so even further down the path of “keep floating point separate”.
This makes sense on its own, but I do not know how to square it with AbstractFloat <: Real. This is a question related to (or exactly about?) co/contra/invariance, and how the subtyping relation shows up, if at all, in figuring out what rules apply when. Is the following a bug? Or is it a necessary consequence of the rule defined on Real?
julia> y = SymbolicUtils.Sym{AbstractFloat}(:y)
y
julia> y - y
0
I would say yes, it’s a missing specialization.
You have subtypes which override defaults of their supertype all of the time. AbstractArray for example has many array types which change underlying behaviors of broadcast, indexing, etc.
But subtypes should not override the defaults (or “fallbacks”) of their supertype “willy nilly”. There should be some kind of consistency with the fallback. e.g.:
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.
That’s a good question. I would assume that for number types you’d want to have those simplifications by default, and choose on a subset of AbstractFloats to disable them. I could see someone going the other way and having @variables x default to SymbolicReal <: Real and define those rules on SymbolicReal. It’s worth a discussion on the SymbolicUtils.jl repo. Open an issue.
There is no problem whatsoever with defining the square root function (sqrt or \sqrt{}) to be the principal square root. This is the standard approach, as described in the Wikipedia article on square roots and as implemented in Julia, Python, Mathematica, Maple, Maxima, etc. Insisting that we define \sqrt{9} to be \pm3 is a pathway to confusion.
I would think that if you want to show that \sqrt{x^2} and x are unequal as expressions, you’d want to point out the example of x = -3. That way the argument does not depend on a non-standard interpretation of the square root symbol.
From a conventional programming language point-of-view, yes,
From a logico-symbolic point-of-view where
y = sqrt(x) means that x == y^2 , y = sqrt(9) means that 9 == y^2
asking Maple to solve(y^2 = 9, y) gives (3, -3)
A good deal of the back-and-forth is frisson between Numerics and Symbolics (both are math).
We are working toward a best (or a very good) case implementation approach that allows us to glean important benefits that each perspective brings, sometimes one better than the other and vice versa.
My point is exactly that the square root function is, in all common usage, distinct from the solve function. If we insist on conflating the two, how do we interpret an expression like \sqrt{2} + \sqrt{3}? It’s got to be a four-element set, I suppose. Here we’re already verging on chaos, wishing we had a concise way to say “the positive square root” when we want to.
I’m comfortable granting the point that there are many interesting loci of interplay between the demands of symbolic and numeric computation, some of which may require creative thinking about how common concepts need to be treated differently. But I remain skeptical about the specific proposal of abandoning the square root function (which is implied by an assertion like \sqrt{9} = \pm 3).
We are not going to abandon any of these print(names(Base.Math)). I understand how that impression may have been conveyed, inadvertently. What is being considered is how to respect the numerical math functions without abandoning the symbolic/theoretical math functions that may carry the same names in ordinary discourse (hence the tangling of notions). It is my guess that we will find a way to distinguish the desired context (symbolic, numeric, some of each) for any occurrence of the named function – there are multiple avenues for this. The simplest way to distinguish the Julia numerical square_root from the more constructively mathematical symbolic square_root would be use the name from Base.Math as the numerical function and to use its titlecase spelling as the symbolic function sqrt(2) !== Sqrt(2). That leaves out a mixed version (maybe Sqrt is mixed, and purely symbolic is a context, as with let, and do blocks).
see next post by @DNF for one of the more familiar approaches within Julia
Doesn’t that run a bit counter to the general Julia philosophy? Could one instead leverage multiple dispatch with sqrt(2) vs sqrt(Symbolic(2)) (just making up names here)?
Disclaimer: I know nothing about CASs in general, and hardly anything about Symbolics.jl, so this is a question rather than a suggestion.
of course. I just wanted to clarify that these two things (the numerical evaluation and the symbolic handling) are not required to be one thing (although synthesis is good if it work).
The principal square root function is a logically, mathematically, symbolically and numerically perfectly fine function. Using the names sqrt/Sqrt to denote anything else is a mistake. Not because other notions of square root aren’t useful too, but because this convention is already firmly established. If you want a function that returns a set, call it something else (all_square_roots(x) or whatever).
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.
This confusion essentially goes away if you build symbolic sets with well-typed quantifiers into the system. If you want to solve x^2 = y, write down S = \{x : x \in \mathbb{C} \land x^2 = y \} or S = \{x : x \in \mathbb{R} \land x \ge 0 \land x^2 = y \} or something else, depending on what you mean. Now various notions of solving can be reduced to computations on the set S: checking its cardinality, extracting an arbitrary element, iterating over the elements, computing some superset, taking the intersection with another set, etc. To keep things convenient for users, you may want unique_solution(), any_solution() and all_solutions() functions or something similar.
There is actually nothing inherently wrong with simplifying sqrt(x^2) → x or solving x^2 = y → x = sqrt(x) as purely formal operations on formal symbolic expressions. You just have to remember that these are formal operations which don’t commute with evaluation for x \in \mathbb{R}, say. When we do mathematics by hand, we often ignore conditions and check that the result makes sense in the end. The traditional CAS design approach mimics human operation but leaves the final checking to the user. This is clearly very useful in many cases, but it makes it seriously difficult to write correct complex programs that compose operations (where the user can’t realistically check each step). I’m convinced that it’s possible to do better; there’s no need to throw out the classical heuristic symbolic manipulation completely, but it should be a complement to and not a substitute for correct domain-aware rewriting. There are some parallels with floating-point versus interval arithmetic here (including the fact that floating-point solutions much like heuristic symbolic solutions sometimes can be checked a posteriori!).
I also completely agree with Oscar Benjamin that automatic evaluation is a mistake. Simple, unevaluated expressions by default are a win-win: easier to implement, easier and more versatile to use (the user can decide precisely which transformations to apply), and more efficient (symbolic expressions are a terrible data structure for algebraic computation). Instead of performing clever and invariably buggy automatic “canonicalization” on expressions, a far better solution is to convert to specialized data structures internally for computations (e.g. polynomial types for expanding or factoring).
Yes, that’s why in the tracer it’s planned to not do automatic evaluation, as described in Symbolic sqrt(2) and sin(pi) - #21 by ChrisRackauckas, with dispatches for alternative constructor-based simplifications all being type based.
That’s a good idea.
That is what we do. Sounds like Symbolics.jl is on the right path then, thanks.