Algoritmo/formula veloce per gamma seriale di moduli di numeri co-prime

Question

Nel mio progetto, una parte del problema è lì. Ma per semplificare, qui il problema viene formulato. Esistono due numeri interi co-prime positivi: a e b, dove a < b. Multipli di a da 1 a b-1 sono elencati seguito da funzionamento modulo da b.

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

Ora, c'è un altro intero, dire n (1 <= n < b). Attraverso i primi numeri n nell'elenco, dobbiamo trovare quanti numeri è inferiore a, ad esempio m (1 <= m < b). Questo può essere fatto in approccio a forza bruta, dando così un O(n).

Un esempio:

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

List è:

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

Questa è una permutazione dei numeri da 1 a 12 per il funzionamento del modulo di qualsiasi due co-primi produce un permutazione dei numeri se includiamo un altro numero, ovvero 0. Se prendiamo lo a= 2, b=13, la lista sarebbe stata 2, 4, 6, 8, 10, 12, 1, 3, 5, 7, 9, 11, che fornisce uno schema. Considerando che se a e sono molto grandi (nel mio progetto possono andare fino a 10^20), quindi non ho idea di come dedurre un modello di numeri così grandi.

Ora tornare all'esempio, si prende la prima n = 8 numeri dalla lista, che dà

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

Applicando l'operatore less-than con m = 6, dà il numero totale dei numeri inferiori a m essere 3 come spiegato di seguito nell'elenco

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

dove 0 si riferisce alle non essere inferiore a m e 1 si riferisce ad essere inferiore a m.

Dal momento che, l'algoritmo di cui sopra è un O(n), che non è accettabile per la gamma di [0, 10^20], così può fare la comunità dare un suggerimento/indizio/suggerimento mi consentano di raggiungere una soluzione O(log n), o meglio ancora la soluzione O(1)?

Answer 1

(Attenzione: ho avuto un po 'di nervosismo sul range di moltiplicatori non essendo [0, n), quindi l'ho regolato. È abbastanza facile compensare)

Ho intenzione di disegnare, con codice Python testato, un'implementazione che viene eseguita nel tempo O (log max {a, b}). Innanzitutto, ecco alcune funzioni di utilità e un'implementazione ingenua.

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)) 

Ora, come possiamo accelerare?Il primo miglioramento consiste nel suddividere i moltiplicatori in intervalli disgiunti in modo che, all'interno di un intervallo, i multipli corrispondenti di a siano compresi tra due multipli di . Conoscendo i valori più bassi e più alti, possiamo contare tramite una divisione del soffitto il numero di multipli inferiore a 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 

Questo non è abbastanza veloce. Il secondo miglioramento consiste nel sostituire la maggior parte del ciclo while con una chiamata ricorsiva. Nel codice sottostante, j è il numero di iterazioni "complete" nel senso che è presente un avvolgimento. term3 account per l'iterazione rimanente, utilizzando la logica che è simile a count2.

Ognuna delle iterazioni complete contribuisce a floor(m/a) o floor(m/a) + 1 residui sotto la soglia m. Se otteniamo il + 1 dipende da ciò che first è per quella iterazione. first inizia da 0 e modifica da a - (b % a) modulo a su ogni iterazione attraverso il ciclo while. Otteniamo il + 1 ogni volta che è sotto una certa soglia e questo conteggio è calcolabile tramite una chiamata ricorsiva.

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 

Il tempo di esecuzione può essere analizzato in modo analogo all'algoritmo GCD euclideo.

Ecco un codice di prova per dimostrare le prove per le mie affermazioni di correttezza. Ricorda di eliminare le asserzioni prima di testare le prestazioni.

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)