using Symbolics, LinearAlgebra
@variables q::Real
iden = collect(I(4) .* 1.0)
k = inv(iden .* q^2)
I get true/q^2 on the diagonal instead of 1/q^2
2 Likes
julia> true == 1
true
true is the smallest value that represents the integer 1
julia> true isa Integer
true
1 Like
This looks like a bug that is worth reporting. Even more extreme is
using Symbolics, LinearAlgebra
@variables q
inv(collect(q^2*2.0*I(4)))
giving
4Ă—4 Matrix{Num}:
true / (2.0(q^2)) 0 0 0
0 true / (2.0(q^2)) 0 0
0 0 true / (2.0(q^2)) 0
0 0 0 true / (2.0(q^2))
1 Like
What do you find problematic? What answer would you have expected, taking into account that true == 1?
1 Like
With more complicated expressions it becomes difficult to read 1 instead of true. I tried “simplify” but it doesn’t convert true to 1. Maybe there is an alternative? I am declaring q as Real and would not want a bool result.
There is only a Bool in the symbolic expression, the result will be a floating-point value.
julia> q = rand();
julia> true/2q^2 === 1/2q^2 # Absolutely identical
true
julia> true/2q^2 # Result is a floating-point value
2.8977550256308775
julia> true/2q^2 |> typeof
Float64
The whole point is that true == 1, just like 1.0 == 1 , and it makes sense to use the smallest version of 1 possible to avoid type promotion if not required.
1 Like
You’re right, but I might want to use Symbolics to get a symbolic result without necessarily wanting to calculate it numerically.
In a symbolic expression I would like numbers, not equivalent quantities like Bool. Thank you.
1 Like
You could modify how Bool is printed, but it may have unintended consequences
julia> Base.show(io::IO, b::Bool) = print(io, (b ? "1" : "0"))
julia> true
1
julia> false
0
1 Like
But they are not equivalent. 1.0 === 1 is false and we can get different answers in subsequent computations
julia> x = sum(1.0 for i in range(1,2^16))
65536.0
julia> 2^48*x
1.8446744073709552e19
julia> x = sum(1 for i in range(1,2^16))
65536
julia> 2^48*x
0
1 Like
If you write the prompt more concisely, the result is more readable, too.
julia> using Symbolics, LinearAlgebra
julia> @variables q
1-element Vector{Num}:
q
julia> inv(collect(q^2*2I(4)))
4Ă—4 Matrix{Num}:
1 / (2(q^2)) 0 0 0
0 1 / (2(q^2)) 0 0
0 0 1 / (2(q^2)) 0
0 0 0 1 / (2(q^2))
1 Like
The computations you show here are unrelated to the expression true/2q^2 which is in fact identical to 1/2q^2. If you intend to make further numerical computations, you can always make sure you have a suitable type for the computations you want to perform.
1 Like
To elaborate on why true is a good choice here, consider the following trivial function, if we define it using false, we see that false is adaptive in the sense that multiplying with false will make the result have the same type as the input, irrespective of what the type is:
julia> multiply_by_zero(x) = false*x
multiply_by_zero (generic function with 1 method)
julia> multiply_by_zero(true)
false
julia> multiply_by_zero(1)
0
julia> multiply_by_zero(1f0)
0.0f0
julia> multiply_by_zero(1.0)
0.0
julia> multiply_by_zero(big"1.0")
0.0
Had we defined the function using a floating-point value, we see that the type of the result is different from the type of the input whenever the input type is smaller than Float64
julia> multiply_by_zero(x) = 0.0*x
multiply_by_zero (generic function with 1 method)
julia> multiply_by_zero(true)
0.0
julia> multiply_by_zero(1)
0.0
julia> multiply_by_zero(1f0)
0.0
julia> multiply_by_zero(1.0)
0.0
julia> multiply_by_zero(big"1.0")
0.0
julia>
3 Likes
With a float coefficient that is unfortunately not the case
julia> using Symbolics, LinearAlgebra
julia> @variables q
1-element Vector{Num}:
q
julia> z = inv(collect(q^2*2.0I(4)))
4Ă—4 Matrix{Num}:
true / (2.0(q^2)) 0 0 0
0 true / (2.0(q^2)) 0 0
0 0 true / (2.0(q^2)) 0
0 0 0 true / (2.0(q^2))
Going further with this example worries me
julia> y = z[1,1]
true / (2.0(q^2))
julia> x = sum(y for i in range(1,2^16))
65536(true / (2.0(q^2)))
julia> 2^48*x
0(true / (2.0(q^2)))
1 Like
simplifying or expanding doesn’t seem to help either:
julia> simplify(y)
true / (2.0(q^2))
julia> simplify(x)
65536(true / (2.0(q^2)))
julia> expand(x)
65536 / (2.0(q^2))
julia> 2^48*expand(x)
0
julia> simplify(expand(x))
65536 / (2.0(q^2))
1 Like
This is simply integer overflow
julia> 2^48
281474976710656
julia> 65536 * 2^48
0
it is completely unrelated to Symbolics
1 Like
qed
julia> 2^48
281474976710656
julia> 65536 * 2^48
0
julia> 65536 * BigInt(2)^48
18446744073709551616
1 Like
Yes, that is how I constructed the example. If you reduce your float to an integer then you run in to the overflow problem in subsequent operations. If you keep it is a float, then that problem disappears, at the expense of the answer not being exactly representable. If you change 2.0 to 2.1 or 2.01063072 in the above examples the problems persist, and in the latter cases you would not expect simple replacements with integral or rational constants.
How does one get symbolics to reduce
true/(2.0*q^2)
to a more consistent form such as
0.5 q^(-2)
or
0.5 / q^2
?
2 Likes
I see what you mean; the overflow is hard to diagnose if it accumulates in different parts of the code.
First option I can think of: Use BigInt everywhere; that should be okay in a symbolic computation.
Second option is to implement a float method which recursively turns all numeric coeffiecients into floating-point numbers. see Symbolics.jl/num.jl at master · JuliaSymbolics/Symbolics.jl · GitHub
1 Like
The Num type has builtin choices for its zero element and its unit element given by
julia> using Symbolics
julia> zero(Num)
0
julia> typeof(zero(Num))
Num
julia> one(Num)
1
julia> typeof(one(Num))
Num
We have for example
julia> one(Num)/2.0
0.5
julia> one(Num)/(2.0*q^2)
1 / (2.0(q^2))
So why is this not used?
2 Likes
I would have expected either
0.5 / q^2
or
0.5 q^(-2)
rather than
true / (2.0(q^2))
1 Like