We present a randomized parallel algorithm that computes the greatest common
divisor of two integers of n bits in length with probability 1-o(1) that takes
O(n loglog n / log n) expected time using n^{6+\epsilon} processors on the EREW
PRAM parallel model of computation. We believe this to be the first randomized
sublinear time algorithm on the EREW PRAM for this problem

Bach

Beame

Bernstein

Borodin

Chor

Crandall

Davida

Greenlaw

Hildebrand

Jonathan P. Sorenson

Kannan

Karp

Meyer Eikenberry

Parsell

Purdom

Schönhage

Sedjelmaci

Sorenson

Tenenbaum

English

arXiv

We present a randomized parallel algorithm that computes the greatest common divisor of two integers of n bits in length with probability 1−o(1) that takes O(n log logn/ logn) time using O(n6 + ) processors for any \u3e 0 on the EREW PRAM parallel model of computation. The algorithm either gives a correct answer or reports failure. We believe this to be the first randomized sublinear time algorithm on the EREW PRAM for this problem

Sorenson, Jonathan P.

Digital Commons @ Butler University

Butler UniversityDigital Commons @ Butler UniversityScholarship and Professional Work - LAS College of Liberal Arts & Sciences2010A Randomized Sublinear Time Parallel GCDAlgorithm for the EREW PRAMJonathan P. SorensonButler University, jsorenso@butler.eduFollow this and additional works at: http://digitalcommons.butler.edu/facsch_papersPart of the Theory and Algorithms CommonsThis Article is brought to you for free and open access by the College of Liberal Arts & Sciences at Digital Commons @ Butler University. It has beenaccepted for inclusion in Scholarship and Professional Work - LAS by an authorized administrator of Digital Commons @ Butler University. For moreinformation, please contact fgaede@butler.edu.Recommended CitationSorenson, Jonathan P., "A Randomized Sublinear Time Parallel GCD Algorithm for the EREW PRAM" Information Processing Letters /(2010): 198-201.Available at http://digitalcommons.butler.edu/facsch_papers/81Permission to post this publication in our archive was granted by the copyright holder, Elsevier (http://www.elsevier.com/wps/find/journaldescription.cws_home/505612/description#description). This copy should be used for educational and research purposes only.  The original publication appeared at:  Sorenson, Jonathan P.  “A Randomized Sublinear Time Parallel GCD Algorithm for the EREW PRAM.”  Information Processing Letters 110 (2010) 198-201.     DOI: Not Available  Author's personal copyInformation Processing Letters 110 (2010) 198–201Contents lists available at ScienceDirectInformation Processing Letterswww.elsevier.com/locate/iplA randomized sublinear time parallel GCD algorithm for the EREW PRAMJonathan P. SorensonComputer Science and Software Engineering, Butler University, Indianapolis, IN, USAa r t i c l e i n f o a b s t r a c tArticle history:Received 1 July 2009Received in revised form 29 October 2009Accepted 19 December 2009Available online 28 December 2009Communicated by M. YamashitaKeywords:Parallel algorithmsRandomized algorithmsAlgorithm analysisGreatest common divisorNumber theoretic algorithmsSmooth numbersWe present a randomized parallel algorithm that computes the greatest common divisor oftwo integers of n bits in length with probability 1−o(1) that takes O (n log logn/ logn) timeusing O (n6+) processors for any  > 0 on the EREW PRAM parallel model of computation.The algorithm either gives a correct answer or reports failure.We believe this to be the ﬁrst randomized sublinear time algorithm on the EREW PRAMfor this problem.© 2009 Elsevier B.V. All rights reserved.1. IntroductionThe parallel complexity of computing integer greatestcommon divisors is an open problem (see [5]), and no newcomplexity results have been published since the early1990s. This problem is not known to be either P-completeor in N C [10,13,16].The ﬁrst sublinear time parallel algorithm that uses apolynomial number of processors is due to Kannan, Miller,and Rudolph [12]. Adleman and Kompella [1] presented arandomized algorithm that runs in polylog time, but usesa superpolynomial, yet subexponential number of pro-cessors. The fastest currently known algorithm is due toChor and Goldreich [6] which takes O (n/ logn) time us-ing O (n1+) processors. See also [22], and Sedjelmaci [19]who showed a clear way to do extended GCDs in the samecomplexity bounds. However, all of these algorithms usethe concurrent-read concurrent-write (CRCW) parallel RAM(PRAM) model of computation.The algorithms of Chor and Goldreich [6] and the au-thor [22] can be readily modiﬁed for the weaker concur-E-mail address: sorenson@butler.edu.URL: http://www.butler.edu/~sorenson.rent-read exclusive-write (CREW) PRAM to obtain runningtimes of O (n log logn/ logn) using a polynomial numberof processors. And of course one can take a CRCW PRAMalgorithm and emulate it on an exclusive-read exclusive-write (EREW) PRAM at a cost of a factor of O (logn) in therunning time, giving linear time algorithms for the EREWPRAM using a polynomial number of processors. For moreinformation, see [2,7,18,20,21].In this paper, we present what we believe is the ﬁrstsublinear time, polynomial processor EREW PRAM algo-rithm for computing greatest common divisors. Note thatthe EREW PRAM is weaker than the CREW or CRCW PRAMmodels of parallel computation. We do make use of ran-dom numbers in a fundamental way.Theorem 1.1. There exists a randomized algorithm to computethe greatest common divisor of two integers of total length n inbinary with probability 1 − o(1) in O (n log logn/ logn) timeusing a polynomial number of processors on the EREW PRAM.In the next section we describe our algorithm, and inSection 3 we prove correctness, give a complexity analysis,and ﬂesh out the details of the algorithm. We conclude inSection 4 with a simple result on the relative density of0020-0190/$ – see front matter © 2009 Elsevier B.V. All rights reserved.doi:10.1016/j.ipl.2009.12.008Author's personal copyJ.P. Sorenson / Information Processing Letters 110 (2010) 198–201 199integers with large polynomially smooth divisors, which isneeded for the analysis of the algorithm.2. Algorithm descriptionDeﬁne the inputs as u, v of total length n in bi-nary. Let B , our small prime bound, be deﬁned as B =B(n) := n2. A larger value for B can be chosen, so long aslog B = o(n), but correctness would be compromised if Bwere signiﬁcantly smaller (see Section 3.1).1. Find a list of primes up to B . Also, for each primep 6 B , compute and save pe for e = 1 . . . n/ log2 p.2. Remove and save common prime factors of u, v thatare 6 B , and let u0, v0 denote these modiﬁed inputs.WLOG we assume u0 > v0.3. Main Loop. Repeat while ui vi = 0. Here i indicates thecurrent loop iteration, starting at i = 0.(a) For j := 1 to 2B logn in parallel do:i. Choose aij uniformly at random from 1 . . . vi−1.ii. Compute ri j := aijui mod vi .iii. Compute si j as ri j with all prime factors p 6 Bremoved.(We elaborate on how to do this below.)(b) Find si := min j{si j}. Let jmin denote the value of jfor which si = si j , and for later reference, let ai =aijmin .(c) ui+1 := vi ; vi+1 := si .4. ui+vi is, with probability 1−o(1), equal to gcd(u0, v0)(as we show below). If we err, it is by including spu-rious factors that do not belong, so verify that ui + vievenly divides both u0, v0, and if not, report an error.Otherwise, include any saved common prime factorsfound in step 2 above, and the algorithm is complete.3. Algorithm analysisIn this section we prove correctness, and compute theparallel complexity of our algorithm from the previous sec-tion.3.1. CorrectnessNote that in Steps 2 and 4 we handle any prime divi-sors 6 B of the gcd(u, v), so WLOG we can assume eithergcd(u, v) = 1 or gcd(u, v) > B .At iteration i of the main loop, we perform the trans-formation(ui, vi) → (vi, si).Since si is equal to aiui mod vi , ignoring factors below B ,this transformation will only fail to preserve the greatestcommon divisor if ai and vi share a common factor. Fur-thermore, this common factor must be composed only ofprimes exceeding B . Since ai is chosen uniformly at ran-dom, the probability ai and vi share a prime factor largerthat B is at most∑p|vi , p>B1p6∑p|vi , p>B1B6 logB viB= O(1n logn).As we will see below, with high probability, the number ofmain loop iterations is o(n). Thus, the probability that anyof the ai values introduces a spurious factor is o(1).Note that in [23], a similar, but not identical, transfor-mation was analyzed. It was observed that in practice, withno removal of small prime divisors, the expected numberof bits contributed by spurious factors was constant permain loop iteration.3.2. Runtime analysisFirst we calculate the number of main loop iterations,and then we describe how each iteration can be computedin O (logn) time using a polynomial number of processors.3.2.1. Main loop iterationsLet W := 0.5(log B)2/ log log B . Then by Theorem 4.2,which we prove in the next section, the length of si j issmaller than ri j < vi by at least Θ(W ) bits with probabil-ity at least 1/B . (Note that we chose 0.5 to get a clean 1/Bprobability – other choices for the constant can be made towork with the right adjustments.)So, the probability all 2B logn choices for j fail to havelog si j 6 log vi − W is(1− 1B)2B logn= O(1n2).So, with probability 1− O (1/n2), log si 6 log vi − W .We remove roughly (log B)2/ log log B bits each mainloop iteration. Thus, the number of main loop iterationsis O (n log log B/(log B)2) = o(n). The probability that anyone loop iteration fails to remove the needed Θ(W ) bits isO (1/n), so the probability we exceed this number of mainloop iterations and terminate without computing an an-swer is o(1).3.2.2. Computation cost and algorithm detailsUnless stated otherwise, cost is given for the EREWPRAM. For a brief overview of the cost of parallel arith-metic, see [22, Section 6.2].Step 1. We can ﬁnd the primes 6 B in O (log B) time usingO (B) processors (see [24]). For each prime p 6 B ande 6 n, we can compute pe in at most O (logn) mul-tiplications, each of which takes O (log B) time usingB1+o(1) processors [17]. See also [16, Theorem 12.2].As there are O (B/ log B) primes, this is O (logn log B)time using nB2+o(1) processors.Step 2. For each prime p and exponent e, we assign agroup of processors to see if pe divides u but pe+1does not. Division takes O (logn) time using n1+ pro-cessors for any  > 0 using the algorithm of Beame,Cook, and Hoover [3], giving a total processor count ofO (n2+ B/ log B).The result is a vector of the form [pekk ] that liststhe primes dividing u with maximal exponents. Sincethere are at most n/ log B integers in the vector > 1,and they total at most n bits (their product is 6 u), theiterated product algorithm of [3] can take their prod-uct in O (logn) time using n1+ processors. Dividing uby this product can be done at the same cost.Author's personal copy200 J.P. Sorenson / Information Processing Letters 110 (2010) 198–201We repeat this for v , and obtain a similar vector.We combine these two vectors using a minimum op-eration, and take the product of the entries, to obtainthe shared prime power divisors of u, v which mustbe saved for Step 4.The total cost of this step is O (logn) time usingO (n2+ B/ log B) processors.See also [8] and references therein.Step 3. Checking for zero takes O (logn) time using O (n)processors.Step 3(a)i. For each j, choosing an n-bit number at ran-dom takes constant time using O (n) processors. Wereduce it modulo vi in O (logn) time using n1+ pro-cessors [3].Step 3(a)ii. This is simply a multiplication and a division,again taking O (logn) time using n1+ processors.Step 3(a)iii. Here we use the same method as describedin Step 2 above. This is O (logn) time using O (n2+ B/log B) processors for each j.Step 3(a). And so, the total cost of this parallel step isO (log B) time using O (n2+(logn)B2/ log B) proces-sors.Step 3(b). This can be done in O (log(B logn)) = O (log B)time using O (B logn) processors.Step 3(c). This takes constant time using O (n) processors.We conclude that the cost of one main loop iterationis O (log B) time using O (n2+(logn)B2/ log B) proces-sors. Step 3(a)iii is the bottleneck.Earlier we showed that the number of iterations isO (n log log B/(log B)2), for a total time of O (n log log B/log B) for all iterations of the main loop.Step 4. This is an addition, a division, and a multiplicationusing the results from Step 2; O (logn) time using n1+processors.Clearly, the bottleneck of the algorithm is Step 3(a). Theoverall complexity isO(n log log Blog B)= O(n log lognlogn)time, andO(n2+(logn) B2log B)= O (n6+) processors,where  > 0,for the EREW PRAM. This completes our proof of Theo-rem 1.1.One could take B to be superpolynomial in n; for ex-ample, if B = exp[√n ] we can obtain a running time ofroughly√n using exp[O (√n )] processors. Similar resultscould be obtained from some of the CRCW PRAM algo-rithms mentioned in the introduction by porting them tothe EREW PRAM and setting parameters appropriately.We can also obtain an O (n/ logn) running time on therandomized CRCW PRAM; see [4] for how to perform thenecessary main loop operations in O (logn/ log logn) timevia the explicit Chinese remainder theorem. See also [8].It would be interesting to see if this algorithm can bemodiﬁed to compute Jacobi symbols quickly in parallel. See[9] and references therein.4. Numbers with smooth divisorsLet P (n) denote the largest prime divisor of n. IfP (n)6 y we say that n is y-smooth. LetΨ (x, y) = #{n6 x: P (n)6 y},the number of integers 6 x that are y-smooth. Letu = u(x, y) := log x/ log y. We will make use of the fol-lowing lemma.Lemma 4.1. (See [11, Corollary 1.3].) Let  > 0 and assumeu < y1− . ThenΨ (x, y) = xu−u(1+o(1)) for x > y > 2.Note that the o(1) here tends to zero for large u,and the implied constant depends on  . Better results areknown, but this suﬃces for our purposes. For additionalreferences see [11,25], and for references on approxima-tion algorithms for Ψ (x, y) see [14].We recall the deﬁnition of Hk , the kth harmonic num-ber asHk =k∑i=11i.It is well known that Hk = logk + γ + O (1/k), whereγ = 0.57721 . . . is Euler’s constant (for example, see [15,4.5.4]).Fix a constant c > 0. Deﬁne B(x) to be a strictly in-creasing function of x, but with log B(x) = o(log x). (We areprimarily interested in B(x) polynomial in log x.) DeﬁneW (x) := c · (log B(x))2log log B(x),F (x) := #{n6 x: n =my, P (m)6 B(x), logm>W (x)}.In other words, F (x) counts integers n 6 x where n has aB(x)-smooth divisor that is > expW (x).Theorem 4.2. Let  > 0. For suﬃciently large x we haveF (x)> xB(x)c(1+).Proof. Choose δ > 0 such that (1 + δ)3 < 1 +  . From thedeﬁnition, we haveF (x) =x/exp[W (x)]∑y=1Ψ(xy, B(x)).First, we limit the range of summation to obtain the lowerboundF (x)>x/exp[W (x)]∑y=x/(exp[(1+δ)W (x)])Ψ(xy, B(x)).Next, we apply Lemma 4.1. We also observe that u−u(1+o(1))> u−(1+δ)u for large u, and for a lower bound, we can ﬁxu at its largest value on the interval of summation, namelyu = u(x) = (1+ δ)W (x)/ log B(x). This gives usAuthor's personal copyJ.P. Sorenson / Information Processing Letters 110 (2010) 198–201 201F (x)>x/exp[W (x)]∑y=x/(exp[(1+δ)W (x)])xy· u−(1+δ)u.Using∑ba 1/t = Hb − Ha > (1 − δ) log(b/a) for suﬃcientlylarge a, we obtain thatF (x)> x · δ(1− δ)W (x) · u−(1+δ)u> x · u−(1+δ)uas W is a strictly increasing function of x for large x. Nextwe plug in for u as follows:log(u−(1+δ)u)= −(1+ δ)u logu= −(1+ δ) (1+ δ)W (x)log B(x)log((1+ δ)W (x)log B(x))= −(1+ δ) (1+ δ) c log B(x)log log B(x)log((1+ δ) c log B(x)log log B(x))>−c (1+ δ)3 log B(x)for x suﬃciently large. We now haveF (x)> x · B(x)−c(1+δ)3for suﬃciently large x. Only a lower bound is needed for our purposes, but onecan obtain an upper bound on F (x) of similar shape usingthe same general methods.AcknowledgementsWe would like to thank the referee for many help-ful comments that improved this paper, and the HolcombAwards Committee at Butler University for their support.References[1] L.M. Adleman, K. Kompella, Using smoothness to achieve parallelism,in: 20th Annual ACM Symposium on Theory of Computing, 1988,pp. 528–538.[2] Eric Bach, Jeffrey O. Shallit, Algorithmic Number Theory, vol. 1, MITPress, 1996.[3] P.W. Beame, S.A. Cook, H.J. Hoover, Log depth circuits for division andrelated problems, SIAM J. Comput. 15 (1986) 994–1003.[4] Daniel J. Bernstein, Jonathan P. Sorenson, Modular exponentiationvia the explicit Chinese remainder theorem, Math. Comp. 76 (257)(2007) 443–454.[5] A. Borodin, J. von zur Gathen, J. Hopcroft, Fast parallel matrix andGCD computations, Inform. and Control 52 (1982) 241–256.[6] B. Chor, O. Goldreich, An improved parallel algorithm for integerGCD, Algorithmica 5 (1990) 1–10.[7] R. Crandall, C. Pomerance, Prime Numbers, a Computational Perspec-tive, Springer, 2001.[8] George Davida, Bruce Litow, Guangwu Xu, Fast arithmetics using Chi-nese remaindering, Inform. Process. Lett. 109 (13) (2009) 660–662.[9] Shawna M. Meyer Eikenberry, Jonathan P. Sorenson, Eﬃcient algo-rithms for computing the Jacobi symbol, J. Symbolic Comput. 26 (4)(1998) 509–523.[10] R. Greenlaw, H.J. Hoover, W.L. Ruzzo, Limits to Parallel Computation,Oxford University Press, 1995.[11] A. Hildebrand, G. Tenenbaum, Integers without large prime factors,J. Théor. Nombres Bordeaux 5 (1993) 411–484.[12] R. Kannan, G. Miller, L. Rudolph, Sublinear parallel algorithm forcomputing the greatest common divisor of two integers, SIAM J.Comput. 16 (1) (1987) 7–16.[13] R. Karp, V. Ramachandran, Parallel algorithms for shared-memorymachines, in: J. van Leeuwen (Ed.), Algorithms and Complexity, in:Handbook of Theoretical Computer Science, vol. A, Elsevier and MITPress, 1990.[14] Scott Parsell, Jonathan P. Sorenson, Fast bounds on the distributionof smooth numbers, in: Florian Hess, Sebastian Pauli, Michael Pohst(Eds.), Proceedings of the 7th International Symposium on Algorith-mic Number Theory (ANTS-VII), in: LNCS, vol. 4076, Springer, Berlin,Germany, ISBN 3-540-36075-1, July 2006, pp. 168–181.[15] P. Purdom Jr., C. Brown, The Analysis of Algorithms, Holt, Rinehart,and Winston, 1985.[16] J.H. Reif (Ed.), Synthesis of Parallel Algorithms, Morgan Kaufman, SanMateo, CA, 1993.[17] A. Schönhage, V. Strassen, Schnelle Multiplikation großer Zahlen,Computing (Arch. Elektron. Rechnen) 7 (1971) 281–292, MR 45#1431.[18] Sidi Mohamed Sedjelmaci, On a parallel extended euclidean algo-rithm, in: AICCSA, IEEE Computer Society, 2001, pp. 235–241.[19] Sidi Mohamed Sedjelmaci, A parallel extended GCD algorithm, J. Dis-crete Algorithms 6 (3) (2008) 526–538.[20] Sidi Mohamed Sedjelmaci, On a parallel Lehmer–Euclid GCD algo-rithm, in: Proceedings of the 2001 International Symposium on Sym-bolic and Algebraic Computation, ACM, New York, 2001, pp. 303–308.[21] Sidi Mohamed Sedjelmaci, A modular reduction for GCD computa-tion, J. Comput. Appl. Math. 162 (1) (2004) 17–31.[22] Jonathan P. Sorenson, Two fast GCD algorithms, J. Algorithms 16(1994) 110–144.[23] Jonathan P. Sorenson, An analysis of the generalized binary gcd al-gorithm, in: Alf van der Poorten, Andreas Stein (Eds.), High Primesand Misdemeanors: Lectures in Honour of the 60th Birthday of HughCowie Williams, Banff, Alberta, Canada, 2004, pp. 327–340.[24] Jonathan P. Sorenson, Ian Parberry, Two fast parallel prime numbersieves, Inform. and Comput. 144 (1) (1994) 115–130.[25] Gérald Tenenbaum, Introduction to Analytic and Probabilistic Num-ber Theory, English edition, Cambridge Stud. Adv. Math., vol. 46,Cambridge University Press, 1995.

A Randomized Sublinear Time Parallel GCD Algorithm for the EREW PRAM

Crossref

A randomized sublinear time parallel GCD algorithm for the EREW PRAM

CiteSeerX

http://digitalcommons.butler.edu/cgi/viewcontent.cgi?article=1081&context=facsch_papers

A Randomized Sublinear Time Parallel GCD Algorithm for the EREW PRAM

Abstract

Similar works

Full text

Available Versions

Digital Commons @ Butler University

Crossref

CiteSeerX