十進BASIC 第2掲示板 過去ログ

!準モンテカルロ法　∫[0,1]f(x)dx=1/N*Σ[i=1,N]f(xi)

DEF f(x)=x^2 !被積分関数 1/3

LET ITER=1000 !発生個数

LET UINT_MAX=2^32-1 !4294967295UL
LET s=0
FOR n=0 TO ITER-1
   LET g2n=0 !超一様分布列　ヴァンデルコープット列（van der Corput）
   FOR i=0 TO 30
      LET g2n=(g2n+REMAINDER(r_sft(n,i),2))*2
   NEXT i
   IF REMAINDER(r_sft(n,31),2)=1 THEN LET g2n=g2n+1

   LET x=g2n/(UINT_MAX+1) !区間[0,1)={0,1/2,1/4,3/4,1/8,5/8,3/8,7/8,1/16,9/16,5/16,13/16,…}
   !!!PRINT x !debug

   LET s=s+f(x) !Σf
   PRINT s/(n+1)
NEXT n
!PRINT s/ITER !結果を表示する

FUNCTION r_sft(n,i) !Ｃ言語 n>>i 等価
   LET r_sft=IP(n/2^i)
END FUNCTION

END

!準モンテカルロ法　∫[0,1]f(x)dx=1/N*Σ[i=1,N]f(xi)

DEF f(x)=x^2 !被積分関数 1/3

LET ITER=1000 !発生個数

LET UINT_MAX=2^32-1 !4294967295UL　※32bit
LET s=0
FOR n=0 TO ITER-1
   LET g2n=0 !超一様分布列　ヴァンデルコープット列（van der Corput）
   LET a=n
   FOR i=1 TO 32 !２進法への変換より、ビット列の並びを反転（逆順に）する
      LET aa=INT(a/2)
      LET g2n=g2n*2+(a-aa*2) !下位桁を上位桁へ
      LET a=aa
   NEXT i
   !!!PRINT n;g2n,right$(REPEAT$("0",32)&BSTR$(g2n,2),32) !debug

   LET x=g2n/(UINT_MAX+1) !区間[0,1)={0,1/2,1/4,3/4,1/8,5/8,3/8,7/8,1/16,9/16,5/16,13/16,…}
   !!!PRINT x !debug

   LET s=s+f(x) !Σf
   !PRINT s/(n+1)
NEXT n
PRINT s/ITER !結果を表示する

END