Square Line Picking

Author

Eric W. Weisstein
March 14, 2017

This notebook downloaded from http://mathworld.wolfram.com/notebooks/ComputationalGeometry/SquareLinePicking.nb.

For more information, see Eric's MathWorld entry http://mathworld.wolfram.com/SquareLinePicking.html.

©2017 Wolfram Research, Inc. except for portions noted otherwise

Additional Authors

Analytic length distribution derived by
Michael Trott
March 12, 2004

Points in Interior

RandomPoint

lines=RandomPoint[Rectangle[],{(10^7),2}];//Timing

Length[lens=EuclideanDistance@@@lines]//Timing

dist[l_]:=({{\[Piecewise], ({{2 l (-4 l+(l^2)+π), 0≤l≤1}, {2 l (-2-(l^2)+4 (Sqrt[-1+(l^2)])+π-4 ArcTan[(Sqrt[-1+(l^2)])]), 1<l≤(Sqrt[2])}})}})

Plot[dist[l],{l,0,(Sqrt[2])}]

With[{bins=100},Show[{Histogram[lens,bins],Plot[Evaluate[(1/bins)Length[lens]dist[s]],{s,0,(Sqrt[2])},PlotStyle->Opacity[1]]}]]

Mean[lens]

Monte Carlo

x=Norm/@Subtract@@@Table[Random[],{10^6},{2},{2}];//Timing

{Mean[x],Median[x],Variance[x]}

<<Graphics`Graphics`

BarChart[Length/@Split[Sort[Floor[100x]]],Ticks->None,BarSpacing->-.2,BarStyle->Red,BarEdges->False]

Length Distribution Derivation

Integrate[DiracDelta[l^2-Abs[(x1-x2)^2+(y1-y2)^2]],{x1,0,1}]

Integrate[DiracDelta[Expand[l-((((x1-x2))^2)+(((y1-y2))^2))]],{x1,0,1}]

Wrong...

Length Distribution

Attempt 1

(Integrate[(1/(2(Sqrt[(x1-x2)^2+(y1-y2)^2])))(DiracDelta[d-(Sqrt[(x1-x2)^2+(y1-y2)^2])]+DiracDelta[d+(Sqrt[(x1-x2)^2+(y1-y2)^2])]),{x2,0,1},Assumptions->0<d<Sqrt[2]&&0<y1<1&&0<y2<1])//Timing

Attempt 2

Integrate[DiracDelta[d^2-((x1-x2)^2+(y1-y2)^2)],{x1,0,1}]

Integrate[DiracDelta[(d^2)-(((x1-x2))^2)-(((y1-y2))^2)],{x1,0,1},Assumptions->!(x2+(Sqrt[(d^2)-(((y1-y2))^2)])<1&&-1<-x2+(Sqrt[(d^2)-(((y1-y2))^2)])<0)]

FullSimplify[!(-1<-x2+(Sqrt[(d^2)-(((y1-y2))^2)])<0&&x2+(Sqrt[(d^2)-(((y1-y2))^2)])≥1)||!(x2+(Sqrt[(d^2)-(((y1-y2))^2)])<1&&-1<-x2+(Sqrt[(d^2)-(((y1-y2))^2)])<0)]

Integrate[Piecewise[{{(1/(Sqrt[Abs[(d^2)-(((y1-y2))^2)]])),x2+(Sqrt[(d^2)-(((y1-y2))^2)])<1&&-1<-x2+(Sqrt[(d^2)-(((y1-y2))^2)])<0}},(1/(Sqrt[(d^2)-(((y1-y2))^2)]))],{x2,0,1},Assumptions->0<d<Sqrt[2]&&0<y1<1&&0<y2<1]//InputForm

Integrate[Piecewise[{{1/Sqrt[d^2-y1^2+2*y1*y2-y2^2],
d-y1+y2≥0&&d+y1-y2≥0}},0],{y1,0,1},Assumptions->0<d<Sqrt[2]&&0<y2<1]

Integrate[%,{y2,0,1},Assumptions->0<d<Sqrt[2]]

FullSimplify[%,0<d<Sqrt[2]]

<<Graphics`Colors`

<<PSDirect`

Plot[If[0<d<1,-2 d+π,2 (-d+(Sqrt[-1+(d^2)])+ArcCot[(Sqrt[-1+(d^2)])])],{d,0,Sqrt[2]},PlotRange->All]

Attempt 3

i=Integrate[DiracDelta[d^2-((x1-x2)^2+(y1-y2)^2)],{x1,0,1},Assumptions->0<d<Sqrt[2]&&0<x2<1&&0<y1<1&&0<y2<1]

i[[3]]

Integrate[(1/(Sqrt[Abs[(d^2)-(((y1-y2))^2)]]))(Boole[x2+(Sqrt[(d^2)-(((y1-y2))^2)])<1&&-1<-x2+(Sqrt[(d^2)-(((y1-y2))^2)])<0]+Boole[-1<-x2+(Sqrt[(d^2)-(((y1-y2))^2)])<0&&x2+(Sqrt[(d^2)-(((y1-y2))^2)])≥1]),{x2,0,1},Assumptions->0<d<Sqrt[2]&&0<y1<1&&0<y2<1]//FullSimplify

Integrate[%,{y1,0,1},Assumptions->0<d<Sqrt[2]&&0<y2<1]//FullSimplify

Integrate[%,{y2,0,1},Assumptions->0<d<Sqrt[2]]//Simplify

Plot[%36,{d,0,Sqrt[2]},PlotRange->All]

Attempt 4

Integrate[Boole[d^2==(x1-x2)^2+(y1-y2)^2],{x1,0,1},Assumptions->0<d<Sqrt[2]&&0<x2<1&&0<y1<1&&0<y2<1]

Integrate[Boole[x==y],{x,0,1}]

Integrate[Piecewise[{{1,x==y}}],{x,0,1}]

Attempt 5: Michael's Version

This is correct for 0≤l≤(Sqrt[2])except l==1

p[λ_]:=UnitStep[λ] UnitStep[1-λ] 2 λ (Pi+λ (λ-4))+UnitStep[λ-1] UnitStep[Sqrt[2]-λ]*(2 λ (Pi-2-λ^2+4 Sqrt[λ^2-1]-4 ArcTan[Sqrt[λ^2-1]]))

PiecewiseExpand[p[l]]

{2 l (-4 l+(l^2)+π),
-2 l (2+(l^2)-4 (Sqrt[-1+(l^2)])-π+4 ArcTan[(Sqrt[-1+(l^2)])]),
-4 (l+2 (l^2)-2 l (Sqrt[-1+(l^2)])-l π+2 l ArcTan[(Sqrt[-1+(l^2)])]),
p[l]}/.l->1//Simplify

Limit[{2 l (-4 l+(l^2)+π),
-2 l (2+(l^2)-4 (Sqrt[-1+(l^2)])-π+4 ArcTan[(Sqrt[-1+(l^2)])]),
-4 (l+2 (l^2)-2 l (Sqrt[-1+(l^2)])-l π+2 l ArcTan[(Sqrt[-1+(l^2)])]),
p[l]},l->1]

p[l_]:=Piecewise[{
{2l((l^2)-4l+π),0≤l≤1},{2l(4(Sqrt[(l^2)-1])-((l^2)+2-π)-4ArcTan[(Sqrt[(l^2)-1])]),1<l<=(Sqrt[2])}}]

p[l]

Integrate[p[l],{l,0,x},Assumptions->0<x<1]

Integrate[p[l],{l,0,x},Assumptions->1<x<(Sqrt[2])]

Expand/@Integrate[p[l],{l,0,x},Assumptions->0<x<(Sqrt[2])]

SquareLinePicking /: PDF[SquareLinePicking,x_]:=(Piecewise[{{2 l (-4 l+(l^2)+π), 0≤l≤1}, {2 l (-2-(l^2)+4 (Sqrt[-1+(l^2)])+π-4 ArcTan[(Sqrt[-1+(l^2)])]), 1<l≤(Sqrt[2])}}])
SquareLinePicking /:CDF[SquareLinePicking,x_]:=(Piecewise[{{π (x^2)-((8 (x^3))/3)+((x^4)/2), 0<x≤1}, {(1/3)-2 (x^2)+π (x^2)-((x^4)/2)+(4/3) (Sqrt[-1+(x^2)])+(8/3) (x^2) (Sqrt[-1+(x^2)])-4 (x^2) ArcTan[(Sqrt[-1+(x^2)])], 1<x<(Sqrt[2])}}])
SquareLinePicking /:Domain[SquareLinePicking]:=Interval[{0,Sqrt[2]}]

Show[GraphicsArray[
Block[{$DisplayFunction=Identity},Plot[#[SquareLinePicking,l],{l,0,Sqrt[2]},PlotStyle->Red,AxesLabel->TraditionalForm/@{l,#[l]}]&/@{PDF,CDF}
]],ImageSize->500]

Raw Moments

General

Integrate[l^n p[l],{l,0,Sqrt[2]}]

Integrate[l^n 2 l (-4 l+(l^2)+π),{l,0,1},Assumptions->n≥0]

Integrate[l^n 2 l (-2-(l^2)+4 (Sqrt[-1+(l^2)])+π-4 ArcTan[(Sqrt[-1+(l^2)])]),{l,1,(Sqrt[2])},Assumptions->n≥0]

Integrate[l^n 2 l (-2-(l^2)+4 (Sqrt[-1+(l^2)])+π),{l,1,(Sqrt[2])},Assumptions->n≥0]

Integrate[2 (l^(1+n)) 4 ArcTan[(Sqrt[-1+(l^2)])],{l,1,(Sqrt[2])},Assumptions->n≥0]

Integrate[2 (l^(1+n)) 4 ArcTan[(Sqrt[-1+(l^2)])],l]

Table[Integrate[l^n p[l],{l,0,Sqrt[2]}],{n,8}]

Normalization

Integrate[p[l],{l,0,Sqrt[2]}]

Mean

μ=With[{n=1},Integrate[l^n p[l],{l,0,Sqrt[2]}]]

Central Moments

Integrate[(l-μ)^n p[l],{l,0,Sqrt[2]}]

Integrate[(l-μ)^n p[l],{l,0,Sqrt[2]},Assumptions->n≥0]

Table

Table[Integrate[(l-μ)^n p[l],{l,0,Sqrt[2]}],{n,4}]

FullSimplify[%]

N[%]

Variance

{#,N[#]}&@Integrate[(l-μ)^2p[l],{l,0,Sqrt[2]}]//FullSimplify

Median

f[x_?NumericQ]:=NIntegrate[p[l],{l,0,x}]

FindRoot[f[x]-1/2,{x,.5},WorkingPrecision->20]

d=Integrate[p[l],{l,0,x}]

FullSimplify[d]

Expand[(1/6) (2-3 (x^4)+8 (Sqrt[-1+(x^2)])+16 (x^2) (Sqrt[-1+(x^2)])+6 (x^2) (-2+π-4 ArcTan[(Sqrt[-1+(x^2)])]))]

Expand[(1/6) (x^2) (6 π+x (-16+3 x))]

FindRoot[π (x^2)-((8 (x^3))/3)+((x^4)/2)==(1/2),{x,1/2},WorkingPrecision->20]

Mean via Direct Integration

((Subscript[x, 1]),(Subscript[y, 1]),(Subscript[x, 2]),(Subscript[y, 2]))

i=Integrate[Sqrt[(x1-x2)^2+(y1-y2)^2],{x1,0,1},{y1,0,1},{x2,0,1},{y2,0,1}];//Timing

i

...probably wrong...

i/.Integrate[x__]:>NIntegrate[Evaluate[x],Method->QuasiMonteCarlo]

...uneavaluated with no messages...

v4.2

Integrate[Sqrt[(x1-x2)^2+(y1-y2)^2],{x1,0,1},{x2,0,1},{y1,0,1},{y2,0,1}]//Timing

Integrate[Sqrt[(x1-x2)^2+(y1-y2)^2],{x1,0,1},{y1,0,1},{x2,0,1},{y2,0,1}]//Timing

Integrate[Sqrt[(x1-x2)^2+(y1-y2)^2],{x1,0,1},{y1,0,1},{x2,0,1},{y2,0,1},GenerateConditions->True]//Timing

v5.0

Integrate[Sqrt[(x1-x2)^2+(y1-y2)^2],{x1,0,1},{x2,0,1},{y1,0,1},{y2,0,1}]//Timing

Integrate[Sqrt[(x1-x2)^2+(y1-y2)^2],{x1,0,1},{y1,0,1},{x2,0,1},{y2,0,1}]//Timing

Integrate[Sqrt[(x1-x2)^2+(y1-y2)^2],{x1,0,1},{y1,0,1},{x2,0,1},{y2,0,1},GenerateConditions->True]//Timing

v5.1

Integrate[Sqrt[(x1-x2)^2+(y1-y2)^2],{x1,0,1},{x2,0,1},{y1,0,1},{y2,0,1}]//Timing

Integrate[Sqrt[(x1-x2)^2+(y1-y2)^2],{x1,0,1},{y1,0,1},{x2,0,1},{y2,0,1}]//Timing

Integrate[Sqrt[(x1-x2)^2+(y1-y2)^2],{x1,0,1},{y1,0,1},{x2,0,1},{y2,0,1},GenerateConditions->True]//Timing

v6 ((Subscript[x, 1]),(Subscript[x, 2]),(Subscript[y, 1]),(Subscript[y, 2])) [CORRECT]

Integrate[Sqrt[(x1-x2)^2+(y1-y2)^2],{x1,0,1},{x2,0,1},{y1,0,1},{y2,0,1}]//Timing

(1/15) (2+(Sqrt[2])+5 ArcSinh[1])//N

NIntegrate[Sqrt[(x1-x2)^2+(y1-y2)^2],{x1,0,1},{x2,0,1},{y1,0,1},{y2,0,1}]//Timing

v6.1 ((Subscript[x, 1]),(Subscript[x, 2]),(Subscript[y, 1]),(Subscript[y, 2])) [CORRECT]

Integrate[Sqrt[(x1-x2)^2+(y1-y2)^2],{x1,0,1},{x2,0,1},{y1,0,1},{y2,0,1}]//Timing

v6.1 ((Subscript[x, 1]),(Subscript[x, 2]),(Subscript[y, 1]),(Subscript[y, 2])) [CORRECT]

Integrate[Sqrt[(x1-x2)^2+(y1-y2)^2],{x1,0,1},{x2,0,1},{y1,0,1},{y2,0,1},GenerateConditions->True]//Timing

N[(1/15) (2+(Sqrt[2])+5 ArcSinh[1])]

v6 ((Subscript[x, 1]),(Subscript[y, 1]),(Subscript[x, 2]),(Subscript[y, 2])) [WRONG]

Integrate[Sqrt[(x1-x2)^2+(y1-y2)^2],{x1,0,1},{y1,0,1},{x2,0,1},{y2,0,1}]//Timing

...wrong...

Timing[Last[%]/.Integrate[x__]:>NIntegrate[x]]

Integrate[Sqrt[(x1-x2)^2+(y1-y2)^2],{x1,0,1},{y1,0,1},{x2,0,1},{y2,0,1},GenerateConditions->True]//Timing

Integrate[Sqrt[(x1-x2)^2+(y1-y2)^2],{x1,0,1},{y1,0,1},{x2,0,1},{y2,0,1},Assumptions->And@@(0<#<1&/@{x1,x2,y1,y2})]//Timing

Mean from -1/2 to 1/2

v6 ((Subscript[x, 1]),(Subscript[y, 1]),(Subscript[x, 2]),(Subscript[y, 2])) [WRONG]

(i1=With[{h=1/2},Integrate[Sqrt[(x1-x2)^2+(y1-y2)^2],{x1,-h,h},{y1,-h,h},{x2,-h,h},{y2,-h,h}]])//Timing

N[i1]

...wrong...

v6.1 ((Subscript[x, 1]),(Subscript[y, 1]),(Subscript[x, 2]),(Subscript[y, 2])) [WRONG]

$Version

(i1=With[{h=1/2},Integrate[Sqrt[(x1-x2)^2+(y1-y2)^2],{x1,-h,h},{y1,-h,h},{x2,-h,h},{y2,-h,h}]])//Timing

N[i1]

v6.1 ((Subscript[x, 1]),(Subscript[y, 1]),(Subscript[x, 2]),(Subscript[y, 2])) WithAssumptions

(i=With[{h=1/2},IntegrateWithAssumptions[Sqrt[(x1-x2)^2+(y1-y2)^2],{x1,-h,h},{y1,-h,h},{x2,-h,h},{y2,-h,h}]])//Timing

v6 ((Subscript[x, 1]),(Subscript[y, 1]),(Subscript[x, 2]),(Subscript[y, 2])) GenerateConditions -> True [UNEVALUATED]

(i1=With[{h=1/2},Integrate[Sqrt[(x1-x2)^2+(y1-y2)^2],{x1,-h,h},{y1,-h,h},{x2,-h,h},{y2,-h,h},GenerateConditions->True]])//Timing

v6 ((Subscript[x, 1]),(Subscript[x, 2]),(Subscript[y, 1]),(Subscript[y, 2])) [UNEVALUATED]

(i1=With[{h=1/2},Integrate[Sqrt[(x1-x2)^2+(y1-y2)^2],{x1,-h,h},{x2,-h,h},{y1,-h,h},{y2,-h,h}]])//Timing

v6.1 ((Subscript[x, 1]),(Subscript[x, 2]),(Subscript[y, 1]),(Subscript[y, 2])) [UNEVALUATED]

(i1=With[{h=1/2},Integrate[Sqrt[(x1-x2)^2+(y1-y2)^2],{x1,-h,h},{x2,-h,h},{y1,-h,h},{y2,-h,h}]])//Timing

v6.1 ((Subscript[x, 1]),(Subscript[x, 2]),(Subscript[y, 1]),(Subscript[y, 2])) WithAssumptions

(i=With[{h=1/2},IntegrateWithAssumptions[Sqrt[(x1-x2)^2+(y1-y2)^2],{x1,-h,h},{x2,-h,h},{y1,-h,h},{y2,-h,h}]])//Timing

v6 ((Subscript[x, 1]),(Subscript[x, 2]),(Subscript[y, 1]),(Subscript[y, 2])) GenerateConditions -> True  [UNEVALUATED]

(i1=With[{h=1/2},Integrate[Sqrt[(x1-x2)^2+(y1-y2)^2],{x1,-h,h},{x2,-h,h},{y1,-h,h},{y2,-h,h},GenerateConditions->True]])//Timing

One at a time

(i=With[{h=1/2},IntegrateWithAssumptions[Sqrt[(x1-x2)^2+(y1-y2)^2],{x1,-h,h},{x2,-h,h},{y1,-h,h},{y2,-h,h},Simplification->FullSimplify]])//Timing

FullSimplify[%7,Abs[x1-x2]≠0]

With[{h=1/2},Integrate[Sqrt[(x1-x2)^2+(y1-y2)^2],{x1,-h,h},Assumptions->-h<x2<h&&-h<y1<h&&-h<y2<h]]

With[{h=1/2},FullSimplify[%,-h<x2<h&&-h<y1<h&&-h<y2<h]]

With[{h=1/2},Integrate[%,{x2,-h,h},Assumptions->-h<y1<h&&-h<y2<h]]

FullSimplify[(1/8)(2/3) (4 (Sqrt[1+(y1^2)-2 y1 y2+(y2^2)])+4 (y1^2) (Sqrt[1+(y1^2)-2 y1 y2+(y2^2)])-12 (((y1-y2))^2) (Sqrt[1+(y1^2)-2 y1 y2+(y2^2)])-8 y1 y2 (Sqrt[1+(y1^2)-2 y1 y2+(y2^2)])+4 (y2^2) (Sqrt[1+(y1^2)-2 y1 y2+(y2^2)])+8 (((y1-y2))^2) Abs[y1-y2]-6 (y1^2) Log[2]+18 (((y1-y2))^2) Log[2]+12 y1 y2 Log[2]-6 (y2^2) Log[2]-6 (((y1-y2))^2) Log[(((y1-y2))^2)]+6 (((y1-y2))^2) Log[(((y1-y2))^2) (((1+(Sqrt[1+(y1^2)-2 y1 y2+(y2^2)])))^2)]-6 (((y1-y2))^2) Log[-1+(Sqrt[1+4 (y1^2)-8 y1 y2+4 (y2^2)])]-6 (((y1-y2))^2) Log[1+(Sqrt[1+4 (y1^2)-8 y1 y2+4 (y2^2)])]-3 (y1^2) Log[Abs[y1-y2]]+3 (((y1-y2))^2) Log[Abs[y1-y2]]+6 y1 y2 Log[Abs[y1-y2]]-3 (y2^2) Log[Abs[y1-y2]]),-(1/2)<y1<(1/2)&&-(1/2)<y2<(1/2)]

?IntegrateWithAssumptions

Integrate[(1/6) (2 (1-2 (((y1-y2))^2)) (Sqrt[1+(((y1-y2))^2)])+4 ((Abs[y1-y2])^3)+3 (((y1-y2))^2) (2 Log[1+(Sqrt[1+(((y1-y2))^2)])]-Log[(((y1-y2))^2)])),{y1,-1/2,1/2},Assumptions->-1/2<y2<1/2]

NIntegrate[(1/288) (6+144 (y2^2)+96 (y2^4)+15 (Sqrt[5-4 y2+4 (y2^2)])-18 y2 (Sqrt[5-4 y2+4 (y2^2)])-36 (y2^2) (Sqrt[5-4 y2+4 (y2^2)])+24 (y2^3) (Sqrt[5-4 y2+4 (y2^2)])+15 (Sqrt[5+4 y2+4 (y2^2)])+18 y2 (Sqrt[5+4 y2+4 (y2^2)])-36 (y2^2) (Sqrt[5+4 y2+4 (y2^2)])-24 (y2^3) (Sqrt[5+4 y2+4 (y2^2)])+24 ArcSinh[(1/2)-y2]+24 ArcSinh[(1/2)+y2]-144 (y2^2) Log[2]+384 (y2^3) Log[2]+72 (y2^2) Log[4]-144 (y2^3) Log[4]-36 (y2^2) Log[16]+72 (y2^3) Log[16]-16 (y2^3) Log[64]-12 Log[1-2 y2]+72 y2 Log[1-2 y2]-192 (y2^3) Log[1-2 y2]-144 (y2^2) Log[(1/2)-y2]+288 (y2^3) Log[(1/2)-y2]-12 Log[1+2 y2]-72 y2 Log[1+2 y2]-144 (y2^2) Log[1+2 y2]-96 (y2^3) Log[1+2 y2]-72 y2 Log[2+(Sqrt[5+4 (-1+y2) y2])]+144 (y2^2) Log[2+(Sqrt[5+4 (-1+y2) y2])]+12 Log[2+(Sqrt[5-4 y2+4 (y2^2)])]-96 (y2^3) Log[2+(Sqrt[5-4 y2+4 (y2^2)])]+12 Log[2+(Sqrt[5+4 y2+4 (y2^2)])]+72 y2 Log[2+(Sqrt[5+4 y2+4 (y2^2)])]+144 (y2^2) Log[2+(Sqrt[5+4 y2+4 (y2^2)])]+96 (y2^3) Log[2+(Sqrt[5+4 y2+4 (y2^2)])]),{y2,-1/2,1/2}]

FullSimplify[(1/288) (6+144 (y2^2)+96 (y2^4)+15 (Sqrt[5-4 y2+4 (y2^2)])-18 y2 (Sqrt[5-4 y2+4 (y2^2)])-36 (y2^2) (Sqrt[5-4 y2+4 (y2^2)])+24 (y2^3) (Sqrt[5-4 y2+4 (y2^2)])+15 (Sqrt[5+4 y2+4 (y2^2)])+18 y2 (Sqrt[5+4 y2+4 (y2^2)])-36 (y2^2) (Sqrt[5+4 y2+4 (y2^2)])-24 (y2^3) (Sqrt[5+4 y2+4 (y2^2)])+24 ArcSinh[(1/2)-y2]+24 ArcSinh[(1/2)+y2]-144 (y2^2) Log[2]+384 (y2^3) Log[2]+72 (y2^2) Log[4]-144 (y2^3) Log[4]-36 (y2^2) Log[16]+72 (y2^3) Log[16]-16 (y2^3) Log[64]-12 Log[1-2 y2]+72 y2 Log[1-2 y2]-192 (y2^3) Log[1-2 y2]-144 (y2^2) Log[(1/2)-y2]+288 (y2^3) Log[(1/2)-y2]-12 Log[1+2 y2]-72 y2 Log[1+2 y2]-144 (y2^2) Log[1+2 y2]-96 (y2^3) Log[1+2 y2]-72 y2 Log[2+(Sqrt[5+4 (-1+y2) y2])]+144 (y2^2) Log[2+(Sqrt[5+4 (-1+y2) y2])]+12 Log[2+(Sqrt[5-4 y2+4 (y2^2)])]-96 (y2^3) Log[2+(Sqrt[5-4 y2+4 (y2^2)])]+12 Log[2+(Sqrt[5+4 y2+4 (y2^2)])]+72 y2 Log[2+(Sqrt[5+4 y2+4 (y2^2)])]+144 (y2^2) Log[2+(Sqrt[5+4 y2+4 (y2^2)])]+96 (y2^3) Log[2+(Sqrt[5+4 y2+4 (y2^2)])]),-1/2<y2<1/2]

Integrate[(1/96) (2+5 (Sqrt[5+4 (-1+y2) y2])+5 (Sqrt[5+4 y2 (1+y2)])+2 y2 (-3 (Sqrt[5+4 (-1+y2) y2])+3 (Sqrt[5+4 y2 (1+y2)])+2 y2 (2 y2 (4 y2+(Sqrt[5+4 (-1+y2) y2])-(Sqrt[5+4 y2 (1+y2)]))-3 (-4+(Sqrt[5+4 (-1+y2) y2])+(Sqrt[5+4 y2 (1+y2)]))))+8 ArcSinh[(1/2)-y2]+8 ArcSinh[(1/2)+y2]+8 y2 (3+4 (y2^2)) Log[-(((-1+2 y2) (2+(Sqrt[5+4 y2 (1+y2)])))/((1+2 y2) (2+(Sqrt[5+4 (-1+y2) y2]))))]+(4+48 (y2^2)) Log[-(((2+(Sqrt[5+4 (-1+y2) y2])) (2+(Sqrt[5+4 y2 (1+y2)])))/(-1+4 (y2^2)))]),{y2,-1/2,1/2}]

N[%]

(i=With[{h=1/2},IntegrateWithAssumptions[Sqrt[(x1-x2)^2+(y1-y2)^2],{x1,-h,h},{x2,-h,h},{y1,-h,h},{y2,-h,h},Simplification->FullSimplify]])//Timing

Points on Different Faces

Mean

TrigToExp[(1/9) (2+(Sqrt[2])+5 ArcSinh[1])]

N[%]

v5.0

2/3Integrate[Sqrt[x^2+y^2],{x,0,1},{y,0,1}]+1/3Integrate[Sqrt[1+(y-u)^2],{u,0,1},{y,0,1}]//FullSimplify//Timing

v5.1

(2/3)Integrate[Sqrt[x^2+y^2],{x,0,1},{y,0,1}]+(1/3)Integrate[Sqrt[1+(y-u)^2],{u,0,1},{y,0,1}]//FullSimplify//Timing

v6

(2/3)Integrate[Sqrt[x^2+y^2],{x,0,1},{y,0,1}]+(1/3)Integrate[Sqrt[1+(y-u)^2],{u,0,1},{y,0,1}]//FullSimplify//Timing

Numerical

FacePoint[1]:={0,Random[Real]}
FacePoint[2]:={1,Random[Real]}
FacePoint[3]:={Random[Real],0}
FacePoint[4]:={Random[Real],1}

<<DiscreteMath`Combinatorica`

f=KSubsets[Range[4],2]

Mean[Norm/@Subtract@@@Map[FacePoint,f[[#]]&/@Table[Random[Integer,{1,6}],{10^6}],{2}]]//Timing