Let’s get right to the 3rd problem on ProjectEuler.net.
Question 3:
The prime factors of 13195 are 5, 7, 13 and 29.
What is the largest prime factor of the number 600851475143 ?
That’s a good question, and one that had me scratching my head for a second or two. I will be the first to admit that my knowledge of math revolving around primes is lacking. As such, my very first Prime algorithm was very bad. I attempted to find all of the primes between 1 and limit/2 (in this case: 300425737571). I then attempted to perform a modulus on each prime (iteratively) to find if it was a factor. Once the loop completed, I would know my largest prime factor. Sound in theory, terrible in practice. 20 minutes into the program running, I was still calculating all of the prime numbers (granted, that calculation could have used some optimization as well). 10 minutes after that, I was finally testing my numbers against my limit, but at that point it was a forgone conclusion that this was not going to work. I am not even going to list that code here as it would just be a waste of time. This leads me to my revised code.
Psuedo-Prime: After doing some research and trying a few things out, I came up with the following:
double Run()
{
limit = 600851475143;
double factor = 2;
while(factor < limit)
{
if(fmod(limit, factor) == 0)
limit /= factor;
else
factor += 1;
}
return limit;
}
Results:
Limit Speed 600851475143 216 microseconds 9512811111 137 milliseconds 600851475111 95 minutes
I realize now how simple this really is. Assuming you start at 2 and divide it out as many times as possible (very important) before moving on, by the time your factor is bigger than your limit, it will be prime (you have removed all of the smaller non-prime factors already). This algorithm finds the solution to the given question in 216 microseconds. It wasn’t until I tried to gauge the efficiency for a larger number that I realized just how tricky that was going to be. You see, a bigger number does not always mean a more difficult number. For example, given this algorithm, the number 600851475143 is found quickly, but the number 9512811111 (a number ~63 times smaller than the first) takes 137 milliseconds. This is because the original number has a relatively small (4 digits) prime factor, however the prime factor of the second number is 4510579. Obviously it takes many more iterations to solve the second number than the first. Even a number like 1*10^300 would finish in 0 milliseconds because the largest prime factor is 5 (that goes for any number in the 10,100,1000,10000…. series). Imagine if I tried the same algorithm on the number 600851475111. It is still a smaller number than our first, but its largest prime factor is 200283825037. Completing that number took 95 minutes!
Optimization 1: When thinking about optimizing this problem, the first thing I considered was recursion. If I could quickly get all of the primes between 1 and limit/2, I could start at the end and work my way down. In fact, if I did that, I could solve the question for my third number (600851475111) in no time. Would this be an improvement in efficiency? The answer is yes and no. Flipping the order doesn’t do much to improve the algorithm, it just flips the order. With recursion, I would find answers to problems that had huge prime factors very quickly, but finding the prime factor of 1*10^300 would take a vast amount of time since the largest prime factor is in fact very small. This excludes the time “tax” that would be involved in getting the primes to begin with (some languages have fast built in methods, some don’t).
With this in mind, I decide to continue with iteration and try to find a way to cut the number of loops to improve performance. The first thing I notice is that we are starting with the prime 2 and incrementing by 1 from there. Since 2 is the only even prime, we can do 2 first and then start with 3 and increment by 2 thereafter. This should theoretically cut the number of iterations by half (a good number in my book). The code:
double Run()
{
limit = 600851475143;
double factor = 2;
while(factor == 2)
{
if(fmod(limit, factor) == 0)
limit /= factor;
else
factor += 1;
}
while(factor < limit)
{
if(fmod(limit, factor) == 0)
limit /= factor;
else
factor += 2;
}
return limit;
}
Results:
Limit Speed 600851475143 107 microseconds 9512811111 69 milliseconds 600851475111 48 minutes
While it is more code, it can be seen to have increased by speed by a ratio of about 1:2 (which is exactly what we guessed).
Optimization 2: The Project Euler website has an optimization mentioned that I touched upon at the beginning of this process. Originally, I was finding all of the prime numbers between 1 and limit / 2 (where limit is the given number). This was in the right track, but not exactly correct. The actual rule is that no number ‘n’ can have more than 1 prime factor higher than Sqrt(n). This in essence reduces the number of iterations required by a square. I believe my code also has some minor improvements over the code proposed by Project Euler in that I do not store the “last factor” and instead directly modify the variable ‘limit’ which will eventually contain the answer. The optimization is minimal, but it saves a couple extra comparisons and variable swaps. Here is my code:
double Run()
{
limit = 600851475143;
double factor = 2;
double maxFactor = sqrt(limit);
while(factor == 2)
{
if(fmod(limit, factor) == 0)
limit /= factor;
else
factor += 1;
}
while(factor < limit && factor < maxFactor)
{
if(fmod(limit, factor) == 0){
limit /= factor;
maxFactor = sqrt(limit);
}
else
factor += 2;
}
return limit;
}
Results:
Limit Speed 600851475143 25 microseconds 9512811111 36 microseconds 600851475111 9 milliseconds
Obviously this method works very well. We shaved a ton of time off of our calculations. Comparing Optimization 2 against our Prime algorithm yields the following speed ratios (in order): 1:9, 1:3805, 1:633333
I think we can call it a day.
Optimization 3: So, it turns out we couldn’t call it a day. For whatever reason, I had just assumed that doubles would be the best variable to use since we were dealing with very large numbers. I neglected the 64-bit integer variable as a possible solution. Even though _int64 needs to be converted to double to use the sqrt function, it allows us to use the much faster % operator (as opposed to fmod). The results speak for themselves. The code:
double Run()
{
_int64 limit = 600851475111;
_int64 factor = 2;
_int64 maxFactor = sqrt((double)limit);
while(factor == 2)
{
if(limit % factor == 0)
limit /= factor;
else
factor += 1;
}
while(factor < limit && factor < maxFactor)
{
if(limit % factor == 0){
limit /= factor;
maxFactor = sqrt((double)limit);
}
else
factor += 2;
}
return limit;
}
Results:
Limit Speed 600851475143 9 microseconds 9512811111 13 microseconds 600851475111 2 milliseconds
As you can see, we have speed things up significantly. Now, dare I say, we can call it a day (again)?
Try looking into “sieve of eratosthenes”
The sieve is a great way to find the primes, I am just not sure that it would be necessarily applicable to this problem. If I attempted to write the algorithm finding the primes first I would think it would take longer. For instance, I would have to find and mark the primes, then iterate through and see if each prime divided the number I am looking for evenly, then, even if it did, I would need to keep going looking for the biggest to do so. I may give it a whirl this weekend. Maybe I will be surprised.
You may like this: http://coding.pressbin.com/70/List-of-primes-as-CSV-for-PHP-etc/
list of the first 1000 primes, comma separated. I did a dirty little hack method, by having an array of prime numbers.
Essentially it would try dividing by a prime. If it didn’t divide, it would move to the next highest prime number in the array. If it did divide, then it would reset the array back to 0, and begin looping again.
This way, I’d ensure it would only divide the number by the smallest prime numbers.
My only issue with this method was that it just flat out will not work if I use the square root of the number.
The code that works:
#include “stdafx.h”
#include
#include
using namespace std;
int _tmain(int argc, _TCHAR* argv[])
{
long aNumber = 600851475143;
int divider = 0;
int primes[1000] = { 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 131, 137, 139, 149, 151, 157, 163, 167, 173, 179, 181, 191, 193, 197, 199, 211, 223, 227, 229, 233, 239, 241, 251, 257, 263, 269, 271, 277, 281, 283, 293, 307, 311, 313, 317, 331, 337, 347, 349, 353, 359, 367, 373, 379, 383, 389, 397, 401, 409, 419, 421, 431, 433, 439, 443, 449, 457, 461, 463, 467, 479, 487, 491, 499, 503, 509, 521, 523, 541, 547, 557, 563, 569, 571, 577, 587, 593, 599, 601, 607, 613, 617, 619, 631, 641, 643, 647, 653, 659, 661, 673, 677, 683, 691, 701, 709, 719, 727, 733, 739, 743, 751, 757, 761, 769, 773, 787, 797, 809, 811, 821, 823, 827, 829, 839, 853, 857, 859, 863, 877, 881, 883, 887, 907, 911, 919, 929, 937, 941, 947, 953, 967, 971, 977, 983, 991, 997, 1009, 1013, 1019, 1021, 1031, 1033, 1039, 1049, 1051, 1061, 1063, 1069, 1087, 1091, 1093, 1097, 1103, 1109, 1117, 1123, 1129, 1151, 1153, 1163, 1171, 1181, 1187, 1193, 1201, 1213, 1217, 1223, 1229, 1231, 1237, 1249, 1259, 1277, 1279, 1283, 1289, 1291, 1297, 1301, 1303, 1307, 1319, 1321, 1327, 1361, 1367, 1373, 1381, 1399, 1409, 1423, 1427, 1429, 1433, 1439, 1447, 1451, 1453, 1459, 1471, 1481, 1483, 1487, 1489, 1493, 1499, 1511, 1523, 1531, 1543, 1549, 1553, 1559, 1567, 1571, 1579, 1583, 1597, 1601, 1607, 1609, 1613, 1619, 1621, 1627, 1637, 1657, 1663, 1667, 1669, 1693, 1697, 1699, 1709, 1721, 1723, 1733, 1741, 1747, 1753, 1759, 1777, 1783, 1787, 1789, 1801, 1811, 1823, 1831, 1847, 1861, 1867, 1871, 1873, 1877, 1879, 1889, 1901, 1907, 1913, 1931, 1933, 1949, 1951, 1973, 1979, 1987, 1993, 1997, 1999, 2003, 2011, 2017, 2027, 2029, 2039, 2053, 2063, 2069, 2081, 2083, 2087, 2089, 2099, 2111, 2113, 2129, 2131, 2137, 2141, 2143, 2153, 2161, 2179, 2203, 2207, 2213, 2221, 2237, 2239, 2243, 2251, 2267, 2269, 2273, 2281, 2287, 2293, 2297, 2309, 2311, 2333, 2339, 2341, 2347, 2351, 2357, 2371, 2377, 2381, 2383, 2389, 2393, 2399, 2411, 2417, 2423, 2437, 2441, 2447, 2459, 2467, 2473, 2477, 2503, 2521, 2531, 2539, 2543, 2549, 2551, 2557, 2579, 2591, 2593, 2609, 2617, 2621, 2633, 2647, 2657, 2659, 2663, 2671, 2677, 2683, 2687, 2689, 2693, 2699, 2707, 2711, 2713, 2719, 2729, 2731, 2741, 2749, 2753, 2767, 2777, 2789, 2791, 2797, 2801, 2803, 2819, 2833, 2837, 2843, 2851, 2857, 2861, 2879, 2887, 2897, 2903, 2909, 2917, 2927, 2939, 2953, 2957, 2963, 2969, 2971, 2999, 3001, 3011, 3019, 3023, 3037, 3041, 3049, 3061, 3067, 3079, 3083, 3089, 3109, 3119, 3121, 3137, 3163, 3167, 3169, 3181, 3187, 3191, 3203, 3209, 3217, 3221, 3229, 3251, 3253, 3257, 3259, 3271, 3299, 3301, 3307, 3313, 3319, 3323, 3329, 3331, 3343, 3347, 3359, 3361, 3371, 3373, 3389, 3391, 3407, 3413, 3433, 3449, 3457, 3461, 3463, 3467, 3469, 3491, 3499, 3511, 3517, 3527, 3529, 3533, 3539, 3541, 3547, 3557, 3559, 3571, 3581, 3583, 3593, 3607, 3613, 3617, 3623, 3631, 3637, 3643, 3659, 3671, 3673, 3677, 3691, 3697, 3701, 3709, 3719, 3727, 3733, 3739, 3761, 3767, 3769, 3779, 3793, 3797, 3803, 3821, 3823, 3833, 3847, 3851, 3853, 3863, 3877, 3881, 3889, 3907, 3911, 3917, 3919, 3923, 3929, 3931, 3943, 3947, 3967, 3989, 4001, 4003, 4007, 4013, 4019, 4021, 4027, 4049, 4051, 4057, 4073, 4079, 4091, 4093, 4099, 4111, 4127, 4129, 4133, 4139, 4153, 4157, 4159, 4177, 4201, 4211, 4217, 4219, 4229, 4231, 4241, 4243, 4253, 4259, 4261, 4271, 4273, 4283, 4289, 4297, 4327, 4337, 4339, 4349, 4357, 4363, 4373, 4391, 4397, 4409, 4421, 4423, 4441, 4447, 4451, 4457, 4463, 4481, 4483, 4493, 4507, 4513, 4517, 4519, 4523, 4547, 4549, 4561, 4567, 4583, 4591, 4597, 4603, 4621, 4637, 4639, 4643, 4649, 4651, 4657, 4663, 4673, 4679, 4691, 4703, 4721, 4723, 4729, 4733, 4751, 4759, 4783, 4787, 4789, 4793, 4799, 4801, 4813, 4817, 4831, 4861, 4871, 4877, 4889, 4903, 4909, 4919, 4931, 4933, 4937, 4943, 4951, 4957, 4967, 4969, 4973, 4987, 4993, 4999, 5003, 5009, 5011, 5021, 5023, 5039, 5051, 5059, 5077, 5081, 5087, 5099, 5101, 5107, 5113, 5119, 5147, 5153, 5167, 5171, 5179, 5189, 5197, 5209, 5227, 5231, 5233, 5237, 5261, 5273, 5279, 5281, 5297, 5303, 5309, 5323, 5333, 5347, 5351, 5381, 5387, 5393, 5399, 5407, 5413, 5417, 5419, 5431, 5437, 5441, 5443, 5449, 5471, 5477, 5479, 5483, 5501, 5503, 5507, 5519, 5521, 5527, 5531, 5557, 5563, 5569, 5573, 5581, 5591, 5623, 5639, 5641, 5647, 5651, 5653, 5657, 5659, 5669, 5683, 5689, 5693, 5701, 5711, 5717, 5737, 5741, 5743, 5749, 5779, 5783, 5791, 5801, 5807, 5813, 5821, 5827, 5839, 5843, 5849, 5851, 5857, 5861, 5867, 5869, 5879, 5881, 5897, 5903, 5923, 5927, 5939, 5953, 5981, 5987, 6007, 6011, 6029, 6037, 6043, 6047, 6053, 6067, 6073, 6079, 6089, 6091, 6101, 6113, 6121, 6131, 6133, 6143, 6151, 6163, 6173, 6197, 6199, 6203, 6211, 6217, 6221, 6229, 6247, 6257, 6263, 6269, 6271, 6277, 6287, 6299, 6301, 6311, 6317, 6323, 6329, 6337, 6343, 6353, 6359, 6361, 6367, 6373, 6379, 6389, 6397, 6421, 6427, 6449, 6451, 6469, 6473, 6481, 6491, 6521, 6529, 6547, 6551, 6553, 6563, 6569, 6571, 6577, 6581, 6599, 6607, 6619, 6637, 6653, 6659, 6661, 6673, 6679, 6689, 6691, 6701, 6703, 6709, 6719, 6733, 6737, 6761, 6763, 6779, 6781, 6791, 6793, 6803, 6823, 6827, 6829, 6833, 6841, 6857, 6863, 6869, 6871, 6883, 6899, 6907, 6911, 6917, 6947, 6949, 6959, 6961, 6967, 6971, 6977, 6983, 6991, 6997, 7001, 7013, 7019, 7027, 7039, 7043, 7057, 7069, 7079, 7103, 7109, 7121, 7127, 7129, 7151, 7159, 7177, 7187, 7193, 7207, 7211, 7213, 7219, 7229, 7237, 7243, 7247, 7253, 7283, 7297, 7307, 7309, 7321, 7331, 7333, 7349, 7351, 7369, 7393, 7411, 7417, 7433, 7451, 7457, 7459, 7477, 7481, 7487, 7489, 7499, 7507, 7517, 7523, 7529, 7537, 7541, 7547, 7549, 7559, 7561, 7573, 7577, 7583, 7589, 7591, 7603, 7607, 7621, 7639, 7643, 7649, 7669, 7673, 7681, 7687, 7691, 7699, 7703, 7717, 7723, 7727, 7741, 7753, 7757, 7759, 7789, 7793, 7817, 7823, 7829, 7841, 7853, 7867, 7873, 7877, 7879, 7883, 7901, 7907, 7919 };
while (divider < 999)
{
if (aNumber % primes[divider] == 0)
{
aNumber /= primes[divider];
divider = 0;
}
else divider++;
cout << "aNumber is " << aNumber << endl;
cout << "Divider is " << divider << endl;
cout << "Primes is " << primes[divider] << endl;
if (aNumber == 1)
break;
};
return 0;
}
Hey I know this is off topic but I was wondering if you knew of any
widgets I could add to my blog that automatically tweet my newest twitter updates.
I’ve been looking for a plug-in like this for
quite some time and was hoping maybe you would have some experience with something like this.
Please let me know if you run into anything. I truly enjoy reading your blog and I look forward to your
new updates.
Thanks for finally talking about > Project Euler Question #3 | | Fix By ProximityFix By Proximity < Liked it!
Hey there! I could have sworn I’ve been to this blog before
but after browsing through some of the post I realized it’s new to me.
Anyhow, I’m definitely delighted I found it
and I’ll be book-marking and checking back frequently!
Undeniably consider that that you stated.
Your favorite justification seemed to be on the internet the easiest factor to be
aware of. I say to you, I definitely get irked whilst
people think about worries that they just don’t
recognize about. You controlled to hit the nail upon the top and outlined out the
whole thing without having side effect , other people could take a signal.
Will probably be again to get more. Thank you
My blog: develop mobile apps for iphone and android (Wm)