Intrinsic-функции и управление точностью

Встроенные функции Fortran работают с целыми массивами разом, а параметр kind даёт точный и переносимый контроль над точностью вещественных чисел.

Intrinsic-функция — встроенная в язык процедура (например sum, sqrt, matmul), часто работающая поэлементно или над массивом целиком. kind — целочисленный параметр типа, задающий разрядность (точность и диапазон) числовых типов переносимым образом.

Почему встроенные функции — это про производительность

Fortran предоставляет богатую библиотеку встроенных функций, и они не просто удобны — они быстры. Многие из них принимают массивы и работают со всеми элементами сразу. sqrt(x) для массива x извлекает корень из каждого элемента; sum(a) складывает весь массив; maxval(a) находит максимум. Это не сахар: компилятор знает семантику встроенных функций и оптимизирует их агрессивно — векторизует, разворачивает циклы, использует специальные инструкции процессора. Ваш ручной цикл может оказаться медленнее, чем sum(a), потому что в встроенную функцию вложены годы оптимизации. Поэтому правило хорошего тона: если есть встроенная функция для задачи — используйте её, а не пишите цикл руками.

Маски и условные редукции

Особая сила встроенных функций — работа с масками. Многие редукции принимают необязательный аргумент mask — логический массив той же формы, и учитывают только элементы, где маска истинна. Это позволяет выражать условные вычисления без явных циклов и ветвлений.

program intrinsics_demo
  implicit none
  real :: a(8) = [3.0, -1.0, 4.0, -1.0, 5.0, -9.0, 2.0, 6.0]

  print *, "Сумма всех:        ", sum(a)
  print *, "Сумма положительных:", sum(a, mask = a > 0.0)
  print *, "Сколько отрицательных:", count(a < 0.0)
  print *, "Максимум:           ", maxval(a)
  print *, "Индекс максимума:   ", maxloc(a)
  print *, "Все ли > -10:       ", all(a > -10.0)
end program intrinsics_demo

Вывод:

 Сумма всех:            9.00000000    
 Сумма положительных:   20.0000000    
 Сколько отрицательных:           3
 Максимум:            6.00000000    
 Индекс максимума:            8
 Все ли > -10:        T

Выражение sum(a, mask = a > 0.0) складывает только положительные элементы — компактно и без цикла. Заметьте: maxloc возвращает массив индексов (для многомерного случая — по одному на измерение), поэтому даже для вектора результат — массив из одного элемента. Маски — основа «массивного» стиля Fortran, где логика выражается операциями над целыми массивами, а компилятор сам генерирует эффективный код.

Точность: проблема kind

Вещественные числа в компьютере приближённы, и от их разрядности зависит точность расчёта. Исторически Fortran предлагал real (одинарная точность, обычно 32 бита, ~7 десятичных цифр) и double precision (двойная, 64 бита, ~15-16 цифр). Но привязка к словам ненадёжна: «одинарная» на разных машинах могла отличаться. Современный Fortran ввёл параметр kind — целое число, точно идентифицирующее разновидность типа. real(kind=8) или просто real(8) у большинства компиляторов означает 64-битное вещественное. Но числовые значения kind не стандартизованы (у одного компилятора double — это kind 8, у другого мог бы быть иной), поэтому опираться на «магическое 8» не вполне переносимо.

Переносимый выбор точности

Правильный, переносимый способ — запросить kind по требованиям к точности и диапазону через встроенную функцию selected_real_kind(p, r): «дай мне тип хотя бы с p десятичными цифрами и диапазоном показателя r». Компилятор вернёт подходящий kind или −1, если такого нет. Этот kind определяют как именованную константу и используют повсюду.

module precision_mod
  implicit none
  ! двойная точность переносимо: >= 15 значащих цифр, диапазон 10^+-300
  integer, parameter :: dp = selected_real_kind(15, 307)
  integer, parameter :: sp = selected_real_kind(6, 37)   ! одинарная
end module precision_mod

program precision_demo
  use precision_mod
  implicit none
  real(dp) :: x
  x = 1.0_dp / 3.0_dp        ! литерал с суффиксом kind!
  print *, "1/3 в double:", x
  print *, "epsilon(dp):", epsilon(x)
  print *, "huge(dp):   ", huge(x)
  print *, "digits:     ", precision(x)
end program precision_demo

Вывод:

 1/3 в double:  0.33333333333333331     
 epsilon(dp):   2.2204460492503131E-016
 huge(dp):      1.7976931348623157E+308
 digits:              15

Критически важная деталь — суффикс kind у литералов: 1.0_dp, а не просто 1.0. Без суффикса литерал 1.0 имеет тип одинарной точности, и выражение 1.0 / 3.0 вычисляется в одинарной, теряя цифры ещё до присваивания в x. Запись 1.0_dp / 3.0_dp заставляет вычислять в двойной точности. Это типичнейшая скрытая ошибка точности: правильный тип переменной, но литералы без суффикса. В качественном численном коде каждая вещественная константа несёт суффикс kind.

Функции запроса свойств типа

Fortran даёт встроенные функции, сообщающие характеристики числового типа: epsilon(x) — машинное эпсилон (наименьшее ε, при котором 1+ε ≠ 1); huge(x) и tiny(x) — наибольшее и наименьшее положительное представимое; precision(x) — число надёжных десятичных цифр; radix, digits, maxexponent — детали внутреннего представления. Они незаменимы для написания численно устойчивого кода: например, сравнивать вещественные на «равенство» нужно не через ==, а через близость на уровне нескольких epsilon. Эти функции делают код самонастраивающимся под точность типа, а не зашивают магические пороги.

Как работает под капотом

Вещественные числа Fortran хранит в формате IEEE-754: одинарная точность — 1 бит знака, 8 бит экспоненты, 23 бита мантиссы (отсюда ~7 десятичных цифр); двойная — 1+11+52 бита (~15-16 цифр). Параметр kind — это компиляторный «тег», по которому выбирается формат и набор машинных инструкций. selected_real_kind(15, 307) компилятор вычисляет на этапе компиляции, перебирая доступные форматы, и возвращает kind двойной точности (на типичной платформе). Поскольку всё разрешается статически, накладных расходов в рантайме нет. Встроенные редукции вроде sum компилятор реализует развёрнутыми векторизованными циклами, а порядок суммирования может отличаться от наивного слева-направо — это влияет на последний бит результата при суммировании чисел разного масштаба, что нормально и обычно желательно (частичные суммы точнее). Понимание, что вещественная арифметика неточна и порядок операций влияет на младшие биты, — основа грамотного численного программирования.

Элементальные функции и whole-array стиль

Глубинная причина, по которой встроенные функции так удобны и быстры, — концепция элементальности. Элементальная функция определена для скаляра, но автоматически применяется к массиву любой формы поэлементно. Все математические встроенные функции — sqrt, sin, exp, abs и десятки других — элементальны. Поэтому sqrt(x) работает одинаково естественно и для одного числа, и для массива из миллиона элементов, и для трёхмерного поля. Это позволяет писать в whole-array стиле — оперируя целыми массивами как едиными величинами, без явных циклов:

real :: temp(100, 100), kelvin(100, 100)
kelvin = temp + 273.15           ! сдвиг всего поля одной строкой
real :: energy(100, 100)
energy = 0.5 * mass * velocity**2 ! формула применена к каждой точке

Такой код не только короче и ближе к математической записи, но и эффективнее: компилятор сам генерирует оптимальный, кэш-дружелюбный и векторизуемый цикл, зная семантику операции над всем массивом. Более того, вы можете объявлять собственные функции элементальными атрибутом elemental — тогда ваша скалярная функция автоматически распространится на массивы. Элементальная функция обязана быть pure (без побочных эффектов) и работать только со скалярными аргументами, но взамен получает «массивную» применимость даром. Это вершина выразительности численного Fortran: математическую модель пишут в терминах скаляра (одна точка, один элемент), а применяется она ко всему полю автоматически и быстро.

Накопление ошибки и численная устойчивость

Понимание неточности вещественной арифметики — не теоретическая тонкость, а практическая необходимость, потому что ошибки округления накапливаются и могут разрушить результат. Каждая операция над числами с плавающей точкой вносит крошечную относительную ошибку порядка машинного эпсилон. В одиночку она ничтожна, но в длинных вычислениях миллионы таких ошибок складываются, а иногда катастрофически усиливаются. Классический пример — потеря значимости (cancellation) при вычитании близких чисел: если a и b почти равны, то a - b теряет большинство верных цифр, потому что совпадающие старшие разряды сокращаются, обнажая накопленный «шум» в младших. Формула, математически верная, может давать численно бессмысленный ответ именно из-за такого вычитания.

! Неустойчиво: при больших x корни близки, вычитание теряет точность
! x1 = (-b + sqrt(b**2 - 4*a*c)) / (2*a)
! Устойчивый вариант для одного из корней использует
! свойство x1*x2 = c/a, чтобы избежать вычитания близких чисел.
real(8) :: a, b, c, disc, q, x1, x2
disc = sqrt(b**2 - 4.0d0*a*c)
q = -0.5d0 * (b + sign(disc, b))   ! sign выбирает сложение, не вычитание
x1 = q / a
x2 = c / q                         ! второй корень без потери значимости

Этот приём решения квадратного уравнения — хрестоматийный: наивная формула теряет точность для одного из корней, а переформулировка через sign и соотношение Виета её сохраняет. Подобные «численно устойчивые» переформулировки — целая дисциплина (численный анализ), и грамотный вычислитель их знает для типовых задач: суммирования (алгоритм Кэхэна компенсирует ошибку), решения систем (выбор главного элемента), вычисления дисперсии (одно-проходные формулы неустойчивы). Ключевые инструменты Fortran здесь — двойная точность по умолчанию для серьёзных расчётов (она отодвигает порог накопления ошибки), функции запроса свойств типа (epsilon для разумных порогов сравнения) и понимание, что порядок операций влияет на результат. Запомните главное: вещественное число в компьютере — это приближение, арифметика над ним неточна, и профессиональный численный код пишется с учётом этого, а не в наивном предположении, что 0.1 + 0.2 в точности равно 0.3 (оно не равно).

Частые ошибки

• Литералы без суффикса kind. x = 1.0 / 3.0 при real(dp) :: x считает деление в одинарной точности. Пишите 1.0_dp / 3.0_dp.
• Опора на «магический» kind 8. Значения kind не стандартизованы. Используйте selected_real_kind и именованную константу dp.
• Сравнение вещественных через ==. Из-за округления точное равенство почти никогда не выполняется. Сравнивайте по близости с порогом порядка epsilon.
• Ручные циклы вместо встроенных редукций. sum, maxval, dot_product обычно быстрее и понятнее ручного цикла — компилятор их отлично оптимизирует.
• Игнорировать переполнение. Промежуточные значения могут превысить huge; масштабируйте данные или используйте более широкий тип, опираясь на huge/tiny.

Итоги

• Встроенные функции работают над массивами целиком и оптимизированы компилятором — предпочитайте их ручным циклам.
• Маски (mask=) позволяют условные редукции без явных циклов: sum(a, mask = a > 0).
• Точность задаётся параметром kind; переносимо его получают через selected_real_kind(p, r).
• Литералы обязаны нести суффикс kind (1.0_dp), иначе вычисления идут в одинарной точности.
• Функции epsilon/huge/tiny/precision делают численный код самонастраивающимся и устойчивым.