|
投稿者：たろさ
投稿日：2016年 7月20日(水)17時42分2秒
|
|
|
|
『無限解析のはじまり』高瀬正仁:著 ちくま学芸文庫
P328
3.34オイラーの公式の証明
を眺めています。数式をプログラムにしてみました。
OPTION ARITHMETIC RATIONAL
LET format$="#." & REPEAT$("#",999)
LET e=1
LET f=1
FOR i=1 TO 450
LET f=f*i
LET e=e+1/f ! 自然対数の底(ネピア数)1000桁
NEXT i
! PRINT USING format$:e
FOR r=3 TO 4
LET n=10^r
LET x=2
LET x1=5
LET m=(((1+x/n)^n)/(E^x))/(((1+x1/n)^n)/(E^x1))
PRINT USING format$:m
NEXT r
END
LET n=10^6 辺りから重いです。LET n=10^7 何日?
1+x/n)^n=E^x
虚数z
(1+z/n)^n=E^z
この数式からからe^(i*pi)=-1に至るまでの過程をプログラムしてみたいと模索しています。
それと、非自明な零点の関係も
OPTION ARITHMETIC COMPLEX !複素数モード
SET WINDOW -40,40,-4,4 !x,y
DRAW grid(10,1)
SET POINT STYLE 6
FOR x=-40 TO 0.99 STEP 0.01
LET m=fnZeta(COMPLEX(1/2,X))
LET R=RE(m)
LET I=IM(m)
! PRINT m
SET POINT COLOR 3
PLOT POINTS:x,r
PLOT
SET POINT COLOR 4
PLOT POINTS:x,i
PLOT
NEXT x
FOR x=1.01 TO 40 STEP 0.01
LET m=fnZeta(COMPLEX(1/2,X))
LET R=RE(m)
LET I=IM(m)
! PRINT m
SET POINT COLOR 2
PLOT POINTS:x,r
PLOT
SET POINT COLOR 7
PLOT POINTS:x,i
PLOT
NEXT x
END
EXTERNAL FUNCTION fnZeta(X)
OPTION ARITHMETIC COMPLEX !複素数モード
OPTION BASE 0 !ゼロベース
LET M=200
DIM C(M)
LET C(0)=1
for K=1 to M
FOR J=K TO 1 STEP -1
LET C(J)=(C(J)+C(J-1))/2
NEXT j
NEXT k
for J=1 to M
LET U=U+((-1)^INT(J-1))*C(J)/J^X
IF C(J)<1E-15 THEN EXIT FOR
NEXT j
LET fnZeta=U/(1-2^(1-X))
END FUNCTION
このグラフの揺らぎと、素数の個数をx単位の間隔で数えた場合の揺らぎは?
どのように変化
私の場合は、予想もできません。
http://blogs.yahoo.co.jp/donald_stinger
|
|
|
|
投稿者：しばっち
投稿日：2016年 7月23日(土)20時39分45秒
|
|
|
|
> No.4103[元記事へ]
たろささんへのお返事です。
> (1+x/n)^n=E^x
e=lim(1+1/n)^n n→∞
n=2^(2^10) とする
2^1024≒10^300 10進300桁程度
OPTION ARITHMETIC DECIMAL_HIGH
LET N=10
LET A=.5
FOR I=1 TO N !'(((((((((.5^2)^2)^2)^2)^2)^2)^2)^2)^2)^2
LET A=A*A
PRINT A
NEXT I
LET A=1+A
FOR I=1 TO 2^N !' ((1+1/2^1024)^2)^1024
LET A=A*A
PRINT A
NEXT I
END
> 虚数z
> (1+z/n)^n=E^z
>
> この数式からe^(i*pi)=-1に至るまでの過程をプログラムしてみたいと模索しています。
exp(z)=lim(1+z/n)^n n→∞
z=πi , n=2^(2^10) として (a+bi)^2=(a*a-b*b)+(2*a*b)i (i^2=-1) を使って計算します
OPTION ARITHMETIC DECIMAL_HIGH
LET N=10
LET A=1 !'A+Bi
LET B=.5
FOR I=1 TO N
LET B=B*B
NEXT I
!'LET B=1
!'FOR I=1 TO 2^N !'(.5)^(2^10)
!' LET B=B/2
!'NEXT I
LET B=B*PI !' B=π/2^1024
FOR I=1 TO 2^N !' ((1+πi/2^1024)^2)^1024
LET X=A*A-B*B
LET Y=2*A*B
LET A=X
LET B=Y
PRINT A
!'PRINT B
NEXT I
END
※
A=COMPLEX(SQR((1+1/SQR(2))/2),SQR((1-1/SQR(2))/2))
↓
A^2=COMPLEX(1/SQR(2),1/SQR(2))
↓
A^4=COMPEX(0,1)
↓
A^8=COMPLEX(-1,0)
|
|
|
|
投稿者：たろさ
投稿日：2016年 7月25日(月)16時47分9秒
|
|
|
|
> No.4105[元記事へ]
しばっちさんへのお返事です。
ありがとうございます。
> > (1+x/n)^n=E^x
>
> e=lim(1+1/n)^n n→∞
> n=2^(2^10) とする
> 2^1024≒10^300 10進300桁程度
私の個人的感想です。 計算の速さに驚きました。
> > 虚数z
> > (1+z/n)^n=E^z
> exp(z)=lim(1+z/n)^n n→∞
> z=πi , n=2^(2^10) として (a+bi)^2=(a*a-b*b)+(2*a*b)i (i^2=-1) を使って計算します
programの計算結果を見ての感想です。
.70710678118654752440084436210484903=SQR(0.5)
1.59906916149560598660344186636935668
-1.00000000000000000000000000000000000
SQR(0.5)は、何度も計算結果に出てくるので、見覚えがあります。
話しは、ぶっ飛びます。
What Is The Factorial Of 1/2? SURPRISING (1/2)! = (√π)/2
(1/2)! = (√π)/2
私が素直にプログラムすると
!n!(n階乗)
LET f=1
LET n=4
FOR i=0 TO n-1
LET f=f*(n-i)
NEXT i
PRINT f
LET f=1
LET n=0.5
FOR i=0 TO n-1
LET f=f*(n-i)
NEXT i
PRINT f
END
---------------------------
計算結果 24 1
(1/2)! = (√π)/2
になりません。
ここはメビュウスの輪
メビウスの帯
オイラー定数
#3801
初めてこのプログラムを見た時
私の個人的感想です。 計算の速さに驚きました。
数字の4とpi
zeta(±0.5*x*i)=0
zeta(sin(pi)^x*i)=-0.5
私の場合は まだ、勉強不足です。
http://blogs.yahoo.co.jp/donald_stinger
|
|
|
|
投稿者：nagram
投稿日：2016年 8月 5日(金)16時12分10秒
|
|
|
|
> No.4103[元記事へ]
たろささんへのお返事です。
> 『無限解析のはじまり』高瀬正仁:著 ちくま学芸文庫
> P328
> 3.34オイラーの公式の証明
>
> を眺めています。数式をプログラムにしてみました。
>
> OPTION ARITHMETIC RATIONAL
> LET format$="#." & REPEAT$("#",999)
> LET e=1
> LET f=1
> FOR i=1 TO 450
> LET f=f*i
> LET e=e+1/f ! 自然対数の底(ネピア数)1000桁
> NEXT i
> ! PRINT USING format$:e
>
> FOR r=3 TO 4
> LET n=10^r
> LET x=2
> LET x1=5
> LET m=(((1+x/n)^n)/(E^x))/(((1+x1/n)^n)/(E^x1))
> PRINT USING format$:m
> NEXT r
> END
>
> LET n=10^6 辺りから重いです。LET n=10^7 何日?
>
1000桁モードで実行すれば、r=18まで高速に計算できます。
OPTION ARITHMETIC DECIMAL_HIGH
DECLARE EXTERNAL FUNCTION EXP ! 外部関数宣言(必須)
LET x=2
LET x1=5
LET ex=EXP(x)
LET ex1=EXP(x1)
FOR r=3 TO 18
LET n=10^r
LET m=((1+x/n)^n/ex)/((1+x1/n)^n/ex1)
PRINT m
NEXT r
END
!
EXTERNAL FUNCTION EXP(x) ! EXP関数 (精度約1000桁 有理数モード不可)
OPTION ARITHMETIC DECIMAL_HIGH
DECLARE NUMERIC a,lim
FUNCTION ex(t,n) ! e^a テイラー展開
IF t>lim THEN LET ex=t+ex(t*a/n,n+1) ELSE LET ex=0
END FUNCTION
LET a=ABS(x) ! 定義域 a<=2321.0057737… (10進モード a<=227.9559…)
LET lim=EPS(2.71828^INT(a))*10^(-9) ! EPS(e^a)の9桁下まで計算し加算
IF x>=0 THEN LET EXP=ex(1,1) ELSE LET EXP=1/ex(1,1)
END FUNCTION
外部関数EXP(x)は、約1000桁の精度があります。テイラー展開は関数の再帰呼出しを利用し、最後の項から足し上げることによって精度を高くしました。xは負数でも可能です。
PRINT EXP(1) で自然対数の底を小数点以下1007桁まで得られます。
【おまけ】静岡理工科大学の菅沼義昇先生が公開しているガンマ関数の近似値を求めるプログラムをBASICに移植しました。ガンマ関数が階乗の拡張と言われる理由がわかります。
!/****************************/
!/* ガンマ関数の計算 */
!/* coded by Y.Suganuma */
!/****************************/
DECLARE FUNCTION gamma
DECLARE NUMERIC x,y
INPUT PROMPT "data?":x !-69.38<=x<56.905…(10進モード)
PRINT "x ="; x
LET y=gamma(x)
PRINT " Γ(x) ="; y
PRINT "x*Γ(x) ="; x*y
PRINT "Γ(x+1) ="; gamma(x+1)
!/****************************************/
!/* Γ（ｘ）の計算（ガンマ関数，近似式） */
!/* ier : =0 : normal */
!/* =-1 : x=-n (n=0,1,2,･･･) */
!/* return : 結果 */
!/****************************************/
FUNCTION gamma(x)
LOCAL ier,err,g,s,t,v,w,y,k
LET ier=0
IF x>5 THEN
LET v=1/x
LET s=((((((-0.000592166437354*v+0.0000697281375837)*v+0.00078403922172)*v &
& -0.000229472093621)*v-0.00268132716049)*v+0.00347222222222)*v+0.0833333333333)*v+1
LET g=2.506628274631001*EXP(-x)*x^(x-0.5)*s
ELSE
LET err=1E-20
LET w=x
LET t=1
IF x<1.5 THEN
IF x<err THEN
LET k=INT(x)
LET y=k-x
IF ABS(y)<err OR ABS(1-y)<err THEN LET ier=-1
END IF
IF ier=0 THEN
DO WHILE w<1.5
LET t=t/w
LET w=w+1
LOOP
END IF
ELSE
IF w>2.5 THEN
DO WHILE w>2.5
LET w=w-1
LET t=t*w
LOOP
END IF
END IF
LET w=w-2
LET g=(((((((0.0021385778*w-0.0034961289)*w+0.0122995771)*w-0.00012513767)*w &
& +0.0740648982)*w+0.0815652323)*w+0.411849671)*w+0.422784604)*w+0.999999926
LET g=g*t
END IF
IF ier=-1 THEN PRINT "[ ier=";ier;"]"
LET gamma=g
END FUNCTION
END
|
|
|
戻る