プログラマ専用SNS ミクプラ

以下のアルゴリズムに従って，2項乱数を生成し，niko.csvというファイルを作成するプログラムを作ったのですが
プログラムを起動してみるといつまでたってもプログラムが終わりません．

どこがおかしいのでしょうか…
n,pの値に指定はないですが，n=1000,p=0.95で計算しました．
よろしくお願いします．

2項乱数B(n,p)生成のための逆関数法によるアルゴリズム
STEP 1: Generate a random number U.
STEP 2: c= p/(1-p), i= 0, pr= (1-p)n, F= pr.
STEP 3: If U< F, set X= iand stop.
STEP 4: pr= [c(n-i)/(i+1)]pr, F= F+ pr, i= i+ 1.
STEP 5: Go to STEP 3.

コード:

 #include<stdio.h>
#include<stdlib.h>
#include<math.h>

#define lambda 1.0

main(){
  FILE *fp;
  double u,p,f,c,pr;
  int n,i,j,x,z;
  
  fp = fopen("niko.csv","wt");
  
  if(fp==NULL){
    fprintf(stderr,"cannot open file !\n");
    exit(1);
  }
  printf("試行回数を入力してください: ");
  scanf("%d", &n);
  printf("成功確率を入力してください: ");
  scanf("%lf", &p);
  
  for(j=0;j<n;j++){
    u=(double)rand()/(RAND_MAX+1.0);
    c = p/(1-p);
    i=0; 
    pr=(1-p);
    for(z=1;i<n;z++){
      pr = pr*(1-p);
    }
    f=pr;
    while(u>=f){
      pr=(c*(n-i)/(i+1))*pr; 
      f=f+pr; 
      i=i+1;
    }
    x=i;
    fprintf(fp,"%d\n",x);
  }
  fclose(fp);
}

見直してみたら28行目のfor文
for(z=1;i<n;z++)
は
for(z=1;z<n;z++)
でした．ですがこの部分を直しても依然プログラムは終わらないままです…

どうかよろしくお願いします．

p=0.95なので、(1-p)=0.05であり、0.05を1000乗すると小数点以下に0がおよそ1300個連続する数になるため、
正の数は10の-325乗程度までしか表現できないIEEE 754 倍精度浮動小数点数では0になってしまいます。
0に何をかけても0であり、fにどれだけ0を足しても初期値pr=0のままです。
よって、u>=fがずっと成り立ち、無限ループになる可能性があります。
試しにn=10, p=0.01で実行したところ、正常終了しました。
また、計算用変数をdoubleからlong doubleにしたところ、ideoneでは正常終了しました。

コード:

#include<stdio.h>
#include<stdlib.h>
#include<math.h>

#define lambda 1.0

main(){
  FILE *fp;
  long double u,p,f,c,pr;
  double pin;
  int n,i,j,x,z;
  
  fp = fopen("niko.csv","wt");
  
  if(fp==NULL){
    fprintf(stderr,"cannot open file !\n");
    exit(1);
  }
  printf("試行回数を入力してください: ");
  scanf("%d", &n);
  printf("成功確率を入力してください: ");
  scanf("%lf", &pin);
  p=(long double)pin;
  
  for(j=0;j<n;j++){
    u=(long double)rand()/(RAND_MAX+1.0);
    c = p/(1-p);
    i=0; 
    pr=(1-p);
    for(z=1;z<n;z++){
      pr = pr*(1-p);
    }
    f=pr;
    while(u>=f){
      pr=(c*(n-i)/(i+1))*pr; 
      f=f+pr; 
      i=i+1;
    }
    x=i;
    fprintf(fp,"%d\n",x);
  }
  fclose(fp);
}

オフトピック

n=1000, p=0.95ので、
なぜかdouble型のままでもideoneで正常終了しました。
しかし、試しにprの値を出力させてみると、TLEになりました。
long doubleだとprの値を出力させても正常終了したので、最適化の影響かもしれません。

みけCATさん

上手くいきました，ありがとうございます！
とてもわかりやすき説明して頂き感謝しています．
表示できる少数の桁数に限界があることをすっかり忘れていました．

また，何かの期会がありましたらよろしくお願いしますｌ．

プログラマ専用SNS ミクプラ

プログラムが終わらない

プログラムが終わらない

Re: プログラムが終わらない

Re: プログラムが終わらない

Re: プログラムが終わらない