Large numbers in magnitudes cause compiler crashes or long compilation times

[Tehforsch]

This is more of a curiosity than an actual issue, so feel free to close this if it doesn't seem relevant.
I was inspired by the nice talk you (@chiphogg ) gave at CppCon to implement vector space magnitudes in my rust compile time units library (not because the library is mature enough to really need them, but because rusts const-exprs like integers a lot better than floats).
While doing so I stumbled over the problem that factorizing the integers in the magnitude in the standard/naive way via trial division runs into performance problems (especially if the user really wants to run into performance problems). In my library, the magnitudes are currently read from the exponent / mantissa of a float, but the largest primes that fit into the 52+1 bit mantissa of a double (9007199254740881 for example), can already take a really long time. While implementing some more sophisticated algorithms, I checked what au does to solve this issue and realized that it also uses the standard factorization algorithm, so I tried to confirm my suspicion and compiled the following code:

struct LargePrimes : decltype(Meters{} * mag<9007199254740881>());

if compiled with clang, this gives

note: constexpr evaluation hit maximum step limit; possible infinite loop?

and with gcc, I get

au.hh:442:5: error: ‘constexpr’ loop iteration count exceeds limit of 262144 (use ‘-fconstexpr-loop-limit=’ to increase the limit)

it eventually compiles with
g++ -fconstexpr-loop-limit=2147483647 -fconstexpr-ops-limit=2147483647
but takes about ~100s on my machine. Putting in a ~64bit prime should take about an hour.

This is probably not much of a real-world concern since I maliciously wrote in large prime numbers, but if somebody, for example, were to autogenerate their unit definitions from a bunch of sensor data, they might get unlucky and accidentally stumble into this. Probably not.

There are much faster algorithms that handle 64 bit integers just fine. If you're at all interested in fixing this, I am by no means knowledgeable in number theory, but I implemented the Miller Rabin test to check for primality first and then the Pollard's Rho algorithm to factor the composite numbers and with some tweaking, that seems to work fine.

[chiphogg]

Very well spotted! This is a known issue, although not previously tracked in this project.

We encountered this in the mp-units project after I contributed vector space magnitudes there (mpusz/mp-units#300). The problem is more practical than you might imagine! Once you start using more "extreme" units, such as those found in subatomic or cosmic scales, you find magnitudes whose smallest prime factor is surprisingly large, and you hit this compiler iteration limit.

In mp-units, I assumed this problem was due to naive trial division, so I learned about wheel factorization and implemented it there. Unfortunately, it seems that for wheel division, the space requirements grow very rapidly while the iteration count reduces only slowly. Ultimately, while it might make find some answers a little faster, it's not a solution to this problem.

Our current plan is to write a paper proposing a new standard library function, perhaps something like std::first_factor. Vector space magnitudes are (as far as we know now) a hard requirement for a standard units library, and reliable prime factorization is a hard requirement for vector space magnitudes. The unix program factor, written in C, can easily factor any 64-bit integer in a tiny fraction of a second --- so we know this can be done. If we had this functionality in C++, it would make all kinds of basic number theory code easier to write.

I'm planning to finish this paper before the Feb. 15 mailing deadline, and have it discussed at the Tokyo meeting in March.

Pollard's rho algorithm sounds fascinating and exciting! I'm a little hesitant to code it up in this project just yet, because of its pseudorandom character. It sounds like it can fail, and I'm not confident I could write test cases that correctly exercise these failure cases, and cover the selection of a new parameter for the retry.

Perhaps we'll see how things develop with the paper. Perhaps we'll bite the bullet and try Pollard's rho in the meantime, ideally with plenty of tests and perhaps also fuzzing to gain enough confidence. Perhaps both!

[chiphogg]

What piqued my interest the most in your post was the idea of implementing vector space magnitudes in rust. In C++, Au and mp-units use variadic templates for the implementation. But I had thought rust didn't support variadic templates! (I've never used rust.) So how can you implement vector space magnitudes?

The reason I'm familiar with this is because I checked out the iliekturtles/uom project when it briefly ascended to the top of the dimensional-analysis star chart. I was scandalized to find that they eagerly convert all quantities of any given dimension to some base unit for that dimension! That's undesirable because it complicates end users' choice: it forces them to change the binary program they produce if they want to get unit safety.

After thinking it through, I speculated that they were more-or-less forced into this design choice because rust's lack of variadic templates prevents them from being able to represent unit magnitudes in the type. (Yes, they could have used some analogue of std::ratio, but this representation doesn't meet the basic requirements for a multidimensional units library.) Thus, my assumption was that you can't actually write a good units library in rust as it stands today. But I could easily be missing something here due to my lack of experience with the language.

[Tehforsch]

We encountered this in the mp-units project after I contributed vector space magnitudes there (mpusz/mp-units#300). The problem is more practical than you might imagine! Once you start using more "extreme" units, such as those found in subatomic or cosmic scales, you find magnitudes whose smallest prime factor is surprisingly large, and you hit this compiler iteration limit.

Our current plan is to write a paper proposing a new standard library function, perhaps something like std::first_factor. Vector space magnitudes are (as far as we know now) a hard requirement for a standard units library, and reliable prime factorization is a hard requirement for vector space magnitudes. The unix program factor, written in C, can easily factor any 64-bit integer in a tiny fraction of a second --- so we know this can be done. If we had this functionality in C++, it would make all kinds of basic number theory code easier to write.

I'm planning to finish this paper before the Feb. 15 mailing deadline, and have it discussed at the Tokyo meeting in March.

That's very interesting to hear, thanks! I was briefly looking at parsec or mole to check if their conversion factors to SI happened to have large prime factors, but then gave up on the search. Out of curiosity, is there a specific one you could point to?

Pollard's rho algorithm sounds fascinating and exciting! I'm a little hesitant to code it up in this project just yet, because of its pseudorandom character. It sounds like it can fail, and I'm not confident I could write test cases that correctly exercise these failure cases, and cover the selection of a new parameter for the retry.
Perhaps we'll see how things develop with the paper. Perhaps we'll bite the bullet and try Pollard's rho in the meantime, ideally with plenty of tests and perhaps also fuzzing to gain enough confidence. Perhaps both!

I was initially under the impression that it would always succeed if retried with (a number of) different starting values, but after some more googling I realize now that that assumption is wrong, so you're right. I think for my purposes, I can "simply" choose to use slow trial division (or throw a compiler error) if Pollard's rho doesn't find a factor but the primality test says that the number is composite. That is not exactly a great solution, but it could be an improvement, since the small primes can still be covered by trial division, so everything that worked before will still work.

Note that the normal Miller-Rabin test is also probabilistic, but apparently it's been shown to always give the correct result for 64 bit ints if a large enough set of witnesses are tested. I won't claim to have double checked that proof though.

I'd be curious to know if fuzzing can help with this problem. Very handwavy again, but I wouldn't be surprised if the unpredictable nature of integer factorization prevents fuzzing from slowly encroaching on a failing test case.

What piqued my interest the most in your post was the idea of implementing vector space magnitudes in rust. In C++, Au and mp-units use variadic templates for the implementation. But I had thought rust didn't support variadic templates! (I've never used rust.) So how can you implement vector space magnitudes?

So, I should preface this by saying that I am just in the process of implementing this and haven't seen all the parts come together, but I think that I have verified that all individual components of my implementation should work.

My library works using const generics. The currently stable version of that feature only allows one to write code that is generic over ints/bool/chars. However, there are also unstable features that allow one to use generic const expressions and to have "arbitrary" structs as const generics, which is what my library uses. These still only allow a tiny subset of what C++ templates can do, but it definitely makes things a little easier.

Now, real variadics don't exist, so my current plan is to express the vector space as a fixed size array of at most 15 integer factors (plus special entries for pi/zero/...). 15 since that should make all 64 bit integers representable (the primorial numbers are the worst case here and the number 32589158477190044730, which is the product of the first 16 primes is already too large for 64 bits).
With the unstable features, an array like that is allowed as a const generic and it shouldnt be too hard to implement multiplication etc. on such an array.

The reason I'm familiar with this is because I checked out the iliekturtles/uom project when it briefly ascended to the top of the dimensional-analysis star chart. I was scandalized to find that they eagerly convert all quantities of any given dimension to some base unit for that dimension! That's undesirable because it complicates end users' choice: it forces them to change the binary program they produce if they want to get unit safety.

After thinking it through, I speculated that they were more-or-less forced into this design choice because rust's lack of variadic templates prevents them from being able to represent unit magnitudes in the type. (Yes, they could have used some analogue of std::ratio, but this representation doesn't meet the basic requirements for a multidimensional units library.) Thus, my assumption was that you can't actually write a good units library in rust as it stands today. But I could easily be missing something here due to my lack of experience with the language.

I think your assumption is mostly correct. I don't know how exactly the decision was made in uom, but since they work on stable Rust, they rely on typenum which performs some wonderful type magic to allow them to represent integers at compile time. I suspect representing the actual unit magnitudes at compile time without const generics would probably be challenging at least.

However, I'm afraid mine and all other Rust dimensional-analysis libraries that I am aware of currently also represent quantities by converting them down to a base representation. For now, the reason I want the compile-time magnitudes is mainly to make constructing quantities easier (again, inspired by your talk - I want to be able to write v = 5.0 * meter / second or l = length.round_in(meters), and this will enable that syntax).

I don't think it is necessarily the lack of variadic templates that is keeping us from doing this properly - with const generics, we could choose to make our type generic over the exponent of every available unit (i.e. instead of { length: 1, time: -1 }, i'd have { meters: 1, seconds: -1, feet: 0, hours: 0, ... }).
This certainly gives ugly error messages, but I think the main reason most don't do this right now is that we currently cannot express even remotely as many things in Rusts trait system as one can with C++ templates. I don't know enough about C++, but it seems that you are able to state things like "perform implicit conversion between two units if the ratio of their magnitudes fulfills some criterion". I wouldn't know how to do this in Rust, which means that the only option if I want my quantities to be based on units instead of just the dimension, is to disallow all operations between quantities with different units and make every conversion explicit.
That might be an option, but for now I preferred the simplicity of the first approach. I used my library to write scientific code where the alternative seemed to be to just convert everything to a base representation anyways, but to cross your fingers and hope that nobody forgets a conversion factor - coming from there, any kind of unit safety seems preferable and having a nice user interface makes it easier to convince others to use it too.

That's undesirable because it complicates end users' choice: it forces them to change the binary program they produce if they want to get unit safety.

I'm curious what exactly you mean by this point. My impression was that the main advantage of proper unit safety was to make it possible to simultaneously represent quantities in units that are not easily convertible into one another (for example because their ratio is non-integral, or it might overflow the underlying storage type, ...), but now I am thinking I might have understood this incorrectly.

[chiphogg]

That's very interesting to hear, thanks! I was briefly looking at parsec or mole to check if their conversion factors to SI happened to have large prime factors, but then gave up on the search. Out of curiosity, is there a specific one you could point to?

Here, have two! I tracked down the motivating units.

• The proton mass definition involves 1,672,621,923,695, one of whose factors is 334,524,384,739, which is prime.
• The unit eV / c^2 involves 17,826,619,216,279, one of whose factors is 225,653,407,801, which is prime.

Now, real variadics don't exist, so my current plan is to express the vector space as a fixed size array of at most 15 integer factors (plus special entries for pi/zero/...). 15 since that should make all 64 bit integers representable (the primorial numbers are the worst case here and the number 32589158477190044730, which is the product of the first 16 primes is already too large for 64 bits). With the unstable features, an array like that is allowed as a const generic and it shouldnt be too hard to implement multiplication etc. on such an array.

Thanks for the interesting explanation! I don't yet have any experience writing rust, so I appreciate the insights.

I wonder if you might want to (conservatively) double that size of 15? After all, you'll very commonly need to represent ratios, and the numerator and denominator won't have any prime factors in common.

Of course, even this wouldn't be enough because you could have fractional exponents, and magnitudes that could be represented in double even though they couldn't be represented as a ratio of two 64-bit integers. However, I expect that in practice, 30 basis vectors would be more than enough.

I used my library to write scientific code where the alternative seemed to be to just convert everything to a base representation anyways, but to cross your fingers and hope that nobody forgets a conversion factor - coming from there, any kind of unit safety seems preferable and having a nice user interface makes it easier to convince others to use it too.

Yes --- definitely! The robust unit safety is IMO clearly worth the cost of forcing users to change their computations in order to get it. Of course, it'd be even better if you could get the benefit without the cost, but that seems pretty challenging for rust at the moment.

That's undesirable because it complicates end users' choice: it forces them to change the binary program they produce if they want to get unit safety.

I'm curious what exactly you mean by this point. My impression was that the main advantage of proper unit safety was to make it possible to simultaneously represent quantities in units that are not easily convertible into one another (for example because their ratio is non-integral, or it might overflow the underlying storage type, ...), but now I am thinking I might have understood this incorrectly.

Let me elaborate. What I had in mind here was the principle that I referred to in my talk as "same program, only safer". Imagine a user writing a program without a units library. Let's suppose that their preferred choice for some particular length variable is to use an int, and represent the result in "mils" (i.e., one thousandth of an inch).

Now let's say they want to use a units library so they can prevent unit errors.

Some libraries have one unit for each dimension. I believe uom is like this. So in switching to this units library, they would have to convert their value to meters. (Obviously, this is majorly lossy in this case. It's not clear to me whether uom would raise a compiler error, or silently truncate.)

Some other libraries have one storage type for all user facing types. The nholthaus C++ library is like this. So in this case, the user can store their value natively in mils, but they will have to use double (or float, if they select that globally) rather than int. This makes the assembly code meaningfully different, which gives the user a tradeoff that must be weighed.

What "same program, only safer" means is that the quantity type should support both the same unit and storage type that the user would have had without any units library. In this case, something like Au's Quantity<Milli<Inches>, int> fulfills these criteria. Once the compiler finishes optimizing out the units types, we'll get the same assembly we would have had with just int ("same program"), but with the added knowledge that if we made a mistake, the compiler would catch it ("only safer").

[Tehforsch]

Here, have two! I tracked down the motivating units.

The proton mass definition involves 1,672,621,923,695, one of whose factors is 334,524,384,739, which is prime.
The unit eV / c^2 involves 17,826,619,216,279, one of whose factors is 225,653,407,801, which is prime.

Thanks for the nice examples! I especially like that I do have the definition of the proton mass in the code for which I wrote my library, so it is quite likely that I would have run into this issue.

I wonder if you might want to (conservatively) double that size of 15? After all, you'll very commonly need to represent ratios, and the numerator and denominator won't have any prime factors in common.

Of course, even this wouldn't be enough because you could have fractional exponents, and magnitudes that could be represented in double even though they couldn't be represented as a ratio of two 64-bit integers. However, I expect that in practice, 30 basis vectors would be more than enough.

That's true, in order to do ratios I will probably have to go to 30 basis vectors.

What "same program, only safer" means is that the quantity type should support both the same unit and storage type that the user would have had without any units library. In this case, something like Au's Quantity<Milli, int> fulfills these criteria. Once the compiler finishes optimizing out the units types, we'll get the same assembly we would have had with just int ("same program"), but with the added knowledge that if we made a mistake, the compiler would catch it ("only safer").

Thank you, I understand your point now. I believe implementing the "same program, only safer"-approach in Rust could be possible, but only at the cost of lots of convenience (and error message quality), which could make it a bit harder to convince people to move to a unit library. This is definitely something to think about though.

[chiphogg]

@Tehforsch, try the newly landed changes on main!

[Tehforsch]

@Tehforsch, try the newly landed changes on main!

Cool, thanks!

[chiphogg]

get_value<uint64_t>(mag<585226005592931977>())

This fails to compile, because it's composite, but the prime factors are too large. (I got the idea by browsing https://miller-rabin.appspot.com/.)

It seems like the next move would be to try Pollard's Rho algorithm for factoring known composites with no small prime factors, using Brent's cycle detection method.

[chiphogg]

I got Pollard's rho working with Brent's cycle detection method. I even used a variant of the batching improvement, which reduced the number of steps.

The bad news is that it's still not enough to factor 585226005592931977 at compile time. The good news is that I was able to trace the problem to our modular multiplication function, which appears to be way more expensive than I thought. (It's an algorithm of my own original design --- go figure.)

The key to solving this was a hacky little script that measures the minimum value of -fconstexpr-steps we need to make clang happy:

#!/bin/bash

# This script finds the minimum necessary value of `-fconstexpr-steps`.

# The command to run.
function command() {
  printf "\n\n\n\nTesting with -fconstexpr-steps=%d\n\n\n\n" $1
  bazel test --copt=-fconstexpr-steps=$1 //au:utility_test
}

# Increase exponentially until we find a value that passes.
I="1048576"
INC="$I"
while ! command $I; do
  I=$((I + INC))
  INC=$((INC * 2))
done
INC=$((INC / 2))

while [ $INC -gt 1 ]; do
  INC=$((INC / 2))
  if command $((I - INC)); then
    I=$((I - INC))
  fi
done

printf "The minimum value of -fconstexpr-steps is %d\n" $I

I did three arrangements. First, the current modular multiplication implementation. Then, a version that calls the code twice every time an external user calls the function --- this is a way to estimate the cost of calling it once (by noting the increase), and put a ceiling on how much speedup we could hope to get in the very best case. Finally, I replaced our algorithm with a simple (but non-portable) 128-bit multiplication, as a way to get that best-case scenario.

Here are the results, for both the "regular" Pollard's rho (using Brent's cycle detection, of course), and the batched version (which is mainly notable for using more multiplications, in order to buy fewer expensive GCD calls):

Multiplication	Steps (regular rho)	Steps (batched rho)
Current	4'802'069	3'959'173
Double cost	7'121'157	7'790'206
128-bit	2'567'401	268'250

Conclusions:

1. The batched version is a meaningful win even with our current modular multiplication implementation.

2. Doubling the cost of multiplication has staggering effects. This is even more pronounced for the batched version, where it is literally almost double. This implies that making multiplication more efficient should remove almost all of the cost!

3. Using 128-bit multiplication indeed removes almost all of the cost, and puts us well under whatever clang's default is (I think it's roughly 1,000,000).

Thus, I foresee roughly four PRs:

1. Change the contract from find_first_factor to find_prime_factor, paving the way for something like Pollard's rho where we don't guarantee we'll get the "first" one.

2. Use "regular" Pollard's rho.

3. Switch to the "batched" version.

4. Find a faster implementation of modular multiplication.

For the latter, I think Montgomery multiplication looks very promising. I got the idea from the README for hurchalla/modular_arithmetic.

[chiphogg]

I think this is done enough for now. Adding Pollard's rho instead of more trial division is already a huge improvement. I filed #328 to track a potential future round of improvements.