Использовать ПОТОЛОК без эффекта ошибки округления

ProgramBox

Использовать ПОТОЛОК без эффекта ошибки округления

Post author:admin
Запись опубликована:7 февраля, 2023
Post category:Вопросы по программированию

#fortran #fortran90 #intrinsics

#fortran #fortran90 #встроенные

Вопрос:

Я пытаюсь использовать встроенную функцию ‘ПОТОЛОК’, но ошибка округления иногда затрудняет получение того, что я хочу. Пример кода очень прост:

 PROGRAM MAIN

  IMPLICIT NONE

  INTEGER, PARAMETER        ::      ppm_kind_double = KIND(1.0D0)
  REAL(ppm_kind_double)     ::      before,after,dx

  before = -0.112
  dx = 0.008
  after = CEILING(before/dx)

  WRITE(*,*) before, dx, before/dx, after

  END

И я получил результаты:

Значение, которое я даю ‘before’ и ‘dx’ в коде, предназначено только для демонстрации. Например, для тех, кто был до /dx = -13,5, я хочу использовать ПОТОЛОК, чтобы получить -13. Но для изображения, которое я показываю, я действительно хочу получить -14. Я рассмотрел возможность использования некоторых аргументов, таких как

  IF(ABS(NINT(before/dx) - before/dx) < 0.001)

Но это просто некрасиво. Есть ли лучший способ сделать это?

Обновить:

Я был удивлен, обнаружив, что проблема не возникнет, если я установлю переменные в константы в ppm_kind_double . Поэтому я предполагаю, что эта «ошибка округления» произойдет только тогда, когда количество цифр для точности округления используемой мной машины больше, чем определено в ppm_kind_double . На самом деле я запускаю свою программу (не этот демонстрационный код) в кластере, о точности которого я не знаю. Так что, может быть, именно квадратичная точность на этой машине приводит к проблеме?

После того, как я установил константы с двойной точностью:

 before = -0.112_ppm_kind_double
dx = 0.008_ppm_kind_double

1. Вы немного сбили меня с толку, желая -14. Мне потребовалось некоторое время, чтобы понять, что -0.112 / 0.008 — это точно -14, но из-за ошибок округления результат немного больше, и CEILING then усиливает эту незначительную ошибку.

2. @chw21, на самом деле я видел ваш предыдущий комментарий и ждал, когда вы его найдете. 🙂

3. Учитывая, что точность здесь является частью проблемы, стоит отметить, что ваши литеральные константы имеют только одинарную точность — например -0.112 , вероятно, должны быть записаны -0.112_ppm_kind_double . Однако это изменение не решит основную проблему.

4. В дополнение к созданию констант согласованного типа данных, трудно догадаться, почему вы не использовали ANINT вместо CEILING .

5. Кстати, ваша интерпретация неверна. На самом деле произошло то, что изначально вы устанавливали для своих переменных неточную одинарную точность, а с суффиксом вида вы используете правильную точность. Я настоятельно рекомендую вам прочитать о числовой точности выражений в Fortran. Числа по умолчанию имеют одинарную точность. И ваша «точность округления машины» не существует.

Ответ №1:

Это немного сложно, потому что вы никогда не знаете, откуда берется ошибка округления. Если dx бы он был чуть больше, чем 0.008 тогда, деление before/dx все равно могло бы быть округлено до того же значения, но теперь -13 это был бы правильный ответ.

Тем не менее, наиболее распространенный метод, который я видел, — это просто немного подтолкнуть предыдущее значение в противоположном направлении. Что-то вроде этого:

 program sign_test
    use iso_fortran_env
    implicit none
    real(kind=real64) :: a, b
    integer(kind=int32) :: c
    a = -0.112
    b = 0.008
    c = my_ceiling(a/b)
    print*, a, b, c
contains
    function my_ceiling(v)
        implicit none
        real(kind=real64), intent(in) :: v
        integer(kind=int32) :: my_ceiling
        my_ceiling = ceiling(v - 1d-6, kind=int32)
    end function my_ceiling
end program sign_test

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

Ответ №2:

обратите внимание, что если ваши значения условно «точны» с указанной точностью, вы можете сделать что-то вроде этого:

  after=nint(1000*before)/nint(1000*dx)

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