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

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

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

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

    function genrand_real()

    real   genrand_real
    real*8 genrand_real2

    genrand_real = real(genrand_real2())

    return
    end

Я использую real и real*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

random rounding fortran mersenne-twister fortran77
person Henrique M. Cezar    schedule 15.06.2016    source источник

Ответы (1)


arrow_upward
1
arrow_downward

Опубликованная вами функция НЕ генерирует случайные числа, она ограничивает только случайные целые числа (от genrand_int32()) до интервала [0,1) путем деления на 2 ^ 32 (что точно равно 4294967296) или сначала добавляет 2 ^ 32, если int отрицательное . 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 и так далее. Впрочем, это скорее академический момент.

person StefanS    schedule 15.06.2016
comment
Большое спасибо! Вы не только решили мою проблему, но и сделали мою программу намного быстрее, так как: mt19937 истек: 11.5120001, пользователь: 11.5120001, sys: 0.00000000 mt19937 одинарная точность истек: 4.83599997, пользователь: 4.83599997, sys: 0.000000 - person Henrique M. Cezar; 15.06.2016
comment
Пожалуйста. Это эффект просто не использования более высокой точности. Однако есть и обратная сторона, см. stackoverflow.com/a/17951021 (это для C#, но концепция та же). Однако, если остальная часть вашей программы в любом случае реальна с одинарной точностью, это вряд ли имеет значение. - person StefanS; 15.06.2016
comment
Еще раз, спасибо за вашу помощь и извините за то, что я был тем парнем, но у меня есть дополнительный вопрос: stackoverflow.com/questions/37859027/ - person Henrique M. Cezar; 16.06.2016
comment
Прочитав комментарии к другому вопросу: да, я должен был подумать о проблеме, заключающейся в том, что стандартные целые числа не могут точно содержаться в вещественных числах с одинарной точностью. Я подумаю об этом. - person StefanS; 16.06.2016
comment
Я не понимаю. Гарантирует ли real(integer) округление до 0 для чисел, которые не могут быть представлены, а real(double) округляет до ближайшего? - person sh1; 18.06.2016