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

chyanog

mathematica 发表于 2013-10-18 11:10

8#的Mathematica代码还有较大优化余地，s2,s3,s5去了也行

1. Clear["Global`*"];
   
2. fun = Compile[{{n0, _Integer}}, Block[{n = n0, s2, s3, s5, out},
   
3.     s2 = s3 = s5 = 0;
   
4.     While[Mod[n, 2] == 0, s2++; n = Quotient[n, 2]];
   
5.     While[Mod[n, 3] == 0, s3++; n = Quotient[n, 3]];
   
6.     While[Mod[n, 5] == 0, s5++; n = Quotient[n, 5]];
   
7.     n
   
8.     ], RuntimeAttributes -> Listable, RuntimeOptions -> "Speed",
   
9.    CompilationTarget -> "C"
   
10.    ];
    
11. 
    
12. Pick[#, fun@#, 1] &@Range[5 10^5]~Take~200 // AbsoluteTiming

复制代码

mathematica

chyanog 发表于 2013-10-21 12:21
8#的Mathematica代码还有较大优化余地，s2,s3,s5去了也行

Compile:ptx: Unknown option RuntimeAttributes in Compile[{{n0,_Integer}},Block[{n=n0,s2,s3,s5,out},s2=s3=s5=0;While[Mod[n,2]==0,s2++;n=Quotient[n,2]];While[Mod[n,3]==0,s3++;n=Quotient[n,3]];While[Mod[n,5]==0,s5++;n=Quotient[n,5]];n],RuntimeAttributes->Listable,RuntimeOptions->Speed,CompilationTarget->C]. >>

Pick::incomp: Expressions {1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,<<499950>>} and Compile[{{n0,_Integer}},Block[{n=n0,s2,s3,s5,out},s2=s3=s5=0;While[Mod[n,2]==0,s2++;n=Quotient[n,2]];While[Mod[n,3]==0,s3++;n=Quotient[n,3]];While[Mod[n,5]==0,s5++;n=Quotient[n,5]];n],RuntimeAttributes->Listable,RuntimeOptions->Speed,CompilationTarget->C][{1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,<<499950>>}] have incompatible shapes. >>

Take::take: Cannot take positions 1 through 200 in Pick[{1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,<<499950>>},Compile[{{n0,_Integer}},Block[{n=n0,s2,s3,s5,out},s2=s3=s5=0;While[Mod[n,2]==0,s2++;n=Quotient[<<2>>]];While[Mod[n,3]==0,s3++;n=Quotient[<<2>>]];While[Mod[n,5]==0,s5++;n=Quotient[<<2>>]];n],RuntimeAttributes->Listable,RuntimeOptions->Speed,CompilationTarget->C][{1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,<<499950>>}],1]. >>

我运行你的代码有这些错误,不知道啥原因.


我的是mathematica 7.0.1.0 winxp系统

chyanog

mathematica 发表于 2013-10-21 13:39
Compile:ptx: Unknown option RuntimeAttributes in Compile[{{n0,_Integer}},Block[{n=n0,s2,s3,s5,ou ...

我的代码在Mathematica8和9上速度很快（不到0.1s），7.0还不支持Compile的新选项。建议升级Mathematica版本，我曾经做过测试，8.0和9.0的Pick和Compile比7.0快得多，如果你一定要用7.0可以运行下面代码（0.7s）

1. Clear["Global`*"];
   
2. fun = Compile[{{n0, _Integer}},
   
3.    Block[{n = n0, s2, s3, s5, out}, s2 = s3 = s5 = 0;
   
4.     While[Mod[n, 2] == 0, s2++; n = Quotient[n, 2]];
   
5.     While[Mod[n, 3] == 0, s3++; n = Quotient[n, 3]];
   
6.     While[Mod[n, 5] == 0, s5++; n = Quotient[n, 5]];
   
7.     n]];
   
8. 
   
9. Pick[#, fun /@ #, 1] &@Range[5 10^5]~Take~200 // AbsoluteTiming

复制代码

mathe

gxqcn 发表于 2013-10-21 11:31
无需排序，也无需任何除法运算，就能直接构造出来的。

请注意到一个事实：

我怎么觉得过去讨论过的。O(n)的时间和空间复杂度即可

mathe

不知道是否有快于O(n)的算法来计算序列中第n个数，或者给定序列中一个数，计算其在数列中的序数

hujunhua

我上周末也在考虑这种算法，这个规律应该不是那么容易揣摸的。所以我采取的是一种“马后炮”式的算法。
考虑这个序列中相邻两项的比值，在前200项数据中统计199个相邻比，可知只有17个不同者，它们从小到大依次列于表1：
表1$={15625/15552, 2048/2025, 81/80, 3125/3072, 128/125, 250/243, 25/24,135/128, 16/15, 27/25, 10/9, 9/8, 6/5, 5/4, 4/3, 3/2, 2}$

用一个三元组(a,b,c)来记录相邻比中的2，3，5的幂指，例如: r(6/5)={1,1,-1}, r(10/9)={1,-2,1}等, 那么上述17个相邻比列于表2：
表2$={{-6,-5,6},{11,-4,-2},{-4,4,-1},{-10,-1,5},{7,0,-3},{1,-5,3},{-3,-1,2},{-7,3,1},{4,-1,-1},{0,3,-2},{1,-2,1},{-3,2,0},{1,1,-1},{-2,0,1},{2,-1,0},{-1,1,0},{1,0,0}}$.

原序列也用这样三元组来表示，前40项列于表3
表3$={{0,0,0},{1,0,0},{0,1,0},{2,0,0},{0,0,1},{1,1,0},{3,0,0},{0,2,0},{1,0,1},{2,1,0},{0,1,1},{4,0,0},{1,2,0},{2,0,1},{3,1,0},{0,0,2},{0,3,0},{1,1,1},{5,0,0},{2,2,0},$
${3,0,1},{0,2,1},{4,1,0},{1,0,2},{1,3,0},{2,1,1},{6,0,0},{3,2,0},{0,1,2},{4,0,1},{0,4,0},{1,2,1},{5,1,0},{2,0,2},{2,3,0},{3,1,1},{0,0,3},{7,0,0},{0,3,1},{4,2,0}}$

在表3中，由前项${r_1,s_1,t_1}$生成后项${r_2,s_2,t_2}$的方法是：${r_2,s_2,t_2}={r_1,s_1,t_1}+{a_i,b_i,c_i}$
其中${a_i,b_i,c_i}$是表2中第1个使得${r_2,s_2,t_2}$各项皆为非负的三元组。

这个算法的问题是，表17如何先验地得到，以及它能够生成后多少项。

gxqcn

刚刚新鲜出炉的代码：

1. #include <stdio.h>
   
2. 
   
3. typedef unsigned int UINT;
   
4. 
   
5. 
   
6. // 快速搜素前 $n$ 个 $2^r*3^s*5^t$ 形式的自然数序列；返回有效数目
   
7. static UINT Search( UINT array[], UINT array_size )
   
8. {
   
9.         UINT *p, *p2, *p3, *p5;
   
10.         UINT count = 1, x2, x3, x5;
    
11. 
    
12.         *( p = p2 = p3 = p5 = &array[0] ) = 1;
    
13.         x2 = *p2 * 2;
    
14.         x3 = *p3 * 3;
    
15.         x5 = *p5 * 5;
    
16. 
    
17.         do
    
18.         {
    
19.                 *++p = ( x2 < x3 ) ? x2 : x3;
    
20.                 if ( x5 < *p ) *p = x5;
    
21.                 if ( 0xFFFFFFFFu == *p ) break;
    
22. 
    
23.                 if ( x2 == *p ){ x2 = *++p2 * 2;        if ( x2 < *p ) x2 = 0xFFFFFFFFu; }
    
24.                 if ( x3 == *p ){ x3 = *++p3 * 3;        if ( x3 < *p ) x3 = 0xFFFFFFFFu; }
    
25.                 if ( x5 == *p ){ x5 = *++p5 * 5;        if ( x5 < *p ) x5 = 0xFFFFFFFFu; }
    
26.         } while ( ++count < array_size );
    
27. 
    
28.         return count;
    
29. }
    
30. 
    
31. 
    
32. // 输出
    
33. static void Print_Output( UINT array[], UINT array_size )
    
34. {
    
35.         UINT k = 0;
    
36. 
    
37.         printf( "{ %u", array[k] );
    
38.         while ( array_size != ++k )
    
39.         {
    
40.                 printf( ", %u", array[k] );
    
41.         }
    
42.         printf( " }\n" );
    
43. }
    
44. 
    
45. 
    
46. int main( void )
    
47. {
    
48. #        define ARRAY_SIZE                2048
    
49. 
    
50.         UINT array[ ARRAY_SIZE ];
    
51.         UINT array_size = ARRAY_SIZE;
    
52. 
    
53.         array_size = Search( array, array_size );
    
54.         Print_Output( array, array_size );
    
55. 
    
56.         return 0;
    
57. }

复制代码

说明：

• 提升难度：搜索全部可用uint32数据类型表达的$2^r*3^s*5^t$形式的数，实际可搜得 1848 项（最大值 = 4,271,484,375 = 3^7 * 5^9）；
• 以上代码兼容 C/C++ 编译；
• 效率非常之高（数据输出反而为主要占时）。

gxqcn

基本原理为：

• 追加项，即新项 *p=min{ x2, x3, x5 };
• 移动待乘项指针：若 x_k==*p，则 ++p_k.  以使新的对应的 x_k>*p 继续成为下一次的候选值。
  

northwolves

VBA代码，速度也能接受：

Sub Main()
Dim num&, k&, n&, h(), x, y, i(2)
num = 6000
x = Array(2, 3, 5)
y = Array(2, 3, 5)
ReDim h(num - 1)
h(0) = 1
For n = 1 To num - 1
h(n) = y(0)
For k = 1 To 2
If h(n) > y(k) Then h(n) = y(k)
Next
For k = 0 To 2
If h(n) = y(k) Then i(k) = i(k) + 1: y(k) = x(k) * h(i(k))
Next
Next
MsgBox h(num - 1)
End Sub

返回408146688000000

northwolves

这个链接http://rosettacode.org/wiki/Hamming_numbers#PARI.2FGP 的代码：

(22:18) gp >
(22:18) gp > allocatemem (10 ^ 8)
  ***   Warning: new stack size = 100000000 (95.367 Mbytes).
(22:18) gp > Hupto(n)={
  my(v=vector(n),x2=2,x3=3,x5=5,i=1,j=1,k=1,t);
  v[1]=1;
  for(m=2,n,
    v[m]=t=min(x2,min(x3,x5));
    if(x2 == t, x2 = v[i++] << 1);
    if(x3 == t, x3 = 3 * v[j++]);
    if(x5 == t, x5 = 5 * v[k++]);
  );
  v
};
(22:18) gp > h(n)=Hupto(n)[n];
(22:18) gp > h (10 ^ 6)
%31 = 51931278044838873608958984375000000000000000000000000000000000000000000000
0000000000