Округление двойной точности до одинарной точности: форсирование верхней границы

Question

Я использую реализацию Mersenne Twister, которая предоставляет мне числа с двойной точностью.

http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/VERSIONS/FORTRAN/fortran.html (реализация в Fortran 77 Цуёси Тада, я использую genrand_real2)

Однако мое приложение нуждается во избежании предупреждений при умножении чисел с разной точностью на случайное число одинарной точности. Итак, я написал небольшую функцию для преобразования между двумя типами данных:

    function genrand_real()

    real   genrand_real
    real*8 genrand_real2

    genrand_real = real(genrand_real2())

    return
    end

Я использую реальные и реальные *8, чтобы соответствовать коду, над которым я работаю. Он прекрасно работает большую часть времени (кроме того факта, что я не уверен, насколько быстрым является real()), однако он меняет верхнюю границу моего ГСЧ, поскольку преобразование изменяет [0,1) на [0,1]. Я никогда не думал об этом, пока у меня не было проблем с этим.

Мой вопрос заключается в том, как я могу эффективно обеспечить верхнюю границу или даже как написать функцию, аналогичную genrand_real2 (оригинальной), которая предоставляет мне реал с одинарной точностью. Я думаю, мне нужно только заменить делитель 4294967296.d0, но я не знаю, по какому количеству

  function genrand_real2()

  double precision genrand_real2,r
  integer genrand_int32
  r=dble(genrand_int32())
  if(r.lt.0.d0)r=r+2.d0**32
  genrand_real2=r/4294967296.d0

  return
  end

Answer 1

Размещенная вами функция НЕ генерирует случайные числа, она ограничивает только случайные целые числа (из genrand_int32()) к интервалу [0,1) путем деления на 2^32 (что в точности равно 4294967296) или добавления сначала 2 ^ 32, если целое число отрицательно. 2^32 - это число значений, которое может содержать стандартное целое число, наполовину отрицательное, наполовину положительное (приблизительно, в положительном конце отсутствует 1) и, следовательно, происходит от функции genrand_int32(),

Представьте, что у вас есть числа от -10 до 10 и вы хотите ограничить их интервалом [0,1]. Самое простое решение состоит в том, чтобы добавить 20 к отрицательным числам (таким образом, положительные остаются 0-10, а отрицательные становятся 10-20), а затем делятся на 20. Это именно то, что делает функция, просто с 2^31 вместо 10.

Если вам интересно, почему интервал для вашей функции равен [0, 1): поскольку для числа 0 также требуется место, а битовое представление может хранить только 2 ^ 32 числа, вы не можете иметь 2^31 отрицательных и 2^31 положительное число И 0. Решение состоит в том, чтобы опустить значение +2^31 (самое высокое положительное) и, следовательно, 1 исключается из вашего интервала.

Итак, чтобы свести все это к единой точности:

function genrand_real2()

real genrand_real2,r
integer genrand_int32
r=real(genrand_int32())
if(r.lt.0)r=r+2**32
genrand_real2=r/4294967296

return
end

Магические числа должны оставаться неизменными, потому что они относятся к целым числам, а не к действительным числам.

Изменить: Вы уже сказали это сами, так что я просто повторяю для других людей: Для переносимости технически не рекомендуется использовать типы по умолчанию без указания точности. Так что вы должны сделать sp = selected_real_kind(6, 37) (sp для одинарной точности) где-то, а затем real(kind=sp)... а также 2.0_sp и так далее. Однако это скорее академический момент.