Вычисляем число π

Сегодня я расскажу об одном интересном числе, а точнее о методах его вычисления. Кто не понял, сегодняшняя статья будет о вычислениях числа π (Пи), делать мы это будем с помощью Icon.

Итак, метод первый — используем сходящийся числовой ряд, найденный математиком С. Раманужаном (формулу можете найти в интернете сами):

procedure ramanujan(n)
local a,k,pi,b,c,s,d,e
a:=(2*sqrt(2))/9801.0
k:=0
s:=0
while k<n do {
	b:=factorial(4*k)
	c:=1103+26390*k
	d:=factorial(k)^4
	e:=396^(4*k)
	s+:=(b*c)/(d*e)
	k+:=1
}
return 1/(s*a)
end

procedure factorial_i(n)
local i,p
if n<2 then return 1 else {
 i:=2
 p:=1
 while i<n do {
 p*:=i
 i+:=1
 }
 return p
}
end

Из процедуры понятно, что сам ряд дает лишь единицу, деленную на π (но я исправил этот недостаток). Процедура требует лишь один аргумент — количество итераций (советую не ставить слишком большим), и приемлемый результат можно увидеть уже после 100 итераций (получаем число — 3.14159273).

Следующий интересный подход к вычислению числа π — это теория вероятности, скомбинированная с теорией чисел.

Из теории чисел известно, что вероятность того, что два наугад взятых числа взаимно простые (т.е. не имеют общих делителей, кроме 1) составляет 6 деленное на π квадрат. Вот этим я и воспользовался:

procedure pi_simple(n)
local x,y,i,m,p
m:=0.0
randomize()
every i:=1 to n do {
   x:=?400
   y:=?1000
   if gcd(x,y)=1 then m+:=1.0
}
p:=m/n
return sqrt(6/p)
end

Стоит сказать пару слов о том, как я до этого додумался… Да в принципе, никак. Нашел интересную статью (ссылка для любопытствующих — О числе Пи) и опробовал, но…

Мне не во всем понравился алгоритм, предложенный автором для этого случая — зачем предлагать сравнивать делители, если можно воспользоваться НОД (Наибольшим Общим Делителем)? Непонятно.

В общем, я усовершенствовал алгоритм, заменив сравнение делителей на сравнение НОД двух чисел (которые выбраны случайно) с 1, так как если НОД двух чисел равен 1, то они являются взаимно простыми.

Поначалу, я написал такую процедуру для вычисления НОД:

procedure gcd(m,n)
if m=0 then return n
if n=0 then return m
if m=1|n=1 then return 1 else {
  if m%2=0 & n%2=1 then return gcd(m/2,n)
  if m%2=1 & n%2=0 then return gcd(m,n/2)
  if m%2=0 & n%2=0 then return 2*gcd(m/2,n/2)
  if m%2=1 & n%2=1 then {
    if n>m then return gcd((n-m)/2,m) else return gcd((m-n)/2,n)
  }
}
end

Но, используя библиотеку IPL random я выпал в осадок, увидев в терминале ошибку «inconsistence redeclaration» (недопустимое переопределение) для процедуры gcd!

Оказывается, с подключением библиотеки random подключается библиотека factors, которая содержит процедуры factorial и gcd. Это позволило отказаться от написанных своими руками факториала и НОД (однако, я все-таки решил сохранить эти функции, но с небольшим изменением имен — уж очень показателен опыт).

Кстати, данный метод после 100 итераций дает результат — 2.828427125. Не густо, как в общем-то и в других вероятностных методах.

P.S: В Icon константа &pi соответствует числу Пи. Не верите? Испытайте команду write(&pi)!

И на всякий случай полный код вычислений:

link random
procedure main()
write(pi_simple(100))
write(ramanujan(100))
end

procedure ramanujan(n)
local a,k,pi,b,c,s,d,e
a:=(2*sqrt(2))/9801.0
k:=0
s:=0
while k<n do {
	b:=factorial(4*k)
	c:=1103+26390*k
	d:=factorial(k)^4
	e:=396^(4*k)
	s+:=(b*c)/(d*e)
	k+:=1
}
return 1/(s*a)
end

procedure factorial_i(n)
local i,p
if n<2 then return 1 else {
 i:=2
 p:=1
 while i<n do {
 p*:=i
 i+:=1
 }
 return p
}
end

procedure pi_simple(n)
local x,y,i,m,p
m:=0.0
randomize()
every i:=1 to n do {
 x:=?400
 y:=?1000
 if gcd(x,y)=1 then m+:=1.0
}
p:=m/n
return sqrt(6/p)
end

procedure gcd_b(m,n)
if m=0 then return n
if n=0 then return m
 	if m=1|n=1 then return 1 else {
 	if m%2=0 & n%2=1 then return gcd(m/2,n)
        if m%2=1 & n%2=0 then return gcd(m,n/2)
        if m%2=0 & n%2=0 then return 2*gcd(m/2,n/2)
        if m%2=1 & n%2=1 then {
        	if n>m then return gcd((n-m)/2,m) else return gcd((m-n)/2,n) 
        	
        }
 	}

end