From 66a4490bf8f9b7b193133364bfd2d8add5d9f437 Mon Sep 17 00:00:00 2001 From: Kagre Date: Tue, 14 Jun 2016 15:57:11 -0600 Subject: [PATCH 1/8] iterate j over the odd numbers up to sqrt(i) - I don't have gcc, so I'm not too sure how well this'll compile - I used McDougall/Wotherspoon Newton-Raphson Method because I wasn't quite sure how to code the inline assembler equivalent ` _asm { fild i fsqrt }` - mostly revised for the odd j, and because i wanted to try and code `sqrt(i)` --- programs/prime.c | 37 +++++++++++++++++++++++++------------ 1 file changed, 25 insertions(+), 12 deletions(-) diff --git a/programs/prime.c b/programs/prime.c index 020bfc8..2f8adce 100644 --- a/programs/prime.c +++ b/programs/prime.c @@ -12,10 +12,10 @@ // 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 #include #include @@ -49,16 +49,29 @@ 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, iRoot=0xFFFFFFFF; // root finding variables - 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 McDougall/Wotherspoon to find square root of i + xx_k = iRoot; // x*_0 = x_0 + // max_number <= max_64 implies sqr(max_number) <= sqr(max_64) + f_val = iRoot * iRoot - i; // f(x_0) + df_val = iRoot * 2; // f'((x_0+x*_0)/2) = f'(x_0) + x_k = iRoot / 2 + i / 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 - i; + 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(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); From 40ad4a2e5b5c2adbbffb05e00fbaa397e969976b Mon Sep 17 00:00:00 2001 From: Kagre Date: Tue, 14 Jun 2016 16:40:15 -0600 Subject: [PATCH 2/8] Update primesmp.c --- programs/primesmp.c | 45 ++++++++++++++++++++++++++++++++++++--------- 1 file changed, 36 insertions(+), 9 deletions(-) diff --git a/programs/primesmp.c b/programs/primesmp.c index a0882bf..c232bb0 100644 --- a/programs/primesmp.c +++ b/programs/primesmp.c @@ -32,7 +32,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_number=0, max_root=0xFFFFFFFF, primes=1, local=0, process_stage=0, processes=0, max_processes=0, singletime=0, k=0; float speedup; time_t start, finish; @@ -81,6 +81,24 @@ 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 for all threads + unsigned long xx_k, x_k, f_val, df_val; // root finding variables + xx_k = max_root; // x*_0 = x_0 + // max_number <= max_64 implies sqr(max_number) <= sqr(max_64) + 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; k2; 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; + // using McDougall/Wotherspoon to find root of max_number for i + xx_k = iRoot; // x*_0 = x_0 + // max_number <= max_64 implies sqr(max_number) <= sqr(max_64) + f_val = iRoot * iRoot - i; // f(x_0) + df_val = iRoot * 2; // f'((x_0+x*_0)/2) = f'(x_0) + x_k = iRoot / 2 + i / 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 - i; + 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(j=3; j<=iRoot && i%j; j+=2); + if(j>iRoot)tprimes++; } // Continue loop up to max number // Add tprimes to primes. From 5dc294d3dd63b8c3b4683bf41a9dac16602311ec Mon Sep 17 00:00:00 2001 From: Kagre Date: Wed, 15 Jun 2016 11:17:24 -0600 Subject: [PATCH 3/8] Update primesmp.c --- programs/primesmp.c | 37 +++++++++++++++++++------------------ 1 file changed, 19 insertions(+), 18 deletions(-) diff --git a/programs/primesmp.c b/programs/primesmp.c index c232bb0..33a604f 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 @@ -32,7 +34,8 @@ 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, max_root=0xFFFFFFFF, primes=1, local=0, process_stage=0, processes=0, max_processes=0, singletime=0, k=0; +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; @@ -84,10 +87,9 @@ int main(int argc, char *argv[]) if(!(max_number&1)) max_number--; // drop max to odd - // using McDougall/Wotherspoon to find root of max_number for all threads + // 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 - // max_number <= max_64 implies sqr(max_number) <= sqr(max_64) 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) @@ -121,7 +123,7 @@ int main(int argc, char *argv[]) #else for (k=0; k2; i-=h) + for(i = max_number - (max_number - iFloor) % h; i>=iFloor; i-=h) { - // using McDougall/Wotherspoon to find root of max_number for i - xx_k = iRoot; // x*_0 = x_0 - // max_number <= max_64 implies sqr(max_number) <= sqr(max_64) - f_val = iRoot * iRoot - i; // f(x_0) - df_val = iRoot * 2; // f'((x_0+x*_0)/2) = f'(x_0) - x_k = iRoot / 2 + i / df_val; // x_1 = x_0 - f(x_0)/f'((x_0+x*_0)/2) = x_0 - f(x_0)/f'(x_0) + // find root of i + xx_k = iRoot; // x*_0 = x_0 + f_val = iRoot * iRoot - i; // f(x_0) + df_val = iRoot * 2; // f'((x_0+x*_0)/2) = f'(x_0) + x_k = iRoot / 2 + i / 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 - i; @@ -194,7 +195,7 @@ void *prime_process(void *param) // end root algo for(j=3; j<=iRoot && i%j; j+=2); if(j>iRoot)tprimes++; - } // Continue loop up to max number + } // Continue loop down from max number // Add tprimes to primes. #ifdef BAREMETAL From 4057002c4e8812fe0bd97b4669e0328fd82189f3 Mon Sep 17 00:00:00 2001 From: Kagre Date: Wed, 15 Jun 2016 11:39:54 -0600 Subject: [PATCH 4/8] clean up tab spacing --- programs/prime.c | 15 ++++++++------- 1 file changed, 8 insertions(+), 7 deletions(-) diff --git a/programs/prime.c b/programs/prime.c index 2f8adce..01c6733 100644 --- a/programs/prime.c +++ b/programs/prime.c @@ -16,6 +16,7 @@ // maxn = 1000000 =0x00000000000f4240 primes = 78498 // maxn = 5000000 =0x00000000004c4b40 primes = 348513 // maxn = 10000000 =0x0000000000989680 primes = 664579 +// maxn = 4294967295 =0x00000000ffffffff primes = 203280221 (max_32) #include #include @@ -51,16 +52,16 @@ int main(int argc, char *argv[]) time(&start); if(!(max_number&1))max_number--; // drop max to odd - unsigned long xx_k, x_k, f_val, df_val, iRoot=0xFFFFFFFF; // root finding variables + 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) for(i=max_number; i>2; i-=2) // reversed i to count down so previous step's iRoot is a good starting point { // using McDougall/Wotherspoon to find square root of i - xx_k = iRoot; // x*_0 = x_0 - // max_number <= max_64 implies sqr(max_number) <= sqr(max_64) - f_val = iRoot * iRoot - i; // f(x_0) - df_val = iRoot * 2; // f'((x_0+x*_0)/2) = f'(x_0) - x_k = iRoot / 2 + i / df_val; // x_1 = x_0 - f(x_0)/f'((x_0+x*_0)/2) = x_0 - f(x_0)/f'(x_0) + xx_k = iRoot; // x*_0 = x_0 + f_val = iRoot * iRoot - i; // f(x_0) + df_val = iRoot * 2; // f'((x_0+x*_0)/2) = f'(x_0) + x_k = iRoot / 2 + i / 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 - i; @@ -68,7 +69,7 @@ int main(int argc, char *argv[]) df_val = x_k + xx_k; x_k = (df_val * x_k - f_val) / df_val; } - // end root algo + // 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 From dcba6787374c4b7aa3df8475127cd9aa0d9b6f19 Mon Sep 17 00:00:00 2001 From: Kagre Date: Wed, 15 Jun 2016 16:12:46 -0600 Subject: [PATCH 5/8] update root finding split McDougall to run once/first, then Raphson thereafter --- programs/prime.c | 30 ++++++++++++++++-------------- 1 file changed, 16 insertions(+), 14 deletions(-) diff --git a/programs/prime.c b/programs/prime.c index 01c6733..4961eee 100644 --- a/programs/prime.c +++ b/programs/prime.c @@ -54,23 +54,25 @@ int main(int argc, char *argv[]) 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=max_number; i>2; i-=2) // reversed i to count down so previous step's iRoot is a good starting point { - // using McDougall/Wotherspoon to find square root of i - xx_k = iRoot; // x*_0 = x_0 - f_val = iRoot * iRoot - i; // f(x_0) - df_val = iRoot * 2; // f'((x_0+x*_0)/2) = f'(x_0) - x_k = iRoot / 2 + i / 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 - i; - 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 - + // using a step of Raphson to find the square root of i from previous iRoot + iRoot = (iRoot * iRoot + i) / (2 * iRoot); + // 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 From bd98800192be45027a674bf6a915f7185be8c885 Mon Sep 17 00:00:00 2001 From: Kagre Date: Wed, 15 Jun 2016 16:47:24 -0600 Subject: [PATCH 6/8] Update primesmp.c --- programs/primesmp.c | 19 +++++-------------- 1 file changed, 5 insertions(+), 14 deletions(-) diff --git a/programs/primesmp.c b/programs/primesmp.c index 33a604f..b8dcb38 100644 --- a/programs/primesmp.c +++ b/programs/primesmp.c @@ -157,7 +157,6 @@ int main(int argc, char *argv[]) void *prime_process(void *param) { register unsigned long h, i, j, tprimes=0, iRoot=max_root, iFloor; - unsigned long xx_k, x_k, f_val, df_val; // root finding variables // Lock process_stage, copy it to local var, subtract 1 from process_stage, unlock it. #ifdef BAREMETAL @@ -180,21 +179,13 @@ void *prime_process(void *param) // Process for(i = max_number - (max_number - iFloor) % h; i>=iFloor; i-=h) { - // find root of i - xx_k = iRoot; // x*_0 = x_0 - f_val = iRoot * iRoot - i; // f(x_0) - df_val = iRoot * 2; // f'((x_0+x*_0)/2) = f'(x_0) - x_k = iRoot / 2 + i / 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 - i; - 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; - } + // use a step of Raphson to find root of i + // I think one step is probably ok until around h >= 64379, then you may need two steps + iRoot = (iRoot * iRoot + i) / (2 * iRoot); // end root algo + for(j=3; j<=iRoot && i%j; j+=2); - if(j>iRoot)tprimes++; + if(j>iRoot) tprimes++; } // Continue loop down from max number // Add tprimes to primes. From 3b79c4e43ed2c6629e78b9b908f70a290debf47d Mon Sep 17 00:00:00 2001 From: Kagre Date: Wed, 15 Jun 2016 17:14:00 -0600 Subject: [PATCH 7/8] switch Raphson to Babylonian less opps, same diff. --- programs/primesmp.c | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/programs/primesmp.c b/programs/primesmp.c index b8dcb38..42a39c3 100644 --- a/programs/primesmp.c +++ b/programs/primesmp.c @@ -179,9 +179,9 @@ void *prime_process(void *param) // Process for(i = max_number - (max_number - iFloor) % h; i>=iFloor; i-=h) { - // use a step of Raphson to find root of i + // 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 * iRoot + i) / (2 * iRoot); + iRoot = (iRoot + i / iRoot) / 2; // end root algo for(j=3; j<=iRoot && i%j; j+=2); From 16eacc87edc67ed09223f67c480187f58f491845 Mon Sep 17 00:00:00 2001 From: Kagre Date: Wed, 15 Jun 2016 17:23:33 -0600 Subject: [PATCH 8/8] switch Raphson to Babylonian --- programs/prime.c | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/programs/prime.c b/programs/prime.c index 4961eee..f28d406 100644 --- a/programs/prime.c +++ b/programs/prime.c @@ -70,9 +70,10 @@ int main(int argc, char *argv[]) for(i=max_number; i>2; i-=2) // reversed i to count down so previous step's iRoot is a good starting point { - // using a step of Raphson to find the square root of i from previous iRoot - iRoot = (iRoot * iRoot + i) / (2 * iRoot); + // 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