论坛 › Python交流 › 小练习：20171113 哈沙德数

冬雪雪冬 发表于 2017-11-12 20:22:20

小练习：20171113 哈沙德数

本帖最后由 冬雪雪冬 于 2017-11-19 14:33 编辑

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

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

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

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

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

这里已经有了500余题。

http://bbs.fishc.com/static/image/hrline/1.gif

要求：

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

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

注重程序效率和创意。

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

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

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

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

http://bbs.fishc.com/static/image/hrline/1.gif387

能够被其各位数字和整除的数被称为哈沙德数或奈文数。
201是一个哈沙德数，因为它能被3整除（其各位数字和是3）。
如果我们截去201的最后一位数字，我们得到20，同样是一个哈沙德数。
如果我们截去20的最后一位数字，我们得到2，仍然是一个哈沙德数。
如果一个哈沙德数不断截去最后一位数字的结果始终是哈沙德数，我们称之为可右截哈沙德数。此外：
201/3=67是一个素数。
如果一个哈沙德数被其各位数字除的结果是一个素数，我们称之为强哈沙德数。现在我们取素数2011。
如果我们截去2011的最后一位数字，我们得到201，一个可右截的强哈沙德数。
我们称这样的素数为可右截强哈沙德素数。已知所有小于10000的可右截强哈沙德素数之和为90619。求所有小于1014的可右截强哈沙德素数之和。

jerryxjr1220 发表于 2017-11-13 11:38:25

这题的难点其实是大数的质数检验，因为10**14的质数判断，用普通方法实在是太慢了，我参考了维基百科上的Miller Rabin质数检验法，可以非常快速的判断超大质数，而且准确率相当高。
https://rosettacode.org/wiki/Miller%E2%80%93Rabin_primality_test#Python
其次，就是利用广度优先搜索生成右截强哈沙德数
import random

N=10**14

def is_prime(n):
    """
    Miller-Rabin primality test.
    A return value of False means n is certainly not prime. A return value of
    True means n is very likely a prime.
    """
    _mrpt_num_trials = 5 # number of bases to test
    if n < 2: return False
    if n == 2: return True
    # ensure n is odd
    if n % 2 == 0:
      return False
    # write n-1 as 2**s * d
    # repeatedly try to divide n-1 by 2
    s = 0
    d = n-1
    while True:
      quotient, remainder = divmod(d, 2)
      if remainder == 1:
            break
      s += 1
      d = quotient
    assert(2**s * d == n-1)

    # test the base a to see whether it is a witness for the compositeness of n
    def try_composite(a):
      if pow(a, d, n) == 1:
            return False
      for i in range(s):
            if pow(a, 2**i * d, n) == n-1:
                return False
      return True # n is definitely composite

    for i in range(_mrpt_num_trials):
      a = random.randrange(2, n)
      if try_composite(a):
            return False

    return True # no base tested showed n as composite

def strong(n, s):
    if n >= N: return False
    if is_prime(n//s): return True
    return False

def harshads(n, s):
    if n>=N or not n%s==0:
      return
    else:
      yield (n,s)
      v = 10*n
      for i in range(10):
            for (hh, ss) in harshads(v+i, s+i):
                yield (hh, ss)

def primeExtensions(n):
    v = 10*n
    for d in :
      if v+d < N and is_prime(v+d):
            yield v+d

ans=0
for i in range(1,10):
    for h, s in harshads(i,i):
      if strong(h,s):
            for p in primeExtensions(h):
                ans += p
print (ans)
696067597313468

gunjang 发表于 2017-11-14 22:19:51

本帖最后由 gunjang 于 2017-11-14 22:22 编辑

696067597313468
没优化大概6s多
import time
import numpy as np

st = time.time()
limit = 10000000 #10^7
primes = np.ones(limit, dtype=np.bool)
for i in range(2, int(limit**0.5)+1):
        if primes: primes = False
primes = False
plist = np.nonzero(primes)

def isprime(n):
        if n < limit:
                return primes
        else:
                sq = int(n**0.5)
                for p in plist:
                        if p > sq:
                                return True
                        if n%p == 0:
                                return False

def generaterRTHarshadNum(digno):
        hnlist = list(range(1, 10)) #digit number==1
        r = 0
        for dn in range(2, digno):
                nextHN = []
                for n in hnlist:
                        d = n*10
                        s = sum()
                        for last in range(10):
                                if d%s==0:
                                        nextHN.append(d)
                                        if isprime(d//s):
                                                #d is strong Harshad number
                                                for step in :
                                                        if isprime(d*10+step):
                                                                r += d*10+step
                                d +=1
                                s +=1
                hnlist = nextHN
        return r

#print(generaterRTHarshadNum(4)) # less than 10000 is 90619.

print(generaterRTHarshadNum(14)) #696067597313468 cost 6.6s
print('cost %fs'%(time.time()-st))

mrchenqi1981 发表于 2017-11-15 16:36:07

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

查看完整版本: 小练习：20171113 哈沙德数