LINUX.ORG.RU

Вычислить число пи с точностью до 1000-го знака


0

2

Предлагаю холивар
В интернете полно реализаций на разных языках
Я чуть позже выложу свой на питоне, в нем всего каких-то 10 строк, решение выводится за минуту на обычном компе, Архимед наверно перевернулся бы в гробу
В основу моего решения положен далеко не самый быстрый алгоритм, описанный еще Антифонтом: сторона квадрата, вписанного в окружность с единичным радиусом, равна корню из двух. Если удвоить число сторон до 8, то сторона будет равна более сложному выражению: нужно вычесть корень квадратный из разницы между двойкой и все тем же корнем из двух. И т.д. для 16,32, 64 .... - угольников

★★★★★

Последнее исправление: kto_tama (всего исправлений: 2)

Ответ на: комментарий от invy

pi - не целое. И вычисления там тоже не целочисленные.

и что? Просто умножить все на 10^1000 как самое тупое решение, хотя в питоне наверное что-то и так сделано.

abs ★★★
()
Ответ на: комментарий от post-factum

Только разница в том, что ты — современник Беллара, а не Ньютона.

посыпаю голову пеплом

kto_tama ★★★★★
() автор топика

давай чтоб под вендой еще работало!

dk-
()
Ответ на: комментарий от abs

подсказка:

from decimal import *
getcontext().prec = 1000
print "%s" % (Decimal(math.pi))
.

kto_tama ★★★★★
() автор топика
Ответ на: комментарий от devl547

Он француз. Так что таки Беллар.

anonymous
()
Ответ на: комментарий от LongLiveUbuntu

Беллар использовал в вычислениях формулу Чудновского, а не Рамануджана, хотя они подобны.

post-factum ★★★★★
()

В интернете полно реализаций на разных языках

https://gmplib.org/download/misc/gmp-chudnovsky.c

Ну и в чём смысл? Никто ведь не будет писать это заново. А если и будет, как проверить?

i-rinat ★★★★★
()
Ответ на: комментарий от Eddy_Em

А, фигушки: монте-карлой нужна будет длинная математика и жутко долго будет.

Тут нужен метод, позволяющий сразу же цифры клепать поочередно.

Eddy_Em ☆☆☆☆☆
()
Ответ на: комментарий от Eddy_Em

Да и обычный генератор ПСЧ — говно:

cat montepi.c
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <stdint.h>

#define MAX (uint64_t)40000000

int main(){
	uint64_t L, In = 0;
	double x,y;
	srand48(time(NULL));
	for(L = 0; L < MAX; L++){
		x = drand48();
		y = drand48();
		if(x*x + y*y < 1) In++;
	}
	printf("%llu\n", In);
}

gcc montepi.c -o montepi && ./montepi 
31419386
./montepi 
31423255
./montepi 
31419588

Eddy_Em ☆☆☆☆☆
()
Ответ на: комментарий от zolden

лехко:

from decimal import *

getcontext().prec = 2000

radical = 0
a=4
for i in range(1,2000):
  if i==1:
    radical = Decimal(2).sqrt()
  else :
    _radical = Decimal(2 + radical).sqrt()
    res = Decimal(2 - _radical).sqrt()
    print "i=%s  π=%s\n" % (i, res * a )  
    radical = _radical
  a = a * 2
kto_tama ★★★★★
() автор топика
Последнее исправление: kto_tama (всего исправлений: 1)
Ответ на: комментарий от Manhunt

А что, похоже:

octave:9> function x=pi_i(N);i=[0:N];x=sum(1./16.^i.*(4./(8*i+1)-2./(8*i+4)-1./(8*i+5)-1./(8*i+6)));endfunction
octave:10> for z=[0:10];printf("%1.20f\n", pi_i([z]));endfor
3.13333333333333330373
3.14142246642246636412
3.14158739034658163192
3.14159245756743565892
3.14159264546033645260
3.14159265322808778365
3.14159265357288086662
3.14159265358897288323
3.14159265358975225979
3.14159265358979133964
3.14159265358979311600

Если педивикия не врет, то первые 20 знаков после запятой: 3.14159265358979323846

К сожалению, последнее число в моем листинге — предел точности double. Дальше выжать из октавы не знаю как.

Eddy_Em ☆☆☆☆☆
()
Ответ на: комментарий от Eddy_Em

Обалдеть! Несчастных 11 итераций дают 15 чисел после запятой!!!

Вот только физический смысл этой формулы совершенно непонятен.

Eddy_Em ☆☆☆☆☆
()

другой алгориттм - по формуле Сриниваса Рамануджана:

k1 = (Decimal(2).sqrt() * 2) 
k2 = Decimal(k1/9801)
#print k2
f5 = Decimal(0) 

for i in range(0,1000):
   f1  = Decimal(math.factorial(4*i))
   f2  = Decimal(1103 + 26390 * i)
   ff1 = Decimal(f1 * f2)
   f3  = Decimal((math.factorial(i))**4)
   f4  = Decimal(396**(4*i))
   ff3 = Decimal(f3 * f4)
   ff4 = Decimal(ff1 / ff3)
   f5  = f5 + ff4 

f5 = f5 * k2
  
print 1/f5

и сюда же в кучу - алгоритм братьев Чудновских:

k1 = Decimal(10005).sqrt()
k2 = 426880 * k1
k3 = Decimal(1/(k2))

s = 0 

for i in range(0,100):
   f1  = Decimal(math.factorial(6*i))
   f2  = Decimal(13591409 + 545140134*i)
   f3  = f1 * f2
   f4  = Decimal(math.factorial(3*i))
   f5  = Decimal((math.factorial(i))**3)
   f6  = Decimal((-640320)**(3*i))
   f7  = f4 * f5 * f6
   f8  = f3 / f7
   s   = s + f8

s = s * k3

print Decimal(1 / s)

беллара мне уже лень писать

kto_tama ★★★★★
() автор топика
Последнее исправление: kto_tama (всего исправлений: 3)
#include <stdio.h>
#include <mpfr.h>
int main(int argc, const char** argv)
{
    mpfr_t pi;
    mpfr_init2(pi, 1024*1024*8);
    mpfr_const_pi(pi, MPFR_RNDD);
    mpfr_out_str (stdout, 10, 0, pi, MPFR_RNDD);
    mpfr_clear(pi);
    mpfr_free_cache();
    return 0;
}
$ clang -I/opt/mpfr/include -I/opt/gmp/include -L/opt/mpfr/lib -L/opt/gmp/lib -Wl,-no_pie -march=core-avx2 -O3 -o main main.c -lgmp -lmpfr
$
$
$ time sh -c './main > pi.txt'

real	0m3.768s
user	0m3.722s
sys	0m0.041s
$
$ stat -c "%s" pi.txt
2525225
slyjoeh ★★★
()
Ответ на: комментарий от kto_tama

Проверка в цикле лишняя. Без неё будет красивее смотреться.

Reset ★★★★★
()
Ответ на: комментарий от zolden

А есть какой-нибудь алгоритм, не требующий использования библиотеки с длинными числами?

Как генерируют тысячные и больше знаки числа пи? Это ж оперативы не хватит!

Eddy_Em ☆☆☆☆☆
()
Ответ на: комментарий от invy

pi - не целое. И вычисления там тоже не целочисленные.

И что? Почему результат не хранить в целых?

ieeya
()

Предлагаю холивар

В военное время по приказу командования π может быть целочисленным 3 или 4 в зависимости от обстановки :)

quickquest ★★★★★
()
Ответ на: комментарий от Eddy_Em

Да и обычный генератор ПСЧ — говно

идиты к Кнуту. У него есть хорошие ГПСЧ.

emulek
()
Ответ на: комментарий от fijiol

а я уж надеялся тут кто-нибудь тейлора начнет советовать.. эх

_что_ ты собрался раскладывать в ряд Тейлора?

emulek
()
Ответ на: комментарий от Eddy_Em

К сожалению, последнее число в моем листинге — предел точности double. Дальше выжать из октавы не знаю как.

man 1 bc

emulek
()
Вы не можете добавлять комментарии в эту тему. Тема перемещена в архив.