Generic functions with irrationals that work with Symbolics / ModelingToolkit

General Usage
, ,
1

I have some functions that basically evaluate polynomials with coefficients given by enormous expressions. I want to support multiple precisions — like Float64 and Double64 — but also symbolic operations. And some of the coefficients involve irrationals. Let’s just take \pi as the example, though I actually need quite a few irrationals that I define. To support the various precisions, I simply put the polynomial inside a let block, at the start of which I redefine the irrationals (and some other things I need). Here’s a simplified version.

function p(x)
    let π=oftype(x, π)
        1 + (2//3 + 4π/5) * x
    end
end

This seems to work great for the floating-point types. But then I want to do symbolic stuff, like get the expression for p when x is a Num. In the let block, π is defined to be a Num type, which seems great, until I actually do anything with it — like multiply, divide, add, subtract. In those cases, π seems to be immediately converted to Float64, which is not what I want. This also seems to be a problem when I pass the result to ModelingToolkit.

Is there some way to hold off on evaluating expressions involving Num(π), somewhat like @register?

2

Open an issue. I don’t think anyone has looked at interactions with Irrationals?

3

New issue is here.

Reading through some of the linked posts, I found a hack that seems to work for me, for now: using things like Term(identity, π).

For the code above, the simplest (and extra hacky) approach is to just add a new method to oftype:

Base.oftype(x::Num, y::Irrational) = Term(identity, y)

With that, p(x) returns the appropriate type, and retains all appropriate digits of \pi for floating-point input, but leaves the term unevaluated for symbolic input.

Is there a downside to this approach?

1 Like
4

That looks correct to me, but I’d have @shashi review it. Make a PR?

5

The oftype code above works for me, but wouldn’t work for most cases, so I’m not sure what the PR should be. A slightly more general/useful version might be

Base.convert(::Type{Num}, y::Irrational) = Symbolics.Term(identity, y)

But I think that’s still not what most people would want. From the discussion here (and maybe here), it looks like people want to be able to do something like Num(π) and get the sensible thing, but surely it would be too much to do Num(y::Irrational) = Term(identity, y), since that wouldn’t actually result in a Num, right?

I’m happy to do a PR if it helps move in a positive direction, but I feel like this problem calls for more thought and more knowledge of the inner workings than I have.

[But if Term(identity, y) is actually used, it might be worth adding a special case to show_term to just show the argument if the operation is identity, because \mathrm{identity}(\pi) is pretty ugly and could be confusing.]