periodic2 Ewald の k≠0 項 exp(kz)erfc(karg0 + adz) を overflow させずに返す。 kz = ±|k|dz と adz = ±alphadz は同じ符号で、karg0 = |k|/(2 alpha) なので kz - (karg0 + adz)^2 = -(karg0^2 + adz^2) が成り立つ。 exp(kz) が binary64 で有限な kz では元の積をそのまま返す(従来の式とビット一致)。 exp(kz) が overflow する kz では y = karg0 + adz >= sqrt(2 kz) > 37 なので erfc は 0 に underflow し、 元の積は inf0 = NaN になる。この範囲だけ exp(-(karg0^2 + adz^2))erfc_scaled(y) で評価する。 真値はこのとき exp(-kz) 以下で、binary64 では 0 か非正規数になる。
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| real(kind=dp), | intent(in) | :: | kz |
指数の引数 ±|k|dz。 |
||
| real(kind=dp), | intent(in) | :: | karg0 |
指数の引数 ±|k|dz。 |k|/(2 alpha)。 |
||
| real(kind=dp), | intent(in) | :: | adz |
指数の引数 ±|k|dz。 |k|/(2 alpha)。 kz と同符号の ±alpha*dz。 |
pure real(dp) function periodic2_ewald_exp_erfc(kz, karg0, adz) result(term) real(dp), intent(in) :: kz, karg0, adz real(dp), parameter :: max_finite_exp_arg = log(huge(1.0_dp)) if (kz <= max_finite_exp_arg) then term = exp(kz)*erfc(karg0 + adz) else term = exp(-(karg0*karg0 + adz*adz))*erfc_scaled(karg0 + adz) end if end function periodic2_ewald_exp_erfc