LoginSignup
0
0

More than 1 year has passed since last update.

fortran総称名の罠: cmplx,real,imagは単精度への型変換

Last updated at Posted at 2023-03-09

https://www.nag-j.co.jp/fortran/FI_10.html
には「総称名を使うことを強く推奨する」といいつつヤバイ例が載せてある。
それは以下のコードです。ヤバイ点はわかります?危険なんです。

program generic_vs_specific
  implicit none
  double precision :: a = 2.0
  complex(kind(0d0)) :: b = (2.0,1.0)
  print *, sqrt(a)
  print *, sqrt(b)
end program generic_vs_specific

complex(kind(0d0)) :: b = (2d0,1d0)としないとまずいです。complex(kind(0d0)) :: b = (2,1)でもいいですけどいつもd0つけると覚えとくのが良い。できるだけ脳の無駄遣いをしなくていいコードを書きましょう。上のコードでa=2.1,b=1.1として

  print *, a
  print *, b

を追加して実行するとわかるけど、a,bは単精度しか出ません。a=2.0,b=1.0だと一見精度が出ているわけでとても危険な例です。
数値計算は「丸め誤差と離散化誤差の狭間での戦い」なので状況によりヤバイことにもなりえます。

上の例は比較的トリビアルだけど、もっと危険な例があります。以下のコードを実行してみてください。

  complex(kind(0d0)) :: ccc = (5.1,3.2),c,cc,cc1,cc2,cc0,cc3
  real(8)::a,b
  write(6,*) 'ccc=',ccc
  a=5.1d0
  b=3.2d0
  cc1=cmplx(a,b)
  write(6,*) 'cc1=',cc1
  cc2=dcmplx(a,b)
  write(6,*) 'cc2=',cc2
  cc3=complex(a,b)
  write(6,*) 'cc3=',cc3
  end

倍精度が保てるのはcc2,cc3です。cc1も単精度ですー「総称名を使いましょう」の罠です。
real(a)と書くと当然単精度への変換となるわけだけど、cmplxだと複素化する機能だからね、と錯覚しちゃうわけです。
とにかく「cmplxは使ったらいかん」ということは記憶しておかないといけない。これは精度を単精度に落とす関数なんです。(realやimagも「実数部分、虚数部分をとるだけではない」ことには注意しとく必要はある)。

gnufortranだとhttps://gcc.gnu.org/onlinedocs/gcc-3.4.6/g77/Complex-Intrinsic.html#Complex-Intrinsicのcomplex
という関数をつかえなくはないですーこの最後の方に説明を書いています。がfortranの「標準実装」ではないようで、ぼくのifortバージョンでは実装されてないです(ifortとgfortranで動けば「標準実装」としていますーnvfortranもやったほうがいいとは思ってはいるけれど)。

丸め誤差と離散化誤差:
数値計算は「丸め誤差と離散化誤差の狭間での戦い」なので状況によりヤバイことにもなりえます。普通に微分を取るとまあ3桁ぐらいは桁落ちするー2階微分なら6桁ぐらいって感じだとすると8桁しかない単精度だと精度が2桁しか残らない。100万回ぐらい計算するとその最悪のやつが悪さをして全体の精度に響くということもあり得るー正定値加算など問題が生じない場合もあるけれど、負符号も取りうる量を足すときには悩まされます。基本として入力に対する出力のスムーズネスをみれば数値的偶然誤差は計測できるけれども数値的系統誤差は制御できない.
また数値的偶然誤差を減らすには、4倍精度計算を入れていって様子をみるしかないーよほど賢くやらないと、なかなかにタフな仕事になる。

台形公式とシンプソン公式ーシンプソン公式なら倍精度は必須(がある種のルールオブサムだと思う)。あるいは丸め誤差が積算してる状況なら台形公式のほうが計算安定性がよい。そういう情報が過去ログになりすぎてるような気がするな。正定値積算なのか、なめらかな関数の積算なのか、などをおよそ想像して公式を選ばないといけない。

駄:ずっとむかし、akaikkrにこのcmplxのバグを発見したときヒトの未来予測能力の限界を感じました。たぶんfortranの当初は総称名(ポリモルフィズム)というアイデアはなかったと思うーほんの10年先のごくごく当たり前なはずの未来が読めない。いまの「スマホ社会」を西暦2000年に予言したヒトはいなかった。OLEDは3000時間持てば実用素材だけど、昔は50000時間の寿命がないと実用にならんと言われていた。Apple最大の発明は電池が入れ替えれないってことでしょう?これも驚天動地でした。まあだからこそ未来が楽しみってことでもあるけれど。

0
0
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
0
0