Click here for the unformatted .cpp file




  1. /*
  2.     5 Prime counting functions implementations based on Linnik's Identity, including the
  3.     star here, an implementation that runs in O(n^2/3 ln n) time and O(n^1/3 ln n) space
  4.     (search for the function countprimes5, below).
  5.  
  6.     Nathan McKenzie, 8-9-2011
  7.     nathan _AT_ icecreambreakfast.com
  8.  
  9.     See http://www.icecreambreakfast.com/primecount/primecounting.html for descriptions
  10.     of the theory behind all of this stuff.
  11.  
  12.     See http://www.icecreambreakfast.com/primecount/primecounting.html#1_1 for a basic
  13.     description  of Linnik's identity,and
  14.     http://www.icecreambreakfast.com/primecount/primecounting.html#ch2  for a bit more
  15.     elaboration how this identity can be used to compute the prime counting function.
  16.  
  17.     As a general note, 4 of the 5 of these algorithms could be sped up considerably by
  18.     making use of a suitably sized wheel - the fourth algorithm (countprimes4) already
  19.     makes use of such a wheel.
  20.  
  21. */
  22.  
  23.  
  24. /* !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  25.  
  26. NOTE!  For easy exploration, change the following line
  27. */
  28.  
  29. #define primeCountingFunction countprimes5
  30.  
  31. /* to any of the following values : countprimes1, countprimes2, countprimes3,
  32.     countprimes4, countprimes5 to try out the different prime counting algorithms.
  33.  
  34.     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  35. */
  36.  
  37.  
  38.  
  39.  
  40.  
  41. #include "stdio.h"
  42. #include "stdlib.h"
  43. #include "math.h"
  44. #include "conio.h"
  45. #include "time.h"
  46.  
  47.  
  48. typedef long long BigInt;
  49. const int g_wheelLargestPrime = 19;
  50. int scaleNum = 10;
  51.  
  52. /* ________________________________________________________________________________________
  53.  
  54.         This is the main timing loop, running the prime counting function for increasingly
  55.         large input.
  56.  
  57.         Change the #define "primeCountingFunction" to the various "countprime" functions to
  58.         run and performance time them.
  59.    ________________________________________________________________________________________
  60. */
  61.  
  62. BigInt countprimes1( BigInt n );
  63. BigInt countprimes2( BigInt n );
  64. BigInt countprimes3( BigInt n );
  65. BigInt countprimes4( BigInt n );
  66. BigInt countprimes5( BigInt n );
  67. void MakeWheel( int );
  68.  
  69. int main(int argc, char* argv[]){
  70.     if( primeCountingFunction == countprimes4 ) MakeWheel( g_wheelLargestPrime );
  71.  
  72.     int oldClock = (int)clock();
  73.     int lastDif = 0;
  74.  
  75.     printf("                                                                    ");
  76.     printf( "Time\n");
  77.     printf("                                                                    ");
  78.     printf( "Increase\n");
  79.     printf("                                                                    ");
  80.     printf( "for x%d\n", scaleNum);
  81.     printf("         ");
  82.     printf( "__ Input Number __   __ Output Number __ _ MSec _ _ Sec _  Input\n");
  83.     printf( "                                                                    \n");
  84.     for( BigInt i = scaleNum; i <= 1000000000000000000; i *= scaleNum ){
  85.         printf( "%17I64d(10^%4.1f): ", i, log( (double)i )/log(10.0) );
  86.         BigInt total = (BigInt)(primeCountingFunction( i )+.00001);
  87.         int newClock = (int)clock();
  88.         printf( " %20I64d %8d : %4d: x%f\n",
  89.             total, newClock - oldClock, ( newClock - oldClock ) / CLK_TCK,
  90.             ( lastDif ) ? (double)( newClock - oldClock ) / (double)lastDif : 0.0 );
  91.         lastDif = newClock - oldClock;
  92.         oldClock = newClock;
  93.     }
  94.  
  95.     getch();
  96.  
  97.     return 0;
  98. }
  99.  
  100.  
  101.  
  102.  
  103.  
  104.  
  105.  
  106.  
  107.  
  108.  
  109.  
  110.  
  111. static BigInt mu[] = { 0, 1, -1, -1, 0, -1, 1, -1, 0, 0, 1, -1, 0, -1, 1, 1, 0, -1, 0,
  112.     -1, 0, 1, 1, -1, 0, 0, 1, 0, 0, -1, -1, -1, 0, 1, 1, 1, 0, -1, 1, 1, 0, -1, -1, -1,
  113.     0, 0, 1, -1, 0, 0, 0, 1, 0, -1, 0, 1, 0, 1, 1, -1, 0, -1, 1, 0, 0, 1, -1, -1, 0, 1,
  114.     -1, -1, 0, -1, 1, 0, 0, 1, -1, -1, 0, 0 };
  115.  
  116. /* ________________________________________________________________________________________
  117.  
  118.     ALGORITHM 1 : countprimes1
  119.  
  120.     This is the first implementation of the prime counting function.  It is not at all
  121.     fast, but it is at least rather simple.  It runs in something like O(n^1.7) time and
  122.     O(episolon) space - which makes it worse, in fact, that trial division.
  123.  
  124.     See http://www.icecreambreakfast.com/primecount/primecounting.html#3_1 for the source
  125.     identity of countdivisors.
  126.     See http://www.icecreambreakfast.com/primecount/primecounting.html#2_4 for the body of
  127.     countprimes1
  128.    ________________________________________________________________________________________
  129. */
  130.  
  131. double countdivisors(BigInt n, BigInt k){
  132.     if (k == 0) return 1;
  133.     if (k == 1) return n - 1;
  134.     BigInt t = 0;
  135.     for (BigInt i = 2; i <= n; i++) t += countdivisors( n / i, k - 1 );
  136.     return t;
  137. }
  138.  
  139. BigInt countprimes1( BigInt n){
  140.     double t = 0.0;
  141.     for (BigInt j = 1; j < log((double)n) / log(2.0); j++)
  142.         for (BigInt k = 1; k < log( pow( n, 1.0 / j ) ) / log(2.0); k++)
  143.             t += pow( -1.0, (double)k+1 ) * countdivisors(pow( n, 1.0 / j ), k)
  144.                 / k / j * mu[j];
  145.     return (t+.1);
  146. }
  147.  
  148. /* ________________________________________________________________________________________
  149.  
  150.     ALGORITHM 2 : countprimes2
  151.  
  152.     This is a recursive implementation of countprimes1.  It's ultimate quite similar,
  153.     including its lousy performance.
  154.  
  155.     See http://www.icecreambreakfast.com/primecount/primecounting.html#3_2 for the
  156.     source identity of recurseprimes.
  157.     See http://www.icecreambreakfast.com/primecount/primecounting.html#2_3 for the body
  158.     of countprimes2
  159.    ________________________________________________________________________________________
  160. */
  161.  
  162. double recurseprimes( BigInt n, BigInt k ){
  163.     double t= 0;
  164.     for( BigInt j = 2; j <= n; j++ )t += 1.0 / k - recurseprimes( n/j, k+1 );
  165.     return t;
  166. }
  167.  
  168. BigInt countprimes2( BigInt n ){
  169.     double t = 0;
  170.     for (BigInt j = 1; j < log((double)n) / log(2.0); j++)
  171.         t += recurseprimes(pow( n, 1.0 / j ), 1) / j * mu[j];
  172.     return t+.1;
  173. }
  174.  
  175. /* ________________________________________________________________________________________
  176.  
  177.     ALGORITHM 3 : countprimes3
  178.  
  179.     This is another, substantially faster, implementation of the prime counting function
  180.     via Linnik's Identity.  It relies on internal symmetries of the number of divisors
  181.     summatory function to be sped up.  For an algorithm that uses O(epsilon) memory, it
  182.     runs pretty quickly, appearing to clock in just a bit below the O(n) mark.
  183.    
  184.     This algorithm can be sped up massively, in constant time terms, by using a suitably
  185.     large wheel.
  186.  
  187.     See http://www.icecreambreakfast.com/primecount/primecounting.html#3_3 for the
  188.     source identity of countdivisorsfast
  189.     See http://www.icecreambreakfast.com/primecount/primecounting.html#2_4 for the
  190.     body of countprimes3
  191.    ________________________________________________________________________________________
  192. */
  193.  
  194. int inversepow( int n, int k) {
  195.     return pow(n, 1.0 / k) + .00000001;
  196. }
  197.  
  198. double factorial(BigInt val){
  199.     double total = 1.0;
  200.     for (int i = 1; i <= val; i++) total *= i;
  201.     return total;
  202. }
  203.  
  204. double binomial_(BigInt val, BigInt div) {
  205.     return factorial(val) / (factorial(div) * factorial(val - div));
  206. }
  207.  
  208. double countdivisorsfast( BigInt n, BigInt k, BigInt a ){
  209.     if( k == 0 )return 1;
  210.     if( k == 1 )return n - a + 1;
  211.     double t = 0;
  212.     for( BigInt m = a; m <= inversepow(n,k); m++ )
  213.         for( BigInt j = 0; j < k; j++ )
  214.             t += countdivisorsfast( n / pow( (double)m, (double)k - j ), j, m+1 )
  215.                 * binomial_( k, j );
  216.     return t;
  217. }
  218.  
  219. BigInt countprimes3(BigInt n){
  220.     double t = 0.0;
  221.     for (BigInt j = 1; j < log((double)n) / log(2.0); j++)
  222.         for (BigInt k = 1; k < log( pow( n, 1.0 / j ) ) / log(2.0); k++)
  223.             t += pow( -1.0, (double)k + 1 ) * countdivisorsfast(inversepow( n, j ), k, 2 )
  224.                 / k / j * mu[j];
  225.     return t + .1;
  226. }
  227.  
  228.  
  229.  
  230.  
  231.  
  232.  
  233.  
  234.  
  235.  
  236. /* ________________________________________________________________________________________
  237.  
  238.     ALGORITHM 4 : countprimes4
  239.  
  240.     Although it might be difficult to tell, the following code implements the same
  241.     algorithm, essentially, as ALGORITHM 3.  So, its core identities are found in
  242.     http://www.icecreambreakfast.com/primecount/primecounting.html#3_3
  243.     and http://www.icecreambreakfast.com/primecount/primecounting.html#2_4
  244.  
  245.     The main difference here is that a lot of integration has happened for performance
  246.     reasons (so, for example, partially computed binomials are kept around as computation
  247.     happens), and, especially, a large wheel has been implemented, which speeds up
  248.     computation drastically.
  249.  
  250.     This is a pretty fast prime counting algorithm considering it uses O(episilon) memory.  
  251.     I've tried, from a number of different angles, to come up with ways of trading off
  252.     memory usage to increase its runtime performance, but I haven't found any such approaches.
  253.  
  254.     I'm not entirely sure what the Big O runtime performance of this algorithm is,
  255.     ultimately... Up to around 10^15, it seems to be running in something like O(n^4/5),
  256.     but it also seems to be taking slightly longer each power of 10.
  257.   ________________________________________________________________________________________
  258. */
  259.  
  260.  
  261. const double EPSILON = .00000000001;
  262.  
  263. typedef long long BigInt;
  264.  
  265. /* A wheel including 19 has 9.6 million entries. */
  266.  
  267. /* Used for the construction of the wheel - include more primes as needed,
  268.    but this is already enough primes to consume over 75 gigs of RAM */
  269.  
  270. const int g_primes[] ={ 2,   3,   5,   7,  11,  13,  17,  19,  23,  29 };
  271.  
  272. int         g_wheelCycleEntries;
  273. int         g_wheelCyclePeriod;
  274. int         g_wheelFirstPrime;
  275. int         g_wheelBasePrimes;
  276. int*        g_wheelTranslation          = 0;
  277. int*        g_wheelOffsets              = 0;
  278.  
  279. BigInt      g_latticepoints;
  280. BigInt      g_minVarValue;
  281. BigInt      g_boundary;
  282.  
  283. BigInt      g_scale;
  284. BigInt      g_divisor;
  285.  
  286. BigInt      g_lastMax;
  287.  
  288. int         g_variablesLeft;
  289. BigInt      g_lastScaleDivisor;
  290. BigInt      g_scaleVal;
  291.  
  292. int g_mu[] = {
  293.   0,   1,  -1,  -1,   0,  -1,   1,  -1,   0,   0,
  294.   1,  -1,   0,  -1,   1,   1,   0,  -1,   0,  -1,
  295.   0,   1,   1,  -1,   0,   0,   1,   0,   0,  -1,
  296.  -1,  -1,   0,   1,   1,   1,   0,  -1,   1,   1,
  297.   0,  -1,  -1,  -1,   0,   0,   1,  -1,   0,   0,
  298.   0,   1,   0,  -1,   0,   1,   0,   1,   1,  -1,
  299.   0,  -1,   1,   0,   0,   1,  -1,  -1,   0,   1,
  300.  -1,  -1,   0,  -1,   1,   0,   0,   1,  -1,  -1,
  301.   0,   0,   1,  -1,   0,   1,   1,   1,   0,  -1,
  302.   0,   1,   0,   1,   1,   1,   0,  -1,   0,   0,
  303.   0,  -1,  -1,  -1,   0,  -1,   1,  -1,   0,  -1,
  304.  -1,   1,   0,  -1,  -1,   1,   0,   0,   1,   1,
  305.   0,   0,   1,   1,   0,   0,   0,  -1,   0,   1,
  306.  -1,  -1,   0,   1,   1,   0,   0,  -1,  -1,  -1,
  307.   0,   1,   1,   1,   0,   1,   1,   0,   0,  -1,
  308.   0,  -1,   0,   0,  -1,   1,   0,  -1,   1,   1,
  309.   0,   1,   0,  -1,   0,  -1,   1,  -1,   0,   0,
  310.  -1,   0,   0,  -1,  -1,   0,   0,   1,   1,  -1,
  311.   0,  -1,  -1,   1,   0,   1,  -1,   1,   0,   0,
  312.  -1,  -1,   0,  -1,   1,  -1,   0,  -1,   0,  -1,
  313.   0,   1,   1,   1,   0,   1,   1,   0,   0,   1,
  314.   1,  -1,   0,   1,   1,   1,   0,   1,   1,   1,
  315.   0,   1,  -1,  -1,   0,   0,   1,  -1,   0,  -1,
  316.  -1,  -1,   0,  -1,   0,   1,   0,   1,  -1,  -1,
  317.   0,  -1,   0,   0,   0,   0,  -1,   1,   0,   1,
  318.   0,  -1,   0,   1,   1,  -1,   0,
  319. };
  320.  
  321. /* Note that with 64 bit ints, we can't go above factorial( 20 ) anyway. */
  322. BigInt g_factorial[] = {
  323.     0, 1, 2, 6, 24, 120, 720, 5040, 40320, 362880, 3628800, 39916800, 479001600,
  324.     6227020800, 87178291200, 1307674368000, 20922789888000, 355687428096000,
  325.     6402373705728000, 121645100408832000, 2432902008176640000
  326. };
  327.  
  328. inline BigInt InversePower( BigInt x, BigInt y ){
  329.     return ( (BigInt)( pow( (double)x + EPSILON, ( 1.0 / (double)y ) ) + EPSILON ) );
  330. }
  331.  
  332. inline int GetwheelCyclePeriod( int cap ){
  333.     int val = 1;
  334.     int i   = 0;
  335.  
  336.     while( g_primes[ i ] <= cap ){
  337.         val *= g_primes[ i ];
  338.         i++;
  339.     }
  340.     return val;
  341. }
  342.  
  343. inline int GetFirstIncludedPrime( int cap ){
  344.     int i = 0;
  345.    
  346.     while( g_primes[ i ] <= cap ){
  347.         i++;
  348.     }
  349.     return g_primes[ i ];
  350. }
  351.  
  352. inline int GetBasePrimes( int cap ){
  353.     int i = 0;
  354.     while( g_primes[ i ] <= cap ){
  355.         i++;
  356.     }
  357.     return i;
  358. }
  359.  
  360. inline void IncrementWheel( int &offset ){
  361.     offset++;
  362.     if( offset >= g_wheelCycleEntries ){
  363.         offset = 0;
  364.     }
  365. }
  366.  
  367. void MakeWheel( int cap ){
  368.     g_wheelBasePrimes       = GetBasePrimes( cap );
  369.     g_wheelCyclePeriod      = GetwheelCyclePeriod( cap );
  370.     g_wheelCycleEntries     = 0;
  371.  
  372.     int cur                 = 0;
  373.     int offset              = -1;
  374.  
  375.     int* wheelBase          = 0;
  376.  
  377.     wheelBase               = ( int* )malloc( g_wheelCyclePeriod * sizeof( int ) );
  378.     g_wheelTranslation      = ( int* )malloc( ( g_wheelCyclePeriod + 1 ) * sizeof( int ) );
  379.     g_wheelOffsets          = ( int * )malloc( g_wheelCyclePeriod * sizeof( int ) );
  380.     g_wheelFirstPrime       = GetFirstIncludedPrime( cap );
  381.  
  382.     for( int i = 0; i < g_wheelCyclePeriod; i++ ){
  383.         wheelBase[ i ] = 1;
  384.         for( int j = 2; j <= cap; j++ ){
  385.             if( !( ( i + 1 ) % j ) ){
  386.                 wheelBase[ i ] = 0;
  387.                 break;
  388.             }
  389.         }
  390.     }
  391.  
  392.     while( cur < g_wheelCyclePeriod ){
  393.         if( wheelBase[ cur ] && cur != 0 ){
  394.             g_wheelOffsets[ g_wheelCycleEntries ] = offset + 1;
  395.             offset = 0;
  396.             g_wheelCycleEntries++;
  397.         }else{
  398.             offset++;
  399.         }
  400.         cur++;
  401.     }
  402.  
  403.     g_wheelOffsets[ g_wheelCycleEntries ] = 2;
  404.     g_wheelCycleEntries++; 
  405.  
  406.     int total = 0;
  407.  
  408.     g_wheelTranslation[ 0 ] = 0;
  409.     for( BigInt i = 0; i < g_wheelCyclePeriod; i++ ){
  410.         if( i && wheelBase[ i - 1 ] ){
  411.             total++;
  412.         }
  413.         g_wheelTranslation [ i + 1 ] = total;
  414.     }
  415.  
  416.     free( wheelBase );
  417. }
  418.  
  419. /* This function calculates how many entries the wheel leaves
  420.    in the range from (rangeStart, rangeStop].*/
  421. inline BigInt CountWheelEntries( BigInt rangeStart, BigInt rangeEnd ){
  422.     rangeEnd++;
  423.  
  424.     int a = rangeStart % g_wheelCyclePeriod;
  425.     int b = rangeEnd % g_wheelCyclePeriod;
  426.  
  427.     return ( rangeEnd - b - rangeStart + a ) / g_wheelCyclePeriod *
  428.         g_wheelCycleEntries + g_wheelTranslation[ b ] - g_wheelTranslation[ a ] ;
  429. }
  430.  
  431. void CountHyperbolaLattice_2_Variables( void ){
  432.     BigInt  finalBoundary   = g_boundary / g_minVarValue;
  433.     BigInt  boundaryRoot    = (BigInt)( sqrt( (double)finalBoundary ) + EPSILON );
  434.  
  435.     /* For if the final two digits happen to be the same. */
  436.     g_latticepoints += g_scale / ( g_divisor * ( g_divisor + 1 ) );
  437.  
  438.     /* Leading digit is same, final digit is not. */
  439.     g_latticepoints +=
  440.         ( CountWheelEntries( g_minVarValue, finalBoundary / g_minVarValue ) - 1 )
  441.         * ( g_scale / g_divisor );
  442.  
  443.     /* For if the final two digits happen to be the same,
  444.         but both differ from the previous. */
  445.     g_latticepoints +=
  446.         ( CountWheelEntries( g_minVarValue, boundaryRoot ) - 1 ) * ( g_scale / 2 );
  447.  
  448.     /* Both digits differ from all other digits - This is the hellish evil
  449.         loop of portentious doom. */
  450.  
  451.     int curWheelOffset      = g_wheelTranslation[ g_minVarValue % g_wheelCyclePeriod ];
  452.  
  453.     BigInt  curLeadingVar   = g_minVarValue + g_wheelOffsets[ curWheelOffset ];
  454.     BigInt  subTotal        = 0;
  455.  
  456.     IncrementWheel( curWheelOffset );
  457.  
  458.     while( curLeadingVar <= boundaryRoot ){
  459.         subTotal += CountWheelEntries( curLeadingVar, finalBoundary / curLeadingVar ) - 1;
  460.         curLeadingVar += g_wheelOffsets[ curWheelOffset ];
  461.         IncrementWheel( curWheelOffset );
  462.     }
  463.  
  464.     g_latticepoints += subTotal * g_scale;
  465. }
  466.  
  467. void CountHyperbolaLattice_3_Variables( BigInt hyperbolaBoundary, BigInt minVarValue ){
  468.     BigInt maxVarValue  = InversePower( hyperbolaBoundary, 3 );
  469.  
  470.     int curWheelOffset  = g_wheelTranslation[ minVarValue % g_wheelCyclePeriod ];
  471.  
  472.     g_boundary          = hyperbolaBoundary;
  473.     g_minVarValue       = minVarValue;
  474.     g_scale             = g_scaleVal / g_lastScaleDivisor;
  475.     g_divisor           = g_lastScaleDivisor + 1;
  476.  
  477.     CountHyperbolaLattice_2_Variables();
  478.  
  479.     g_minVarValue += g_wheelOffsets[ curWheelOffset ];
  480.     IncrementWheel( curWheelOffset );
  481.  
  482.     g_scale     = g_scaleVal;
  483.     g_divisor   = 2;
  484.     while( g_minVarValue <= maxVarValue ){
  485.         CountHyperbolaLattice_2_Variables();
  486.         g_minVarValue += g_wheelOffsets[ curWheelOffset ];
  487.         IncrementWheel( curWheelOffset );
  488.     }
  489. }
  490.  
  491. void CountHyperbolaLattice_X_Variables( BigInt hyperbolaBoundary, BigInt minVarValue ){
  492.     BigInt         maxVarValue      = InversePower( hyperbolaBoundary, g_variablesLeft );
  493.  
  494.     /* Save global variables that will be restored at end of function */
  495.  
  496.     BigInt      scaleVal                = g_scaleVal;
  497.     BigInt      lastScaleDivisor        = g_lastScaleDivisor;
  498.  
  499.     int curWheelOffset = g_wheelTranslation[ minVarValue % g_wheelCyclePeriod ];
  500.  
  501.     g_variablesLeft--;
  502.     g_lastScaleDivisor                  = lastScaleDivisor + 1;
  503.     g_scaleVal                          = scaleVal / lastScaleDivisor;
  504.  
  505.     if( g_variablesLeft == 3 ){
  506.         CountHyperbolaLattice_3_Variables( hyperbolaBoundary / minVarValue, minVarValue );
  507.     }else{
  508.         CountHyperbolaLattice_X_Variables( hyperbolaBoundary / minVarValue, minVarValue );
  509.     }
  510.  
  511.     g_lastScaleDivisor                  = 2;
  512.     g_scaleVal                          = scaleVal;
  513.     minVarValue += g_wheelOffsets[ curWheelOffset ];
  514.     IncrementWheel( curWheelOffset );
  515.  
  516.     while( minVarValue <= maxVarValue ){
  517.         if( g_variablesLeft == 3 ){
  518.             CountHyperbolaLattice_3_Variables( hyperbolaBoundary / minVarValue, minVarValue );
  519.         }else{
  520.             CountHyperbolaLattice_X_Variables( hyperbolaBoundary / minVarValue, minVarValue );
  521.         }
  522.         minVarValue += g_wheelOffsets[ curWheelOffset ];
  523.         IncrementWheel( curWheelOffset );
  524.     }
  525.  
  526.     /* Restore global variables */
  527.  
  528.     g_lastScaleDivisor = lastScaleDivisor;
  529.     g_variablesLeft++;
  530. }
  531.  
  532. BigInt CountHyperbolaLattice( BigInt hyperbolaBoundary, int hyperbolaVariables ){
  533.     g_latticepoints     = 0;
  534.     g_variablesLeft     = hyperbolaVariables;
  535.     g_lastScaleDivisor  = 1;
  536.     g_scaleVal          = g_factorial[ hyperbolaVariables ];
  537.  
  538.     if( hyperbolaBoundary <
  539.             (BigInt)pow( (double)g_wheelFirstPrime, (double)hyperbolaVariables ) ){
  540.         return 0;
  541.     }
  542.  
  543.     switch( hyperbolaVariables ){
  544.     case 1:
  545.         g_latticepoints = CountWheelEntries( g_wheelFirstPrime, hyperbolaBoundary );
  546.         break;
  547.     case 2:
  548.         /* CountHyperbolaLattice_2_Variables expects a number of global variables
  549.            to be initialized when it is called, which generally happens in
  550.            CountHyperbolaLattice_3_Variables.  We have to do it manually here. */
  551.         g_minVarValue   = g_wheelFirstPrime;
  552.         g_boundary      = g_wheelFirstPrime * hyperbolaBoundary;
  553.         g_scale         = 2;
  554.         g_divisor       = 1;
  555.  
  556.         CountHyperbolaLattice_2_Variables();
  557.         break;
  558.     case 3:
  559.         CountHyperbolaLattice_3_Variables( hyperbolaBoundary, g_wheelFirstPrime );
  560.         break;
  561.     default:
  562.         CountHyperbolaLattice_X_Variables( hyperbolaBoundary, g_wheelFirstPrime );
  563.         break;
  564.     }
  565.     return g_latticepoints;
  566. }
  567.  
  568. BigInt countprimes4( BigInt n ){
  569.     int maxPower = ( int )( log( ( double )n + EPSILON )
  570.         / log ( ( double )g_wheelFirstPrime + EPSILON ) + EPSILON ) + 1;
  571.     double  total           = 0.0;
  572.  
  573.     int     oldClock        = (int)clock();
  574.     int     totalTime       = 0;
  575.    
  576.     for( int curPower = 1; curPower < maxPower; curPower++ ){
  577.         if( !g_mu[ curPower ] ){
  578.             continue;
  579.         }
  580.  
  581.         BigInt  curMax              = InversePower( n, curPower );
  582.         double  subTotal            = 0.0;
  583.         BigInt  hyperbolaEntries    = 1;
  584.         double  sign                = 1;
  585.  
  586.         while( 1 ){
  587.             double temp     = sign / hyperbolaEntries;
  588.             sign            *= -1;
  589.  
  590.             double v2       = (double)CountHyperbolaLattice( curMax, hyperbolaEntries );
  591.             temp            *= v2;
  592.  
  593.             if( temp == 0.0 ){
  594.                 break;
  595.             }
  596.  
  597.             subTotal        += temp;
  598.             hyperbolaEntries++;
  599.  
  600.             int newClock    = (int)clock();
  601.             totalTime       += newClock - oldClock;
  602.             oldClock        = newClock;
  603.  
  604.         }
  605.         subTotal    /= curPower * g_mu[ curPower ];
  606.         total       += subTotal;
  607.     }
  608.     total += g_wheelBasePrimes;
  609.  
  610.     /* the .5 is to prevent truncation errors - but it's clearly sloppy */
  611.     return (BigInt)( total + 0.5 );
  612. }
  613.  
  614.  
  615.  
  616.  
  617.  
  618.  
  619.  
  620.  
  621.  
  622.  
  623.  
  624.  
  625.  
  626. /* ________________________________________________________________________________________
  627.  
  628.     ALGORITHM 5
  629.  
  630.     This is the fastest, in terms of O() runtime, and by far most complex, implementation
  631.     of a prime counting algorithm based on Linnik's identity that I have found.
  632.  
  633.     Its performance is none too shabby - it runs in roughly O( n^2/3 ln n ) time and
  634.     O( n^1/3 ln n ) space, which puts it in spitting distance of some of the faster methods
  635.     for counting primes.
  636.  
  637.     See http://www.icecreambreakfast.com/primecount/primecounting.html#ch4 ,
  638.         http://www.icecreambreakfast.com/primecount/primecounting.html#ch5 , and
  639.         http://www.icecreambreakfast.com/primecount/primecounting.html#ch6
  640.     for general descriptions of the algorithm implemented here.
  641.  
  642.     As currently implemented here, the code stops working for relatively low input values
  643.     (say, 10^11 ish) - this is a consequence of precision and overflow issues, not anything
  644.     to do with the underlying math.  Obviously those concerns would have to be addressed to
  645.     use it for higher numbers.
  646.  
  647.     As repeatedly mentioned above, this algorithm should be sped up considerably if a wheel
  648.     is implemeneted, as was done in ALGORITHM 4 / countprimes4.  To date, I haven't written
  649.     such an implementation.  A suitable wheel will come with the side bit of making precision
  650.     issues somewhat more tractable.
  651.  
  652.     I haven't gone out of my way to do much in the way of constant time performance
  653.     optimization here - it's already unreadable enough without fiddling with the code
  654.     thusly, and integrating a wheel ought to come first anyway.  But there are plenty
  655.     of opportunities.
  656.  
  657.     I'm really curious how fast this approach can be sped up.
  658.  
  659.     The fuzzy transitions between double and BigInt are a cause of imprecision and
  660.     consternation.
  661.    ________________________________________________________________________________________
  662. */
  663.  
  664. static BigInt* binomials;       /* This is used as a doubly subscripted array, 128x128.  
  665.                                 Indexing is done manually.*/
  666. static BigInt nToTheThird;
  667. static BigInt logn;
  668.  
  669. static BigInt numPrimes;
  670. static BigInt* primes;
  671.  
  672. static BigInt* factorsMultiplied;
  673. static BigInt* totalFactors;
  674. static BigInt* factors;         /* This is used as a doubly subscripted array, n^1/3 x ln n.
  675.                                 Indexing is done manually.*/
  676. static BigInt* numPrimeBases;
  677.  
  678. static BigInt* DPrime;          /*This is used as a doubly subscripted array, n^1/3 x ln n.
  679.                                 Indexing is done manually.*/
  680.  
  681. static BigInt curBlockBase;
  682.  
  683. static double t;
  684.  
  685. static BigInt nToTheHalf;
  686. static BigInt numDPowers;
  687. static double* dPrime;
  688.  
  689. static BigInt S1Val;
  690. static BigInt S1Mode;
  691. static BigInt* S3Vals;
  692. static BigInt* S3Modes;
  693.  
  694. static bool ended;
  695. static BigInt maxSieveValue;
  696.  
  697. static BigInt ceilval;
  698.  
  699. static BigInt n;
  700.  
  701. BigInt binomial( double n, int k ){
  702.     double t = 1;
  703.     for( int i = 1; i <= k; i++ ){
  704.         t *= ( n - ( k - i ) ) / i;
  705.     }
  706.     return BigInt( t + .1 );
  707. }
  708.  
  709. static BigInt invpow(double n, double k) {
  710.     return (BigInt)(pow(n, 1.0 / k) + .00000001);
  711. }
  712.  
  713. /* See http://www.icecreambreakfast.com/primecount/primecounting.html#ch5 for a description
  714.     of calculating d_k'(n) from a complete factorization of a number n. */
  715. static BigInt d1(BigInt* a, BigInt o, BigInt k, BigInt l){
  716.     BigInt t = 1;
  717.     for (BigInt j = 0; j < l; j++) t *= binomials[(a[o*logn+ j] - 1 + k)*128 + a[o*logn+ j]];
  718.     return t;
  719. }
  720.  
  721. /* See http://www.icecreambreakfast.com/primecount/primecounting.html#ch5 for a description
  722.     of calculating d_k'(n) from a complete factorization of a number n.*/
  723. static BigInt d2(BigInt* a, BigInt o, BigInt k, BigInt l, BigInt numfacts ){
  724.     if (numfacts < k) return 0;
  725.     BigInt t = 0;
  726.     for (BigInt j = 1; j <= k; j++) t += ( ( k - j ) % 2 == 1 ? -1:1 )
  727.         * binomials[k * 128 + j] * d1(a, o, j, l);
  728.     return (BigInt)t;
  729. }
  730.  
  731. static void allocPools( BigInt n ){
  732.     nToTheThird = (BigInt)pow(n, 1.0 / 3);
  733.  
  734.     logn = (BigInt)(log(pow(n, 2.00001 / 3)) / log(2.0)) + 1;
  735.     factorsMultiplied = new BigInt[nToTheThird];
  736.     totalFactors = new BigInt[nToTheThird];
  737.     factors = new BigInt[nToTheThird * logn];
  738.     numPrimeBases = new BigInt[nToTheThird];
  739.     DPrime = new BigInt[(nToTheThird + 1) * logn];
  740.     binomials = new BigInt[128*128+ 128];
  741.     for (BigInt j = 0; j < 128; j++)
  742.         for (BigInt k = 0; k <= j; k++)
  743.             binomials[j * 128 + k]= binomial(j, k);
  744.     for (BigInt j = 0; j < logn; j++) DPrime[j] = 0;
  745.     curBlockBase = 0;
  746.  
  747.     t = n - 1;
  748.  
  749.     nToTheHalf = (BigInt)pow(n, 1.0 / 2);
  750.     numDPowers = (BigInt)(log(pow(n, 2.00001 / 3)) / log(2.0)) + 1;
  751.     dPrime = new double[(nToTheThird + 1) * (numDPowers + 1)];
  752.  
  753.     S1Val = 1;
  754.     S1Mode = 0;
  755.     S3Vals = new BigInt[nToTheThird + 1];
  756.     S3Modes = new BigInt[nToTheThird + 1];
  757.  
  758.     ended = false;
  759.     maxSieveValue = (BigInt)(pow(n, 2.00001 / 3));
  760.  
  761.     for (BigInt j = 2; j < nToTheThird + 1; j++){
  762.         S3Modes[j] = 0;
  763.         S3Vals[j] = 1;
  764.     }
  765. }
  766.  
  767. static void deallocPools(){
  768.     delete factorsMultiplied;
  769.     delete totalFactors;
  770.     delete factors;
  771.     delete numPrimeBases;
  772.     delete DPrime;
  773.     delete binomials;
  774.     delete dPrime;
  775.     delete S3Vals;
  776.     delete S3Modes;
  777.     delete primes;
  778. }
  779.  
  780. /* This finds all the primes less than n^1/3, which will be used for sieving and
  781.     generating complete factorizations of numbers up to n^2/3*/
  782. static void fillPrimes(){
  783.     BigInt* primesieve = new BigInt[nToTheThird + 1];
  784.     primes = new BigInt[nToTheThird + 1];
  785.     numPrimes = 0;
  786.     for (BigInt j = 0; j <= nToTheThird; j++) primesieve[j] = 1;
  787.     for (BigInt k = 2; k <= nToTheThird; k++){
  788.         BigInt cur = k;
  789.         if (primesieve[k] == 1){
  790.             primes[numPrimes] = k;
  791.             numPrimes++;
  792.             while (cur <= nToTheThird){
  793.                 primesieve[cur] = 0;
  794.                 cur += k;
  795.             }
  796.         }
  797.     }
  798.     delete primesieve;
  799. }
  800.  
  801. /*This resets some state used for the sieving and factoring process.*/
  802. static void clearPools(){
  803.     for (BigInt j = 0; j < nToTheThird; j++){
  804.         numPrimeBases[j] = -1;
  805.         factorsMultiplied[j] = 1;
  806.         totalFactors[j] = 0;
  807.     }
  808. }
  809.  
  810. /* We can use sieving on our current n^1/3 sized block of numbers to
  811.  get their complete prime factorization signatures, with which we can then
  812.  quickly compute d_k' values.*/
  813. static void factorRange(){
  814.     for (BigInt j = 0; j < numPrimes; j++){
  815.         /*mark everything divided by each prime, adding a new entry.*/
  816.         BigInt curPrime = primes[j];
  817.         if (curPrime * curPrime > curBlockBase + nToTheThird) break;
  818.         BigInt curEntry = ( curBlockBase % curPrime == 0 ) ? 0:curPrime
  819.             - (curBlockBase % curPrime);
  820.         while (curEntry < nToTheThird){
  821.             if( curEntry+curBlockBase != 0 ){
  822.                 factorsMultiplied[curEntry] *= curPrime;
  823.                 totalFactors[curEntry]++;
  824.                 numPrimeBases[curEntry]++;
  825.                 factors[curEntry*logn+ numPrimeBases[curEntry]] = 1;
  826.             }
  827.             curEntry += curPrime;
  828.         }
  829.         /*mark everything divided by each prime power*/
  830.         BigInt cap = (BigInt)( log((double)(nToTheThird+curBlockBase))
  831.             / log((double)curPrime) + 1 );
  832.         BigInt curbase = curPrime;
  833.         for (BigInt k = 2; k < cap; k++){
  834.             curPrime *= curbase;
  835.             curEntry = (curBlockBase % curPrime == 0) ? 0 : curPrime
  836.                 - (curBlockBase % curPrime);
  837.             while (curEntry < nToTheThird){
  838.                 factorsMultiplied[curEntry] *= curbase;
  839.                 totalFactors[curEntry]++;
  840.                 if (curEntry + curBlockBase != 0)
  841.                     factors[curEntry*logn+ numPrimeBases[curEntry]] = k;
  842.                 curEntry += curPrime;
  843.             }
  844.         }
  845.     }
  846.     /*account for prime factors > n^1/3*/
  847.     for (BigInt j = 0; j < nToTheThird; j++){
  848.         if (factorsMultiplied[j] < j+curBlockBase){
  849.             numPrimeBases[j]++;
  850.             totalFactors[j]++;
  851.             factors[j*logn+ numPrimeBases[j]] = 1;
  852.         }
  853.     }
  854. }
  855.  
  856. /* By this point, we have already factored, through sieving, all the numbers in the
  857.     current n^1/3 sized block we are looking at.  With a complete factorization, we
  858.     can calculate d_k'(n) for a number.  Then, D_k'(n) = d_k'(n) + D_k'(n-1).*/
  859. static void buildDivisorSums(){
  860.     for (BigInt j = 1; j < nToTheThird+1; j++){
  861.         if (j + curBlockBase == 1 || j + curBlockBase == 2) continue;
  862.         for (BigInt k = 0; k < logn; k++){
  863.             DPrime[j * logn + k] = DPrime[(j - 1) * logn + k] +
  864.                 d2(factors, j - 1, k, numPrimeBases[j - 1] + 1, totalFactors[j - 1]);
  865.         }
  866.     }
  867.     for (BigInt j = 0; j < logn; j++) DPrime[j] = DPrime[nToTheThird*logn+ j];
  868. }
  869.  
  870. /*  This general algorithm relies on values of D_k' <= n^2/3 and d_k' <= n^1/3.  
  871.     This function calculates those values of d_k'.*/
  872. static void find_dVals(){
  873.     curBlockBase = 1;
  874.     clearPools();
  875.     factorRange();
  876.     buildDivisorSums();
  877.  
  878.     for (BigInt j = 2; j <= nToTheThird; j++){
  879.         for (BigInt m = 1; m < numDPowers; m++){
  880.             double s = 0;
  881.             for (BigInt r = 1; r < numDPowers; r++)
  882.                 s += pow(-1.0, (double)( r + m )) * (1.0 / (r + m + 1)) *
  883.                     (DPrime[j * logn + r] - DPrime[(j - 1) * logn + r]);
  884.             dPrime[j*(numDPowers + 1)+ m] = s;
  885.         }
  886.     }
  887. }
  888.  
  889. static void resetDPrimeVals(){
  890.     curBlockBase = 0;
  891.     for (BigInt k = 0; k < nToTheThird + 1; k++)
  892.         for (BigInt j = 0; j < logn; j++)
  893.             DPrime[k * logn + j] = 0;
  894. }
  895.  
  896. /* This function is calculating the first two sums of
  897.     http://www.icecreambreakfast.com/primecount/primecounting.html#4_4  
  898.     It is written to rely on values of D_k' from smallest to greatest,
  899.     to use the segmented sieve.*/
  900. static void calcS1(){
  901.     if (S1Mode == 0){
  902.         while (S1Val <= ceilval){
  903.             BigInt cnt = (n / S1Val - n / (S1Val + 1));
  904.             for (BigInt m = 1; m < numDPowers; m++)
  905.                 t += cnt * (m % 2 == 1 ? -1 : 1) * (1.0 / (m + 1)) *
  906.                     DPrime[(S1Val - curBlockBase + 1) * logn + m];
  907.             S1Val++;
  908.             if (S1Val >= n / nToTheHalf){
  909.                 S1Mode = 1;
  910.                 S1Val = nToTheHalf;
  911.                 break;
  912.             }
  913.         }
  914.     }
  915.     if (S1Mode == 1){
  916.         while (n / S1Val <= ceilval){
  917.             for (BigInt m = 1; m < numDPowers; m++)
  918.                 t += (m % 2 == 1 ? -1 : 1) * (1.0 / (m + 1)) *
  919.                     DPrime[(n / S1Val - curBlockBase + 1) * logn + m];
  920.             S1Val--;
  921.             if (S1Val < nToTheThird + 1){
  922.                 S1Mode = 2;
  923.                 break;
  924.             }
  925.         }
  926.     }
  927. }
  928.  
  929. /* This loop is calculating the 3rd term that runs from 2 to n^1/3 in
  930.     http://www.icecreambreakfast.com/primecount/primecounting.html#4_4 */
  931. static void calcS2(){
  932.     for (BigInt j = 2; j <= nToTheThird; j++)
  933.         for (BigInt k = 1; k < numDPowers; k++)
  934.             t += (n / j - 1) * pow(-1.0, (double)k) * (1.0 / (k + 1)) *
  935.                 (DPrime[j * logn + k] - DPrime[(j - 1) * logn + k]);
  936. }
  937.  
  938. /* This loop is calculating the two double sums in
  939.     http://www.icecreambreakfast.com/primecount/primecounting.html#4_4
  940.     It is written to rely on values of D_k' from smallest to greatest,
  941.     to use the segmented sieve.*/
  942. static void calcS3(){
  943.     for (BigInt j = 2; j <= nToTheThird; j++){
  944.         if (S3Modes[j] == 0){
  945.             BigInt endsq = (BigInt)(pow(n / j, .5));
  946.             BigInt endVal = (n / j) / endsq;
  947.             while (S3Vals[j] <= ceilval){
  948.                 BigInt cnt = (n / (j * S3Vals[j]) - n / (j * (S3Vals[j] + 1)));
  949.                 for (BigInt m = 1; m < numDPowers; m++)
  950.                     t += cnt * DPrime[(S3Vals[j] - curBlockBase + 1) * logn + m] *
  951.                         dPrime[j*(numDPowers + 1)+ m];
  952.                 S3Vals[j]++;
  953.                 if (S3Vals[j] >= endVal){
  954.                     S3Modes[j] = 1;
  955.                     S3Vals[j] = endsq;
  956.                     break;
  957.                 }
  958.             }
  959.         }
  960.         if (S3Modes[j] == 1){
  961.             while (n / (j * S3Vals[j]) <= ceilval){
  962.                 for (BigInt m = 1; m < numDPowers; m++)
  963.                     t += DPrime[(n / (j * S3Vals[j]) - curBlockBase + 1) * logn + m] *
  964.                         dPrime[j * (numDPowers + 1) + m];
  965.                 S3Vals[j]--;
  966.                 if (S3Vals[j] < nToTheThird / j + 1){
  967.                     S3Modes[j] = 2;
  968.                     break;
  969.                 }
  970.             }
  971.         }
  972.     }
  973. }
  974.  
  975. /*  This is the most important function here. How it works:
  976.      first we allocate our n^1/3 ln n sized pools and other variables.
  977.      Then we go ahead and sieve to get our primes up to n^1/3
  978.      We also calculate, through one pass of sieving, values of d_k'(n) up to n^1/3
  979.      Then we go ahead and calculate the loop S2 (check the description of the algorithm),
  980.         which only requires values of d_k'(n) up to n^1/3, which we already have.
  981.      Now we're ready for the main loop.
  982.      We do the following roughly n^1/3 times.
  983.      First we clear our sieving variables.
  984.      Then we factor, entirely, all of the numbers in the current block sized n^1/3 that
  985.         we're looking at.
  986.      Using our factorization information, we calculate the values for d_k'(n) for the
  987.         entire range we're looking, and then sum those together to have a rolling set
  988.         of D_k'(n) values
  989.      Now we have values for D_k'(n) for this block sized n^1/3
  990.      First we see if any of the values of S1 that we need to compute are in this block. We
  991.         can do this by (see the paper) walking through the two S1 loops backwards, which
  992.         will use the D_k'(n) values in order from smallest to greatest
  993.      We then do the same thing will all of the S3 values
  994.      Once we have completed this loop, we will have calculated the prime power function for n.
  995.    
  996.     This loop is essentially calculating
  997.     http://www.icecreambreakfast.com/primecount/primecounting.html#4_4
  998. */
  999.  
  1000. static double calcPrimePowerCount(BigInt nVal){
  1001.     n = nVal;
  1002.     allocPools(n);
  1003.     fillPrimes();
  1004.     find_dVals();
  1005.     calcS2();
  1006.     resetDPrimeVals();
  1007.  
  1008.     for (curBlockBase = 0; curBlockBase <= maxSieveValue; curBlockBase += nToTheThird ){
  1009.         clearPools();
  1010.         factorRange();
  1011.         buildDivisorSums();
  1012.  
  1013.         ceilval = curBlockBase + nToTheThird - 1;
  1014.         if (ceilval > maxSieveValue) {
  1015.             ceilval = maxSieveValue;
  1016.             ended = true;
  1017.         }
  1018.  
  1019.         calcS1();
  1020.         calcS3();
  1021.         if (ended) break;
  1022.     }
  1023.  
  1024.     deallocPools();
  1025.  
  1026.     return t;
  1027. }
  1028.  
  1029. /*Make sure to read the header comments labeled "ALGORITHM 5", above*/
  1030. static BigInt countprimes5(BigInt num) {
  1031.     double total = 0.0;
  1032.     for (BigInt i = 1; i < log((double)num) / log(2.0); i++) {
  1033.         double val = calcPrimePowerCount( invpow(num, i)) / (double)i * mu[i];
  1034.         total += val;
  1035.     }
  1036.     return total+.1;
  1037. }


Share |





(c) 2011 icecreambreakfast.com | Contact