On the Factor Refinement Principle and its Implementation on - - PowerPoint PPT Presentation
On the Factor Refinement Principle and its Implementation on Multicore Architectures Masters Public Lecture Presented By: Md. Mohsin Ali Supervisor: Dr. Marc Moreno Maza Dept. of Computer Science (ORCCA Lab) The University of Western Ontario,
The University of Western Ontario, London, ON, Canada
1 / 65
2 / 65
3 / 65
◮ Let D be a UFD and m1, m2, . . . , mr be elements of D. ◮ Let m be the product of m1, m2, . . . , mr. ◮ We say that elements n1, n2, . . . , ns of D form a GCD-free basis
◮ Let e1, e2, . . . , es be positive integers. ◮ We say that the pairs (n1, e1), (n2, e2), . . . , (ns, es) form a
1≤j≤s nfj j = mi,
1≤i≤s nei i = m.
4 / 65
5 / 65
◮ Simplifying systems of polynomial equations and inequations,
S1 S2 S3 S1 S2 S3 A B C D E F G = ⇒ ,
◮ consolidation of independent factorizations, ◮ etc.
6 / 65
7 / 65
◮ Given a partial factorization of an integer m, say m = m1m2,
◮ This process is continued until all the factors are pairwise
◮ This is also used for the general case of more than two inputs,
8 / 65
9 / 65
◮ Basic idea same as before but organizing the computations
◮ The trick is to keep tracks of the pairs (nj, nk) in an ordered
10 / 65
11 / 65
◮ Divide the input into sub-problems until a base case is
◮ Conquer the sub-problems from leaves to the root applying
2 d), where ◮ d is the sum of the degrees of the input polynomials, ◮ M(d) is a multiplication time of two univariate polynomials of
12 / 65
◮ by Leiserson, Li, Moreno Maza, and Xie in 2010 [ELLMX10]. ◮ This is a multithreaded (fork-join parallelism) approach which
◮ Its Cilk++ shows nearly linear speed-up on 32-core processors
13 / 65
14 / 65
◮ A multi-core processor is a single computing component with
◮ They also share the same bus and memory controller; thus
◮ In order to maintain memory consistency, synchorization is
15 / 65
◮ This model represents the execution of a multithreaded
◮ Assuming unit cost of execution for all threads, the number of
◮ The maximum length of a path from the root to a leaf is the
◮ The paralleisim is the ratio work to span (= average amount
16 / 65
17 / 65
◮ The processor can only refer to words that reside in the cache
◮ If the referenced line of a word is not in cache, the
◮ Cache complexity analyzes algorithms in terms of cache
18 / 65
◮ Cilk (resp. Cilk++) is an extension of C (resp. C++)
◮ A Cilk (resp. Cilk++) program has the same semantics as
◮ for almost all parallel steps, there are at least p units of work
◮ each processor is either working or stealing work, ◮ each thread executes in unit time.
19 / 65
◮ In practice, the observed speedup factor may be less
◮ Many factors explain that fact: simplification assumptions of
◮ Cilkview is a perforance analyzer which caclulates the work,
◮ Cilkview also estimates the running time Tp on p processors
20 / 65
◮ Parallel algorithm based on the naive refinement principle
◮ Parallel algorithm based on the augment refinement principle
◮ Parallel algorithm based on subproduct tree [MORE
21 / 65
22 / 65
2, 6, 7, 10, 15, 21, 22, 26 2, 6, 7, 10 15, 21, 22, 26 2, 6 7, 10 15, 21 22, 26 2 6 7 10 15 21 22 26 21 61 71 101 151 211 221 261 31, 22 71, 101 51, 71, 32 111, 131, 22 31, 71, 51, 23 51, 71, 32, 111, 131, 22 111, 131, 33, 72, 52, 25 Input Output Expand Merge Done in parallel Done in parallel
23 / 65
24 / 65
1:
2:
3:
4:
5:
6:
7:
8:
25 / 65
Input: Arrays of square-free positive integers A = l1, . . . , lk, E = e1, . . . , ek, B = m1, . . . , mr and F = f1, . . . , fr, where A, E (resp. B, F) is regarded as a factor refinement (l1, e1), . . . , (lk, ek) (resp. (m1, f1), . . . , (mr, fr)). Output: A factor refinement of the concatenation of A, E and B, F. Allocate G and H of (k × r) row-major arrays, A′ an k-array and B′ an r-array; 1: parallel for (i, j) ∈ {1, . . . , k} × {1, . . . , r} do 2: Hi,j ← 0; G ← GcdOfAllPairs(A, B) ; 3: parallel for i = 1 to k do 4: pi ←
1≤j≤r Gi,j ;
5: A′
i ← li/pi ;
6: for j = 1 to r do 7: if Gi,j = 1 then 8: Hi,j ← Hi,j + ei; 9: For each column, do similar work, from Line 4 to 9 ; 10: Write the non-one entries of A′ (for each row), B′ (for each column), and G to C; 11: Write the corresponding exponents from E, F, and H to D; 12: return C, D ; 13:
26 / 65
Input: 1-D arrays of positive integers A = l1, l2, . . . , lk, B = m1, m2, . . . , mr and a 2-D array G = [Gi,j | 1 ≤ i ≤ k, 1 ≤ j ≤ r]. Output: All possible pairs of gcd(li, mj) for 1 ≤ i ≤ k, 1 ≤ j ≤ r with gcd(li, mj) stored in Gi,j. comment C is a global variable equal to a base-case threshold, say 16. Assume C ≥ 2. if k ≤ C and r ≤ C then
1:
for (i, j) ∈ {1, . . . , k} × {1, . . . , r} do
2:
Gi,j ← gcd(Ai, Bj); [Place of INEFFICIENCY for cache complexity]
3:
else if k > C and r ≤ C then
4:
Divide A into two halves A1 = l1, . . . , lk/2, A2 = lk/2+1, . . . , lk ;
5:
GcdOfAllPairsInner(A1, B, G);
6:
GcdOfAllPairsInner(A2, B, G);
7:
else if k ≤ C and r > C then
8:
Divide B into two halves B1 = m1, . . . , mr/2, B2 = mr/2+1, . . . , mr ;
9:
GcdOfAllPairsInner(A, B1, G);
10:
GcdOfAllPairsInner(A, B2, G);
11:
else
12:
Divide A and B arrays into two halves
13:
A1 = l1, . . . , lk/2, A2 = lk/2+1, . . . , lk; B1 = m1, . . . , mr/2, B2 = mr/2+1, . . . , mr ; spawn GcdOfAllPairsInner(A1, B1, G);
14:
spawn GcdOfAllPairsInner(A2, B2, G);
15:
spawn GcdOfAllPairsInner(A1, B2, G);
16:
spawn GcdOfAllPairsInner(A2, B1, G);
17:
27 / 65
◮ Work = O(n2), ◮ Span = O(n), and ◮ Parallelism = O(n). GOOD! :-)
28 / 65
◮ there exists a positive integer w such that each integer
◮ each input array A or B is packed, that is, its successive
29 / 65
◮ The speedup estimates are obtained by Cilkview feature of
◮ Timings are obtained by running in an Intel(R) Core(TM) i7
◮ Timing results are obtained with the average value of 5
30 / 65
2 4 6 8 10 12 14 16 2 4 6 8 10 12 14 16 Speedup Cores Parallelism = 2395.73432, Ideal Speedup Lower Performance Bound Measured Speedup
31 / 65
10 20 30 40 50 60 100 200 300 400 500 600 Time (Sec) Input size Algorithm: Cores = 1 Algorithm: Cores = 8 Maple
32 / 65
33 / 65
1:
2:
3:
4:
5:
6:
7:
8:
34 / 65
Input: Two sequences of square-free pairwise coprime polynomials A = (a1, . . . , an) and B = (b1, b2, . . . , br) ∈ K[x] together with two sequences of positive integers E = (e1, e2, . . . , en) and F = (f1, f2, . . . , fr). Output: A factor refinement of (a1, e1), . . . , (an, en), (b1, f1), . . . , (br, fr). if n ≤ BASESIZE or r ≤ BASESIZE then 1: return MergeRefineTwoSeq(A, E, B, F); 2: else 3: Divide A, E, B, and F into two halves called A1, A2, E1, E2, B1, B2, and F1, F2, 4: respectively ; (L1, M1, Q1, R1, S1, T1) ← spawn MergeRefinementDNC(A1, E1, B1, F1) ; 5: (L2, M2, Q2, R2, S2, T2) ← spawn MergeRefinementDNC(A2, E2, B2, F2) ; 6: sync ; 7: (L3, M3, Q3, R3, S3, T3) ← spawn MergeRefinementDNC(L1, M1, S2, T2) ; 8: (L4, M4, Q4, R4, S4, T4) ← spawn MergeRefinementDNC(L2, M2, S1, T1) ; 9: sync ; 10: return 11: {L3 + L4, M3 + M4, Q1 + Q2 + Q3 + Q4, R1 + R2 + R3 + R4, S3 + S4, T3 + T4} ;
35 / 65
Input: Two sequences of square-free pairwise coprime polynomials A = (a1, . . . , an) and B = (b1, b2, . . . , br) of K[x] together with two sequences of positive integers E = (e1, e2, . . . , en) and F = (f1, f2, . . . , fr). Output: A factor refinement of (a1, e1), . . . , (an, en), (b1, f1), . . . , (br, fr). L ← ∅; 1: M ← ∅; 2: Q ← ∅; 3: R ← ∅; 4: S0 ← B; 5: T0 ← F; 6: for i from 1 to n do 7: (ℓi, mi, Qi, Ri, Si, Ti) ← MergeRefinePolySeq(ai, ei, Si−1, Ti−1) ; 8: Q ← Q + Qi ; 9: R ← R + Ri ; 10: if ℓi = 1 then 11: L ← L + (ℓi) ; 12: M ← M + (mi) ; 13: return (L, M, Q, R, Sn, Tn); 14:
36 / 65
Input: A square-free polynomial a ∈ K[x], a positive integer e, a sequence of square-free pairwise coprime polynomials B = (b1, b2, . . . , bn) of K[x] and a sequence of positive integers F = (f1, f2, . . . , fn). Output: A factor refinement of ae, b1f1, . . . , bnfn. ℓ0 ← a; 1: m0 ← e; 2: Q ← ∅; 3: R ← ∅; 4: S ← ∅; 5: T ← ∅; 6: for i from 1 to n do 7: (ℓi, mi, Gi, Vi, di, wi) ← PolyRefine(ℓi−1, mi−1, bi, fi) ; 8: Q ← Q + Gi ; 9: R ← R + Vi ; 10: if di = 1 then 11: S ← S + (di) ; 12: T ← T + (wi) ; 13: return (ℓn, mn, Q, R, S, T); 14:
37 / 65
Input: Two square-free univariate polynomials a, b ∈ K[x] for a field K and two positive integers e, f . Output: A factor refinement of (a, e), (b, f ) . g ← gcd(a, b); 1: a′ ← a quotient g; 2: b′ ← b quotient g; 3: if g = 1 then 4: return (a, e, ∅, ∅, b, f ) ; // Here ∅ designates the empty sequence 5: else if a = b then 6: return (1, 1, (a), (e + f ), 1, 1) 7: else 8: (ℓ1, e1, G1, V1, r1, f1) ← PolyRefine(a′, e, g, e + f ); 9: (ℓ2, e2, G2, V2, r2, f2) ← PolyRefine(r1, f1, b′, f ); 10: if ℓ2 = 1 then 11: G2 ← G2 + (ℓ2) ; // Here + designates sequence concatenation 12: V2 ← V2 + (e2); 13: return (ℓ1, e1, G1 + G2, V1 + V2, r2, f2) ; 14:
38 / 65
◮ Work = O(n2), ◮ Span = O(Cn), and ◮ Parallelism = O(n/C). NOT BAD! :-)
39 / 65
◮ each input or output sequence in Algorithms 8, 7, 6 or 5 is
◮ each element of the field K can be stored in one machine
40 / 65
2 4 6 8 10 12 14 16 2 4 6 8 10 12 14 16 Speedup Cores Parallelism = 170.733083, Ideal Speedup Lower Performance Bound Measured Speedup
41 / 65
5 10 15 20 25 30 35 40 45 50 5000 5200 5400 5600 5800 6000 Time (Sec) Input size Algorithm: Cores = 1 Algorithm: Cores = 8 Maple
42 / 65
2 4 6 8 10 12 14 16 2 4 6 8 10 12 14 16 Speedup Cores Parallelism = 265.204981, Ideal Speedup Lower Performance Bound Measured Speedup
43 / 65
20 40 60 80 100 120 140 100 200 300 400 500 600 Time (Sec) Input size Algorithm: Cores = 1 Algorithm: Cores = 8 Maple
44 / 65
2 4 6 8 10 12 14 16 2 4 6 8 10 12 14 16 Speedup Cores Parallelism = 136.931277, Ideal Speedup Lower Performance Bound Measured Speedup
45 / 65
5 10 15 20 25 30 35 40 45 50 200 300 400 500 600 Time (Sec) Input size Algorithm: Cores = 1 Algorithm: Cores = 8 Maple
46 / 65
47 / 65
◮ Useful construction to devise fast algorithms with univariate polynomials, ◮ If m0, m1, . . . , mn−1 be monic, non-constant, polynomials in K[x], then its subproduct tree is as follows [cost: O(M(d) log2(d)), d = n−1
i=0 degree(mi)] m0.m1 . . . mn−1 m0.m1 . . . mn/2−1 mn/2.mn/2+1 . . . mn−1 m0.m1 m2.m3 mn−2.mn−1 m0 m1 m2 m3 mn−2 mn−1 . . . . . . . . . . . . . . . . . . i = k i = k − 1 i = 1 i = 0
48 / 65
◮ Let m0, m2, . . . , mn−1 be as before and assume their subproduct tree M has been computed. ◮ To compute f rem m0, f rem m1, . . . , f rem mn−1, one reduces f by the root M, then by its children (leading to two remainders) and so on. ◮ deg(f ) < n, this amounts to O(M(d) log2(d)) field operations.
m0.m1 . . . mn−1 m0.m1 . . . mn/2−1 mn/2.mn/2+1 . . . mn−1 m0.m1 m2.m3 mn−2.mn−1 m0 m1 m2 m3 mn−2 mn−1 . . . . . . . . . . . . . . . . . . i = k i = k − 1 i = 1 i = 0
49 / 65
◮ Suppose we have a polynomial f ∈ K[x] and a sequence of square-free polynomials A = m0, m1, . . . , mn−1 in K[x] with their subproduct tree. ◮ We compute multiGcd(f , A), that is, all gcd(f , ai) for 0 ≤ i ≤ (n − 1) as follows:
a0.a1 . . . an−1 a0.a1 . . . an/2−1 an/2.an/2+1 . . . an−1 a0.a1 a2.a3 an−2.an−1 a0 a1 a2 a3 an−2 an−1 . . . . . . . . . . . . . . . mk,0 = f mod (a0.a1 . . . an−1) mk−1,0 = mk,0 mod (a0.a1 . . . an/2−1) mk−1,1 = mk,0 mod (an/2.an/2+1 . . . an−1) m1,0 = m2,0 mod (a0.a1) m0,0 = m1,0 mod a0 return gcd(m0,0, a0) m0,1 = m1,0 mod a1 return gcd(m0,1, a1)
50 / 65
◮ Suppose we have two sequences of square-free polynomials
◮ We compute parallelPairsOfGcd(A, B), that is,
◮ This is done using the subproduct tree of b0, b1, . . . , bs−1 and
51 / 65
◮ Once we have a routine parallelPairsOfGcd(A, B), we easily
◮ Let us call it parallelMergeGcdFreeBases(f1, f2).
52 / 65
1:
2:
3:
4:
5:
6:
7:
8:
53 / 65
◮ Work = O(M(d) log3
2 d), GOOD thanks to fast arithmetic
◮ Span = O( d
n log2 2(n)), and
◮ Parallelism = O(n
log4
2 d
log2
2(n)), NOT GOOD, because almost similar to
i=1deg(inputi),
2 d log2 n).
54 / 65
◮ Better work than that of quadratic. ◮ Almost similar parallelism that of quadratic. ◮ But, implementation on multicore is more challenging,
◮ a negative impact on memory consumption and ◮ no cache-friendly algorithms in multicore is known. 55 / 65
56 / 65
◮ We have proposed and analyzed parallel algorithms for factor
◮ We have observed that cache complexity has major effects on
◮ Algorithm 5 with cache complexity O(n2/ZL) is more efficient
◮ Asymptotically fast polynomial arithmetic does not always pay
◮ Work in progress:
◮ Handling multi-precision (= big) integers efficiently. ◮ Supporting the multivariate case in the context of polynomial
57 / 65
58 / 65
59 / 65
Let W1(n), W2(n), W3(n) (resp. S1(n), S2(n), S3(n)) the work (resp. span) of Algorithms 1, 2 and 3 on input data of order n. Thus we have: W1(n) = 2W1(n/2) + W2(n) and S1(n) = S1(n/2) + S2(n). (1) W2(n) = W3(n) + O(n2) and S2(n) = S3(n) + O(n). (2) W3(n) = 4W3(n/2) + Θ(1) and S3(n) = S3(n/2) + Θ(1). (3) From Equation 3, we have: W3(n) = Θ(n2) and S3(n) = Θ(log2(n)). From Equation 2, we obtain W2(n) = Θ(n2) and S2(n) = Θ(n). Finally, from Equation 1, we deduce W1(n) = Θ(n2) and S1(n) = Θ(n).
60 / 65
Let Q(n) is the number of cache misses in Algorithm 3. Under the ideal cache model, there exists a constant α such that we have, Q(n) ≤
for n < αZ 4Q(n/2) + Θ(1)
The solution of the above recurrence is as follows: Q(n) ≤ 4Q(n/2) + Θ(1) ≤ 4[4Q(n/4) + Θ(1)] + Θ(1) . . . ≤ 4k Q(n/2k ) +
k−1
4j Θ(1) where 2k = n/αZ in base case. Under this base case, Q(n) ≤ (n/αZ)2[(αZ)2/L + αZ] + Θ((n/αZ)2) ≤ Θ(n2/L + n2/Z + n2/Z2) 61 / 65
Let us denote by W4(n), W5(n), W6(n), W7(n), W8(n) (resp. S4(n), S5(n), S6(n), S7(n), S8(n)) the work (resp. span) of Algorithms 4, 5, 6, 7 and 8 on input data of order n. Thus we have: W4(n) ≤ 2W4(n/2) + W5(n) and S4(n) ≤ S4(n/2) + S5(n). (4) Algorithm 5 also proceeds in a divide-and-conquer manner, dividing the input data into two parts and performing two recursive calls. W5(n) ≤
for n < C 4W5(n/2) + Θ(1)
(5) and S5(n) ≤
for n < C 2S5(n/2) + Θ(1)
(6) It is observed that W6(n), W7(n), W8(n) fit within Θ(n2), Θ(n), Θ(1), respectively. Moreover, since Algorithms 6, 7 and 8 are serial, we also have S6(n), S7(n), S8(n) within Θ(n2), Θ(n), Θ(1), respectively. Let k = ⌈log2(n/C)⌉. Then we have W5(n) ≤ O(4kC 2) = O(n2) and S5(n) ≤ O(2kC 2) = O(Cn). Therefore, from Relation (4), we deduce W4(n) ∈ O(n2) and S4(n) ∈ O(Cn).
62 / 65
The recursive structure of Algorithm 5 presents that there exists a positive constant α such that Q(n) satisfies the following relation: Q(n) ≤
for n < αZ 4Q(n/2) + Θ(1)
provided C < αZ. The above recurrence leads to the following inequality for all n ≥ 2: Q(n) ≤ 4Q(n/2) + Θ(1) ≤ 4[4Q(n/4) + Θ(1)] + Θ(1) . . . ≤ 4k Q(n/2k ) +
k−1
4j Θ(1) where k = ⌈log2(n/αZ)⌉. Since n/2k ≤ αZ, we deduce: Q(n) ≤ (n/αZ)2(αZ/L + 1) + Θ((n/αZ)2) = O(n2/ZL + n2/Z2). . 63 / 65
64 / 65
65 / 65