题目69：找出一百万以下的n中使得n/φ(n)取到最大的n

欧拉计划

Totient maximum

Euler's Totient function, φ(n) [sometimes called the phi function], is used to determine the number of numbers less than n which are relatively prime to n. For example, as 1, 2, 4, 5, 7, and 8, are all less than nine and relatively prime to nine, φ(9)=6.



It can be seen that n=6 produces a maximum n/φ(n) for n ≤ 10.

Find the value of n ≤ 1,000,000 for which n/φ(n) is a maximum.

题目：

欧拉函数φ(n)（有时也叫做 phi 函数）可以用来计算小于 n 的数字中与 n 互质的数字的个数。例如，因为 1,2,4,5,7,8 全部小于 9 并且与 9 互质，所以φ(9)=6。



可以看出对于n ≤ 10，n=6 时的 n/φ(n) 取到最大值。

找出 n ≤ 1,000,000 的 n 中使得 n/φ(n) 取到最大的 n 值。

jerryxjr1220

Mark一下：
解题的关键：
欧拉函数phi(m)：当m>1是，phi(m)表示比m小且与m互质的正整数个数
1、费马定理：a的p-1次方mod p余1。（其中p是素数，a是不能被p整除的正整数。
2、欧拉定理
    2.1 欧拉函数（RSA的证明用到）
          定义：欧拉函数phi(m)：当m>1是，phi(m)表示比m小且与m互质的正整数个数
               如：phi(24)=8    （1，5，7，11，13，17，19，23）
          性质：（1）当m为素数时，phi(m)=m-1
                    （2）当m=pq，且p和q是互异的素数，
                         则有：phi(m)=phi(p)*phi(p)=(p-1)(q-1) （这很好用）
                    （3）m＝p^e，且p为素数，e为正整数，则
                         phi(m)=p^e-p^(e-1)=(p^(e-1))*(p-1)
          定理：若m=p1^e1*p2^e2....pt^et则：
                phi(m)=m(1-1/p1)(1-1/p2)....(1-1/pt) （pi是素数）
     2.2 欧拉定理
            a^phi(n)=1 mod n
          注：[1] n=p时候，有a^(p-1)=1 mod p，为费马定理
                [2] a^(phi(n)+1)=a mod n
                [3] 若n=pq，p与q为相异素数，取大于0的m,n互质数，有
                    m^(phi(n)+1)=m mod n; m^((p-1)(q-1)+1)=m mod n
欧拉函数的定义:E(k)=([1,n-1]中与n互质的整数个数).
   
     因为任意正整数都可以唯一表示成如下形式:
                     k=p1^a1*p2^a2*……*pi^ai;(即分解质因数形式)
    可以推出:E(k)=(p1-1)(p2-1)……(pi-1)*(p1^(a1-1))(p2^(a2-1))……(pi^(ai-1))
               =k*(p1-1)(p2-1)……(pi-1)/(p1*p2*……pi);
               =k*(1-1/p1)*(1-1/p2)....(1-1/pk)
     ps：在程序中利用欧拉函数如下性质，可以快速求出欧拉函数的值(a为N的质因素)
      若(N%a==0 && (N/a)%a==0) 则有:E(N)=E(N/a)*a;
若(N%a==0 && (N/a)%a!=0) 则有:E(N)=E(N/a)*(a-1);
第一次写欧拉函数的题，琢磨的半天，最后还是只能按照最开始的想法写......
      欧拉函数PHI(n)表示的是比n小，并且与n互质的正整数的个数(包括1)。比如：
      PHI(1) = 1; PHI(2) = 1; PHI(3) = 2; PHI(4) = 2; ... PHI(9) = 6; ...
要计算一个正整数n的欧拉函数的方法如下：
1. 将n表示成素数的乘积: n = p1 ^ k1 * p2 ^ k2 * ... * pn ^ kn（这里p1, p2, ..., pn是素数）
2. PHI(n) = (p1 ^ k1 - p1 ^ (k1 - 1)) * (p2 ^ k2 - p2 ^ (k2 - 1)) * ... * (pn ^ kn - pn ^ (kn - 1))
              = Mult { pi ^ ki - pi ^ (ki -1) }
证明过程如下：
    1. 容易想到：当n为素数时，PHI(n) = n - 1。因为每个比n小的正整数都和n互素。当n为素数p的k次方时，PHI(n) = p ^ k - p ^ (k - 1)。因为在1到n之间的正整数只有p的倍数和n不互素，这样的数有(p ^ k / p)个（为什么？）。
    2. 如果m和n互素，即GCD(m, n) = 1，那么PHI(m * n) = PHI(m) * PHI(n)。用中国剩余定理可以证明，证明的思路是建立这样一种一一对应的关系(a, b) <-> x，其中正整数a小于m并且gcd(a, m) = 1，正整数b小于n并且gcd(b, n) = 1，正整数x小于m*n并且gcd(m*n, x) = 1。证明过程如下：
    1)根据中国剩余定理，如果m和n互素，那么关于未知量x的方程组x % m = a, x % n = b(0 <= a < m, 0 <= b < n)，当0 <= x < m * n时存在并且仅存在一个解。容易证明，如果两个这样的方程组有相同的m, n但是a, b不同，那么他们的解x一定不同。
    2)首先用反正法证明：gcd(m, a) = 1且gcd(n, b) = 1是gcd(m*n, x) = 1的必要条件：假设gcd(a, m) = k > 1，由此可得：a = a' * k; m = m' * k => x = k' * m + a = k' * k * m' + k * a' = k * (k' * m' + a'); 所以gcd(x, m) = k > 1。同理可证，如果gcd(b, n) > 1， 那么gcd(x, n) > 1。所以x和m * n互素的必要条件是a和m互诉且b和n互素。
    3)接下来我们证明充分性：由x % m = a 可以得到x = k * m + a；由欧几里德算法求最大公约数的过程（就不证明了，呵呵，还得想）可以知道gcd(x, m) = gcd(m, a) = 1；同理可得，如果gcd(n, b) = 1那么gcd(x, n) = 1。接下来很容易得到：gcd(m*n, x) = 1。从而证明了充分性。
    4)上面三步的结论表明，数对(a, b)是可以和x建立起一一对应的关系的，所以有多少个不同的(a, b)，就有多少个不同的x。
    3.将n分解成素数乘积后，显然对于任意的i, j(i != j)都满足 pi ^ ki和pj ^ kj是互素的，于是可以的到上面的公式。
跟据上面的公式，可以得到关于欧拉函数的递推关系：
假设素数p能整除n，那么
如果p还能整除n / p, PHI(n) = PHI(n / p) * p;
如果p不能整除n / p, PHI(n) = PHI(n / p) * (p - 1);

jerryxjr1220

解题思路：
用标记法求解，先生成PHI = [0]*1000000，然后依次用质数（或质数的倍数）替换PHI中的值，直至全部完成。 最后求解max(N/PHI)，可以用enumerate函数辅助。

jerryxjr1220

答案：
510510 5.539388020833333
[Finished in 2.9 sec]

def getPrime(n):
    primes = [True]*n
    primes[0], primes[1] = False, False
    for (i, prime) in enumerate(primes):
        if prime:
            for j in range(i*i,n,i):
                primes[j] = False
    primelist = []
    for (k, trueprime) in enumerate(primes):
        if trueprime:
            primelist.append(k)
    return primelist

def philist(n):
    primelist = getPrime(n)
    PHI = [0]*(n+1)
    PHI[0], PHI[1] = -1, -1
    for eachprime in primelist:
        PHI[eachprime] = eachprime - 1
    for i in range(len(primelist)):
        p = primelist[i]
        for k in range(p*p,n+1,p):
                PHI[k] = PHI[int(k/p)]*(p-1)
        for j in range(p*p,n+1,p*p):
            PHI[j] = PHI[int(j/p)]*p
    return PHI
plist = philist(1000001)

maxn, maxnphi = 0, 0
for i in range(2,1000001):
    if i/plist[i] > maxnphi:
        maxnphi = i/plist[i]
        maxn = i
print (maxn,maxnphi)

芒果加黄桃

jerryxjr1220 发表于 2016-10-16 08:17
答案：
510510 5.539388020833333
[Finished in 2.9 sec]

学习 感谢分享

najin

用的matlab

ans =

      510510

时间已过 0.001139 秒。
>>

debuggerzh

510510

Process returned 0 (0x0)   execution time : 3.888 s
Press any key to continue.
2#大大的递推思路很硬核，只是我没实现出来……
不论如何涨姿势了
转化为倒数的最小值，利用两个结论预先打部分表
1.p为质数时，有

                               
登录/注册后可看大图

2.正整数m，n互质时，有

                               
登录/注册后可看大图

#include<iostream>
#include<cmath>
#include<set>
#include<algorithm>
using namespace std;

const int M = 1e6;
const int M_RANGE = 1e6;
int ct = 0;
bool prime[M];
int factor[M];
int phi[M] = {0};

void euler_sieve(){
  prime[1] = false;
  for (int i = 2;i <= M;i++) prime[i] = true;

  for (int i = 2;i <= M;i++){
    if (prime[i]) {factor[++ct] = i;  phi[i] = i-1;}
    for (int j = 1;j <= ct && i*factor[j] < M;j++){
      prime[i*factor[j]] = false;
      if (i % factor[j] == 0) break;
    }
  }
}

void ini(){
  for (int i = 2;i <= sqrt(M);i++){
    if (prime[i]){
      for (int j = i;j*i <= M;j*=i){
        phi[j*i] = phi[j]*i;
      }
      for (int j = i;j*i <= M;j++){
        if (prime[j]) phi[j*i] = phi[j]*phi[i];
      }
    }
  }
}

void solve(set<int> & p,int n){
  for (int i = 1;n != 1;i++){
    int t = factor[i];
    while(n % t == 0) {p.insert(t); n /= t;}
  }
}

double quot(const set<int> & p){
  double res = 1;

  for (set<int>::iterator it = p.begin();it != p.end();++it){
    res *= *it - 1;
    res /= *it;
  }
  return res;
}

int main(){
  euler_sieve();
  ini();
  double mn = 1;//min
  int ans;

  for (int i = 2;i <= M;i++){
    set<int> p;
    double t;

    if (phi[i]) t = phi[i] / (double)i;
    else{
      solve(p,i);
      t = quot(p);
    }
    if (t < mn) {mn = t;  ans = i;}
    //cout << i << " " << t << endl;
  }
  cout << ans << endl;
  return 0;
}

debuggerzh

补图

                               
登录/注册后可看大图

                               
登录/注册后可看大图

永恒的蓝色梦想

C++ 12ms

#include<iostream>
#include<vector>
#include<cstring>
using namespace std;



template<size_t size, typename T>
void euler_sieve(bool (&is_prime)[size], vector<T>& primes, T (&phi)[size]) {
    T i, k;
    memset(is_prime, -1, size);
    is_prime[0] = is_prime[1] = 0;

    for (i = 2; i < size; i++) {
        if (is_prime[i]) {
            primes.push_back(i);
            phi[i] = i - 1;
        }

        for (T j : primes) {
            k = i * j;

            if (k >= size) {
                break;
            }

            is_prime[k] = false;

            if (!(i % j)) {
                phi[k] = phi[i] * j;
                break;
            }

            phi[k] = phi[i] * (j - 1);
        }
    }
}



int main() {
    ios::sync_with_stdio(false);

    constexpr unsigned int UPPER_BOUND = 1000000;
    static bool is_prime[UPPER_BOUND];
    static unsigned int phi[UPPER_BOUND];
    vector<unsigned int> primes;
    double max_div = 0, div;
    unsigned int max_n = 0, n;

    euler_sieve(is_prime, primes, phi);


    for (n = 2; n < UPPER_BOUND; n++) {
        div = (double)n / phi[n];

        if (div > max_div) {
            max_div = div;
            max_n = n;
        }
    }


    cout << max_n << endl;
    return 0;
}

guosl

我们可以应用筛法来批量求出各数的欧拉函数值。

/*
答案：510510
耗时：0.0148415秒
*/
#include <iostream>
#include <algorithm>
#include <cstring>
using namespace std;

const int nN = 1000001;
int f[nN];//记录φ(i)

int main(void)
{
  double d = 0.0;
  int m = 0;
  //计算φ(i)
  for (int i = 1; i < nN; ++i)
  {
    f[i] = i;
    if ((i & 1) == 0)
      f[i] /= 2;
  }
  for (int i = 3; i < nN; i += 2)
  {
    if (f[i] < i)
      continue;
    for (int j = i; j < nN; j += i)
      f[j] = f[j] / i * (i - 1);
  }
  for (int i = 2; i <= 1000000; ++i)
  {
    double td = double(i) / double(f[i]);
    if (td > d)
    {
      m = i;
      d = td;
    }
  }
  cout << m << endl;
  return 0;
}