sharkdp / numbat

A statically typed programming language for scientific computations with first class support for physical dimensions and units
https://numbat.dev
Apache License 2.0
1.3k stars 53 forks source link

unit_of is unsound #521

Open sharkdp opened 4 months ago

sharkdp commented 4 months ago

The unit_of function

fn unit_of<T: Dim>(x: T) -> T

which returns e.g. 1 km/h for an input of 30 km/h can lead to soundness bugs when the argument is zero.

Consider this case, for example:

>>> fn normalize_length(x: Length) = unit_of(x)
>>> normalize_length(0) + 1 m
error: runtime error
 = Conversion error: unit 'm' can not be converted to ''

The type checker doesn't complain since the signature of unit_of promises to return a Length for this case. But unit_of(0) == 1, which leads to the subsequent conversion error.

It can even lead to crashes, like in this example:

fn diff(f, x) = f(x + unit_of(x)) / unit_of(x)
fn amp(t) = cos(1 Hz × t)

diff(amp, 0)
rben01 commented 2 months ago

Just to have an MRE, cos(1Hz * unit_of(0)) is sufficient to crash numbat. Also, an issue that is surely related is the following:

>>> fn f(x: Length) = x  # ok

  fn f(x: Length) -> Length = x

>>> f(1Hz * 1)  # ok
error: while type checking
  ┌─ <input:18>:1:6
  │
1 │ fn f(x: Length) = x
  │      ^ Length
  │
  ┌─ <input:19>:1:3
  │
1 │ f(1Hz * 1)
  │ - ^^^^^^^ Activity or Frequency
  │ │  
  │ incompatible dimensions in argument 1 of function call to 'f'
  │
  = parameter type: Length
     argument type: Time⁻¹    [= Activity, Frequency]

>>> f(1Hz * 0)  # ok

  f(1 hertz × 0)

    = 0    [Length]

>>> unit_of(0) == 1  # ok

  unit_of(0) == 1

    = true    [Bool]

>>> f(1Hz * unit_of(0))  # huh??

  f(1 hertz × unit_of(0))

    = 1 Hz    [Length]  # ???
rben01 commented 2 months ago

Also it appears only to affect the literal 0, not something like 1-1. Neither cos(unit_of((1-1)m)*Hz) nor cos(unit_of(1-1)*Hz) type check, but neither crashes Numbat. Could be something to do with the constraint solver treating literal 0 differently?

If you could reach into a unit_of call and replace any literal 0s with 1-1, that would be an interesting hack.

rben01 commented 2 months ago

On second thought is the underlying issue that unit_of is unsound, or that handling of a literal zero obliterating its unit is special cased in the type checker? During elaboration (turning an act::Statement into a typed_ast::Statement), special care is taken to make sure scalars equal to 0 become polymorphic. But this is not extended to values equal to 0, such as 1-1. Therefore you end up being able to compute cos(0m) but not cos((1-1)m, as the second’s zero value is not polymorphic.

The solution probably has to be applied at a lower level than elaboration — something at runtime that makes a unitful zero coercible to any other unit. Then at this stage unit_of can just read off the unit currently attached to its zero. Of course, this would drastically change how unit_of is implemented, as it could no longer be an ordinary function; it would have to special-cased into the runtime.

sharkdp commented 2 months ago

Also it appears only to affect the literal 0, not something like 1-1.

scalars equal to 0 become polymorphic. But this is not extended to values equal to 0, such as 1-1

This is the intended behavior. There is no way we could infer a polymorphic type for arbitrary expressions that evaluate to zero. Just think of something like rand_int(0, 1). And even for expressions where we could do it (like the 1 - 1 example which could be constant-folded to 0), we don't want to infer a polymorphic type.

Why do we want a polymorphic type for 0 anyway? Because it's very convenient for things like this:

fn failing_sqrt<D: Dim>(x: D^2) -> D =
  if x >= 0
    then sqrt(x)
    else error("Argument to failing_sqrt is negative")

This generic function can be used with dimensionful arguments like failing_sqrt(2 m²). In order to be able to compare arbitrary dimensionful quantities to zero, it needs to be polymorphic. Otherwise, we would need to resort to a comparison like value_of(x) >= 0 or x >= 0 · unit_of(x), which is unexpected and inconvenient.

On the other hand, even if 1 - 1 is mathematically equivalent to 0, it's clearly an expression involving plain numbers, that just happens to have a value of zero. It should have a type of Scalar. And an expression like f(0) where f(t) = 30 cm × sin(20 Hz × t) should have a type of Length, even if f(0) = 0 m = 0. Because it's clearly a calculation involving a length — and not a polymorphic expression.

I understand that this is a strange concept, but I am quite confident that this is the right behavior. It's very similar to how the empty list [] is typed. The plain expression [] has a polymorphic type forall A. List<A>. But something like take(0, ["foo", "bar"]) has a concrete type List<String> even though it evaluates to []. This is not something special to Numbat. You can also see this in other languages like Haskell:

ghci> let x = []
ghci> let y = take 0 ["foo", "bar"]
ghci> :t x
x :: [a]
ghci> :t y
y :: [String]
ghci> x
[]
ghci> y
[]

The solution probably has to be applied at a lower level than elaboration — something at runtime that makes a unitful zero coercible to any other unit. Then at this stage unit_of can just read off the unit currently attached to its zero. Of course, this would drastically change how unit_of is implemented, as it could no longer be an ordinary function; it would have to special-cased into the runtime.

Yes. I think the "unit attached to zero" thing might work. I need to think about this. The ideal solution for me would be to get rid of unit_of and value_of. Even if we can fix this somehow, these functions call for trouble. I need to look into the existing use cases of these functions. Maybe there is a better way.

rben01 commented 2 months ago

After tinkering around a bit, I agree with you that polymorphic (literal) zero is 100% the correct decision. (I tried removing it and was met with immediate pain.) I guess my question is where it should be applied: permanently during elaboration, or temporarily everywhere it's needed to make the types check? With the second option, 0 would be polymorphic during type checking and evaluation, but would live in the VM itself as a concrete unitful value. This seems like it may be overkill though if all it does is enable a sound unit_of.

sharkdp commented 2 months ago

I think the issue goes deeper than this. The mathematical equivalence 0 m == 0 should always hold true. So we need to be able to (actually or mentally) replace 0 m with 0 in a Numbat program and get the same result.

Even if we change the runtime representation of 0 m to somehow retain the m unit and make unit_of(0 m) = m work somehow, I don't see a way how we can uphold this property above. How could unit_of(0) possibly return m where every other unit is an equally possible answer?

I quickly thought about run-time type information or something similar. If we have a "problematic" function like

fn f(x: Length) = unit_of(x)

in the program, we do infer a type of Length -> Length for the unit_of call. If we proceed to call this function with zero, we actually have something like unit_of<Length>(0). And I first thought that we could use this fact somehow by returning some canonical unit of type Length (e.g. the base unit meter). This would result in sound behavior (in a type-level sense). But would still be problematic in practical applications because unit_of(0.001 foot) would be foot, but unit_of(0.000 ft) would be meter. So there would be a discontinuity at zero.


I think it's useful to think about alternatives to unit_of. One example where we need it is numerical differentiation. The scalar version of this has no problem:

fn diff(f: Fn[(Scalar) -> Scalar], x: Scalar) -> Scalar =
  (f(x + Δx) - f(x - Δx)) / 2 Δx
  where
    Δx = 1e-10

But the fully polymorphic version somehow needs to create a Δx of the right unit/type.

fn diff<X: Dim, Y: Dim>(f: Fn[(X) -> Y], x: X) -> Y / X =
  (f(x + Δx) - f(x - Δx)) / 2 Δx
  where
    Δx = 1e-10 × unit_of(x)

One simple (and not totally unreasonable) way to fix this would be to ask the user to provide Δx:

fn diff<X: Dim, Y: Dim>(f: Fn[(X) -> Y], x: X, Δx: X) -> Y / X =
  (f(x + Δx) - f(x - Δx)) / 2 Δx

This would also get rid of the "magic" 1e-10 value that can never be correct for all applications.

But we need to look at other instances of unit_of/value_of to see if it's reasonable to get rid of them.

For now, I'll merge a hotfix that patches the crash: #615