鱼C论坛

 找回密码
 立即注册
查看: 2559|回复: 2

[技术交流] 小练习:20170911 巨大的阶乘

[复制链接]
发表于 2017-9-10 18:20:45 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能^_^

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

x
本帖最后由 冬雪雪冬 于 2017-9-20 19:55 编辑

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

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

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

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

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

这里已经有了500余题。


                               
登录/注册后可看大图


要求:

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

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

注重程序效率和创意。

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

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

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

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


                               
登录/注册后可看大图
288

题目:

1.jpg



2.jpg

想知道小甲鱼最近在做啥?请访问 -> ilovefishc.com
回复

使用道具 举报

发表于 2017-9-13 05:24:44 | 显示全部楼层
好吧,出差也不忘记解个题

这题一开始看题目觉得很难,其实仔细推敲起来还可以,下面说说我的解题思路。

首先,当然要把复杂的概念化简,我们先看看NF(3,10000)的前几项是什么?

  1. T0=1
  2. N(3,0)=T0
  3. Nfac(3,0)=N(3,0)!=1!
  4. Nf(3,0)=0

  5. T1=1
  6. N(3,1)=T0+T1*3**1
  7. Nfac(3,1)=N(3,1)!=(T0+3*T1)!=4!=4*3*2*1
  8. Nf(3,1)=1

  9. T2=1
  10. N(3,2)=T0+T1*3**1+T2*3**2
  11. Nfac(3,2)=N(3,2)!=(T0+T1*3+3*T2*3)!=13!
  12. Nf(3,2)=13//3+4//3=5

  13. T3=2
  14. Nfac(3,3)=(T0+T1*3+3*3*T2+3*3*3*T3)!=(13+54)!=(67)!
  15. Nf(3,3)=67//3+22//3+7//3=22+7+2=31=(T1+3*T2+3*3*T3)+(T1+3*T2+3*3*T3)%3+((T1+3*T2+3*3*T3)%3)%3
复制代码


这题的难点其实是在算n阶乘中质数p的个数。从上面的展开项可以看到,n阶乘中p的个数可以用n反复整除p来计算。

于是,求解Nf(3,10000)就可以写成:
  1. S=290797
  2. T=[S%3]
  3. for i in range(10000):
  4.         S=(S*S)%50515093
  5.         T.append(S%3)
  6. q=0
  7. for i in range(1,10001):
  8.         q += T[i]*3**(i-1)
  9. s = q % 3**20
  10. while q>=3:
  11.         q //= 3
  12.         s = (s+q) % 3**20
  13. print(s)  # 624955285
复制代码


答案也是正确的,但是如果要求Nf(61,10**7)的话,上述的运算效率就比较慢了,不过我们可以借鉴上面这个的思路。

既然
Nf(3,3)=67//3+22//3+7//3=22+7+2=31=(T1+3*T2+3*3*T3)+(T1+3*T2+3*3*T3)%3+((T1+3*T2+3*3*T3)%3)%3


那么我们是不是可以从内层开始逐渐往外层推呢?然后每层都可以取模运算,减小运算量。

于是Nf(61,10**7)就可以写成:
  1. S=290797
  2. T=[S%61]
  3. for i in range(10**7):
  4.         S=(S*S)%50515093
  5.         T.append(S%61)
  6. q=0
  7. s=0
  8. for i in range(10**7,0,-1):
  9.         q = (q*61+T[i])%61**10
  10.         s = (s+q)%61**10
  11. print(s)  # 605857431263981935
复制代码


差不多9s可以出答案。

评分

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

查看全部评分

想知道小甲鱼最近在做啥?请访问 -> ilovefishc.com
回复 支持 反对

使用道具 举报

发表于 2017-9-16 23:47:21 | 显示全部楼层
本帖最后由 gunjang 于 2017-9-17 08:38 编辑

终于赶上了。。。605857431263981935
  1. import numpy as np
  2. import time
  3. #S0 = 290797
  4. #Sn+1 = Sn^2 mod 50515093
  5. #Tn = Sn mod p

  6. '''
  7. N(p,q) = t0 + t1*p^1+t2*p^2+...+tn*p^n
  8. NF(p,q) = npq//p + npq//p^2 + ...npq//p^n
  9. because t0...tn < p
  10. NF(p,q) div p   = t1 + t2*p^1 + t3*p^2 + ...+tn*p^(n-1)
  11. NF(p,q) div p^2 = t2 + t3*p^1 + ...+tn*p^(n-2)
  12. NF(p,q) div p^n = tn
  13. sum of NF(p,q) div p^(1..n) = p^0*(t1+t2+..+tn) + p^1*(t2+t3+..tn) + p^(n-1)*(tn)

  14. let item(x) = p^x*(t(x+1)+t(x+2)+..tn)
  15. if x >= k, item(x) mod p^k == 0
  16. so NF(p,q) mod p^k == p^0*(t1+t2+..+tn) + p^1*(t2+t3+..tn)+ ...p^(k-1)(tk+tk+1+..tn)
  17. '''

  18. def getNFpqModpK(p, q, k):
  19.         tn = np.zeros((k+1), dtype=np.int16)
  20.         s = 290797
  21.         for i in range(1, k+1):
  22.                 s = (s*s) % 50515093
  23.                 tn[i] = s % p

  24.         sumtk1ton = 0 #sum of T(k+1)..T(n)
  25.         for i in range(k+1, q+1):
  26.                 s = (s*s) % 50515093
  27.                 sumtk1ton += s%p
  28.        
  29.         r = 0
  30.         ppi = 1
  31.         pn = np.zeros((k), dtype=np.int64) #p^0...p^k-1
  32.         for i in range(k):
  33.                 pn[i] = (sum(tn[i+1:]) + sumtk1ton) % p**(k-i)
  34.                 #r += (sum(tn[i+1:]) + sumtk1ton) * ppi # *p**i  #may RuntimeWarning: overfow encountered in long_scalars
  35.                 r += pn[i] * ppi
  36.                 ppi *= p

  37.         return r % ppi  #(p**k)

  38. st = time.time()
  39. print(getNFpqModpK(3, 10000, 20)) #624955285
  40. print(getNFpqModpK(61, 10**7, 10)) #605857431263981935
  41. print('cost {0} s'.format(round(time.time()-st, 2))) #6.31s
复制代码

评分

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

查看全部评分

想知道小甲鱼最近在做啥?请访问 -> ilovefishc.com
回复 支持 反对

使用道具 举报

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

本版积分规则

小黑屋|手机版|Archiver|鱼C工作室 ( 粤ICP备18085999号-1 | 粤公网安备 44051102000585号)

GMT+8, 2024-3-29 15:07

Powered by Discuz! X3.4

© 2001-2023 Discuz! Team.

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