Сегодня я расскажу об одном интересном числе, а точнее о методах его вычисления. Кто не понял, сегодняшняя статья будет о вычислениях числа π (Пи), делать мы это будем с помощью 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