小练习：20170828 分数的弹性

冬雪雪冬

从现在开始我们要开展一批欧拉计划的习题练习。

其实在我们论坛中已有欧拉计划的板块，可能有些鱼油还没注意到。

什么是欧拉计划：http://bbs.fishc.com/thread-60405-1-1.html

我们欧拉板块现已给出了200余题，这批练习将从欧拉计划中选题。其实用python语言完成有很多的优势，可以更简洁更方便的实现。

如果大家有兴趣也可浏览欧拉的英文网站：https://projecteuler.net/archives

这里已经有了500余题。

                               
登录/注册后可看大图

要求：

以python语言完成，如果是python2请注明。

程序以代码文字格式发帖。

注重程序效率和创意。

答题在一周内完成，即9.4 10：00之前，其后将公开大家的答案，并评比成绩。

另程序和答案可以在网上搜到，希望大家独立完成。题目不难，大家看看谁的效率高。

-- 回帖需写明解题思路，鼓励在程序中加上注释 --

-- 一些鱼油反映题目有些过难，为此略过一部分偏难的题目 --

                               
登录/注册后可看大图

243

jerryxjr1220

欧拉计划的第243题，可以解锁一个新成就


这题我写了2种解法，第一种是Bruce Force，利用欧拉函数硬算。
欧拉函数在“欧拉计划”的第69题有说明：传送门
R（d）=Phi（d）/（d-1）

利用numba的jit加速运算，不过还是很慢，525s才出结果，毕竟算了10亿以下的所有数字了。

1. import numpy as np
   
2. from numba import jit
   
3. @jit
   
4. def euler243(target, limit):
   
5.         primes = np.ones(limit, dtype='int8')
   
6.         for i in range(2, int(limit**0.5+1)):
   
7.                 primes[i*i::i] = 0
   
8.         plist=np.nonzero(primes)[0][2:].tolist()
   
9.         PHI = np.zeros(limit+1, dtype='int32')
   
10.         PHI[:2:1] = -1
    
11.         for p in plist:
    
12.                 PHI[p] = p - 1
    
13.         for p in plist:
    
14.                 for k in range(p*p, limit+1, p):
    
15.                         PHI[k] = PHI[k//p]*(p-1)
    
16.                 for j in range(p*p, limit+1, p*p):
    
17.                         PHI[j] = PHI[j//p]*p
    
18.         x = 2
    
19.         while (PHI[x]/(x-1)) >= target:
    
20.                 x += 1
    
21.         else:
    
22.                 print(x)
    
23.                 return
    
24. euler243(15499/94744, 1000000000)

复制代码

另外一种，找规律，我利用上述的欧拉函数列出了1000以下所有使R（d）减小的d。

1. 4 0.666666666667 #2的倍数
   
2. 6 0.4 #2*3的倍数
   
3. 12 0.363636363636
   
4. 18 0.352941176471
   
5. 24 0.347826086957
   
6. 30 0.275862068966 #2*3*5的倍数
   
7. 60 0.271186440678
   
8. 90 0.269662921348
   
9. 120 0.268907563025
   
10. 150 0.268456375839
    
11. 180 0.268156424581
    
12. 210 0.22966507177 #2*3*5*7的倍数
    
13. 420 0.229116945107
    
14. 630 0.22893481717
    
15. 840 0.22884386174

复制代码

从规律中，可以发现所有使R（d）减小的d都是“2*3*5*..."的倍数，而这个倍数i必然小于下一个质数。
利用这个性质：

1. def solve(limit):
   
2.     def R(d):
   
3.         num = d
   
4.         for prime in [x for x in primes if not d % x]:
   
5.             num -= num // prime
   
6.         return num / (d - 1)
   
7.     primes = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37]
   
8.     d = 1
   
9.     for prime in primes:
   
10.         d *= prime
    
11.         for i in range(1, prime):
    
12.             if R(d*i) < limit:
    
13.                 return d*i
    
14. print (solve(15499 / 94744))

复制代码

892371480

gunjang

1. from math import sqrt
   
2. 
   
3. '''
   
4. 欧拉φ函数：φ(n)是所有小于n的正整数里，和n互素的整数的个数。n是一个正整数。
   
5. 欧拉证明了下面这个式子：
   
6. 如果n的标准素因子分解式是p1^a1*p2^a2*……*pm^am，其中众pj(j=1,2,……,m)都是素数，而且两两不等。则有
   
7. φ(n)=n(1-1/p1)(1-1/p2)……(1-1/pm)2)....(1-1/pt)
   
8. #Rd = n/(n-1)*(1-1/p1)(1-1/p2)....(1-1/pt)
   
9. '''
   
10. def genPrime():
    
11.         for n in [2,3,5,7]:
    
12.                 yield n
    
13.         n = 11
    
14.         while True:
    
15.                 for i in range(3, int(sqrt(n))+1, 2):
    
16.                         if n % i == 0:
    
17.                                 break
    
18.                 if n % i != 0:
    
19.                         yield n
    
20.                 n += 2
    
21. 
    
22. n = 1
    
23. r = 1
    
24. for p in genPrime():
    
25.         n *= p
    
26.         r *= (1-1/p)
    
27.         if r*(n/(n-1)) < 15499/94744:
    
28.                 print(p, n) #29 6469693230
    
29.                 break
    

复制代码

6469693230

chunchun2017

答案是： 892371480，对应的R(d)=0.16358819553891188<15499/94744(=0.1635881955585578)
=============
思路：
R(n)=弹性分数个数/真分数个数=弹性分数个数/(n-1)=小于n且与n互质的数的个数/(n-1)
网上找了一下，求小于给定的n且与互质的数的个数，可以通过对欧拉函数算法稍加修改实现
初步运行看了一下100以内的数的R(n)，再加上题目中的12的提示，发现了一个规律：
0<n<5时，最小的R(n)为R(4)
0<n<10时，最小的R(n)为R(6)
0<n<15时，最小的R(n)为R(12)
0<n<20时，最小的R(n)为R(18)
0<n<25时，最小的R(n)为R(24)
0<n<40时，最小的R(n)为R(30)
0<n<80时，最小的R(n)为R(60)
0<n<100时，最小的R(n)为R(90)
而
4=2*2
6=2*3
12=2^2*3
18=2*3^2
24=2^3*3
30=2*3*5
60=2^2*3*5
90=2*3^2*5
可见，R(n)最小值都是出现在当n等于[2,3,5,7,11...]这些素数从小到大累积相乘时
为了快速锁定可能满足条件的区间，循环时，不考虑同一素数的幂，而直接采用2*3*5*7这样的方式递增，第一次发现素数的乘积n对应的R(n)小于15499/94744时，马上跳出循环。此时的n(n=p1*p2*p3..*pn）虽然满足R(n)<15499/94744,但不一定是最小值，因为n1=p1*p2*p3...*p(n-1)肯定是不符合，而n=p1*p2*p3...*p(n)肯定是符合条件的,而对于(n1,n2]这个区间中，可能会有符合条件的，也可能不会有符合条件的，为止，需要进一步确定：
(n1,n2]这个区间范围仍然很大，逐一计算耗时太长，根据前面观察到的规律可以分析出，（n1,n2]这个区间中如果有数n0满足R(n0)<15499/94744,这个数必定也是p1*p2*p3...（若干个素数的乘积，可能还会有幂的形式）因此，需要将前面求出的n值对应的素数前面的所有素数与n1滚动相乘，放到一个列表list1中，然后将其从小到大排序，只要找到第一个符合R(n0)<15499/94744的n0，马上跳出循环，此时的n0就是符合条件的最小n0（因为后面即使有ni满足R(ni)<15499/94744,ni也会比n0大）
代码如下：
=====================================

1. #count表示小于n的，与n互质的数的个数（包含1),同时count也表示n的弹性分数的个数
   
2. #R(n)=count/(n-1)
   
3. def PrintPrime(n): #输出[2,n]范围内所有的质数
   
4.     list0=[]   
   
5.     for i in range(2,n+1):
   
6.           for j in range(2,int(math.sqrt(i))+1):
   
7.             if(i%j==0):
   
8.                break
   
9.           else:
   
10.             list0.append(i)
    
11.     return list0
    
12. 
    
13. def euler(n): #count表示[1,n)范围内，不包含N，包含1，与互质的数的个数，返回值为R(n)
    
14.    n0 = n   
    
15.    sqrt_n = int(math.sqrt(n))
    
16.    count = n
    
17.    for i in range(2,sqrt_n+1):
    
18.       if(n%i==0):
    
19.          count = int(count/i*(i-1))
    
20.          while(n%i == 0):
    
21.             n/=i
    
22.    if(n>1):
    
23.        count = int(count/n*(n-1))
    
24.    return count/(n0-1)
    
25. 
    
26. import math
    
27. min0 = 15499/94744
    
28. list0 = PrintPrime(30)  #用于存放从2到30范围内的所有素数
    
29. sum = 1
    
30. for i in list0:   #先找到满足条件的N的上界
    
31.    sum*=i
    
32.    if(euler(sum)<min0):
    
33.        n0 = sum
    
34.        min = euler(sum)
    
35.        break;
    
36. list1=[int(sum/i),]#list1用于存放[sum/i,sum]范围内可能满足R(d)<15499/94744的数,list1[0]本来是不满足R(d)要求的，但为了计算方便，也放入list1
    
37. for k in list0[0:list0.index(i)]:
    
38.     for each in list1:
    
39.         sum1 = k*each
    
40.         if(sum1>=n0):
    
41.            break
    
42.         else:
    
43.            list1.append(sum1)
    
44.     list1.sort()
    
45. for n in list1[1:]: #跳过list1[0]
    
46.    if(euler(n)<min0):
    
47.       break
    
48. print('符合条件的最小N为: %d' % n)
    
49. print('%d对应的R(d)为:' %  n,euler(n))
    

复制代码

运行结果：
符合条件的最小N为: 892371480
892371480对应的R(d)为: 0.16358819553891188

gunjang

chunchun2017 发表于 2017-9-1 15:33
答案是： 892371480，对应的R(d)=0.16358819553891188

还是兄考虑周全，我没有考虑冥的情况

chunchun2017

gunjang 发表于 2017-9-5 11:12
还是兄考虑周全，我没有考虑冥的情况

虽然考虑到了，但是计算过程有点繁琐，比jerryxjr1220 大神的差远了