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

liangbch

39#(我写的代码）只是用来介绍基本算法的，代码比较简单。我实际中用的代码比这速度更快，当然代码也更复杂些。

mathematica

(*计算小于等于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个这个数的大小

mathematica

mathematica 发表于 2013-10-23 18:12
(*计算小于等于m的Pn的光滑数的个数的子函数*)
Clear["Global`*"];(*Clear all variables*)
prs=Table;(* ...

可以利用筛法加加减减,然后较快地确定个数,但是个数多了也算得比较慢

hujunhua

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

我想再加一个bag,用来装第1个bag中相邻2项互质的对，写成 bag1=Bag[]行么？

mathematica

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
理论上可以把所有的都计算出来

mathematica

前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

mathematica

(*计算小于等于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

mathematica

前六个素数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

chyanog

hujunhua 发表于 2013-10-24 10:41
我想再加一个bag,用来装第1个bag中相邻2项互质的对，写成 bag1=Bag[]行么？

可以的

mathematica

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越大,结果越精确!