Request for some numerics related features

#1

I’ve been experimenting with toy implementations of popular libraries that rely on “template magic” and/or operator overloading, to test the boundaries of Rust’s capabilities.

An interesting set of such libraries are the various “extended numbers”. Examples:

  • Dual numbers, which can be used to implement Automatic Differentiation
  • Interval Arithmetic, which can be used to implement more robust implementations of many numeric algorithms.
  • Probability Distributions, which are used for Bayesian inference, often using the Markov Chain Monte Carlo algorithm.

Using Traits and Templates, these can be combined in interesting ways, such as Dual numbers of Intervals, which can be used to differentiate functions while simultaneously reporting the accumulated rounding errors incurred by the computation.

What I’ve found to be missing:

  • Rounding modes – The “round()” functions are insufficient to implement Interval Arithmetic. There is a need to set the rounding mode for all computations, not just casts. This is in the standard C and C++ APIs.
  • Trait Constants – the only available constant via a Trait is Default, which is typically the numeric value “zero”, or the equivalent. There is no elegant way to specify a Trait for “a type with a multiplicative identity” (one) or the various inverses. The closest is From<i32>, which allows the use of From::from( 1 ). That may not make sense for some types, such as matrices or tensors.
  • Integers in Templates – I know this has been requested by many people, but I’d like to throw my hat in the ring as well. Dual numbers for example are typically used with multiple variables, and could benefit from small constant sized arrays specified via template arguments.

While I recognize that Rust is targetted at “systems” and “web” programming, not scientific mathematics, the extended numbers described above are highly relevant to many common applications, such as machine learning, image filtering, robotics and sensor data processing, etc…

Interval Arithmetic in particular is a good match for a “robust” language like Rust!

#2

I looked into this very briefly a while back, but couldn’t find any concrete, well-defined support for it in LLVM (not an expert, could be wrong). It seems to be something that LLVM, at least, sort of ignores; you can go behind its back to change the rounding modes, but what happens when you do is ambiguous. I do recall that D had a period where this “worked”, but was entirely unspecified until it became an aspect of the ABI.

RFC issue #1062: Constants that depend on type parameters in generic code.

Relevant note in RFC #195.

There are associated constants (an unstable, nightly-only feature), but my understanding is that they greatly complicate things once generics get involved. So this looks to be a case of “in progress and wanted, but could be a while”.

RFC issue #1038: Parameterize types over numerics / type level integers.

Yes, this is wanted, but doesn’t seem as easy as it sounds.

#3

I’m pretty keen on Rust for scientific programming, and there’s certainly interest in this space from the team (e.g. I spent the summer doing an internship on Rust, working on SIMD).

Also, FWIW, Rust isn’t really targeted at ‘web’ programming: there is Servo as a driving force, but the demands of a rendering engine are closer to systems programming than to other sorts of web programming. It’s mostly a happy accident that a ton of people are using Rust in this space/developing libraries to allow Rust to work for servers etc… I guess it’s a reflection of the modern world & what a lot of vocal early adopters like doing. :slight_smile:

You’re (unfortunately) not, e.g. LLVM mishandles constant evaluation in the presence of fesetround. Although, ignoring that rather large hole, FFI calls to fenv.h functions should work fine, e.g.:

extern {
    fn fesetround(x: i32);
}

fn main() {
    unsafe {
        let x = 1.23_f32;
        // (implicitly) to nearest
        let a = x.sqrt();

        fesetround(0x400); // to negative infinity
        let b = x.sqrt();

        fesetround(0x800); // to positive infinity
        let c = x.sqrt();

        println!("{}\n{}\n{}", a, b, c);
    }
}

1.1090536
1.1090536
1.1090537

This example also demonstrates the constant evaluation problems with rustc’s use of LLVM: compiling with optimisations makes c have value 1.090536 like the others, since the x.sqrt() is constant folded with round-to-nearest.

I think @peter_bertok (correct me if I’m wrong) is talking more about traits like Zero and One than associated constants.

On that point, there’ve been various experiments with numeric trait towers in std in the past, but none really stuck.

#4

I’m interested where that sentiment comes from. I never had the impression that Rust cares a lot about web programming that much.

#5

…the demands of a rendering engine are closer to systems programming than to other sorts of web programming…

That’s actually a good point: a common issue with graphics rendering engines is that floating point rounding can create 1-pixel cracks between shapes that should otherwise form a continuous area of colour, resulting in ugly visual artefacts.

If I remember correctly, one of the driving forces for the adoption of IEEE 754 floating point arithmetic in GPUs was to ensure at least some level of consistency to help developers avoid issues exactly like this.

Also, your point about constant folding is a good one! If the language doesn’t have built-in rounding mode support, FFI calls to set the rounding mode will leave too many pits behind to fall into.

Performance is another issue. One of the C++ Interval Arithmetic library developers mentioned that switching rounding directions was expensive, but the library would (naively) require switching the mode at least twice for just about every basic operation! His solution was to set the entire program to round in one direction consistently, then he used what amounts to “double negations” to efficiently simulate rounding in the other direction. However, that kind of code is likely to be optimised away by the compiler!

In essence, FP arithmetic needs a set of special sections vaguely similar to the C# unchecked { } block:

fn foo() -> f64 {
    // The compiler is free to replace this constant with '1.0'
    let x = 0.5 + 0.5;
    // Set rounding mode and prevent the compiler from
    // constant folding with any other mode by strictly
    // following IEEE 754 rules. It can still do folding,
    // but only rounding away from zero instead of whatever
    // the default is.
    strict round(Round::away_from_zero) {
        x / 10.0
    }
}
#6

Just some additional context: this is something I recall D having to deal with: the ABI specifies that the rounding mode must be restored before leaving a function, or all floating point operations beyond that point are undefined. I know Walter Bright (and especially Don Clugston) care about floating point, so my assumption has always been that “fast”, “correct”, and “rounding control” are a “pick two” situation: if there was an obvious solution, I expect D would have adopted it.

From what I recall, I don’t believe block control would really work; I mean, changing the rounding mode fundamentally changes what the code means. It’d be like a viral, implicit generic parameter on all code that touches any FP anywhere. Or, a bit like a purity annotation.

Again, not an expert, just trying to recall possibly salient details from barely remembered forum threads. :slight_smile:

#7

I was initially thinking that this could be vaguely solved with something like variants of the f32/f64 types, such as f64_round_up or something, but thinking about the problem a little more made me realise that it’s more akin to collations in SQL, where it’s the operators that have modifiers, not the types. It makes no sense to multiply an f64_round_up type with an f64_round_down type, but a mul<round_up> operator makes sense, much like a sort using a particular collation.

However, it sort-of makes sense to allow the rounding direction to be passed as a template parameter, because there are specialisations available for different rounding modes. E.g.: if the rounding mode is already correct, perform the computation directly. If it’s the exact reverse of what is required, perform a double-negation of the computation. If it’s any other mode, then set the required mode for the duration and then undo.

#8

This, Posit Arithmetic John L. Gustafson 21 February 2017, may be helpful.
I’d really like to do it the Rust way.
Not unum 1 or 2 … and definitely not IEEE 754. Very recent insight by John; posit.

Beating Floating Point at its Own Game: Posit Arithmetic
Stanford Seminar: Beyond Floating Point: Next Generation Computer Arithmetic

Isaac Yonemoto is doing interesting work with juila lang and posit arithmetic.
and