[Ether]

Greg, try this for the internal points problem: taking advantage of symmetry, generate 2 random points in Quadrant1, then reflect one of the points into Quadrants 2 thru 4 so you get 4 lengths per iteration. Then you can use 4 times fewer iterations. Also, doing this allows you to significantly optimize the code to eliminate a lot of repetitious floating point operations. I think you'll find that approach to be much faster. I was able to get 1 ppm accuracy in 7 minutes of runtime on a 10-year-old Pentium D desktop machine (using 32 bit compiled code). Here's the pseudocode for it:

pio2=pi/2; // initialize constant

// iterate:
r1=random; r2=random; dq=pio2*(random-random); 
tmp1=r1+r2; tmp2=2*sqrt(r1*r2);
cosdq=cos(dq); sindq=sin(dq);
sum +=
  sqrt(tmp1-tmp2*cosdq)	
+ sqrt(tmp1+tmp2*sindq)	
+ sqrt(tmp1+tmp2*cosdq)	
+ sqrt(tmp1-tmp2*sindq);	

// don't forget to divide by 4*iterations when done iterating

Using symmetry for the chords problem, you can reduce the length computation to:

sum += sin(pi*abs(random-random));

... and when the iterations are complete, multiply the sum by 2.