I am not a professional mathematics theorist. These are just my personal notes. There are probably errors here. You’ve been warned!
Resources

This trignometry web ap is very cool.

A nice discussion of eigenvectors with active graphical examples.

PDF of "Data Structures and Algorithm Analysis" (C version) by Clifford Shaffer.

Not exactly math, but a nice discussion of consistency models.

A nice linear algebra mini textbook.
Vector Math
I feel vaguely guilty building a transformation matrix by concatenating rotation / transform / scale instead of directly.
— John Carmack
Frankly, I find it excessively difficult and pointless to follow the gobbledygook notations of people like Pascal, LaPlace, Euler, and their pals. I can do it, but it’s a lot of work. Why not cut the ambiguity and go right to code? For all practical purposes, this will suffice. I keep needing to do geometry things and I think about looking up how to do that and writing code. Well, I’ve already written this kind of code a dozen times so just having it in one convenient place is what I really need. Here is one of my vector libraries simplified a bit; this won’t compile, but it should be clear what it does and how. It is worth noting that this is optimized for clarity, not performance. I usually have two functions for most things, one to act on the object itself and one to return an object; I mostly left out the former functions (referred to with "set" by the other code). If you need performance or just intense fanciness, why not check out third party libraries like…

Eigen C++ template library, etc. Recommended by SDCarND. (
aptget install libeigen3dev libeigen3doc
)
Vector(){ x = 0; y = 0; z = 0; }
Vector(float X, float Y, float Z){ x=X; y=Y; z=Z; }
~Vector(){ }
float valX(){ return x; }
float valY(){ return y; }
float valZ(){ return z; }
Vector inverse() { Vector I; I.x=x; I.y=y; I.z=z; return I; }
float magnitude() { float m=sqrt(x*x + y*y + z*z); return m;}
float distance(Vector P2) { Vector P=(*thisP2); return P.magnitude(); }
Vector direction() {
float mag;
Vector unit;
mag = magnitude(); // this>magnitude() works too...
if (mag > 0) {
unit.x = valX() / mag;
unit.y = valY() / mag;
unit.z = valZ() / mag; }
else { unit.x = 0; unit.y = 0; unit.z =0;}
return unit; }
Vector add(Vector offset) {
Vector sum;
sum.x = x + offset.x;
sum.y = y + offset.y;
sum.z = z + offset.z;
return sum; }
Vector subtract(Vector offset) {
Vector difference;
difference.x = x  offset.x;
difference.y = y  offset.y;
difference.z = z  offset.z;
return difference; }
Vector rotate(float &rads) {
Vector rotated;
float s, c;
s=sin(rads); c=cos(rads);
rotated.x = x*c  y*s;
rotated.y= y*c + x*s;
rotated.z= z;
return rotated; }
Vector rotate(float &rads, Vector &basept) {
Vector rotated;
float y2, s, c;
rotated = this>subtract(basept);
s=sin(rads); c=cos(rads);
y2= rotated.y*c + rotated.x*s;
rotated.x = rotated.x*c  rotated.y*s;
rotated.y = y2;
rotated.set2_add(basept);
return rotated; }
/* rotation in XZ axis */
Vector tilt(float &rads) {
Vector rotated;
float s, c;
s=sin(rads); c=cos(rads);
rotated.x = x*c  y*s;
rotated.z= z*c + z*s;
rotated.y= y;
return rotated; }
Vector tilt(float &rads, Vector &basept) {
Vector rotated;
float z2, s, c;
rotated = this>subtract(basept);
s=sin(rads); c=cos(rads);
z2= rotated.z*c + rotated.x*s;
rotated.x = rotated.x*c  rotated.z*s;
rotated.z = z2;
rotated.set2_add(basept);
return rotated; }
Vector viewvector(Vector &VV) {
if ((VV.x != 0)  (VV.y != 0)) {
Vector T, Xplus(1,0,0);
float r, t;
T.x=VV.x; T.y=VV.y; T.z=0;
r= ( T.y<0 ? T.angle(Xplus) : T.angle(Xplus) );
t = (M_PI/2)  atan(VV.z / sqrt(VV.x * VV.x + VV.y * VV.y) );
T = rotate(r);
T.set2_tilt(t);
r=T.x; T.x=T.y; T.y=r; // turn 90deg Z always points up
return T;}
else return *this; }
Vector perspective(float &lens, float &camdist) {
Vector T; float tz;
tz = camdist  z;
if (tz) {
T.x= (x * lens) / tz;
T.y= (y * lens) / tz;}
return T; }
Vector scale(float &fact) {
Vector scaled;
scaled.x = x*fact;
scaled.y= y*fact;
scaled.z= z*fact;
return scaled; }
Vector scale(float &fact, Vector &basept) {
Vector scaled;
scaled = this>subtract(basept);
scaled.x = scaled.x*fact;
scaled.y = scaled.y*fact;
scaled.z = scaled.z*fact;
scaled.set2_add(basept);
return scaled; }
Vector scale(Vector &fact3d) {
Vector scaled;
scaled.x = x * fact3d.x;
scaled.y = y * fact3d.y;
scaled.z = z * fact3d.z;
return scaled; }
Vector scale(Vector &fact3d, Vector &basept) {
Vector scaled;
scaled = this>subtract(basept);
scaled.x = scaled.x * fact3d.x;
scaled.y = scaled.y * fact3d.y;
scaled.z = scaled.z * fact3d.z;
scaled.set2_add(basept);
return scaled; }
float dot(Vector &B) {
float dotprod;
dotprod = x * B.x + y * B.y + z * B.z;
return dotprod; }
float angle(Vector &B) {
float ang;
ang = magnitude() * B.magnitude();
if ( ang == 0 ) ang = 0;
else ang = acos( dot(B) / ang ) ;
if (!(ang == ang)) return 0; /* nan check!!!
else return ang; }
Vector midpoint(Vector &B) {
return (((*this  B) * .5) + B); }
float angle(Vector B, Vector basept) {
set2_subtract(basept); B.set2_subtract(basept);
return angle(B); }
Vector cross(Vector &B) {
Vector C;
C.x = y * B.z  z * B.y;
C.y = z * B.x  x * B.z;
C.z = x * B.y  y * B.x;
return C; }
Geometry
Miscellaneous geometry notes. Possibly erroneous! Double check!
Point On A Line Perpendicular To Another Point
Points A and B are line endpoints, Point P is the point to get a perpendicular from off the line, Point N is the point on the line that is where P lines up at 90deg.
NA dot NP = 0
NB dot NP = 0
AN dot BN = 1
(AxNx)(PxNx)+(AyNy)(PyNy)+(AzNz)(PzNz)=0
(BxNx)(PxNx)+(ByNy)(PyNy)+(BzNz)(PzNz)=0
(AxNx)(BxNx)+(AyNy)(BYNy)+(AzNz)(BzNz)=1
(AxNx)(PxNx)+(AyNy)(PyNy)+(AzNz)(PzNz)=0
AxPxNxPxNxAx+Nx^2 + AyPyNyPyNyAy+Ny^2 + AzPzNzPzNzAz+Nz^2 = 0
BxPxNxPxNxBx+Nx^2 + ByPyNyPyNyBy+Ny^2 + BzPzNzPzNzBz+Nz^2 = 0
BxAxNxAxNxBx+Nx^2 + ByAyNyAyNyBy+Ny^2 + BzAzNzAzNzBz+Nz^2 = 1
AxPx+NxPx+NxAxNx^2 = AyPyNyPyNyAy+Ny^2 + AzPzNzPzNzAz+Nz^2
Nx^2 = AyPyNyPyNyAy+Ny^2 + AzPzNzPzNzAz+Nz^2 + AxPx+ NxPx +NxAx
Nx^2 = AyPy+NyPy+NyAyNy^2  AzPz+NzPz+NzAzNz^2  AxPx NxPx NxAx
+AxPx NxPx NxAx AyPy +NyPy +NyAy Ny^2  AzPz +NzPz +NzAz Nz^2 AxPx NxPx NxAx +AyPy NyPy NyAy +Ny^2 +AzPz NzPz NzAz +Nz^2 = 0
+AxPx AyPy +NyPy +NyAy Ny^2  AzPz +NzPz +NzAz Nz^2 AxPx +AyPy NyPy NyAy +Ny^2 +AzPz NzPz NzAz +Nz^2 = NxPx +NxAx +NxPx +NxAx
2(Px +Ax)Nx = +AxPx AyPy +NyPy +NyAy Ny^2  AzPz +NzPz +NzAz Nz^2 AxPx +AyPy NyPy NyAy +Ny^2 +AzPz NzPz NzAz +Nz^2
Nx = [(+AxPx AyPy +NyPy +NyAy Ny^2  AzPz +NzPz +NzAz Nz^2 AxPx +AyPy NyPy NyAy +Ny^2 +AzPz NzPz NzAz +Nz^2)/2(Px+Ax)]
BxPx[(+AxPx AyPy +NyPy +NyAy Ny^2  AzPz +NzPz +NzAz Nz^2 AxPx +AyPy NyPy NyAy +Ny^2 +AzPz NzPz NzAz +Nz^2)/2(Px+Ax)]Px[(+AxPx AyPy +NyPy +NyAy Ny^2  AzPz +NzPz +NzAz Nz^2 AxPx +AyPy NyPy NyAy +Ny^2 +AzPz NzPz NzAz +Nz^2)/2(Px+Ax)]Bx+[(+AxPx AyPy +NyPy +NyAy Ny^2  AzPz +NzPz +NzAz Nz^2 AxPx +AyPy NyPy NyAy +Ny^2 +AzPz NzPz NzAz +Nz^2)/2(Px+Ax)]^2 + ByPyNyPyNyBy+Ny^2 + BzPzNzPzNzBz+Nz^2 = 0
Midpoint Of A Line
Line with endpoint vectors A and B:
midpoint= ((A  B) * .5) + B
Intersection Of Line And Plane

N=(Nx,Ny,Nz) normal vector of plane

A=(Ax,Ay,Az) one endpoint of the line in question

B=(Bx,By,Bz) other endpoint of the line in question

P=(Px,Py,Pz) the point of intersection on plane along line

Px=(BxAx)alpha+Ax alpha=ratio of P along line AB

alpha=(PxAx)/(BxAx)=(PyAy)/(ByAy)=(PzAz)/(BzAz)
Three Equations For Three Variables Px,Py,Pz
1. NxPx + NyPy + NzPz = 0 or N dot P=0
2. Py=(PxAx)(ByAy)/(BxAx) + Ay
3. Pz=(PxAx)(BzAz)/(BxAx) + Az
0=NxPx+Ny((PxAx)(ByAy)/(BxAx)+Ay)+Nz((PxAx)(BzAz)/(BxAx)+Az)
<after much fussy and punctilious rearranging>
Ny(AyBxAxBy) + Nz(AzBxAxBz)
Px = 
Nx(AxBx) + Ny(AyBy) + Nz(AzBz)
Intersection Of Two 2d Line Segments
Line 1 is (Ax,Ay) to (Bx,By)
Line 2 is (Cx,Cy) to (Dx,Dy)
y1 = m1*x1+b1
y2 = m2*x2+b2
Where the lines intersect, x1=x2 and y1=y2, so:
m1*x+b1 = m2*x+b2
(m1m2)*x = b2b1
x = (b2b1)/(m1m2)
Ah, but what is intercept b and slope m?
m1 = (ByAy)/(BxAx)
m2 = (DyCy)/(DxCx)
And now the intercept b?
Ay = m1*Ax + b1
Cy = m2*Cx + b2
b1 = Ay  (m1*Ax)
b2 = Cy  (m2*Cx)
b1 = Ay  ((ByAy)/(BxAx)*Ax)
b2 = Cy  ((DyCy)/(DxCx)*Cx)
Recall:
x = (b2b1)/(m1m2)
y = m1*x + b1 (or (y = m2*x + b2) at intersection)
Fill in b, m and the intersection point is x,y:
x = {[Cy((DyCy)/(DxCx)*Cx)][Ay((ByAy)/(BxAx)*Ax)]} /
{[(ByAy)/(BxAx)][(DyCy)/(DxCx)]}
y = [(ByAy)/(BxAx)]*x + [Ay  ((ByAy)/(BxAx)*Ax)]
Of course this should be set up to catch special cases like parallel lines.
Intersection Of A Plane And A Z Vector At Arbitrary X,Y
Note

UNTESTED 
This is a special case of the line/plane intersection. It is a simplified case with a particular use in Hidden Line Removal (HLR). When the view is lined up along the Z axis, and a point’s X,Y is within an opaque triangle’s X,Y borders, then it is important to know whether the point is nearer to the viewer than the opaque triange (visable) or farther (hidden). The Z, or nearness to viewer, of the point is known. How near the point on the triangle is at X,Y is not automatically known; the triangle could be quite oblique to the view.
Setup:
Coordinates of the triangular face:
(T1x,T1y,T1z) (T2x,T2yT2z) (T3x,T3y,T3z)
Coordinates of the intersecting point:
(Ax,Ay,Az) Ax and Ay are known constants, Az is the variable
First move the whole mess so that a triangle vertex is 0,0,0:
T2x = T1x; T2y = T1y; T2z = T1z (i.e. T2x=T2xT1x, etc.)
T3x = T1x; T3y = T1y; T3z = T1z
Ax = T1x; Ay = T1y
T1 now is effectively (0,0,0)
The triangle needs a useful handle in the form of a normal vector at a vertex (the new origin vertex T1 is a good choice):
N = T2 crossproduct T3
Nx= T2y*T3z  T2z*T3y
Ny= T2x*T3z  T2z*T3x
Nz= T2x*T3y  T2y*T3x
If N is normal to the tri and A is on the tri, then they must be perpendicular and satisfy this equation:
A dotproduct N = 0
And that's it. One unknown, one equation. Onward:
Nx*Ax + Ny*Ay + Nz*Az = 0
Nz*Az = Nx*Ax  Ny*Ay
Az = (Nx*Ax  Ny*Ay)/Nz
Now put the thing back where it was (only Z is necessary):
Az += T1z;
So, the complete deal with everything included should be:
Az = {[[(T2yT1y)*(T3zT1z)  (T2zT1z)*(T3yT1y)]*(AxT1x) 
[(T2xT1x)*(T3zT1z)  (T2zT1z)*(T3xT1x)]*(AyT1y)] /
[(T2xT1x)*(T3yT1y)  (T2yT1y)*(T3xT1x)]} + T1z
Intersection Of A Line And An Intersecting Z Vector
Calculates Z for a point of known X,Y on a 3d line with known endpoints. The goal of this is to determine which of 2 lines is in front of the other. If two lines are found to intersect in xy dimensions (i.e. a Z aligned view) at X,Y, then what are the points for Z for each of the lines. This approach looks at each line one at a time. So there is the xy screen intersection point, (Px,Py), and a line to find the depth of contact of, (Ax,Ay)→(Bx,By). The object is to find Pz for this line. Then find the Z at this point for the other line and compare.
There seem to be two simple approaches to this  one slightly simpler than the other. Since AB is a line, things are linear. This means that the Z value will change with respect to X or Y. So one approach is to drop X or Y since only one is needed to relate to the Z. Here’s that approach:
Ignoring Y Ignoring X
AxBx/AzBz=AxPx/AzPz AyBy/AzBz=AyPy/AzPz
AzPz=(AxPx)*(AzBz)/AxBx AzPz=(AyPy)*(AzBz)/AyBy
Pz=Az(AxPx)*(AzBz)/AxBx Pz=Az(AyPy)*(AzBz)/AyBy
The problem with this approach is that if the Y is ignored, and the line is mostly along the Y axis, then the chance of precision errors gets big. Of course if the line is completely on the Y axis, there is a division by zero problem. One approach is to wisely choose which of the above formulas to use by this heuristic:
if (AxBx) > (AyBy) then (Ignore Y) else (Ignore X)
The grand complete solution is to not use progress along the line in some specific axis, but rather the 2d distance from the starting point as the independent variable. This means that when the point has travelled halfway along the line in XY, it has travelled half of the way in Z. This is opposed to saying that when it has travelled half of the X displacement, then it has also moved half of the Z. This may be true, but the X displacement might be 0. Here’s the full formula:
sqrt[(AxBx)^2+(AyBy)^2]/(AzBz) = sqrt[(AxPx)^2+(AyPy)^2]/(AzPz)
Pz=Az{sqrt[(AxPx)^2+(AyPy)^2]*(AzBz)/sqrt[(AxBx)^2+(AyBy)^2]}
Yuck…
Graph Threory
Linear Algebra and Matrix Operations
/* Input: */ /* A= matrix of dimension x by y. */ /* B= matrix of dimension y by z. */ /* Output: */ /* C= matrix of dimension x by z /* where C[i][j] is dot product of ith row of A and jth column of B */ for (i=1;i<=x;i++) for (j=1;j<=z;j++) { C[i][j]= 0; for (k=1;k<=y;k++) C[i][j]+= A[i][k] * B[k][j]; }
Here is a small comprehensible 2d matrix library from Sebastian Thrun (and Ernesto P. Adorio for the Cholesky functions). I find this far easier to understand than the obfuscatory typography gloop found in math books. Not only does this clearly explain how matrix operations work in unambiguous terms, you can actually use this library directly to do useful things.
#!/usr/bin/python3 from math import * class matrix: def __init__(self, value): self.value = value self.dimx = len(value) self.dimy = len(value[0]) if value == [[]]: self.dimx = 0 def zero(self, dimx, dimy): assert dimx >= 1 or dimy >= 1, "Invalid size of matrix" self.dimx = dimx self.dimy = dimy self.value = [[0 for row in range(dimy)] for col in range(dimx)] def identity(self, dim): assert dim >= 1, "Invalid size of matrix" self.dimx = dim self.dimy = dim self.value = [[0 for row in range(dim)] for col in range(dim)] for i in range(dim): self.value[i][i] = 1 def show(self): for i in range(self.dimx): print(self.value[i]) print(' ') def __add__(self, other): assert self.dimx == other.dimx or self.dimy == other.dimy, \ "Matrices must be of equal dimensions to add" res = matrix([[]]) res.zero(self.dimx, self.dimy) for i in range(self.dimx): for j in range(self.dimy): res.value[i][j] = self.value[i][j] + other.value[i][j] return res def __sub__(self, other): assert self.dimx == other.dimx or self.dimy == other.dimy, \ "Matrices must be of equal dimensions to subtract" res = matrix([[]]) res.zero(self.dimx, self.dimy) for i in range(self.dimx): for j in range(self.dimy): res.value[i][j] = self.value[i][j]  other.value[i][j] return res def __mul__(self, other): assert self.dimy == other.dimx, \ "Matrices must be m*n and n*p to multiply" res = matrix([[]]) res.zero(self.dimx, other.dimy) for i in range(self.dimx): for j in range(other.dimy): for k in range(self.dimy): res.value[i][j] += self.value[i][k] * other.value[k][j] return res def transpose(self): res = matrix([[]]) res.zero(self.dimy, self.dimx) for i in range(self.dimx): for j in range(self.dimy): res.value[j][i] = self.value[i][j] return res def Cholesky(self, ztol=1.0e5): # Computes the upper triangular Cholesky factorization of # a positive definite matrix. res = matrix([[]]) res.zero(self.dimx, self.dimx) for i in range(self.dimx): S = sum([(res.value[k][i])**2 for k in range(i)]) d = self.value[i][i]  S if abs(d) < ztol: res.value[i][i] = 0.0 else: assert d >= 0.0, "Matrix not positivedefinite" res.value[i][i] = sqrt(d) for j in range(i+1, self.dimx): S = sum([res.value[k][i] * res.value[k][j] for k in range(self.dimx)]) if abs(S) < ztol: S = 0.0 res.value[i][j] = (self.value[i][j]  S)/res.value[i][i] return res def CholeskyInverse(self): # Computes inverse of matrix given its Cholesky upper Triangular # decomposition of matrix. res = matrix([[]]) res.zero(self.dimx, self.dimx) # Backward step for inverse. for j in reversed(range(self.dimx)): tjj = self.value[j][j] S = sum([self.value[j][k]*res.value[j][k] for k in range(j+1, self.dimx)]) res.value[j][j] = 1.0/tjj**2  S/tjj for i in reversed(range(j)): res.value[j][i] = res.value[i][j] = sum([self.value[i][k]*res.value[k][j] for k in range(i+1, self.dimx)])/self.value[i][i] return res def inverse(self): aux = self.Cholesky() res = aux.CholeskyInverse() return res def __repr__(self): return repr(self.value)