PARI/GP script to find the center of a circle, given 3 points on its circumference

[randomfox]

PARI is a high-precision math library and computer algebra system. GP is a scripting interface to that library.

Coordinates should be converted to UTM for use with the script.

To run a script from GP, enter \r filename on the command line. If the script fails to run because some of the symbols already exist, enter kill(circle), for example, before running the script, or start a fresh GP session.

dist(a, b) = sqrt(sqr(a[1]-b[1]) + sqr(a[2] - b[2]));

circle(b, c, d) = 
{
    local(temp, bc, cd, det, circ);

    temp = sqr(c[1]) + sqr(c[2]);
    bc = (sqr(b[1]) + sqr(b[2]) - temp) / 2.;
    cd = (temp - sqr(d[1]) - sqr(d[2])) / 2.;
    det = (b[1]-c[1]) * (c[2]-d[2]) - (c[1]-d[1]) * (b[2]-c[2]);

    circ = [
	(bc * (c[2]-d[2]) - cd * (b[2]-c[2])) / det,
	((b[1]-c[1]) * cd - (c[1]-d[1]) * bc) / det
    ];

    [ circ, [ dist(circ, b), dist(circ, c), dist(circ, d) ] ];
}

b = [515939, 4413282];
c = [519054, 4405649];
d = [516437, 4399896];

circle(b, c, d)