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

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

[复制链接]
发表于 2013-10-23 17:48:53 | 显示全部楼层
本帖最后由 chyanog 于 2013-10-23 17:54 编辑
hujunhua 发表于 2013-10-23 02:37
是由楼上链结中的程序改编的,效率当然也很高喽。

当list很大的时候,AppendTo会越来越慢,不如用Internal`Bag
  1. Module[{list, X, dp, checkb, pointer},
  2.    list = {};
  3.    X = {2, 3, 5};
  4.    checkb = {1, 1, 1};
  5.    pointer = {0, 0, 0};
  6.    Do[AppendTo[list, Min@checkb];
  7.     dp = 1 - Sign@(checkb - Last@list);
  8.     pointer += dp;
  9.     checkb += (list[[pointer]]*dp*X - checkb*dp)
  10.     , {2 10^4}];
  11.    list
  12.    ] // Hash // AbsoluteTiming
复制代码
{5.515315, 1434246505}
  1. \$ContextPath = Union@Append[\$ContextPath, "Internal`"];
  2. Module[{bag, X, dp, checkb, pointer},
  3.    bag = Bag[];
  4.    X = {2, 3, 5};
  5.    checkb = {1, 1, 1};
  6.    pointer = {0, 0, 0};
  7.    Do[
  8.     StuffBag[bag, Min@checkb];
  9.     dp = 1 - Sign@(checkb - BagPart[bag, -1]);
  10.     pointer += dp;
  11.     checkb += (BagPart[bag, pointer]*dp*X - checkb*dp),
  12.     {2 10^4}];
  13.    bag~BagPart~All
  14.    ] // Hash // AbsoluteTiming
复制代码
{0.905052, 1434246505}

点评

36#的Thread果然多余,去掉就更简洁了。内部包,果然强大,还没搞懂,有待研习。谢谢!  发表于 2013-10-23 18:35
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 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: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: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[]行么?

可以的
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
您需要登录后才可以回帖 登录 | 欢迎注册

本版积分规则

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

GMT+8, 2024-12-22 11:05 , Processed in 0.034792 second(s), 18 queries .

Powered by Discuz! X3.5

© 2001-2024 Discuz! Team.

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