challenge 2^i * 3^j * 5^k なる整数

2^i * 3^j * 5^k の形で表される整数を小さい方から順に 100 個列挙するプログラムを書いてください。 i, j, k は 0 以上の整数です。アルゴリズムのオーダーについても考えてみてください。

例えば最初の 10 個は次のようになります:

 1 = 2^0 * 3^0 * 5^0
 2 = 2^1 * 3^0 * 5^0
 3 = 2^0 * 3^1 * 5^0
 4 = 2^2 * 3^0 * 5^0
 5 = 2^0 * 3^0 * 5^1
 6 = 2^1 * 3^1 * 5^0
 8 = 2^3 * 3^0 * 5^0
 9 = 2^0 * 3^2 * 5^0
10 = 2^1 * 3^0 * 5^1
12 = 2^2 * 3^1 * 5^0

※解答では i, j, k の各値を示す必要はありません。

Posted feedbacks - Flatten

Nested Hidden
何も考えずに。
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
#include <stdio.h>

int main(){
    int m=1;
    int n=0;
    int i,j,k;
    int v;
    
    do{
        v=m;
        i=0;
        j=0;
        k=0;
        while(v%2==0)v/=2,i++;
        while(v%3==0)v/=3,j++;
        while(v%5==0)v/=5,k++;
        if(v==1){
            printf("%d = 2^%d * 3^%d * 5^%d\n",m,i,j,k);
            n++;
        }
        m++;
    }while(n<100);
    return 0;
}

priority queue で。

 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
#include <cstdio>
#include <queue>
using namespace std;

const int DNUM = 100;

int main(void)
{
    int pre = 0, i = 0;
    priority_queue<int, vector<int>, greater<int> > pq;
    pq.push(1);
    
    while(i < DNUM)
    {
        int a = pq.top();
        pq.pop();
        if(a == pre) continue;
        
        printf("%d\n", a);
        pq.push(a * 2);
        pq.push(a * 3);
        pq.push(a * 5);
        
        pre = a;
        i++;
    }
    
    return 0;
}

おなじアルゴリズムになってしまいました。
100個目は1536。
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
public class Int235 {
    public static void main(String[] a) {
        int limit = 100;
        for (int n = 0, i = 1; n < limit; i++) {
            int tmp = i;
            while (tmp % 2 == 0) tmp /= 2;
            while (tmp % 3 == 0) tmp /= 3;
            while (tmp % 5 == 0) tmp /= 5;
            if (tmp == 1) {
                System.out.println(i);
                n++;
            }
        }
    }
}

軽く候補の絞込み(2と3と5の倍数しか回答は無いので)を入れてみました。
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
list max b n
    | (b^n) > max   =   []
    | otherwise = n : list max b (n+1)

scan n = [(n,i,j,k)|i<-list n 2 0,j<-list n 3 0,k<-list n 5 0,n == (2^i*3^j*5^k)]

result = take 100 $ concatMap (scan) $ filter (check) [1..]

check n
    | n == 1            =   True
    | (n `mod` 2) == 0  =   True
    | (n `mod` 3) == 0  =   True
    | (n `mod` 5) == 0  =   True
    | otherwise         =   False

main=do
    mapM (\x-> putStrLn $ format x) result
    where
        format (n,i,j,k) = concat [show n," = ",fmt 2 i," * ",fmt 3 j," * ",fmt 5 k]
        fmt b n = show b ++ "^" ++ show n

i,j,k は示さなくてもいいそうなので、簡単に。
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
#include <iostream>

int main()
{
  int n = 0;
  for ( int i = 1; n < 100; ++i )
    if ( i == 1 || i % 2 == 0 || i % 3 == 0 || i % 5 == 0 ) {
      std::cout << i << ", ";
      ++n;
    }
  std::cout << "\n";
  return 0;
}

#7638 さんのアルゴリズムを借用。
priority queueの代わりに、標準ライブラリのヒープキューを使います。
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
#!/usr/bin/python

from heapq import heappush, heappop

def doukaku206(num_answers):
    answers = []
    heap = [1]
    prev_answer = 0
    while len(answers) < num_answers:
        a = heappop(heap)
        if a == prev_answer:
            continue
        answers.append(a)
        heappush(heap, a * 2)
        heappush(heap, a * 3)
        heappush(heap, a * 5)
        prev_answer = a
    return answers

print doukaku206(100)

# eof

しまった、これでは12の次が15(3*5)じゃなくて14(2*7)になってしまう。
反省して--self。
そして、まともに動くものを。
 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
#include <iostream>
#include <map>

std::map<int,int> factor(int n)
{
  std::map<int, int> r;
  r[2] = r[3] = r[5] = 0;
  while ( n%2==0 ) { n/=2; ++r[2]; }
  while ( n%3==0 ) { n/=3; ++r[3]; }
  while ( n%5==0 ) { n/=5; ++r[5]; }

  if ( n != 1 ) r.clear();
  return r;
}

int main()
{
  for ( int i = 1, n = 0; n < 100; ++i ) {
    std::map<int,int> f = factor(i);

    if ( !f.empty() ) {
      std::cout << i << " = 2^" << f[2] <<
                        " * 3^" << f[3] <<
                        " * 5^" << f[5] << "\n";
      ++n;
    }
  }
  return 0;
}

Squeak Smalltalk で。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
| count n |
count := 0.
n := 0.
World findATranscript: nil.
[n := n + 1. count < 100] whileTrue: [
    (#(2 3 5) inject: n into: [:quo :each |
            [quo isDivisibleBy: each] whileTrue: [quo := quo / each]. quo]) = 1
         ifTrue: [
            count := count + 1.
            Transcript cr; show: n]]

適当に生成しながら小さい順に並べてます。動けばいい的な作りですが。考え方は 84q さんのと同じでしょうか。

計算量は時間 O(N^2) 空間 O(N) かと思いましたが、実際に試してみた感じだともっと小さいかもしれません。また balanced tree をつかうなど真面目に効率化をやればもっと速くなると思います。

1
2
3
4
5
6
7
8
9
(defun add (n list)
  (if (find n list) list (merge 'list list (list n) #'<)))

(defun h (n)
  (let ((a (list 1)) (c 0))
    (loop (let ((m (pop a)))
            (print m)
            (setf a (add (* m 2) (add (* m 3) (add (* m 5) a))))
            (if (= (incf c) n) (return))))))

別のアルゴリズムで解いてみました。

数が大きくなれば有利かな?

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
public class Sample206 {
    public static void main(String[] args) {
        OUTER: for (int i = 1, n = 0; n < 100; i++) {
            for (int p = 7; p <= i; p += 2) {
                if (p % 3 == 0) continue;
                if (p % 5 == 0) continue;
                if (i % p == 0) {
                    continue OUTER;
                }
            }
            System.out.println(i);
            n++;
        }
    }
}

Javascript@Firebug。 n2,n3,n5の最大値は試行錯誤で決めてます

1
2
3
4
5
6
7
8
var o={};
for(var n2=0; n2<=11; n2++)
  for(var n3=0; n3<=8; n3++)
    for(var n5=0; n5<=7; n5++)
      o[Math.pow(2,n2)*Math.pow(3,n3)*Math.pow(5,n5)] = [n2,n3,n5];
for(var i=0,c=0;c<100;i++)
  if(o[i]!=null)
    console.debug(++c, ". ", i, "= 2^", o[i][0], " + 3^", o[i][1], " + 5^", o[i][2]);

one-linerというかgolfというか...
アルゴリズムは他の皆さんと同じです。
1
2
3
1==shift@z&&++$n&&printf"%d = 2^%d * 3^%d *%s 5^%d\n",@z
while$n<100&&(@z=(++$i,$i,0,0,'',0))&&
map{$z[0]/=$_,++$z[$_]until$z[0]%$_}2,3,5

ちゃんと計算してないけど、時間計算量はO(n logn)、空間計算量はO(1)ぐらい? 再帰してるから空間計算量はO(logn)なのかな。

 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
class P
{
  static void Main(string[] args)
  {
    for (int i = 1, c = 0; c < 100; i++)
    {
      int[] f = F(i);
      if (f != null)
      {
        ++c;
        System.Console.WriteLine
          ("{0} = 2^{1} * 3^{2} * 5^{3}", i, f[2], f[3], f[5]);
      }
    }
  }

  static int[] F(int n)
  {
    if (n == 1) return new int[6];
    int[] a = { 2, 3, 5 };
    foreach (int d in a)
      if (n % d == 0)
      {
        int[] r = F(n / d);
        if (r == null) return null;
        ++r[d];
        return r;
      }
    return null;
  }
}

アプローチとしては#7638と同じでしょうか。配列。

1
2
3
4
5
6
7
8
n=1
for ((i = 0; i < 100; i++));do
    echo $n
    ary[$((n*2))]=1
    ary[$((n*3))]=1
    ary[$((n*5))]=1
    while [ "${ary[$((++n))]}" != 1 ];do :;done
done

ハミング数ですね。
anarchy golf の Hamming Numbers問題
http://golf.shinh.org/p.rb?Hamming+Numbers
のときに投稿した golf用の富豪コードです。

30のx乗 ≡ 0 (mod x) となる x がハミング数に
なることを使っています。(30 は 2,3,5 の最小公倍数)

富豪だから、け、計算量なんて気にしてないんだからね!
1
2
3
4
5
6
7
8
9
main = do
  n <- readLn
  mapM_ print $ take n [x| x <- [1..], mod (30^x) x == 0]

{- anarchy golf に投稿したコードはこちら

main=readLn>>=mapM print.(`take`[x|x<-[1..],mod(30^x)x<1])

-}

素因数分解していくバージョンで
 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
:new
:let s:n=1
:let s:m=0
:while s:m < 100
: if s:n != 1 && s:n % 2 != 0 && s:n % 3 != 0 && s:n % 5 != 0
:   let s:n = s:n + 1
:   continue
: endif
: let s:_ = s:n
: let s:i = 0
: let s:j = 0
: let s:k = 0
: while s:_ % 2 == 0
:   let s:_ = s:_ / 2
:   let s:i = s:i + 1
: endwhile
: while s:_ % 3 == 0
:   let s:_ = s:_ / 3
:   let s:j = s:j + 1
: endwhile
: while s:_ % 5 == 0
:   let s:_ = s:_ / 5
:   let s:k = s:k + 1
: endwhile
: if s:_ == 1
:   call append(s:m, printf("%d = 2^%d * 3^%d * 5^%d", s:n, s:i, s:j, s:k))
:   let s:m = s:m + 1
: endif
: let s:n = s:n + 1
:endwhile
:unlet! s:n
:unlet! s:m
:unlet! s:i
:unlet! s:j
:unlet! s:k
:unlet! s:_

配列の添字が冗長だったので微修正します。

1
2
3
4
5
6
7
8
n=1
for ((i = 0; i < 100; i++));do
    echo $n
    ary[n*2]=1
    ary[n*3]=1
    ary[n*5]=1
    while [ "${ary[++n]}" != 1 ];do :;done
done

 #7652を参考に書いてみました。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
class HummingNumbers
    def self.get(c)
        (30**c % c == 0) ? c : get(c+1)
    end
    def self.take(n,c=1)
        (n==0) ? [] : ((c = get(c)) && take(n-1,c+1).unshift(c))
    end
end

puts HummingNumbers.take(ARGV.length == 1 ? ARGV[0].to_i : 100).join("\n")

 [1..100]>>=penさんの投稿を参考に書いてみました。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
class CHummingNumbers {
    def next(c:Int):Int = (BigInt(30).pow(c) % c).intValue match {
            case 0 => c
            case _ => next(c+1)
        }
    def take(n:Int,c:Int):List[Int] = n match {
            case 0 => List()
            case _ => next(c) match { case v => v::take(n-1,v+1) }
        }
    def take(n:Int):List[Int] = take(n,1)
}
object HummingNumbers {
    def main(args:Array[String]):Unit = {
        try {
            val    n:Int = args.length match {
                case 1 => args(0).toInt
                case _ => 100
            }
            println((new CHummingNumbers).take(n).mkString("\n"))
        } catch {
            case e => e.printStackTrace
        }
    }
}

Haskellらしく素直に無限リストで :)

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
main :: IO ()
main = print $ take 100 hamming

hamming :: [Integer]
hamming = 1 : foldr1 (#) (zipWith map (map (*) [2,3,5]) (repeat hamming))

(#) :: Ord a => [a] -> [a] -> [a]
xxs@(x:xs) # yys@(y:ys) | x < y     = x : xs  # yys
                        | y < x     = y : xxs # ys
                        | otherwise = x : xs  # ys

よく考えたら2で割るのにループはいりませんね。
i,j,kも省いてこんな感じかな。
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
#include <stdio.h>

int main(){
    int m,n,v;
    
    for(n=0,m=1;n<100;m++){
        v=m/(m&-m);
        while(v%3==0)v/=3;
        while(v%5==0)v/=5;
        if(v==1){
            printf("%d\n",m);
            n++;
        }
    }
    return 0;
}

25までのものが出てくるようにDNUMを指定したときに16が出てくるのでしょうか?

(注:2 ^ 4 = 16, 5 ^ 2 = 25)

25までの数を表示するようにDNUMを指定すると、heapに入っている数はi + j + k <= 3を満たすものしかない気がします。しかし、16も指定の数です。


dictでmemo化

本来は配列でmemo化したほうがよいのだが、pythonでListのreserveを適切に行う方法がわからない。index errorをexceptするのもだるいし・・・。

appendするならdictを使うほうがマシでしょう。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
class IsN235WithMemo(object):
  def __init__(self):
    self.memo = dict()
  def __call__(self, n):
    r = self.memo.get(n, None)
    if r is not None:
      return r
    if n == 1:
      return True
    d, m = divmod(n, 2)
    if d and m == 0:
      return is_n235(d)
    d, m = divmod(n, 3)
    if d and m == 0:
      return is_n235(d)
    d, m = divmod(n, 5)
    if d and m == 0:
      return is_n235(d)
    return False

is_n235 = IsN235WithMemo()

for i in range(10):
  print i, is_n235(i)

元になったコードへのコメントを確認してください。


DNUM=16 としてみました。25 まで正しく列挙されているように思いますが…… http://codepad.org/QQfKglM7


勘違いです。失礼しました。 8が先にqueueから抜かれますね。orz


8を取り出したときに, 16, 24, 40が入れられるので, 25の前に16出ますよ.

あとこのアルゴリズムの正当性も証明できます. n (> 1) 回目が終わったとき出てなかったものの最小値をAとします. このAがn回目に出た数より小さいと仮定します. A = 2^i 3^j 5^kで, queのことを考えると, A/2またはA/3またはA/5のどれかも出ていません. なので最小性に反します. よって矛盾.


出題時に考えていた答。平衡木を使うので、求める数の個数を n としたとき、時間計算量は O(log n)、空間計算量は O(n) のはず。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
(define (main args)
  (let ((tm (alist->tree-map '((1 . #t)) = <)))
    (let loop ((n (string->number (cadr args)))
               (rs '()))
      (cond
       ((zero? n)
        (print (reverse rs))
        0)
       (else
        (let ((m (car (tree-map-pop-min! tm))))
          (for-each (cut tree-map-put! tm <> #t)
                    (map (cute * m <>) '(2 3 5)))
          (loop (- n 1) (cons m rs))))))))

無限リストのつもりです。
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
class Hamming < Array
  def hamming(&p)
    x = 1
    loop {
      p.call x if 30**x%x == 0
      x += 1
    }
  end
end

p Enumerable::Enumerator.new(Hamming.new, :hamming).take(100)

求める数の個数 n なのに 時間計算量 O(log n) というのがよく判らないんです。
素人考えだと、最低 O(n)はかかるような気がするんですけど、解説をお願いできますでしょうか

VBA for Excel (2003)

 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
Sub main()
  Dim r, s, t
  Range("B1").Value = "2^x"
  Range("C1").Value = "3^y"
  Range("D1").Value = "5^z"
  Range("B2:D2").Value = 0
  Set r = Range("A1")
  n = 1
  c = 0
  While c < 100
    Set s = r.Offset(n)
    s.Value = n
    If Not IsEmpty(s.Offset(0, 1)) Then
      c = c + 1
      For i = 0 To 2
        Set t = r.Offset(n * ((i * (i + 1) / 2) + 2))
        For j = 0 To 2
          t.Offset(0, j + 1).Value = s.Offset(0, j + 1).Value + IIf(i = j, 1, 0)
        Next
      Next
    End If
    n = n + 1
  Wend
  Cells.AutoFilter Field:=1, Criteria1:="<>"
  Cells.AutoFilter Field:=2, Criteria1:="<>"
End Sub

確認しました。問題ないということでよろしいですね。

#7642 にマイナス評価が付いてるけど、 ジェネレイタにしなかった罰だと思うことにしますw


しらみつぶしに計算してからsortする方式ですが、その境界の目処はたてるようにしてみました。 http://en.wikipedia.org/wiki/Image:Regular_divisibility_lattice.svg の図を、3次元空間のある象限に存在する三角錐とみて、それを含む直方体の中の格子点を取り出します。最低でも、およそ5/6が無駄なデータなのですが、うまく三角錐の部分だけの格子点を巡回する順序を決めるアルゴリズムが思い浮かばなかったので、そのままにしてあります。

$ erlc スクリプトのファイル名

$ erl -noshell -s len main 100 -s init stop

として実行します。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
-module(len).
-export([main/1]).

main([N|_]) ->
        Limit = list_to_integer(atom_to_list(N)),
        Factor = math:pow(Limit * math:log(5) * math:log(2)/ (math:log(2) * math:log(3)), 1 / 3),
        L = lists:sort([round(math:pow(2, X)*math:pow(3,Y)*math:pow(5,Z)) ||
                 X <- lists:seq(0,round(1 + Factor * math:log(5) / math:log(2)) ),
                 Y <- lists:seq(0,round(1 + Factor * math:log(5) / math:log(3)) ),
                 Z <- lists:seq(0,round(1 + Factor))]),
        io:format("~p~n",[lists:sublist(L,Limit)]).

もりっとArrowで
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
module Main where

import Control.Arrow

f = Left >>> f' 2 >>> f' 3 >>> f' 5 >>> const False ||| const True                                                           
  where f' n = g n ||| Right

g n = loop ((snd &&& fst >>> app) &&& (snd >>> g'))
  where g' f m | n > m = Right ()
               | True  = case m `divMod` n of (x,0) -> f x; _ -> Left m

main = runKleisli func [1..]
  where func = arr (filter f) >>> arr (take 100) >>> Kleisli (mapM_ print)

有名なエラトステネスのふるいを変形して。Haskellはあんまりやったことがないから勝手が分かりませんが。

1
2
3
4
5
sieve (x:xs) | divide235 x = x:sieve xs
             | otherwise = sieve [a | a <- xs, a `mod` x /= 0]
    where divide235 x = x `mod` 2 == 0 || x `mod` 3 == 0 || x `mod` 5 == 0

take 100 $ [1..6] ++ sieve [7..]

あー、またログインし忘れ…orz


(雑談です)
この v=m/(m&-m) というのは私には逆立ちしても思いつきそうにないんですが、有名なテクニックなんですか?
いまだになぜ上手くいくのかよくわかってないけど、ビットを書き出してみると確かにそうなるぽいですが・・・教えてえらいひと

m&-mの部分は一応知られたコードのはずです。
下位ビットから見て最初に1になるビットを検出するというコードとして知ったと思います。(この数字の意味は今回はじめて気づきましたが^^;;)

-mって言うのはmと足したとき0になる数字ですので、2進表記した場合の下位ビットから見て最初に1が出て来る場所まではmと同じ、それ以上は反転となる数字です。
-m=~m+1でもありますね。
4ビットでの例:
 2=0010
-2=1110

 3=0011
-3=1101

ここでm&-mをとれば、下から見て一番最初の1のみが残ることがわかります。
証明:
m=1*A+2*B+4*C+8*Dとすると
~m=1*~A+2*~B+4*~C+8*~D

-m=~m+1=1*A'+2*B'+4*C'+8*D'とすると
ビット加算はsum=a^b,Carry=a&bなので
A'=~A^1  = A
B'=~B^(~A&1) = ~B^~A =B^C
C'=~C^(~B&(~A&1)) =~C^(~B&~A)=C^(B|A)
D'=~D^(~C&(~B&(~A&1))) = ~D^(~C&~B&~A)=D^(C|B|A)

m&-m=1*A''+2*B''+4*C''+8*D''

a&(a^b)=a&(a&~b|~a&b)=a&~bより

A''=A&A'=A
B''=B&B'=B&(B^A) = B&~A
C''=C&C'=C&(C^(B|A)) = C&~(B|A) = C&~B&~A
D''=D&D'=D&(D^(C|B|A))=D&~(C|B|A)=D&~C&~B&~A

下位ビットがすべて0のときのみビットがたつことがわかる。
~~証明ここまで

さて、2で割るということは2進数で考えれば、右シフトです。2で割った余りというのはm&1となるわけなので、最下位ビットが0なら2の倍数です。
ここから考えれば最下位ビットから見て0の数だけ2で割れるわけです。

つまり先ほど求めたm&-mというのは、設問で言う2^iとなっていることがわかると思います。

というのでどうでしょう?

解説ありがとうございます。何をやっているかがようやく理解できました。
(式変形はこれからじっくり解読しようと思います。。)

#7671のExcelの解が面白かったので、GNU Emacsに移植してみました。ネタ実装です。

できるだけプログラミング言語というよりエディターのマクロって感じの書きかたを目指しました。カラムの表現は面倒なのでS式にしちゃいましたが。

 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
(defun hamming (num)
  (interactive (list (read-minibuffer "Number: " "100")))
  (let ((n 1) (c 0))
    (switch-to-buffer "*Hamming*")
    (erase-buffer)
    (insert "(nil 0 0 0)")
    (beginning-of-buffer)
    (while (< c num)
      (unless (eolp)
        (save-excursion
          (let ((s (read (current-buffer))))
            (setq c (1+ c))
            (backward-kill-sexp)
            (prin1 (cons n (cdr s)) (current-buffer))
            (mapc (lambda (i)
                    (hamming::goto-line-force (* n i))
                    (if (eolp)
                        (prin1 (cons nil
                                     (mapcar (lambda (j)
                                               (+ (nth (/ (1+ j) 2) s)
                                                  (if (= i j) 1 0) ))
                                             '(2 3 5) ))
                               (current-buffer) )))
                  '(2 3 5) ))))
      (setq n (1+ n))
      (forward-line) ))
  (setq outline-regexp "([0-9]")
  (outline-mode)
  (hide-body) )

(defun hamming::goto-line-force (n)
  (let ((r (goto-line n)))
    (unless (bolp) (insert ?\n))
    (while (< 0 r)
      (insert ?\n)
      (setq r (1- r) ))))

1 2 3 5 のリストを掛け合わせ uniq して、を繰り返す。
最後にソートして n 個とりだす。
1
wd"0(100){./:~~.,1 2 3 5*/^:10[1

#7639を参考に移植。for文がないのがきつい・・・

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
limit=100
n=0
i=1
(n<limit)の間
    tmp=i
    (tmp%2=0)の間,tmp=tmp/2
    (tmp%3=0)の間,tmp=tmp/3
    (tmp%5=0)の間,tmp=tmp/5
    もし(tmp=1)ならば
        iを表示
        n=n+1
    i=i+1

SQLiteで。2,3,5の11乗まで計算し、その直積をとるだけ。
SQLiteで100行だけ取り出す方法がわからなかったので結果テーブルに一度格納してます。Oracleだとwhere rowno<=100が使えたような
 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
-- 数字テーブル
create table d (
  n   integer,
  p2  integer,
  p3  integer,
  p5  integer
);
insert into d values (0,1,1,1);
insert into d select max(n)+1, max(p2)*2, max(p3)*3, max(p5)*5 from d;
insert into d select max(n)+1, max(p2)*2, max(p3)*3, max(p5)*5 from d;
insert into d select max(n)+1, max(p2)*2, max(p3)*3, max(p5)*5 from d;
insert into d select max(n)+1, max(p2)*2, max(p3)*3, max(p5)*5 from d;
insert into d select max(n)+1, max(p2)*2, max(p3)*3, max(p5)*5 from d;
insert into d select max(n)+1, max(p2)*2, max(p3)*3, max(p5)*5 from d;
insert into d select max(n)+1, max(p2)*2, max(p3)*3, max(p5)*5 from d;
insert into d select max(n)+1, max(p2)*2, max(p3)*3, max(p5)*5 from d;
insert into d select max(n)+1, max(p2)*2, max(p3)*3, max(p5)*5 from d;
insert into d select max(n)+1, max(p2)*2, max(p3)*3, max(p5)*5 from d;
-- 結果テーブル
create table result (
  id  integer primary key,
  i   integer,
  j   integer,
  k   integer,
  val integer
);
insert into result (i, j, k, val)
  select a.n, b.n, c.n, a.p2 * b.p3 * c.p5 val from d a, d b, d c order by val;
select val,'2^'||i||' * 3^'||j||' * 5^'||k from result where id <= 100;

バージョン 2.7.3の頃に使ったのでちょっと自信がありませんが、先頭から 100件という
ことであれば、

  e.g.
    select ... from ... limit 100;

のように limit句で、Oracleのrownumと同等のことを実現できたと思います。

# 先頭の位置(オフセット)はoffset句で調節します。

 Streamを使って書いてみました。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
object HammingNumbers {
    def from(n:Int):Stream[Int] = Stream.cons(n,from(n+1))
    def hammings(s:Stream[Int]):Stream[Int] = {
        def divide(n:Int,b:Int):Int = (n % b == 0) match {
                case false => n
                case _ => divide(n/b,b)
            }
        Stream.cons(s.head,hammings(s.tail.filter { n => List(2,3,5).foldLeft(n) { (v,b) => divide(v,b) } == 1 }))
    }
    def hammings():Stream[Int] = hammings(from(1))
    def main(args:Array[String]):Unit = {
        try {
            val    n:Int = args.length match {
                case 1 => args(0).toInt
                case _ => 100
            }
            println(hammings().take(n).mkString("\n"))
        } catch {
            case e => e.printStackTrace
        }
    }
}

 #7638と#7671のやり方で。

 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
import    scala.collection.immutable.SortedSet
import    scala.collection.immutable.TreeSet
abstract class CHammingNumbers {
    def take(n:Int):List[Int]
}
class CHammingNumbersG extends CHammingNumbers {
    def next(s:List[Int],h:List[Int]):Tuple2[List[Int],List[Int]] = {
        val    m:List[Int] = s.zip(List(2,3,5)).map { v => h.apply(v._1)*v._2 }
        val    n:Int = m.sort { (a,b) => a < b }.head
        (s.zip(m).map { v => (v._2 == n) match { case true => v._1 + 1; case _ => v._1 } },h+n)
    }
    def take(n:Int,s:List[Int],h:List[Int]):Tuple2[List[Int],List[Int]] = n match {
            case 0 => (s,h)
            case _ => next(s,h) match { case v => take(n-1,v._1,v._2) }
        }
    def take(n:Int):List[Int] = take(n-1,List(0,0,0),List(1))._2
}
class CHammingNumbersS extends CHammingNumbers {
    def next(c:SortedSet[Int]):Tuple2[Int,SortedSet[Int]] = (c.firstKey,(TreeSet(c.firstKey*2,c.firstKey*3,c.firstKey*5)++(c-c.firstKey)))
    def take(n:Int,c:SortedSet[Int]):List[Int] = n match {
            case 0 => List()
            case _ => next(c) match { case v => v._1::take(n-1,v._2) }
        }
    def take(n:Int):List[Int] = take(n,TreeSet(1))
}
object HammingNumbers {
    def main(args:Array[String]):Unit = {
        try {
            val    n:Int = args.length match {
                case 1 => args(0).toInt
                case _ => 100
            }
            List(new CHammingNumbersG,new CHammingNumbersS).foreach { h => println(h.take(n).mkString("\n")) }
        } catch {
            case e => e.printStackTrace
        }
    }
}

グッドです。以下のように書いて結果テーブルは不要になりました。
バージョン 2.8.17 と 3.6.1 で動くことを確認してあります。
1
select (a.p2*b.p3*c.p5) val, a.n, b.n, c.n from d a, d b, d c order by val limit 100;

sedだと割り算のコストが高すぎるので、他のアルゴリズムでは難しそうです。
乗算+ソートのために内部は2進で計算しています。
 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
#!/bin/sed
s/.*/lol/
: loop
    s/^\([^l]*\)l*o\([lo]*\)\(\n.*\)\?$/\1a\2\3\nx\2o\nx\2o+\2=\nx\2oo+\2=/
    : 2->10
        s/[f-j]/#&/g
        s/[a-j][#l]/\U&/g
        y/AbBcCdDeEfFgGhHiIjJ/bcdefghijabcdefghij/
        s/^#/b/m
        s/#//g
        s/\([a-j]\)[Lo]/\1/
    /[a-j][lo]/ b 2->10
    y/abcdefghij/0123456789/
    : add
        s/\([xo]l*\)l+\([lo]*\)l=/\U\1\E+\2=o/g
        s/X/xl/g
        y/LO/ol/
        s/l+\([lo]*\)o=/+\1=l/g
        s/o+\([lo]*\)\([lo]\)=/+\1=\2/g
        s/+=//g
    t add
    s/^x\([lo]*\)$/\U\1\Eo\1/gm
    y/LO/ll/
    : sort
        s/^\([lo]*\)\n\1$/\1/gm
        s/^\(\([lo]*\)l[lo]*\)\n\(\2o[lo]*\)$/\3\n\1/gm
    t sort
    s/^\(\(\w*\n\)\{99\}\w*\)\n.*$/\1/
/[lo]/ b loop

因数分解する方針で。

 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
import std.stdio;

struct HammingNumbers {
    static bool isHammingNumber(uint n) {
        static bool[uint] memo;
        
        if(n == 1) return true;
        if(auto p = n in memo) return *p;
        
        if(n % 2 == 0)
            return memo[n] = isHammingNumber(n / 2);
        if(n % 3 == 0)
            return memo[n] = isHammingNumber(n / 3);
        if(n % 5 == 0)
            return memo[n] = isHammingNumber(n / 5);
        return memo[n] = false;
    }
    
    static int opApply(int delegate(ref const(size_t), ref const(uint)) dg) {
        size_t i = 0;
        uint n = 1;
        while(true) {
            if(isHammingNumber(n)) {
                if(auto result = dg(i, n))
                    return result;
                i++;
            }
            n++;
        }
    }
}

void main(){
    foreach(i, n; HammingNumbers) {
        if(i == 100) break;
        writeln(i, ": ", n);
    }
}

itertoolsを使って。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
import itertools

def numbers():
    for n in itertools.count(1):
        t = n
        for i in (2, 3, 5):
            while t != 1 and t % i == 0:
                t //= i
        if t == 1:
            yield n

def main():
    for n in itertools.islice(numbers(), 100):
        print n

if __name__ == '__main__':
    main()

2,3,5以外の素数の組合せ(3,7,11)にも対応。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
def hamming(n ,a = [1 ,2 ,3 ,5])
  r = []
  h = {}
  a.each{|x|h[x] = true}
  1.upto(n){
    r << cnt = h.keys.sort[0]
    a.each{|x|h[x * cnt] = true}
    h.delete(cnt)
  }
  r
end

p "n = "
p hamming(gets.chomp.to_i)

Forth の実績追加のため、実装してみました。

 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
: DIV_CHECK ( n d -- n/d^m )
BEGIN
  2DUP MOD 0=
WHILE
  DUP -ROT / SWAP
REPEAT
DROP ;

: HummingNumbers ( n -- )
1
BEGIN
  OVER 0>
WHILE
  DUP
  2 DIV_CHECK
  3 DIV_CHECK
  5 DIV_CHECK
  1 = IF
    DUP .
    SWAP 1- SWAP
  THEN
  1+
REPEAT
2DROP
;

オーダーは考えていませんが。

1
2
import List
main = print $ take 100 $ sort $ [2^i*3^j*5^k|lm<-[0..99],i<-[0..lm],j<-[0..lm-i],let k = lm-i-j]

SICP Prloblem 3.56を参考にストリームで。 ただし、OCamlではvalue recursionは容易ではないので、ストリーム構築時に一部黒魔術(Objモジュール)を使っています。

時間効率性はO(n) (のはず)

 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
type 'a t = Nil | Cons of 'a * 'a t lazy_t

let (!) = Lazy.force

let rec ( ** ) factor = function
    | Nil -> Nil
    | Cons (a, b) -> Cons (a * factor, lazy (factor ** !b))

let rec ( ++ ) s1 s2 = match s1, s2 with
    | Nil,          _            -> s2
    | _,            Nil          -> s1
    | Cons(h1, t1), Cons(h2, t2) ->
        if h1 < h2 then Cons(h1, lazy (!t1 ++  s2)) else
        if h1 > h2 then Cons(h2, lazy ( s1 ++ !t2)) else
                        Cons(h1, lazy (!t1 ++ !t2))

let hamming =
    let s = Cons(1, lazy Nil) in
    let _ = Obj.set_field (Obj.repr s) 1
        (Obj.repr (lazy (2**s ++ 3**s ++ 5**s) )) in
    s

let rec print n s = if n > 0 then match s with
    | Nil -> ()
    | Cons (a, b) -> Printf.printf "%d\n" a; print (n-1) !b

let _ = print 100 hamming

C++でO(n)のがなさげだったので.

 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
#include <iostream>
#include <list>

int main(int,char**)
{
  std::list<int> q;
  q.push_back(1);
  std::cout << "1\n";

  std::list<int>::iterator i2,i3,i5;
  i2 = i3 = i5 = q.begin();

  for(int i=0;i<100-1;++i)
  {
    int n = std::min(2**i2,std::min(3**i3,5**i5));
    q.push_back(n);
    if(2**i2==n)++i2;
    if(3**i3==n)++i3;
    if(5**i5==n)++i5;

    if(q.front()<std::min(*i2,std::min(*i3,*i5)))
      q.pop_front();

    std::cout << n << std::endl;
  }
  return 0;
}

素因数分解用の関数があればもっと楽になるかな

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
def 素因数分解(num){
    def map = [:]
    for( int i = 2; i <= num ;i++ ){
        if( num%i == 0 ){
            map[i] = (map[i]?:0) + 1
            num = (int)(num/i)
            i--
        }
    }
    map
}

def list = [:]
def oklist = [2, 3, 5]
for(def i = 1;list.size()<100;i++){
    def map = 素因数分解(i)
    if(map.keySet() - oklist)
        continue
    list[i] = map
}
list.each{ key, value ->
    println("${key}".padLeft(4) + " = " + oklist.collect{ "${it}^${value[it]?:0}" }.join(" * "))
}

関数名に困ったので、こういう数の集合を「アンダロート数」と呼ぶことにします。
 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
51
52
53
54
55
56
57
58
59
60
61
62
63
64
#module

#define true  1
#define false 0

// 可能な限り割り続ける
#defcfunc div_force int p1, int p2, var count, local val
    val   = p1
    count = 0
    
    repeat
        if ( (val \ p2) == 0 ) {
            val /= p2
            count ++
        } else {
            break
        }
    loop
    
    return val
    
// アンダロート数の集合かどうか
#defcfunc IsSetOfAndarote int p1, array condition, array result, local num, local count
    num = p1
    dim result, length(condition)
    foreach condition
        num = div_force(num, condition(cnt), result(cnt))
    loop
    return ( num == 1 )
    
#global

#define MAX_COUNT 100

*main
    dim condition, 3
        condition = 2, 3, 5
    dim result, 3
    sdim buf, 32000
    
    buf = "集合 { "
    foreach condition
        if ( cnt != 0 ) { buf += ", " }
        buf += condition(cnt)
    loop
    buf += " } に対するアンダロート数を"+ MAX_COUNT +"個列挙します。\n----------\n"
    
    count = 0
    repeat , 1
        if ( IsSetOfAndarote(cnt, condition, result) ) {
            buf += strf("%5d = ", cnt)
            foreach condition
                buf += "("+ condition(cnt) +" ^"+ strf("%2d", result(cnt)) +")"
            loop
            buf += "\n"
            count ++
            if ( count == MAX_COUNT ) { break }
        }
        await 0
    loop
    
    objmode 2
    mesbox buf, ginfo(12), ginfo(13)
    stop

十進BASICで。
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
LET n = 100
DO WHILE(n > 0)
   LET i = i + 1
   LET t = i
   LET h = 2
   DO WHILE(h < 6)
      DO WHILE(MOD(t,h) = 0)
         LET t = t / h
      LOOP
      LET h = h * 2 - 1
   LOOP
   IF t = 1 THEN
      PRINT i
      LET n = n - 1
   END IF
LOOP
END

Index

Feed

Other

Link

Pathtraq

loading...