QQ登录

只需一步,快速开始

登录 | 立即注册 | 找回密码

主题

帖子

荣誉

版主

Rank: 20Rank: 20Rank: 20Rank: 20

技术值
查看: 353|回复: 4

[活动专区] 小练习:20171113 哈沙德数

[复制链接]
最佳答案
495 
累计签到:715 天
连续签到:14 天
冬雪雪冬 发表于 2017-11-12 20:22:20 3534 | 显示全部楼层 |阅读模式

马上注册加入鱼C,享用更多服务吧^_^

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

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

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

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

什么是欧拉计划: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。
求所有小于1014的可右截强哈沙德素数之和。


楼层
跳转到指定楼层
最佳答案
54 
累计签到:351 天
连续签到:2 天
jerryxjr1220 发表于 2017-11-13 11:38:25 | 显示全部楼层
这题的难点其实是大数的质数检验,因为10**14的质数判断,用普通方法实在是太慢了,我参考了维基百科上的Miller Rabin质数检验法,可以非常快速的判断超大质数,而且准确率相当高。
https://rosettacode.org/wiki/Miller%E2%80%93Rabin_primality_test#Python

其次,就是利用广度优先搜索生成右截强哈沙德数
  1. import random

  2. N=10**14

  3. def is_prime(n):
  4.     """
  5.     Miller-Rabin primality test.
  6.     A return value of False means n is certainly not prime. A return value of
  7.     True means n is very likely a prime.
  8.     """
  9.     _mrpt_num_trials = 5 # number of bases to test
  10.     if n < 2: return False
  11.     if n == 2: return True
  12.     # ensure n is odd
  13.     if n % 2 == 0:
  14.         return False
  15.     # write n-1 as 2**s * d
  16.     # repeatedly try to divide n-1 by 2
  17.     s = 0
  18.     d = n-1
  19.     while True:
  20.         quotient, remainder = divmod(d, 2)
  21.         if remainder == 1:
  22.             break
  23.         s += 1
  24.         d = quotient
  25.     assert(2**s * d == n-1)

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

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

  38.     return True # no base tested showed n as composite

  39. def strong(n, s):
  40.     if n >= N: return False
  41.     if is_prime(n//s): return True
  42.     return False

  43. def harshads(n, s):
  44.     if n>=N or not n%s==0:
  45.         return
  46.     else:
  47.         yield (n,s)
  48.         v = 10*n
  49.         for i in range(10):
  50.             for (hh, ss) in harshads(v+i, s+i):
  51.                 yield (hh, ss)

  52. def primeExtensions(n):
  53.     v = 10*n
  54.     for d in [1,3,7,9]:
  55.         if v+d < N and is_prime(v+d):
  56.             yield v+d

  57. ans=0
  58. for i in range(1,10):
  59.     for h, s in harshads(i,i):
  60.         if strong(h,s):
  61.             for p in primeExtensions(h):
  62.                 ans += p
  63. print (ans)
复制代码

696067597313468
[Finished in 0.4s]

评分

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

查看全部评分

最佳答案
1 

尚未签到

gunjang 发表于 2017-11-14 22:19:51 | 显示全部楼层
本帖最后由 gunjang 于 2017-11-14 22:22 编辑

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

  3. st = time.time()
  4. limit = 10000000 #10^7
  5. primes = np.ones(limit, dtype=np.bool)
  6. for i in range(2, int(limit**0.5)+1):
  7.         if primes[i]: primes[i*i::i] = False
  8. primes[0:2] = False
  9. plist = np.nonzero(primes)[0]

  10. def isprime(n):
  11.         if n < limit:
  12.                 return primes[n]
  13.         else:
  14.                 sq = int(n**0.5)
  15.                 for p in plist:
  16.                         if p > sq:
  17.                                 return True
  18.                         if n%p == 0:
  19.                                 return False

  20. def generaterRTHarshadNum(digno):
  21.         hnlist = list(range(1, 10)) #digit number==1
  22.         r = 0
  23.         for dn in range(2, digno):
  24.                 nextHN = []
  25.                 for n in hnlist:
  26.                         d = n*10
  27.                         s = sum([int(x) for x in str(n)])
  28.                         for last in range(10):
  29.                                 if d%s==0:
  30.                                         nextHN.append(d)
  31.                                         if isprime(d//s):
  32.                                                 #d is strong Harshad number
  33.                                                 for step in [1,3,7,9]:
  34.                                                         if isprime(d*10+step):
  35.                                                                 r += d*10+step
  36.                                 d +=1
  37.                                 s +=1
  38.                 hnlist = nextHN
  39.         return r

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

  41. print(generaterRTHarshadNum(14)) #696067597313468 cost 6.6s
  42. print('cost %fs'%(time.time()-st))
复制代码

评分

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

查看全部评分

最佳答案
0 

尚未签到

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

点评

不,我坚决不同意楼主的看法!: 5.0
不,我坚决不同意楼主的看法!: 5
你确定10**14的数量级用穷举法能算得出来?  发表于 2017-11-16 12:10

发表回复

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

本版积分规则

关闭

小甲鱼强烈推荐 上一条 /1 下一条

    移动客户端下载(未启用)
    微信公众号

小黑屋|手机版|Archiver|鱼C工作室 ( 粤ICP备11014136号

Copyright 2018 鱼C论坛 版权所有 All Rights Reserved.

Powered by Discuz! X3.1 Copyright
© 2001-2018 Comsenz Inc.    All Rights Reserved.

小黑屋|手机版|Archiver|鱼C工作室 ( 粤公网安备 44051102000370号 | 粤ICP备11014136号

GMT+8, 2017-11-25 19:13

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