fortran & C：計算結果が毎回異なる

Question

fortranでcos(x)をYからXの範囲で台形積分するプログラムをC言語のプログラムから呼び出すコードを書いています。
今回の質問はそのfortranの部分についてで、毎回異なった計算結果が表示されるのですが、fortran初心者すぎて何が理由かわかりません。以下がそのfortranのコードです。また、このコードを呼び出すCプログラムを一番下に置いておきます。

入力はプログラムの実行時にコマンドライン引数で任意の２数を指定して、出力は入力した２数の範囲でcos(x)を台形積分した結果として正しい値が帰ってくることを期待しています。
実際に実行したところ、実行する度に違う結果が出ました。(以下は入出力例です。)

./JikkouFile_Name 1.5707963 0

b=1.570796
a=0.000000
   1.5707963000000000        0.0000000000000000        4.59149455E-41
  -469083744  -469083744.00000000     
-469083744.000000

      DOUBLE PRECISION FUNCTION SEKIBUN(X,Y)
      DOUBLE PRECISION X
      DOUBLE PRECISION Y
      INTEGER I
      DOUBLE PRECISION WIDTH=(X-Y)/10000000
      DOUBLE PRECISION MENSEKI=0.0
      PRINT *,X,Y,WIDTH
      DO I=0,10000000
       MENSEKI=MENSEKI+DCOS(I*WIDTH+Y)*WIDTH
      END DO
      SEKIBUN=MENSEKI
      PRINT *,MENSEKI,SEKIBUN
      RETURN
      END

#include <stdio.h>
#include <stdlib.h>

extern double sekibun_(double *, double *);

int main(int argc, char **argv){
    double a,b,menseki;
    b = atof(argv[1]);
    a = atof(argv[2]);
    printf("b=%lf\na=%lf\n",b,a);    //(*1*)
    menseki = sekibun_(&b, &a);
    printf("%lf\n",menseki);         //(*2*)
    return 0;
}

Comments

[otn]

「こういう入力に対して、こういう結果を期待したが、こう出力された」という風に書いてください。

[izmktr]

一度、FORTRAN側の処理を渡した引数を表示するだけとか、
FORTRAN側が固定値を返してCで表示するだけとか、そういう処理で確認してみるのはどうでしょうか
毎回変わるなら、Cがポインタで引数を渡しているのが原因なんじゃないかと思いますが…

[wanwannyaan]

FORTRANの出力系を調べてもよくわからなくて、やっていなかったのですが、再度調べ直したところ理解できたのでやってみました。
やってみた結果、C→fortran, fortran→Cのデータの受け渡しで数値が変わっていることはなかったので、fortran側に間違いがあるようです。
とりあえず、これから質問内容に反映します。

Answer 1

MENSEKI が save 変数になっているので、前回の値が残っていて、どんどん加算されてゆく。

DOUBLE PRECISION MENSEKI=0.0

短冊形の積分をするのに 10000000 分割するのも理解しがたい。しかも端点の処理がされていないので、いくら細かく分割したところで期待している値には収束しないだろう。中点積分の方がまし。

もっと少ない分割数で、分割数を変えて計算して、収束の度合いを見る方がはるかに有益だろう。あと２進法コンピュータなら 2 のべき乗で分割しておいた方が色々都合がいいことが多い。

integer, parameter :: n = 128 

初心者ということで、仕方ないとは思うが FORTRAN77 以前の古い規格の書き方をしつつ、五月雨式に Fortran90 以降の機能を使っているので分かり難い。

まず implicit none はつけたほうがいい。

77 までの正規の名称ではあるが、DOUBLE PRECISION を使うのは避ける方がどちらかといえば賢明。倍精度だからと言って DCOS を使うのもよくない。77 以降では総称名の COS の方が一般的。

FORTRAN77 以前を主たる書き方にするならは、暗黙の型を守るのは絶対なので MENSEKI という変数名もよくない。鼻の穴から牛乳飲むレベルの異常感。

C 言語から呼ぶなら Fortran 2003 以降の機能になるが、bind 属性をつければ、そこで指定した名前で C 側から見える。アンダースコアなどは、ベンダー毎に代わっている。

DOUBLE PRECISION FUNCTION SEKIBUN(X,Y) bind(c, name = 'sekibun')

Answer 2

WIDTHの表示値がおかしいですね。
変数宣言と、代入を分けてみては？

Comments

[wanwannyaan]

WIDTHの部分、
DOUBLE PRECISION WIDTH
WIDTH=(X-Y)/10000000
のように修正したところ、出力が次のようになり、WIDTHの値は正しく計算されているようです。また、積分値は直りませんでしたが（出力1.000000を期待している）、何度やっても同じ結果が表示されるようにはなりました。ありがとうございます。
b=1.570796
a=0.000000
1.5707963000000000 0.0000000000000000 1.5707963000000001E-007
32766 32766.000000000000
32766.000000

[otn]

規格を確認していませんが、宣言文の式右辺の評価は、呼び出し前に行われるのでは？