[Symbolics.jl] Defining √2 as an (unevaluated) symbol

General Usage
1

Hello,

with a colleague of mine, we are trying to evaluate Symbolics.jl (vs. sympy/SymPy.jl vs. maxima vs. maple).

A feature that we would really need would be to define √2 as a symbol, unevaluated. What we have in mind is something like

In [1]: import sympy

In [2]: sympy.init_printing(use_unicode=True)

In [3]: sqrt_2 = sympy.sqrt(2)

In [4]: sqrt_2
Out[4]: √2

in python/sympy, or sqrt(Sym(2)) in julia/SymPy.

Is it possible with Symbolics.jl? I am under the impression that Symbolics.jl was maybe not designed with such applications in mind.

I am aware of the post Symbolic sqrt(2) and sin(pi). However, the solution proposed is not satisfactory for our use cases. Indeed

julia> using Symbolics

julia> sqrt_2 = Symbolics.Term(sqrt, [2])
sqrt(2)

julia> sqrt_2^2
sqrt(2)^2

What we really want is sqrt_2^2 to return 2 (as SymPy would). Is this feature already existing? If not, is it planned, or completely out of scope?

Thanks for your answers!

4 Likes
2

I don’t know why this rule isn’t part of the standard set (probably because of some floating point shenanigans with complex numbers), but if you truly desire to simplify like that, here’s how to do that. You just have to add that rewriting rule and make it available to simplify:

julia> using Symbolics, SymbolicUtils  
                                       
julia> r = @rule sqrt(~x)^2 => ~x      
sqrt(~x) ^ 2 => ~x                     
                                       
julia> sqrt_2 = Symbolics.Term(sqrt, 2)
sqrt(2)                                
                                       
julia> simplify(sqrt_2^2, RuleSet([r]))
2                                      

# also works with other numbers here, though again,
# in real applications beware of floating point shenanigans
julia> sqrt_2 = Symbolics.Term(sqrt, 1325)
sqrt(1325)                                
                                          
julia> simplify(sqrt_2^2, RuleSet([r]))   
1325

I encourage you to check out the excellent documentation of SymbolicUtils and how those rules work (and how you can define your own!).

2 Likes
3

Another option would be checking out Metatheory.jl (check out a talk introducing it from this years JuliaCon here), which is basically very fast term rewriting/simplification (also has support for Symbolics.jl). Not 100% sure whether they have that rewriting rule defined by default though.

4

Is there a way to make an expression like \sqrt{x} and e.g. substitute 2 for x and get the exact value \sqrt{2}?

julia> using Symbolics, SymbolicUtils

julia> @variables x
1-element Vector{Num}:
 x

julia> e = sqrt(x)
sqrt(x)

julia> substitute(e, (Dict(x => 2)))
1.4142135623730951
2 Likes
5

That would be awesome. Have there been discussions about symbolic computations like \sqrt x|_{x=2}\rightarrow \sqrt 2, @YingboMa @shashi?

1 Like
6

There could probably be a keyword argument for “substitute but don’t interpret”. Can you open an issue for that?

2 Likes
7

No need to submit an issue. The feature already exists.

julia> using Symbolics

julia> @variables x
1-element Vector{Num}:
 x

julia> substitute(sqrt(x), Dict(x=>2), fold=false)
sqrt(2)
3 Likes
8

TBF, undocumented.

9

This is quite exciting indeed. I will watch this video very closely, thanks for the suggestion.

10

I encourage you to check out the excellent documentation of SymbolicUtils and how those rules work (and how you can define your own!).

This indeed seems to solve my problem. I will study SymbolcUtils more closely, then. Best

11

Hello,

First, let me thank everyone which posted on this thread. Just for the record, here is a slightly more complex example (I agree it is quite trivial in this present form, but it might apply to more complex situations)

julia> using Symbolic
SymbolicUtils Symbolics
julia> using SymbolicUtils, Symbolics

julia> r = @rule sqrt((~x)^2) => abs(~x)
sqrt((~x) ^ 2) => abs(~x)

julia> @variables x
1-element Vector{Num}:
 x

julia> sqr_minus_2 = substitute(x^2, Dict(x=>-2), fold=false)
(-2)^2

julia> simplify(sqrt(sqr_minus_2), RuleSet([r]))
2
1 Like
12

I’m not an expert in terms of Symbolics.jl or CAS, but from my understanding, substitute should always be the last step, no? I.e. first simplify the expression given as far as possible, under all known Rules/legal transforms, substitute a variable and continue with evaluation.

As far as I understand from your questions, you basically want to avoid seeing “ugly” results/results containing literal floating point values?

13

You’re right. substitute is, as far as I understand, a “trick” to define sqrt(2) (or in the present case, (-2)^2) as an unevaluated symbol.
Note that we need to be able to define sqrt(2) as a symbol, because we use a matrix representation of second- and fourth- rank, tensors with minor symmetries which use the sqrt(2) weighting.
And we really need that sqrt(2)² should simplify as 2.

14

You mean something like this?

julia> f = [ sqrt(x)^(y*(i+j)) for i in 1:5, j in 1:5 ]               
5×5 Matrix{Num}:                                                      
 sqrt(x)^(2y)  sqrt(x)^(3y)  sqrt(x)^(4y)  sqrt(x)^(5y)   sqrt(x)^(6y)
 sqrt(x)^(3y)  sqrt(x)^(4y)  sqrt(x)^(5y)  sqrt(x)^(6y)   sqrt(x)^(7y)
 sqrt(x)^(4y)  sqrt(x)^(5y)  sqrt(x)^(6y)  sqrt(x)^(7y)   sqrt(x)^(8y)
 sqrt(x)^(5y)  sqrt(x)^(6y)  sqrt(x)^(7y)  sqrt(x)^(8y)   sqrt(x)^(9y)
 sqrt(x)^(6y)  sqrt(x)^(7y)  sqrt(x)^(8y)  sqrt(x)^(9y)  sqrt(x)^(10y)
                                                                      
julia> r_3 = @rule ^(sqrt(~x),~y) => ^(~x, ~y/2)                      
sqrt(~x) ^ ~y => (~x) ^ (~y / 2)                                      
                                                                      
julia> simplify.(f, Ref(RuleSet([r_3])))                              
5×5 Matrix{Num}:                                                      
          x^y  x^((3//2)*y)  x^((2//1)*y)  x^((5//2)*y)  x^((3//1)*y) 
 x^((3//2)*y)  x^((2//1)*y)  x^((5//2)*y)  x^((3//1)*y)  x^((7//2)*y) 
 x^((2//1)*y)  x^((5//2)*y)  x^((3//1)*y)  x^((7//2)*y)  x^((4//1)*y) 
 x^((5//2)*y)  x^((3//1)*y)  x^((7//2)*y)  x^((4//1)*y)  x^((9//2)*y) 
 x^((3//1)*y)  x^((7//2)*y)  x^((4//1)*y)  x^((9//2)*y)  x^((5//1)*y)

May not be exactly what you’re looking for, since I have a Matrix of symbolic variables that I simplify via broadcasting. Afterwards I can just substitute:

# just y
julia> substitute.(res, Ref(Dict(y => 1)))                           
5×5 Matrix{Num}:                                                     
        x  x^(3//2)  x^(2//1)  x^(5//2)  x^(3//1)                    
 x^(3//2)  x^(2//1)  x^(5//2)  x^(3//1)  x^(7//2)                    
 x^(2//1)  x^(5//2)  x^(3//1)  x^(7//2)  x^(4//1)                    
 x^(5//2)  x^(3//1)  x^(7//2)  x^(4//1)  x^(9//2)                    
 x^(3//1)  x^(7//2)  x^(4//1)  x^(9//2)  x^(5//1)                    
           
# just x                                                          
julia> substitute.(res, Ref(Dict(x => 2)))                           
5×5 Matrix{Num}:                                                     
          2^y  2^((3//2)*y)  2^((2//1)*y)  2^((5//2)*y)  2^((3//1)*y)
 2^((3//2)*y)  2^((2//1)*y)  2^((5//2)*y)  2^((3//1)*y)  2^((7//2)*y)
 2^((2//1)*y)  2^((5//2)*y)  2^((3//1)*y)  2^((7//2)*y)  2^((4//1)*y)
 2^((5//2)*y)  2^((3//1)*y)  2^((7//2)*y)  2^((4//1)*y)  2^((9//2)*y)
 2^((3//1)*y)  2^((7//2)*y)  2^((4//1)*y)  2^((9//2)*y)  2^((5//1)*y)
           
# both, which causes all symbolics terms to collapse to equivalent values                                                          
julia> substitute.(res, Ref(Dict(x => 2, y => 1)))                   
5×5 Matrix{Num}:                                                     
 2         2.82843   4.0       5.65685   8.0                         
 2.82843   4.0       5.65685   8.0      11.3137                      
 4.0       5.65685   8.0      11.3137   16.0                         
 5.65685   8.0      11.3137   16.0      22.6274                      
 8.0      11.3137   16.0      22.6274   32.0

The last example collapses into direct values since under the given rules (+ default rules), the numerical value is equivalent to the uncomputed symbolic value. The full floating point value is still there by the way, it’s just not printed:

julia> ans[1,2] |> dump        
Num                            
  val: Float64 2.82842712474619

Again, I’m not a CAS expert, I’m just messing around with broadcasting and what I can find in the docs of Symbolics.jl :slight_smile:

15

Hm, I guess my example isn’t quite what you’re thinking of, but it’s a generalized version taking any sqrt(x) weighting you’d like :thinking:

16

How would I go about defining a cleverer rule that can handle powers other than 2?

julia> using Symbolics, SymbolicUtils
  
julia> @variables x, y
2-element Vector{Num}:
 x
 y

julia> expand((1 + sqrt(y))^10)
1 + sqrt(y)^10 + 10sqrt(y) + 45(sqrt(y)^2) + 120(sqrt(y)^3) + 210(sqrt(y)^4) + 252(sqrt(y)^5) + 210(sqrt(y)^6) + 120(sqrt(y)^7) + 45(sqrt(y)^8) + 10(sqrt(y)^9)

julia> expand((1 + sqrt(y))^6)
1 + sqrt(y)^6 + 6sqrt(y) + 15(sqrt(y)^2) + 20(sqrt(y)^3) + 15(sqrt(y)^4) + 6(sqrt(y)^5)

julia> r = @rule sqrt(~x)^2 => ~x
sqrt(~x) ^ 2 => ~x

julia> simplify(expand((1 + sqrt(y))^6), RuleSet([r]))
1 + 15y + sqrt(y)^6 + 6sqrt(y) + 20(sqrt(y)^3) + 15(sqrt(y)^4) + 6(sqrt(y)^5)

It should be possible to reduce this to an expression of the form a + b\sqrt{y} where a and b are integers.

17

My guess would be going the other way around - define

r = @rule sqrt(~x) => (~x)^(1//2)

instead and all those nested powers should collapse nicely:

julia> r = @rule sqrt(~x) => (~x)^(1//2)                                    
sqrt(~x) => (~x) ^ (1 // 2)                                                 
                                                                            
julia> simplify(expand((1+ sqrt(y))^6), RuleSet([r]))                       
1 + y^(3//1) + 15y + 6(y^(1//2)) + 20(y^(3//2)) + 15(y^(2//1)) + 6(y^(5//2))

That rule alone won’t be enough to make the “ugly” powers go away and bring it into the form you think of, but then again I have to stress that I’m not an expert in CAS :slight_smile: I may be missing something with writing rules though and matching those non-variable powers for some pattern and extracting something else from them may be possible, so me not knowing how to do this should not be taking as an indicator that Symbolics.jl can’t do it.

18

I wonder if Metatheory.jl might support that better.

1 Like
19

Actually that’s not true!

I intended for the y to be 2 so it should be an expression in \sqrt{2} rather than \sqrt{y}. In that case all even powers of \sqrt{2} reduce to rational. Something like \sqrt{y}^3 does not though because it would simplify as y\sqrt{y} and you’d be left with polynomials in y and some \sqrt{y} factors.

To correct the above the expression below should be writeable in the form a + b \sqrt{2}:

julia> f = substitute(1 + sqrt(y), Dict(y=>2), fold=false)
1 + sqrt(2)

julia> f
1 + sqrt(2)

julia> expand(f^6)
1 + sqrt(2)^6 + 6sqrt(2) + 15(sqrt(2)^2) + 20(sqrt(2)^3) + 15(sqrt(2)^4) + 6(sqrt(2)^5)

julia> r = @rule sqrt(~x) => (~x)^(1//2)
sqrt(~x) => (~x) ^ (1 // 2)

julia> simplify(expand(f^6), RuleSet([r]))
189.9949493661167 + 1.4142135623730951^6

julia> simplify(expand(f^6), RuleSet([r]), fold=false)
ERROR:

It doesn’t seem like simplify has the fold=false argument…