The Balloon Popping Problem 2

The Balloon Popping Problemの続き．

昔書いたHaskellのコードを解読したところ，
風船が少ない時について，2通りの戦略で空気を入れたときの最大期待値を求めるものであった．
風船が少ない時についての最適値を求める訳ではない．
ただし，DPにすればいいのになぜか再帰だけで書いていたのでjavaで書き直してみた．

import static java.lang.Math.*;

public class Balloon {
    public static void main(String[] args) {
        int num = Integer.parseInt(args[0]);
        Rational res = new Balloon().solve(num);
        System.out.println(res);
        System.out.println(res.doubleValue());
        for(int i=0;i<num;++i){
            res = res.subtract(new Rational("1",""+(i+1)*(i+1)));
        }
        System.out.println(res);
        System.out.println(res.doubleValue()+(PI*PI/6));
    }
    Rational solve(int num) {
        int size = 1 << num;
        Rational[] rs = new Rational[num+1];
        Rational[] vols = new Rational[num];
        for(int i=0;i<=num;++i)rs[i] = new Rational(""+i);
        for(int i=0;i<num;++i)vols[i] = new Rational("1",""+(i+1));
        Rational[] exp = new Rational[size];
        for(int i=1;i<size;++i){
            int bals = 0;
            int min = 0;
            int max = num-1;
            for(int j=0;j<num;++j)if((i&(1<<j))>0){
                    bals++;
                    if(vols[min].compareTo(vols[j])>0)min=j;
                    if(vols[max].compareTo(vols[j])<0)max=j;
            }
            if(bals==1){
                exp[i]=vols[min];
            }
            else{
                Rational negative = vols[min].multiply(rs[bals]);
                Rational positive = vols[max];
                for(int j=0;j<num;++j)if((i&(1<<j))>0){
                        positive = positive.add(exp[i-(1<<j)]);
                }
                positive = positive.divide(rs[bals]);
                exp[i] = negative.max(positive);
            }
        }
        return exp[size-1];
    }
}

ただし，Rationalはこれ．

実行してみたところ4分ぐらいで20個の場合が解けて期待値が

9620098467715313103701957/5899598819922718310400000

であり，nが十分大きい場合には

1.6794

が下界であることを示した．
2^20程度なのに遅いのは有理数クラスが遅い所為．

計算誤差が怖いが高速化のためdoubleで組み直してみた．
ついでにC++に変更．

#include <iostream>
#include <vector>
#include <cmath>
#include <cstdlib>
using namespace std;

int main(int argc, char *argv[]){
  int num = atoi(argv[1]);
  int size = 1<<num;
  vector<double> vols;
  for(int i=num;i>0;--i)vols.push_back(1.0/i);
  vector<double>  exp(size);
  for(int i=1;i<size;++i){
    int bals = 0;
    int argmin = num-1;
    int argmax = 0;
    for(int j=0;j<num;++j)if((i&(1<<j))>0){
        bals++;
        if(vols[argmin]>vols[j])argmin=j;
        if(vols[argmax]<vols[j])argmax=j;
    }
    if(bals==1){
      exp[i]=vols[argmin];
    }
    else{
      double negative = vols[argmin]*bals;
      double positive = 0;
      for(int j=0;j<num;++j)if((i&(1<<j))>0){
          positive = positive+exp[i-(1<<j)];
      }
      positive += vols[argmax];
      positive = positive/bals;
      exp[i] = max(negative,positive);
    }
  }
  cout << exp[size-1] << endl;
  double tmp = 0.0;
  for(int i=num;i>0;--i)tmp+=1.0/(i*i);
  cout << exp[size-1]-tmp << endl;
  cout << exp[size-1]-tmp+(M_PI*M_PI/6) << endl;
}

このプログラムだと28個の場合を2分ほどで解き，期待値は

1.64549

nが十分大きいときの下界は

1.68057

となった．
これ以上のサイズではメモリが足りない．

論文の下界は1.659なのでなかなか良い改善だと思う．