小练习：20160711 找出为连续数字产生最多质数的二次公式

冬雪雪冬

从现在开始我们要开展一批欧拉计划的习题练习。
其实在我们论坛中已有欧拉计划的板块，可能有些鱼油还没注意到。
什么是欧拉计划：http://bbs.fishc.com/thread-60405-1-1.html
我们欧拉板块现已给出了81道题，这批练习将从欧拉计划中选题。其实用python语言完成有很多的优势，可以更简洁更方便的实现。
如果大家有兴趣也可浏览欧拉的英文网站：https://projecteuler.net/archives
这里已经有了500余题。

                               
登录/注册后可看大图

好了言归正传，我们开始做小练习。




题目要求：
以python语言完成，如果是python2请注明。
程序以代码文字格式发帖。
题目比较简单，注重程序效率和创意。
答题在一周内完成，即7.18 10：00之前，其后将公开大家的答案，并评比成绩。
另程序和答案可以在网上搜到，希望大家独立完成。

----除列出程序外，请给出输出的结果。----


题目：

欧拉曾发表过一个著名的二次公式：

                               
登录/注册后可看大图

这个公式对于0到39的连续数字能够产生 40 个质数。但是当 n = 40 时，

                               
登录/注册后可看大图

能够被 41 整除。当 n = 41 时,

                               
登录/注册后可看大图

显然也能被 41 整除。

利用计算机，人们发现了一个惊人的公式：

                               
登录/注册后可看大图

。这个公式对于 n = 0 到 79 能够产生 80 个质数。这个公式的系数，−79 和 1601 的乘积是 −126479。

考虑如下形式的二次公式：

                               
登录/注册后可看大图

其中 |n| 表示 n 的绝对值。

例如， |11| = 11， |−4| = 4

对于能够为从 0 开始的连续的 n 产生最多数量的质数的二次公式，找出该公式的系数a, b和它们的乘积。


我试了一下，大约用时6s多，看看大家有什么巧妙的方法。


奖励：
对所有完成程序并得出正确答案的将给予加分奖励，优秀的将额外加分。
在完成一批题目后，将根据每期的完成情况总评评出最佳，会有神秘大奖。

冬雪雪冬

1. import math
   
2. def iscomposit(n):
   
3.     for i in range(2, int(math.sqrt(n)) + 1):
   
4.         if n % i == 0:
   
5.             return True
   
6.     return False
   
7. 
   
8. maxprime = 0
   
9. for a in range(-999, 1000):
   
10.     for b in range(-999, 1000):
    
11.         n = 0
    
12.         while True:
    
13.             x = n * n  + a * n + b
    
14.             if x < 2 or iscomposit(x):
    
15.                 break
    
16.             n += 1
    
17.         if maxprime < n:
    
18.             maxprime = n
    
19.             maxa = a
    
20.             maxb = b
    
21. print('a= %d, b = %d, a * b = %d'%(maxa, maxb, maxa *maxb))

复制代码

wangzhenas

考虑 n=0时多项式结果为质数可知c是1000以内的质数
考虑 n=1时，可得1+b+c为质数，因此b>-c,而且b应为奇数（排除质数为2的情况，因为当n=2时多项式因为正数）

1. def primes_sieve(limit):
   
2.     limitn,not_prime,primes = limit+1,set(),list()
   
3.     for i in range(2, limitn):
   
4.         if i in not_prime:
   
5.             continue
   
6.         for f in range(i*2, limitn, i):
   
7.             not_prime.add(f)
   
8.         primes.append(i)
   
9.     return primes
   
10.    
    
11. 
    
12. import time
    
13. 
    
14. t,maxi,prime_l = time.time(),0,list(reversed(primes_sieve(1000)))   
    
15. prime_s = set(primes_sieve(10000))
    
16. 
    
17. for c in prime_l:
    
18.     for b in range(-c+2,1000,2):
    
19.         for i in range(1,100):
    
20.             if (i**2 + b*i + c) not in prime_s:
    
21.                 if i > maxi:
    
22.                     maxi,result = i,b*c
    
23.                 break
    
24.             
    
25. 
    
26. print('time:',time.time()-t,'result:',result)

复制代码

在我的机器上（i5-3210M）：
time: 0.23101305961608887 result: -59231

mather

1. import time
   
2. on = time.time()
   
3. 
   
4. 
   
5. from math import sqrt
   
6. def is_prime(num):
   
7.     if num<=1:
   
8.         return False
   
9.     t = 2
   
10.     res = True
    
11.     while t <= sqrt(num):
    
12.         if num%t==0:
    
13.             res = False
    
14.             break
    
15.         else:
    
16.             t += 1
    
17.     return res
    
18. 
    
19. max_primes = 41#要利用上已知的条件
    
20. res_a_b = 1,41
    
21. 
    
22. for a in range(-999,1000,2):#如果要1+a+b是质数，b是质数的前提下，a必然是奇数
    
23.     for b in range(2,1000):#b只能从2开始取,且b只能是质数
    
24.         if not is_prime(b):#如果b不是质数下一个
    
25.             continue
    
26.         elif b<max_primes:#如果比最多的个数值小下一个
    
27.             continue
    
28.         elif a<0 and not is_prime((max_primes-1+a)*(max_primes-1)+b):#在a是负数的时候很容易产生负数值，所以过滤一下
    
29.             continue
    
30.         else:
    
31.             for j in range(b):#最多从0到b-1，b个质数
    
32.                 if not is_prime((j+a)*j+b):#验证j是不是质数，如果是下一个
    
33.                     break
    
34.             if j > max_primes:#找到更大的，记录一下
    
35.                 max_primes = j
    
36.                 res_a_b = a,b
    
37. print("最多个数的系数组合是a=%d,b=%d，乘积是%d" % res_a_b,res_a_b[0]*res_a_b[1])
    
38. 
    
39. print(time.time()-on)

复制代码

结果：最多个数的系数组合是a=-61,b=971，乘积是 -59231
运行时间：4.0335609912872314

bacon6581

1. from time import time
   
2. start=time()
   
3. 
   
4. def zhishu(n):
   
5.     if n==2:
   
6.         return True
   
7.     if n**0.5==int(n**0.5):
   
8.         return False
   
9.     m=int(n**0.5)
   
10.     for i in range(2,m+1):
    
11.         if n%i==0:
    
12.             return False
    
13.     return True
    
14. 
    
15. a_result=0
    
16. b_result=0
    
17. x_max=0
    
18. 
    
19. a=-1001
    
20. while a<=1000:
    
21.     a+=1
    
22.     b=-1001
    
23. 
    
24.     while b<=1000:
    
25.         b+=1
    
26. 
    
27.         x=-1
    
28.         while 1:
    
29.             x+=1
    
30.             if x*x + a*x +b<=2:
    
31.                 break
    
32.             if zhishu(x*x + a*x +b)==False:
    
33.                 break
    
34. 
    
35.         if x>x_max:
    
36.             x_max=x
    
37.             a_result=a
    
38.             b_result=b
    
39. print('a*b is: ',a_result*b_result)
    
40. print('a is: ',a_result,'b is: ',b_result)
    
41. print('Used time: ',time()-start)

复制代码

huomqh

1. def isPrime(x):
   
2.     i=2
   
3.     if x<=1:
   
4.         return 0
   
5.    
   
6.     while i**2<=x:
   
7.         if x%i==0:
   
8.             return 0
   
9.         i+=1
   
10.     return 1
    
11. 
    
12. def f(x,y,z):
    
13.     return x*x+y*x+z
    
14. 
    
15. def g(x,y):
    
16.     if x<0:
    
17.         return x*x-4*y
    
18.     else:
    
19.         return -1
    
20.    
    
21. import time
    
22. tt=time.time()
    
23. b=3
    
24. count=0
    
25. counta=0
    
26. countb=0
    
27. while b<=999:
    
28.     if isPrime(b)==0:
    
29.         b+=2
    
30.         continue
    
31.     a=999
    
32.     while (a>=-999) and (g(a,b)<0):
    
33.         i=1
    
34.         while isPrime(f(i,a,b))==1:
    
35.             i+=1
    
36.         if i-1>count:
    
37.             count,counta,countb=i-1,a,b
    
38.         a-=2
    
39.     b+=2
    
40. 
    
41. print('a=%d,b=%d,a*b=%d,n=%d'%(counta,countb,counta*countb,count))
    
42. print('用时：%.4f s'%(time.time() - tt))
    

复制代码

结果：

1. a=-61,b=971,a*b=-59231,n=70
   
2. 用时：1.3211 s

复制代码

Spicebush

1. #答案为a = -61, b = 971, a * b = -59231
   
2. import math
   
3. def zs(n):
   
4.     if n>0:
   
5.         half=math.sqrt(n)
   
6.         for i in range(2,int(math.sqrt(n))+1):
   
7.             if n%i==0:
   
8.                 return False
   
9.                 break
   
10.         return True
    
11.     else:
    
12.         return False
    
13. 
    
14. js_l = [x for x in range(-999, 1000) if x%2 != 0]
    
15. zs_l = [x for x in range(-999, 1000) if zs(abs(x))]
    
16. 
    
17. z = 1
    
18. for a in js_l:
    
19.     for b in zs_l:
    
20.         n = 1
    
21.         while zs(n*n + a*n + b):
    
22.             n += 1
    
23.         if n > z:
    
24.             z, A, B = n, a, b
    
25.             
    
26. print('a =', A, '\nb =', B, '\na * b =', A*B)

复制代码

shichaoufo

a为-61,b为971，乘积为-59231

1. import math
   
2. import time
   
3. 
   
4. st = time.time()
   
5. 
   
6. jishu = []
   
7. fujishu = []
   
8. zhishua = []
   
9. zuida = 40
   
10. zuidaa = 1
    
11. zuidab = 41
    
12. 
    
13. def zspd(n):
    
14.     if n <= 1:
    
15.         return False
    
16.     y = int(math.sqrt(n))
    
17.     for i in range(2,y+1):
    
18.         if n % i == 0:
    
19.             return False
    
20.     return True
    
21. 
    
22. for i in range(999,1,-1):
    
23.     if zspd(i):
    
24.         zhishua.append(i)
    
25. 
    
26. for i in range(1,1000):
    
27.     if i % 2 != 0:
    
28.         jishu.append(i)
    
29.         
    
30. for i in jishu:
    
31.     fujishu.append(-i)
    
32. 
    
33. fanweia = fujishu + jishu
    
34. 
    
35. for a in fanweia:
    
36.     for b in zhishua:
    
37.         if zspd(zuida*zuida + a*zuida +b):
    
38.             n = 1
    
39.             x = a+b+1
    
40.             while zspd(x):
    
41.                 x += (2 * n + 1 + a)
    
42.                 n += 1
    
43.             if n>zuida:
    
44.                 zuida = n
    
45.                 zuidaa = a
    
46.                 zuidab = b
    
47.         else:
    
48.             next
    
49. 
    
50. cj = zuidaa * zuidab                  
    
51. print("a为%d,b为%d" % (zuidaa,zuidab))
    
52. print("乘积为%d" % cj)
    
53. leng = time.time() - st
    
54. print("用时%f" % leng)

复制代码

无符号整形

应该是20160711 吧？？？

冬雪雪冬

无符号整形 发表于 2016-7-15 17:36
应该是20160711 吧？？？

晕！犯了个低级错误。