Output first \$n\$ digits of \$\pi^{1/\pi}\$

Question

This challenge is to produce the shortest code for the constant \$\pi^{1/\pi}\$. Your code must output the first \$n\$ consecutive digits of \$\pi^{1/\pi}\$, where \$n\$ is given in the input. Alternatively, your program may accept no input, and output digits indefinitely

This is code golf, so the shortest submission (in bytes) wins except that it must output the 1000 digits for \$n = 1000\$ in less than 10 seconds on a reasonable PC.

You may not use a built-in for \$\pi\$, the gamma function or any trigonometic functions.

If the input is \$n = 1000\$ then the output should be:

1.439619495847590688336490804973755678698296474456640982233160641890243439489175847819775046598413042034429435933431518691836732951984722119433079301671110102682697604070399426193641233251599869541114696602206159806187886346672286578563675195251197506612632994951137598102148536309069039370150219658973192371016509932874597628176884399500240909395655677358944562363853007642537831293187261467105359527145042168086291313192180728142873804095866545892671050160887161247841159801201214341008864056957496443082304938474506943426362186681975818486755037531745030281824945517346721827951758462729435404003567168481147219101725582782513814164998627344159575816112484930170874421446660340640765307896240719748580796264678391758752657775416033924968325379821891397966642420127047241749775945321987325546370476643421107677192511404620404347945533236034496338423479342775816430095564073314320146986193356277171415551195734877053835347930464808548912190359710144741916773023586353716026553462614686341975282183643

Note the 1000 decimal places includes the first \$1\$.

Output format

Your code can output in any format you wish.

This related question tackles a simpler variant but with the same restrictions.

Notes

If it's helpful, you can assume \$\pi^{1/\pi}\$ is irrational and in fact you can even assume it is a Normal number.

Answer 1

Python, 149 bytes

Saved one byte due to @H.PWiz.

k=n=int(input())+1
i=j=8*n
e=l=p=z=10**n
while i:p=2*z-i//2*p//~i;i-=2
q=z-z*z//p
while j:l=q//j+q*l//z;j-=1
while k:e=z+e*l//k//p;k-=1
print(e//100)

Try it online!

Input n = 1000 finishes in less than 0.5s.

\$\sqrt[\pi]{\pi}\$ is calculated as the \$e^\frac{\ln(\pi)}{\pi}\$.\$\pi\$ is calculated with the usual Euler-Leibniz: \$\pi=\sum_{n=0}\limits^{\infty}{\frac{n!}{(2n+1)!!}}\$.

\$e^x\$ and \$\ln(x)\$ are both computed using a Taylor series: \$e^x=\sum\limits_{n=0}^{\infty}{\frac{x^n}{n!}}\$ and \$\ln(1-x)=-\sum\limits_{n=1}^{\infty}{\frac{x^n}{n}}\$. Because the iteration for \$\ln(x)\$ converges only when \$|x|<1\$, this is instead calculated as \$\ln(\pi)=-\ln(\frac{1}{\pi})\$.

Given that Python uses Karasuba Multiplication, the overall runtime complexity is \$\mathcal{O}(n^{1+\log_2(3)})\$ - in other words, twice as many digits will take approximately 6 times as long.

---

Subquadratic Complexity

import sys
from gmpy2 import isqrt, mpz

def piks(a, b):
  if a == b:
    if a == 0:
      return (1, 1, 1123)
    p = a*(a*(32*a-48)+22)-3
    q = a*a*a*24893568
    t = 21460*a+1123
    return (p, -q, p*t)
  m = (a+b) >> 1
  p1, q1, t1 = piks(a, m)
  p2, q2, t2 = piks(m+1, b)
  return (p1*p2, q1*q2, q2*t1 + p1*t2)

n = int(sys.argv[1])-1
m = n*20//3
z = mpz(10)**n

# n / log(777924, 10)
pi_terms = mpz(n*0.16975227728583067)

pp, pq, pt = piks(0, pi_terms)
pq *= 3528

pi2m = (pq << m) // pt

a, b = 2 << m, 8
while a != b:
  a, b = (a + b) >> 1, isqrt(a*b)

mlog2_pi = (z << m) // a

a, b = 2*pi2m, 8
while a != b:
  a, b = (a + b) >> 1, isqrt(a*b)
logpi_pi = z * pi2m // a - mlog2_pi

mlog2 = mlog2_pi * pq // pt

d = e = (15044673 << m) // 10450451
pt //= z
while d:
  a, b = 2*e, 8
  while a != b:
    a, b = (a + b) >> 1, isqrt(a*b)
  lnx = (pq * e) // (pt * a) - mlog2
  d = e * (lnx - logpi_pi) // z
  e -= d

print(e * z >> m)

Try it online!

Input n = 20000 finishes in less than one second.

\$\pi\$ can be computed in subquadratic time by use of Karatsuba splitting a.k.a. Fast E-function Evaluation, which reduces a summation to n terms to a single rational value p/q, splitting in binary descent. I've chosen to use Ramanujan #39, which is the fastest converging series of its kind that doesn't require an arbitrary precision square root, to my knowledge. Explicitly, this is computed as:

\$\large{\left.{3528}\middle/{\sum\limits_{n=0}^{\infty}\frac{(-1)^n (4n)! (1123+21460n)}{(n!)^4 14112^{2n}}}\right.}\$

For \$\ln(x)\$ and \$e^x\$, using the same technique wouldn't provide any benefit, because both are computed as a power series of an arbitrary precision variable. Fortunately, Gauss has gifted us with an elegant quadratically converging formula for \$\ln(x)\$, based on the Arithmetic-Geometric Mean:

\$\DeclareMathOperator{\AGM}{AGM}\large\ln(x)=\lim\limits_{n\rightarrow\infty}\frac{\pi x^n}{2n\AGM(x^n,4)}\$

or, in cases when large powers of x are inconvenient, this can also be computed as:

\$\large\ln(x)\approx\frac{\pi x 2^m}{2\AGM(x 2^m,4)}-m\log(2)\$

Conveniently, division by \$\pi\$ can be acheived simply by not multiplying through by \$\pi\$.

With the natural logarithm thusly defined, \$e^\frac{\ln(\pi)}{\pi}\$ can be computed via Newton's method on \$x_{n+1}=x_n-x_n(\ln(x_n)-\frac{\ln(\pi)}{\pi})\$. The overall complexity is then \$\mathcal{O}(n^*\log^3(n))\$, where \$n^*\$ will vary with the complexity of the multiplication algorithm GMP is using for any given bit length.

Answer 2

Wolfram Language (Mathematica), 62 bytes

In my first try I used Zeta function but this one uses the imaginary part of ln(-1) for pi. (@someone)
Prints more than 4000 digits in the first 10 seconds,

q=Im@Log@-1;Do[Print@RealDigits[N[q^(1/q),2k]][[1,k]],{k,∞}]

Try it online!

9 bytes saved from @someone

Answer 3

Wolfram Language (Mathematica), 23 22 bytes

(p=Log@-1/I)^p^-1~N~#&

Try it online!

^{-1 byte by @someone — the use and reuse of p beats the pure function}

Taking a lot of inspiration from @J42161217 and displaying n digits instead of an infinite stream.

Note that I'm using ToString in Tio to suppress the trailing precision specifier in the output, which does not appear when this code is executed in Mathematica.

Maybe using 180° for pi would work; but maybe that's too close to being trigonometric and is disallowed.

Answer 4

AXIOM, 221 bytes

m(k)==(v:=8.*k;(4/(v+1)-2/(v+4)-1/(v+5)-1/(v+6))*16^-k)
p(n:PI):String==(d:=digits(n+9);e:=10.^-digits();i:=s:=0;repeat(k:=m(i);k<e=>break;s:=s+k;i:=i+1);r:=concat split((s^(1/s))::String,char " ");digits(d);r.(1..(n+1)))

test and ungolf:

(3) -> p 1000
   (3)
  "1.43961949584759068833649080497375567869829647445664098223316064189024343948
  91758478197750465984130420344294359334315186918367329519847221194330793016711
  10102682697604070399426193641233251599869541114696602206159806187886346672286
  57856367519525119750661263299495113759810214853630906903937015021965897319237
  10165099328745976281768843995002409093956556773589445623638530076425378312931
  87261467105359527145042168086291313192180728142873804095866545892671050160887
  16124784115980120121434100886405695749644308230493847450694342636218668197581
  84867550375317450302818249455173467218279517584627294354040035671684811472191
  01725582782513814164998627344159575816112484930170874421446660340640765307896
  24071974858079626467839175875265777541603392496832537982189139796664242012704
  72417497759453219873255463704766434211076771925114046204043479455332360344963
  38423479342775816430095564073314320146986193356277171415551195734877053835347
  93046480854891219035971014474191677302358635371602655346261468634197528218364
  3"
                                                             Type: String
           Time: 0.03 (IN) + 0.53 (EV) + 0.25 (OT) + 0.15 (GC) = 0.97 sec


-- Bailey-Borwein-Plouffe formula for pi
-- https://en.m.wikipedia.org/wiki/Bailey%E2%80%93Borwein%E2%80%93Plouffe_formula   
--       +oo
--      -----          --                                      --
--      \          1   |     4         2         1         1    |
--   pi= |      -------|  ------- - ------- - ------- - ------- |
--      /            k |   8*k+1     8*k+4     8*k+5     8*k+6  |
--      ----- k   16   --                                      --
--        0
m(k)==(v:=8.*k;(4/(v+1)-2/(v+4)-1/(v+5)-1/(v+6))*16^-k)
p(n:PI):String==
      d:=digits(n+9); e:=10.^-digits(); i:=s:=0
      repeat
          k:=m(i)
          k<e=>break
          s:=s+k;i:=i+1
      r:=concat split((s^(1/s))::String,char " ")
      digits(d)
      r.(1..(n+1))

the function m calculate the term of the sum; the p() function loop, sum them in the variable s until the term is < than min float value (one can see as epsilon); p() function return one string of that number pi^(1/pi); it seems the required digits are returned in less than 1 second for input to p() function 1000.