数学研发论坛

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

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

[复制链接]
发表于 2013-10-21 12:21:11 | 显示全部楼层


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. Pick[#, fun@#, 1] &@Range[5 10^5]~Take~200 // AbsoluteTiming
复制代码

点评

@mathematica,看这里,http://reference.wolfram.com/mathematica/ref/CompilationTarget.html  发表于 2013-10-21 21:12
CompilationTarget -> "C"  发表于 2013-10-21 19:10
这个是编译成C语言的意思吗?  发表于 2013-10-21 19:10
我明白了,&与前面的#构成纯函数  发表于 2013-10-21 15:47
代码第12行的&@是啥意思?  发表于 2013-10-21 13:46
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2013-10-21 13:39:16 | 显示全部楼层
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系统
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2013-10-21 15:55:11 | 显示全部楼层
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. Pick[#, fun /@ #, 1] &@Range[5 10^5]~Take~200 // AbsoluteTiming
复制代码
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2013-10-21 17:34:13 | 显示全部楼层
gxqcn 发表于 2013-10-21 11:31
无需排序,也无需任何除法运算,就能直接构造出来的。

请注意到一个事实:

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

点评

确实,在CSDN上n年前曾讨论过前100个$2^r*3^s$形式的自然数序列。  发表于 2013-10-21 20:55
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2013-10-21 17:37:15 | 显示全部楼层
不知道是否有快于O(n)的算法来计算序列中第n个数,或者给定序列中一个数,计算其在数列中的序数

点评

肯定有。计算序数也许只需$O(k^{1/3})$的时间,计算第$n$个数也许只需$O(n^{1/3}*\log n)$的时间。  发表于 2013-10-21 18:54
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2013-10-21 17:58:49 | 显示全部楼层

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如何先验地得到,以及它能够生成后多少项。

点评

我觉得与素数定理有关系,应该不是那么简单的  发表于 2013-10-21 19:21
对于前n个自然数,算出符合条件的个数f(n),然后运用曲线拟合的技术,看看f(n)的表达式如何,你觉得呢?  发表于 2013-10-21 19:20
先把代码递上来!然后再考虑别的问题  发表于 2013-10-21 19:09
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2013-10-21 20:56:43 | 显示全部楼层
刚刚新鲜出炉的代码:
  1. #include <stdio.h>

  2. typedef unsigned int UINT;


  3. // 快速搜素前 $n$ 个 $2^r*3^s*5^t$ 形式的自然数序列;返回有效数目
  4. static UINT Search( UINT array[], UINT array_size )
  5. {
  6.         UINT *p, *p2, *p3, *p5;
  7.         UINT count = 1, x2, x3, x5;

  8.         *( p = p2 = p3 = p5 = &array[0] ) = 1;
  9.         x2 = *p2 * 2;
  10.         x3 = *p3 * 3;
  11.         x5 = *p5 * 5;

  12.         do
  13.         {
  14.                 *++p = ( x2 < x3 ) ? x2 : x3;
  15.                 if ( x5 < *p ) *p = x5;
  16.                 if ( 0xFFFFFFFFu == *p ) break;

  17.                 if ( x2 == *p ){ x2 = *++p2 * 2;        if ( x2 < *p ) x2 = 0xFFFFFFFFu; }
  18.                 if ( x3 == *p ){ x3 = *++p3 * 3;        if ( x3 < *p ) x3 = 0xFFFFFFFFu; }
  19.                 if ( x5 == *p ){ x5 = *++p5 * 5;        if ( x5 < *p ) x5 = 0xFFFFFFFFu; }
  20.         } while ( ++count < array_size );

  21.         return count;
  22. }


  23. // 输出
  24. static void Print_Output( UINT array[], UINT array_size )
  25. {
  26.         UINT k = 0;

  27.         printf( "{ %u", array[k] );
  28.         while ( array_size != ++k )
  29.         {
  30.                 printf( ", %u", array[k] );
  31.         }
  32.         printf( " }\n" );
  33. }


  34. int main( void )
  35. {
  36. #        define ARRAY_SIZE                2048

  37.         UINT array[ ARRAY_SIZE ];
  38.         UINT array_size = ARRAY_SIZE;

  39.         array_size = Search( array, array_size );
  40.         Print_Output( array, array_size );

  41.         return 0;
  42. }
复制代码
说明:
  • 提升难度:搜索全部可用uint32数据类型表达的$2^r*3^s*5^t$形式的数,实际可搜得 1848 项(最大值 = 4,271,484,375 = 3^7 * 5^9);
  • 以上代码兼容 C/C++ 编译;
  • 效率非常之高(数据输出反而为主要占时)。

点评

注:如果不需检测整型溢出的问题,则代码中含“0xFFFFFFFFu”语句均可删除。  发表于 2013-10-22 16:24
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2013-10-21 21:09:20 | 显示全部楼层
基本原理为:
  • 追加项,即新项 *p=min{ x2, x3, x5 };
  • 移动待乘项指针:若 xk==*p,则 ++pk.  以使新的对应的 xk>*p 继续成为下一次的候选值。
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2013-10-21 22:06:25 | 显示全部楼层
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
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2013-10-21 22:21:40 | 显示全部楼层
这个链接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
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
您需要登录后才可以回帖 登录 | 欢迎注册

本版积分规则

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

GMT+8, 2019-7-17 04:28 , Processed in 0.054815 second(s), 22 queries .

Powered by Discuz! X3.4

© 2001-2017 Comsenz Inc.

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