论坛 › 编程擂台 › 求形如 2^r*3^s*5^t 的整数序列

mathematica 发表于 2013-10-22 09:21:44

http://oeis.org/A051037

gxqcn 发表于 2013-10-22 09:39:16

mathematica 发表于 2013-10-22 09:23
讨论的是别人剩下来的东西

请注意你的措词！
什么是“别人剩下来的东西”？！

首先，问题本身是个纯数学问题，不可能是什么新问题；但作为算法题，大家来PK还是很不错的（这才是关键）！
其次，发这个问题之前，我并不知道别处是否有过类似讨论，因为我是将数年前在CSDN上看到的题目再推广的。

发这个帖子的目的是为了追求极致的算法，即如何快速解决问题的过程，而非解决问题本身。
所以，何来“别人剩下来的东西”之说？

mathematica 发表于 2013-10-22 09:44:18

算法求解: 2^i * 3^j * 5^k i,j,k >=0. 如何用程序打印出下列的数 1 2 3 4 5 6 8 9 10 12 15 16 18 20 ..

http://zhidao.baidu.com/question/299302012

mathematica 发表于 2013-10-22 14:54:17

http://acm.fjnu.edu.cn/problem/1426
#include<stdio.h>
#include<iostream.h>
int min(int a,int b)
{
        return (a>b?b:a);
}
int getmin(int a,int b,int c)
{
        return min(min(a,b),c);
}
void f()
{
        int i2,i3,i5,i,a2,a3,a5,n;
        while(scanf("%d%d%d%d",&a2,&a3,&a5,&n)!=EOF)
        {
                int *a=new int;
                a=1;
                i2=i5=i3=0;
                for(i=1;i<=n;i++)
                {
                  a[ i ]=getmin(a2*a,a3*a,a5*a);
                  if(a[ i ]==a2*a)
                          i2++;
                  if(a[ i ]==a3*a)
                          i3++;
                  if(a[ i ]==a5*a)
                          i5++;
                }
                printf("%d\n",a);
                delete a;
        }
}
int main()
{
        f();
        return 1;
}我把代码复制过来,相当于备份一下
http://zhidao.baidu.com/question/299302012

wayne 发表于 2013-10-22 21:25:16

northwolves 发表于 2013-10-21 22:21
这个链接http://rosettacode.org/wiki/Hamming_numbers#PARI.2FGP 的代码：

(22:18) gp >


神链接。#include <iostream>
#include <vector>
// Hamming like sequences Generator
//
// Nigel Galloway. August 13th., 2012
//
class Ham {
private:
        std::vector<unsigned int> _H, _hp, _hv, _x;
public:
        bool operator!=(const Ham& other) const {return true;}
        Ham begin() const {return *this;}
      Ham end() const {return *this;}
        unsigned int operator*() const {return _x.back();}
        Ham(const std::vector<unsigned int> &pfs):_H(pfs),_hp(pfs.size(),0),_hv({pfs}),_x({1}){}
        const Ham& operator++() {
          for (int i=0; i<_H.size(); i++) for (;_hv<=_x.back();_hv=_x[++_hp]*_H);
          _x.push_back(_hv);
          for (int i=1; i<_H.size(); i++) if (_hv<_x.back()) _x.back()=_hv;
          return *this;
        }
};

int main() {
int count = 1;
for (unsigned int i : Ham({2,3,5})) {
    if (count <= 62) std::cout << i << ' ';
    if (count++ == 1691) {
      std::cout << "\nThe one thousand six hundred and ninety first Hamming Number is " << i << std::endl;
      break;
    }
}
return 0;
}

hujunhua 发表于 2013-10-23 02:37:33

Mathematica程序，不显含 if 语句。

本帖最后由 hujunhua 于 2013-10-23 16:22 编辑

list={};(*预备一个空表来装要求的序列*)
X={2,3,5};(*倍乘因子表*)
checkb={1,1,1};(*备选框*)
pointer={0,0,0};(*倍乘因子指针表，指向已得的表list中的位置*)
Do;(*将备选框的最小数追加到表list末*)
dp=1-Sign@(checkb-Last@list);(*用1标记上一步所选最小值在备选框中的位置，其它位置标为0，作为指针的移动表*)
pointer+=dp;(*按指针移动表移动指针*)
checkb+=(Thread]*dp*X]-Thread),(*将有移动的指针所指的表list中的项乘以对应的倍乘因子，以所得结果按标记更新备选框*)
{200}](*返回从备选框中选最小数入表的那一步，直至表长达到200*)
Print是由楼上链结中的程序改编的，效率当然也很高喽。

mathematica 发表于 2013-10-23 09:29:59

这个,每一项必然是前面几项中某一项乘以2或者3或者5而得到的

liangbch 发表于 2013-10-23 17:10:06

这个题目的另一个说法是“求5-smooth数",百度知道(34#)给出的代码简单而有效，但有其局限性。用此代码求更大范围或者更大的质数的smooth数就不行了。其主要局限性是，
1.这个算法必须从1开始求所有smooth数，而不能求某一区间的smooth数。
2.这个算法所有求得的smooth数必须保存在内存。更大范围或者更大的质数的smooth数的个数会很多，可能导致内存不够用。不得不使用内存映射文件，或文件来存储,这会导致程序很复杂。

我正好需要计算大量的smooth数，鉴于磁盘和内存的限制，我不得不使用普通的做法，和素数的分段筛很类似，可分段求smooth数，不过求得需要的smooth数后需要排序操作，另外，我的这个算法好处是，便于分布式计算。

许多算法看上去很美，但很有局限性。和此算法一个类似的例子是《一个线性O(N)时间复杂度求素数筛法》(http://www.cnblogs.com/lzhenf/archive/2011/12/24/2300271.html),但使用后发现，虽然复杂度确实很小，但实际性能不理想，因为程序中对内存的访问局部性很差，导致cache命中率很低，限于这个算法特点，所有的所有素数必须保存在内存，故不能求很大范围内的素数。

下表是我2^64以内的p1到p32-smooth数的个数。


n
Pn
2^64以内的
Pn-smooth数的个数

1
2
64

2
3
1344

3
5
13282

4
7
85348

5
11
378554

6
13
1381206

7
17
4164316

8
19
11169655

9
23
26587554

10
29
56643654

11
31
113399428

12
37
210499995

13
41
370466869

14
43
627484292

15
47
1020181037

16
53
1591363196

17
59
2395982361

18
61
3527728592

19
67
5050259345

20
71
7077350886

21
73
9759368423

22
79
13183354386

23
83
17528499627

24
89
22918160232

25
97
29457339259

26
101
37424970530

27
103
47121505151

28
107
58735924658

29
109
72644250191

30
113
89058878497

31
127
107598623619

32
131
129070232683

mathematica 发表于 2013-10-23 17:42:17

#include <stdio.h>
#include <limits.h>
#include <string.h>
#include <malloc.h>

#define MAX_LEVEL                32

typedef unsigned long DWORD;

#if defined(_MSC_VER)
        typedef unsigned __int64 UINT64;
        #define UINT64_FMT"%I64u"

#elif defined(__GNUC__)
        typedef unsigned long long UINT64;
        #define UINT64_FMT"%llu"
#else
        #error don't support this compiler
#endif

const DWORD primes[]=
{
       2,3,5,7, 11, 13, 17, 19,
        23, 29, 31, 37, 41, 43, 47, 53,
        59, 61, 67, 71, 73, 79, 83, 89,
        97,101,103,107,109,113,127,131
};

typedefstruct _searchContent_tag
{
        struct search_tag
        {
                UINT64 PI;                // product of some numbers
                UINT64 max;
                int idx;                // the index in array pArr
                int maxIdx;                // the index of last element
                UINT64* pArr;        // the address of primes power table
        }arr;
        UINT64 min;
        UINT64 max;
        intmaxLevel;                //level
}SRCH_CONTENT_ST;

void initSearchContent(int level, UINT64 min, UINT64 max, SRCH_CONTENT_ST* p)
{
        UINT64 prime,limit,pp;
        int i,c;
       
        memset(p,0,sizeof(SRCH_CONTENT_ST));
       
        p->maxLevel=level-1;
        p->max=max;         p->min=min;
        p->arr.PI=1;
        p->arr.max=max;
        for (i=0;i<level;i++)
        {
                prime=primes;
                limit=ULLONG_MAX/prime;
               
                c=1;pp=1;
                while (pp<=limit)        //pp is power of prime
                { pp*=prime;c++; }
               
                p->arr.pArr=(UINT64 *)malloc( sizeof(UINT64) * (c+1));
               
                p->arr.pArr=pp=1;        c=1;
                while (pp<=limit)
                { pp*=prime;p->arr.pArr=pp; }
               
                p->arr.maxIdx= c-1;
        }
}

void freeSearchContent( SRCH_CONTENT_ST* p)
{
        int i;
        for (i=0;i<=p->maxLevel;i++)
        {
                if ( p->arr.pArr!=NULL)
                        free(p->arr.pArr);
        }
        memset(p,0,sizeof(SRCH_CONTENT_ST));
}

//recursive form
void searchSmoothNums1( SRCH_CONTENT_ST* p, int order)
{
        int    idx;
        /*
          recurrence formula
                p->arr.PI= 1;
                p->arr.PI    = arr.pArr[ arr.idx] * arr.PI;       
                p->arr.max= p->max / p->arr.PI;
        */

        if (order==0)
        {
                UINT64 min;
                UINT64 prod;

                min= (p->min + p->arr.PI -1)/p->arr.PI;// min >= p->min/PI( when p->min % PI==0, min= p->min/PI)
               
                if ( min > p->arr.max)
                        return;

                idx=0;
                while (idx < p->arr.maxIdx && p->arr.pArr<min )
                        idx++;
               
                while ( p->arr.pArr <= p->arr.max && idx <= p->arr.maxIdx )
                {
                        prod= p->arr.PI * p->arr.pArr;
                        idx++;
                        printf(UINT64_FMT"\n",prod);
                }
        }
        else
        {
                if ( p->arr.max==1 )
                {
                        printf(UINT64_FMT"\n",p->arr.PI);
                        return ;
                }
               
                for (idx=0;
                        idx<= p->arr.maxIdx&& p->arr.pArr <= p->arr.max;
                        idx++ )
                {
                        p->arr.idx=idx;
                        p->arr.PI = p->arr.PI * p->arr.pArr;
                        p->arr.max= p->max / p->arr.PI;
                        searchSmoothNums1(p,order-1);   //call itselft
                }
               
        }
}


#define M_FILL   0    //fill mode
#define M_INC    1    //increase mode

// iterative form
// search all p smooth number which between p->min and p->max
// For exmaple
//   order=0, search 2-smooth numbers
//   order=1, search 3-smooth numbers
//   order=2, search 5-smooth number
void searchSmoothNums2( SRCH_CONTENT_ST* p, int order)
{
        UINT64 prod;
        UINT64 PI;
        UINT64 min;
        int idx;

        int i=order;
        int mode=M_FILL;
        /*
          recurrence formula
                p->arr.PI= 1;
                p->arr.PI   = arr.pArr[ arr.idx] * arr.PI;       
                p->arr.max    = p->max / arr.PI;       
        */
       
        p->arr.PI=1;
        p->arr.max= p->max;

        while (i<=order)
        {
               
                if (mode==M_FILL)//fill mode
                {
                        PI = p->arr.PI;
                        if (i>0)
                        {
                                p->arr.idx=0;
                                p->arr.PI = p->arr.PI;        //p->arr.PI= p->arr.PI;
                                p->arr.max= p->arr.max;
                                i--;                                                //forward to next level
                        }
                        else// i==0
                        {
                                min= (p->min + PI-1)/PI;        // min >= p->min/PI( when p->min % PI==0, min= p->min/PI)
                       
                                if ( min > p->arr.max)
                                {
                                        i++;                                        // backtrace
                                        mode=M_INC;
                                }
                                else
                                {
                                        idx=0;
                                        while (idx < p->arr.maxIdx && p->arr.pArr<min )
                                                idx++;
                                        if (p->arr.pArr > p->arr.max || idx > p->arr.maxIdx)
                                                i++;                                        // backtrace
                                        else   //in this case, varible i don't change
                                        {
                                                p->arr.idx=idx;
                                                printf(UINT64_FMT"\n",PI * p->arr.pArr);

                                        }
                                        mode=M_INC;                //change to increase mode
                                }
                        }
                }
                else//is increase mode
                {
                        idx= p->arr.idx + 1;
                        if (p->arr.pArr > p->arr.max || idx> p->arr.maxIdx)
                                i++;                // backtrace, mode don't change, still is increase mode
                        else
                        {
                                p->arr.idx=idx;
                                prod= p->arr.PI * p->arr.pArr;

                                if (i==0)
                                        printf(UINT64_FMT"\n",prod);
                                else
                                {
                                        p->arr.max= p->max / prod;
                                        if ( p->arr.max ==1)
                                        {
                                                printf(UINT64_FMT"\n",prod);
                                                i++;                        // backtrace, mode don't change, still is increase mode
                                        }
                                        else
                                        {
                                                p->arr.PI=prod;
                                                mode=M_FILL;        // change mode to fill mode
                                                i--;                        // forward to next level
                                        }
                                }
                        }
                }
        }
}

void ref_searchSmoothNums(int order,UINT64 min, UINT64 max)
{
        UINT64 i,m;
        int j,isSmooth;

        for (i=min;i<=max;i++)
        {
                for (m=i,isSmooth=0,j=0;j<order;j++)
                {
                        while (m % primes==0) m/=primes;
                        if (m==1)
                        {        isSmooth=1;        break;        }
                }
                if (isSmooth)
                        printf(UINT64_FMT"\n",i);
        }
}

void test()
{
        SRCH_CONTENT_ST searchContent;
        UINT64 nMin,nMax;
        int cs;//case

        for (cs=0;cs<=1;cs++)
        {
                if (cs==0)
                { nMin=1;nMax=100; }
                else if (cs==1)
                { nMin=100;        nMax=200;        }

                initSearchContent(4,nMin,nMax,&searchContent);                       
               
                printf("\n\nSearch smooth number ref_searchSmoothNums-------\n");
                ref_searchSmoothNums(4,nMin,nMax);
               
                printf("\n\nSearch smooth number searchSmoothNums1-------\n");
                searchSmoothNums1( &searchContent,4-1);
               
                printf("\n\nSearch smooth number searchSmoothNums2-------\n");
                searchSmoothNums2( &searchContent,4-1);

                freeSearchContent( &searchContent);
        }
}

int main(int argc, char* argv[])
{
        test();
        return 0;
}1.什么是Smooth 数

如果一个整数的所有素因子都不大于B，我们称这个整数为B-Smooth数。用符号Pn表示表示第n个素数，例P1=2,P2=3,P3=5. 我们称素因子不大于Pn的Smooth数叫Pn-Smooth数。关于Smooth数更多的信息，请请参阅 Wiki词条 Smooth number http://en.wikipedia.org/wiki/Smooth_number.P1 到 P9-Smooth 数可参阅下面的数列，这些数列在On-Line Encycolpedia of Integer Sequences http://oeis.org 给出

    2-smooth numbers: A000079 (2^i, i>=0)

    3-smooth numbers: A003586 (2^i ×3^j,i,j>=0)

    5-smooth numbers: A051037 (2^i ×3^j×5^k, i,j,k>=0)

    7-smooth numbers: A002473 (2^i×3^j×5^k×7^l, i,j,k,l>=0)

    11-smooth numbers: A051038 (etc...)

    13-smooth numbers: A080197

    17-smooth numbers: A080681

    19-smooth numbers: A080682

    23-smooth numbers: A080683



2. 如何得到smooth数

基本上，我们有2种方法来生成一定范围内的smooth数。

1.   穷举法。基本方法是对指定范围内的每个数做分解素因数，如果其所有素因子都小于等于Pn,那个这个数就是Pn-smooth数。这种方法的优点是简单，缺点是速度较慢，只适合检查小范围内的整数。



2．构造法。其基本方法是用一个数组表示各个素因子的指数。例arr是一个数组，arr表示素数p1（2）的指数，arr表示素数p2（3）的指数，arr表示素数p3（5）的指数. 欲求min和max之间的所有pn-smooth数, 则需要列举枚举这个数组的各个元素的取值。数组的各个元素的值需满足如下2个条件。

1. p[ i]>=0, 且p[ i] ^ arr[ i] <= max.（p[ i]表示第i个素数，i>=1,i<=n）

2. min <=p^arr * p^arr * p^arr * …p*arr <=max

对于一个包含n个元素的数组arr, 每个元素arr[ i] （1<=i<=n）各取一个数，就得到一个n维向量。如arr=0, arr=0, arr=0可用3维向量(0,0,0)来表示，而arr=1, arr=0, arr=0 则用3维向量(1,0,0)来表示。为了使得取得的所有n维向量既不重复也不遗漏。我们这里定义了比较2个n维向量的大小的规则

1.如果2个n维向量a,b的所有分量都是相同的，我们说这两个n维向量是相同的。

2.对于2个n维向量a,b（a的各个分量为a,a, …, a）,如果a>b我们说a>b.

3.如果存在某个i，使得（1）对于所有的j (i<j<=n)都有a=b. (2) a[ i]>b[ i],那么我们说n维向量a>b.

我们列举每个n维向量时，总是从小到到的进行，先是(0,0,0), 然后是(1,0,0), 然后是(2,0,0)…，然后是(0,1,0), (1,1,0), (2,1,0).   

未完待续

http://blog.csdn.net/liangbch/article/details/12186099

chyanog 发表于 2013-10-23 17:48:53

本帖最后由 chyanog 于 2013-10-23 17:54 编辑

hujunhua 发表于 2013-10-23 02:37
是由楼上链结中的程序改编的，效率当然也很高喽。
当list很大的时候，AppendTo会越来越慢，不如用Internal`BagModule[{list, X, dp, checkb, pointer},
   list = {};
   X = {2, 3, 5};
   checkb = {1, 1, 1};
   pointer = {0, 0, 0};
   Do;
    dp = 1 - Sign@(checkb - Last@list);
    pointer += dp;
    checkb += (list[]*dp*X - checkb*dp)
    , {2 10^4}];
   list
   ] // Hash // AbsoluteTiming
{5.515315, 1434246505}\$ContextPath = Union@Append[\$ContextPath, "Internal`"];
Module[{bag, X, dp, checkb, pointer},
   bag = Bag[];
   X = {2, 3, 5};
   checkb = {1, 1, 1};
   pointer = {0, 0, 0};
   Do[
    StuffBag;
    dp = 1 - Sign@(checkb - BagPart);
    pointer += dp;
    checkb += (BagPart*dp*X - checkb*dp),
    {2 10^4}];
   bag~BagPart~All
   ] // Hash // AbsoluteTiming{0.905052, 1434246505}

查看完整版本: 求形如 2^r*3^s*5^t 的整数序列