commit c0f599e563d8b7fcba944c832dffaf7e6cd98e5f
Author: Kyle Isom <kyle@imap.cc>
Date:   Mon Mar 7 23:47:39 2022 -0800

    Initial import.

diff --git a/LICENSE b/LICENSE
new file mode 100644
index 0000000..fafc51e
--- /dev/null
+++ b/LICENSE
@@ -0,0 +1,15 @@
+ISC License (ISC)
+
+Copyright 2022 Kyle Isom <kyle@imap.cc>
+
+Permission to use, copy, modify, and/or distribute this software for any
+purpose with or without fee is hereby granted, provided that the above
+copyright notice and this permission notice appear in all copies.
+
+THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
+WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
+MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
+ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
+WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
+ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
+OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
diff --git a/README.rst b/README.rst
new file mode 100644
index 0000000..a7ca55f
--- /dev/null
+++ b/README.rst
@@ -0,0 +1,5 @@
+kmath
+=====
+
+While learning maths, I've found it help to write code to explore the ideas. If I can implement them, it helps me understand them.
+
diff --git a/kmath/__init__.py b/kmath/__init__.py
new file mode 100644
index 0000000..e69de29
diff --git a/kmath/linea.py b/kmath/linea.py
new file mode 100644
index 0000000..b1e4377
--- /dev/null
+++ b/kmath/linea.py
@@ -0,0 +1,308 @@
+# -*- coding: utf-8 -*-
+"""
+``linea.vector``
+
+Arbitrarily-sized vectors built on numpy's ndarray.
+
+Components:
+
++ class Vector
++ exception NonConformantVectors
++ function: dot
++ function: angle
++ function: parallel
++ function: orthogonal
+"""
+# pylint: disable=C0103
+import math
+
+import numpy
+
+EQUALITY_TOLERANCE = 0.001
+
+
+class NonConformantVectors(Exception):
+    """
+    A NonConformantVector is thrown when attempting to do operations
+    with vectors of differing sizes; generally, the size of the
+    /right-hand/ argument doesn't conform to the size of the /left-hand/
+    argument.
+    """
+
+    def __init__(self, expected, actual):
+        """
+        Initialize a new NonConformantVectors exception.
+        :param expected: length of the first vector
+        :param actual: length of the second vector
+        """
+        self.expected = expected
+        self.actual = actual
+        Exception.__init__(self)
+
+    def __str__(self):
+        msg = "Expected the right-hand vector to have dimension {}, "
+        msg += "but it has a dimension of {}."
+        return msg.format(self.expected, self.actual)
+
+
+class Vector:
+    """
+    A vector is a one-dimensional vector of some arbitrary size. This size
+    is fixed and can't be changed later in the Vector's life.
+    """
+
+    def __init__(self, a=None, *args):
+        """
+        Initialise a vector, either using an iterable passed in or as a sequence
+        of values.
+        >>> print(Vector(1, 2, 3))
+        [1; 2; 3]
+        >>> print(Vector([1, 2, 3]))
+        [1; 2; 3]
+        """
+        if len(args) > 0:
+            a = [a]
+            a.extend(args)
+            self.v = numpy.array(a)
+        else:
+            if isinstance(a, numpy.ndarray):
+                self.v = a
+            else:
+                self.v = numpy.array(a)
+
+    def __str__(self):
+        s = '['
+        for i in range(len(self.v)):
+            s += str(self.v[i])
+            if i < len(self.v) - 1:
+                s += '; '
+        s += ']'
+        return s
+
+    def __len__(self):
+        return len(self.v)
+
+    def __iter__(self):
+        return (x for x in self.v)
+
+    def __mul__(self, other):
+        o = list(map(lambda x: x * other, self.v))
+        return Vector(a=o)
+
+    def __rmul__(self, other):
+        return self * other
+
+    def __add__(self, other):
+        if len(self) != len(other):
+            raise NonConformantVectors(len(self), len(other))
+        return Vector(self.v + other.v)
+
+    def __radd__(self, other):
+        return self + other
+
+    def __sub__(self, other):
+        if len(self) != len(other):
+            raise NonConformantVectors(len(self), len(other))
+        return Vector(self.v - other.v)
+
+    def __eq__(self, other):
+        if not isinstance(other, Vector):
+            raise ValueError
+        if len(self) != len(other):
+            raise NonConformantVectors(len(self), len(other))
+        eq = numpy.isclose(self.v, other.v, EQUALITY_TOLERANCE)
+        if isinstance(eq, bool):
+            return eq
+        return eq.all()
+
+    def __repr__(self):
+        return 'Vector[{}]'.format(len(self))
+
+    def __getitem__(self, item):
+        return self.v[item]
+
+    def magnitude(self):
+        """
+        Return the magnitude of the vector.
+        """
+        return math.sqrt(sum(map(lambda x: x * x, self.v)))
+
+    def is_zero(self, tolerance=EQUALITY_TOLERANCE):
+        """
+        Return True if the vector is a zero vector (within some tolerance).
+        """
+        return numpy.isclose(self.magnitude(), 0, tolerance)
+
+    def unit(self):
+        """
+        Return the unit vector of this vector. If this method is called on
+        a zero vector (i.e. is_zero returns True), a ValueError will be thrown.
+        """
+        mag = self.magnitude()
+        if self.is_zero():
+            raise ValueError("cannot normalise the zero vector")
+        return self * (1 / mag)
+
+    def dot(self, other):
+        """
+        Compute the dot product between this vector and the other vector.
+        """
+        return dot(self, other)
+
+    def angle_with(self, other, in_degrees=False):
+        """
+        Compute the angle between this vector and the other vector. The
+        default is to return the angle in radians. If in_degrees is True,
+        returns the angle in degrees.
+        """
+        return angle(self, other, in_degrees=in_degrees)
+
+    def parallel_to(self, other):
+        """Return True if the vector other is parallel to this vector."""
+        return parallel(self, other)
+
+    def orthogonal_to(self, other):
+        """Return True if the vector other is orthogonal to this vector."""
+        return orthogonal(self, other)
+
+    def project_parallel(self, basis):
+        """
+        Return the projection of this vector onto the given basis vector.
+        """
+        unit_basis = basis.unit()
+        return Vector(dot(self, unit_basis) * unit_basis)
+
+    def project_orthogonal(self, basis):
+        """
+        Compute the orthogonal projection of the vector from the given basis
+        vector.
+        """
+        spar = self.project_parallel(basis)
+        return self - spar
+
+
+# The dot (or inner) product determines the angle between two vectors.
+def dot(v, w):
+    """
+    Return the dot product of vectors v and w.
+    :rtype: int
+    """
+    assert isinstance(v, Vector)
+    assert isinstance(v, Vector)
+    inner = sum([a * b for (a, b) in zip(v.v, w.v)])
+
+    # Cauchy-Schwartz inequality
+    assert abs(inner) <= (v.magnitude() * w.magnitude())
+    return inner
+
+
+def angle(v, w, in_degrees=False, tolerance=EQUALITY_TOLERANCE):
+    """
+    Return the angle between vectors v and w in radians. If in_degrees is
+    True, return the answer in degrees.
+    """
+    try:
+        # check for floating point problems resulting in domain errors
+        inner = dot(v.unit(), w.unit())
+        if inner > 1:
+            inner = clamp_if_close(inner, 1.0, tolerance)
+        if inner < -1:
+            inner = clamp_if_close(inner, -1.0, tolerance)
+    except NonConformantVectors:
+        raise ValueError('Cannot determine the angle between the zero vector and another vector.')
+    except:
+        # Deliberately pass through any other exceptions.
+        raise
+    theta = math.acos(inner)
+    if in_degrees:
+        theta = r2d(theta)
+    return theta
+
+
+def parallel(v, w, tolerance=EQUALITY_TOLERANCE):
+    """Return True if vectors v and w are parallel."""
+    if len(v) != len(w):
+        raise NonConformantVectors(len(v), len(w))
+    if v.is_zero() or w.is_zero():
+        return True
+
+    # HERE LIES THE RUIN OF A DUMB
+    # I'd originally tried to code this by creating a list of v_i / w_i, and then reducing
+    # the list via comparison (e.g. are all the values in the list numpy.isclose?). I got it
+    # right, but... the video showed a better way. Lesson learned, think about things instead
+    # of blindly coding through.
+    theta = angle(v, w)
+    if numpy.isclose(theta, 0, tolerance):
+        return True
+    return numpy.isclose(theta, math.pi, tolerance)
+
+
+def orthogonal(v, w, tolerance=EQUALITY_TOLERANCE):
+    """
+    Return True if vectors v and w are orthogonal.
+    """
+    if len(v) != len(w):
+        raise NonConformantVectors(len(v), len(w))
+    if v.is_zero() or w.is_zero():
+        return True
+    return numpy.isclose(dot(v, w), 0, tolerance)
+
+
+def cross(v, w):
+    """Return the cross product of the 3D vectors v and w."""
+    if len(v) != 3:
+        raise NonConformantVectors(3, len(v))
+    if len(w) != 3:
+        raise NonConformantVectors(3, len(w))
+
+    (x1, y1, z1) = v.v
+    (x2, y2, z2) = w.v
+    xc = (y1 * z2) - (y2 * z1)
+    # - (x1 * z2) - (x2 * z1)
+    yc = (x1 * z2) - (x2 * z1)
+    # x1 * y2 - x2 * y1
+    zc = (x1 * y2) - (x2 * y1)
+    return Vector(xc, -yc, zc)
+
+
+def area_parallelogram(v, w):
+    """Return the area of a parallelogram formed by Vectors v and w."""
+    u = cross(v, w)
+    return u.magnitude()
+
+
+def area_triangle(v, w):
+    """Return the area of a triangle formed by Vectors v and w."""
+    return 0.5 * area_parallelogram(v, w)
+
+
+def r2d(rval):
+    """
+    Convert the integer or floating point radian value to degrees. The
+    radian value must be between :math:`-2*\pi <= rval <= 2*\pi`.
+
+    >>> r2d(math.pi)
+    180.0
+    >>> r2d(3 * math.pi / 4)
+    135.0
+    >>> r2d(3 * math.pi)
+    ---------------------------------------------------------------------------
+    AssertionError                            Traceback (most recent call last)
+    """
+    assert abs(rval) <= (2 * math.pi)
+    return rval * 180 / math.pi
+
+
+def clamp_if_close(value, clamped, tolerance=EQUALITY_TOLERANCE):
+    """
+    Return clamped if value is close to clamped (within tolerance), or return
+    value otherwise.
+
+    >>> clamp_if_close(0.99, 1.0, tolerance=0.1)
+    1.0
+    >>> clamp_if_close(0.99, 1.0, tolerance=0.001)
+    0.99
+    """
+    if numpy.isclose(value, clamped, tolerance):
+        return clamped
+    return value
diff --git a/kmath/rat.py b/kmath/rat.py
new file mode 100644
index 0000000..226dd6d
--- /dev/null
+++ b/kmath/rat.py
@@ -0,0 +1,132 @@
+import math
+
+
+class Rational:
+    """
+    Rational implements a rational number type.
+    """
+
+    def __init__(self, n, d):
+        self.n = int(n)
+        self.d = int(d)
+
+        # hack to avoid DBZ
+        if d == 0:
+            self.n = 1
+            self.d = 0
+
+    def __float__(self) -> float:
+        return self.n / self.d
+
+    def __int__(self) -> int:
+        """
+        The integer version of a Rational is the rounded float representation
+        of the number.
+        """
+        return round(self.n / self.d)
+
+    def __add__(self, other):
+        lcm = int(math.lcm(self.d, other.d))
+        ks = lcm / self.d
+        ko = lcm / other.d
+
+        return Rational((self.n * ks) + (other.n * ko), lcm)
+
+    def __mul__(self, other):
+        return Rational(self.n * other.n, self.d * other.d)
+
+    def __truediv__(self, other):
+        return Rational(self.n * other.d, self.d * other.n)
+
+    def lcm(self, other):
+        return int(math.lcm(self.d, other.d))
+
+    def normalize(self, other):
+        lcm = self.lcm(other)
+        k = lcm / self.d
+        return self.scale(k)
+
+    def __normpair__(self, other):
+        rs = self.reduce()
+        ro = other.reduce()
+        ns = rs.normalize(ro)
+        no = ro.normalize(rs)
+        return ns, no
+
+    def reduce(self):
+        gcd = math.gcd(self.n, self.d)
+        return Rational(self.n / gcd, self.d / gcd)
+
+    def __eq__(self, other):
+        rs = self.reduce()
+        ro = other.reduce()
+        return (rs.n == ro.n) and (rs.d == ro.d)
+
+    def __gt__(self, other):
+        ns, no = self.__normpair__(other)
+        return ns.n > no.n
+
+    def __ge__(self, other):
+        ns, no = self.__normpair__(other)
+        return ns.n >= no.n
+
+    def __lt__(self, other):
+        ns, no = self.__normpair__(other)
+        return ns.n < no.n
+
+    def __le__(self, other):
+        ns, no = self.__normpair__(other)
+        return ns.n <= no.n
+
+    def scale(self, k):
+        return Rational(self.n * k, self.d * k)
+
+    def __sub__(self, other):
+        lcm = self.lcm(other)
+        ks = lcm / self.d
+        ko = lcm / other.d
+
+        return Rational((self.n * ks) - (other.n * ko), lcm)
+
+    def __pow__(self, power):
+        return Rational(self.n ** power, self.d ** power)
+
+    def __repr__(self):
+        return f'Rational({self.n}, {self.d})'
+
+
+def Int(n) -> Rational:
+    """
+    Int produces an integer version of a Rational (e.g. denominator is 1).
+    """
+    return Rational(n, 1)
+
+
+def ipow(n, m):
+    x = n
+
+    while m > 1:
+        x = x * n
+
+    return x
+
+
+def expansion(v):
+    top = Int(1)
+    bottom = Int(2) + v
+
+    return top / bottom
+
+
+def recurse(steps, start=None):
+    x = start
+    if x is None:
+        x = Rational(1, 2)
+    print(f'{x}')
+
+    while steps > 0:
+        x = expansion(x)
+        print(f'{x}')
+        steps -= 1
+
+    return x + Int(1)
diff --git a/requirements.txt b/requirements.txt
new file mode 100644
index 0000000..31ef97d
--- /dev/null
+++ b/requirements.txt
@@ -0,0 +1,2 @@
+numpy
+pytest
\ No newline at end of file
diff --git a/setup.py b/setup.py
new file mode 100644
index 0000000..989bd02
--- /dev/null
+++ b/setup.py
@@ -0,0 +1,23 @@
+# -*- coding: utf-8 -*-
+
+# Learn more: https://github.com/kennethreitz/setup.py
+
+from setuptools import setup, find_packages
+
+with open('README.rst') as f:
+    readme = f.read()
+
+with open('LICENSE') as f:
+    license_file = f.read()
+
+setup(
+    name='kmath',
+    version='0.1.0',
+    description='Maths code for experimentation',
+    long_description=readme,
+    author='Kyle Isom',
+    author_email='kyle@imap.cc',
+    url='https://git.wntrmute.dev/kyle/pymath',
+    license=license_file,
+    packages=find_packages(exclude=('tests', 'docs'))
+)
diff --git a/tests/__init__.py b/tests/__init__.py
new file mode 100644
index 0000000..e69de29
diff --git a/tests/linea.py b/tests/linea.py
new file mode 100644
index 0000000..1f64dce
--- /dev/null
+++ b/tests/linea.py
@@ -0,0 +1,167 @@
+import numpy
+import pytest
+
+import kmath.linea as vec
+
+
+def fequal(a, b):
+    eq = numpy.isclose(a, b, vec.EQUALITY_TOLERANCE)
+    if isinstance(eq, numpy.ndarray):
+        return eq.all()
+    return eq
+
+
+def test_equality():
+    v1 = vec.Vector(1, 2, 3)
+    v2 = vec.Vector(1, 2)
+    v3 = vec.Vector(1, 2, 3)
+    v4 = vec.Vector(3, 4, 5)
+
+    assert (v1 == v3)
+    with pytest.raises(vec.NonConformantVectors):
+        assert v1 != v2
+    assert (v1 != v4)
+
+
+# Video 4
+def test_basic_operations():
+    v1 = vec.Vector(8.218, -9.341)
+    v2 = vec.Vector(-1.129, 2.111)
+    v3 = vec.Vector(7.089, -7.230)
+    assert (v1 + v2 == v3)
+
+    v1 = vec.Vector(7.119, 8.215)
+    v2 = vec.Vector(-8.223, 0.878)
+    v3 = vec.Vector(15.3420, 7.3370)
+    assert (v1 - v2 == v3)
+
+    v1 = vec.Vector(1.671, -1.012, -0.318)
+    k = 7.41
+    v2 = vec.Vector(12.3821, -7.4989, -2.3564)
+    assert (k * v1 == v2)
+    assert (v1 * k == v2)
+
+
+# Video 6
+def test_magnitude():
+    v1 = vec.Vector(-0.221, 7.437)
+    assert (fequal(v1.magnitude(), 7.4403))
+
+    v2 = vec.Vector(5.581, -2.136)
+    unit = vec.Vector(0.933935214087, -0.357442325262)
+    assert (v2.unit() == unit)
+
+    v3 = vec.Vector(8.813, -1.331, -6.247)
+    assert (fequal(v3.magnitude(), 10.8842))
+
+    v4 = vec.Vector(1.996, 3.108, -4.554)
+    unit = vec.Vector(.340401295943, 0.530043701298, -0.776647044953)
+    assert (v4.unit() == unit)
+
+    v5 = vec.Vector(0, 0, 0)
+    with pytest.raises(ValueError):
+        v5.unit()
+
+
+# Video 8
+def test_dot_product():
+    v1 = vec.Vector(1, 2, -1)
+    v2 = vec.Vector(3, 1, 0)
+    assert (fequal(vec.dot(v1, v2), 5))
+    assert (fequal(v1.dot(v2), 5))
+
+    v1 = vec.Vector(7.887, 4.138)
+    v2 = vec.Vector(-8.802, 6.776)
+    assert (fequal(vec.dot(v1, v2), -41.382))
+
+    v3 = vec.Vector(-5.955, -4.904, -1.874)
+    v4 = vec.Vector(-4.496, -8.755, 7.103)
+    assert (fequal(vec.dot(v3, v4), 56.3971))
+
+    v5 = vec.Vector(3.183, -7.627)
+    v6 = vec.Vector(-2.668, 5.319)
+    assert (fequal(vec.angle(v5, v6), 3.072))
+
+    v7 = vec.Vector(7.35, 0.221, 5.188)
+    v8 = vec.Vector(2.751, 8.259, 3.985)
+    theta = vec.angle(v7, v8, in_degrees=True)
+    assert (fequal(theta, 60.2758))
+
+    v9 = vec.Vector(0, 0)
+    with pytest.raises(ValueError):
+        v1.angle_with(v9)
+
+
+# Video 10.
+def test_parallel_orthogonal():
+    v1 = vec.Vector(-7.579, -7.88)
+    v2 = vec.Vector(22.737, 23.64)
+    assert (v1.parallel_to(v2))
+    assert (not v1.orthogonal_to(v2))
+
+    v3 = vec.Vector(-2.029, 9.97, 4.172)
+    v4 = vec.Vector(-9.231, -6.639, -7.245)
+    assert (not v3.parallel_to(v4))
+    assert (not v3.orthogonal_to(v4))
+
+    v5 = vec.Vector(-2.328, -7.284, -1.214)
+    v6 = vec.Vector(-1.821, 1.072, -2.94)
+    assert (not v5.parallel_to(v6))
+    assert (v5.orthogonal_to(v6))
+
+    v7 = vec.Vector(2.118, 4.827)
+    v8 = vec.Vector(0, 0)
+    assert v7.parallel_to(v8)
+    assert v7.orthogonal_to(v8)
+
+
+# Video 12
+def test_projection():
+    # sanity check
+    v1 = vec.Vector(1, 3)
+    v2 = vec.Vector(3, 3)
+    v3 = v1.project_parallel(v2)
+    v4 = v1.project_orthogonal(v2)
+    assert v3 + v4 == v1
+
+    v1 = vec.Vector(3.039, 1.879)
+    v2 = vec.Vector(0.825, 2.036)
+    v3 = vec.Vector(1.0826, 2.6717)
+    assert v1.project_parallel(v2) == v3
+
+    v4 = vec.Vector(-9.88, -3.264, -8.159)
+    v5 = vec.Vector(-2.155, -9.353, -9.473)
+    v6 = vec.Vector(-8.350, 3.376, -1.434)
+    assert v4.project_orthogonal(v5) == v6
+
+    v7 = vec.Vector(3.009, -6.172, 3.692, -2.510)
+    v8 = vec.Vector(6.404, -9.144, 2.759, 8.718)
+    v9 = vec.Vector(1.969, -2.811, 0.848, 2.680)
+    assert v7.project_parallel(v8) == v9
+    v10 = vec.Vector(1.040, -3.361, 2.844, -5.190)
+    assert v7.project_orthogonal(v8) == v10
+
+    assert (v9 + v10) == v7
+
+
+def test_cross_product():
+    # sanity check
+    v1 = vec.Vector(5, 3, -2)
+    v2 = vec.Vector(-1, 0, 3)
+    v3 = vec.Vector(9, -13, 3)
+    assert vec.cross(v1, v2) == v3
+
+    v1 = vec.Vector(8.462, 7.893, -8.187)
+    v2 = vec.Vector(6.984, -5.975, 4.778)
+    v3 = vec.Vector(-11.205, -97.609, -105.685)
+    assert vec.cross(v1, v2) == v3
+
+    v4 = vec.Vector(-8.987, -9.838, 5.031)
+    v5 = vec.Vector(-4.268, -1.861, -8.866)
+    area = 142.122
+    assert fequal(vec.area_parallelogram(v4, v5), area)
+
+    v6 = vec.Vector(1.500, 9.547, 3.691)
+    v7 = vec.Vector(-6.007, 0.124, 5.772)
+    area = 42.565
+    assert fequal(vec.area_triangle(v6, v7), area)