小练习：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多，看看大家有什么巧妙的方法。


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

冬雪雪冬

import math
def iscomposit(n):
    for i in range(2, int(math.sqrt(n)) + 1):
        if n % i == 0:
            return True
    return False

maxprime = 0
for a in range(-999, 1000):
    for b in range(-999, 1000):
        n = 0
        while True:
            x = n * n  + a * n + b
            if x < 2 or iscomposit(x):
                break
            n += 1
        if maxprime < n:
            maxprime = n
            maxa = a
            maxb = b
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时多项式因为正数）

def primes_sieve(limit):
    limitn,not_prime,primes = limit+1,set(),list()
    for i in range(2, limitn):
        if i in not_prime:
            continue
        for f in range(i*2, limitn, i):
            not_prime.add(f)
        primes.append(i)
    return primes
    
 
import time

t,maxi,prime_l = time.time(),0,list(reversed(primes_sieve(1000)))   
prime_s = set(primes_sieve(10000))

for c in prime_l:
    for b in range(-c+2,1000,2):
        for i in range(1,100):
            if (i**2 + b*i + c) not in prime_s:
                if i > maxi:
                    maxi,result = i,b*c
                break
            

print('time:',time.time()-t,'result:',result)

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

mather

import time
on = time.time()


from math import sqrt
def is_prime(num):
    if num<=1:
        return False
    t = 2
    res = True
    while t <= sqrt(num):
        if num%t==0:
            res = False
            break
        else:
            t += 1
    return res

max_primes = 41#要利用上已知的条件
res_a_b = 1,41

for a in range(-999,1000,2):#如果要1+a+b是质数，b是质数的前提下，a必然是奇数
    for b in range(2,1000):#b只能从2开始取,且b只能是质数
        if not is_prime(b):#如果b不是质数下一个
            continue
        elif b<max_primes:#如果比最多的个数值小下一个
            continue
        elif a<0 and not is_prime((max_primes-1+a)*(max_primes-1)+b):#在a是负数的时候很容易产生负数值，所以过滤一下
            continue
        else:
            for j in range(b):#最多从0到b-1，b个质数
                if not is_prime((j+a)*j+b):#验证j是不是质数，如果是下一个
                    break
            if j > max_primes:#找到更大的，记录一下
                max_primes = j
                res_a_b = a,b
print("最多个数的系数组合是a=%d,b=%d，乘积是%d" % res_a_b,res_a_b[0]*res_a_b[1])

print(time.time()-on)

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

bacon6581

from time import time
start=time()

def zhishu(n):
    if n==2:
        return True
    if n**0.5==int(n**0.5):
        return False
    m=int(n**0.5)
    for i in range(2,m+1):
        if n%i==0:
            return False
    return True

a_result=0
b_result=0
x_max=0

a=-1001
while a<=1000:
    a+=1
    b=-1001

    while b<=1000:
        b+=1

        x=-1
        while 1:
            x+=1
            if x*x + a*x +b<=2:
                break
            if zhishu(x*x + a*x +b)==False:
                break

        if x>x_max:
            x_max=x
            a_result=a
            b_result=b
print('a*b is: ',a_result*b_result)
print('a is: ',a_result,'b is: ',b_result)
print('Used time: ',time()-start)

huomqh

def isPrime(x):
    i=2
    if x<=1:
        return 0
    
    while i**2<=x:
        if x%i==0:
            return 0
        i+=1
    return 1

def f(x,y,z):
    return x*x+y*x+z

def g(x,y):
    if x<0:
        return x*x-4*y
    else:
        return -1
    
import time
tt=time.time()
b=3
count=0
counta=0
countb=0
while b<=999:
    if isPrime(b)==0:
        b+=2
        continue
    a=999
    while (a>=-999) and (g(a,b)<0):
        i=1
        while isPrime(f(i,a,b))==1:
            i+=1
        if i-1>count:
            count,counta,countb=i-1,a,b
        a-=2
    b+=2

print('a=%d,b=%d,a*b=%d,n=%d'%(counta,countb,counta*countb,count))
print('用时：%.4f s'%(time.time() - tt))

结果：

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

Spicebush

#答案为a = -61, b = 971, a * b = -59231
import math
def zs(n):
    if n>0:
        half=math.sqrt(n)
        for i in range(2,int(math.sqrt(n))+1):
            if n%i==0:
                return False
                break
        return True
    else:
        return False

js_l = [x for x in range(-999, 1000) if x%2 != 0]
zs_l = [x for x in range(-999, 1000) if zs(abs(x))]

z = 1
for a in js_l:
    for b in zs_l:
        n = 1
        while zs(n*n + a*n + b):
            n += 1
        if n > z:
            z, A, B = n, a, b
            
print('a =', A, '\nb =', B, '\na * b =', A*B)

shichaoufo

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

import math
import time

st = time.time()

jishu = []
fujishu = []
zhishua = []
zuida = 40
zuidaa = 1
zuidab = 41

def zspd(n):
    if n <= 1:
        return False
    y = int(math.sqrt(n))
    for i in range(2,y+1):
        if n % i == 0:
            return False
    return True

for i in range(999,1,-1):
    if zspd(i):
        zhishua.append(i)

for i in range(1,1000):
    if i % 2 != 0:
        jishu.append(i)
        
for i in jishu:
    fujishu.append(-i)

fanweia = fujishu + jishu

for a in fanweia:
    for b in zhishua:
        if zspd(zuida*zuida + a*zuida +b):
            n = 1
            x = a+b+1
            while zspd(x):
                x += (2 * n + 1 + a)
                n += 1
            if n>zuida:
                zuida = n
                zuidaa = a
                zuidab = b
        else:
            next

cj = zuidaa * zuidab                   
print("a为%d,b为%d" % (zuidaa,zuidab))
print("乘积为%d" % cj)
leng = time.time() - st
print("用时%f" % leng)

无符号整形

应该是20160711 吧？？？

冬雪雪冬

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

晕！犯了个低级错误。