chore: Unify math types, utils and functions (#8389)

Co-authored-by: dwelle <5153846+dwelle@users.noreply.github.com>
This commit is contained in:
Márk Tolmács 2024-09-03 00:23:38 +02:00 committed by GitHub
parent e3d1dee9d0
commit f4dd23fc31
No known key found for this signature in database
GPG key ID: B5690EEEBB952194
98 changed files with 4291 additions and 3661 deletions

View file

@ -0,0 +1,70 @@
import * as GA from "./ga";
import { point, toString, direction, offset } from "./ga";
import * as GAPoint from "./gapoints";
import * as GALine from "./galines";
import * as GATransform from "./gatransforms";
describe("geometric algebra", () => {
describe("points", () => {
it("distanceToLine", () => {
const point = GA.point(3, 3);
const line = GALine.equation(0, 1, -1);
expect(GAPoint.distanceToLine(point, line)).toEqual(2);
});
it("distanceToLine neg", () => {
const point = GA.point(-3, -3);
const line = GALine.equation(0, 1, -1);
expect(GAPoint.distanceToLine(point, line)).toEqual(-4);
});
});
describe("lines", () => {
it("through", () => {
const a = GA.point(0, 0);
const b = GA.point(2, 0);
expect(toString(GALine.through(a, b))).toEqual(
toString(GALine.equation(0, 2, 0)),
);
});
it("parallel", () => {
const point = GA.point(3, 3);
const line = GALine.equation(0, 1, -1);
const parallel = GALine.parallel(line, 2);
expect(GAPoint.distanceToLine(point, parallel)).toEqual(0);
});
});
describe("translation", () => {
it("points", () => {
const start = point(2, 2);
const move = GATransform.translation(direction(0, 1));
const end = GATransform.apply(move, start);
expect(toString(end)).toEqual(toString(point(2, 3)));
});
it("points 2", () => {
const start = point(2, 2);
const move = GATransform.translation(offset(3, 4));
const end = GATransform.apply(move, start);
expect(toString(end)).toEqual(toString(point(5, 6)));
});
it("lines", () => {
const original = GALine.through(point(2, 2), point(3, 4));
const move = GATransform.translation(offset(3, 4));
const parallel = GATransform.apply(move, original);
expect(toString(parallel)).toEqual(
toString(GALine.through(point(5, 6), point(6, 8))),
);
});
});
describe("rotation", () => {
it("points", () => {
const start = point(2, 2);
const pivot = point(1, 1);
const rotate = GATransform.rotation(pivot, Math.PI / 2);
const end = GATransform.apply(rotate, start);
expect(toString(end)).toEqual(toString(point(2, 0)));
});
});
});

317
packages/math/ga/ga.ts Normal file
View file

@ -0,0 +1,317 @@
/**
* This is a 2D Projective Geometric Algebra implementation.
*
* For wider context on geometric algebra visit see https://bivector.net.
*
* For this specific algebra see cheatsheet https://bivector.net/2DPGA.pdf.
*
* Converted from generator written by enki, with a ton of added on top.
*
* This library uses 8-vectors to represent points, directions and lines
* in 2D space.
*
* An array `[a, b, c, d, e, f, g, h]` represents a n(8)vector:
* a + b*e0 + c*e1 + d*e2 + e*e01 + f*e20 + g*e12 + h*e012
*
* See GAPoint, GALine, GADirection and GATransform modules for common
* operations.
*/
export type Point = NVector;
export type Direction = NVector;
export type Line = NVector;
export type Transform = NVector;
export const point = (x: number, y: number): Point => [0, 0, 0, 0, y, x, 1, 0];
export const origin = (): Point => [0, 0, 0, 0, 0, 0, 1, 0];
export const direction = (x: number, y: number): Direction => {
const norm = Math.hypot(x, y); // same as `inorm(direction(x, y))`
return [0, 0, 0, 0, y / norm, x / norm, 0, 0];
};
export const offset = (x: number, y: number): Direction => [
0,
0,
0,
0,
y,
x,
0,
0,
];
/// This is the "implementation" part of the library
type NVector = readonly [
number,
number,
number,
number,
number,
number,
number,
number,
];
// These are labels for what each number in an nvector represents
const NVECTOR_BASE = ["1", "e0", "e1", "e2", "e01", "e20", "e12", "e012"];
// Used to represent points, lines and transformations
export const nvector = (value: number = 0, index: number = 0): NVector => {
const result = [0, 0, 0, 0, 0, 0, 0, 0];
if (index < 0 || index > 7) {
throw new Error(`Expected \`index\` between 0 and 7, got \`${index}\``);
}
if (value !== 0) {
result[index] = value;
}
return result as unknown as NVector;
};
const STRING_EPSILON = 0.000001;
export const toString = (nvector: NVector): string => {
const result = nvector
.map((value, index) =>
Math.abs(value) > STRING_EPSILON
? value.toFixed(7).replace(/(\.|0+)$/, "") +
(index > 0 ? NVECTOR_BASE[index] : "")
: null,
)
.filter((representation) => representation != null)
.join(" + ");
return result === "" ? "0" : result;
};
// Reverse the order of the basis blades.
export const reverse = (nvector: NVector): NVector => [
nvector[0],
nvector[1],
nvector[2],
nvector[3],
-nvector[4],
-nvector[5],
-nvector[6],
-nvector[7],
];
// Poincare duality operator.
export const dual = (nvector: NVector): NVector => [
nvector[7],
nvector[6],
nvector[5],
nvector[4],
nvector[3],
nvector[2],
nvector[1],
nvector[0],
];
// Clifford Conjugation
export const conjugate = (nvector: NVector): NVector => [
nvector[0],
-nvector[1],
-nvector[2],
-nvector[3],
-nvector[4],
-nvector[5],
-nvector[6],
nvector[7],
];
// Main involution
export const involute = (nvector: NVector): NVector => [
nvector[0],
-nvector[1],
-nvector[2],
-nvector[3],
nvector[4],
nvector[5],
nvector[6],
-nvector[7],
];
// Multivector addition
export const add = (a: NVector, b: NVector | number): NVector => {
if (isNumber(b)) {
return [a[0] + b, a[1], a[2], a[3], a[4], a[5], a[6], a[7]];
}
return [
a[0] + b[0],
a[1] + b[1],
a[2] + b[2],
a[3] + b[3],
a[4] + b[4],
a[5] + b[5],
a[6] + b[6],
a[7] + b[7],
];
};
// Multivector subtraction
export const sub = (a: NVector, b: NVector | number): NVector => {
if (isNumber(b)) {
return [a[0] - b, a[1], a[2], a[3], a[4], a[5], a[6], a[7]];
}
return [
a[0] - b[0],
a[1] - b[1],
a[2] - b[2],
a[3] - b[3],
a[4] - b[4],
a[5] - b[5],
a[6] - b[6],
a[7] - b[7],
];
};
// The geometric product.
export const mul = (a: NVector, b: NVector | number): NVector => {
if (isNumber(b)) {
return [
a[0] * b,
a[1] * b,
a[2] * b,
a[3] * b,
a[4] * b,
a[5] * b,
a[6] * b,
a[7] * b,
];
}
return [
mulScalar(a, b),
b[1] * a[0] +
b[0] * a[1] -
b[4] * a[2] +
b[5] * a[3] +
b[2] * a[4] -
b[3] * a[5] -
b[7] * a[6] -
b[6] * a[7],
b[2] * a[0] + b[0] * a[2] - b[6] * a[3] + b[3] * a[6],
b[3] * a[0] + b[6] * a[2] + b[0] * a[3] - b[2] * a[6],
b[4] * a[0] +
b[2] * a[1] -
b[1] * a[2] +
b[7] * a[3] +
b[0] * a[4] +
b[6] * a[5] -
b[5] * a[6] +
b[3] * a[7],
b[5] * a[0] -
b[3] * a[1] +
b[7] * a[2] +
b[1] * a[3] -
b[6] * a[4] +
b[0] * a[5] +
b[4] * a[6] +
b[2] * a[7],
b[6] * a[0] + b[3] * a[2] - b[2] * a[3] + b[0] * a[6],
b[7] * a[0] +
b[6] * a[1] +
b[5] * a[2] +
b[4] * a[3] +
b[3] * a[4] +
b[2] * a[5] +
b[1] * a[6] +
b[0] * a[7],
];
};
export const mulScalar = (a: NVector, b: NVector): number =>
b[0] * a[0] + b[2] * a[2] + b[3] * a[3] - b[6] * a[6];
// The outer/exterior/wedge product.
export const meet = (a: NVector, b: NVector): NVector => [
b[0] * a[0],
b[1] * a[0] + b[0] * a[1],
b[2] * a[0] + b[0] * a[2],
b[3] * a[0] + b[0] * a[3],
b[4] * a[0] + b[2] * a[1] - b[1] * a[2] + b[0] * a[4],
b[5] * a[0] - b[3] * a[1] + b[1] * a[3] + b[0] * a[5],
b[6] * a[0] + b[3] * a[2] - b[2] * a[3] + b[0] * a[6],
b[7] * a[0] +
b[6] * a[1] +
b[5] * a[2] +
b[4] * a[3] +
b[3] * a[4] +
b[2] * a[5] +
b[1] * a[6],
];
// The regressive product.
export const join = (a: NVector, b: NVector): NVector => [
joinScalar(a, b),
a[1] * b[7] + a[4] * b[5] - a[5] * b[4] + a[7] * b[1],
a[2] * b[7] - a[4] * b[6] + a[6] * b[4] + a[7] * b[2],
a[3] * b[7] + a[5] * b[6] - a[6] * b[5] + a[7] * b[3],
a[4] * b[7] + a[7] * b[4],
a[5] * b[7] + a[7] * b[5],
a[6] * b[7] + a[7] * b[6],
a[7] * b[7],
];
export const joinScalar = (a: NVector, b: NVector): number =>
a[0] * b[7] +
a[1] * b[6] +
a[2] * b[5] +
a[3] * b[4] +
a[4] * b[3] +
a[5] * b[2] +
a[6] * b[1] +
a[7] * b[0];
// The inner product.
export const dot = (a: NVector, b: NVector): NVector => [
b[0] * a[0] + b[2] * a[2] + b[3] * a[3] - b[6] * a[6],
b[1] * a[0] +
b[0] * a[1] -
b[4] * a[2] +
b[5] * a[3] +
b[2] * a[4] -
b[3] * a[5] -
b[7] * a[6] -
b[6] * a[7],
b[2] * a[0] + b[0] * a[2] - b[6] * a[3] + b[3] * a[6],
b[3] * a[0] + b[6] * a[2] + b[0] * a[3] - b[2] * a[6],
b[4] * a[0] + b[7] * a[3] + b[0] * a[4] + b[3] * a[7],
b[5] * a[0] + b[7] * a[2] + b[0] * a[5] + b[2] * a[7],
b[6] * a[0] + b[0] * a[6],
b[7] * a[0] + b[0] * a[7],
];
export const norm = (a: NVector): number =>
Math.sqrt(Math.abs(a[0] * a[0] - a[2] * a[2] - a[3] * a[3] + a[6] * a[6]));
export const inorm = (a: NVector): number =>
Math.sqrt(Math.abs(a[7] * a[7] - a[5] * a[5] - a[4] * a[4] + a[1] * a[1]));
export const normalized = (a: NVector): NVector => {
const n = norm(a);
if (n === 0 || n === 1) {
return a;
}
const sign = a[6] < 0 ? -1 : 1;
return mul(a, sign / n);
};
export const inormalized = (a: NVector): NVector => {
const n = inorm(a);
if (n === 0 || n === 1) {
return a;
}
return mul(a, 1 / n);
};
const isNumber = (a: any): a is number => typeof a === "number";
export const E0: NVector = nvector(1, 1);
export const E1: NVector = nvector(1, 2);
export const E2: NVector = nvector(1, 3);
export const E01: NVector = nvector(1, 4);
export const E20: NVector = nvector(1, 5);
export const E12: NVector = nvector(1, 6);
export const E012: NVector = nvector(1, 7);
export const I = E012;

View file

@ -0,0 +1,26 @@
import * as GA from "./ga";
import type { Line, Direction, Point } from "./ga";
/**
* A direction is stored as an array `[0, 0, 0, 0, y, x, 0, 0]` representing
* vector `(x, y)`.
*/
export const from = (point: Point): Point => [
0,
0,
0,
0,
point[4],
point[5],
0,
0,
];
export const fromTo = (from: Point, to: Point): Direction =>
GA.inormalized([0, 0, 0, 0, to[4] - from[4], to[5] - from[5], 0, 0]);
export const orthogonal = (direction: Direction): Direction =>
GA.inormalized([0, 0, 0, 0, -direction[5], direction[4], 0, 0]);
export const orthogonalToLine = (line: Line): Direction => GA.mul(line, GA.I);

View file

@ -0,0 +1,52 @@
import * as GA from "./ga";
import type { Line, Point } from "./ga";
/**
* A line is stored as an array `[0, c, a, b, 0, 0, 0, 0]` representing:
* c * e0 + a * e1 + b*e2
*
* This maps to a standard formula `a * x + b * y + c`.
*
* `(-b, a)` corresponds to a 2D vector parallel to the line. The lines
* have a natural orientation, corresponding to that vector.
*
* The magnitude ("norm") of the line is `sqrt(a ^ 2 + b ^ 2)`.
* `c / norm(line)` is the oriented distance from line to origin.
*/
// Returns line with direction (x, y) through origin
export const vector = (x: number, y: number): Line =>
GA.normalized([0, 0, -y, x, 0, 0, 0, 0]);
// For equation ax + by + c = 0.
export const equation = (a: number, b: number, c: number): Line =>
GA.normalized([0, c, a, b, 0, 0, 0, 0]);
export const through = (from: Point, to: Point): Line =>
GA.normalized(GA.join(to, from));
export const orthogonal = (line: Line, point: Point): Line =>
GA.dot(line, point);
// Returns a line perpendicular to the line through `against` and `intersection`
// going through `intersection`.
export const orthogonalThrough = (against: Point, intersection: Point): Line =>
orthogonal(through(against, intersection), intersection);
export const parallel = (line: Line, distance: number): Line => {
const result = line.slice();
result[1] -= distance;
return result as unknown as Line;
};
export const parallelThrough = (line: Line, point: Point): Line =>
orthogonal(orthogonal(point, line), point);
export const distance = (line1: Line, line2: Line): number =>
GA.inorm(GA.meet(line1, line2));
export const angle = (line1: Line, line2: Line): number =>
Math.acos(GA.dot(line1, line2)[0]);
// The orientation of the line
export const sign = (line: Line): number => Math.sign(line[1]);

View file

@ -0,0 +1,42 @@
import * as GA from "./ga";
import * as GALine from "./galines";
import type { Point, Line } from "./ga";
import { join } from "./ga";
export const from = ([x, y]: readonly [number, number]): Point => [
0,
0,
0,
0,
y,
x,
1,
0,
];
export const toTuple = (point: Point): [number, number] => [point[5], point[4]];
export const abs = (point: Point): Point => [
0,
0,
0,
0,
Math.abs(point[4]),
Math.abs(point[5]),
1,
0,
];
export const intersect = (line1: Line, line2: Line): Point =>
GA.normalized(GA.meet(line1, line2));
// Projects `point` onto the `line`.
// The returned point is the closest point on the `line` to the `point`.
export const project = (point: Point, line: Line): Point =>
intersect(GALine.orthogonal(line, point), line);
export const distance = (point1: Point, point2: Point): number =>
GA.norm(join(point1, point2));
export const distanceToLine = (point: Point, line: Line): number =>
GA.joinScalar(point, line);

View file

@ -0,0 +1,41 @@
import * as GA from "./ga";
import type { Line, Direction, Point, Transform } from "./ga";
import * as GADirection from "./gadirections";
/**
* TODO: docs
*/
export const rotation = (pivot: Point, angle: number): Transform =>
GA.add(GA.mul(pivot, Math.sin(angle / 2)), Math.cos(angle / 2));
export const translation = (direction: Direction): Transform => [
1,
0,
0,
0,
-(0.5 * direction[5]),
0.5 * direction[4],
0,
0,
];
export const translationOrthogonal = (
direction: Direction,
distance: number,
): Transform => {
const scale = 0.5 * distance;
return [1, 0, 0, 0, scale * direction[4], scale * direction[5], 0, 0];
};
export const translationAlong = (line: Line, distance: number): Transform =>
GA.add(GA.mul(GADirection.orthogonalToLine(line), 0.5 * distance), 1);
export const compose = (motor1: Transform, motor2: Transform): Transform =>
GA.mul(motor2, motor1);
export const apply = (
motor: Transform,
nvector: Point | Direction | Line,
): Point | Direction | Line =>
GA.normalized(GA.mul(GA.mul(motor, nvector), GA.reverse(motor)));