※ 引述《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
. ●● / ● ● ● \╲ *
(..)●/ /● .
. O ●●
--
※ 發信站: 批踢踢實業坊(ptt.m8.ntu.edu.tw)
◆ From: Heat.m3.ntu.edu