数学研发论坛

 找回密码
 欢迎注册
楼主: gxqcn

[擂台] 求形如 2^r*3^s*5^t 的整数序列

[复制链接]
发表于 2013-10-23 18:11:35 | 显示全部楼层
39#(我写的代码)只是用来介绍基本算法的,代码比较简单。我实际中用的代码比这速度更快,当然代码也更复杂些。
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2013-10-23 18:12:41 | 显示全部楼层
(*计算小于等于m的Pn的光滑数的个数的子函数*)
Clear["Global`*"];(*Clear all variables*)
prs=Table[Prime@n,{n,125}];(*产生前125个素数*)
fun[n0_,m0_]:=Module[{n=n0,m=m0,pr,npr,out},
                 pr=prs[[1;;n]];(*前n个素数*)
                 npr=Floor@Log[pr,m]+1;(*<=m的素数次幂+1*)
                 out=Apply[Times,npr](*乘法得到小于m的pn光滑数的个数*)
]

fun[3, 2^32]=9702对于2 3 5来说,小于2^32的个数是9702
我觉得用二分法,应该可以很快找到第x个这个数的大小

点评

代码有错误!  发表于 2013-10-23 18:21
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2013-10-23 18:26:33 | 显示全部楼层
mathematica 发表于 2013-10-23 18:12
(*计算小于等于m的Pn的光滑数的个数的子函数*)
Clear["Global`*"];(*Clear all variables*)
prs=Table;(* ...

可以利用筛法加加减减,然后较快地确定个数,但是个数多了也算得比较慢
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2013-10-24 10:41:21 | 显示全部楼层

RE: 如果我想再加一个bag怎么办?

chyanog 发表于 2013-10-23 17:48
当list很大的时候,AppendTo会越来越慢,不如用Internal`Bag{5.515315, 1434246505}{0.905052, 143424650 ...

我想再加一个bag,用来装第1个bag中相邻2项互质的对,写成 bag1=Bag[]行么?
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2013-10-24 12:50:36 | 显示全部楼层
liangbch 发表于 2013-10-23 17:10
这个题目的另一个说法是“求5-smooth数",百度知道(34#)给出的代码简单而有效,但有其局限性。用此代码求更 ...

(*计算小于等于m的Pn的光滑数的个数的子函数*)
Clear["Global`*"];(*Clear all variables*)
count=0; minnum=0; maxnum=2^64;
Do[num=
       (Prime@1)^n1
      *(Prime@2)^n2
      *(Prime@3)^n3
      *(Prime@4)^n4
  ;If[num<=maxnum,count=count+1]
  ,{n1,minnum,Log[Prime@1,maxnum]}
  ,{n2,minnum,Log[Prime@2,maxnum]}
  ,{n3,minnum,Log[Prime@3,maxnum]}
  ,{n4,minnum,Log[Prime@4,maxnum]}
];count

前四个素数构成的光滑数,个数是85349
理论上可以把所有的都计算出来
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2013-10-24 12:51:37 | 显示全部楼层
前31个素数构成的光滑数,计算代码如下:
(*计算小于等于m的Pn的光滑数的个数的子函数*)
Clear["Global`*"];(*Clear all variables*)
count=0; minnum=0; maxnum=2^64;
Do[num=
       (Prime@1)^n1
      *(Prime@2)^n2
      *(Prime@3)^n3
      *(Prime@4)^n4
      *(Prime@5)^n5
      *(Prime@6)^n6
      *(Prime@7)^n7
      *(Prime@8)^n8
      *(Prime@9)^n9
      *(Prime@10)^n10
      *(Prime@11)^n11
      *(Prime@12)^n12
      *(Prime@13)^n13
      *(Prime@14)^n14
      *(Prime@15)^n15
      *(Prime@16)^n16
      *(Prime@17)^n17
      *(Prime@18)^n18
      *(Prime@19)^n19
      *(Prime@20)^n20
      *(Prime@21)^n21
      *(Prime@22)^n22
      *(Prime@23)^n23
      *(Prime@24)^n24
      *(Prime@25)^n25
      *(Prime@26)^n26
      *(Prime@27)^n27
      *(Prime@28)^n28
      *(Prime@29)^n29
      *(Prime@30)^n30
      *(Prime@31)^n31
  ;If[num<=maxnum,count=count+1]
  ,{n1,minnum,Log[Prime@1,maxnum]}
  ,{n2,minnum,Log[Prime@2,maxnum]}
  ,{n3,minnum,Log[Prime@3,maxnum]}
  ,{n4,minnum,Log[Prime@4,maxnum]}
  ,{n5,minnum,Log[Prime@5,maxnum]}
  ,{n6,minnum,Log[Prime@6,maxnum]}
  ,{n7,minnum,Log[Prime@7,maxnum]}
  ,{n8,minnum,Log[Prime@8,maxnum]}
  ,{n9,minnum,Log[Prime@9,maxnum]}
  ,{n10,minnum,Log[Prime@10,maxnum]}
  ,{n11,minnum,Log[Prime@11,maxnum]}
  ,{n12,minnum,Log[Prime@12,maxnum]}
  ,{n13,minnum,Log[Prime@13,maxnum]}
  ,{n14,minnum,Log[Prime@14,maxnum]}
  ,{n15,minnum,Log[Prime@15,maxnum]}
  ,{n16,minnum,Log[Prime@16,maxnum]}
  ,{n17,minnum,Log[Prime@17,maxnum]}
  ,{n18,minnum,Log[Prime@18,maxnum]}
  ,{n19,minnum,Log[Prime@19,maxnum]}
  ,{n20,minnum,Log[Prime@20,maxnum]}
  ,{n21,minnum,Log[Prime@21,maxnum]}
  ,{n22,minnum,Log[Prime@22,maxnum]}
  ,{n23,minnum,Log[Prime@23,maxnum]}
  ,{n24,minnum,Log[Prime@24,maxnum]}
  ,{n25,minnum,Log[Prime@25,maxnum]}
  ,{n26,minnum,Log[Prime@26,maxnum]}
  ,{n27,minnum,Log[Prime@27,maxnum]}
  ,{n28,minnum,Log[Prime@28,maxnum]}
  ,{n29,minnum,Log[Prime@29,maxnum]}
  ,{n30,minnum,Log[Prime@30,maxnum]}
  ,{n31,minnum,Log[Prime@31,maxnum]}
];count

点评

运行时间可能不是一般的长,请耐心等待!  发表于 2013-10-24 12:52
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2013-10-24 12:57:49 | 显示全部楼层
(*计算小于等于m的Pn的光滑数的个数的子函数*)
Clear["Global`*"];(*Clear all variables*)
count=0; minnum=0; maxnum=2^64;
Do[num=
       (Prime@1)^n1
      *(Prime@2)^n2
      *(Prime@3)^n3
      *(Prime@4)^n4
      *(Prime@5)^n5
  ;If[num<=maxnum,count=count+1]
  ,{n1,minnum,Log[Prime@1,maxnum]}
  ,{n2,minnum,Log[Prime@2,maxnum/(Prime@1)^n1]}
  ,{n3,minnum,Log[Prime@3,maxnum/(Prime@2)^n2/(Prime@1)^n1]}
  ,{n4,minnum,Log[Prime@4,maxnum/(Prime@3)^n3/(Prime@2)^n2/(Prime@1)^n1]}
  ,{n5,minnum,Log[Prime@5,maxnum/(Prime@4)^n4/(Prime@3)^n3/(Prime@2)^n2/(Prime@1)^n1]}
];count

前五个素数,代码效率稍微高些
378555
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2013-10-24 12:59:54 | 显示全部楼层
前六个素数1381207
(*计算小于等于m的Pn的光滑数的个数的子函数*)
Clear["Global`*"];(*Clear all variables*)
count=0; minnum=0; maxnum=2^64;
Do[num=
       (Prime@1)^n1
      *(Prime@2)^n2
      *(Prime@3)^n3
      *(Prime@4)^n4
      *(Prime@5)^n5
      *(Prime@6)^n6
  ;If[num<=maxnum,count=count+1]
  ,{n1,minnum,Log[Prime@1,maxnum]}
  ,{n2,minnum,Log[Prime@2,maxnum/(Prime@1)^n1]}
  ,{n3,minnum,Log[Prime@3,maxnum/(Prime@2)^n2/(Prime@1)^n1]}
  ,{n4,minnum,Log[Prime@4,maxnum/(Prime@3)^n3/(Prime@2)^n2/(Prime@1)^n1]}
  ,{n5,minnum,Log[Prime@5,maxnum/(Prime@4)^n4/(Prime@3)^n3/(Prime@2)^n2/(Prime@1)^n1]}
  ,{n6,minnum,Log[Prime@6,maxnum/(Prime@5)^n5/(Prime@4)^n4/(Prime@3)^n3/(Prime@2)^n2/(Prime@1)^n1]}
];count

毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2013-10-24 13:03:33 | 显示全部楼层
hujunhua 发表于 2013-10-24 10:41
我想再加一个bag,用来装第1个bag中相邻2项互质的对,写成 bag1=Bag[]行么?

可以的
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2013-10-24 13:33:04 | 显示全部楼层
  1. (*我猜测的通项公式*)
  2. (*前n个素数,小于等于m的光滑数的个数*)
  3. Clear["Global`*"];(*Clear all variables*)
  4. fun[n0_,m0_]:=Module[{n=n0,m=m0,pn,pn2,pn3,out},
  5.                  pn=Table[Prime@k,{k,1,n}];(*前n个素数表*)
  6.                  pn2=Apply[Times,pn];(*前n个素数的乘积*)
  7.                  pn3=Apply[Times,Log[pn]];(*前n个素数的对数的乘积*)
  8.                  out=Ceiling[(Log[m*Sqrt[pn2]])^n/(n!*pn3)]
  9. ]
  10. maxnum=2^64;
  11. fun[1,maxnum](*结果精确命中*)
  12. fun[2,maxnum](*结果精确命中*)
  13. fun[3,maxnum]
复制代码
估计m越大,结果越精确!
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
您需要登录后才可以回帖 登录 | 欢迎注册

本版积分规则

小黑屋|手机版|数学研发网 ( 苏ICP备07505100号 )

GMT+8, 2019-3-26 08:15 , Processed in 0.047470 second(s), 15 queries .

Powered by Discuz! X3.4

© 2001-2017 Comsenz Inc.

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