В этом вопросе есть нечто большее, чем кажется на первый взгляд. Почему 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 n b n-1 .. .b 0 x 2 m . Если обратная тригонометрическая функция придумала наилучшее числовое двоичное приближение для своего аргумента, мы не хотим терять точность при умножении. Для этого нам нужно умножить только на степени 2.
π = 4 arctan(1) = 2 arcsin(1)
= 2 arccos(0)
= -4 arctan(-1) = arccos(-1) = -2 arcsin(-1)
примечание: это видно в представлении IEEE-754 binary64 (наиболее распространенная форма 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)
нечеткий аргумент реализации: с этого момента все немного зависит от реализации обратных тригонометрических функций. Иногда arccos и arcsin являются производными от atan2 и atan2 как
ACOS(x) = ATAN2(SQRT(1-x*x),1)
ASIN(x) = ATAN2(1,SQRT(1-x*x))
или более конкретно с числовой точки зрения:
ACOS(x) = ATAN2(SQRT((1+x)*(1-x)),1)
ASIN(x) = ATAN2(1,SQRT((1+x)*(1-x)))
Кроме того, atan2 является частью x86 Инструкции, установленной как FPATAN, в то время как другие - нет. С этой целью я бы поспорил использование:
π = 4 arctan(1)
над всеми остальными.
Примечание: это нечеткий аргумент. Я уверен, что есть люди с лучшим мнением по этому поводу.
Аргумент Фортрана: почему мы должны аппроксимировать π как:
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
Ответ лежит в стандарте Фортрана . Стандарт никогда гласит, что 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