精華區beta GoodNews 關於我們 聯絡資訊
※ 引述《Action (真的...)》之銘言: : 我把 Pi 當平均值 (0.5), 不過算起來好像怪怪的... 我也check了一下那各 x^2-test: 答案是 probably random 也不知道有沒有問題 :~ 哀 JaVa啊 我感覺不到你 只好先用 C 打嚕 /*by heat*/ #include <stdio.h> #include <math.h> #define N 10000 #define rmax 10000 /*外部變數*/ long module[55]; long fx[rmax],fy[rmax]; static long x,y; /*******************************************************/ /*This the value of parameter to user for linear method*/ /*******************************************************/ long m=10000; long m1=100; long b=1621; int c=1; int j = 54; long a=523; /*********************************/ /* Avoid overflow */ /*********************************/ long mul(long r,long s) { long r0,r1,s0,s1; r0=r/m1;r1=r%m1; s0=s/m1;s1=s%m1; return (((r0*s1+s0*r1)%m1)*m1+r1*s1)%m; } /********************************/ /* Linear Congruential method */ /********************************/ void lrand() { int i=0; for(;i<=54;i++){ a=(mul(a,b)+c)%m; module[i]=a; } } /******************************/ /*Additive Congruenteal method*/ /******************************/ long arand() { j = ((j+1)%55); module[j]=(module[j]+module[(j+31)%55])%m; return module[j]; } /************************/ /*檢查是否亂數為probably*/ /************************/ void forcheck() { int i=0; long tx,ty; double squax,squay; for(;i<N;i++){ tx=fx[i]*fx[i]; ty=fy[i]*fy[i]; } squax = ((m*tx/N)-N); squay = ((m*ty/N)-N); if(squax-m<=2*sqrt(m)) puts("the random numbers x is probably random\n"); else puts("not random\n"); if(squay-m<=2*sqrt(m)) puts("the random numbers y is probably random\n"); else puts("not random\n"); } /*********************************/ /* 找出落在函數裡面的點 */ /*********************************/ void func(int *sum1,int *sum2) { double xc,yc; x=arand();xc=(double)x/m;fx[x]++; y=arand();yc=(double)y/m;fy[y]++; if(xc<=1.0&&yc<=1.0&&pow(xc,2)+pow(yc,2)<=1.0) *sum1+=1; if(xc<=1.0&&yc<=1.0&&yc<=pow(xc,2)) *sum2+=1; } main() { double squax,squay; int i,sum1=0,sum2=0; lrand(); for(i=0;i<=N;i++){ fx[i]=0; fy[i]=0; } for(i=0;i<=N;i++) func(&sum1,&sum2); printf("\npi => %f\n",(double)sum1/N*4.0); printf("\nx^2 => %f\n",(double)sum2/N); forcheck(); } -- . / / . . \▕ ▕ / ● ●/ /● ● /╱* Heat ¯╲\∪/╱¯ . / / Heat . ●● / ● ● ● \╲ * (..)/ / . . ●● -- ※ 發信站: 批踢踢實業坊(ptt.m8.ntu.edu.tw) ◆ From: Heat.m3.ntu.edu