ICPC Probrem C 多倍長で解いてみた。

　一昨日の問題を多倍長で解いてみました。
　
　問題はここから。ここから先は解説なので、読んでみて自分で考えてみようと思う人はURLより下を見ないように。
http://www.logos.ic.i.u-tokyo.ac.jp/acm-icpc/C.jp/C.html
　
　考え方は簡単で、まず1つ目の単位分数を1/1から順番に1/2,1/3,1/4……を前項までの合計値に対して加える。加えた値に対して与えられたp/qと大小を判断する。ここで、
　

• p/q=合計値となった場合、分数が一致したということで正解とする。
• p/q<合計値となった場合、探索した数字では大きすぎるので加算する単位分数の分母を大きくして、より小さい数字にしてから再びp/qと大小の判断を行う。ただし、分母の積が与えられたaよりも大きくなった場合、ループを脱出して再帰を脱出する。
• p/q>合計値となった場合、更に単位分数を加算してみて同様に大小を判断する。ただし、次に加算する単位分数は分母が前項に使用した分母以上の数字を用いる。また、単位分数の数が与えられたnより大きくなった場合、再帰を脱出する。

　
という条件で、再帰関数を書く。ここまでは試合中までに書けた。ただし、その場合以下の問題が発生した。
　

• p/qと単位分数同士の加算の比較時に丸め誤差が大量に発生して、右辺と左辺の一致が取れない。試合中は、適当に整数倍したものをint型に変換するなど、かなり無茶なことをやっていた。

　
　そもそもの話をすると、私が多倍長プログラムを書くのが面倒でやらなかっただけなのだが、それが大きく裏目に出てしまったようだ。そこで、家に帰ってから（その後かなり自己嫌悪に陥っていましたが）、多倍長プログラムに書き直した。それが下のプログラムである。
　
　桁数の精度は4*25=100桁となっている。配列の先頭は整数部分なので、小数第96位まで出ることになる。

#include 

#include 

#include 

#include 

#include 
#include 

#include 
#pragma warning(disable:4786)
using namespace std;
class FileTokenizer

{

public:

　　typedef vector STRARRAY;

　　typedef vector STR2DARRAY;

　　FileTokenizer(const char* psz, const char* sep=" ")

　　　　: doc_()

　　{

　　　　fstream fin(psz, ios::in);

　　　　if(fin.fail())

　　　　　　throw std::runtime_error*1;

　　　　while(!fin.eof())

　　　　{

　　　　　　static char buf[256];

　　　　　　fin.getline(buf, 256);

　　　　　　STRARRAY vct;
　　　　　　char* tp;

　　　　　　tp = strtok(buf, sep);

　　　　　　if(tp != NULL)

　　　　　　　　vct.push_back(tp);

　　　　　　while(tp != NULL)

　　　　　　{

　　　　　　　　tp = strtok(NULL, sep);

　　　　　　　　if(tp != NULL)

　　　　　　　　　　vct.push_back(tp);

　　　　　　}

　　　　　　doc_.push_back(vct);

　　　　}

　　　　fin.close();

　　}

　　STR2DARRAY& get()

　　{

　　　　return doc_;

　　}

private:

　　STR2DARRAY doc_;

};
int ans;
// 初期化

void init_array(int array,int len){

　　for(int i=0;i
　　　　array[i]=0;

　　}

}
// 多倍長分数の生成

void init_taniarray(int array,int undiv,int div,int len){

　　int d;

　　int taihi;

　　for(int i=0;i
　　　　array[i]=0;

　　}

　　array[0]=undiv;

　　d=0;

　　for(i=0;i
　　　　taihi=(array[i]+d*10000)/div;

　　　　d=(array[i]+d*10000)%div;

　　　　array[i]=taihi;

　　}

}
// 多倍長演算同士の足し算

void addarray(int array0,int array1,int array2,int len){

　　int d;

　　d=0;

　　array2[24]=array2[24]+1;

　　for(int i=len-1;i>=0;i--){

　　　　array0[i]=array1[i]+array2[i]+d;

　　　　if(array0[i]>=10000){

　　　　　　d=1;

　　　　　　array0[i]=array0[i]%10000;

　　　　}else{

　　　　　　d=0;

　　　　}

　　}

}
// 多倍長演算同士の評価

// array1が大きい 1

// array2が大きい -1

// 同じ 0

int arraycmp(int array1,int array2,int len){

　　for(int i=0;i
　　　　if(array1[i]
　　　　if(array1[i]>array2[i]) return 1;

　　}

　　return 0;

}
// 再帰による解の発見

// pdivq: p/q

// a: 超えてはならない分母の積

// xi: 現時点での分母の積

// n: 分数の項数

// i: 探索を開始する数字

// total: 現時点での右辺の分数の合計値
void recursive(int pdivq,int a,int xi,int n,int i,int total[]){

　　int n_total[25];

　　int tanni[25];

　　if(n==0) return;

　　init_array(n_total,25);

　　init_array(tanni,25);

　　for(int j=i;xi*j<=a;j++){

　　　　// 演算

　　　　init_taniarray(tanni,1,j,25);

　　　　addarray(n_total,total,tanni,25);

　　　　if(arraycmp(n_total,pdivq,24)==0){

　　　　　　// printf("解の発見\n");

　　　　　　ans++;

　　　　}else if(arraycmp(n_total,pdivq,24)==1) continue;

　　　　else if(arraycmp(n_total,pdivq,24)==-1){

　　　　　　// printf("更に奥への探索\n");

　　　　　　recursive(pdivq,a,xi*j,n-1,j,n_total);

　　　　}

　　}

}
int main(void){

　　int p,q,a,n;

　　int initvalue[25];

　　int dummy[25];

　　FileTokenizer file("in.txt");

　　FileTokenizer::STR2DARRAY array=file.get();

　　init_array(dummy,25);

　　// p,q,a,n 一行目の

　　p=atoi(array.at(0).at(0).c_str());

　　q=atoi(array.at(0).at(1).c_str());

　　a=atoi(array.at(0).at(2).c_str());

　　n=atoi(array.at(0).at(3).c_str());

　　for(int i=1;!(p==0 && q==0 && a==0 && n==0);i++){

　　　　ans=0;

　　　　init_taniarray(initvalue,p,q,25);

　　　　recursive(initvalue,a,1,n,1,dummy);

　　　　printf("%d\n",ans);

　　　　// 次のp,q,a,n

　　　　p=atoi(array.at(i).at(0).c_str());

　　　　q=atoi(array.at(i).at(1).c_str());

　　　　a=atoi(array.at(i).at(2).c_str());

　　　　n=atoi(array.at(i).at(3).c_str());

　　}

　　return 0;

}

　これならば、かなりの精度で比較が出来る。しかし、このプログラムを作って演算させた結果、サンプルと2点違う答えが出てきた。（明らかに少ない）
　
　更に検討した結果、これでも不備があることが分かった。
　例えば、(p,q,a,n)=(2 4 54 2)という値のときに、答えとしては、以下のようなパターンが想定される。

　2/4=1/4+1/4

　2/4=1/4+1/8+1/8

　2/4=1/3+1/6

　上二つの答えは普通に検出できる。
　しかし、一番下の式に対しては多倍長演算を用いた場合でも、

　左辺　1/2=0.5

　右辺　0.33333333333（略）+0.1666666666（略）=0.499999999999999（略）

　と、これでは一致が取れない。
　そこで、少々反則技であるが、加算時に、小数第96位に1を加えた状態で加算する。ただし、比較時には多少精度を緩めないと、0.00000（中略）0001分の誤差が生じるため、小数第80位までの比較に留めておく。これによって、丸め誤差を回避することが出来る。このプログラムによって、入力データに対する出力結果は………
　
　長いのでやめておきましょう。＿|￣|○
　
　私もこの答えがあっているかどうかはわかりませんが、サンプル解答では正解が出ているので、多分大丈夫でしょう。
　しかし、Pentium3-500MHzのマシンでプログラムを動かしてみたら、妙に遅かった……普段はPen4-3.06GHzのマシンを使っているので、全然気がつかなかったのですが、やっぱり普段から速いマシンを使ってると良くないですねぇ……。
　
　初めから多倍長で解けば行けたかもしれない、と考えてみたが、果たして試合中の緊迫した空気の中で、後者のバグが見つけられたか？はちょっと微妙なところ。素直に他の問題を解いていれば良かったかも。……と、考えたらちょっと気分が楽になった。

*1:std::string("File Not Found : ")+psz).c_str(