Calcular pi até 5 casas decimais

15

Isso vem de http://programmers.blogoverflow.com/2012/08/20-controversial-programming-opinions/

"Dado que o Pi pode ser estimado usando a função 4 * (1 - 1/3 + 1/5 - 1/7 + ...) com mais termos dando maior precisão, escreva uma função que calcule o Pi com uma precisão de 5 casas decimais. "

  • Observe que a estimativa deve ser feita calculando a sequência fornecida acima.
n
fonte
8
Você provavelmente deve adicionar mais algumas regras, caso contrário você terá respostas como (Python)p=lambda:3.14159
Matt
1
Você já viu codegolf.stackexchange.com/questions/506/… , que é muito semelhante? No mínimo, funções trigonométricas devem ser banidas para esse problema, pois permitem soluções triviais como este programa QBASIC: INT (4E5 * ATN (1)) /
1E5
Eu acho que você deve exigir que o algoritmo seja de aproximação sucessiva: quanto mais você computa, mais perto fica do pi.
21412
@DavidCarraher, embora isso seja matematicamente inevitável ao usar esta série, do ponto de vista analítico numérico é altamente dúbio. Uma série alternada de convergência lenta é um filho da propaganda por perda de significado.
Peter Taylor
2
Dupe, mas é tão antigo que não está aqui: stackoverflow.com/q/407518/12274
JB

Respostas:

10

JavaScript, 46 58 56 45 bytes

Atualização do ES6 : Acontece que há mais recursos disponíveis agora que cinco anos se passaram.

let f=(i=0,a=0)=>i>1e6?a:f(i+4,a+8/-~i/(i+3))

Esta versão ( 45 bytes; sim, leté necessária) funciona no modo estrito do ES6 em teoria . Na prática, você pode executá-lo na V8 (por exemplo, com nó) com --use-strict --harmony-tailcalls; o recurso Adequado Tailcalls ainda não foi amplamente implementado, infelizmente. No entanto, é um comportamento especificado, portanto tudo deve estar bem.

Se quisermos manter o que é amplamente implementado, e não exigir o modo estrito, podemos simplesmente usar a sintaxe ES6 da seta de fatura para funções, mas manter a mesma implementação de antes (sugerida por Brian H) por um custo de 48 bytes.

a=>{for(a=i=0;i<1e6;a+=8/++i/~-(i+=3));return a}

A escolha do nome para o parâmetro único não importa realmente , mas podemos escolher um dos nomes que usamos para minimizar a poluição no escopo global.


function(){for(a=i=0;i<1e6;a+=8/++i/~-(i+=3));return a}

Esta versão é uma expressão de função; adicione dois caracteres (por exemplo, " f") se desejar que ele seja nomeado. Esta versão derruba os globais ae i; isso pode ser evitado se adicionarmos " a,i" à lista de parâmetros.

Utiliza uma versão reformulada do algoritmo para contornar a necessidade de subtração.

 1/1 - 1/3  +   1/5 - 1/7   +    1/9 - 1/11  + ...
(3/3 - 1/3) + (7/35 - 5/35) + (11/99 - 9/99) + ...
    2/3     +      2/35     +       2/99     + ...
  2/(1*3)   +    2/(5*7)    +     2/(9*11)   + ...

Aqui está uma versão "simples" sem esse ajuste:

function(){for(a=0,i=1;i<1e6;i+=2)a+=[,4,,-4][i%4]/i;return a}

com 64 62 caracteres.

Obrigado a @ardnew pela sugestão de se livrar do 4*antes do return.


História

function(){for(a=i=0;i<1e6;a+=8/++i/~-(i+=3));return a}     // got rid of `i+=4`; restructured
// Old versions below.
function(){for(a=0,i=1;i<1e6;i+=4)a+=8/i/-~-~i;return a}    // got rid of `4*`
function(){for(a=0,i=1;i<1e6;i+=4)a+=2/i/-~-~i;return 4*a}
FireFly
fonte
oO trabalho muito bom, fatorando a subtração.
Acolito
1
excelente trabalho, mas precisa ser escrito como uma função adequada
ardnew
@ardnew: Obrigado, devo ter perdido esse detalhe quando li a descrição do problema. Eu atualizei e agora é uma expressão de função que pode ser chamada (lambda); não tenho certeza se isso é permitido ou se deve receber um nome. Se for esse o caso, são apenas mais dois caracteres.
FireFly
1
@FireFly, você também pode cortar 2 caracteres, alterando a+=2/i/-~-~i;return 4*aparaa+=8/i/-~-~i;return a
ardnew
@ardnew: oh, incrível; não pensei nisso. : D
FireFly
8

Python 59 bytes

print reduce(lambda x,p:p/2*x/p+2*10**999,range(6637,1,-2))

Isso imprime 1000 dígitos; um pouco mais do que o necessário 5. Em vez de usar a iteração prescrita, ele usa o seguinte:

pi = 2 + 1/3*(2 + 2/5*(2 + 3/7*(2 + 4/9*(2 + 5/11*(2 + ...)))))

O 6637(o denominador mais interno) pode ser formulado como:

dígitos * 2 * log 2 (10)

Isso implica uma convergência linear. Cada iteração mais profunda produzirá mais um bit binário de pi .

Se , no entanto, você insistir em usar aidentidade tan -1 , uma convergência semelhante poderá ser alcançada, se você não se importar em abordar o problema de maneira ligeiramente diferente. Analisando as somas parciais:

4.0, 2.66667, 3.46667, 2.89524, 3.33968, 2.97605, 3.28374, ...

é evidente que cada termo pula para frente e para trás em ambos os lados do ponto de convergência; a série tem convergência alternada. Além disso, cada termo está mais próximo do ponto de convergência do que o termo anterior; é absolutamente monotônico em relação ao seu ponto de convergência. A combinação dessas duas propriedades implica que a média aritmética de dois termos vizinhos está mais próxima do ponto de convergência do que qualquer um dos termos em si. Para ter uma idéia melhor do que quero dizer, considere a seguinte imagem:

Somas Parciais

A série externa é a original e a série interna é encontrada calculando a média de cada um dos termos vizinhos. Uma diferença notável. Mas o que é realmente notável é que essa nova série também possui convergência alternada e é absolutamente monotônica em relação ao seu ponto de convergência. Isso significa que esse processo pode ser aplicado repetidamente, ad nauseum.

Está bem. Mas como?

Algumas definições formais. Deixe- P 1 (n) ser o n th termo da primeira sequência, P 2 (n) ser o n th termo da segunda sequência, e semelhante P k (n) o N ° termo do k th sequência tal como acima definido .

P 1 = [P 1 (1), P 1 (2), P 1 (3), P 1 (4), P 1 (5), ...]

P 2 = [(P 1 (1) + P 1 (2)) / 2, (P 1 (2) + P 1 (3)) / 2, (P 1 (3) + P 1 (4)) / 2, (P 1 (4) + P 1 (5)) / 2, ...]

P 3 = [(P 1 (1) + 2P 1 (2) + P 1 (3)) / 4, (P 1 (2) + 2P 1 (3) + P 1 (4)) / 4, (P 1 (3) + 2P 1 (4) + P 1 (5)) / 4, ...]

P 4 = [(P 1 (1) + 3P 1 (2) + 3P 1 (3) + P 1 (4)) / 8, (P 1 (2) + 3P 1 (3) + 3P 1 (4) + P 1 (5)) / 8, ...]

Não surpreendentemente, esses coeficientes seguem exatamente os coeficientes binomiais e podem ser expressos como uma única linha do triângulo de Pascal. Como uma linha arbitrária do Triângulo de Pascal é trivial de calcular, uma série arbitrariamente 'profunda' pode ser encontrada, simplesmente tomando as primeiras n somas parciais, multiplicando cada pelo termo correspondente na k- ésima linha do Triângulo de Pascal e dividindo por 2 k-1 .

Dessa forma, a precisão total do ponto flutuante de 32 bits (~ 14 casas decimais) pode ser alcançada com apenas 36 iterações, momento em que as somas parciais nem convergiram para a segunda casa decimal. Obviamente, isso não é um jogo de golfe:

# used for pascal's triangle
t = 36; v = 1.0/(1<<t-1); e = 1
# used for the partial sums of pi
p = 4; d = 3; s = -4.0

x = 0
while t:
  t -= 1
  p += s/d; d += 2; s *= -1
  x += p*v
  v = v*t/e; e += 1

print "%.14f"%x

Se você queria precisão arbitrária, isso pode ser alcançado com uma pequena modificação. Aqui, mais uma vez, calculando 1000 dígitos:

# used for pascal's triangle
f = t = 3318; v = 1; e = 1
# used for the partial sums of pi
p = 4096*10**999; d = 3; s = -p

x = 0
while t:
  t -= 1
  p += s/d; d += 2; s *= -1
  x += p*v
  v = v*t/e; e += 1

print x>>f+9

O valor inicial de p começa 2 10 maior, para neutralizar os efeitos da divisão inteira de s / d à medida que d se torna maior, fazendo com que os últimos dígitos não convergam. Observe aqui novamente que 3318também é:

dígitos * log 2 (10)

O mesmo número de iterações que o primeiro algoritmo (dividido pela metade porque t diminui em 1 em vez de 2 em cada iteração). Mais uma vez, isso indica uma convergência linear: um bit binário de pi por iteração. Em ambos os casos, são necessárias 3318 iterações para calcular 1000 dígitos de pi , como uma cota ligeiramente melhor que 1 milhão de iterações para calcular 5.

primo
fonte
Isso é muito melhor do que a minha solução:4 * sum(1/(1+i*2) if not i%2 else -1/(1+i*2) for i in xrange(places*10**(places)))
Aaron Hall
1
Isso é muito semelhante à minha abordagem , que por acaso é uma forma diferente da sua. No meu, como k → ∞, se f(-1,k)aproxima da sua soma de Euler.
Simply Beautiful Art
1
Muito legal; análise e explicação impressionantes, obrigado.
Jeremy radcliff
Apenas uma coisa pequena. Você não quis dizer depois de P_1 = ..., P_2 = ..., P_3 = ..., P_4 = ..."... multiplique cada um pelo termo correspondente na kthlinha do Triângulo de Pascal e divida por 2^{k-1}.", Em vez de nthlinha e 2^{n-1}?.
precisa saber é o seguinte
@jeremyradcliff eu fiz, sim. Obrigado pela correção.
Primo
5

Mathematica 42 39 34 33 31 26 32

Abordagem de Arquimedes 26 caracteres

N@#*Sin[180 Degree/#]&

Isso atinge o critério quando a entrada é 822.

Pergunta: Alguém sabe como ele calculou o pecado de 180 graus? Eu não.


Abordagem de Leibniz (série de Gregory) 32 caracteres

Esta é a mesma função que o poser do problema deu como exemplo. Atinge o critério em aproximadamente meio milhão de iterações.

N@4Sum[(-1)^k/(2k+1),{k,0,10^6}]

Abordagem Madhava-Leibniz 37 caracteres

Essa variação usa mais alguns caracteres, mas converge para o critério em apenas 9 iterações!

N@Sqrt@12 Sum[(-1/3)^k/(2k+1),{k,0,9}]
DavidC
fonte
todos os computam pelo algoritmo fornecido na definição do problema?
Acolito
A abordagem de @acolyte Leibniz (agora a primeira listada) é de fato a mencionada na descrição do problema. É muito lento para convergir. Uma ligeira variação nele (Madhava-Leibniz) converge muito rapidamente.
DavidC
Seno de 180 ° é bem fácil. É 180 ° / N que pode ficar complicado fora dos suspeitos do costume para N.
JB
Please explain, @J.B. Tricky to measure?
DavidC
This entry should state "32" because only Leibniz' approach fulfils the requirements (counting the characters in the code as given, I get 34, but both spaces may be safely removed, giving indeed a length of 32).
celtschk
4

APL (14)

4×-/÷1-⍨2×⍳1e6

 

marinus
fonte
1
13, --/4÷1-2×⍳1e6
Timtech
4

Java (67 chars)

float r(){float p=0,s=4,i=1E6f;while(--i>0)p+=(s=-s)/i--;return p;}

Note that this avoids loss of significance by adding the numbers up in the correct order.

Peter Taylor
fonte
this is fully compliant C code too. if posted as C, you could change while(--i>0) to while(i--) and save 2 chars
ardnew
1
@ardnew, true, but with C there are much more interesting tricks to play...
Peter Taylor
4

Haskell, 32

foldr(\k->(4/(2*k+1)-))0[0..8^7]

GHCi> foldr(\k->(4/(2*k+1)-))0[0..8^7]
3.141593130426724

Counting a function name it's

34

π=foldr(\k->(4/(2*k+1)-))0[0..8^7]
ceased to turn counterclockwis
fonte
3

R - 25 chars

sum(c(4,-4)/seq(1,1e6,2))
flodel
fonte
3

C (GCC) (44 chars)

float p(i){return i<1E6?4./++i-p(++i):0;}

That's 41 chars, but it also has to be compiled with -O2 to get the optimiser to eliminate the tail recursion. This also relies on undefined behaviour with respect to the order in which the ++ are executed; thanks to ugoren for pointing this out. I've tested with gcc 4.4.3 under 64-bit Linux .

Note that unless the optimiser also reorders the sum, it will add from the smallest number, so it avoids loss of significance.

Call as p().

Peter Taylor
fonte
Your recursive call is q(), not p(). And I don't think -O2 should be counted (but if you do count it, it's 4 chars because of the required space).
ugoren
Also: 1. gcc 4.1.1 doesn't optimize the recursion (and I don't see how it could), so the stack overflows. 2. it should be called as p(0). 3. Save a char by return++i.... 4. Two ++i makes undefined behavior.
ugoren
@ugoren, thanks for your comments. In order: q - that'll teach me to double-check after renaming. I think I'm following normal practice in counting -O2 as 3 chars, but we can open it up on meta if you want; meta.codegolf.stackexchange.com/questions/19 is the only relevant discussion I can find. I've added the version of gcc which I'm using, and which allows me to call it as p(). Saving the char stops the optimiser and gives segfault. I will clarify that I'm using undefined behaviour, as per meta.codegolf.stackexchange.com/questions/21
Peter Taylor
I added an answer to the meta question about flags. About p() - are you sure calling p() from any context would work? Or is it just what happened to be on the stack in your test?
ugoren
@ugoren, maybe I got lucky consistently. Even if I call it twice in a row, the second one still returns the correct value. gcc does seem to produce slightly different code for p() vs p(0), but I don't know what behaviour it documents and I'm not really a C programmer.
Peter Taylor
3

J, 26 chars

+/+/_2((4 _4)&%)>:+:i.100

Moved from 100 items of sequence to 1e6 items. Also now it's a code tagged and could be copypasted from browser to the console without errors.

+/+/_2((4 _4)&%)\>:+:i.1e6
fftw
fonte
3
-/4%>:2*i.1e6 -- 13 characters. (Thanks to b_jonas in #jsoftware for making me realise that -/ works to compute a sum with alternating sign. [This is since all operators in J are of equal precedence and right-associative, so -/ 1 2 3 4 <=> 1 - (2 - (3 - 4)) <=> 1 - 2 + 3 - 4.])
FireFly
that's neat and twice as awesome. Or even 2^10 more awesome!
fftw
@FireFly that is beautiful
Jonah
2

Javascript - 33 Characters

p=x=>4*(1-(x&2))/x+(x>1?p(x-2):0)

Call p passing a positive odd number x and it will calculate Pi with (x-1)/2 terms.

MT0
fonte
2

Ruby - 82 chars

def f(n,k=n)k>0?(f(n,k-1)+f(n+1,k-1))/2:n<0?0:f(n-1,0)+(-1)**n/(2*n+1.0)end;4*f(9)

Try it : https://repl.it/LQ8w

The approach uses the given series indirectly using a numerical acceleration approach. The resulting output is

pi ≈ 3.14159265161

vs.

pi = 3.14159265359

It starts with

f(n,0) = 1/1 - 1/3 + 1/5 - ... + ((-1)**n)/(2*n+1)

And then, since this is alternating, we can accelerate the convergence using

f(n,1) = (f(n,0) + f(n+1,0))/2

And it repeatedly applies this:

f(n,k) = (f(n,k-1) + f(n+1,k-1))/2

And for simplicity, f(n) = f(n,n).


Ruby - 50 chars

If you don't mind running for a really long while, then you could simply use

def f(n)n<0?0:f(n-1)+(-1)**n/(2*n+1.0)end;4*f(1e7)

or

a=0;for k in 0..1e7 do a+=(-1)**k/(2*k+1.0)end;4*a
Simply Beautiful Art
fonte
1

C, 69 chars

float p,b;void main(a){b++<9e6?p+=a/b++,main(-a):printf("%f\n",4*p);}
  • Run with no command line parameters (so a is initialized to 1).
  • Must be compiled with optimization.
  • void main is strange and non-standard, but makes things work. Without it, the recursion is implemented as a real call, leading to a stack overflow. An alternative is adding return.
  • Two characters 4* can be saved, if running with three command line parameters.
ugoren
fonte
You could shorten that to int main(a) or even main(a), GCC only gives a warning. And it will give a warning for void main anyway, and maybe even because you have only one argument to main.
nyuszika7h
1

Clojure - 79 chars

(fn [](* 4(apply +(map #(*(Math/pow -1 %1)(/ 1.0(+ 1 %1 %1)))(range 377000)))))

This creates a function of no arguments which will calculate a float which approximates pi correctly to five decimal places. Note that this does not bind the function to a name such as pi, so this code must either be evaluated in place with eval as (<code>) or bound to a name in which case the solution is

(defn p[](* 4(apply +(map #(*(Math/pow -1 %1)(/ 1.0(+ 1 %1 %1)))(range 377000)))))

for 82 chars

About

(defn nth-term-of-pi [n] (* (Math/pow -1 n) (/ 1.0 (+ 1 n n))))
(defn pi [c] (* 4 (apply + (map nth-term-of-pi (range c)))))
(def  pi-accuracy-constant (loop [c 1000] (if (< (pi c) 3.14159) (recur (inc c)) c)))
; (pi pi-accuracy-constant) is then the value of pi to the accuracy of five decimal places
arrdem
fonte
1

PHP - 56 55 chars

<?for($j=$i=-1;1e6>$j;){$p+=($i=-$i)/($j+=2);}echo$p*4;

I don't know that I can get it much smaller without breaking the algorithm rule.

TwoScoopsofPig
fonte
1
How about this for 45? <?for(;1e6>$j;)$p+=($i=-$i|4)/~-$j+=2;echo$p;
primo
I was trying to come up with that, but couldn't get the bitwise ops to work. Thanks for the suggestion!
TwoScoopsofPig
You can remove the last semicolon to save 1 character.
nyuszika7h
1

Perl - 43 39 chars

not sure the rules on anonymous subroutines, but here's another implementation using @FireFly's series construction

sub{$s+=8/((4*$_+2)**2-1)for 0..1e6;$s}

sub p{$s+=(-1)**$_*4/(2*$_+1)for 0..1e6;$s}

ardnew
fonte
0

Java - 92 84 chars

I cannot beat by far Peter Taylor's result, but here is mine:

double d(){float n=0,k=0,x;while(n<9E5){x=1/(1+2*n++);k+=(n%2==0)?-x:x;}return 4*k;}

Ungolfed version:

double d() {
    float n = 0, k = 0, x;
    while (n < 9E5) {
        x = 1 / (1 + 2 * n++);
        k += (n % 2 == 0) ? -x : x;
    }
    return 4 * k;
}

Edit: Saved a few characters using ternary operator.

Averroes
fonte
0

Python - 56 chars

Meh, My python-fu is not strong enough. I couldn't see any more shortcuts but maybe a more experienced golfer could find something to trim here?

t=s=0
k=i=1
while t<1e6:t,s,i,k=t+1,k*4./i+s,i+2,-k
chucksmash
fonte
You could use Python 3 to save one byte for the float division (4. -> 4). In other news, I just found a case where Python 3 actually beats Python 2 in code golf!
nyuszika7h
0

Ruby - 54 chars

def a()p=0;1000000.times{|i|p+=8/(4*i*(4*i+2))};p;end;

My first try on console

def a()i=1;p=0;while i<2**100 do p+=8/(i*(i+2));i+=4;end;p;end;

63 chars.

Perello
fonte
You can save a byte by using def a; instead of def a().
nyuszika7h
Another one by removing the last semicolon.
nyuszika7h
0

Perl (76 chars)

$y=1e4;for$x(0..1e4-1){$y--while sqrt($x**2+$y**2)>1e4;$a+=$y}print 4*$a/1e8

(Result: 3.14159052)

Not the shortest possible solution, but maybe interesting. It's a geometrical one. I calculate the area under a circle.

I got another funny approach, but it's really slow. It counts the number of discrete points in a square that are below a quarter circle and calculates pi from it:

$i=shift;for$x(0..$i){for$y(0..$i){$h++if sqrt($x**2+$y**2)<$i}}print$h*4/$i**2

It expects the number of iterations as command line argument. Here you can see how run time relates to accuracy. ;)

$ time perl -e '$i=shift;for$x(0..$i){for$y(0..$i){$h++if sqrt($x**2+$y**2)<$i}}print$h*4/$i**2' 100
3.1796
real    0m0.011s
user    0m0.005s
sys 0m0.003s

$ time perl -e '$i=shift;for$x(0..$i){for$y(0..$i){$h++if sqrt($x**2+$y**2)<$i}}print$h*4/$i**2' 1000
3.14552
real    0m0.354s
user    0m0.340s
sys 0m0.004s

$ time perl -e '$i=shift;for$x(0..$i){for$y(0..$i){$h++if sqrt($x**2+$y**2)<$i}}print$h*4/$i**2' 10000
3.14199016
real    0m34.941s
user    0m33.757s
sys 0m0.097s
memowe
fonte
0

k (25 chars)

4*+/%(i#1 -1)'1+2!i:1000000

Slightly shorter:

+/(i#4 -4)%1+2*!i:1000000
skeevey
fonte
0

Python (49)

print 4*sum((-1)**i/(2*i+1.)for i in range(9**6))
3.14159453527
arshajii
fonte
0

CJam - 21

1e6{WI#4d*I2*)/}fI]:+

Pretty straightforward calculation of the given series.
CJam is http://sf.net/p/cjam

aditsu
fonte
0

Julia - 30 characters

sum(4./[1:4:1e6] - 4./[3:4:1e6])
Milktrader
fonte
0

SQL, 253 bytes

DECLARE @B int=3, @A varchar(max), @C varchar(max)='1'
WHILE @B<100000
BEGIN
SELECT @C=@C+(select case when (@B-1)%4=0 then'+'else'-'end)+
(SELECT cast(cast(1.0/@B as decimal(9,8)) as varchar(max)))
SELECT @B=@B+2
END
EXECUTE('SELECT 4*('+@C+')')

I would provide a SQL Fiddle, but this goes too many loops deep finding the 1/3 1/5 1/7 etc. fractions and gives errors lol. However, if you change @B<100000 to 1000 then it runs (obviously not to the same number of digits of accuracy).

phroureo
fonte
0

Befunge, 129 bytes

p08p109p^v*86%+55:<$$$<
\$\>:#,_@>+\55+/:#^_"."
v>p"~"/:"~"%08p"~"/00p:2\4%-*"(}"
8^90%"~":+2:+g90*+g80*<
>*:**\/+>"~~"00g:"~"`!|

Try it online!

In case anyone is wondering, it's an elephant.

Cartoon elephant

James Holderness
fonte