8#的Mathematica代码还有较大优化余地，s2,s3,s5去了也行Clear["Global`*"];
fun = Compile[{{n0, _Integer}}, Block[{n = n0, s2, s3, s5, out},
s2 = s3 = s5 = 0;
While == 0, s2++; n = Quotient];
While == 0, s3++; n = Quotient];
While == 0, s5++; n = Quotient];
n
], RuntimeAttributes -> Listable, RuntimeOptions -> "Speed",
CompilationTarget -> "C"
];
Pick[#, fun@#, 1] &@Range~Take~200 // AbsoluteTiming chyanog 发表于 2013-10-21 12:21
8#的Mathematica代码还有较大优化余地，s2,s3,s5去了也行
Compile::optx: Unknown option RuntimeAttributes in Compile[{{n0,_Integer}},Block[{n=n0,s2,s3,s5,out},s2=s3=s5=0;While==0,s2++;n=Quotient];While==0,s3++;n=Quotient];While==0,s5++;n=Quotient];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==0,s2++;n=Quotient];While==0,s3++;n=Quotient];While==0,s5++;n=Quotient];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==0,s2++;n=Quotient[<<2>>]];While==0,s3++;n=Quotient[<<2>>]];While==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系统 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）Clear["Global`*"];
fun = Compile[{{n0, _Integer}},
Block[{n = n0, s2, s3, s5, out}, s2 = s3 = s5 = 0;
While == 0, s2++; n = Quotient];
While == 0, s3++; n = Quotient];
While == 0, s5++; n = Quotient];
n]];
Pick[#, fun /@ #, 1] &@Range~Take~200 // AbsoluteTiming gxqcn 发表于 2013-10-21 11:31
无需排序，也无需任何除法运算，就能直接构造出来的。
请注意到一个事实：
我怎么觉得过去讨论过的。O(n)的时间和空间复杂度即可 不知道是否有快于O(n)的算法来计算序列中第n个数，或者给定序列中一个数，计算其在数列中的序数
RE: 20#gxqcn
我上周末也在考虑这种算法，这个规律应该不是那么容易揣摸的。所以我采取的是一种“马后炮”式的算法。考虑这个序列中相邻两项的比值，在前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如何先验地得到，以及它能够生成后多少项。 刚刚新鲜出炉的代码：#include <stdio.h>
typedef unsigned int UINT;
// 快速搜素前 $n$ 个 $2^r*3^s*5^t$ 形式的自然数序列；返回有效数目
static UINT Search( UINT array[], UINT array_size )
{
UINT *p, *p2, *p3, *p5;
UINT count = 1, x2, x3, x5;
*( p = p2 = p3 = p5 = &array ) = 1;
x2 = *p2 * 2;
x3 = *p3 * 3;
x5 = *p5 * 5;
do
{
*++p = ( x2 < x3 ) ? x2 : x3;
if ( x5 < *p ) *p = x5;
if ( 0xFFFFFFFFu == *p ) break;
if ( x2 == *p ){ x2 = *++p2 * 2; if ( x2 < *p ) x2 = 0xFFFFFFFFu; }
if ( x3 == *p ){ x3 = *++p3 * 3; if ( x3 < *p ) x3 = 0xFFFFFFFFu; }
if ( x5 == *p ){ x5 = *++p5 * 5; if ( x5 < *p ) x5 = 0xFFFFFFFFu; }
} while ( ++count < array_size );
return count;
}
// 输出
static void Print_Output( UINT array[], UINT array_size )
{
UINT k = 0;
printf( "{ %u", array );
while ( array_size != ++k )
{
printf( ", %u", array );
}
printf( " }\n" );
}
int main( void )
{
# define ARRAY_SIZE 2048
UINT array[ ARRAY_SIZE ];
UINT array_size = ARRAY_SIZE;
array_size = Search( array, array_size );
Print_Output( array, array_size );
return 0;
}说明：
[*]提升难度：搜索全部可用uint32数据类型表达的2^r*3^s*5^t形式的数，实际可搜得 1848 项（最大值 = 4,271,484,375 = 3^7 * 5^9）；
[*]以上代码兼容 C/C++ 编译；
[*]效率非常之高（数据输出反而为主要占时）。 基本原理为：
[*]追加项，即新项 *p=min{ x2, x3, x5 };
[*]移动待乘项指针：若 xk==*p，则 ++pk.以使新的对应的 xk>*p 继续成为下一次的候选值。
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
这个链接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;
for(m=2,n,
v=t=min(x2,min(x3,x5));
if(x2 == t, x2 = v << 1);
if(x3 == t, x3 = 3 * v);
if(x5 == t, x5 = 5 * v);
);
v
};
(22:18) gp > h(n)=Hupto(n);
(22:18) gp > h (10 ^ 6)
%31 = 51931278044838873608958984375000000000000000000000000000000000000000000000
0000000000