​P​i​ =​= ​3​.​2​

Inspired by this video of Infinite Series.

Introduction

Pi is defined as the ratio of the circumference to the diameter of a circle. But how is a circle defined? Usually a circle is defined as the points with constant distance to the centerpoint (let us assume that the center is at (0,0)). The next question would be: How do we define the distance? In the following we are considering different notions of distances (induced by the Lp-norms):

Given a norm (=something that measures a length) we can easily construct a distance (=distance between two points) as follows:

dist(A,B) := norm (A-B)


The euclidean norm is given by:

norm((x,y)) = (x^2 + y^2)^(1/2)


This is also called the L2-norm. The other Lp-norms are constructed by replacing the 2 in the above formula by other values between 1 and infinity:

norm_p((x,y)) = (|x|^p + |y|^p)^(1/p)


The unit circles for those different norms have quite distinct shapes:

Challenge

Given a p >= 1, calculate the ratio of circumference to diameter of a Lp-circle with respect to the Lp-norm with an accuracy of four significant figures.

Testcases

We can use that for p,q with 1 = 1/p + 1/q we get the same ratio for the Lp as well as the Lq norm. Furthermore for p = q = 2 the ratio is minimal, and for p = 1, q = infinity we get a ratio of 4, so the ratios are always between pi and 4.

p   or  q            ratio
1       infinity     4
2       2            3.141592
1.623   2.60513      3.200
1.5     3            3.25976
4       1.33333      3.39693

• The shapes are known as Lamé curves or superellipses and exist for 0 < p < 1 too, even though the norm itself doesn't (because it violates the triangle inequality). The Wikipedia article for the superellipse includes a closed form for the area. – Neil Jan 7 '17 at 23:31
• @Neil We do however need to consider the circumference, not the area, which - as far as I know - can only be calculated via a arc length integral. – flawr Jan 7 '17 at 23:33
• Sorry, by the time I'd finished reading up on them I'd forgotten what the question had asked for. – Neil Jan 7 '17 at 23:51
• Lovely challenge! – Luis Mendo Jan 8 '17 at 4:16
• It's interesting to note that the area formula (A = πr²) does not hold for p ≠ 2 – Mego Jan 8 '17 at 13:08

Python + scipy, 92 bytes

from scipy.integrate import*


Formula is from this math.SE question.

• When testing an implementation with this appraoch I had trouble with the convergence of that approach, due to the singularity at x=1, how does your submission do? – flawr Jan 7 '17 at 23:12
• Scipy is not part of the Python standard library. Maybe switch to Sage? – busukxuan Jan 7 '17 at 23:13
• @busukxuan There is no requirement on PPCG that allows you to only use standard libraries. But I'll mention it in the title anyways. – orlp Jan 7 '17 at 23:14
• @ChristianSievers I did my own integration to avoid feeling bad for using someone else's closed formula :-P – Luis Mendo Jan 8 '17 at 4:14
• @ChristianSievers I actually also included another formla in the sandbox, in case you're interested=) – flawr Jan 8 '17 at 10:43

MATL, 31 bytes

0:1e-3:1lyG^-lG/^v!d|G^!slG/^sE


Explanation

This generates the x,y coordinates of one quarter of the unit circle sampled at 1001 points with step 0.001 in x. The length of the quarter of circle is approximated by that of the polygonal line that passes through those points; that is, the sum of the lengths of the 1000 segments. Length is of course computed according to p-norm. Multiplying the result by 2 gives the approximate lenght of half a circle, that is, pi.

0:1e-3:1   % Push [0 0.001 0.002 ... 0.999 1]. These are the x coordinates of
% the vertices of the polygonal line that will approximate a quarter
% of the unit circle
l          % Push 1
y          % Duplicate [0 0.001 0.002 ... 0.999 1] onto the top of the stack.
G          % Push input, p
^          % Element-wise power: gives [0^p 0.001^p ... 1^p]
-          % Element-wise subtract from 1: gives [1-0^p 1-0.001^p ... 1-1^p]
lG/        % Push 1, push p, divide: gives 1/p
^          % Element-wise power: gives [(1-0^p)^(1/p) (1-0.001^p)^(1/p) ...
% ... (1-1^p)^(1/p)]. These are the y coordinates of the vertices
% of the polygonal line
v          % Concatenate vertically into a 2×1001 matrix. The first row contains
% the x coordinates and the second row contains the y coordinates
!          % Transpose
d|         % Compute consecutive differences down each column. This gives a
% 1000×2 matrix with the x and y increments of each segment. These
% increments will be referred to as Δx, Δy
G          % Push p
^          % Element-wise power
!          % Transpose
s          % Sum of each column. This gives a 1×1000 vector containing
% (Δx)^p+(Δy)^p for each segment
lG/        % Push 1/p
^          % Element-wise power. This gives a 1×1000 vector containing
% ((Δx)^p+(Δy)^p)^(1/p) for each segment, that is, the length of
% each segment according to p-norm
s          % Sum the lenghts of all segments. This approximates the length of
% a quarter of the unit circle
E          % Multiply by 2. This gives the length of half unit circle, that is,
% pi. Implicitly display


Mathematica, 49 46 bytes

3 bytes saved due to alephalpha.

2NIntegrate[(1+(a^-#-1)^(1-#))^(1/#),{a,0,1}]&


Anonymous function. Takes a number as input and returns a number as output.

• 2NIntegrate[(1+(a^-#-1)^(1-#))^(1/#),{a,0,1}]& – alephalpha Jan 8 '17 at 2:55

PARI/GP, 48 43 bytes

It's easy after @orlp has found the formula, and @alephalpha's version saves 5 bytes:

p->2*intnum(u=0,1,(1+(u^-p-1)^(1-p))^(1/p))


To add something slightly useful, let's calculate the p for which we get 3.2:

? f=p->2*intnum(u=0,1,(1+(u^-p-1)^(1-p))^(1/p));
? solve(p=1,2,f(p)-3.2)
%2 = 1.623002382384469009676324702


Correct usage

While the code gives results that are much more exact than the challenge demands, it can easily be improved a lot: if we replace the upper integration limit 1 with [1,1/p-1] (giving what the manual calls the singularity exponent) then all shown digits of f(2) agree with Pi. This is still true if we increase the precision to 100 (type \p100).

However, after that change the solve computation no longer worked. I changed the inner term to explicitely handle the case u=0 and also changed to a different computer with a newer PARI version and 64 bit (which implies a higher default precision).

Here is the improved calculation of the p value for Pi=3.2, and let's also have a look at the real Pi:

? f=p->2*intnum(u=0,[1,1/p-1],if(u,(1+(u^-p-1)^(1-p))^(1/p),0));
? f(2)
%2 = 3.1415926535897932384626433832795028842
? Pi
%3 = 3.1415926535897932384626433832795028842
? solve(p=1,2,f(p)-3.2)
%4 = 1.6230023823844690096763253745604419761

• p->2*intnum(u=0,1,(1+(u^-p-1)^(1-p))^(1/p)) – alephalpha Jan 8 '17 at 2:57

JavaScript (ES7), 80 bytes

Based on orlp's answer. This JS implementation is quite slow. You may want to try i=1e-7 (or even higher) for a faster approximation.

Note: This is basically intended for Chrome and Edge only. An equivalent ES6 version using Math.pow() on Firefox 50.1 seems to be much slower.

Edit: According to Neil, this should also work fine on Firefox 52.

f=
p=>{for(i=5e-8,s=x=0;(x+=i)<1;)s+=i*(x**(1-p)+(1-x)**(1-p))**(1/p);return 2/p*s}

console.log(f(1).toFixed(3))
console.log(f(2).toFixed(3))
console.log(f(1.623).toFixed(3))

• The ES7 version seemed pretty sprightly when I tried it using Firefox 52 (I didn't measure it scientifically, but it felt about the same speed as Chrome; Edge froze on me). – Neil Jan 15 '17 at 11:06
• @Neil Thanks for your feedback. Updated accordingly. – Arnauld Jan 15 '17 at 11:23