Use modulo with numbers greater than 64bit integer in numpy/numba

ueznem 0
Viewed 69

I‘m trying to implement a prime factorization algorithm leveraging the GPU/CUDA for parallelization as a pet project.

I‘m using python with numpy and numba for the parallelization part.

My problem is now that I hit the 64 bit integer boundary quite fast and I am searching for solutions to work around this. Since numpy/numba only supports integers up to 64 bit (compared to arbitrary large numbers in python itself), this is where I‘m stuck at the moment.

In essence the part on the GPU mainly uses the modulo operation and a bit of iteration.

I found that I can use bigger dividends on the modulo by splitting the modulo into multiple operations with a smaller dividend.

Example:

Let’s assume we cannot use numbers bigger than 10^3.

I would like do do the following operation on the GPU:

1019 % 17 = 16

I can do so by splitting the dividend 1019 into multiple arbitary sizes chunks, for example:

1019 = 500 + 519

And then calculating the modulo on all chunks separately and taking the modulo on the sum of the results again.

((500 % 17) + (519 % 17)) % 17 = 16

My question is now: Is there a similar operation I can perform to work with a bigger divisor and quotient (for example 2037 % 1019)? Or even better, a way to use arbitrary sized numbers in numpy/numba without precision loss?

Bear with me if I didn’t use proper math slang.

math numpy numba cuda python
1 Answers

You can reduce the dividend much more efficiently using modulus thanks to basic congruent identity rules. Indeed:

Assuming:
  a ≡ x [d]
  b ≡ y [d]
  c ≡ z [d]

Then:
  a * b + c ≡ x * y + z [d]

This is an interesting property if d is small. Additionally, all number can be decomposed in a sum of power of two values and it is easy to pre-compute the remainders of 2**n modulus d using modular exponentiation.

For example:

15428794586587458 ≡ 3592296 * 2**32 + 749035842   [2019]
                  ≡ 495 * 1090 + 975              [2019]
                  ≡    477     + 975              [2019]
                  ≡         1452                  [2019]

Unfortunately, when it comes to the divider, AFAIK there is no simple efficient (ie. polynomial) method to reduce with special assumption on both the divider and the dividend (especially the former). I think prime factorization would be much simpler to solve with such a reduction. Assuming the divisor d is a composite number, then the Chinese remainder theorem should certainly help to split the division in simpler pieces but if d is a prime number then AFAIK this is a hard problem.

In that case, using variable precision numbers (or simply large integers in multiple parts) is certainly the simplest solution. A naive solution is to use a binary search since a ≡ b [d] is equivalent to solve a = k * d + b. Indeed, k can be multiplied by two over and over until k * d is bigger than a and then a binary search on k can be used so a ≤ k * d < a + d. The multiplication can be efficiently computed using the Karatsuba algorithm (there are faster method for very large numbers but this one is faster for reasonably large ones certainly used in your case). AFAIK, a Newton's method can be used to archive a faster reduction (both are algorithmically fast but the later converge faster).

Jérôme Richard0
Related