diff --git a/programs/prime.c b/programs/prime.c index 020bfc8..f28d406 100644 --- a/programs/prime.c +++ b/programs/prime.c @@ -12,10 +12,11 @@ // Linux/BSD compile using GCC (Tested with 4.7) // gcc primesmp.c -o primesmp // -// maxn = 500000 primes = 41538 -// maxn = 1000000 primes = 78498 -// maxn = 5000000 primes = 348513 -// maxn = 10000000 primes = 664579 +// maxn = 500000 =0x000000000007a120 primes = 41538 +// maxn = 1000000 =0x00000000000f4240 primes = 78498 +// maxn = 5000000 =0x00000000004c4b40 primes = 348513 +// maxn = 10000000 =0x0000000000989680 primes = 664579 +// maxn = 4294967295 =0x00000000ffffffff primes = 203280221 (max_32) #include #include @@ -49,16 +50,32 @@ int main(int argc, char *argv[]) printf("Prime v1.5 - Searching up to %ld.\nProcessing...\n", max_number); time(&start); + + if(!(max_number&1))max_number--; // drop max to odd + unsigned long xx_k, x_k, f_val, df_val; // root finding variables + unsigned long iRoot=0xFFFFFFFF; // max_number <= max_64 implies sqr(max_number) <= sqr(max_64) + // using McDougall/Wotherspoon to find square root of max_number + xx_k = iRoot; // x*_0 = x_0 + f_val = iRoot * iRoot - max_number; // f(x_0) + df_val = iRoot * 2; // f'((x_0+x*_0)/2) = f'(x_0) + x_k = iRoot / 2 + max_number / df_val;// x_1 = x_0 - f(x_0)/f'((x_0+x*_0)/2) = x_0 - f(x_0)/f'(x_0) + for(j=1; j<30 && (iRoot - x_k); j++){ + iRoot = x_k; + f_val = x_k * x_k - max_number; + xx_k = (df_val * x_k - f_val) / df_val; + df_val = x_k + xx_k; + x_k = (df_val * x_k - f_val) / df_val; + } + // end root algo - for(i=3; i<=max_number; i+=2) + for(i=max_number; i>2; i-=2) // reversed i to count down so previous step's iRoot is a good starting point { - for(j=2; j*j<=i; j++) - { - if(i%j==0) break; //Number is divisble by some other number. So break out - } - if(j*j>i) - primes++; - + // using a step of Babylonian to find the square root of i from previous iRoot + iRoot = (iRoot + i / iRoot) / 2; + // end root algo + + for(j=3; j<=iRoot && i%j; j+=2); // test i for divisibility by j + if(j>iRoot)primes++; // count as prime when not divisible } //Continue loop up to max number time(&finish); diff --git a/programs/primesmp.c b/programs/primesmp.c index a0882bf..42a39c3 100644 --- a/programs/primesmp.c +++ b/programs/primesmp.c @@ -17,6 +17,8 @@ // maxn = 1000000 primes = 78498 // maxn = 5000000 primes = 348513 // maxn = 10000000 primes = 664579 +// maxn = 4294967295 primes = 203280221 +// maxn = 18446744073709551615 primes aprox 4.15829ee17 #include #include @@ -33,6 +35,7 @@ void *prime_process(void *param); // primes is set to 1 since we don't calculate for '2' as it is a known prime number unsigned long max_number=0, primes=1, local=0, process_stage=0, processes=0, max_processes=0, singletime=0, k=0; +unsigned long max_root=0xFFFFFFFF; // max_number <= max_64 implies sqr(max_number) <= sqr(max_64) float speedup; time_t start, finish; @@ -67,8 +70,8 @@ int main(int argc, char *argv[]) printf("Using a maximum of %ld process(es). Searching up to %ld.\n", max_processes, max_number); - for (processes=1; processes <= max_processes; processes*=2) - { + for (processes=1; processes <= max_processes; processes++) // changed to ++ from *=2 so processes matches comments below + { primes = 1; process_stage = processes; @@ -81,6 +84,23 @@ int main(int argc, char *argv[]) printf("Processing with %ld process(es)...\n", processes); time(&start); // Grab the starting time + + if(!(max_number&1)) max_number--; // drop max to odd + + // using McDougall/Wotherspoon to find root of max_number as starting point for all threads + unsigned long xx_k, x_k, f_val, df_val; // root finding variables + xx_k = max_root; // x*_0 = x_0 + f_val = max_root * max_root - max_number; // f(x_0) + df_val = max_root * 2; // f'((x_0+x*_0)/2) = f'(x_0) + x_k = max_root / 2 + max_number / df_val; // x_1 = x_0 - f(x_0)/f'((x_0+x*_0)/2) = x_0 - f(x_0)/f'(x_0) + for(int i=1; i<30 && (max_root - x_k); i++){ + max_root= x_k; + f_val = x_k * x_k - max_number; + xx_k = (df_val * x_k - f_val) / df_val; + df_val = x_k + xx_k; + x_k = (df_val * x_k - f_val) / df_val; + } + // end root algo // Spawn the worker processes for (k=0; k=iFloor; i-=h) { - for(j=2; j*j<=i; j++) - { - if(i%j==0) break; // Number is divisible by some other number. So break out - } - if(j*j>i) - { - tprimes = tprimes + 1; - } - } // Continue loop up to max number + // use a step of Babylonian to find root of i + // I think one step is probably ok until around h >= 64379, then you may need two steps + iRoot = (iRoot + i / iRoot) / 2; + // end root algo + + for(j=3; j<=iRoot && i%j; j+=2); + if(j>iRoot) tprimes++; + } // Continue loop down from max number // Add tprimes to primes. #ifdef BAREMETAL