Is there anyway to parallelize the PowerMod function?
Here is my Left-To-Right modular exponentation:
AfshinPowerMod[a_, b_, m_] := (Output = 1;
Do[If[n == 1, Output = Mod[Output*Output*a, m],
Output = PowerMod[Output, 2, m]], {n, IntegerDigits[b, 2]}];
Return[Output])
It computes roughly in the same time as original PowerMod for 10K digits numbers.
Update
(Since for prime testing, 2 is the lowest base, I've used it for faster computation)
Timing Result (Mathematica 8):
In[1]:= Total[DigitCount[2^100000 + 1]]
Out[1]= 30103
In[2]:= AbsoluteTiming[AfshinPowerMod[2, 2^100000 + 1, 2^100000 + 1]]
Out[2]= {91.839778,}
In[3]:= AbsoluteTiming[PowerMod[2, 2^100000 + 1, 2^100000 + 1]]
Out[3]= {92.9851312,}
As per J.M. request for RandomPrime (Based on the generated random prime timing differs, but yet it is roughly as fast as PowerMod)
In[41]:= prime = RandomPrime[10^2000]
In[43]:= AbsoluteTiming[PowerMod[2, prime, prime]]
Out[43]= {0.1406312, 2}
In[42]:= AbsoluteTiming[AfshinPowerMod[2, prime, prime]]
Out[42]= {0.1562399, 2}
Fold[]instead ofDo[]in your implementation:pmod[a_, b_, m_] := Fold[Mod[a^#2 #1^2, m] &, 1, IntegerDigits[b, 2]], though I'm not sure if this approach is still any better than the built-in one... – J. M.'s missing motivation Nov 22 '12 at 16:50ParallelMaporParallelTable to runPowerModon your 100k numbers over all your cores. – s0rce Nov 22 '12 at 19:10RandomPrime[]as the modulus and base... – J. M.'s missing motivation Nov 23 '12 at 10:03RandomPrime[10^1000]. PowerMod remains faster than AfshinPowerMod. Finding the random primes takes almost all the CPU time, so I calculated them first then wrappedAbsoluteTimingaround the two mod functions. – KennyColnago Nov 24 '12 at 02:50