鱼C论坛

 找回密码
 立即注册
查看: 4718|回复: 9

[技术交流] 小练习:20171023 能被两个素数整除的最大整数

[复制链接]
发表于 2017-10-22 19:52:44 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能^_^

您需要 登录 才可以下载或查看,没有账号?立即注册

x
本帖最后由 冬雪雪冬 于 2017-10-31 20:17 编辑

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

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

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

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

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

这里已经有了500余题。


                               
登录/注册后可看大图


要求:

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

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

注重程序效率和创意。

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

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

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

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


                               
登录/注册后可看大图
347

在≤ 100的整数中,能整除的素数只有2和3的最大整数是96,因为96=32*3=25*3。对于两个不同的素数p和q,在≤N的正整数中,能整除的素数只有p和q的最大整数记为M(p,q,N);如果这样的正整数不存在,则记M(p,q,N)=0。
例如,M(2,3,100)=96。
M(3,5,100)=75而非90,因为90能被2、3、5整除。
此外M(2,73,100)=0,因为不存在≤ 100的正整数使得能整除的素数只有2和73。
记S(N)为所有不同的M(p,q,N)之和。S(100)=2262。
求S(10 000 000)。

想知道小甲鱼最近在做啥?请访问 -> ilovefishc.com
回复

使用道具 举报

发表于 2017-10-22 22:21:33 | 显示全部楼层
本帖最后由 gunjang 于 2017-10-23 20:40 编辑

答案11109800204052
原来花费25s,现在7s左右,本质还是暴力法,需要优化

  1. import time
  2. import numpy as np
  3. from math import log

  4. st = time.time()
  5. def solve(limit):
  6.         primes = np.ones(limit, dtype=np.bool)
  7.         for i in range(2, limit//2):
  8.                 if primes[i]: primes[i*i::i] = False
  9.         primes[0], primes[1] = False, False
  10.         plist = np.nonzero(primes)[0]

  11.         r = 0
  12.         for i in range(len(plist)-1):
  13.                 if plist[i] * plist[i+1] > limit:
  14.                         break
  15.                 for j in range(i+1, len(plist)):
  16.                         if plist[i] * plist[j] > limit:
  17.                                 break
  18.                         m = 0
  19.                         pa = plist[i] #p^a
  20.                         qb = plist[j]**int(log(limit/pa, plist[j]))        #q^b
  21.                         while True:
  22.                                 t = pa * qb
  23.                                 if t > limit:
  24.                                         qb //= plist[j]
  25.                                         if qb == 1:
  26.                                                 break
  27.                                         t = pa * qb
  28.                                 if t > m:
  29.                                         m = t
  30.                                 pa *= plist[i]
  31.                         r += m
  32.         return r #overflow in 32bit

  33. print(solve(10000000)) #11109800204052
  34. print('cost %fs'%(time.time()-st)) #10m cost 25.9s => 6.7s
复制代码

评分

参与人数 1荣誉 +10 鱼币 +10 收起 理由
冬雪雪冬 + 10 + 10

查看全部评分

想知道小甲鱼最近在做啥?请访问 -> ilovefishc.com
回复 支持 反对

使用道具 举报

发表于 2017-10-23 09:50:50 | 显示全部楼层
本帖最后由 jerryxjr1220 于 2017-10-23 15:01 编辑

一开始写的程序,就是完全按照题目逻辑来,暴力解法,1000万的话,差不多15秒出结果。
  1. import math
  2. def gen_primes(limit):
  3.     primes = [1]*limit
  4.     for i in xrange(2,int(limit**0.5+1)):
  5.         if primes[i]:
  6.             for j in xrange(i*i,limit,i):
  7.                 primes[j]=0
  8.     return [k for k,v in enumerate(primes) if v][2:]
  9. def get_max(p,q,n):
  10.     mp = int(math.log(n,p)+1)
  11.     mq = int(math.log(n,q)+1)
  12.     maxi = 0
  13.     for i in range(1,mp):
  14.         for j in range(1,mq):
  15.             if p**i*q**j>n: break
  16.             maxi = max(maxi, p**i*q**j)
  17.     return maxi
  18. def solve(N):
  19.     result = 0
  20.     plist=gen_primes(N//2)
  21.     for i in xrange(len(plist)-1):
  22.         for j in xrange(i+1,len(plist)):
  23.             gm = get_max(plist[i],plist[j],N)
  24.             if gm==0: break
  25.             result += gm
  26.     return result
  27. print solve(10000000)
复制代码

11109800204052
[Finished in 15.5s]
程序很容易写,逻辑清晰。

对这个程序做一下优化,get_max()函数分情况讨论,不用穷举所有组合,可以进一步缩短运算时间。
  1. import math
  2. def gen_primes(limit):
  3.     primes = [1]*limit
  4.     for i in xrange(2,int(limit**0.5+1)):
  5.         if primes[i]:
  6.             for j in xrange(i*i,limit,i):
  7.                 primes[j]=0
  8.     return [k for k,v in enumerate(primes) if v][2:]
  9. def get_max(p,q,n):
  10.     mp = int(math.log(n,p))
  11.     mq = int(math.log(n,q))
  12.     if p*q>n:
  13.         maxi = 0
  14.     elif q*q>n:
  15.         maxi = p**mp*q
  16.         while maxi>n:
  17.             maxi//=p
  18.     elif n**(1/3)<q<n**0.5 and (n/q)**0.5<p<q:
  19.         maxi = p*q
  20.     else:
  21.         maxi = 0
  22.         for i in range(1,mp+1):
  23.             for j in range(1,mq+1):
  24.                 if p**i*q**j>n: break
  25.                 maxi = max(maxi, p**i*q**j)
  26.     return maxi
  27. def solve(N):
  28.     result = 0
  29.     plist=gen_primes(N//2)
  30.     for i in xrange(len(plist)-1):
  31.         for j in xrange(i+1,len(plist)):
  32.             gm = get_max(plist[i],plist[j],N)
  33.             if gm==0: break
  34.             result += gm
  35.     return result
  36. print solve(10000000)
复制代码

11109800204052
[Finished in 7.6s]

用arrdio改写,并做个界面:
  1. import win.ui;
  2. import time;

  3. /*DSG{{*/
  4. var winform = win.form(text="Project Euler Calculation";right=460;bottom=278;border="thin";max=false)
  5. winform.add(
  6. button={cls="button";text="Calculate";left=120;top=182;right=334;bottom=233;font=LOGFONT(h=-27;weight=700);z=3};
  7. edit={cls="edit";left=135;top=117;right=423;bottom=158;edge=1;font=LOGFONT(h=-23;weight=700);multiline=1;z=2};
  8. limit={cls="edit";left=136;top=69;right=423;bottom=110;edge=1;font=LOGFONT(h=-23;weight=700);multiline=1;z=7};
  9. static={cls="static";text="Result";left=7;top=117;right=133;bottom=153;align="center";font=LOGFONT(h=-23;weight=700);transparent=1;z=1};
  10. static2={cls="static";text="Euler 347";left=21;top=21;right=216;bottom=58;font=LOGFONT(h=-27;weight=700;italic=255);transparent=1;z=4};
  11. static3={cls="static";text="Limit";left=7;top=72;right=133;bottom=108;align="center";font=LOGFONT(h=-23;weight=700);transparent=1;z=6};
  12. time_calc={cls="static";left=10;top=249;right=438;bottom=276;font=LOGFONT(h=-17);transparent=1;z=5}
  13. )
  14. /*}}*/

  15. var calc = function (limit) {
  16.         var gen_primes = function (n) {
  17.                 var primes = {};
  18.                 for i=1;n table.push(primes,1);
  19.                 for i=2;n**0.5 {
  20.                         if primes[i] {
  21.                                 for j=i*i;n;i primes[j]=0;
  22.                         }
  23.                 }
  24.                 var plist = {};
  25.                 for k,v in primes if v table.push(plist,k);
  26.                 table.remove(plist,1);
  27.                 return plist;
  28.         }

  29.         var get_max = function (p,q,n) {
  30.                 var mp=math.floor(math.log10(n)/math.log10(p));
  31.                 var mq=math.floor(math.log10(n)/math.log10(q));
  32.                 if p*q>n return 0;
  33.                 elseif q*q>n {
  34.                         max = p**mp*q;
  35.                         while max>n max/=p;
  36.                         return max;
  37.                 }
  38.                 elseif n**(1/3)<q and q<n**0.5 and (n/q)**0.5<p and p<q return p*q;
  39.                 else {
  40.                         max = 0;
  41.                         for i=1;mp {
  42.                                 for j=1;mq {
  43.                                         if p**i*q**j>n break;
  44.                                         max = math.max(max,p**i*q**j);
  45.                                 }
  46.                         }
  47.                         return max;
  48.                 }
  49.         }
  50.         var result = 0;
  51.         var plst = gen_primes(limit/2);
  52.         for x=1;#plst-1 {
  53.                 for y=x+1;#plst {
  54.                         var tmp = get_max(plst[x],plst[y],limit);
  55.                         if not tmp break;
  56.                         result += tmp;
  57.                 }
  58.         }
  59.         return result;
  60. }
  61. winform.button.oncommand = function(id,event){
  62.         //winform.msgbox( winform.button.text );
  63.         var st = time.tick();
  64.        
  65.         winform.edit.text = calc ( eval(winform.limit.text) );
  66.         winform.time_calc.text = string.format('Total time used: %d ms',time.tick()-st);
  67. }
  68. winform.show()
  69. win.loopMessage();
复制代码

Untitled.png
3.5秒出结果,同样的逻辑,果然还是python的运算效率拼不过啊。

换一种逻辑,不是依次测试每个质数对,而是改用标记法,把每个质数的倍数都标记出来。
  1. N=10**7
  2. def gen_primes(limit):
  3.     primes = [1]*limit
  4.     for i in xrange(2,int(limit**0.5+1)):
  5.         if primes[i]:
  6.             for j in xrange(i*i,limit,i):
  7.                 primes[j]=0
  8.     return [k for k,v in enumerate(primes) if v][2:]
  9. plist=gen_primes(N)
  10. multiprime=[[i] for i in xrange(N+1)]
  11. for p in plist:
  12.     for q in xrange(p,N+1,p):
  13.         if not multiprime[q]: continue
  14.         if len(multiprime[q]) >= 3:
  15.             multiprime[q] = False
  16.         else: multiprime[q].append(p)
  17. used = set()
  18. result = 0
  19. for i in xrange(N,-1,-1):
  20.     if not multiprime[i]: continue
  21.     if len(multiprime[i]) == 3:
  22.         mark = (multiprime[i][1], multiprime[i][2])
  23.         if mark not in used:
  24.             used.add(mark)
  25.             result += i
  26. print result
复制代码

差不多20多秒出结果。

评分

参与人数 1荣誉 +10 鱼币 +10 收起 理由
冬雪雪冬 + 10 + 10

查看全部评分

想知道小甲鱼最近在做啥?请访问 -> ilovefishc.com
回复 支持 反对

使用道具 举报

发表于 2017-10-23 22:55:13 | 显示全部楼层
  1. import math

  2. #素数判断函数
  3. def is_prime(n):
  4.   if n < 2:
  5.     return False
  6.   for i in range(2, int(n**0.5+1)):
  7.     if n%i == 0:
  8.       return False
  9.   return True

  10. #对于两个不同的素数p和q,在≤N的正整数中,
  11. #能整除的素数只有p和q的最大整数记为M(p,q,N);
  12. #如果这样的正整数不存在,则记M(p,q,N)=0
  13. def M(p,q,N):
  14.   lst = []
  15.   if p * q > N:
  16.     return 0        #如果这样的正整数不存在,则记M(p,q,N)=0
  17.   else:
  18.     for i in range(1,int(math.log(N,p))+1):
  19.       for j in range(1,int(math.log(N/p**i,q))+1):
  20.         lst.append(p**i * q**j)   #生成p,q两个素数正整数幂运算,小于N的所有结果
  21.     return max(lst)

  22. # S(N)为所有不同的M(p,q,N)之和
  23. def S(N):
  24.     prime_list = [x for x in range(1,N//2) if is_prime(x)] #通过N//2优化掉一部分pq组合
  25.     total = 0
  26.     for m in range(len(prime_list)-1):
  27.         for n in range(m+1,len(prime_list)):
  28.             if prime_list[m] * prime_list[n] > N: #通过p*q不能大于N再优化掉一部分pq组合
  29.                 break
  30.             else:
  31.                 total += M(prime_list[m],prime_list[n],N)
  32.     return total
复制代码


运行结果:
>>> S(100)
2262
>>> S(10000000)
11109800204052

评分

参与人数 1荣誉 +10 鱼币 +10 收起 理由
冬雪雪冬 + 10 + 10

查看全部评分

想知道小甲鱼最近在做啥?请访问 -> ilovefishc.com
回复 支持 反对

使用道具 举报

发表于 2017-10-24 13:17:47 | 显示全部楼层
jerryxjr1220 发表于 2017-10-23 09:50
一开始写的程序,就是完全按照题目逻辑来,暴力解法,1000万的话,差不多15秒出结果。

11109800204052
3.5秒出结果,同样的逻辑,果然还是python的运算效率拼不过啊。


lua 和 python的底层核心原理是不同的,这是它们之间效率差别的根源。
lua 刚开始是堆栈虚拟机,后来改为 寄存器虚拟机,用保守的c89标准写c程序。因为多数系统都支持c89,lua在源码级别就实现了跨平台。
py自始至终都是堆栈虚拟机,所有的字节码都是在操作栈结构,目的就是用简单的数据结构达到跨平台。
堆栈虚拟机天然地比 寄存器虚拟机慢,栈结构里,所有的操作都是串行的。

点评

我很赞同!: 1.0
我很赞同!: 1
昨天我用lua和arrdio做了个对比测试,都是循环10亿次累加,lua环境用时30多秒,arrdio用时才10秒,不明白为啥lua反而慢那么多?用python跑花了200秒...  发表于 2017-10-26 08:45

评分

参与人数 1荣誉 +5 鱼币 +5 贡献 +5 收起 理由
jerryxjr1220 + 5 + 5 + 5 学习了!

查看全部评分

想知道小甲鱼最近在做啥?请访问 -> ilovefishc.com
回复 支持 反对

使用道具 举报

发表于 2017-10-24 14:59:48 | 显示全部楼层
本帖最后由 JonTargaryen 于 2017-10-24 15:52 编辑

最后没有得到正确结果,测试100可以得到2262
10000000得到的结果是11109740320590

1. 得到n = 10000000以内的全部素数列表
2. 遍历之(大于sqrt(n)可以退出)
3. 得到每两个素数组合的结果
4. 求和

求指导

代码
  1. import math

  2. # get the list of all the primes that smaller than n
  3. def all_prime(n):
  4.     primes = [0] * (n+1)
  5.     primes[2] = 1
  6.     primes[3] = 1
  7.     primes[5] = 1

  8.     result = [2]

  9.     for num in range(1, n + 1, 6):
  10.         primes[num] = 1;

  11.     for num in range(5, n + 1, 6):
  12.         primes[num] = 1;

  13.     for num in range(3, int(math.sqrt(n) + 1), 2):
  14.         if primes[num]:
  15.             for j in range(num + num, n + 1, num):
  16.                 primes[j] = 0

  17.     for num in range(3, n + 1, 2):
  18.         if primes[num]:
  19.             result.append(num);

  20.     return result


  21. def largest_divisible(p, q, n) :
  22.     num = p * q  # the answer must contain both p and q
  23.     if num > n:
  24.         return 0
  25.     elif num == n:
  26.         return n

  27.     goal = n // num  # what we want is to approach the "goal" as much as possible

  28.     # if q >= sqrt(n), then n contains q just once
  29.     if (q >= math.sqrt(n)):
  30.         temp = int(math.log(goal, p))

  31.         return num * (p ** temp)

  32.     diff_min = goal  # minimize the diff and approach the goal
  33.     result = 1

  34.     for i in range(int(math.log(goal, p)) + 1):
  35.         temp_p = p ** i

  36.         for j in range(int(math.log(goal, q)) + 1):
  37.             temp_q = q ** j
  38.             temp = temp_p * temp_q

  39.             if temp * p < goal:  # too small
  40.                 continue

  41.             if temp > goal:
  42.                 break
  43.             elif temp == goal:
  44.                 return temp * num
  45.             else:
  46.                 diff_temp = goal - temp
  47.                 if diff_temp < diff_min:
  48.                     result = temp
  49.                     diff_min = diff_temp
  50.     return num * result


  51. def main():
  52.     n = int(input("please input a number:"))
  53.     result = 0

  54.     primes = all_prime(n // 2)  # get the list of all the primes smaller than n // 2 + 1

  55.     # p is the smaller prime, so it shoule be smaller than sqrt(n)
  56.     for p in primes:
  57.         if p > math.sqrt(n):
  58.             break

  59.         for q in primes:
  60.             if q > p:  # q is the larger prime
  61.                 result += largest_divisible(p, q, n)

  62.     print(result)


  63. if __name__ == "__main__":
  64.     main()
复制代码
想知道小甲鱼最近在做啥?请访问 -> ilovefishc.com
回复 支持 反对

使用道具 举报

发表于 2017-10-25 16:18:31 | 显示全部楼层
本帖最后由 小Q学Python 于 2017-10-25 16:22 编辑
  1. def primeslist(max):
  2.     a = [True]*(max+1)
  3.     a[0],a[1]=False,False
  4.     for index in range(2, len(a)):
  5.         if a[index]==True:
  6.             for i in range(index+index,len(a),index):
  7.                 a[i]=False
  8.     return [i for i, j in enumerate(a) if j==True]

  9. n = 10000000
  10. pl = primeslist(n//2)
  11. l = []
  12. for i in range(len(pl)-1):
  13.     for j in range(i+1, len(pl)):
  14.         if pl[i]*pl[j] > n:
  15.             break
  16.         x = [pl[i]*pl[j]]
  17.         for k in x:
  18.             if k*pl[i] <= n:
  19.                 x.append(k*pl[i])
  20.         for k in x:
  21.             if k*pl[j] <= n:
  22.                 x.append(k*pl[j])
  23.         l.append(max(x))
  24. print(sum(l))
复制代码



11109800204052

评分

参与人数 1荣誉 +10 鱼币 +10 收起 理由
冬雪雪冬 + 10 + 10

查看全部评分

想知道小甲鱼最近在做啥?请访问 -> ilovefishc.com
回复 支持 反对

使用道具 举报

发表于 2017-10-25 16:54:20 | 显示全部楼层
代码写完了,前排占个坑,性能很差。。。。
  1. import time


  2. #判断X是不是质数,zs为小于x的质数列表
  3. def FindZS(x,zs):#一个数(N)不能被小于根号N的数整除,那么这个数就是质数。
  4.     for j in zs:
  5.         if j > x**0.5:#所以当除数大于根号N的数时,跳出循环
  6.             break
  7.         if x%j == 0: #余数为0时,这个数就是合数返回False
  8.             return False
  9.     return True

  10. #返回X以内的质数列表,6N±1法则
  11. def zhisu(x):
  12.     z = [2,3]
  13.     for i in range(6,x,6):
  14.         if FindZS(i-1,z) is True:
  15.             z.append(i-1)
  16.         if FindZS(i+1,z) is True:
  17.             z.append(i+1)
  18.     return z

  19. #分解质因数
  20. def F(n):
  21.     global zs_half
  22.     if n <= 1 :
  23.         return
  24.     f = set()  #存放质因数的无重复列表
  25.     while n > 1:
  26.         x = n
  27.         for i in zs_half:
  28.             if n % i == 0:  #若能整除
  29.                 f.add(i)  #则表示i是质因数
  30.                 n = n//i  #N/i返回N重新进入判断
  31.                 break  #找到1个质因数后马上break,跳出当前for循环
  32.         if x == n:#判断X是不是质数
  33.             break
  34.     if len(f) == 0:
  35.         return
  36.     return f


  37. def M(q,p,N):#计算q、p在N以内的最大整数
  38.     global zs_all
  39.     if q == p or p*q>N:#如果p等于q,或者q和p的积大于N那么返回0
  40.         return 0
  41.     else:
  42.         while N:#循环递减N
  43.             if N not in zs_all:
  44.                 a = F(N)
  45.                 if a != None:
  46.                     if set([q,p]) == a:
  47.                         return N
  48.             N-=1
  49.         return 0

  50. def S(n):
  51.     global zs_half
  52.     x = 0
  53.     for i in zs_half:
  54.         for j in zs_half:
  55.             if i > j:
  56.                 x =x+M(i,j,n) #i大于j进行计算累加
  57.             else:
  58.                 break  #当i<=j的时候直接跳出J循环,进入下一个i循环。
  59.     return x




  60. if __name__ == '__main__':
  61.     n=10000
  62.     start = time.time()
  63.     zs_all=zhisu(n)
  64.     zs_half=zhisu(n//2)
  65.     zs = zhisu(int(n**0.5))
  66.     start1=time.time()
  67.     print('质数创建耗时:',start1-start,'秒','n =',n)
  68.     print(S(n))
  69.     print('完成计算耗时:',time.time()-start1,'秒','n =',n)
复制代码
想知道小甲鱼最近在做啥?请访问 -> ilovefishc.com
回复 支持 反对

使用道具 举报

发表于 2017-10-26 11:30:10 | 显示全部楼层
@jerryxjr1220
你试试 luajit
想知道小甲鱼最近在做啥?请访问 -> ilovefishc.com
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

小黑屋|手机版|Archiver|鱼C工作室 ( 粤ICP备18085999号-1 | 粤公网安备 44051102000585号)

GMT+8, 2024-3-29 14:22

Powered by Discuz! X3.4

© 2001-2023 Discuz! Team.

快速回复 返回顶部 返回列表