nbtheory_8cpp-source.html

著名的密码库Crypto++的文档 C++语言的杰作。程序员必备。

<a name="l00364"></a>00364         {
<a name="l00365"></a>00365                 assert(<a class="code" href="class_prime_sieve.html#e4a118cd164def0137049b041a40df5a">m_step</a>%2==0);
<a name="l00366"></a>00366                 Integer qFirst = (<a class="code" href="class_prime_sieve.html#505052e929448429944305dbbf8e3d9b">m_first</a>-<a class="code" href="class_prime_sieve.html#b088fe1cbe49ca396ffee53f94da0706">m_delta</a>) &gt;&gt; 1;
<a name="l00367"></a>00367                 Integer halfStep = <a class="code" href="class_prime_sieve.html#e4a118cd164def0137049b041a40df5a">m_step</a> &gt;&gt; 1;
<a name="l00368"></a>00368                 <span class="keywordflow">for</span> (<span class="keywordtype">unsigned</span> <span class="keywordtype">int</span> i = 0; i &lt; primeTableSize; ++i)
<a name="l00369"></a>00369                 {
<a name="l00370"></a>00370                         word16 p = primeTable[i];
<a name="l00371"></a>00371                         word16 stepInv = (word16)<a class="code" href="class_prime_sieve.html#e4a118cd164def0137049b041a40df5a">m_step</a>.<a class="code" href="class_integer.html#881f9c714ee42f35718725a43d4d7db3" title="calculate multiplicative inverse of *this mod n">InverseMod</a>(p);
<a name="l00372"></a>00372                         <a class="code" href="class_prime_sieve.html#ef7632203a02890568748c22506d6247">SieveSingle</a>(<a class="code" href="class_prime_sieve.html#43d578612bed4a6c9ef353987c0407f8">m_sieve</a>, p, <a class="code" href="class_prime_sieve.html#505052e929448429944305dbbf8e3d9b">m_first</a>, <a class="code" href="class_prime_sieve.html#e4a118cd164def0137049b041a40df5a">m_step</a>, stepInv);
<a name="l00373"></a>00373 
<a name="l00374"></a>00374                         word16 halfStepInv = 2*stepInv &lt; p ? 2*stepInv : 2*stepInv-p;
<a name="l00375"></a>00375                         <a class="code" href="class_prime_sieve.html#ef7632203a02890568748c22506d6247">SieveSingle</a>(<a class="code" href="class_prime_sieve.html#43d578612bed4a6c9ef353987c0407f8">m_sieve</a>, p, qFirst, halfStep, halfStepInv);
<a name="l00376"></a>00376                 }
<a name="l00377"></a>00377         }
<a name="l00378"></a>00378 }
<a name="l00379"></a>00379 
<a name="l00380"></a>00380 <span class="keywordtype">bool</span> FirstPrime(Integer &amp;p, <span class="keyword">const</span> Integer &amp;max, <span class="keyword">const</span> Integer &amp;equiv, <span class="keyword">const</span> Integer &amp;mod, <span class="keyword">const</span> <a class="code" href="class_prime_selector.html">PrimeSelector</a> *pSelector)
<a name="l00381"></a>00381 {
<a name="l00382"></a>00382         assert(!equiv.<a class="code" href="class_integer.html#d767ae81c89be3804da8785e132d2d1f">IsNegative</a>() &amp;&amp; equiv &lt; mod);
<a name="l00383"></a>00383 
<a name="l00384"></a>00384         Integer gcd = GCD(equiv, mod);
<a name="l00385"></a>00385         <span class="keywordflow">if</span> (gcd != <a class="code" href="class_integer.html#8c070592581bf6c2f928c72bfa1c1638" title="avoid calling constructors for these frequently used integers">Integer::One</a>())
<a name="l00386"></a>00386         {
<a name="l00387"></a>00387                 <span class="comment">// the only possible prime p such that p%mod==equiv where GCD(mod,equiv)!=1 is GCD(mod,equiv)</span>
<a name="l00388"></a>00388                 <span class="keywordflow">if</span> (p &lt;= gcd &amp;&amp; gcd &lt;= max &amp;&amp; IsPrime(gcd) &amp;&amp; (!pSelector || pSelector-&gt;<a class="code" href="class_prime_selector.html#44f47be2ff2e8ce1813cb0781d725181">IsAcceptable</a>(gcd)))
<a name="l00389"></a>00389                 {
<a name="l00390"></a>00390                         p = gcd;
<a name="l00391"></a>00391                         <span class="keywordflow">return</span> <span class="keyword">true</span>;
<a name="l00392"></a>00392                 }
<a name="l00393"></a>00393                 <span class="keywordflow">else</span>
<a name="l00394"></a>00394                         <span class="keywordflow">return</span> <span class="keyword">false</span>;
<a name="l00395"></a>00395         }
<a name="l00396"></a>00396 
<a name="l00397"></a>00397         <span class="keywordtype">unsigned</span> <span class="keywordtype">int</span> primeTableSize;
<a name="l00398"></a>00398         <span class="keyword">const</span> word16 * primeTable = GetPrimeTable(primeTableSize);
<a name="l00399"></a>00399 
<a name="l00400"></a>00400         <span class="keywordflow">if</span> (p &lt;= primeTable[primeTableSize-1])
<a name="l00401"></a>00401         {
<a name="l00402"></a>00402                 <span class="keyword">const</span> word16 *pItr;
<a name="l00403"></a>00403 
<a name="l00404"></a>00404                 --p;
<a name="l00405"></a>00405                 <span class="keywordflow">if</span> (p.<a class="code" href="class_integer.html#13ddbfd8e9729932c2a99b0dff530978">IsPositive</a>())
<a name="l00406"></a>00406                         pItr = std::upper_bound(primeTable, primeTable+primeTableSize, (word)p.<a class="code" href="class_integer.html#2e90d8f4c5a13e203b94f9abc24d733f" title="return equivalent signed long if possible, otherwise undefined">ConvertToLong</a>());
<a name="l00407"></a>00407                 <span class="keywordflow">else</span>
<a name="l00408"></a>00408                         pItr = primeTable;
<a name="l00409"></a>00409 
<a name="l00410"></a>00410                 <span class="keywordflow">while</span> (pItr &lt; primeTable+primeTableSize &amp;&amp; !(*pItr%mod == equiv &amp;&amp; (!pSelector || pSelector-&gt;<a class="code" href="class_prime_selector.html#44f47be2ff2e8ce1813cb0781d725181">IsAcceptable</a>(*pItr))))
<a name="l00411"></a>00411                         ++pItr;
<a name="l00412"></a>00412 
<a name="l00413"></a>00413                 <span class="keywordflow">if</span> (pItr &lt; primeTable+primeTableSize)
<a name="l00414"></a>00414                 {
<a name="l00415"></a>00415                         p = *pItr;
<a name="l00416"></a>00416                         <span class="keywordflow">return</span> p &lt;= max;
<a name="l00417"></a>00417                 }
<a name="l00418"></a>00418 
<a name="l00419"></a>00419                 p = primeTable[primeTableSize-1]+1;
<a name="l00420"></a>00420         }
<a name="l00421"></a>00421 
<a name="l00422"></a>00422         assert(p &gt; primeTable[primeTableSize-1]);
<a name="l00423"></a>00423 
<a name="l00424"></a>00424         <span class="keywordflow">if</span> (mod.<a class="code" href="class_integer.html#ed4bb7208a18b986ef3e1a7d92e06d1d">IsOdd</a>())
<a name="l00425"></a>00425                 <span class="keywordflow">return</span> FirstPrime(p, max, CRT(equiv, mod, 1, 2, 1), mod&lt;&lt;1, pSelector);
<a name="l00426"></a>00426 
<a name="l00427"></a>00427         p += (equiv-p)%mod;
<a name="l00428"></a>00428 
<a name="l00429"></a>00429         <span class="keywordflow">if</span> (p&gt;max)
<a name="l00430"></a>00430                 <span class="keywordflow">return</span> <span class="keyword">false</span>;
<a name="l00431"></a>00431 
<a name="l00432"></a>00432         <a class="code" href="class_prime_sieve.html">PrimeSieve</a> sieve(p, max, mod);
<a name="l00433"></a>00433 
<a name="l00434"></a>00434         <span class="keywordflow">while</span> (sieve.NextCandidate(p))
<a name="l00435"></a>00435         {
<a name="l00436"></a>00436                 <span class="keywordflow">if</span> ((!pSelector || pSelector-&gt;<a class="code" href="class_prime_selector.html#44f47be2ff2e8ce1813cb0781d725181">IsAcceptable</a>(p)) &amp;&amp; FastProbablePrimeTest(p) &amp;&amp; IsPrime(p))
<a name="l00437"></a>00437                         <span class="keywordflow">return</span> <span class="keyword">true</span>;
<a name="l00438"></a>00438         }
<a name="l00439"></a>00439 
<a name="l00440"></a>00440         <span class="keywordflow">return</span> <span class="keyword">false</span>;
<a name="l00441"></a>00441 }
<a name="l00442"></a>00442 
<a name="l00443"></a>00443 <span class="comment">// the following two functions are based on code and comments provided by Preda Mihailescu</span>
<a name="l00444"></a>00444 <span class="keyword">static</span> <span class="keywordtype">bool</span> ProvePrime(<span class="keyword">const</span> Integer &amp;p, <span class="keyword">const</span> Integer &amp;q)
<a name="l00445"></a>00445 {
<a name="l00446"></a>00446         assert(p &lt; q*q*q);
<a name="l00447"></a>00447         assert(p % q == 1);
<a name="l00448"></a>00448 
<a name="l00449"></a>00449 <span class="comment">// this is the Quisquater test. Numbers p having passed the Lucas - Lehmer test</span>
<a name="l00450"></a>00450 <span class="comment">// for q and verifying p &lt; q^3 can only be built up of two factors, both = 1 mod q,</span>
<a name="l00451"></a>00451 <span class="comment">// or be prime. The next two lines build the discriminant of a quadratic equation</span>
<a name="l00452"></a>00452 <span class="comment">// which holds iff p is built up of two factors (excercise ... )</span>
<a name="l00453"></a>00453 
<a name="l00454"></a>00454         Integer r = (p-1)/q;
<a name="l00455"></a>00455         <span class="keywordflow">if</span> (((r%q).Squared()-4*(r/q)).IsSquare())
<a name="l00456"></a>00456                 <span class="keywordflow">return</span> <span class="keyword">false</span>;
<a name="l00457"></a>00457 
<a name="l00458"></a>00458         <span class="keywordtype">unsigned</span> <span class="keywordtype">int</span> primeTableSize;
<a name="l00459"></a>00459         <span class="keyword">const</span> word16 * primeTable = GetPrimeTable(primeTableSize);
<a name="l00460"></a>00460 
<a name="l00461"></a>00461         assert(primeTableSize &gt;= 50);
<a name="l00462"></a>00462         <span class="keywordflow">for</span> (<span class="keywordtype">int</span> i=0; i&lt;50; i++) 
<a name="l00463"></a>00463         {
<a name="l00464"></a>00464                 Integer b = a_exp_b_mod_c(primeTable[i], r, p);
<a name="l00465"></a>00465                 <span class="keywordflow">if</span> (b != 1) 
<a name="l00466"></a>00466                         <span class="keywordflow">return</span> a_exp_b_mod_c(b, q, p) == 1;
<a name="l00467"></a>00467         }
<a name="l00468"></a>00468         <span class="keywordflow">return</span> <span class="keyword">false</span>;
<a name="l00469"></a>00469 }
<a name="l00470"></a>00470 
<a name="l00471"></a>00471 Integer MihailescuProvablePrime(<a class="code" href="class_random_number_generator.html" title="interface for random number generators">RandomNumberGenerator</a> &amp;rng, <span class="keywordtype">unsigned</span> <span class="keywordtype">int</span> pbits)
<a name="l00472"></a>00472 {
<a name="l00473"></a>00473         Integer p;
<a name="l00474"></a>00474         Integer minP = <a class="code" href="class_integer.html#de53248f5dbb520273a70856b975417c" title="return the integer 2**e">Integer::Power2</a>(pbits-1);
<a name="l00475"></a>00475         Integer maxP = <a class="code" href="class_integer.html#de53248f5dbb520273a70856b975417c" title="return the integer 2**e">Integer::Power2</a>(pbits) - 1;
<a name="l00476"></a>00476 
<a name="l00477"></a>00477         <span class="keywordflow">if</span> (maxP &lt;= Integer(s_lastSmallPrime).Squared())
<a name="l00478"></a>00478         {
<a name="l00479"></a>00479                 <span class="comment">// Randomize() will generate a prime provable by trial division</span>
<a name="l00480"></a>00480                 p.<a class="code" href="class_integer.html#0f0574b9cae3cddf62c155da93085f0d">Randomize</a>(rng, minP, maxP, <a class="code" href="class_integer.html#9b4088ac01abf76b9ba60060abccb7a3fe686f55e5b6768b20009a12522bd0d9">Integer::PRIME</a>);
<a name="l00481"></a>00481                 <span class="keywordflow">return</span> p;
<a name="l00482"></a>00482         }
<a name="l00483"></a>00483 
<a name="l00484"></a>00484         <span class="keywordtype">unsigned</span> <span class="keywordtype">int</span> qbits = (pbits+2)/3 + 1 + rng.<a class="code" href="class_random_number_generator.html#c06c9f5e66f35d643e4192f231b8a4cd" title="generate a random 32 bit word in the range min to max, inclusive">GenerateWord32</a>(0, pbits/36);
<a name="l00485"></a>00485         Integer q = MihailescuProvablePrime(rng, qbits);
<a name="l00486"></a>00486         Integer q2 = q&lt;&lt;1;
<a name="l00487"></a>00487 
<a name="l00488"></a>00488         <span class="keywordflow">while</span> (<span class="keyword">true</span>)
<a name="l00489"></a>00489         {
<a name="l00490"></a>00490                 <span class="comment">// this initializes the sieve to search in the arithmetic</span>
<a name="l00491"></a>00491                 <span class="comment">// progression p = p_0 + \lambda * q2 = p_0 + 2 * \lambda * q,</span>