современный фортран , Бортеньев. О. В. Бартеньев Современный Фортран
 Скачать 2.24 Mb.
|
|
6.8. Математические элементные функции Фортран содержит математические функции вычисления корня, логарифмов, экспоненты и тригонометрических функций. Тип и параметр разновидности типа результата такие же, как у первого аргумента. В разделе приведены формы вызова функций с родовыми именами. Специфические имена функций даны в разд. 8.18.2. Когда параметрами логарифмических и тригонометрических функций являются комплексные числа, то функции возвращают комплексное число, аргумент ϑ которого равен главному значению аргумента комплексного числа в радианах (- π < ϑ ≤ π). 6.8.1. Экспоненциальная, логарифмическая функции и квадратный корень EXP(x) - возвращает e x = e**x для вещественного или комплексного x. В случае комплексного x = (a, b) результат равен e**a*(cos(b) + i sin(b)). LOG(x) - возвращает значение натурального логарифма для вещественного или комплексного x. В случае вещественного аргумента значение x должно быть больше нуля. В случае комплексного аргумента x не должен быть нулем. Если х комплексного типа, то действительный компонент результата равен натуральному логарифму модуля x, мнимый компонент - главному значению аргумента x в радианах, т. е. если x = (a, b), то LOG(x) = = (LOG(SQRT(a**2 + b**2)), ATAN2(a, b)). LOG10(x) - возвращает десятичный логарифм вещественного аргумента. Значение х должно быть больше нуля. 179 О. В. Бартеньев. Современный ФОРТРАН SQRT(x) - возвращает квадратный корень для вещественного или комплексного аргумента x. В случае вещественного аргумента значение x должно быть больше нуля. В случае комплексного х функция возвращает число, модуль которого равен корню квадратному из модуля x и угол которого равен половине угла x. Так, если x = (a, b), то SQRT(x) = SQRT(a**2 + b**2) * * e -j * 0.5 * atan(a / b) Извлечь корень можно также, применив операцию возведения в степень: SQRT(x) = x**0.5. Однако применять следует функцию SQRT, поскольку SQRT(x) выполняется быстрее, чем x**0.5. 6.8.2. Тригонометрические функции В Фортране существуют как тригонометрические функции, в которых аргумент должен быть задан в радианах, например SIN(x), так и функции, аргумент которых задается в градусах, например SIND(x). Синус и арксинус SIN(x) - возвращает синус вещественного или комплексного аргумента x, который интерпретируется как значение в радианах. SIND(x) - возвращает синус вещественного или комплексного аргумента x, который интерпретируется как значение в градусах. ASIN(x) - возвращает арксинус вещественного аргумента x (| x | ≤ 1), выраженный в радианах в интервале - π/2 ≤ ASIN(x) ≤ π/2. ASIND(x) - возвращает арксинус вещественного аргумента x (| x | ≤ 1), выраженный в градусах в интервале -90 ≤ ASIN(x) ≤ 90. Косинус и арккосинус COS(x) - возвращает косинус вещественного или комплексного аргумента x, который интерпретируется как значение в радианах. COSD(x) - возвращает косинус вещественного или комплексного аргумента x, который интерпретируется как значение в градусах. ACOS(x) - возвращает арккосинус вещественного аргумента x (| x | ≤ 1), выраженный в радианах в интервале 0 ≤ ACOS(x) ≤ π. ACOSD(x) - возвращает арккосинус вещественного аргумента x (| x | ≤ 1), выраженный в градусах в интервале 0 ≤ ACOS(x) ≤ 180. Тангенс, котангенс и арктангенс TAN(x) - возвращает тангенс вещественного аргумента x, который интерпретируется как значение в радианах. TAND(x) - возвращает тангенс вещественного аргумента x, который интерпретируется как значение в градусах. COTAN(x) - возвращает котангенс вещественного аргумента x (x ≠0), который интерпретируется как значение в радианах. 180 6. Встроенные процедуры ATAN(x) - возвращает арктангенс вещественного аргумента x, выраженный в радианах в интервале - π/2 < ATAN(x) < π/2. ATAND(x) - возвращает арктангенс вещественного аргумента x, выраженный в градусах в интервале -90 < ATAN(x) < 90. ATAN2(y, x) - возвращает арктангенс (y/x), выраженный в радианах в интервале - π ≤ ATAN2(y, x) ≤ π. Аргументы y и x должны быть вещественного типа с одинаковым значением параметра разновидности и не могут одновременно равняться нулю. ATAN2D(y, x) - возвращает арктангенс (y/x), выраженный в градусах в интервале -180 ≤ ATAN2D(y, x) ≤ 180. Аргументы y и x должны быть вещественного типа с одинаковым значением параметра разновидности и не могут одновременно равняться нулю. Диапазон результатов функций ATAN2 и ATAN2D: Аргументы Результат (в рад ) Результат (в град. ) y > 0 Результат > 0 Результат > 0 y = 0 и x > 0 Результат = 0 Результат = 0 y = 0 и x < 0 Результат = π Результат = 180 y < 0 Результат < 0 Результат < 0 x = 0 и y > 0 Результат = π/2 Результат = 90 x = 0 и y < 0 Результат = - π/2 Результат = -90 Гиперболические тригонометрические функции SINH(x) - гиперболический синус для выраженного в радианах вещественного аргумента x. COSH(x) - гиперболический косинус для выраженного в радианах вещественного аргумента x. TANH(x) - гиперболический тангенс для выраженного в радианах вещественного аргумента x. 6.9. Функции для массивов Встроенные функции для работы с массивами позволяют выполнять вычисления в массивах, получать справочные данные о массиве и преобразовывать массивы. Встроенные функции обработки массивов рассмотрены в разд. 4.12. Помимо встроенных CVF и FPS содержат дополнительные подпрограмму SORTQQ сортировки одномерного массива и целочисленную функцию BSEARCHQQ бинарного поиска в отсортированном массиве. Для их вызова необходимо сослаться на модуль MSFLIB. CALL SORTQQ(adrarray, count, size) - сортирует одномерный массив, адрес которого равен adrarray. Для вычисления адреса применяется 181 О. В. Бартеньев. Современный ФОРТРАН функция LOC. Сортируемый массив не должен быть производного типа. Параметр count имеет вид связи INOUT и при вызове равен числу элементов массива, подлежащих сортировке, а на выходе - числу реально отсортированных элементов. Тип параметров adrarray, count - стандартный целый. Параметр size является положительной константой стандартного целого типа (size < 32,767), задающей тип и разновидность типа сортируемого массива. В файле msflib.f90 определены следующие константы: Константа Типы массива SRT$INTEGER1 INTEGER(1) SRT$INTEGER2 INTEGER(2) или эквивалентный SRT$INTEGER4 INTEGER(4) или эквивалентный SRT$REAL4 REAL(4) или эквивалентный SRT$REAL8 REAL(8) или эквивалентный Если величина size не является именованной константой приведенной таблицы и меньше чем 32'767, то предполагается, что задан символьный массив, длина элемента которого равна size. Чтобы убедиться в том, что сортировка выполнена успешно, следует сравнить значения параметра count до и после сортировки. При положительном результате они совпадают. Предупреждение. Адрес сортируемого массива должен быть вычислен функцией LOC. Значения параметров count и size должны точно описывать характеристики массива. Если же подпрограмма SORTQQ получила неверные параметры, то будет выполнена попытка сортировки некоторой области памяти. Если память принадлежит текущему процессу, то сортировка будет выполнена, иначе операционная система выполнит функции защиты памяти и остановит вычисления. BSEARCHQQ(adrkey, adrarray, length, size) - выполняет бинарный поиск значения, которое содержится в переменной, расположенной по адресу adrkey. Поиск выполняется в отсортированном одномерном массиве, первый элемент которого имеет адрес adrarray. Функция возвращает индекс искомого элемента или 0, если элемент не найден. Тип результата и параметров adrkey, adrarray, length и size - стандартный целый. Элементы массива не могут быть производного типа. Параметр length равен числу элементов массива. Смысл параметра size пояснен при рассмотрении подпрограммы SORTQQ. До выполнения поиска массив должен быть отсортирован по возрастанию значений его элементов. 182 6. Встроенные процедуры Предупреждение. Адреса сортируемого массива и элемента должны быть вычислены функцией LOC. Значения параметров count и size должны точно описывать характеристики массива. К тому же искомый элемент должен иметь тот же тип и разновидность типа, что и массив, в котором выполняется поиск. Эти характеристики задаются параметром size. Если же функция BSEARCHQQ получила неверные параметры, то, если память принадлежит текущему процессу, будет выполнена попытка поиска в некоторой области памяти, иначе операционная система выполнит функции защиты памяти и остановит вычисления. Пример. Найти в массиве все равные заданному значению элементы. Для решения предварительно отсортируем массив, а далее выполним поиск, учитывая, что равные элементы отсортированного массива следуют подряд. use msflib integer(4) a(20000), n, n2, ada, i, k integer(4) :: ke = 234 ! Вычислим, сколько раз в массиве a real(4) rv ! встречается число ke n = size(a); n2 = n ! Запомним размер массива a do i = 1, n ! Заполним массив a случайным образом call random(rv) a(i) = int(rv * 1000.0) end do ada = loc(a) ! Адрес массива a call sortqq(ada, n, SRT$INTEGER4) ! Сортировка массива a if(n .ne. n2 ) stop 'Sorting error' k = bsearchqq(loc(ke), ada, n, SRT$INTEGER4) if(k == 0) then print *, 'Элемент ke = ', ke, ' не найден' stop end if i = k ! Поиск всех равных ke элементов do while(a(i) == ke .and. i <= n) print *, 'Элемент с индексом i = ', i, ' равен ke; ke = ', ke i = i + 1 end do print *, 'Всего найдено элементов: ', i - k ! 21 end На самом деле результат неверен. Действительно, вычислим nke - число равных ke элементов в цикле nke = 0 do i = 1, size(a) if(a(i) == ke) nke = nke + 1 end do print *, nke ! 28 183 О. В. Бартеньев. Современный ФОРТРАН Оказывается, что nke = 28, а не 21, как это было найдено выше. Дело в том, что BSEARCHQQ в общем случае находит в упорядоченном векторе не первый равный заданному числу элемент, а один из последующих. Чтобы вернуться к первому искомому элементу, в программу после вызова BSEARCHQQ нужно добавить код if(k > 1) then ! Если в векторе есть элемент, равный ke, it = k ! выполним перемещение назад, корректируя do k = it, 1, -1 ! величину k - позицию первого равного ke элемента if(a(k) /= ke) exit end do k = k + 1 end if 6.10. Справочные функции для любых типов ALLOCATED(array) - возвращает .TRUE., если память выделена под массив array, .FALSE. - если не выделена. Параметр array должен иметь атрибут ALLOCATABLE. Результат будет неопределенным, если не определен статус массива array. Результат имеет стандартный логический тип. Пример: real, allocatable :: a(:) if(.not. allocated(a)) allocate(a(10)) ASSOCIATED(pointer [, target]). Параметр pointer должен быть ссылкой. Состояние привязки pointer не должно быть неопределенным. Если параметр target опущен, то функция ASSOCIATED возвращает .TRUE., если ссылка pointer прикреплена к какому-либо адресату. Если параметр target задан и имеет атрибут TARGET, то результат равен .TRUE., если ссылка pointer прикреплена к target. Если target задан и имеет атрибут POINTER, то результат равен .TRUE., если как pointer, так и target прикреплены к одному адресату. При этом состояние привязки ссылки tar- get не должно быть неопределенным. Во всех других случаях результат равен .FALSE. Результат имеет стандартный логический тип. В случае массивов .TRUE. возвращается, если совпадают формы аргументов и если элементы ссылки в порядке их следования прикреплены к соответствующим элементам адресата. Пример: real, pointer :: a(:), b(:), c(:) real, target :: e(10) a => e ! Назначение ссылки b => e c => e(1:10:2) ! Ссылка прикрепляется к сечению массива print *, associated(a, e) ! T 184 6. Встроенные процедуры print *, associated(a, b) ! T print *, associated(c) ! T print *, associated(c, e) ! F print *, associated(c, e(1:10:2)) ! T PRESENT(a) - определяет, задан ли необязательный формальный параметр a при вызове процедуры. Функция PRESENT может быть вызвана только в процедуре с необязательными параметрами. Функция возвращает .TRUE., если в вызове процедуры присутствует фактический параметр, соответствующий формальному параметру a. В противном случае PRE- SENT возвращает .FALSE.. Результат имеет стандартный логический тип. Пример: call who(1, 2) ! Напечатает: a present ! b present call who(1) ! Напечатает: a present call who(b = 2) ! Напечатает: b present call who( ) ! Напечатает: No one contains subroutine who(a, b) integer(4), optional :: a, b if(present(a)) print *, 'a present' if(present(b)) print *, 'b present' if(.not. present(a) .and. .not. present(b)) print *, 'No one' end subroutine who end KIND(x) - возвращает стандартное целое, равное значению параметра разновидности аргумента x. Пример: real(kind(1e0)), parameter :: one = 1.0 ! Вещественные константа one real(kind(1e0)) :: err ! и переменная err одинарной точности 6.11. Числовые справочные и преобразовывающие функции 6.11.1. Модели данных целого и вещественного типа Каждая разновидность целого и вещественного типа содержит конечное множество чисел. Так, тип INTEGER(2) представляет все целые числа из диапазона от -32,768 до +32,767. Каждое такое множество чисел может быть описано моделью. Данные о параметрах модели заданной разновидности типа и о конкретных характеристиках числа в задающей его модели позволяют получать встроенные числовые справочные и преобразовывающие функции, а также функции IMACH, AMACH и DMACH библиотеки IMSL. Двоичное представление целого числа i задается формулой 185 О. В. Бартеньев. Современный ФОРТРАН 1 1 0 ) 1 ( − − = Q знак s b b b i … E мантисса P знак s b b b b x 2 ) 1 ( 1 2 1 0 × − = − … , где s - это 0 или 1 (+ или -); b i - двоичное число (0 или 1); Q - число цифр в целом числе по основанию 2. Вещественные числа с плавающей точкой представляются в CVF и FPS в близком соответствии со стандартом IEEE для арифметики с плавающей точкой (ANSI/IEEE Std 754-1985, IEEE Standard for Binary Floating-Point Arithmetic, 1985). Фортран поддерживает форматы одинарной точности - тип REAL(4), двойной точности - тип REAL(8) и повышенной точности, используемой для выполнения промежуточных операций. Например, в следующем коде: real(4) :: a, b, c, d, f ! Задана опция /Op компилятора FPS b = 0.0; c = 1.0e30; d = 1.0e30; f = 1.0e-30 a = (b + (c * d)) / 2.0 * f print *, a ! 5.000000E+29 промежуточные вычисления, если задана опция компилятора /Op, FPS выполнит с повышенной точностью. Если же при компиляции задана опция /Oxp, предусматривающая полную оптимизацию скорости вычислений и проверку ошибок, то на этапе компиляции возникнут предупреждения вида (по причине переполнения в результате умножения c * d = 1.0e30 * 1.0e30) warning F4756: overflow in constant arithmetic а результатом вычислений будет машинная бесконечность - 1#INF........ Двоичное представление вещественного числа x с плавающей точкой задается формулой , где s - это 0 или 1 (+ или -); b i - двоичное число (0 или 1); P - число цифр в мантиссе нормализованного представления вещественного числа по основанию 2; E - целое число, называемое (двоичным) порядком, из отрезка E min ≤ E ≤ E max . В табл. 6.1 приводятся значения параметров модели вещественных чисел для одинарной, двойной и повышенной точности. Таблица 6.1. Параметры модели чисел стандарта IEEE Std 754 Параметр Точность одинарная двойная повышенная Число бит для знака 1 1 1 P 24 53 64 E max +128 +1024 +16384 186 6. Встроенные процедуры E min -125 -1021 -16381 Смещение двоичного порядка +126 +1022 +16382 Число бит для двоичного порядка 8 11 15 Число бит для числа 32 64 80 Стандарт требует, чтобы числа одинарной и двойной точности представлялись в нормализованном виде, поэтому всегда b 0 = 1 и, следовательно, для представления мантиссы чисел одинарной точности требуется 23 бита, а двойной - 52. Числа повышенной точности не нормализуются, поэтому для мантиссы требуется 64 бита. При записи порядка к нему с целью повышения скорости вычислений добавляется смещение, поэтому порядок всегда представляется в виде положительного числа e. Реально, однако, порядок E = e - смещение. Структура формата вещественных чисел двойной точности с плавающей точкой приведена на рис. 6.1. 1 11 52 Знак Порядок + смещение Мантисса Рис. 6.1. Структура IEEE-формата вещественных чисел двойной точности |