#!/usr/local/src/cmd/calc/calc -q -f/* * 4dsphere - determine if 6 points lie on the surface of a sphere in R^4 * * usage: *      4dsphere x0 y0 z0 w0 x1 y1 z1 w1 ... x5 y5 z5 w5 * *      x0 y0 z0 w0	point 0 in R^4 *      x1 y1 z1 w1	point 1 in R^4 *	...		... *      x5 y5 z5 w5	point 5 in R^4 * * Copyright (C) 2001  Landon Curt Noll * * Calc is open software; you can redistribute it and/or modify it under * the terms of the version 2.1 of the GNU Lesser General Public License * as published by the Free Software Foundation. * * Calc is distributed in the hope that it will be useful, but WITHOUT * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY * or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General * Public License for more details. * * A copy of version 2.1 of the GNU Lesser General Public License is * distributed with calc under the filename COPYING-LGPL.  You should have * received a copy with calc; if not, write to Free Software Foundation, Inc. * 59 Temple Place, Suite 330, Boston, MA  02111-1307, USA. * * @(#) $Revision: 1.4 $ * @(#) $Id: 4dsphere.calc,v 1.4 2007/02/07 00:34:04 chongo Exp $ * @(#) $Source: /usr/local/src/cmd/calc/cscript/RCS/4dsphere.calc,v $ * * Under source code control:   2001/05/03 19:02:03 * File existed as early as:    2001 * * chongo <was here> /\oo/\     http://www.isthe.com/chongo/ * Share and enjoy!  :-)        http://www.isthe.com/chongo/tech/comp/calc/ *//* * parse args */argc = argv();if (argc != 25) {    fprintf(files(2), "usage: %s x0 y0 z0 w0 x1 y1 z1 w1 ... x5 y5 z5 w5\n",	argv(0));    exit;}x0 = eval(argv(1));y0 = eval(argv(2));z0 = eval(argv(3));w0 = eval(argv(4));x1 = eval(argv(5));y1 = eval(argv(6));z1 = eval(argv(7));w1 = eval(argv(8));x2 = eval(argv(9));y2 = eval(argv(10));z2 = eval(argv(11));w2 = eval(argv(12));x3 = eval(argv(13));y3 = eval(argv(14));z3 = eval(argv(15));w3 = eval(argv(16));x4 = eval(argv(17));y4 = eval(argv(18));z4 = eval(argv(19));w4 = eval(argv(20));x5 = eval(argv(21));y5 = eval(argv(22));z5 = eval(argv(23));w5 = eval(argv(24));/* * verbose output setup */print "(":x0:",":y0:",":z0:",":w0:") ":;print "(":x1:",":y1:",":z1:",":w1:") ":;print "(":x2:",":y2:",":z2:",":w2:") ":;print "(":x3:",":y3:",":z3:",":w3:") ":;print "(":x4:",":y4:",":z4:",":w4:") ":;print "(":x5:",":y5:",":z5:",":w5:") ":;/* * * Given the 5 points: * * 	(x0,y1,z1,w1) * 	(x1,y1,z1,w1) * 	(x2,y2,z2,w2) * 	(x3,y3,z3,w3) * 	(x4,y4,z4,w4) * 	(x5,y5,z5,w5) * * we can determine if they lie in the surface of 4D sphere in R^4 if the * following matrix is 0: * * 	| x0^2+y0^2+z0^2+w0^2	x0	y0	z0	w0	1 | * 	| x1^2+y1^2+z1^2+w1^2	x1	y1	z1	w1	1 | * 	| x2^2+y2^2+z2^2+w2^2	x2	y2	z2	w2	1 |	= 0 * 	| x3^2+y3^2+z3^2+w3^2	x3	y3	z3	w3	1 | * 	| x4^2+y4^2+z4^2+w4^2	x4	y4	z4	w4	1 | * 	| x5^2+y5^2+z5^2+w5^2	x5	y5	z5	w5	1 | */if ((w0*(-x1*(-y2*(-z4*(z5^2+y5^2+x5^2+w5^2)    -z3*(-z5^2+z4^2-y5^2+y4^2-x5^2+x4^2-w5^2+w4^2) +    (z4^2+y4^2+x4^2+w4^2)*z5+ (z3^2+y3^2+x3^2+w3^2)*(z4-z5))    +y3*((z4^2+y4^2+x4^2+w4^2)*z5-z4*(z5^2+y5^2+x5^2+w5^2))    +z2*(-y4*(z5^2+y5^2+x5^2+w5^2)    -y3*(-z5^2+z4^2-y5^2+y4^2-x5^2+x4^2-w5^2+w4^2)    +y5*(z4^2+y4^2+x4^2+w4^2) + (y4-y5)*(z3^2+y3^2+x3^2+w3^2))    -z3*(y5*(z4^2+y4^2+x4^2+w4^2) -y4*(z5^2+y5^2+x5^2+w5^2)) +    (z2^2+y2^2+x2^2+w2^2)*(y4*z5+y3*(z4-z5) -y5*z4- (y4-y5)*z3) -    (z3^2+y3^2+x3^2+w3^2)*(y4*z5-y5*z4)) +y1*(-x2*(-z4*(z5^2+y5^2+x5^2+w5^2)    -z3*(-z5^2+z4^2-y5^2+y4^2-x5^2+x4^2-w5^2+w4^2) +    (z4^2+y4^2+x4^2+w4^2)*z5+ (z3^2+y3^2+x3^2+w3^2)*(z4-z5))    +x3*((z4^2+y4^2+x4^2+w4^2)*z5-z4*(z5^2+y5^2+x5^2+w5^2))    +z2*(-x4*(z5^2+y5^2+x5^2+w5^2)    -x3*(-z5^2+z4^2-y5^2+y4^2-x5^2+x4^2-w5^2+w4^2) +x5*(z4^2+y4^2+x4^2+w4^2)    + (x4-x5)*(z3^2+y3^2+x3^2+w3^2)) -z3*(x5*(z4^2+y4^2+x4^2+w4^2)    -x4*(z5^2+y5^2+x5^2+w5^2)) + (z2^2+y2^2+x2^2+w2^2)*(x4*z5+x3*(z4-z5)    -x5*z4- (x4-x5)*z3) - (z3^2+y3^2+x3^2+w3^2)*(x4*z5-x5*z4))    -x2*(-y3*((z4^2+y4^2+x4^2+w4^2)*z5-z4*(z5^2+y5^2+x5^2+w5^2))    +z3*(y5*(z4^2+y4^2+x4^2+w4^2) -y4*(z5^2+y5^2+x5^2+w5^2))    + (z3^2+y3^2+x3^2+w3^2)*(y4*z5-y5*z4))    +y2*(-x3*((z4^2+y4^2+x4^2+w4^2)*z5-z4*(z5^2+y5^2+x5^2+w5^2))    +z3*(x5*(z4^2+y4^2+x4^2+w4^2) -x4*(z5^2+y5^2+x5^2+w5^2)) +    (z3^2+y3^2+x3^2+w3^2)*(x4*z5-x5*z4)) -z1*(-x2*(-y4*(z5^2+y5^2+x5^2+w5^2)    -y3*(-z5^2+z4^2-y5^2+y4^2-x5^2+x4^2-w5^2+w4^2) +y5*(z4^2+y4^2+x4^2+w4^2)    + (y4-y5)*(z3^2+y3^2+x3^2+w3^2)) +x3*(y5*(z4^2+y4^2+x4^2+w4^2)    -y4*(z5^2+y5^2+x5^2+w5^2)) +y2*(-x4*(z5^2+y5^2+x5^2+w5^2)    -x3*(-z5^2+z4^2-y5^2+y4^2-x5^2+x4^2-w5^2+w4^2)    +x5*(z4^2+y4^2+x4^2+w4^2) + (x4-x5)*(z3^2+y3^2+x3^2+w3^2))    -y3*(x5*(z4^2+y4^2+x4^2+w4^2) -x4*(z5^2+y5^2+x5^2+w5^2)) -    (x4*y5-x5*y4)*(z3^2+y3^2+x3^2+w3^2) + (x4*y5+x3*(y4-y5) -x5*y4-    (x4-x5)*y3)*(z2^2+y2^2+x2^2+w2^2)) -z2*(-x3*(y5*(z4^2+y4^2+x4^2+w4^2)    -y4*(z5^2+y5^2+x5^2+w5^2)) +y3*(x5*(z4^2+y4^2+x4^2+w4^2)    -x4*(z5^2+y5^2+x5^2+w5^2)) + (x4*y5-x5*y4)*(z3^2+y3^2+x3^2+w3^2))    + (z1^2+y1^2+x1^2+w1^2)*(x2*(y4*z5+y3*(z4-z5) -y5*z4-    (y4-y5)*z3) -x3*(y4*z5-y5*z4) -y2*(x4*z5+x3*(z4-z5) -x5*z4-    (x4-x5)*z3) +y3*(x4*z5-x5*z4) - (x4*y5-x5*y4)*z3+    (x4*y5+x3*(y4-y5) -x5*y4- (x4-x5)*y3)*z2) +    (z2^2+y2^2+x2^2+w2^2)*(x3*(y4*z5-y5*z4) -y3*(x4*z5-x5*z4) +    (x4*y5-x5*y4)*z3)) -x0*(-w1*(-y2*(-z4*(z5^2+y5^2+x5^2+w5^2)    -z3*(-z5^2+z4^2-y5^2+y4^2-x5^2+x4^2-w5^2+w4^2) +    (z4^2+y4^2+x4^2+w4^2)*z5+ (z3^2+y3^2+x3^2+w3^2)*(z4-z5))    +y3*((z4^2+y4^2+x4^2+w4^2)*z5-z4*(z5^2+y5^2+x5^2+w5^2))    +z2*(-y4*(z5^2+y5^2+x5^2+w5^2)    -y3*(-z5^2+z4^2-y5^2+y4^2-x5^2+x4^2-w5^2+w4^2)    +y5*(z4^2+y4^2+x4^2+w4^2) + (y4-y5)*(z3^2+y3^2+x3^2+w3^2))    -z3*(y5*(z4^2+y4^2+x4^2+w4^2) -y4*(z5^2+y5^2+x5^2+w5^2)) +    (z2^2+y2^2+x2^2+w2^2)*(y4*z5+y3*(z4-z5) -y5*z4- (y4-y5)*z3) -    (z3^2+y3^2+x3^2+w3^2)*(y4*z5-y5*z4)) -y1*(w2*(-z4*(z5^2+y5^2+x5^2+w5^2)    -z3*(-z5^2+z4^2-y5^2+y4^2-x5^2+x4^2-w5^2+w4^2) +    (z4^2+y4^2+x4^2+w4^2)*z5+ (z3^2+y3^2+x3^2+w3^2)*(z4-z5))    -w3*((z4^2+y4^2+x4^2+w4^2)*z5-z4*(z5^2+y5^2+x5^2+w5^2))    -z2*(-w4*(z5^2+y5^2+x5^2+w5^2)    -w3*(-z5^2+z4^2-y5^2+y4^2-x5^2+x4^2-w5^2+w4^2) +w5*(z4^2+y4^2+x4^2+w4^2)    + (w4-w5)*(z3^2+y3^2+x3^2+w3^2)) +z3*(w5*(z4^2+y4^2+x4^2+w4^2)    -w4*(z5^2+y5^2+x5^2+w5^2)) + (z2^2+y2^2+x2^2+w2^2)*(-w4*z5-w3*(z4-z5)    +w5*z4+ (w4-w5)*z3) - (z3^2+y3^2+x3^2+w3^2)*(w5*z4-w4*z5))