小练习：20171113 哈沙德数

冬雪雪冬

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

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

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

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

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

这里已经有了500余题。

                               
登录/注册后可看大图

要求：

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

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

注重程序效率和创意。

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

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

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

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

                               
登录/注册后可看大图

387

能够被其各位数字和整除的数被称为哈沙德数或奈文数。
201是一个哈沙德数，因为它能被3整除（其各位数字和是3）。
如果我们截去201的最后一位数字，我们得到20，同样是一个哈沙德数。
如果我们截去20的最后一位数字，我们得到2，仍然是一个哈沙德数。
如果一个哈沙德数不断截去最后一位数字的结果始终是哈沙德数，我们称之为可右截哈沙德数。

此外：
201/3=67是一个素数。
如果一个哈沙德数被其各位数字除的结果是一个素数，我们称之为强哈沙德数。

现在我们取素数2011。
如果我们截去2011的最后一位数字，我们得到201，一个可右截的强哈沙德数。
我们称这样的素数为可右截强哈沙德素数。

已知所有小于10000的可右截强哈沙德素数之和为90619。

求所有小于10¹⁴的可右截强哈沙德素数之和。

jerryxjr1220

这题的难点其实是大数的质数检验，因为10**14的质数判断，用普通方法实在是太慢了，我参考了维基百科上的Miller Rabin质数检验法，可以非常快速的判断超大质数，而且准确率相当高。

https://rosettacode.org/wiki/Miller%E2%80%93Rabin_primality_test#Python

其次，就是利用广度优先搜索生成右截强哈沙德数

1. import random
   
2. 
   
3. N=10**14
   
4. 
   
5. def is_prime(n):
   
6.     """
   
7.     Miller-Rabin primality test.
   
8.     A return value of False means n is certainly not prime. A return value of
   
9.     True means n is very likely a prime.
   
10.     """
    
11.     _mrpt_num_trials = 5 # number of bases to test
    
12.     if n < 2: return False
    
13.     if n == 2: return True
    
14.     # ensure n is odd
    
15.     if n % 2 == 0:
    
16.         return False
    
17.     # write n-1 as 2**s * d
    
18.     # repeatedly try to divide n-1 by 2
    
19.     s = 0
    
20.     d = n-1
    
21.     while True:
    
22.         quotient, remainder = divmod(d, 2)
    
23.         if remainder == 1:
    
24.             break
    
25.         s += 1
    
26.         d = quotient
    
27.     assert(2**s * d == n-1)
    
28. 
    
29.     # test the base a to see whether it is a witness for the compositeness of n
    
30.     def try_composite(a):
    
31.         if pow(a, d, n) == 1:
    
32.             return False
    
33.         for i in range(s):
    
34.             if pow(a, 2**i * d, n) == n-1:
    
35.                 return False
    
36.         return True # n is definitely composite
    
37. 
    
38.     for i in range(_mrpt_num_trials):
    
39.         a = random.randrange(2, n)
    
40.         if try_composite(a):
    
41.             return False
    
42. 
    
43.     return True # no base tested showed n as composite
    
44. 
    
45. def strong(n, s):
    
46.     if n >= N: return False
    
47.     if is_prime(n//s): return True
    
48.     return False
    
49. 
    
50. def harshads(n, s):
    
51.     if n>=N or not n%s==0:
    
52.         return
    
53.     else:
    
54.         yield (n,s)
    
55.         v = 10*n
    
56.         for i in range(10):
    
57.             for (hh, ss) in harshads(v+i, s+i):
    
58.                 yield (hh, ss)
    
59. 
    
60. def primeExtensions(n):
    
61.     v = 10*n
    
62.     for d in [1,3,7,9]:
    
63.         if v+d < N and is_prime(v+d):
    
64.             yield v+d
    
65. 
    
66. ans=0
    
67. for i in range(1,10):
    
68.     for h, s in harshads(i,i):
    
69.         if strong(h,s):
    
70.             for p in primeExtensions(h):
    
71.                 ans += p
    
72. print (ans)

复制代码

696067597313468
[Finished in 0.4s]

gunjang

696067597313468
没优化大概6s多

1. import time
   
2. import numpy as np
   
3. 
   
4. st = time.time()
   
5. limit = 10000000 #10^7
   
6. primes = np.ones(limit, dtype=np.bool)
   
7. for i in range(2, int(limit**0.5)+1):
   
8.         if primes[i]: primes[i*i::i] = False
   
9. primes[0:2] = False
   
10. plist = np.nonzero(primes)[0]
    
11. 
    
12. def isprime(n):
    
13.         if n < limit:
    
14.                 return primes[n]
    
15.         else:
    
16.                 sq = int(n**0.5)
    
17.                 for p in plist:
    
18.                         if p > sq:
    
19.                                 return True
    
20.                         if n%p == 0:
    
21.                                 return False
    
22. 
    
23. def generaterRTHarshadNum(digno):
    
24.         hnlist = list(range(1, 10)) #digit number==1
    
25.         r = 0
    
26.         for dn in range(2, digno):
    
27.                 nextHN = []
    
28.                 for n in hnlist:
    
29.                         d = n*10
    
30.                         s = sum([int(x) for x in str(n)])
    
31.                         for last in range(10):
    
32.                                 if d%s==0:
    
33.                                         nextHN.append(d)
    
34.                                         if isprime(d//s):
    
35.                                                 #d is strong Harshad number
    
36.                                                 for step in [1,3,7,9]:
    
37.                                                         if isprime(d*10+step):
    
38.                                                                 r += d*10+step
    
39.                                 d +=1
    
40.                                 s +=1
    
41.                 hnlist = nextHN
    
42.         return r
    
43. 
    
44. #print(generaterRTHarshadNum(4)) # less than 10000 is 90619.
    
45. 
    
46. print(generaterRTHarshadNum(14)) #696067597313468 cost 6.6s
    
47. print('cost %fs'%(time.time()-st))

复制代码

mrchenqi1981

import math

def isHarshadNumber(num):
    s = 0
    n = num
    while (num != 0):
        s += num % 10
        num = int(num / 10)

    if (s == 0):
        return False
   
    return (n % s == 0)

def isRightTruncatedHarshadNumber(num):
    while (num != 0):
        if (not isHarshadNumber(num)):
            return False

        num = int(num / 10)

    return True

def isPrime(num):
    if (num == 1):
        return False
   
    for i in range(2, int(math.sqrt(num)) + 1):
        if (num % i == 0):
            return False

    return True

def isStrongHarshadNumber(num):
    s = 0
    n = num
    while (num != 0):
        s += num % 10
        num = int(num / 10)

    if (s == 0):
        return False
   
    if (isHarshadNumber(n) and isPrime(int(n / s))):
        return True

    return False

def isRightTruncatedStrongHarshadPrime(num):
    return (isPrime(num) and isRightTruncatedHarshadNumber(int(num / 10)) and isStrongHarshadNumber(int(num / 10)))

def isRightTruncatedHarshadPrime(num):
    return (isRightTruncatedHarshadNumber(num) and isPrime(num))

def mysum(n):
    s = 0
    for i in range (1, n):
        if (isRightTruncatedStrongHarshadPrime(i)):
            print (i)
            s += i

    return s