Fast algorithm / formula for the serial range modulo co-primes

Question

There is one part of the problem in my project. But for simplicity, a problem is formulated here. There are two positive integers: a and b , where a < b . The set of values a from 1 to b-1 displayed, followed by the operation of module b .

a mod b , 2*a mod b , 3*a mod b , ..., (b-1)*a mod b

Now there is another integer, say n ( 1 <= n < b) . Through the first n numbers in the list, we need to find how many numbers are less, for example m ( 1 <= m < b ). This can be done using brute force, thereby giving O(n) .

Example:

a=6, b=13, n=8, m=6

List:

6, 12, 5, 11, 4, 10, 3, 9, 2, 8, 1, 7

This is a permutation of numbers from 1 to 12, because the operation of the module of any two accompanying numbers causes a permutation of numbers if we include another number, that is, 0 . If we take a= 2, b=13 , then the list would be 2, 4, 6, 8, 10, 12, 1, 3, 5, 7, 9, 11 , which gives a template. If a and b very large (in my project they can go up to 10 ^ 20), then I have no idea how to derive a pattern of such large numbers.

Now, returning to the example, we take the first n = 8 numbers from the list that gives

6, 12, 5, 11, 4, 10, 3, 9

Using the less-than operator with m = 6 , it gives a total number of numbers less than m equal to 3, as described below in the list

0, 0, 1, 0, 1, 0, 1, 0

where 0 means at least m , and 1 means that it is less than m .

Since the above algorithm is O(n) , which is unacceptable for the range [0, 10^20] , so can the community give a hint / hint / hint so that I can reach the solution O(log n ) or even better O(1) decision?

Answer 1

(Warning: I jerked a bit when the range of factors was not [0, n), so I adjusted it. It is easy to compensate.)

I am going to sketch out with tested Python code an implementation that runs on time O (log max {a, b}) . Firstly, there are some utility functions and a naive implementation.

 from fractions import gcd from random import randrange def coprime(a, b): return gcd(a, b) == 1 def floordiv(a, b): return a // b def ceildiv(a, b): return floordiv(a + b - 1, b) def count1(a, b, n, m): assert 1 <= a < b assert coprime(a, b) assert 0 <= n < b + 1 assert 0 <= m < b + 1 return sum(k * a % b < m for k in range(n)) 

Now, how can we speed it up? The first improvement is to split the multipliers into disjoint ranges so that within the range the corresponding multiples of a are between two multiples of b . Knowing the lowest and highest values, we can calculate with the help of ceiling division the number of multiples less than m .

 def count2(a, b, n, m): assert 1 <= a < b assert coprime(a, b) assert 0 <= n < b + 1 assert 0 <= m < b + 1 count = 0 first = 0 while 0 < n: count += min(ceildiv(m - first, a), n) k = ceildiv(b - first, a) n -= k first = first + k * a - b return count 

This is not fast enough. The second improvement is to replace most of the while loop with a recursive call. The code j below shows the number of iterations that are “completed” in the sense that there is a workaround. term3 takes into account the remaining iteration using logic similar to count2 .

Each of the complete iterations contributes floor(m / a) or floor(m / a) + 1 under the threshold m . Whether we get + 1 depends on what first for this iteration. first starts at 0 and changes to a - (b % a) modulo a at each iteration through a while loop. We get + 1 whenever it is under some threshold, and this score is calculated through a recursive call.

 def count3(a, b, n, m): assert 1 <= a < b assert coprime(a, b) assert 0 <= n < b + 1 assert 0 <= m < b + 1 if 1 == a: return min(n, m) j = floordiv(n * a, b) term1 = j * floordiv(m, a) term2 = count3(a - b % a, a, j, m % a) last = n * a % b first = last % a term3 = min(ceildiv(m - first, a), (last - first) // a) return term1 + term2 + term3 

Runtime can be analyzed similarly to the Euclidean GCD algorithm.

Here is some test code that confirms my statements about correctness. Remember to remove claims before performance testing.

 def test(p, f1, f2): assert 3 <= p for t in range(100): while True: b = randrange(2, p) a = randrange(1, b) if coprime(a, b): break for n in range(b + 1): for m in range(b + 1): args = (a, b, n, m) print(args) assert f1(*args) == f2(*args) if __name__ == '__main__': test(25, count1, count2) test(25, count1, count3)