2014年11月16日日曜日

「Genuine」なエラトステネスのふるいで無限素数リストを作る仕組み

Haskell の Data.Numbers.Primes の 素数の無限リストってどうやって作ってるのだろうと思って ソースを読んでたら、インスパイアされたという論文『The Genuine Sieve of Eratosthenes』へのリンクがソースコメントに記載してあった。

Haskell 遅延評価の例としてよく見かける、
primes = sieve [2..]
sieve (p : xs) = p : sieve [x | x <− xs, x ‘mod‘ p > 0]
というアルゴリズムとは異なる「真のエラトステネスのふるい」を、関数型らしく遅延評価と無限リストで書いてみるという趣向で、読み応えがある。

なにしろ無限リストなので、オリジナルの「ふるい」のように上限を決めて合成数をキャンセルしていくわけにはいかないので、素数候補を順次チェックしていくために必要な分だけ合成数を保持しておく必要がある。

そのコードを引用したものが以下のもの(元の論文では、PriorityQueue を使った発展形がこれの後に記述されているが、簡単なのでこっちで考えてみる)
sieve xs = sieve’ xs Map.empty
  where
    sieve’ [] table = []
    sieve’ (x:xs) table =
      case Map.lookup x table of
          Nothing −> x : sieve’ xs (Map.insert (x*x) [x] table)
          Just facts −> sieve’ xs (foldl reinsert (Map.delete x table) facts)
        where
          reinsert table prime = Map.insertWith (++) (x+prime) [prime] table        
このコードをもとに sieve' (x:xs) table が再帰するにつれて、どんな感じで素数が得られたり(あるいは得られなかったり)、合成数のテーブルが更新されていくのか、ワンステップずつ様子が観察できるような Javascript コードを書いてみた。

        (function(prime01, $, undefined) {
        var map = null;
        var wheel = null;
        var wheelIndex = 0;
        var currentCandidate = null;
        var nextCandidate = null;

        var $table = null;

        function nextStep() {
          var result = wheel[wheelIndex];
          wheelIndex = rotateWheel(wheelIndex);
          return result;
        }
        function rotateWheel(wheelIndex_) {
          ++wheelIndex_;
          return (wheelIndex_ == wheel.length) ? 0: wheelIndex_;
        }
        prime01.shiftCandidate = function () {
          var found = search(currentCandidate.p, 0, map.length-1);
          var $lastRow = $table.find('tr.calc-step').last().find('td').first();
          if (found || found==0) {
            var elem = map[found];
            map.splice(found, 1);
            _.each(elem.v, function(prime) {
              insert(elem.k + prime.p * wheel[prime.w], prime, 0, map.length-1);
              prime.w = rotateWheel(prime.w);
              $lastRow.css('color', '#FF0000');
            });
          } else {
            insertSqr(map, currentCandidate);
            appendPrimeToDiv(currentCandidate.p);
            $lastRow.css('color', '#00FF00');
          }
          currentCandidate = nextCandidate;
          nextCandidate = { p: nextCandidate.p+nextStep(), w:wheelIndex};

          addStepRow();
        }
        function appendPrimeToDiv(prime) {
          var $primeList = $('#prime-list');
          var currentList = $primeList.text();
          $primeList.text(currentList + (currentList=="" ? "" : ", ") + prime);
        }
        function addStepRow() {
          updateRow(createStepRow().appendTo($table));
        }
        function keyAt(index) {
          var elem = map[index];
          return elem ? elem.k : null;
        }
        function search(num, min, max) {
          if      (max < min)         return null;
          else if (keyAt(min) == num) return min;
          else if (keyAt(max) == num) return max;

          var mid    = Math.floor((min+max)/2);
          var midKey = keyAt(mid);

          if      (midKey == num) return mid;
          else if (midKey <  num) return search(num, mid+1, max-1);
          else                    return search(num, min+1, mid-1);
        }
        function insertSqr(map, num) {
          insert(num.p * num.p, num, 0, map.length-1);
        }
        function insertAt(index, element) { map.splice(index, 0, element); }
        function pushAt(index, value) { map[index].v.push(value); }
        function insert(key, value, min, max) {
          var e = {k:key, v:[value]};
          var insertInBetween = function(min_, max_) {
            if      (keyAt(min_) > key) return insertAt(min_,   e);
            else if (keyAt(max_) < key) return insertAt(max_+1, e);
            else                        return insertAt(max_,   e);
          };
          if      (map.length==0)    map.push({k:key, v:[value]});
          else if (keyAt(min)== key) pushAt(min, value);
          else if (keyAt(max)== key) pushAt(max, value);
          else if (keyAt(max) < key) insertAt(max+1, e);
          else if (keyAt(min) > key) insertAt(min,   e);
          else if (min+1 == max) insertAt(max, e);
          else {
            var mid = Math.floor((min+max)/2);
            if        (keyAt(mid) == key) pushAt(mid, value);
            else   if (keyAt(mid) <  key) {
              if (mid+1==max) insertInBetween(mid, max);
              else insert(key, value, mid+1, max-1);
            } else if (keyAt(mid) > key) {
              if (min+1==mid) insertInBetween(min, mid);
              else insert(key, value, min+1, mid-1);
            }
          }
        }
        function mapToString2(map) {
          return _.reduce(map, function(acc, e) {
               var facts = _.map(e.v, function(a) { return a.p; } ).join();
              return acc+(acc=="" ? "": ", ") + e.k +":[" + facts+"]";
          }, "");
        }
        function updateRow(row) {
          row.find('td')
            .first().text(currentCandidate.p)
            .next().text(currentCandidate.p + "+" + wheel[currentCandidate.w] + ",,, ")
            .next().text(mapToString2(map))

        }
          function initializeWheel(primes_, current_, wheel_) {
            map = [];
            wheel = wheel_;
            wheelIndex = 0;
            currentCandidate = {p:current_, w:0};
            $('#prime-list').text("");
            _.each(primes_, function(p) { appendPrimeToDiv(p); });
            nextCandidate = {p:currentCandidate.p + nextStep(), w:wheelIndex};
            updateRow($table.find('tr'));
          }
          var wheelSetups = {
            0: function() {initializeWheel([],      2, [1]);},
            1: function() {initializeWheel([2],     3, [2]);},
            2: function() {initializeWheel([2,3],   5, [2,4]);},
            3: function() {initializeWheel([2,3,5], 7, [4,2,4,2,4,6,2,6]);},
          };
          prime01.prepareWheel = function (wheelNumber) {
            wheelSetups[wheelNumber]();
          }
          function createStepRow() {
            return $('<tr class="calc-step">')
              .append('<td style="text-align:right;padding:1px 5px 1px 5px"/>')
              .append('<td style="padding:1px 5px 1px 5px;text-align:center;"/>')
              .append('<td/></tr>');
          }
          prime01.initializeWheelSelection = function () {
            $table.find('tr.calc-step').remove();
            $table.append(createStepRow());
            $('input:radio.wheel-selection-radio:checked').removeAttr('checked');
            $('input:radio[value=1]').prop('checked', true).trigger('change');
          }
          prime01.initialize = function() {
            $table = $('table#calculation-steps');
            $('input:radio.wheel-selection-radio').change(function(){
              prime01.prepareWheel(parseInt($(this).attr('value')));
              $('button#next').removeAttr('disabled');
            });
            $("button#next").click(function(){
              $('input.wheel-selection-radio').attr('disabled', 'disabled');
              prime01.shiftCandidate();
            });
            $("button#clear").click(function(){
              $('input.wheel-selection-radio').removeAttr('disabled');
              prime01.initializeWheelSelection();
            });
            prime01.initializeWheelSelection();
          }
        } (window.prime01 = window.prime01 || {}, jQuery));

        $(function() {
          prime01.initialize();
        });
sieve' (x:xs) table
x xs table
wheel 0: [1]
wheel 1: [2]
wheel 2: [2, 4]
wheel 3: [4,2,4,2,4,6,2,6]
PRIMES


next ボタンを押せば、sieve' の再帰が進み、それぞれの「ふるいにかけられる素数候補」と「合成数のテーブル」を示す行が、表の末尾に追加されていく。clear ボタンを押すと行が削除され初期状態にもどる。素数が得られたら、ボタンの上のリストに追加されていく。

またホイールも選べるようにしてみた。

wheel 0 はホイール無しで、2から上の自然数を全てふるいにかける。最初はこれを選択して観察すると仕組みが分かりやすい。デフォルトは wheel 1 にした。これは、3以上の奇数をふるいにかけるもの。wheel 2, 3 の意味も、実際に選択して next を何回か押下してみればわかると思う。

wheel n の n が増えると、車輪の円周が大きくなっていき、予め除外される合成数の割合が増えてチェックする回数が減っていくが、増えすぎてもメモリを圧迫して効果が悪くなるらしい。Haskell の Data.Numbers.Primes では 6 くらいが調度良いという。

仕組みがわかれば、Haskell 以外の遅延評価がデフォルトじゃなかったり、何なら手続き型の'不純な'言語でも、同じような感じで素数の無限リストを生成するコードが書けると思う。

0 件のコメント:

コメントを投稿