This question of determining whether points form a circle or not can be resolved by solving two sub-problems:

1. Express the points with just two coordinates rather than $n=3$.

2. Test whether a collection of points in a plane is concyclic.

We might as well solve the first sub-problem for dimensions $n\ge 3$, because it all works the same: we find an orthonormal basis of $n$-space where the first two basis elements (appropriately translated) span the plane containing the points. Expressing the points in that basis gives the desired 2D coordinates.

An elegant way to identify concyclic planar points is to view them as complex numbers and check them in groups of four each to see whether they have a real cross ratio: if so, each group lies on a common generalized circle (that is, they are either collinear or concyclic). We may pick any three distinct points and, adjoining each of the remaining points in turn, check each of these groups of four. If they are all found to have real cross-ratios, then they all must lie on the (unique) generalized circle determined by those three distinct points.

This plan is executed by concyclicQ, which includes the checks for degenerate situations (such as all points collinear). Some care is needed in the use of Orthogonalize: just a little bit of noise in the data orthogonal to the plane can cause it to conclude the points are not coplanar. To force the points to be considered coplanar, just remove the check If[Length[f] > 2.... (Notice the uses of Chop to reduce the effects of floating point imprecision.)

crossRatio[z1_, z2_, z3_, z4_] := (z1 - z3) (z2 - z4) / ((z2 - z3) (z1 - z4));
concyclicQ[pts_List] /; Length[pts] >= 3 := Module[{points, k, f, g, q, p, test, tol=10^-6},
   (* Step 1: project into the Complex plane *)
   points = Union[pts];                             (* Eliminate any duplicates *)
   q = # - Mean[points] & /@ points;                (* Center the points *)
   f = Select[Chop[Orthogonalize[q, Tolerance -> tol]], Norm[#] > 0 &]; (* Find a basis *)
   If[Length[f] == 0, Return[True]];                (* All are coincident *)
   If[Length[f] == 1, Return[Length[points] <= 2]]; (* Collinear *)
   If[Length[f] > 2, Return[False]];                (* Not coplanar *)
   f = Select[Chop[Orthogonalize[f~Join~IdentityMatrix[Dimensions[points][[2]]]]], 
     Norm[#] > 0 &];                                (* Find an adapted basis *)
   g = ( Inverse[f] )[[All, 1 ;; 2]];               (* Change-of-basis matrix *)
   p = Complex @@@ (q . g);                         (* Change the basis *)

   (* Step 2: test for concyclicality *)
   test = Append[p[[1 ;; 3]], #] & /@ p[[4 ;;]];    (* Test last n-3 against the first 3 *)
   Min[Boole[# \[Element] Reals] & /@ Chop[crossRatio @@@ test, tol]] == 1
   ];
concyclicQ[points_List] := True (* Three points or fewer *)

This code works in greater than three dimensions, too.

Examples:

With[{n = 600, e = Orthogonalize[RandomReal[NormalDistribution[0, 1], {3, 3}]]},
 points = Join[#, ConstantArray[1,  3 - 2]] & /@ (Through[{Cos, Sin}[ 
         RandomReal[{0, 2 \[Pi]}, n]]]\[Transpose]) . e;
]; (* Many random concyclic points *)
concyclicQ[points]

True

q = Append[points, {0, 0, 0}]; (* Throw in a "bad" point *)
concyclicQ[q]

False

concyclicQ[RandomReal[{0, 1}, {4, 3}]] (* Four points in general position *)

False

concyclicQ[{RandomReal[{0, 1}, 4]}\[Transpose] . {{1, 1, 1}}] (* Collinear *)

False

concyclicQ[ConstantArray[ {0, 0, 0}, 4]] (* Coincident *)

True