| | |
| Small primes: |
| Batch trial division |
| rho |
| p-1 |
| ECM |
| |
|
Batch trial division
Given
a sequence P of primes
and a sequence X of nonzero integers,
the batch-trial-division algorithm
figures out
- which of the primes in P divide X[0],
- which of the primes in P divide X[1],
- which of the primes in P divide X[2],
- etc.
Standard trial division does the same job
by handling each X[i] separately,
dividing each X[i] by each P[j].
Batch trial division is much faster.
How it works
Batch trial division has four steps:
- Use a product tree
to efficiently compute the product
x = X[0]X[1]X[2]...
- Use a remainder tree
to efficiently figure out which of the primes divide x.
Remove the other primes from P:
those primes can't divide any of the original integers.
- If X has only a single element, print the new P and stop.
- Recursively apply batch trial division
to the first and second halves of X, in both cases using the new P.
The Sage code has three layers.
The first layer, primesin,
reads a sequence P of primes and an integer x,
and returns the list of primes in P that divide x:
def primesin(P,x):
result = remainders(x,P)
return [p for p,r in zip(P,result) if r == 0]
The second layer, primesinproduct,
reads a sequence P of primes and a sequence X,
and returns the list of primes in P that divide the product of the elements of X:
def primesinproduct(P,X):
return primesin(P,product(X))
The third layer, primesineach,
is the complete batch-trial-division algorithm,
returning for each X[i]
the list of primes in P that divide X[i]:
def primesineach(P,X):
n = len(X)
if n == 0: return []
P = primesinproduct(P,X)
if n == 1: return [P]
return primesineach(P,X[:n//2]) + primesineach(P,X[n//2:])
# example:
P = list(primes(2,10))
X = [50,157,266,377,490,605]
print primesineach(P,X)
# output: [[2, 5], [], [2, 7], [], [2, 5, 7], [5]]
# double check:
print [[p for p in P if x%p==0] for x in X]
# output: [[2, 5], [], [2, 7], [], [2, 5, 7], [5]]
Version:
This is version 2012.12.27 of the batchtrial.html web page.
|