Handling of rewrites with respect to types

[HarrisonGrodin]

How much should we care about the types of equivalent expressions? For example:

• Given a negative Int value of a, should (3 ^ a) ^ a be rewritten to 3 ^ (a^2), even though this would normally throw an error if a isn't converted to a float?

julia> x = -2;

julia> (3 ^ x) ^ x
ERROR: DomainError with -2:
Cannot raise an integer x to a negative power -2.
Make x a float by adding a zero decimal (e.g., 2.0^-2 instead of 2^-2), or write 1/x^2, float(x)^-2, or (x//1)^-2
Stacktrace:
 [1] throw_domerr_powbysq(::Int64, ::Int64) at ./intfuncs.jl:176
 [2] power_by_squaring(::Int64, ::Int64) at ./intfuncs.jl:196
 [3] ^(::Int64, ::Int64) at ./intfuncs.jl:220
 [4] top-level scope at none:0

julia> 3 ^ (x^2)
81

• Should A * 0.0 be rewritten to zero(A), even though the resulting value may be of a different type?

julia> A = [1 2; 3 4];

julia> A * 0.0
2×2 Array{Float64,2}:
 0.0  0.0
 0.0  0.0

julia> zero(A)
2×2 Array{Int64,2}:
 0  0
 0  0

[MasonProtter]

Related to this is your comment on issue #19 with regards to associativity. It's well known that floating point addition is not associative in general.

However, I tend to be of the opinion that the differences between floating point numbers and real numbers are a bug and not a feature so I don't see much point in trying to mimic them. I think that in general, if Rewrite can do better than floating point by taking advantage of known algebraic structure, then we should and the differences that arise should be treated as a shortcoming of floating point arithmetic, not the symbolic manipulations.

In practice, I think my opinion is that floating point addition should be treated as associative, and (3 ^ a) ^ a should be treated as 3^(a^2).

Also, in the case of 3^a I think that we should automatically promote 3 -> float(3) if a < 0.

[HarrisonGrodin]

I would absolutely agree that all imprecision with floating point calculations should be considered more of a bug than a feature.

julia> f(x) = sin(x)^2 + cos(x)^2;  # expected: always 1.0

julia> f(2)
1.0

julia> f(3)
0.9999999999999999

Depending on use cases, though, types may or may not be important to retain. For example, if we run Rewrite as a Cassette @pass and end up rewriting 5.0 - 5.0 to 0::Int, this could lead to confusion, slowdowns, or dispatch failures. On the other hand, if the user is dealing only with symbolic mathematics, I would argue that having any difference between 5 and 5.0 would probably be counterintuitive.

Maybe it's worth having two different (but intersecting) sets of standard rules: one which strictly retains types and one which is more lenient about types but with stronger rules.

---

After thinking about this more, I'd like to build on it even further. It seems like that it could be very helpful (if not beneficial) to have a custom numerical type for storing "infinite-precision" decimals using Rationals. (The usefulness can be seen by the counterintuitivity of having 0.1 + 0.2 != 0.3 alone, from the perspective of purely symbolic mathematics.) Rather than storing 0.12345 as a floating-point decimal, we would translate to the rational 12345 // 100000 and operate on it as such.

Here are a few side-effects of this approach that may be worth considering:

• We would not be able to interact well with functions which return purely floating-point results (e.g. sin). However, this shouldn't be an issue in general; as far as pure symbolic algebra goes, we would never want to evaluate sin(x / y) to a numerical value anyways (except for special cases, like sin(π / 2) => 1 or sin(π / 6) => 1 / 2).
• We would have the potential to display repeating terms in their precise form. (Given ThisNewType(1064//495), we can display 2.14̅9̅.)
• Working with these custom numbers could get cumbersome in standard Julia code. However, this should hopefully be mitigated by custom parsing of numerical types inside of @term expressions.
• It seems to make sense that the type would internally parameterize Rational on BigInt, to avoid confusion with fixed-bit integers and such. However, this isn't set in stone. Should we consider parameterizing the entire type? Is it worth using Complex{BigInt}, even though it has some internal "normalization" logic? (I suppose the same could be said about using Rational in the first place, since it automatically converts 2 // 4 to 1 // 2.)

[MasonProtter]

I like the idea of using rationals everywhere that's possible. That's actually what Mathematica does, come to think of it. We can just do rationals everywhere its possible to do so and then fall back on BigFloat or Float64.

I'm not sure it's a good idea to use Complex unless we're actually in a complex domain, though we should definitely try and support it.

By the way, does the rules system have a way of parameterizing rules on type? It'd solve a lot of problems to be able to do something like

@term RULES [
    sin(x{T})^2 + cos(x{T})^2 where {T<:AbstractFloat} = T(1.0)
    sin(x{T})^2 + cos(x{T})^2 where {T<:AbstractInt} = T(1)
]

[HarrisonGrodin]

Ah, fantastic! Glad there's a precedence for that. I'd rather not fall back to any floating-point arithmetic, since that's accepting an approximation; instead, it seems more pure to just retain the full term (e.g. sin(3)). We can definitely allow for evaluation of terms that converts this new type to BigFloats.

---

The rule system has something close to this feature, but not using the syntax you showed (yet?). It does not allow for usage T in the output pattern. (Note that this example will change slightly with #45.)

x = Variable(:x, TypeSet(AbstractFloat))
y = Variable(:y, TypeSet(Int))

@term RULES [
    sin($x)^2 + cos($x)^2 => 1.0
    sin($y)^2 + cos($y)^2 => 1
end

[MasonProtter]

We should maybe make some sort of rewrite specific type promote function to be applied to heterogeneous types in a term.