Implement `powi` without using MPFR

[Chris00]

I've implemented powi without the use of MPFR, using repeated squaring (distinguishing cases according to the sign of the interval bounds). It is 4.5–5 times faster than the MPFR version — something that may matter when powi is a building block of code that will be run millions of times. However, the returned intervals are a few ulp larger so the tests no longer pass. One example is tests/itf1788_tests/libieeep1788_elem.rs:3567:5:

    assert_eq2!(
        n2i(13.1, 13.1).powi(8),
        n2i(867302034.6900622, 867302034.6900623)
    );

where I get (with the new Debug printing):

thread 'itf1788_tests::libieeep1788_elem::minimal_pown_test' panicked at tests/itf1788_tests/libieeep1788_elem.rs:3567:5:
assertion failed: `(left == right)`
  left: `Interval [867302034.6900619, 867302034.6900629]`,
 right: `Interval [867302034.6900622, 867302034.6900623]`

How important is it for the intervals to be tight? If it is, one may locally use double-doudle arithmetic to get the rounding correct. There is the qd crate that can serve as an inspiration. I think however it is better to implement the operations (we need very few of them) ourselves to

• get upward rounding;
• use __m256d to process two double-double simultaneously.

Before doing the work, I'd like to have your opinion about this endeavor.

[Chris00]

P.S. One may also have a feature flag fast-powi to choose the slightly less accurate but faster version.

[unageek]

That sounds promising! How about introducing a module that serves non-standard or non-tightest operations, including the fast powi?

As for double-double arithmetic, https://github.com/mskashi/kv could be worth looking into, although it is written in C++.

[Chris00]

How about introducing a module that serves non-standard or non-tightest operations, including the fast powi?

Unfortunately, with such a scheme, one needs to rewrite the code to change the implementation (one may imagine using the precise functions first and wanting to switch to faster but less precise ones if the code is too slow for the problem.

Another possibility would be to implement these functions as various traits, the one in scope saying which version is used. This has the drawback that a prelude is desirable: if f is implemented as a method then it will be used without warning instead of the one of the trait — even if the trait is in scope. Thus all implementations of, say, transcendental functions must be in traits. If there is no prelude and one needs transcendental functions, the trait must be present in use inari::{...}; (which may be Ok too).

Alternatively, none of the transcendental functions are methods and one needs to use inari::...::* (or use inari::... as ... if one does not want to bring these functions in scope) the appropriate module to access them.

[unageek]

@Chris00 Could you share some code snippets that illustrate your design ideas?

[Chris00]

Here is an example. In the inari library, on would define

pub trait TranscendentalMPFR {
   fn cos(self) -> Self;
}

pub trait TranscendentalFast {
   fn cos(self) -> Self;
}

and implement these traits for Interval.

In the "client" code, one can write

use inari::TranscendentalMPFR;

or

use inari::TranscendentalFast;

to have access to the appropriate cos method. In any case, the fully qualified name such as in inari::TranscendentalFast::cos(x) is always accessible. One can also use modules to "hide" the trait imported by super and bring another trait into scope. Example:

use inari::{Interval as I, TranscendentalMPFR};

mod Fast {
    use super::I;
    use inari::TranscendentalFast;

    fn f(x: I) -> I {
        x.cos()
    }
}

[unageek]

Thank you for the clarification. What if you have two or more implementations of the cosine function that vary in terms of speed and tightness trade-offs? How do we arrange such functions in traits?

[Chris00]

Each implementation must be in its own trait. To minimize the number of traits, these should group functions with similar characteristics. This is not much different than what you would do with modules — it's the calling (method v.s. function) that differs.

[Chris00]

Another possibility may be to use a phantom type to mark the kind of operations supported by the intervals. Here is some proof of concept code:

struct Itv<F> {
    // ...
    marker: PhantomData<F>,
}

enum A {}
enum B {}

impl Itv<A> {
    fn cos(self) -> Self { todo!() }
}

impl Itv<B> {
    fn cos(self) -> Self { todo!() }
}

In order to switch to a different implementation a zero cost conversion function Itv<A> → Itv<B> is necessary.