nbtheory_8cpp-source.html

Crypto++是一个非常强大的密码学库,主要是功能全

00396         {
00397                 <span class="keywordflow">for</span> (<span class="keywordtype">unsigned</span> <span class="keywordtype">int</span> i = 0; i &lt; primeTableSize; ++i)
00398                         SieveSingle(m_sieve, primeTable[i], m_first, m_step, m_step.InverseMod(primeTable[i]));
00399         }
00400         <span class="keywordflow">else</span>
00401         {
00402                 assert(m_step%2==0);
00403                 Integer qFirst = (m_first-m_delta) &gt;&gt; 1;
00404                 Integer halfStep = m_step &gt;&gt; 1;
00405                 <span class="keywordflow">for</span> (<span class="keywordtype">unsigned</span> <span class="keywordtype">int</span> i = 0; i &lt; primeTableSize; ++i)
00406                 {
00407                         word p = primeTable[i];
00408                         word stepInv = m_step.<a class="code" href="class_integer.html#_integerz49_7">InverseMod</a>(p);
00409                         SieveSingle(m_sieve, p, m_first, m_step, stepInv);
00410 
00411                         word halfStepInv = 2*stepInv &lt; p ? 2*stepInv : 2*stepInv-p;
00412                         SieveSingle(m_sieve, p, qFirst, halfStep, halfStepInv);
00413                 }
00414         }
00415 }
00416 
00417 <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> PrimeSelector *pSelector)
00418 {
00419         assert(!equiv.<a class="code" href="class_integer.html#_integerz41_10">IsNegative</a>() &amp;&amp; equiv &lt; mod);
00420 
00421         Integer gcd = GCD(equiv, mod);
00422         <span class="keywordflow">if</span> (gcd != <a class="code" href="class_integer.html#_integerz37_11">Integer::One</a>())
00423         {
00424                 <span class="comment">// the only possible prime p such that p%mod==equiv where GCD(mod,equiv)!=1 is GCD(mod,equiv)</span>
00425                 <span class="keywordflow">if</span> (p &lt;= gcd &amp;&amp; gcd &lt;= max &amp;&amp; IsPrime(gcd))
00426                 {
00427                         p = gcd;
00428                         <span class="keywordflow">return</span> <span class="keyword">true</span>;
00429                 }
00430                 <span class="keywordflow">else</span>
00431                         <span class="keywordflow">return</span> <span class="keyword">false</span>;
00432         }
00433 
00434         BuildPrimeTable();
00435 
00436         <span class="keywordflow">if</span> (p &lt;= primeTable[primeTableSize-1])
00437         {
00438                 word *pItr;
00439 
00440                 --p;
00441                 <span class="keywordflow">if</span> (p.<a class="code" href="class_integer.html#_integerz41_12">IsPositive</a>())
00442                         pItr = std::upper_bound(primeTable, primeTable+primeTableSize, (word)p.<a class="code" href="class_integer.html#_integerz41_1">ConvertToLong</a>());
00443                 <span class="keywordflow">else</span>
00444                         pItr = primeTable;
00445 
00446                 <span class="keywordflow">while</span> (pItr &lt; primeTable+primeTableSize &amp;&amp; *pItr%mod != equiv)
00447                         ++pItr;
00448 
00449                 <span class="keywordflow">if</span> (pItr &lt; primeTable+primeTableSize)
00450                 {
00451                         p = *pItr;
00452                         <span class="keywordflow">return</span> p &lt;= max;
00453                 }
00454 
00455                 p = primeTable[primeTableSize-1]+1;
00456         }
00457 
00458         assert(p &gt; primeTable[primeTableSize-1]);
00459 
00460         <span class="keywordflow">if</span> (mod.<a class="code" href="class_integer.html#_integerz41_15">IsOdd</a>())
00461                 <span class="keywordflow">return</span> FirstPrime(p, max, CRT(equiv, mod, 1, 2, 1), mod&lt;&lt;1, pSelector);
00462 
00463         p += (equiv-p)%mod;
00464 
00465         <span class="keywordflow">if</span> (p&gt;max)
00466                 <span class="keywordflow">return</span> <span class="keyword">false</span>;
00467 
00468         PrimeSieve sieve(p, max, mod);
00469 
00470         <span class="keywordflow">while</span> (sieve.NextCandidate(p))
00471         {
00472                 <span class="keywordflow">if</span> ((!pSelector || pSelector-&gt;IsAcceptable(p)) &amp;&amp; FastProbablePrimeTest(p) &amp;&amp; IsPrime(p))
00473                         <span class="keywordflow">return</span> <span class="keyword">true</span>;
00474         }
00475 
00476         <span class="keywordflow">return</span> <span class="keyword">false</span>;
00477 }
00478 
00479 <span class="comment">// the following two functions are based on code and comments provided by Preda Mihailescu</span>
00480 <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)
00481 {
00482         assert(p &lt; q*q*q);
00483         assert(p % q == 1);
00484 
00485 <span class="comment">// this is the Quisquater test. Numbers p having passed the Lucas - Lehmer test</span>
00486 <span class="comment">// for q and verifying p &lt; q^3 can only be built up of two factors, both = 1 mod q,</span>
00487 <span class="comment">// or be prime. The next two lines build the discriminant of a quadratic equation</span>
00488 <span class="comment">// which holds iff p is built up of two factors (excercise ... )</span>
00489 
00490         Integer r = (p-1)/q;
00491         <span class="keywordflow">if</span> (((r%q).Squared()-4*(r/q)).IsSquare())
00492                 <span class="keywordflow">return</span> <span class="keyword">false</span>;
00493 
00494         assert(primeTableSize &gt;= 50);
00495         <span class="keywordflow">for</span> (<span class="keywordtype">int</span> i=0; i&lt;50; i++) 
00496         {
00497                 Integer b = a_exp_b_mod_c(primeTable[i], r, p);
00498                 <span class="keywordflow">if</span> (b != 1) 
00499                         <span class="keywordflow">return</span> a_exp_b_mod_c(b, q, p) == 1;
00500         }
00501         <span class="keywordflow">return</span> <span class="keyword">false</span>;
00502 }
00503 
00504 Integer MihailescuProvablePrime(<a class="code" href="class_random_number_generator.html">RandomNumberGenerator</a> &amp;rng, <span class="keywordtype">unsigned</span> <span class="keywordtype">int</span> pbits)
00505 {
00506         Integer p;
00507         Integer minP = <a class="code" href="class_integer.html#_integerz37_13">Integer::Power2</a>(pbits-1);
00508         Integer maxP = <a class="code" href="class_integer.html#_integerz37_13">Integer::Power2</a>(pbits) - 1;
00509 
00510         <span class="keywordflow">if</span> (maxP &lt;= Integer(lastSmallPrime).Squared())
00511         {
00512                 <span class="comment">// Randomize() will generate a prime provable by trial division</span>
00513                 p.<a class="code" href="class_integer.html#_integerz43_10">Randomize</a>(rng, minP, maxP, Integer::PRIME);
00514                 <span class="keywordflow">return</span> p;
00515         }
00516 
00517         <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#_x917_r_n_ga3">GenerateWord32</a>(0, pbits/36);
00518         Integer q = MihailescuProvablePrime(rng, qbits);
00519         Integer q2 = q&lt;&lt;1;
00520 
00521         <span class="keywordflow">while</span> (<span class="keyword">true</span>)
00522         {
00523                 <span class="comment">// this initializes the sieve to search in the arithmetic</span>
00524                 <span class="comment">// progression p = p_0 + \lambda * q2 = p_0 + 2 * \lambda * q,</span>
00525                 <span class="comment">// with q the recursively generated prime above. We will be able</span>
00526                 <span class="comment">// to use Lucas tets for proving primality. A trick of Quisquater</span>
00527                 <span class="comment">// allows taking q &gt; cubic_root(p) rather then square_root: this</span>
00528                 <span class="comment">// decreases the recursion.</span>
00529 
00530                 p.<a class="code" href="class_integer.html#_integerz43_10">Randomize</a>(rng, minP, maxP, Integer::ANY, 1, q2);
00531                 PrimeSieve sieve(p, STDMIN(p+PrimeSearchInterval(maxP)*q2, maxP), q2);
00532 
00533                 <span class="keywordflow">while</span> (sieve.NextCandidate(p))
00534                 {
00535                         <span class="keywordflow">if</span> (FastProbablePrimeTest(p) &amp;&amp; ProvePrime(p, q))
00536                                 <span class="keywordflow">return</span> p;
00537                 }
00538         }
00539 
00540         <span class="comment">// not reached</span>
00541         <span class="keywordflow">return</span> p;
00542 }
00543 
00544 Integer MaurerProvablePrime(<a class="code" href="class_random_number_generator.html">RandomNumberGenerator</a> &amp;rng, <span class="keywordtype">unsigned</span> <span class="keywordtype">int</span> bits)
00545 {
00546         <span class="keyword">const</span> <span class="keywordtype">unsigned</span> smallPrimeBound = 29, c_opt=10;
00547         Integer p;
00548 
00549         BuildPrimeTable();
00550         <span class="keywordflow">if</span> (bits &lt; smallPrimeBound)
00551         {
00552                 <span class="keywordflow">do</span>
00553                         p.<a class="code" href="class_integer.html#_integerz43_10">Randomize</a>(rng, Integer::Power2(bits-1), Integer::Power2(bits)-1, Integer::ANY, 1, 2);
00554                 <span class="keywordflow">while</span> (TrialDivision(p, 1 &lt;&lt; ((bits+1)/2)));
00555         }
00556         <span class="keywordflow">else</span>
00557         {
00558                 <span class="keyword">const</span> <span class="keywordtype">unsigned</span> margin = bits &gt; 50 ? 20 : (bits-10)/2;
00559                 <span class="keywordtype">double</span> relativeSize;
00560                 <span class="keywordflow">do</span>
00561                         relativeSize = pow(2.0, <span class="keywordtype">double</span>(rng.<a class="code" href="class_random_number_generator.html#_x917_r_n_ga3">GenerateWord32</a>())/0xffffffff - 1);
00562                 <span class="keywordflow">while</span> (bits * relativeSize &gt;= bits - margin);
00563 
00564                 Integer a,b;
00565                 Integer q = MaurerProvablePrime(rng, <span class="keywordtype">unsigned</span>(bits*relativeSize));
00566                 Integer I = <a class="code" href="class_integer.html#_integerz37_13">Integer::Power2</a>(bits-2)/q;
00567                 Integer I2 = I &lt;&lt; 1;
00568                 <span class="keywordtype">unsigned</span> <span class="keywordtype">int</span> trialDivisorBound = (<span class="keywordtype">unsigned</span> <span class="keywordtype">int</span>)STDMIN((<span class="keywordtype">unsigned</span> <span class="keywordtype">long</span>)primeTable[primeTableSize-1], (<span class="keywordtype">unsigned</span> <span class="keywordtype">long</span>)bits*bits/c_opt);
00569                 <span class="keywordtype">bool</span> success = <span class="keyword">false</span>;
00570                 <span class="keywordflow">while</span> (!success)
00571                 {
00572                         p.<a class="code" href="class_integer.html#_integerz43_10">Randomize</a>(rng, I, I2, Integer::ANY);
00573                         p *= q; p &lt;&lt;= 1; ++p;
00574                         <span class="keywordflow">if</span> (!TrialDivision(p, trialDivisorBound))
00575                         {
00576                                 a.Randomize(rng, 2, p-1, Integer::ANY);
00577                                 b = a_exp_b_mod_c(a, (p-1)/q, p);
00578                                 success = (GCD(b-1, p) == 1) &amp;&amp; (a_exp_b_mod_c(b, q, p) == 1);
00579                         }
00580                 }
00581         }
00582         <span class="keywordflow">return</span> p;
00583 }
00584 
00585 Integer CRT(<span class="keyword">const</span> Integer &amp;xp, <span class="keyword">const</span> Integer &amp;p, <span class="keyword">const</span> Integer &amp;xq, <span class="keyword">const</span> Integer &amp;q, <span class="keyword">const</span> Integer &amp;u)
00586 {
00587         <span class="comment">// isn't operator overloading great?</span>
00588         <span class="keywordflow">return</span> p * (u * (xq-xp) % q) + xp;
00589 }
00590 
00591 Integer CRT(<span class="keyword">const</span> Integer &amp;xp, <span class="keyword">const</span> Integer &amp;p, <span class="keyword">const</span> Integer &amp;xq, <span class="keyword">const</span> Integer &amp;q)
00592 {
00593         <span class="keywordflow">return</span> CRT(xp, p, xq, q, EuclideanMultiplicativeInverse(p, q));
00594 }
00595 
00596 Integer ModularSquareRoot(<span class="keyword">const</span> Integer &amp;a, <span class="keyword">const</span> Integer &amp;p)