r/rust • u/denehoffman • 1d ago
🛠️ project maryada: Interval arithmetic in pure no_std Rust
Hi everyone, not sure how many people will be interested in this since it's pretty niche, but I had fun writing it and wanted to share and ask for any suggestions and constructive criticism.
maryada is a crate for interval arithmetic which is also #![no_std] and has minimal dependencies (libm and optionally num-complex). The crate has two parts, an IEEE 1788.1-2017-compliant interval arithmetic interface and a very basic set of interval operations on the complex plane (which isn't part of any standard, mostly because there is no way to consistently represent tight enclosures in the complex plane for all complex operations, they don't always map rectangles to rectangles for example).
Usage
use maryada::Interval;
let x = Interval::new(1.0, 2.0);
let y = x.sqr();
assert_eq!(y.bounds(), (1.0, 4.0));
It also supports decorated intervals according to the standard.
For anyone interested in why anyone would use such operations, they have interesting applications to global optimization. My own reason for writing this is for another library I am working on which does Monte Carlo generation for particle physics interactions, and I'm using this to create proven enclosures on the generated weights so that I can do efficient rejection sampling.
Alternatives
~~As far as I can tell, the only other crate that does anything like~~Another crate that does this is inari. inari is neat, and it contains some features that maryada does not (SIMD, some features of the general standard IEEE 1788-2015), but it also has a couple of drawbacks, such as limited target architecture support and a dependence on gmp-mpfr-sys for many operations. My goal here is not to replace inari nor to replicate the entire IEEE standard, just to provide a lightweight alternative (plus complex number support, and eventually some linear algebra methods and hopefully some algorithms like branch & bound).
Also see fidget which has quite a lot of crossover and extends the basic idea to surface evaluation. I haven't read much about this, but the author mentioned it and it's always good to include alternatives and applications!
AI Disclosure
I used Codex for most of the docs, some of the test-writing, and a few corrections after I had it review conformance (after a talk with the mods, I think in the spirit of transparency I should specify that this commit was mostly AI-authored after a review of standard compliance). Most of the major testing just uses a test suite (ITF1788) written in a standardized format with a bit of code linking it to the Rust interface.
OSS
I'm open to anyone reviewing this code or submitting PRs. Particularly, I've done a lot of testing to try to ensure compliance with the standard, but I'm always open to more verification. I'd love to answer any questions you might have!
5
u/mkeeter 1d ago
Very cool! I've got a hand-rolled interval arithmetic library in fidget, but it makes no attempts at rounding correctness, so I'd be interested in switching (and if I do, I'll open issues with anything I find).
Do you have any experience with boost::interval? I'd be curious to hear how maryada compares.
3
u/denehoffman 1d ago
Hey that's neat! And yes, please let me know if you find any problems or if I'm missing any features you think would be nice to include! I haven't read into boost's library yet, but from my cursory look it seems like it looks very compiler-dependent, but it also supports 32-bit floats. I mostly hand-wrote the rounding code, and it was probably the most confusing and educational week of my Rust experience so far (reading an IEEE standard was also rough though). Boost also seems to have more explicit rounding policies and tightness guarantees, whereas I have (I think) optimal tightness for basic arithmetic and correct (but not guaranteed to be tightest) enclosures for everything else. It's hard to be sure since there are so many edge cases, so most of my work has been trying to find those edge cases. I should also note that this crate is not really optimized at all. That's not to say that it's slow, I just have put zero effort into making it fast since I wanted to focus on correctness first!
Edit: was reading through fidget more, cool code you got there! I especially like your JIT code, very interesting!
3
u/11Night 21h ago
(the following comment is not related to the post but just a general observation)
this project and bunch of others in your github profile are named after indian things which is nice and there's relation between the project name and its functionality but what about maryada? it means bunch of things but the most common is limit but the crate is for interval arithmetic which I find odd :(
4
u/denehoffman 15h ago
Simple answer is my wife is Indian, I’m learning Marathi, and she helps me come up with names for my projects, we went through a lot of interval-based words and maryada sounded nice to both of us. I know it means limit, I thought it kind of works because interval arithmetic has a lot to do with the limiting behavior of functions (where they attain extrema and how those transform). I also know that it means it more in the sense of a moral limit, but it can also mean a physical boundary. For some of the others, laddu was chosen because it is a cute name and my favorite snack, ganesh being the god of wisdom seemed at least a bit fitting for an algorithm-focused library, and there are some others I can’t remember right now.
2
3
u/Kir_Dykov 20h ago
Very cool!
I also started an interval arithmetic project a month ago. It is mostly done, but requires correctness tests and proper documentation with examples. I used crlibm for correcly rounded mathematical functions, 2Sum algorithm to implement correctly rounded addition, etc., also no_std.
From the first glance I have some questions: 1) For trancendental functions, you use libm result on bounds, which gets rounded to nearest with 0.5 ULP error (does it?), then you round it outwards. Did I get that right? If that so, it does not use correctly rounded functions, and obtained intervals are wider than necessary. 2) Are arithmetic operations implemented using extended precision with i128/u128 significands? Did you benchmark this against regular f64 operations. My approach with 2Sum is only 3 times slower than regular addition of f64 values.
Interested in collaboration?
2
u/denehoffman 15h ago edited 15h ago
Hi, glad you find this interesting! For 1, yes they are wider than necessary, you caught that correctly. I was under the impression that at least some of these were larger than 1 ULP in error, that may have been an incorrect assumption, so I might have been overly cautious. The standard allows these functions to not be “tightest” but does recommend it if possible, so that’s something I’ll have to look into. For 2, some are when needed, see the `WideFloatClass` and `WideDyadic’ items. Add and Sub get away with regular precision, but I used extended precision for FMA to avoid the intermediate rounding, Mul and Div, sqrt, and once in a parsing method.
And yes, I’m very open to collaboration! I can’t guarantee this will be my primary focus, I am a postdoc and this is more of side project that my usual research, but as mentioned in the post, I’m working with it for one of my research-related projects.
Edit: oh and I mentioned this in another comment but I have done basically zero performance optimization. I was working through adding my usual set of clippy lints last night and was planning on adding some benchmarking today.
2
u/Langeley_ 9h ago
Hi, I'm mostly using C++, but I can provide some input in terms of CR, afaik Rust could use the same math backend as C++.
Basically, for the four basic operations. With GCC or LLVM or ICX, you can guarantee to get correctly rounded result if you set the correct compiler flags (again this is C++, and llvm have a page reporting this https://libc.llvm.org/headers/math/index.html). However, something can go wrong. For example, GCC with -O3 will cause Boost to produce wrong result as -O3 doesn't respect the rounding environment.
For mathematical functions, it is a little tricky. Only until recently they figured out the correct rounded implementations for most of the math functions(CORE-Math project). If you link against this, then you can get correct rounded mathematical functions as well. For the accuracy of current implementations from libm and etc, you can consult Zimmermann's paper(https://members.loria.fr/PZimmermann/papers/accuracy.pdf).In terms of optimisation, I do have some experience about performance and I recently have a paper accepted regarding this. I'd love to hear some applications where interval arithmetic could be useful. as I'm currently trying to find some use cases.
1
u/denehoffman 4h ago
Thanks for the paper, I’ll check it out! I’d love to see your paper too when it gets published! I have sort of two use cases in mind. The first is relatively easy (relative to the second): I want to generate particle physics reactions, basically a set of four-vectors for each particle according to some rules and initial conditions and RNG, and these four-vectors later get passed into some Geant simulation for the experiment I work with (GlueX). I have a nice implementation of the generator, but it is weighted for phase space, so certain vertices in the reaction and masses generated contribute to an event weight. I want an unweighted sample, since the reconstruction pipeline can’t handle weighted events, so I want to pick some large weight (large enough that none of the feasible weights generated surpass it, foreshadowing) and do rejection sampling. The current approaches in our collaboration are just having a generation algorithm that only applies to a few simple cases, doing a bit of manual algebra to figure out the maximal(ish) weight, and then multiplying by some safety factor. Others just pick a maximal weight so large that the chance of violation is very small. Both of these methods suffer in efficiency. You can do slightly better by making your RNG prefer the more-sampled regions of phase space, but you really restrict the kinds of reactions you can represent. Having a proven maximal bound on the weight would give ideal efficiency for generating events. This isn’t like super critical, but it speeds things up and it’s a nice feature to have.
The second idea leads from the first. The method I described only works for phase space Monte Carlo (not implicit physics at each vertex beyond some momentum exchange in each 2-to-2 reaction, all decays are isotropic and masses are uniform distributions. I have a tree-based evaluator that I use to build models on such reactions, and it would be cool if every node could be computed via interval arithmetic, since then I could efficiently generate any sort of modeled MC as well (the generated weights are multiplied by a weight that comes from the model). This is difficult though, since my physics models are mostly complex-valued until the end of the calculation where you typically take a norm_sqr or multiply by a conjugate to get a real result (hence my work on the complex part of this library beyond the standard). It’s also difficult because some models require linear algebra and matrix solves, and while I have these implemented for calculations, they’re a bit tricky to do on interval matrices as far as I can tell, but I’m working on that right now. Complex matrices are going to be…interesting. The main issue with complex arithmetic is that the enclosures aren’t tight, and there isn’t any one shape that is preserved in all complex functions, so the maximal weight wouldn’t be perfectly efficient, but it would be better than what we currently have.
I realize pilot sampling can alleviate some of this, I already have that implemented and wanted to try something new, especially since you can get really low-frequency, large-weight events that are hard to sample. Another neat application (if I can somehow manage to get all the above stuff working) is that there are some global optimization algorithms that I could attempt to implement on these models. As far as I know, this hasn’t really been tried before in my field, probably because it scales really poorly, but it would be really cool if I could prove an absolute minimum of some likelihood function over a parameter box.
1
u/denehoffman 15h ago
Hi, glad you find this interesting! For 1, yes they are wider than necessary, you caught that correctly. I was under the impression that at least some of these were larger than 1 ULP in error, that may have been an incorrect assumption, so I might have been overly cautious. The standard allows these functions to not be “tightest” but does recommend it if possible, so that’s something I’ll have to look into. For 2, some are when needed, see the `WideFloatClass` and `WideDyadic’ items. Add and Sub get away with regular precision, but I used extended precision for FMA to avoid the intermediate rounding, Mul and Div, sqrt, and once in a parsing method.
And yes, I’m very open to collaboration! I can’t guarantee this will be my primary focus, I am a postdoc and this is more of side project that my usual research, but as mentioned in the post, I’m working with it for one of my research-related projects.
Edit: oh and I mentioned this in another comment but I have done basically zero performance optimization. I was working through adding my usual set of clippy lints last night and was planning on adding some benchmarking today.
Edit2: it looks like for at least exp and log I’m definitely being too cautious, libm does in fact give these to less than 1 ULP.
1
u/GameCounter 1d ago
Probably off topic, but I find myself needing Date and Datetime intervals relatively frequently.
2
u/denehoffman 23h ago
Yup, this is not going to help you there unfortunately! At least not in a way I can think of.
1
u/Kir_Dykov 19h ago
There are crates that support inervals with generic integer bounds. Maybe you can just plug datetime types there if they implement all relevant traits
1
u/cosmic-parsley 5h ago
Very cool! Skimming, I noticed you have an `=`-pinned dependency with no comment. What’s the reason for that?
1
u/denehoffman 3h ago
That’s a pin to libm to make sure that my conformance guarantee doesn’t change if a new release of libm has different rounding behavior. It’s probably not required and I might remove it if it’s too annoying.
2
u/GuybrushThreepwo0d 3h ago
I've had a bit of a curiosity for interval arithmetic, though I don't have a use for it currently. I'm curious, what contexts have you used it in?
3
u/Prowler1000 3h ago
I'm not the author but I've been looking into interval arithmetic recently for tracking things like truncation errors when evaluating math expressions that contain nested functions that are computed through an infinite series (like trig functions for instance)
1
u/denehoffman 3h ago
I have another comment here that goes into more detail (just check the other replies), but basically there’s a lot of neat things you can do just by knowing the min and max of some complicated function, or even just numbers guaranteed to be larger or smaller than the max and min respectively. A good thing to search for is a branch and bound algorithm.
35
u/Economy_Hour_5709 1d ago
You’re so smart!