くだすれFORTRAN（超初心者用）

くだすれFORTRAN（超� ..

596:デフォルトの名無しさん
06/07/05 05:08:37
さぁテポドンの弾道をFORTRANで計算しろ！

597:デフォルトの名無しさん
06/07/05 07:31:54
>>596
do i=1.10
宇宙の果てまで、ちゅぼーん

598:デフォルトの名無しさん
06/07/05 14:07:32
>>597
ほんとに１０発撃ったらしいぞｗ　４発は不発だったらしいがｗｗｗ

599:デフォルトの名無しさん
06/07/05 19:01:15
問題：L行M列とM行N列の行列の積を計算するサブルーチン副プログラムを作れ。
サブルーチン副プログラムを呼び出す主プログラムも作成し、結果を出力すること。
主プログラムでは二つの行列のデータ読み込みと結果の出力を行うこととする。
(この後プログラムが正しく動作することを確認するための例題が続きます
途中の文章をいれたプリント文はそれを入力させるためのものです)
real function SEKI(A,B,C,J,K,L)
real A(2,2),B(2,3),C(2,3)
c
do 30 I=1,L
do 20 J=1,N
C(I,J)=0
do 10 K=1,M
C(I,J)=C(I,J)+A(I,K)*B(K,J)
10 continue
20 continue
30 continue
c
return
end

600:デフォルトの名無しさん
06/07/05 19:02:58
c
do 70 K=1,N
print*,C(1,K)
70 continue
do 80 K=1,N
print*,C(2,K)
80 continue
end
c
print*,'1,-2,3,0,4,-1,2,1,0,-3を入力して下さい'
do 50 I=1,L
do 40 K=1,M
read*,A(I,K)
40 continue
50 continue

601:デフォルトの名無しさん
06/07/05 19:03:53
do 70 K=1,M
do 60 J=1,N
read*,B(K,J)
60 continue
70 continue
c
end

エラーが出るんですがもうどこが間違ってるのか全然わかりません。お手上げです。
誰か分かる方、ご教授願います。

602:デフォルトの名無しさん
06/07/05 19:26:34
>>599
function ではなく　SUBROUTINE だべｗ
配列引数のサイズを、固定長で与えては一般性のあるプログラムにならない。
（これを次元の宣言と割り切る、はみ出し無視のFORTRAN６６時代の猛者なら別だがｗｗｗ）

あとは、まぁまぁ大丈夫でないの？　
７０番が2回出てくる上に、妙なところにend文が割り込んで理解に苦しむが、
こぴぺ時の間違いとみなしておくｗ

603:デフォルトの名無しさん
06/07/05 20:05:21
できたおｗ

program hoge
implicit none
real, allocatable :: A(:,:), B(:,:), C(:,:)
integer :: i,j, l,m,n
read *, l, m, n
allocate(A(l,m), B(m,n), C(l,n))
read *, ((A(i,j), j=1,m), i=1,l)
read *, ((B(i,j), j=1,n), i=1,m)
call seki(A, B, C, l, m, n)
do i=1,l
print *, (C(i,j), j=1,n)
end do
deallocate(A, B, C)
stop
contains
subroutine seki(A,B,C, l,m,n)
integer, intent(in) :: l,m,n
real, intent(in) :: A(1:l,1:m), B(1:m,1:n)
real, intent(out) :: C(1:l,1:n)
C = MATMUL(A, B)
end subroutine seki
end program hoge

604:デフォルトの名無しさん
06/07/05 20:19:03
>>603
おまｗ　
９０でやったらわざわざMATMUL相当のサブルーチン作る意味無いｗ

ついでに７７でも作ってやれｗ

605:デフォルトの名無しさん
06/07/05 20:53:04
77だとこんな感じだお

SUBROUTINE SEKI(A,B,C,L,M,N)
REAL A(L,M), B(M,N), C(L,N)
DO 30 I=1,L
DO 20 J=1,N
C(I,J) = 0.
DO 10 K=1,M
C(I,J) = C(I,J) + A(I,K)*B(K,J)
10 CONTINUE
20 CONTINUE
30 CONTINUE
RETURN
END
C
PROGRAM HOGE
REAL A(3,3), B(3,2), C(3,2)
L=3
M=3
N=2
READ(*,*) ((A(I,J), J=1,M), I=1,L)
READ(*,*) ((B(I,J), J=1,N), I=1,M)
CALL SEKI(A,B,C, L,M,N)
DO 10 I=1,L
WRITE(*,*) (C(I,J), J=1,N)
10 CONTINUE
STOP
END

606:学生（アホ１６歳）
06/07/05 23:57:31
以下の課題についてですが…
【自然対数 e は，1 + (1)-1 + (2*1)-1 + (3*2*1)-1 + ・・・ + (10*9*8*・・・・*2*1)-1 というように近似することができます．この値をディスプレイに表示するプログラムを作成しなさい．】

ちなみに７７で書かないと点くれないそうで…（苦笑

７７で階乗計算って…どうすればいいのですか？

607:デフォルトの名無しさん
06/07/06 00:30:23
>>606
1! = 0!*1 = 1*1 = 1
2! = 1!*2 = 1*2 = 2
3! = 2!*3 = 2*3 = 6
4! = 3!*4 = 6*4 = 24
　　：
n! = (n-1)! * n

k = 1
do 10 n=1, 10
k = n * k
write(*,*) 'n = ',n, ', n! =', k
10 continue
end

608:デフォルトの名無しさん
06/07/06 01:18:11
>>606
階乗は>>607でいいが、１項毎に定義どうりに計算するのではなく
一つ前の項に違いだけ掛ければいい。
y=1, y = y*x/1, y = y*x/2, y = y * x/3 ..............

が、さらに順を変えて尻から進めるほうがもっといい。
1 + x (1 + x(1 + x/3) / 2 ) / 1

609:学生（アホ１６歳）
06/07/06 03:22:58
>>607-608
ありがとうございます！
出来ました。

610:デフォルトの名無しさん
06/07/07 11:11:35
fortranでC言語の#defineみたいに、プリプロセッサで定数指定できるのってありますか？

611:デフォルトの名無しさん
06/07/07 13:02:41
拡張子を .F にしたらプリプロセッサを解析してくれるコンパイラも多い。

612:610
06/07/07 13:35:41
>>611
返信どうもっす。じゃあ、fortranの標準では#defineみたいなのないんですね？

613:デフォルトの名無しさん
06/07/07 13:56:35
>>611
標準では存在しない。

#define は、便利な反面、批判も多いのでたぶん意図的に新しい規格でも不採用とおもわれ。
ただ問題点を無くした条件コンパイルなどのプリプロセッサへの対応は論じつづけられている。

614:610
06/07/07 18:17:48
>>613

なるほど、納得いたしました。どうもありがとうございます！

615:デフォルトの名無しさん
06/07/07 20:37:06
定数は parameter だな。

616:デフォルトの名無しさん
06/07/11 19:31:21
学校の課題をやっていただきたいのですが、、、
次の数列の和を求めよ。（DO文）
1/1+x=1-x+x**2-x**3+x**4-x**5・・・
打ち切り条件　1.0*10**-6
提出明日です。よろしくおねがいします。

617:デフォルトの名無しさん
06/07/11 22:20:20
割り込みすいません…
分数の足し算プログラムですが、
どこが違うかわからないので教えてくださいませんか？
PROGRAM FRACTION
INTEGER I,J,K,L,M,N
READ(*,*) I,J,K,L
CALL TUUBUN(I,J,K,L)
WRITE(*,*)M,N
CALL REDUCTION(M,N)
M=M/A
N=N/A
WRITE(*,*)M,'/',N
END

SUBROUTINE TUUBUN(A,B,C,D)
INTEGER A,B,C,D
M=A*D+B*C
N=C*D
RETURN
END

618:デフォルトの名無しさん
06/07/11 22:24:36
>>616
FORTRAN７７で適当に作った。適当に直してくれ。
数学的には　-1<x<1 で収束する式だが、１に近づくにつれて収束が遅くなる。
つまり、たくさんの項の和を取らねばならない。しかし、延々とやると数値誤差が
たまるので意味がなくなる。ゆえに適当なところであきらめて打ち切る。

PROGRAM VIP
PARAMETER (eps = 1.0e-6)
WRITE(6, *) 'INPUT x (-1 < x < 1)'
READ(5, *) x
s = 1.0
t = 1.0
DO 10 i = 1, 10000
t = - t * x
s = s + t
IF (ABS(t) .LT. eps) GOTO 999
10 CONTINUE
STOP 'NOT CONVERGED'
999 WRITE(6, *) ' sum =', s, ' iteration', i
WRITE(6, *) ' 1/(1+x) =', 1.0 / (1.0 + x)
STOP
END

619:デフォルトの名無しさん
06/07/11 22:32:24
>>617
サブルーチンの変数の範囲について根本的に誤解している。

メインプログラムとサブルーチンでは変数はお互いに見えない。

ゆえに、通分した分子と分母はM、Nはサブルーチンの引数で返す必要がある。

CALL TUUBUN(I,J,K,L, m, n)

SUBROUTINE TUUBUN(iA,iB,iC,iD, m, n)
M=iA*iD+iB*iC
N=iC*iD
RETURN
END

あと、暗黙の型は死んでも守るべし。周りがなんと言おうと守るべし。
かっぱに尻子玉を抜かれそうになっても寧ろ暗黙の型守るべし。

620:617
06/07/12 10:25:30
>>619
途中からなくなっていたので、
もう一度直して
PROGRAM FRACTION
INTEGER I,J,K,L,M,N
READ(*,*) I,J,K,L
CALL TUUBUN(I,J,K,L)
CALL REDUCTION(M,N)
M=M/A N=N/A
WRITE(*,*)M,'/',N
END

SUBROUTINE TUUBUN(iA,iB,iC,iD)
INTEGER iA,iB,iC,iD
M=iA*iD+iB*iC
N=iC*iD
RETURN
END

621:617
06/07/12 10:26:19
>>620
SUBROUTINE REDUCTION(P,Q)
INTEGER P,Q,X,Y
X=P
Y=Q
10 IF(X.EQ.Y) THEN GO TO 20
ELSE IF(X.GT.Y)
THEN X=X-Y
ELSE Y=Y-X
END IF
END IF GO TO 10
20 CONTINUE
A=X
RETURN
END

622:デフォルトの名無しさん
06/07/12 14:25:40
>>621
PROGRAM FRACTION
READ(*,*) I, J, K, L
CALL TUUBUN(I, J, K, L, M, N)
CALL REDUCTION(M, N, IA)
M = M / IA
N = N / IA
WRITE(*,*) M, '/', N
END
SUBROUTINE TUUBUN(N1, N2, ID1, ID2, NUMERA, IDENOM)
NUMERA = N1 * ID2 + N2 * ID1
IDENOM = ID1 * ID2
RETURN
END
SUBROUTINE REDUCTION(NUMERA, IDENOM, IGCD)
IX = NUMERA
IY = IDENOM
10 IF(IX .EQ. IY) THEN
GO TO 20
ELSE IF (IX .GT. IY) THEN
IX = IX - IY
ELSE
IY = IY - IX
END IF
GO TO 10
20 CONTINUE
IGCD = IX
RETURN
END
>>622

623:さおり
06/07/12 16:17:47
すいません！無理を承知でカキコします。
私は今理系で某国立大学に通ってるんですけど、水商売のバイトにハマってしまって全然授業についていってない状態なんです。
それでプログラミングの授業があるんですけど課題が出て「放物の運動（空気抵抗を考慮する）」のプログラムをfortranで作れって課題が出たんです。。
初速度と投げる角度を入力して一定時間ごと（１秒ごととか）の物体の位置をｘとｙで出すらしいんです。
私、授業にほとんど出てなくて本当にまったくわかってない状態なんです。
少し真面目な大学なんで私みたいな子は浮いてしまって、学校には助けてくれる友達もいないんです。
誰かプログラム作ってください。お願いします。この授業を落としてしまうと本当に留年なんです。おねがいします！！

624:デフォルトの名無しさん
06/07/12 16:43:42
>>623
まずはおっぱいうｐ汁

625:デフォルトの名無しさん
06/07/12 21:23:23
>>623
担当教官をお店で接待すれば優もらえるおｗ

626:デフォルトの名無しさん
06/07/13 00:14:44
金があるなら留年しろ。

627:デフォルトの名無しさん
06/07/13 01:04:04
俺は少年なら助けてやるのだが・・・・女は興味ないｗｗｗ

628:デフォルトの名無しさん
06/07/13 07:18:28
>>623
課題だけ書く方が答えが得られるすれのようです。

629:さおり
06/07/13 10:47:04
本当にお願いします！
いろいろ事情があって留年できないんです。

630:デフォルトの名無しさん
06/07/13 12:22:06
>>629
課題は省略しないで出題されたまま提示するべき
ヒントになるかどうか知らん
URLﾘﾝｸ(www12.plala.or.jp)

631:デフォルトの名無しさん
06/07/13 12:26:25
>>623
空気定数って物体のサイズによって二通りの計算方法があるけど、どっちを使うの？

632:631
06/07/13 12:29:35
>>629
>>630 のヒントの方の計算方法だと楽でいいよね。さあどっちだろうか？

633:631
06/07/13 16:25:24
>>631
×空気定数→○空気抵抗

634:さおり
06/07/13 21:43:35
えーっと空気抵抗が速度の２乗に比例する場合らしいです。本当にありがとうございます。

635:デフォルトの名無しさん
06/07/13 22:08:31
解決したみたいだから、次の宿題。

636:デフォルトの名無しさん
06/07/13 22:56:22
これから宿題をストーリーつきで頼む人は、

『ボク宿題が解けなくて困ってるんだ。
お兄ちゃん様に一緒に解いてほしいなぁ～。云々』

という、もっとやる気を出させるストーリーを考えること！

637:デフォルトの名無しさん
06/07/13 23:00:37
>>634
ああ、めんどくさい方ね。
それ、解析的に解が求まらないから、微分方程式出して数値解析するしかないかな？
一応、微分方程式を立ててみるけど、チェックしてみて
vx, vy を速度の x,y 成分として、物体の質量を m 、重力定数を g、空気抵抗の
係数を k とすると

m dvx/dt = -k vx √(vx^2+vy^2)
m dvy/dt = -k vy √(vx^2+vy^2) - mg

となるかな？あとは、授業で習った数値解析で vx, vy をまず求めるのだが……

638:デフォルトの名無しさん
06/07/13 23:19:43
>>636
兄貴ぃ～！
俺、宿題分かんないッス！！
徹夜で手取り足取り教えて欲しいッス！！

639:デフォルトの名無しさん
06/07/14 00:15:18
>>638
サブ！　手取り腰取り教えてやるぜ！！

640:デフォルトの名無しさん
06/07/14 01:05:54
>>637
最後まで解けよｗ

641:デフォルトの名無しさん
06/07/14 02:10:46
>>637
言いだしっぺのh(ry

642:641
06/07/14 02:11:58
間違った
637でなくて>>640

643:637
06/07/14 10:13:37
>>640
そうじゃなくて、学校で習った数値解析法で解かないと
合格になる気がしないんですけど。XX法とか言ってもらえれば
続きも出来るのですが……

644:デフォルトの名無しさん
06/07/14 11:08:41
>>634の書き込みを見ると、「らしい」なんて書いているから、課題の正確な内容を
把握してなさそう。仮に折角>>637が解いたところで「やっぱり違った」とかで別解を
要求してきそうだ。

645:デフォルトの名無しさん
06/07/14 13:55:36
C++スレ見たいに宿題依頼のテンプレが必要なんでね？

見たところ、宿題はいまだF77中心のようだし。ｆ９０かｆ７７かも選択してもらわないと。

646:さおり
06/07/14 16:17:43
そんなこと絶対しません。よろしくおねがいします。

647:デフォルトの名無しさん
06/07/14 18:27:20
つーことで、さおりタンの仕事は
1. 微分方程式の確認
2. 授業でならった微分方程式の数値解析法の(せめて)名前の確認
だね。

648:デフォルトの名無しさん
06/07/14 19:44:38
>>646
「そんなことしません」じゃなくて、レポートとして期待されているプログラムを
作成するためには情報が足りねーの。
課題を入手してきなさい。あと、おっぱいうｐな。

649:デフォルトの名無しさん
06/07/14 20:48:58
>>637に従って解いたｗ
微分方程式は、一番原始的な方法で積分したｗ　
大学１，２生ならこれでも合格ラインを超えるだろう。３，４年生ならルンゲ・クッタを使うべきだろう。
原点（０，０）から初速ｖｘ、ｖｙで玉を放り投げる。地面より下に座標がきたら（ｙ＜０．０）終了する。

PROGRAM yamajun
PARAMETER( g = 9.8, fric = 0.3 )
t = 0.0
dt = 1.0 / 8.0
amass = 1.0
x = 0.0
y = 0.0
vx = 10.0
vy = 10.0
ax = 0.0
ay = 0.0
1 CONTINUE
WRITE(9, *) x, y, t !, vx, vy !, ax, ay
IF ( y < 0.0 ) GOTO 999
t = t + dt
x = x + dt * vx
y = y + dt * vy
ax = - fric * vx * SQRT(vx**2 + vy**2) / amass
ay = -g - fric * vy * SQRT(vx**2 + vy**2) / amass
vx = vx + dt * ax
vy = vy + dt * ay
GOTO 1
999 STOP
END

FORTRANは、うほっの世界であるべきなので、さおりは次回からは
カヲルという少年に性転換すること。

650:637
06/07/14 21:17:58
>>649
thx.
いやー。数値解析の教科書とか見て気が重かったのだが、ありがとう。
さおりタン満足してくれるといいね。

651:649
06/07/14 21:40:27
>>650
>>649はあまりに手抜きｗ　誰かもっと良いの作ってやれｗ

２元２階の微分方程式なので、本のルンゲ・クッタルーチンにそのまま放り込めないｗ

ちょっと脳みそ使う必要がある。数独１回分くらいｗ　

今日暑いから考えらんねー　冷房ほしいよー

652:デフォルトの名無しさん
06/07/14 21:43:43
2,3,5.6,2.7,2,-,9,3
という「-」のような記号の混じった数列データが入ったファイルがあって、
それを下のように読み込みたいのですが、
real*8 a,b,c,d,e,f,g,h
read(10,*) a,b,c,d,e,f,g,h
これで「-」を読み込むときエラーが出て終了します。
「-」が出てきた場合には-1.0d-10などの任意の数字を代入しておく
ような方法はないでしょうか？
よろしくお願いします。

653:デフォルトの名無しさん
06/07/14 23:37:38
>>652　ほい。　F９０で書いたが、主要なところはF７７でもおｋ。　質問あればおｋ。
やり方は色々あるが、基本的にはまずデータを１行文字列として読み込んで、
カンマで区切って１こづつ数値として読み取るのが標準だろう。

ここでは読み取りエラーが起きてから例外処理的に”－”を処理しているが、
”－”が多いならはじめから一度文字列で読み込んで処理するほうがすっきりする。
PROGRAM vip
IMPLICIT NONE
CHARACTER (LEN = 136) :: text
REAL(KIND = 8) :: x(8)
INTEGER :: i, ipos, io
DO
READ(9, *, IOSTAT = io) x
IF (io == -1) EXIT
IF (io > 0) THEN
BACKSPACE(9)
READ(9, '(A)') text
DO i = 1, 8
ipos = INDEX(text, ',') - 1
IF (ipos == -1) ipos = LEN(text) ! 8th number
IF (TEXT(ipos:ipos) == '-') THEN
x(i) = -1.0d-10
ELSE
READ(TEXT(1:ipos), *) x(i)
END IF
text = text(ipos + 2:)
END DO
END IF
PRINT '(8g9.2)', x
END DO
STOP
END PROGRAM vip

654:６５３
06/07/14 23:39:00
■入力データ

2,3,5.6,2.7,2,1,9,3
2,3,5.6,2.7,2,-,9,3
2,3,5.6,2.7,2,2,9,3

■出力
2.0 3.0 5.6 2.7 2.0 1.0 9.0 3.0
2.0 3.0 5.6 2.7 2.0 -0.10E-09 9.0 3.0
2.0 3.0 5.6 2.7 2.0 2.0 9.0 3.0
Press any key to continue

655:デフォルトの名無しさん
06/07/15 03:04:03
　　　　　 |　|/ノ二__‐─ｧ　　ヽﾆニ二二二ヾ } ,'⌒ヽ
　　　　 /⌒!|　 =彳ｏ｡ﾄ￣ヽ　　　　　'´ !o_ｼ｀ヾ　| i/ ヽ !　　　　ウホッ！　いい男・・・
　　　　 ! ハ!|　　ー─　'　　i　　!　　 `' 　 '' " 　 ||ヽ　l |
　　　　| | /ヽ!　　　　　　　　|　　　　　　　　　　　 |ヽ i　!
　

656:652
06/07/15 09:43:09
>>653
ウホッ！

657:デフォルトの名無しさん
06/07/15 09:44:53
>>623,>>629,>>634,>>646

>>637をもとに>>649が解いたものは、m、k、初期位置、時間の刻みに
適当な値を代入して解いている。本来はあんたが提示すべき情報。

>>630,>>647,>>648を良く読め。

658:デフォルトの名無しさん
06/07/15 10:57:31
教えて！goo にも同じような質問が・・・

659:637
06/07/15 12:20:32
>>658
じゃあ我々の役目も終ったのかな？
みんなお疲れ様。

660:デフォルトの名無しさん
06/07/15 17:35:26
女が来るとふいんきが乱れるから、ホモOnlyスレにしないか？

661:デフォルトの名無しさん
06/07/15 18:59:43
プログラムの実行中に一瞬ポーズする関数ってありますか？

662:デフォルトの名無しさん
06/07/15 19:17:53
>>661
PAUSE

663:デフォルトの名無しさん
06/07/15 19:30:41
>>662
ただPAUSEはFORTRAN2003あたりで廃止勧告が出されている。
TSS利用で便利だから漏れは好きなんだけどなｗ

664:デフォルトの名無しさん
06/07/15 20:28:29
write(*,*) abs(-16.8-(6.8))>10.0
を出力すると「T」になってしまいます。
助けてください。

665:664
06/07/15 20:32:23
間違えました、
real*8 a,b
a=-16.80000000000000
b=- 6.80000000000000
write(*,*) abs(a-b)>10
とやったら「T」と出てしまいます。
何が原因でしょうか？

666:664
06/07/15 20:53:41
更に状況を書かせて頂きますと、
読み取りファイルに
-16.8,-6.8
とあって、
real*8 a,b
open(10,file="test")
read(10,*)a,b
write(*,*) abs(a-b)>10
とやったら「T」と出力されるのです。

667:デフォルトの名無しさん
06/07/15 20:59:25
FORTRAN77は大文字でプログラム書くのがデフォルトなのか

668:デフォルトの名無しさん
06/07/15 21:04:15
>>664
それは良くあることｗ　気にするな。　浮動小数点之罪。

669:デフォルトの名無しさん
06/07/15 21:06:04
>>667
本来のFORTRAN77の規格では大文字のみ。

そもそも、昔は小文字が表示できない端末とか良くあったｗ

670:デフォルトの名無しさん
06/07/15 21:22:32
>>664
参考までに、INTEL　FORTRANの結果を見てくれ。
PROGRAM unko
IMPLICIT NONE
REAL (KIND = 8) :: a, b
a = -16.8
b = -6.8
PRINT *, a, b, ABS(a - b) < 10
a = -16.8d0
b = -6.8d0
PRINT *, a, b, ABS(a - b) < 10.0d0
STOP
END PROGRAM unko

■結果
-16.7999992370605 -6.80000019073486 T
-16.8000000000000 -6.80000000000000 F
Press any key to continue

この場合は倍精度（８バイト）の変数に単精度（４バイト）の常数を与えているので、
長さの足りない４バイトにゴミが入っている。

■
とにかく浮動小数点の誤差とか桁落ちとか実数の比較とかでググレば説明が出るはず。
この辺は昔は電子計算機利用の一番最初に注意されて来た事柄だ。

671:664
06/07/15 21:45:51
>>670
どうもありがとうございます。
一応計算中の変数を画面に表示させてチェックしたのですが、
aは-16.8000000000000
bは -6.8000000000000
と表示されていました。
この上でabs(a-b)>10をやるとTになりました。

672:デフォルトの名無しさん
06/07/15 22:03:21
>>671
(T_T)
30桁くらい出力してごらん。
write ( *, '(2f35.30)' ) a, b

673:670
06/07/16 00:28:13
>>671
漏れの書き方が誤解を招いたらしい。すまんこ。

数直線のうちFORTRANであらわされるのは、０を中心として正負にほぼ同じだけの
範囲の中の、有限個の有理数だけである。数直線の上に、定規の目のように切った点が
あって、定規の目に挟まれた部分は近くの目で代表することになる。

（浮動小数点のばあい、この定規の目は体育で使う巻尺のように、０の近傍は細かく
区切ってあるが、遠くへ行くほど荒く区切るようになっている。）

コンパイラにも依存するのだが、基本的には単精度や倍精度の浮動小数点は、
２の（負を含む）べき乗の和であらわされるような有理数になっている。
このため、１０進表記で書いた有限小数が、かならずしも数直線の上に切った目に
一致するとは限らない。

つまり
１．５＝１＋０．５＝２＾０＋２＾－１
０．７５＝０．５＋０．２５＝２＾－１＋２＾－２
これらの数は、きっかり割り切れる（すなわち、定規の目の上にのっている）

ところが
０．８＝０．５＋０．２５＋０．０６２５＋・・・・＝２＾－１＋２＾－２＋２＾－４＋・・・・・
となって、数直線の上に切った目の上に載っていない。
つまり、近くの別の数に四捨五入されている。
>>672のように、FORMATで小数点いかの数を長く出すと、以下のようになる。
-16.7999992370605 -6.80000019073486 T
-16.8000000000000 -6.80000000000000 F
-16.800000000000000710542735760100 -6.799999999999999822364316059975
つまり、倍精度の－１６．８は、実際は-16.800000000000000710542735760100・・・・・・
という数に四捨五入されている。－６．８も同様。
（ただし、有効桁からすれば、＊フォーマットで出力される分しか制度が無い。
つまり＊フォーマットで出力される数の最後の桁程度の範囲で四捨五入がされている。）

674:デフォルトの名無しさん
06/07/16 14:38:42
　積分区間の下限Aと上限B、分割数N（偶数）を読み込んで、
Simpson公式を用いて、定積分I=(A,B)∫f(x)dxを計算するプログラムを作れ。
　ただし、数値積分の部分のプログラムは、任意の関数に対して使用できる
サブルーチン・サブプログラムを用いよ。その場合関数f(x)の計算は関数サブプログラムで行う。

　さらに上記の分割数Nが奇数に対しても使えるように拡張せよ。
Simpson公式を用いて、Nが奇数のときは1区画だけ台形公式を用いるようにすればよい。

手も足も出ません…、お願いできますか？

675:デフォルトの名無しさん
06/07/16 16:37:20
>>674
PROGRAM HOGE
READ(*,*) A, B, N
IF (MOD(N,2) .NE. 0) GOTO 9999　… ※１
C
DX = (B - A) / N
S = 0.
DO 10 I = 1, N/2
X1 = A + 2*(I-1)*DX
X2 = A + (2*I - 1)*DX
X3 = A + 2*I*DX
S = S + (F(X1) + 4*F(X2) + F(X3))
10 CONTINUE
S = S * DX / 3.
C IF (MOD(N,2) .NE. 0) THEN　… ※２
C X1 = B - DX
C X2 = B
C S = S + (F(X1) + F(X2)) * DX / 2
C END IF
WRITE(*,*) S
9999 CONTINUE
STOP
END

前半部分はそのまま。（※２の所はいらない）
後半部分は※１のIF文を消して※２のコメントアウトされてるIFブロックを追加する。

676:デフォルトの名無しさん
06/07/16 21:55:10
ファイルを読み込んだり作ったりするとき、
実行ファイルと同じディレクトリの中に
別のディレクトリを作って保存したりするにはどうしたら良いでしょうか。

677:デフォルトの名無しさん
06/07/16 22:07:35
>>676
既出の質問だが・・・・お答えしよう。

FORTRANの標準命令ではできない。

コンパイラメーカー固有の拡張命令でシステムの関数を呼び出すものがあるだろうから
それらを通じてやるしかない。マニュアルの後ろのほうの章を探して味噌。

678:676
06/07/16 23:32:02
>>677
ｻﾝｸｽｺ

679:デフォルトの名無しさん
06/07/17 12:21:13
番組の途中ですが・・・・・・・・・・・・

盧武鉉（ノ・ムヒョン）大統領が今月11日に行われた与党ヨルリン・ウリ党の指導部
および国会の統一外交通商委員会に所属する議員らとの晩さん懇談会で行った
発言が波紋を呼んでいる。

一部新聞は懇談会出席者の証言を引用し、盧大統領は「ブッシュ米大統領が
北朝鮮問題を善と悪の対立概念で見ているため、説得が難しくなっている。
米国は友邦なので厳しく責めることは出来ないが、
日本とは対決しなければならない」と語ったという。
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

ソース:朝鮮日報
URLﾘﾝｸ(japanese.chosun.com)

680:デフォルトの名無しさん
06/07/18 05:15:28
多変数の4次のルンゲクッタのプログラムを書いているのですが、どうもうまくいきません。

今、厳密解を

y(1) = exp(x) + exp(-x)
y(2) = exp(x) - exp(-x)

とおいて、その微分方程式

dydx(1) = exp(x) - exp(-x)
dydx(2) = exp(x) + exp(-x)

を数値的に解いて、厳密解と比較しています。

サブルーチンのderive中で、

dydx(1) = exp(x) - exp(-x)
dydx(2) = exp(x) + exp(-x)

とおいた場合は、きちんと数値解と厳密解が一致するのですが、

dydx(1) = y(2)
dydx(2) = y(1)

とおくと、両者は一致しません。これはなぜでしょうか？
考えているのですが、一向にわかりません。
プログラムは以下になります。よろしくお願いします。

681:デフォルトの名無しさん
06/07/18 05:17:34
program rk4
implicit none
double precision x0, y0, x, h, y(2), dydx(2), yout(2)
external derive
h = 1.0d-2
x0 = 0.0d0
y(1) = 2.0d0
y(2) = 0.0d0
do x = x0, 2.0d0, h
write(*, *) x, dabs((dexp(x) + dexp(-x)) - y(1))
write(*, *) x, dabs((dexp(x) - dexp(-x)) - y(2))
write(*, *)
call derive(x, y, dydx)
call rk4(y, dydx, x, h, yout, derive)
y(1) = yout(1)
y(2) = yout(2)
end do
stop
end

subroutine derive(x, y, dydx)
implicit none
double precision x, y(2), dydx(2)
c dydx(1) = y(2)
c dydx(2) = y(1)
dydx(1) = dexp(x) - dexp(-x)
dydx(2) = dexp(x) + dexp(-x)
return
end

682:デフォルトの名無しさん
06/07/18 05:18:07
subroutine rk4(y, dydx, x, h, yout, derive)
implicit none
integer i
double precision h, x, dydx(2), y(2), yout(2)
external derive
double precision h6, hh, xh, dym(2), dyt(2), yt(2)
hh = h*0.5d0
h6 = h/6d0
xh = x + hh
do i = 1, 2
yt(i) = y(i) + hh*dydx(i)
call derive(xh, yt, dyt)
yt(i) = y(i) + hh*dyt(i)
call derive(xh, yt, dym)
yt(i) = y(i) + h*dym(i)
dym(i) = dyt(i) + dym(i)
call derive(x+h, yt, dyt)
yout(i) = y(i) + h6*(dydx(i) + dyt(i) + 2.0d0*dym(i))
end do
return
end

683:デフォルトの名無しさん
06/07/18 15:03:58
>>680
うまく行かないのは、subroutine rk4 に誤りがあるから。

以下のループ部分は、i　=　1 について積分し、次にi　=　2について積分しているが、
これはおかしい。i=1,2は独立な量ではなく、相互に依存している。
i=1,2を同時に積分して、ちょっとづつ前進させなければならない。

do i = 1, 2
yt(i) = y(i) + hh*dydx(i)
call derive(xh, yt, dyt)
yt(i) = y(i) + hh*dyt(i)
call derive(xh, yt, dym)
yt(i) = y(i) + h*dym(i)
dym(i) = dyt(i) + dym(i)
call derive(x+h, yt, dyt)
yout(i) = y(i) + h6*(dydx(i) + dyt(i) + 2.0d0*dym(i))
end do

こう直せばおｋ
yt = y + hh*dydx
call derive(xh, yt, dyt)
yt = y + hh*dyt
call derive(xh, yt, dym)
yt = y + h*dym
dym = dyt + dym
call derive(x+h, yt, dyt)
yout = y + h6*(dydx + dyt + 2.0d0*dym)

数値誤差は１０＾－１０程度になるが、これは数値微分を使っているのだからまぁおｋ

684:デフォルトの名無しさん
06/07/18 21:45:16
683さん、お答え頂きありがとうございます！
おっしゃるとおり、サブルーチンを変えたところうまく回りました！
どうもありがとうございました！

いまいち、サブルーチンと配列の組み合わせが理解できないでいます。
例えば、yを配列として(例えばy(3)などと)定義した場合、
yt = y + hh*dydx
で全てのyを計算するのと、
do i = 1, 2
yt(i) = y(i) + hh*dydx(i)
で、ひとつずつyを計算するのは違うのでしょうか？
というのも、ルンゲクッタギルというサブルーチンをひろってきたのですが、
そのコードでは、iを1からnまで計算していています。これを参考に、
間違いを指摘して頂いたコードを作ったのですが・・・。
一体何が違うのでしょうか？お手を煩わせてすみません。
よろしくお願いします。

685:デフォルトの名無しさん
06/07/18 21:46:09
subroutine rkg(derive, x0, y, n, ndim, h, nx, f)
implicit real*8 (a-h, o-z)
external derive
dimension y(ndim, 0:*), f(ndim, 4)
data fp / 1. 70710 67811 86547 52440 d0 /
data fm / 0. 29289 32188 13452 47560 d0 /
data fn / 0. 16666 66666 66666 66667 d0 / ! 1/6
data fp2 / 3. 41421 35623 73095 04880 d0 / ! 2 + sqrt(2)
data fm2 / 0. 58578 64376 26904 95119 d0 /
data fpp / -4. 12132 03435 59642 57320 d0 /
data fmm / 0. 12132 03435 59642 57320 d0 /
hh = 0.5d0*h
fmh = fm*h
fph = fp*h
fnh = fn*h

686:デフォルトの名無しさん
06/07/18 21:47:18
do 100 j = 1, nx
x = dble(j-1)*h + x0
call derivs(x, y(1, j-1), f(1, 3))
do i = 1, n
f(i, 2) = y(i, j-1) + hh*f(i, 3)
end do
call derivs(x + hh, f(1, 2), f(1, 4))
do i = 1, n
f(i, 1) = f(i, 2) + fmh*(f(i, 4)-f(i, 3))
f(i, 2) = fm2*f(i, 4) + fmm*f(i, 3)
end do
call derivs(x + hh, f(1, 1), f(1, 3))
do i = 1, n
f(i, 4) = f(i, 1) + fph*(f(i, 3)-f(i, 2))
f(i, 1) = fp2*f(i, 3) + fpp*f(i, 2)
end do
call derivs(x+h, f(1, 4), f(1, 2))
do i = 1, n
y(i, j) = f(i, 4) + fnh*(f(i, 2)-2.0d0*f(i, 1))
end do
100 continue
return
end

687:デフォルトの名無しさん
06/07/18 22:58:38
>>684
DO　LOOPを消したのは、めんどくさかったからで、あらわに書けば>>685のようになる。

>>685のコードと、元の>>682のコードを見比べてみるとわかるが、
>>685ではｄｅｒｉｖｅを前ごとにループを廻しているのに対して、>>682では全体に対して
ループをまわしている。これでは意味が違ってしまう。

そもそも、もっとも素朴なオイラー法ではy(x+h)=y(x)+hy'(x)なのだが、誤差が大きい。
そこでルンゲ・クッタでは、一気にｈだけ進まないで、まずh/2進んで、そこでのｙをもとめて、
これを元にｙ’をもとめて、さらにまたh/2進むと言う感じになっている。

イメージ的には、オイラーは両足そろえてｈだけジャンプするが、ルンゲクッタは
右足１歩左足一歩、もう一回右足一歩、左足一歩だ。
ところが、>>682のように一気にDOLOOPをまわすと右足ケンケン２回、左足ケンケン２回
的な状況になっていて、なんだかおかしいのだ。
なっている。
（あくまでイメージなので、真面目に考えてくださいｗ）

688:687
06/07/18 23:06:43
日本語狂いまくりんぐｗ

do i = 1, 2
　yt(i) = y(i) + hh*dydx(i)
end do

call derive(xh, yt, dyt)

do i = 1, 2
　yt(i) = y(i) + hh*dyt(i)
end do

call derive(xh, yt, dym)

do i = 1, 2
　yt(i) = y(i) + h*dym(i)
end do

do i = 1, 2
　dym(i) = dyt(i) + dym(i)
end do

call derive(x+h, yt, dyt)

do i = 1, 2
　yout(i) = y(i) + h6*(dydx(i) + dyt(i) + 2.0d0*dym(i))
end do

こういう風に、個々にDoLoopになってればおｋ！　F90なら、DoLoop省略できる。

689:デフォルトの名無しさん
06/07/19 07:56:54
こんにちは、fortran90での課題が手に負えないので助けてください。
「n×nの実数正方行列Ａの行列式｜Ａ｜を再帰手続きでもとめよ」というもの
ですが、再帰手続きの使い方、使いどころが見当もつきません。３×３までは
できたんですが再帰手続きを使っておらず、それ以上に拡張できません。
長いこと詰まっています。どうか教えてください。

690:デフォルトの名無しさん
06/07/19 09:40:39
再帰手続きでってことなら、
一番泥臭い方法を使えって事かな。

先ず、行列式は
　　　　　　　　　 p　　n
｜Ａ｜ = Σ(-1)　P Πa_ii
　　　　　　P　　　　　 i

と書ける。
多分、線形代数の教科書見たら書いてあると思う。

ここで、Π_i^n a_ii は a_11 a_22 ... a_nn という積を表している。
で、P は添字の置換演算子。
a_11 a_22 ... a_nn を a_21 a_12 ... a_nn とか a_11 a_n2 ... a_2n とか、
行列の添字の一方に注目して、それを置換する演算子。
Σ_P は、あらゆる置換を考慮して、全部足せ、ということを表している。
で、p は偶置換だと 0 で、奇置換だと 1 になる。

この「あらゆる置換を考慮して」というところを、
再帰を使って実現すればいい。

691:デフォルトの名無しさん
06/07/19 10:48:05
690さん、ありがとうございます。アドバイスをもとに考えてみます。
fortranをはじめたばかりの初心者なので、不躾なんですがよろしければ
なにか具体的な例を示していただければ幸いです。行列を読み込んで添字に
着目して、部分行列を作ってそれから再帰手続きとしようとしているんですが、
再帰手続きの手順がこんがらがりよく理解できません。どうかお願いします。

692:デフォルトの名無しさん
06/07/19 11:19:36
>>689
690 ではない者ですが、余因子展開すると n×n 行列の行列式を
n-1×n-1 行列の行列式で表せると思います。

693:デフォルトの名無しさん
06/07/19 15:39:04
>>689

MODULE mod_det
IMPLICIT NONE
CONTAINS
RECURSIVE FUNCTION det(a) RESULT(res)
IMPLICIT NONE
REAL (KIND = 8), INTENT(IN) :: a(:, :)
REAL (KIND = 8) :: res, wk(SIZE(a, 1) - 1, SIZE(a, 1) - 1)
INTEGER :: i, n
n = SIZE(a, 1)
res = 0.0d0
IF (n == 1) THEN
res = a(1, 1)
ELSE
DO i = 1, n
wk(:, 1:i - 1) = a(2:, 1:i - 1)
wk(:, i:n - 1) = a(2:, i + 1: n)
res = res + REAL((-1)**(i - 1), KIND = 8) * a(1, i) * det(wk)
END DO
END IF
RETURN
END FUNCTION det
END MODULE mod_det

>>693

694:デフォルトの名無しさん
06/07/19 15:40:57
PROGRAM vip
USE mod_det
IMPLICIT NONE
INTEGER, PARAMETER :: n = 7
INTEGER :: i, j
REAL (KIND = 8) :: chk, d, a(n, n)
DO i = 1, n
DO j = 1, n
a(i, j) = 1.0d0 / REAL(i + j, KIND = 8)
END DO
END DO
d = det(a)
chk = 1.0d0
DO i = 1, n
DO j = 1, n
IF ( i > j ) chk = chk * REAL(i - j, KIND = 8) * REAL(i - j, KIND = 8)
chk = chk / REAL( i + j, KIND = 8 )
END DO
END DO
PRINT *, d, d - chk
STOP
END PROGRAM vip

>>694

余因子展開で解いた。　上の１行を元に展開している。
メインルーチンでは、Ｃａｕｃｈｙの定理でチェックしている。多分大丈夫？

695:デフォルトの名無しさん
06/07/19 15:45:01
＞＞６９３の式では、小行列式を作るのに、あらわにワーク行列を作っているが、
ＦＯＲＴＲＡＮ９０なら、マスキングとかを作って、ワークをあらわに経由しないで、
小行列を引数として再帰関数を呼び出せるような気がする。

手元に紙マニュアルがないので、調べる気がしない。だれか頼む。

でもどっちにしろmask行列を作らないと駄目かな？
意味ないかもｗ

696:デフォルトの名無しさん
06/07/19 15:51:12
687さん、早速お答え頂きどうもありがとうございます！
全体と個々のstepとでloopを回すのでは全く違いますね。
ありがとうございました。

697:デフォルトの名無しさん
06/07/19 16:24:45
692さん、693さん、694さん、695さん、具体例およびアドバイスありがとうございます。
本当に助かります。早速これらのプログラムを参考に自分なりに考えてみます。

698:デフォルトの名無しさん
06/07/20 00:14:15
すいません、fortran77で、
「sin(x)をマクローリン展開し、第ｎ項まで計算したものを、
ｘを0から180度まで30度ずつ変えて表示するプログラムを作成せよ。
但しsin(x)の計算はサブルーチンを使うこと」
という課題が出されたのですが、うまく作れません。
どうかヒントだけでも教えてもらえないでしょうか。

699:デフォルトの名無しさん
06/07/20 03:03:56
>>698
多分マクローリン展開が間違ってるんじゃないかな？

700:デフォルトの名無しさん
06/07/20 03:24:06
>>698
PROGRAM oppai
pi = 4.0 * ATAN(1.0)
DO 10 i = 0, 180, 30
rad = pi * REAL(i) /180.0
CALL sine(10, rad, s)
WRITE(6, *) ' SIN(', i, ') =', s, SIN(rad)
10 CONTINUE
STOP
END
C
SUBROUTINE sine(n, rad, s)
s = 0.0
x = rad
DO 10 i = 1, n
s = s + x
x = - x * rad * rad / REAL( 2 * i * (2 * i + 1) )
10 CONTINUE
RETURN
END

701:700
06/07/20 11:34:13
補足すると、SINをテイラー展開で求めるサブルーチンは、
１番目の引数が第何項まで足すかと言うもの、２番めはSIN（Θ）のΘ、３番目が結果の値。

貼り付けたソースでは１０項まで足しているが、単精度なのでこんなに項を足しても意味無い。
有効桁からより適切な項数を編み出してくれ。

SIN(x)=x-x^3/3!+x^5/5!-x^7/7!+........(-1)^((n - 1)/2)x^n/n! (n odd)

702:デフォルトの名無しさん
06/07/20 12:58:47
>>700
SIN(x)=x-x^3/3!+x^5/5!-x^7/7!
これは偶数次の展開が0になって消えているから7次(または8次)の項まで
展開したという理解で良かったっけ?

703:700
06/07/20 13:50:19
>>702
ういうい。
sin は奇関数だから、偶べきの項は出てこない。ゆえに、偶数次の項は消える。

704:困ってます
06/07/20 16:08:00
2次元のﾎﾟｱｿﾝ方程式

方程式 d2f/dx2 + d2f/dy2 = c

境界値　　f = 0

範囲　 x= -1 ～ 1 , y= 0 ～ 1

定数 c= -1 (x<0) , c = 0 (x=0) , c = 1 (x>0)

離散式 -4*f(i,j)+f(i-1,j)+f(i+1,j)+f(i,j-1)+f(i,j+1)= dx*dx*c

これをvisual Fortranを使ってやるのですが全くわかりません。。。
誰か出来る人いませんか？？？？
あとDOSで開いた後立体的に見たいんですけど
よかったら貼り付けてください！！！

705:デフォルトの名無しさん
06/07/20 16:15:40
>>704
こんなの平均値を求めるのと変わらんぞｗ
本質は１行だｗ
あっちのスレで解いてやったが、今後はこっちで継続せよ。

706:困ってます
06/07/20 16:50:26
あれでやってみたんですが全くできません。。。
なんとか馬鹿でも分かるようにおしえてください！！！！

707:デフォルトの名無しさん
06/07/20 16:57:12
>>706
配列ｆの中に、答えが入っている。
これを出力してグラフを書けばいい。

　＿
/　　\
　　　　\＿/
⇒ｘ軸

　\＿/
⇒ｙ軸　　ｘ＜０

　＿
/　　\
⇒ｙ軸　　ｘ＞０

こんな感じになっている。

708:困ってます
06/07/20 17:10:13
ありがとうございます！！！
でもそれでもいまいち分かってない自分がいます。。。
プログラムをそのままコンパイル出来る形で教えていただけないでしょうか？？？？

709:困ってます
06/07/20 17:13:06
出力したいんですがやり方がやっぱわかないです！！！！！

710:デフォルトの名無しさん
06/07/20 17:56:17
>>709
ごめ、元の奴括弧の範囲を継続行に直したところで、ずらしてしまっていた。
結果は大して変わらんが。訂正しておくｗ
PROGRAM unko
IMPLICIT NONE
INTEGER, PARAMETER :: nx = 19, ny = 9
REAL :: f(-nx - 1:nx + 1, 0:ny + 1), ff(-nx:nx, ny), c(-nx:nx, ny), dx, dy
INTEGER :: i, ix, iy
dx = 1.0 / REAL(nx + 1)
dy = 1.0 / REAL(ny + 1)
c(-nx:-1, :) = -1.0
c( 0, :) = 0.0
c( 1:nx , :) = 1.0
f = 0.0
f(-nx - 1, :) = 0.0
f( nx + 1, :) = 0.0
f(:, 0) = 0.0
f(:, ny + 1) = 0.0
DO i = 1, 100
DO ix = -nx, nx
DO iy = 1, ny
ff(ix, iy) = ( -c(ix, iy) * dx * dy &
+ f(ix - 1, iy - 1) + f(ix - 1, iy + 1) + f(ix + 1, iy - 1) + f(ix + 1, iy + 1) ) / 4.0
END DO
END DO
f(-nx:nx, 1:ny) = ff
END DO
DO ix = -nx - 1, nx + 1
WRITE(9, '(1000es11.3)') dx * ix , f(ix, :)
END DO
STOP
END PROGRAM unko

711:デフォルトの名無しさん
06/07/21 00:18:35
教えてください！！

program oppai
implicit none integer:: saikou,saitei,mark,n
n=0
saikou=0
saitei=100
do
read *,mark
if (mark<0) exit
n=n+1
if (mark<=saitei) then
saitei=mark
end if
if (mark>=saikou) then
saikou=mark
end if
end do
print *,'最高点',saikou
print *,'最低点',saitei
pause
end

このプログラムは、生徒のテストの点数を次々に読み込み、
負の数が入力されたら入力を終了し、その最高点、最低点を出力する
プログラムなんですが、100点満点でないテストでは、うまく動きません。
テストが100点満点かどうかわかっていない場合でも、動くように改良してほしいの
です。。。

お願いします～

712:デフォルトの名無しさん
06/07/21 00:22:43
自然対数の底 e を100桁計算するプログラムを教えてください、
e = 1/0!+1/1! + 1/2! +1/3!+… となります
適当なMを決めて
e = 1/0!+1/1! + 1/2! +1/3!+…+1/M! を計算してやればいいみたいです

ヒントがあります。

これを次のように変形する。

e=(.....(1/M+1)/(M-1)+1)/(M-2)+1)....)/2+1)/1+1
M=80で充分
80! > 10^110
(n! ～ (n/e)^n*sqrt(2*pi*n))
おおまかな求める手順 1.e=1でスタート
2.k=M,M-1,...,1について以下を繰り返す
e←e/k+1 とする

100桁求めるにはM=80で十分

100桁の筆算を計算機上でおこなう
e(0) に整数部分
e(1) に小数第１桁
e(i) に小数大 i 桁

と考え、割り算のプログラムを作る
この100桁の数にたいして kで割るということができればよい。

713:デフォルトの名無しさん
06/07/21 00:23:15

1.234 / 4
-> 1 / 4 --> 0 あまり 1
-> (1*10+2)/4 --> 3 あまり 0
-> (0*10+3)/4 --> 0 あまり 3
-> (3*10+4)/4 --> 8 あまり 2

----> .308 あまり 0.002

e(0) に整数部分
e(1) に小数第１桁
e(i) に小数大 i 桁

e/k (kは1～M=100)

e=(e(0).e(1)e(2)....) をk でわる
amari=0
do i=0,100
jo=amari*10+e(i)
e(i)=jo/k
amari=mod(jo,k)
end do

計算終了後、求めたeの値は次のようにして出力する。

print *, e(0),"."
print '(" ",50I1)',(e(i),i=1,100)
pause
end

お願いします！

714:デフォルトの名無しさん
06/07/21 10:04:01
>>711

saitei=100

これがあるせいで、１００点以上にならない。

saitei = 2**30

こうしておけばおｋ！

715:デフォルトの名無しさん
06/07/21 14:26:10
>>712
これだけヒントをもらっているんだから、もう少しがんばれｗ
かったるいがF77で書いてやった。
１０分で作ったから間違ってても知らんｗ

PROGRAM EXPO70
PARAMETER (ne = 100, n = 80)
INTEGER ie(ne)
ie(1) = 1
DO 10 i = 2, ne
ie(i) = 0
10 CONTINUE
DO 20 m = n, 1, -1
k = ie(1)
DO 30 i = 1, ne - 1
ie(i) = k / m
k = 10 * MOD(k, m) + ie(i + 1)
30 CONTINUE
ie(1) = ie(1) + 1
20 CONTINUE
WRITE(*, '(i1, a1, 99i1)') ie(1), '.', (ie(i), i = 2, ne)
STOP
END

URLﾘﾝｸ(www.math.utah.edu)
値はここで比較しろ。最後の桁が違って居るが、これはたぶんおｋ。

716:715
06/07/21 14:39:13
微妙に改良＋ワーク領域を１０１こ取れば１００桁目まで正しくなる。

PROGRAM EXPO70
PARAMETER (ne = 100, n = 80)
INTEGER ie(ne + 1)
DO 10 i = 1, ne + 1
ie(i) = 0
10 CONTINUE
DO 20 m = n, 1, -1
ie(1) = ie(1) + 1
k = ie(1)
DO 30 i = 1, ne
ie(i) = k / m
k = 10 * MOD(k, m) + ie(i + 1)
30 CONTINUE
20 CONTINUE
ie(1) = ie(1) + 1
WRITE(*, '(i1, a1, 99i1)') ie(1), '.', (ie(i), i = 2, ne)
STOP
END

717:デフォルトの名無しさん
06/07/21 15:42:29
>>716
うは、ne=10000 n=3250で１００００桁までもとまりんぐｗ

718:デフォルトの名無しさん
06/07/22 03:43:44
横35、縦35のマスを考える。そのマス目なかにグループAの人が何人か、グループBの人が何人か、いるものとする。
（各マス目の中には、最大1人しか入ることができない。また、誰もいないマス目もある。）
各人は自分の回りに同じグループの人が多くいた方が幸福だと思っている。そこで「幸福度」を次のように定義する。
自分の場所の上下左右、斜め上下左右、合計8箇所を「自分の回り」とする。この中にいる自分と同じグループにいる仲間の数をMとする。
また、「自分の回り」にいる全ての人の数をNとする。この時、幸福度を
・N=0　なら幸福度 0
・N>0 なら幸福度 M/N
と定義する。

さて、各人は自分の幸福度がある値（以後、しきい値と呼ぶ）より小さいと、現状に対する不満から、ランダムに空いているマス目に移動をしようとする。
ただし、移動しない可能性もある。例えば、上図のように、ある人の「自分の回り」のマス目のうち、
上、下、右斜め下、左斜め上の4箇所が空いており、他のマス目には他の人がいるとする。
この時、この人はそれぞれ確率1/5で
・上に移動
・下に移動
・右斜め下に移動
・左斜め上に移動
・移動しない
という行動をとる。35×35のマスの中には多くの人がいるが、この移動は
・まず、一番上のマスの一番左の人が（もしいれば）移動を試みる。（もちろん、この人の幸福度がしきい値以上なら、移動しない）
・次に一番上のマスの左から2番目の人が移動を試みる。以降、順にその右の人が移動を試みる。
・一番上のマスにいる人の移動が終了したら、次は上から二番目のマスの一番左の人が移動を試みる。
・以下同様。

次ページ