I decided to do problem #69. I initially just thought I'd just go the naive route and brute force using the Euclidean algorithm for the greatest common denominator:
This of course turned out to be way too slow. But upon further reading it turns out the totient function can be calculated using Euler's product formula, which calculates the totient based on a number's prime factorization. Since I am going over the range [1,1000000] that means I only have to check <= pi(sqrt(1000000)) = 168 primes to factor an integer by trial division (pi(n) being the amount of primes <= n). Here's the code:
Which gave me the answer in about ~15 seconds. After solving this, I googled and apparently it can be solved in a much better way.
#include <stdio.h>
int
gcd (int a, int b)
{
return a == 0 ? b : gcd(b % a, a);
}
int
phi (int n)
{
int i, o = 0;
for (i = 1; i <= n; i++)
{
if (gcd(n, i) == 1)
{
o++;
}
}
return o;
}
int
main (void)
{
float max;
float cur;
int i;
int k;
max = 0.0;
for (i = 1; i <= 1000000; i++)
{
cur = (float)i/phi(i);
if (cur > max)
{
max = cur;
k = i;
}
printf("%d: %.6f\n", i, cur);
}
printf("\nMax: %d: %.6f\n", k, max);
return 0;
}
This of course turned out to be way too slow. But upon further reading it turns out the totient function can be calculated using Euler's product formula, which calculates the totient based on a number's prime factorization. Since I am going over the range [1,1000000] that means I only have to check <= pi(sqrt(1000000)) = 168 primes to factor an integer by trial division (pi(n) being the amount of primes <= n). Here's the code:
#include <stdio.h>
#define MIN(a,b) (((a)<(b))?(a):(b))
/* All primes under sqrt(n), n = 1,000,000 */
int primes[] = {
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
};
int
phi (int n)
{
int o = n;
int i;
for (i = 0; i*i <= MIN(n, 168); i++)
{
if (n % primes[i] == 0)
{
o -= o / primes[i];
}
}
return o;
}
int
main (void)
{
float max;
float cur;
int i;
int k;
max = 0.0;
for (i = 1; i <= 1000000; i++)
{
cur = (float)i/phi(i);
if (cur > max)
{
max = cur;
k = i;
}
printf("%d: %.6f\n", i, cur);
}
printf("\nMax: %d: %.6f\n", k, max);
return 0;
}
Which gave me the answer in about ~15 seconds. After solving this, I googled and apparently it can be solved in a much better way.