Entries

スポンサーサイト (この記事を編集する[管理者用])

上記の広告は1ヶ月以上更新のないブログに表示されています。
新しい記事を書く事で広告が消せます。

コメント

コメントの投稿

コメントの投稿
管理者にだけ表示を許可する

トラックバック

トラックバック URL
http://rsujskf.blog32.fc2.com/tb.php/1895-f284b367
この記事にトラックバックする(FC2ブログユーザー)

February 2011 Challenge 1問目 - Bogosort (この記事を編集する[管理者用])

Source

codechef February 2011 Challenge 1問目
Problem Statement
February 2011 Challengeの参加記録

問題概要

長さNの配列で,最初1~Nの置換で,k番目の要素はkではないものが与えられる.
それを以下のアルゴリズムでソートする場合,必要ステップ数 (以下の2が実行される回数) の期待値を分数で厳密に求める問題.
 1. 考えている区間 (未ソートの部分) を[1, N]とする.
 2. 考えている区間の要素を全てランダムに並び替える
 3. 考えている区間の左端の要素が,考えている区間内の最小値である限り,考えている区間の左端を+1する
 4. 考えている区間の右端の要素が,考えている区間内の最大値である限り,考えている区間の右端を-1する
 5. 考えている区間が空なら終了.そうでなければ2~4を繰り返す.
テストケースの数は150以下,Nは2以上150以下.ソースコード制限は10000バイト以下.

解法

多倍長を使って適当にDPする.
DPで工夫するか,多倍長の計算で工夫するかしないとTLEになる.
以下は素朴なDPだが,多倍長の計算で約分をできるだけ排除し,さらに約分するときのあまり大きい約数で分子分母ともに割れることがないのを利用して高速化している.

C言語のスパゲッティなコード
#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#define REP(i,a,b) for(i=a;i<b;i++)
#define rep(i,n) REP(i,0,n)

int getPrime(int n,int p[]){int i,j,n2=n/2;rep(i,n2)p[i]=1;for(i=3;i*i<n;i+=2)if(p[i>>1])for(j=(i*i)>>1;j<n2;j+=i)p[j]=0;j=1;p[0]=2;REP(i,1,n2)if(p[i])p[j++]=i*2+1;return j;}
int ps, p[1200];

#define ll long long

#define BIG_INT_SIZE 80
#define BIG_INT_BASE 100000000LL
#define BIG_INT_DIGITS 8
#define BIG_INT_CHAR_SIZE 4000 /* BIG_INT_SIZE*BIG_INT_DIGITS+2以上 */

typedef struct rational_number{ll a,b;}rational;
typedef struct rational_struct_point{rational x,y;}rpnt;
typedef struct big_integer{ll a[BIG_INT_SIZE];}bigInt;
typedef struct big_rational{bigInt a,b;}bigRational;
typedef struct big_rational_struct_point{bigRational x,y;}brpnt;

bigInt bigIntGlobalTmp;

int bigIntToChar(bigInt a,char ret[]);
void printBigInt(bigInt a);
void putBigInt(bigInt a);


ll llGCD(ll a,ll b){
  if(!a) return b; return llGCD(b%a,a);
}

bigInt bigIntZero(){
  bigInt a; int i;
  rep(i,BIG_INT_SIZE) a.a[i]=0;
  return a;
}

bigInt bigIntOne(){
  bigInt a; int i;
  REP(i,1,BIG_INT_SIZE) a.a[i]=0; a.a[0]=1;
  return a;
}

bigInt bigIntOrder(bigInt a){
  int i; ll k;
  REP(i,1,BIG_INT_SIZE) if(a.a[i-1]<0 || a.a[i-1]>=BIG_INT_BASE){
    k=a.a[i-1]/BIG_INT_BASE; a.a[i-1]-=k*BIG_INT_BASE;
    if(a.a[i-1]<0) k--, a.a[i-1]+=BIG_INT_BASE; a.a[i]+=k;
  }
  if(a.a[BIG_INT_SIZE-1]<0){
    rep(i,BIG_INT_SIZE) a.a[i]=-a.a[i]; a=bigIntOrder(a);
    rep(i,BIG_INT_SIZE) a.a[i]=-a.a[i];
  }
  return a;
}

void bigIntOrderPt(bigInt *a){
  int i; ll k;
  REP(i,1,BIG_INT_SIZE) if(a->a[i-1]<0 || a->a[i-1]>=BIG_INT_BASE){
    k=a->a[i-1]/BIG_INT_BASE; a->a[i-1]-=k*BIG_INT_BASE;
    if(a->a[i-1]<0) k--, a->a[i-1]+=BIG_INT_BASE; a->a[i]+=k;
  }
  if(a->a[BIG_INT_SIZE-1]<0){
    rep(i,BIG_INT_SIZE) a->a[i]=-a->a[i]; bigIntOrderPt(a);
    rep(i,BIG_INT_SIZE) a->a[i]=-a->a[i];
  }
}

bigInt llToBigInt(ll a){
  bigInt c; int i;
  REP(i,1,BIG_INT_SIZE) c.a[i]=0; c.a[0]=a;
  bigIntOrderPt(&c); return c;
}

int bigIntGreaterThan(bigInt a,bigInt b){
  int i;
  for(i=BIG_INT_SIZE-1;i>=0;i--){
    if(a.a[i]>b.a[i]) return 1;
    if(a.a[i]<b.a[i]) return 0;
  }
  return 0;
}

int bigIntIsZero(bigInt a){
  int i; rep(i,BIG_INT_SIZE) if(a.a[i]) return 0; return 1;
}

bigInt bigIntPlus(bigInt a,bigInt b){
  int i; bigInt c;
  rep(i,BIG_INT_SIZE) c.a[i]=a.a[i]+b.a[i];
  bigIntOrderPt(&c); return c;
}

bigInt bigIntPlusEqual(bigInt *a,bigInt b){
  int i;
  rep(i,BIG_INT_SIZE) a->a[i]+=b.a[i];
  bigIntOrderPt(a);
}

bigInt bigIntMinus(bigInt a,bigInt b){
  int i; bigInt c;
  rep(i,BIG_INT_SIZE) c.a[i]=a.a[i]-b.a[i];
  bigIntOrderPt(&c); return c;
}

bigInt bigIntMultipleLL(bigInt a,ll b){
  int i; rep(i,BIG_INT_SIZE) a.a[i]*=b;
  bigIntOrderPt(&a); return a;
}

bigInt bigIntMultiple(bigInt a,bigInt b){
  int i,j,ii,jj; bigInt c;
  for(ii=BIG_INT_SIZE-1;ii;ii--) if(a.a[ii]) break; ii++;
  if(ii==1) return bigIntMultipleLL(b,a.a[0]);
  for(jj=BIG_INT_SIZE-1;jj;jj--) if(b.a[jj]) break; jj++;
  if(jj==1) return bigIntMultipleLL(a,b.a[0]);
  rep(i,BIG_INT_SIZE) c.a[i]=0;
  rep(i,ii)if(a.a[i])for(j=0;j<jj&&i+j+1<BIG_INT_SIZE;j++) c.a[i+j]+=a.a[i]*b.a[j];
  bigIntOrderPt(&c); return c;
}

void bigIntDivisionsLL(bigInt a,ll b,bigInt *c,ll *d){
  int i;
  rep(i,BIG_INT_SIZE) c->a[i]=a.a[i];
  for(i=BIG_INT_SIZE-1;i;i--)
    c->a[i-1]+=(c->a[i]%b)*BIG_INT_BASE, c->a[i]/=b;
  *d = c->a[0]%b; c->a[0]/=b;
}

/* c=a/b, d=a%b */
void bigIntDivisions(bigInt a,bigInt b,bigInt *c,bigInt *d){
  int i,j,s,sa,sb; ll ma,mb,mc; bigInt tmp;
  sa=bigIntSign(a); sb=bigIntSign(b);
  if(sa==-1) a=bigIntMultipleLL(a,-1);
  if(sb==-1) b=bigIntMultipleLL(b,-1);
  
  for(j=BIG_INT_SIZE-1;j;j--) if(b.a[j]) break;
  if(!j){
    REP(i,1,BIG_INT_SIZE) d->a[i]=0;
    bigIntDivisionsLL(a,b.a[0],c,&(d->a[0]));
  }else{
    for(i=BIG_INT_SIZE-1;i;i--) if(a.a[i]) break;
    s=i-j; if(s<0) s=0;
    rep(i,BIG_INT_SIZE) c->a[i]=0;
    while(s>=0){
      ma=0; mb=BIG_INT_BASE-1;
      while(ma!=mb){
        mc = (ma+mb)/2 + (ma+mb)%2;
        c->a[s]=mc; tmp=bigIntMultiple(*c,b);
        if(bigIntGreaterThan(tmp,a)) mb=mc-1; else ma=mc;
      }
      c->a[s]=ma; s--;
    }
    tmp = bigIntMultiple(b,*c);
    *d = bigIntMinus(a,tmp);
  }

  if(sa==-1 && sb==-1){
    *d=bigIntMultipleLL(*d,-1);
  } else if(sa==-1 && sb!=-1){
    *c=bigIntMultipleLL(*c,-1);
    *d=bigIntMultipleLL(*d,-1);
  } else if(sa!=-1 && sb==-1){
    *c=bigIntMultipleLL(*c,-1);
  }
}

bigInt bigIntDivision(bigInt a,bigInt b){
  bigInt c,d;
  bigIntDivisions(a,b,&c,&d);
  return c;
}

bigInt bigIntModular(bigInt a,bigInt b){
  bigInt c,d;
  bigIntDivisions(a,b,&c,&d);
  return d;
}

int bigIntSign(bigInt a){
  int i;
  for(i=BIG_INT_SIZE-1;i>=0;i--) if(a.a[i]){
    if(a.a[i]<0) return -1; else return 1;
  }
  return 0;
}

bigInt bigIntAbs(bigInt a){
  if(bigIntSign(a)==-1) return bigIntMultipleLL(a,-1LL); return a;
}

bigInt bigIntGCD(bigInt a,bigInt b){
  if(bigIntSign(a)==-1) a=bigIntMultipleLL(a,-1);
  if(bigIntSign(b)==-1) b=bigIntMultipleLL(b,-1);
  if(bigIntIsZero(a)) return b;
  return bigIntGCD(bigIntModular(b,a),a);
}

bigInt bigIntLCM(bigInt a,bigInt b){
  return bigIntMultiple(bigIntDivision(a,bigIntGCD(a,b)),b);
}

int bigIntToChar(bigInt a,char ret[]){
  int i,j,s=0,len=0; char ct[BIG_INT_CHAR_SIZE]; ll lt;
  if(bigIntSign(a)==-1){
    ret[0]='-'; len=bigIntToChar(bigIntMultipleLL(a,-1LL),ret+1); return len+1;
  }
  rep(i,BIG_INT_SIZE){
    lt=a.a[i]; rep(j,BIG_INT_DIGITS) ct[s++]=lt%10, lt/=10;
  }
  j=0;
  while(s--){
    if(ct[s]) j=1;
    if(j) ret[len++]=ct[s]+'0';
  }
  if(!len) ret[len++]='0';
  ret[len]='\0'; return len;
}

void printBigInt(bigInt a){
  int i,k; char tmp[BIG_INT_CHAR_SIZE];
  k=bigIntToChar(a,tmp); rep(i,k) putchar(tmp[i]);
}

void putBigInt(bigInt a){
  char tmp[BIG_INT_CHAR_SIZE];
  bigIntToChar(a,tmp); puts(tmp);
}

bigInt bigIntDivisionLL(bigInt a,ll b){
  bigInt res; ll tmp;
  bigIntDivisionsLL(a, b, &res, &tmp);
  return res;
}



bigRational bigRationalReduction(bigRational c){
  bigInt k=bigIntGCD(c.a,c.b); k=bigIntAbs(k);
  if(bigIntGreaterThan(k,bigIntOne()))
    c.a=bigIntDivision(c.a,k), c.b=bigIntDivision(c.b,k);
  return c;
}

int bigRationalIsZero(bigRational a){
  return bigIntIsZero(a.a);
}

bigRational bigRationalZero(){
  bigRational c; c.a=bigIntZero(); c.b=bigIntOne(); return c;
}

bigRational llToBigRational(ll a){
  bigRational c; c.a=llToBigInt(a); c.b=bigIntOne(); return c;
}

bigRational rationalToBigRational(rational a){
  bigRational c; c.a=llToBigInt(a.a); c.b=llToBigInt(a.b); return c;
}

bigRational twoLlToBigRational(ll a,ll b){
  bigRational c; c.a=llToBigInt(a); c.b=llToBigInt(b); return bigRationalReduction(c);
}

bigRational twoBigIntToBigRational(bigInt a,bigInt b){
  bigRational c; c.a=a; c.b=b; return bigRationalReduction(c);
}

void printBigRational(bigRational c){
  printBigInt(c.a);
  if(c.b.a[0]==1 && c.b.a[1]==0) return;
  putchar('/'); printBigInt(c.b);
}



#define N 151

bigInt fact[N], nowLCM, dpLCM[N], tmp, tmpA, tmpB;
bigRational dp[N];

int main(){
  int i,j,k,l,m,n,a;
  ll modA, modB;
  bigInt used, go, mul, sum, add;
  bigRational per, one = llToBigRational(1);
  int size;

  ps = getPrime(250,p);

  fact[0] = bigIntOne();
  REP(i,1,N) fact[i] = bigIntMultipleLL(fact[i-1],i);

  dp[0] = dp[1] = bigRationalZero();
  dpLCM[0] = dpLCM[1] = bigIntZero(); nowLCM = bigIntOne();
  
  REP(i,2,N){
    used = bigIntZero();
    sum = bigIntZero();
    rep(j,i){
      if(j==1) continue;
      if(j==0) go = bigIntOne();
      else {
        k = j - 2;
        a = i - k - 2;
        go = bigIntMultipleLL(fact[k], (a+1)*(k*k+k+1));
      }
      bigIntPlusEqual(&used, go);
      add = bigIntMultiple(dpLCM[j], go);
      bigIntPlusEqual(&sum, add);
    }
    dp[i].a = sum; dp[i].b = nowLCM;
    dp[i].a = bigIntPlus(dp[i].a, bigIntMultiple(dp[i].b, fact[i]));
    dp[i].b = bigIntMultiple(dp[i].b, used);

    dp[i+1] = dp[i];
    rep(k,ps){
      for(;;){
        bigIntDivisionsLL(dp[i].a, p[k], &tmpA, &modA);
        if(modA) break;
        bigIntDivisionsLL(dp[i].b, p[k], &tmpB, &modB);
        if(modB) break;
        dp[i].a = tmpA;
        dp[i].b = tmpB;
      }
    }

    tmp = nowLCM;
    nowLCM = bigIntLCM(nowLCM, dp[i].b);
    mul = bigIntDivision(nowLCM, tmp);
    rep(j,i) dpLCM[j] = bigIntMultiple(dpLCM[j], mul);
    mul = bigIntDivision(nowLCM, dp[i].b);
    dpLCM[i] = bigIntMultiple(dp[i].a, mul);
  }

  scanf("%d",&size);
  while(size--){
    scanf("%d",&n);
    printBigRational(dp[n]); puts("");
  }
  
  return 0;
}

コメント

コメントの投稿

コメントの投稿
管理者にだけ表示を許可する

トラックバック

トラックバック URL
http://rsujskf.blog32.fc2.com/tb.php/1895-f284b367
この記事にトラックバックする(FC2ブログユーザー)

Appendix

Recent Articles

ブログ内検索

Ads


(プライバシーポリシー)
上記広告は1ヶ月以上更新のないブログに表示されています。新しい記事を書くことで広告を消せます。