4dsphere.calc

Calc Software Package for Number Calc

#!/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))