Численные рецепты в предположении целочисленной модели fortran 90 ran_init

Question

ran_init Подпрограмма Числовых Рецептов содержит следующие строки:

INTEGER(K4B) :: new,j,hgt
...                                                                                            
hgt=hg 
...
if (hgt+1 /= hgng)    call nrerror('ran_init: arith assump 3 fails')

куда K4B, hgng а также hg глобально объявлены в модуле через:

INTEGER, PARAMETER :: K4B=selected_int_kind(9) 
INTEGER(K4B), PARAMETER :: hg=huge(1_K4B), hgm=-hg, hgng=hgm-1 

Проблема в том, что в одном конкретном компьютере (но не в других) я получаю ошибку ran_init: arith assump 3 fails, Единственное, что я получаю из документации об этой ошибке:

Немного грязного белья здесь! Мы проверяем, переходит ли самое положительное целое число hg в самое отрицательное целое число hgng, когда к нему добавляется 1. Мы не можем просто написать hg+1, так как некоторые компиляторы оценивают это во время компиляции и возвращают сообщение об ошибке переполнения. Если ваш компилятор видит сквозь шараду временную переменную hgt, вам придется найти другой способ ее обмануть!

Как я могу обмануть это?

Answer 1

Непосредственная причина аварии:

компилятор на этапе оптимизации может легко доказать, что hgt+1 а также hgng не может быть таким же, и поэтому он оптимизирует все условия для .false., Добавление двух положительных целых чисел не может создать отрицательное целое число.

Вы могли бы избежать этого, используя меньший уровень оптимизации или какой-то маскарад с временными переменными, которые были бы более сложными, чем то, что они делали в старых Числовых Рецептах. Вы можете попробовать использовать -fwrapv или же -fno-strict-overflow флаг. Но это крайне неуверенно и будет работать только в данной версии компилятора с определенными флагами компилятора.

Самое простое "исправление" - это, вероятно, просто удалить оскорбительную проверку. Вполне вероятно, что это сработает во многих обстоятельствах (и ужасно провалится в других)

Объяснение и решение:

То, что делают Числовые Рецепты, противоречит стандарту Фортрана. Они предполагают, что если вы добавите 1 к наибольшему целому числу, вы получите самое низкое.

Это не разрешено для целых чисел Фортрана и не разрешено для целых чисел со знаком в C. Это разрешено только для целых чисел без знака в C, и Fortran не имеет их.

Таким образом, компилятор может с уверенностью предположить, что добавление двух положительных целых чисел никогда не приведет к отрицательному целому числу. См. https://gcc.gnu.org/bugzilla/show_bug.cgi?id=30475 для подробного обсуждения этих оптимизаций.

Один из вариантов - переписать операции, которые могут привести к переполнению в C, и использовать преобразование в целые числа без знака. Я использую это в генераторе случайных чисел, который я использую (основанный на http://www.cmiss.org/openCMISS/wiki/RandomNumberGenerationWithOpenMP):

#include <stdint.h>
#include <memory.h>

int32_t sum_and_overflow (int32_t a, int32_t b)
{
  uint32_t au, bu, su;
  int32_t s;

  memcpy(&au, &a, sizeof(a));
  memcpy(&bu, &b, sizeof(b));
  su = au + bu;
  memcpy(&s, &su, sizeof(s));
  return s;
} 

с интерфейсом Fortran

interface
  function sum_and_overflow(a, b) result(res) bind(C, name="sum_and_overflow")
    use, intrinsic :: iso_c_binding
    integer(c_int32_t) :: res
    integer(c_int32_t), value :: a, b
  end function
end interface

и вместо

c = a + b

Я звоню

c =  sum_and_overflow(a, b)

Так что в вашем случае

if (sum_and_overflow(hgt,1) /= hgng) ...

но не только там, вам потребуется хотя бы еще одно место, где это предположение используется в генераторе. В используемом мной генераторе есть только одна такая строка.

---

Есть много других хаков, которые могут привести к временному успеху, но позже потерпят неудачу с некоторыми другими параметрами компилятора. Например, неопределенная очистка поведения в GCC не любит переполнение в целых числах со знаком в Fortran и C и завершит вашу программу, если вы это сделаете.

Вот почему я пытаюсь показать решение, которое является более сложным, но работает не по стандарту, а по нему.