// forked from makc3d's Parabola through 3 points
// forked from makc3d's Solving a system of linear equations
package {
import flash.display.Sprite;
import flash.events.Event;
import flash.geom.Point;
import flash.geom.Rectangle;
/**
* Parabola through three points.
* @author makc
*/
public class Test1 extends Sprite {
public var a1:Anchor;
public var a2:Anchor;
public var a3:Anchor;
public var solver:LinearSolver;
public function Test1 () {
addChild (a1 = new Anchor (r, r));
addChild (a2 = new Anchor (r, r));
addChild (a3 = new Anchor (r, r));
solver = new LinearSolver;
addEventListener (Event.ENTER_FRAME, findParabola);
}
public function get r ():int {
return 10 + Math.random () * 445;
}
public function findParabola (whatever:*= null):void {
var x1:Number = a1.x;
var y1:Number = a1.y;
var x2:Number = a2.x;
var y2:Number = a2.y;
var x3:Number = a3.x;
var y3:Number = a3.y;
solver.prepare (3);
// we have three equations A * x² + B * x + C = y
solver.A [0] [0] = x1 * x1; solver.A [0] [1] = x1; solver.A [0] [2] = 1;
solver.A [1] [0] = x2 * x2; solver.A [1] [1] = x2; solver.A [1] [2] = 1;
solver.A [2] [0] = x3 * x3; solver.A [2] [1] = x3; solver.A [2] [2] = 1;
solver.b [0] = y1;
solver.b [1] = y2;
solver.b [2] = y3;
solver.solve ();
graphics.clear ();
graphics.lineStyle (1, 255 * 256 * 256);
for (var i:int = 0; i < 465; i += 3) {
graphics [(i < 1) ? "moveTo" : "lineTo"] (i,
solver.x [0] * i * i + solver.x [1] * i + solver.x [2]
);
}
}
}
}
import flash.display.Sprite;
import flash.events.MouseEvent;
class Anchor extends Sprite {
public function Anchor (x0:int, y0:int) {
x = x0; y = y0;
graphics.beginFill (0xFF7F00, 1);
graphics.drawRect (-4, -4, 8, 8);
graphics.drawRect (-2, -2, 4, 4);
graphics.beginFill (0xFF7F00, 0);
graphics.drawRect (-2, -2, 4, 4);
useHandCursor = buttonMode = true;
addEventListener (MouseEvent.MOUSE_DOWN, startDragMe);
addEventListener (MouseEvent.MOUSE_UP, stopDragMe);
}
public function startDragMe (e:MouseEvent):void {
startDrag ();
}
public function stopDragMe (e:MouseEvent):void {
stopDrag ();
}
}
/**
* Analytically solves the linear system Ax = b.
*
* <p>A is first decomposed into a lower triangular matrix, L and an
* upper triangular matrix U:
*
* <pre>
* Ax = b
* LUx = b
*
* Then, using forward substitution it solves the equation:
*
* Ly = b
*
* Finally, taking the solution to this system (y), we can solve for x
* with back substitution:
*
* Ux = y
* </pre>
*
* @author Drew Cummins, http://lab.generalrelativity.org/file/MatrixND.as
* @author makc, distilled the solver, hacked L[i][i] to be != 0
*/
class LinearSolver {
public var n:int = -1;
public var A:Vector.<Vector.<Number>>;
public var L:Vector.<Vector.<Number>>;
public var U:Vector.<Vector.<Number>>;
public var b:Vector.<Number>;
public var x:Vector.<Number>;
public function prepare (numOfVars:int):void {
if (n != numOfVars) {
n = numOfVars;
A = new Vector.<Vector.<Number>> (n, true);
L = new Vector.<Vector.<Number>> (n, true);
U = new Vector.<Vector.<Number>> (n, true);
b = new Vector.<Number> (n, true); // zeroes
x = new Vector.<Number> (n, true); // zeroes
} else {
// zero b, x
for (var i:int = 0; i < n; i++) {
b [i] = 0; x [i] = 0;
}
}
for (i = 0; i < n; i++) {
var Ai:Vector.<Number> = A [i];
if (Ai == null) {
A [i] = new Vector.<Number> (n, true); // zeroes
L [i] = new Vector.<Number> (n, true); // zeroes
U [i] = new Vector.<Number> (n, true); U [i] [i] = 1; // identity in U
} else {
var Li:Vector.<Number> = L [i];
var Ui:Vector.<Number> = U [i];
for (var j:int = 0; j < n; j++) {
// zero matrices A, L, identity in U
Ai [j] = 0; Li [j] = 0; Ui [j] = (i == j) ? 1 : 0;
}
}
}
}
public function solve ():void {
var dot:Number;
// non-pivoting Crout algorithm to decompose A into LU
for (var k:int = 0; k < n; k++) {
for (var i:int = k; i < n; i++) {
dot = 0;
var Li:Vector.<Number> = L [i];
for (var p:int = 0; p < k; p++) {
dot += Li [p] * U [p] [k];
}
Li [k] = A [i] [k] - dot;
// a hack to work around non-pivoting algorithm limitation:
// http://mymathlib.webtrellis.net/matrices/linearsystems/crout.html
// "If the matrix A is positive definite symmetric or if the matrix is
// diagonally dominant, then pivoting is not necessary; otherwise the
// version using pivoting should be used"
if ((i == k) && (Math.abs (Li [k]) < 1e-13)) Li [k] = 1e-13;
}
var Ak:Vector.<Number> = A [k];
var Lk:Vector.<Number> = L [k];
var Uk:Vector.<Number> = U [k];
for (var j:int = k + 1; j < n; j++) {
dot = 0;
for (p = 0; p < k; p++) {
dot += Lk [p] * U [p] [j];
}
Uk [j] = ( Ak [j] - dot ) / Lk [k];
}
}
// forward substitution to solve Ly = b
// (x is used to hold y and x)
for (i = 0; i < n; i++) {
dot = 0;
Li = L [i];
for (j = 0; j < i; j++) {
dot += Li [j] * x [j];
}
x [i] = ( b [i] - dot ) / Li [i];
}
// backward substitution to solve for Ux = y
for (i = n - 1; i >= 0; i--) {
dot = 0;
var Ui:Vector.<Number> = U [i];
for (j = i + 1; j < n; j++) {
dot += Ui [j] * x [j];
}
// y vector only accessed before solution is set at that index: (x [i] - dot)
// this is why we don't need another structure to hold y
x [i] = ( x [i] - dot ) / Ui [i];
}
}
}