matlab — Fortran COMPLEX вычисляет отличия от переполнения стека

Question

matlab — Fortran COMPLEX вычисляет отличия от переполнения стека

Я завершил порт с Fortran на C ++, но обнаружил некоторые различия в типе COMPLEX. Рассмотрим следующие коды:

PROGRAM CMPLX
COMPLEX*16 c
REAL*8 a
c = (1.23456789, 3.45678901)
a = AIMAG(1.0 / c)
WRITE (*, *) a
END

И C ++:

#include <complex>
#include <iostream>
#include <iomanip>

int main()
{
std::complex<double> c(1.23456789, 3.45678901);
double a = (1.0 / c).imag();

std::cout << std::setprecision(15) << " " << a << std::endl;
}

Компилируя версию C ++ с помощью clang ++ или g ++, я получаю вывод: -0.256561150444368
Однако компиляция версии на Фортране дает мне: -0.25656115049876993

Я имею в виду, разве оба языка не следуют IEEE 754? Если я запускаю следующее в Octave (Matlab):

octave:1> c=1.23456789+ 3.45678901i
c =  1.2346 + 3.4568i
octave:2> c
c =  1.2346 + 3.4568i
octave:3> output_precision(15)
octave:4> c
c =  1.23456789000000e+00 + 3.45678901000000e+00i
octave:5> 1 / c
ans =  9.16290109820952e-02 - 2.56561150444368e-01i

Я получаю так же, как версия C ++. Что случилось с типом Fortran COMPLEX? Я пропускаю некоторые флаги компилятора? -fast-math ничего не меняет. Я хочу создать точно такие же 15 десятичных знаков в C ++ и Fortran, чтобы легче было различать портирование.

Есть ли гуру Фортрана? Спасибо!

2

c++fortran matlab numerical-computing

Решение

Другие решения

Я бы предложил изменить все константы Фортрана на DP

1.23456789_8 (или 1.23456789D00) и т. Д.

и использовать DIMAG вместо AIMAG

0

Источник

Accepted Answer

В коде Фортрана заменить

c = (1.23456789, 3.45678901)

с

c = (1.23456789d0, 3.45678901d0)

Без kind реальные литералы, которые вы используете для rhs, — это, скорее всего, 32-битные реалы, и вам, вероятно, нужны 64-битные реалы. Суффикс d0 заставляет компилятор создавать 64-битные реалы, наиболее близкие к указанным вами значениям. Я вкратце остановился на некоторых деталях этого, и есть другие (возможно, лучшие) способы указания типа литерала действительного числа, но этот подход должен работать нормально на любом текущем компиляторе Фортрана.

Я не очень хорошо знаю C ++, я не уверен, что у кода C ++ такая же проблема.

Если я правильно прочитал ваш вопрос, два кода дают одинаковый ответ на 8sf, предел одинарной точности.

Что касается соответствия IEEE-754, этот стандарт, насколько мне известно, не охватывает вопросы сложной арифметики. Я ожидаю, что арифметика f-p, используемая за кулисами, в большинстве случаев дает результаты по комплексным числам в пределах ожидаемых границ ошибок, но я не знаю, что они гарантированы, как и границы ошибок по арифметике f-p.

5