Зачем определять PI = 4*ATAN(1.d0)

Этот стиль гарантирует, что при присвоении значения PI используется максимальная точность, доступная для ЛЮБОЙ архитектуры.

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

Они эквивалентны, и иногда вы их тоже будете видеть:

PI=DACOS(-1.D0)
PI=2.D0*DASIN(1.D0)

Я полагаю, это потому, что это самая короткая серия на число Пи. Это также означает, что он САМЫЙ ТОЧНЫЙ.

Ряд Грегори-Лейбница (4/1 - 4/3 + 4/5 - 4/7...) равен числу пи.

atan(x) = x^1/1 - x^3/3 + x^5/5 - x^7/7...

Итак, атан(1) = 1/1 - 1/3 + 1/5 - 1/7 + 1/9... 4 * атан(1) = 4/1 - 4/3 + 4/5 - 4/ 7 + 4/9...

Это равно ряду Грегори-Лейбница и, следовательно, равно пи, приблизительно 3,1415926535 8979323846 2643383279 5028841971 69399373510.

Другой способ использовать атан и найти число пи:

pi = 16*atan(1/5) - 4*atan(1/239), но я думаю, что это сложнее.

Надеюсь, это поможет!

(Честно говоря, я думаю, что ряд Грегори-Лейбница был основан на атане, а не 4*атан(1) на основе ряда Грегори-Лейбница. Другими словами, НАСТОЯЩЕЕ доказательство таково:

sin^2 x + cos^2 x = 1 [Теорема] Если x = pi/4 радиан, то sin^2 x = cos^2 x или sin^2 x = cos^2 x = 1/2.

Тогда sin x = cos x = 1/(корень 2). тангенс х (sin х / cos х) = 1, атан х (1 / тангенс х) = 1.

Итак, если atan(x) = 1, x = pi/4 и atan(1) = pi/4. Наконец, 4*atan(1) = pi.)

Пожалуйста, не загружайте меня комментариями — я еще несовершеннолетний.

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

Напротив, указание pi в качестве константы обеспечивает точно такую ​​же точность, как и изначально заданная, что может не подходить для высоконаучных или математических приложений (как часто используется Fortran).

В этом вопросе есть нечто большее, чем кажется на первый взгляд. Почему 4 arctan(1)? Почему бы не любое другое представление, такое как 3 arccos(1/2)?

Это попытается найти ответ путем исключения.

Введение в математику. При использовании обратных тригонометрических функций, таких как arccos, arcsin и arctan, можно легко вычислить число π в Различные пути:

π = 4 arctan(1) = arccos(-1) = 2 arcsin(1) = 3 arccos(1/2) = 6 arcsin(1/2)
  = 3 arcsin(sqrt(3)/2) = 4 arcsin(sqrt(2)/2) = ...

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

Аргумент с плавающей запятой 1: хорошо известно, что конечное двоичное представление с плавающей запятой не может представлять все действительные числа. Некоторые примеры таких чисел: 1/3, 0.97, π, sqrt(2), .... С этой целью мы должны исключить любое математическое вычисление π, где аргумент обратных тригонометрических функций не может быть представлен численно. Это оставляет нам аргументы -1,-1/2,0,1/2 и 1.

π = 4 arctan(1) = 2 arcsin(1)
   = 3 arccos(1/2) = 6 arcsin(1/2)
   = 2 arccos(0)
   = 3/2 arccos(-1/2) = -6 arcsin(-1/2)
   = -4 arctan(-1) = arccos(-1) = -2 arcsin(-1)

Аргумент с плавающей точкой 2: в двоичном представлении число представляется как 0.b_nb_n-1.. .b₀ x 2^м. Если обратная тригонометрическая функция дает наилучшее числовое бинарное приближение для своего аргумента, мы не хотим терять точность при умножении. Для этого мы должны умножать только на степени 2.

π = 4 arctan(1) = 2 arcsin(1)
  = 2 arccos(0)
  = -4 arctan(-1) = arccos(-1) = -2 arcsin(-1)

Примечание. это видно в двоичном файле IEEE-75464 представление (наиболее распространенная форма DOUBLE PRECISION или kind=REAL64). Там у нас есть

write(*,'(F26.20)') 4.0d0*atan(1.0d0) -> "    3.14159265358979311600"
write(*,'(F26.20)') 3.0d0*acos(0.5d0) -> "    3.14159265358979356009"

Этой разницы нет в двоичном формате IEEE-75432 (наиболее распространенная форма REAL или kind=REAL32) и двоичный файл IEEE-754128 (наиболее распространенная форма kind=REAL128 )

Аргумент реализации: на ЦП Intel atan2 является частью инструкции x86. устанавливается как FPATAN, а другая обратная тригонометрическая функция получается из atan2. Потенциальным выводом может быть:

          mathematically         numerically
ACOS(x) = ATAN2(SQRT(1-x*x),1) = ATAN2(SQRT((1+x)*(1-x)),1)
ASIN(x) = ATAN2(1,SQRT(1-x*x)) = ATAN2(1,SQRT((1+x)*(1-x)))

Это можно увидеть в ассемблере этих инструкций (см. здесь). С этой целью я бы поспорил с использованием:

π = 4 arctan(1)

Примечание. это нечеткий аргумент. Я уверен, что есть люди с лучшим мнением по этому поводу.
Интересно почитать о FPATAN: Как реализован ли arctan?, тригонометрические инструкции x87

Аргумент Fortran: почему мы должны аппроксимировать π следующим образом:

integer, parameter :: sp = selected_real_kind(6, 37)
integer, parameter :: dp = selected_real_kind(15, 307)
integer, parameter :: qp = selected_real_kind(33, 4931)

real(kind=sp), parameter :: pi_sp = 4.0_sp*atan2(1.0_sp,1.0_sp)
real(kind=dp), parameter :: pi_dp = 4.0_dp*atan2(1.0_dp,1.0_dp)
real(kind=qp), parameter :: pi_qp = 4.0_qp*atan2(1.0_qp,1.0_qp)

и не :

real(kind=sp), parameter :: pi_sp = 3.14159265358979323846264338327950288_sp
real(kind=dp), parameter :: pi_dp = 3.14159265358979323846264338327950288_dp
real(kind=qp), parameter :: pi_qp = 3.14159265358979323846264338327950288_qp

Ответ содержится в стандарте Fortran. В стандарте никогда указано, что REAL любого вида должно представлять число с плавающей запятой IEEE-754. Представление REAL зависит от процессора. Это означает, что я могу запросить selected_real_kind(33, 4931) и ожидать получения двоичного128 числа с плавающей запятой, но я мог бы получить kind, который представляет число с плавающей запятой с гораздо большей точностью. Может быть, 100 цифр, кто знает. В этом случае моя приведенная выше строка чисел слишком короткая! Нельзя использовать это просто для уверенности? Даже этот файл может быть слишком коротким!

Интересный факт: sin(pi) is never zero

write(*,'(F17.11)') sin(pi_sp) => "   -0.00000008742"
write(*,'(F26.20)') sin(pi_dp) => "    0.00000000000000012246"
write(*,'(F44.38)') sin(pi_qp) => "    0.00000000000000000000000000000000008672"

что понимается как:

pi = 4 ATAN2(1,1) = π + δ
SIN(pi) = SIN(pi - π) = SIN(δ) ≈ δ

program print_pi
! use iso_fortran_env, sp=>real32, dp=>real64, qp=>real128

  integer, parameter :: sp = selected_real_kind(6, 37)
  integer, parameter :: dp = selected_real_kind(15, 307)
  integer, parameter :: qp = selected_real_kind(33, 4931)

  real(kind=sp), parameter :: pi_sp = 3.14159265358979323846264338327950288_sp
  real(kind=dp), parameter :: pi_dp = 3.14159265358979323846264338327950288_dp
  real(kind=qp), parameter :: pi_qp = 3.14159265358979323846264338327950288_qp
  
  write(*,'("SP "A17)') "3.14159265358..."
  write(*,'(F17.11)') pi_sp
  write(*,'(F17.11)')        acos(-1.0_sp)
  write(*,'(F17.11)') 2.0_sp*asin( 1.0_sp)
  write(*,'(F17.11)') 4.0_sp*atan2(1.0_sp,1.0_sp)
  write(*,'(F17.11)') 3.0_sp*acos(0.5_sp)
  write(*,'(F17.11)') 6.0_sp*asin(0.5_sp)

  write(*,'("DP "A26)') "3.14159265358979323846..."
  write(*,'(F26.20)') pi_dp
  write(*,'(F26.20)')        acos(-1.0_dp)
  write(*,'(F26.20)') 2.0_dp*asin( 1.0_dp)
  write(*,'(F26.20)') 4.0_dp*atan2(1.0_dp,1.0_dp)
  write(*,'(F26.20)') 3.0_dp*acos(0.5_dp)
  write(*,'(F26.20)') 6.0_dp*asin(0.5_dp)

  write(*,'("QP "A44)') "3.14159265358979323846264338327950288419..."
  write(*,'(F44.38)') pi_qp
  write(*,'(F44.38)')        acos(-1.0_qp)
  write(*,'(F44.38)') 2.0_qp*asin( 1.0_qp)
  write(*,'(F44.38)') 4.0_qp*atan2(1.0_qp,1.0_qp)
  write(*,'(F44.38)') 3.0_qp*acos(0.5_qp)
  write(*,'(F44.38)') 6.0_qp*asin(0.5_qp)

  write(*,'(F17.11)') sin(pi_sp)
  write(*,'(F26.20)') sin(pi_dp)
  write(*,'(F44.38)') sin(pi_qp)

end program print_pi

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