mpusz / mp-units

The quantities and units library for C++
https://mpusz.github.io/mp-units/
MIT License
1.08k stars 86 forks source link

Faster than lightspeed Constants #169

Closed TobiasWallner closed 1 year ago

TobiasWallner commented 4 years ago

Would it be possible to make a type for every constant so that the value stored is always 1. In that way addition, subtractions multiplications, and divisions will always be the fastest - compiled away or done in out-of-order execution.

so for example the electric vacuum permitivity could be implemented with the template arguments auto eps_0 = unit<... ratio(8.854, 1000), exp = -12> {1}; and the permeability as auto mu_0 = unit< ... ratio(8.854, 1000), exp = -12> {1};

If you would do that you could calculate for example the light speed as: auto c_0 = 1 / sqrt(eps_0 mu_0) // the compiler and CPU at runtime would see 1 / sqrt(1 1) = 1;

then you could go even further since everything is 1 at runtime you could use integral types instead of floating-point types giving you an extra speed advantage. The only issue would be that you lose information if were to divide by a constant that is larger than one.

for example: eps_0 / 4; // ... would return zero.

We could prevent that by using a templated constant wrapper that stores the constant inside the template ratio template auto tw_const = unit<... ratio(a,1)>{1};

then we can write eps_0 / tw_const<4> // will be 1 / 1 = 1 at runtime

I know at some point in the code, for example when printing the value we need to unwrap the whole information and multiply with all the accumulated constants. But the key is that we only perform one multiplication which will take care of all constants from all functions this value went trough.

mpusz commented 4 years ago

Hi @TobiasWallner. I generally like your idea. The biggest problem with this approach is a finite precision of a compile-time ratio arithmetics. We already had issues reported for the cases where a cubic version of a unit overflown ratio or similar cases (e.g. #55, #60). Even though we have the third Exp parameter in ratio it helps only with SI prefixes which are the magnitude of 10. The more strange ratios we will throw into the compiler the faster we will end up with problems.

If anyone has an idea of how to solve the limited precision of ratio please let me know.

TobiasWallner commented 4 years ago

You do not need to limit yourself. you just need to calculate new template parameters without overflow. You can use constexpr functions in combination with auto to calculate the return type of a function, within that function, that is returning that type. .... this definitely sounds more complicated than it is.

Example in Compiler Explorer: https://godbolt.org/z/64Wzxf

Is your problem solved?

TobiasWallner commented 4 years ago

^^note that in my example I divide by 10 but decrement the bit count by 3 which is equivalent to 8, so i am kinda lazy here but since 10 > 8 i am save.

TobiasWallner commented 4 years ago

I just noticed that my algorithm for handling ratio overflows only works for positive numbers. Dou you guys need negative ratios too?

mpusz commented 4 years ago

Using floating-point for ratio is not an option. We cannot afford even one ULP difference for the same value/ratio because it will not result in the same type/unit.

TobiasWallner commented 4 years ago

I am not using floating point types!

TobiasWallner commented 4 years ago

Oh! I see the confusion. I named the struct floatlong. This is because of the template argument 'exp' the long is floating. The value


template <typename ratio, long exp> //parameters are long
struct floatlong{
    long value;   //<--- is a long type 

    double getFactor() { //<-- just used in the example to print to the console
        static constexpr double factor = (ratio::num / (double)ratio::den) * std::pow(10,exp);
        return factor;
    }
};

struct FactorPower {
        long long num; //<-- long long
        long long den; // <-- long long
        long exp; // <-- long 
};

I am using std::ratio and the struct Factor power can be compared to what you call a ratio in your project.

mpusz commented 4 years ago

I was refering to getFactor() returning double.

TobiasWallner commented 4 years ago

yes you don't need that. I already did a push request adding my implementation of multiplications without overflow to your ratio. you will see that no floating point types were used

mpusz commented 4 years ago

Great, thank you!

Could you also tell if this solves the following issue:

#include <units/physical/si/base/length.h>

using namespace units::physical::si;

constexpr auto cubic_au = 1_q_au * 1_q_au * 1_q_au;

Even if it does, imagine multiplying it by constants with huge ratios. How to solve that?

TobiasWallner commented 4 years ago

It is fixing it in the way that no overflow occurs but the result differs in the last digits. So you won't get the precise result.

I once did a large number library maybe we could use it instead of long long to arbitrary precision results? Unfortunately, I made it some years ago in C. I would need to port it to modern C++.

TobiasWallner commented 4 years ago

The relative error is approx 9E-9. ... I have to lose information at some point when dealing with those numbers Seems small but when dealing with E17 or E50 ratios?

I just tested: auto test = cubic_au * cubic_au and it still holds

TobiasWallner commented 4 years ago

oh just noticed there was a logic errror. needed an xnor instead of a xor operation

TobiasWallner commented 4 years ago

the only error that can occur now is if you were to multipy two ratios where the exponent would overflow ... i cannot do anything against that. Except an static_assert maybe? Or an arbitrary precision integer.

TobiasWallner commented 4 years ago

I managed to make it a little bit more precise. now the relative precision is around 8E-10 when doing constexpr auto cubic_au = 1_q_au * 1_q_au * 1_q_au;

mpusz commented 4 years ago

The overflow should be detected at compile-time like it was until now and like it is in case of std::ratio.

mpusz commented 4 years ago

The relative error is approx 9E-9. ... I have to lose information at some point when dealing with those numbers Seems small but when dealing with E17 or E50 ratios?

The human will not care about such an error but the compiler will yield a different unit type for it and this is a no go :-(

No matter what calculations (how many multiplications and divisions) are involved in the way we need to always end up with exactly the same unit type every single time.

TobiasWallner commented 4 years ago

In that case you will need a larger integer type. In the example for au I am already pulling off all tricks and forcing denominators to try reducing the number 149597870700 to 1495978707 which can be squared without information loss. The result 22379522917973918490000 can be reduced to 2237952291797391849 but 9 is not divisible by 5 or 2 so i cannot steal a 10 from the exponent and the third au ratio also has no denominator so I cannot reduce using the other ratio. At this point I have to round up to 2237952291797391850 which can then be reduced to 44759045835947837 / den. And here again I need to lose infromation.

So you will need a larger datatype to preserve all the information.

mpusz commented 4 years ago

Yes, I know.

As I said at the beginning of this thread we already have such problems when ratio is only applied to units. When we apply ratio also to constants we will have even more issues. This is why I think we should not do it as long as we will not solve the ratio resolution problem.

One solution would be to make ratio a class template taking the representation type as a parameter. However, this might produce interop issues between different types?

Another solution would be to use some flexible integer implementation but this is not standardized yet and it will be hard to standardize units with such an internal type (if such type will not reach the standard first).

TobiasWallner commented 3 years ago

Yes, I know.

As I said at the beginning of this thread we already have such problems when ratio is only applied to units. When we apply ratio also to constants we will have even more issues. This is why I think we should not do it as long as we will not solve the ratio resolution problem.

One solution would be to make ratio a class template taking the representation type as a parameter. However, this might produce interop issues between different types?

Another solution would be to use some flexible integer implementation but this is not standardized yet and it will be hard to standardize units with such an internal type (if such type will not reach the standard first).

Hey, its been a long time.

I kind of solved the integer overflow problem. I made a fully constexpr variable length integer called aint, where n is the number of bits representing the integer. This would allow to calculate the size needed to represent new integer values. So if you still wanna go for integer representations of ratios and won't use some kind of lamdas you may do it like so:

template<unsigned long num_bits, unsigned long den_bits, aint<num_bits> _N, aint<den_bits> _D>
struct const_aratio{
    static constexpr auto N = _N;
    static constexpr auto D = _D;

template<unsigned long rn_bits, unsigned long rd_bits, aint<rn_bits> Num, aint<rd_bits> Den>
auto operator*(const const_aratio<rn_bits, rd_bits, Num, Den>& rhs)
{
    static constexpr auto NN = mul_s(this->N, rhs.N);
    static constexpr auto DD = mul_s(this->D, rhs.D);
    static constexpr const_aratio<NN.Bits, DD.Bits, NN, DD> RR;
    return RR;
}
};

int main() {
    bool test = true;

    //test const_aratio

    static constexpr aint<64> n = 446744073709551616;
    static constexpr aint<64> d = 856879542347894362;
    static constexpr const_aratio<n.Bits,d.Bits,n,d> r1;
    static constexpr auto r2 = r1 * r1;
}

The only thing I am struggling with is, that "static constexpr auto r2 = r1 * r1" this deduction fails sometimes.

mpusz commented 3 years ago

Thanks. Unfortunately, I am still extremely busy in the next 10 days. I should have more time to look into this after Dec 3.

I never said that I do not like a lambda approach but I do not see how we can easily replace current ratio logic with lambdas. If someone could create a short POC it would be great. I will also try to spend more time on it in December.

mpusz commented 3 years ago

In fact, I believe angles and ratio issues are the last ones to solve before we could start an ISO standardization effort.

mpusz commented 1 year ago

Done in V2.