diff --git a/README.md b/README.md
index a72c1df..f9c8222 100644
--- a/README.md
+++ b/README.md
@@ -1,6 +1,6 @@
 # Apatite
 
-Apatite is a fundimental package for scientific computing in Crystal. If that sounds like a modified version of the first line from the NumPy homepage, that's because it is. Apatite has (ok, will have) all of the goodness of NumPy sitting atop the blazing speed and beautiful syntax of Crystal.
+Apatite is meant to be a collecion of mathmatical and scientific computing algorithms for the Crystal programming language. I don't expect it to ever reach the level of completeness as numpy, but hopefully it can save some people the trouble of implementing these methods on their own.
 
 ## Installation
 
@@ -25,25 +25,16 @@ TODO: Write usage instructions here, but first write the library...
 ## Roadmap
 
 - [ ] Apetite
-	- [ ] Array Objects
-		- [ ] NArray (see [numpy.ndarray](https://docs.scipy.org/doc/numpy/reference/arrays.ndarray.html))
-		- [ ] DType  (see [numpy.dtype](https://docs.scipy.org/doc/numpy/reference/arrays.dtypes.html))
-		- [ ] Scalars
-		- [ ] Indexing
-	- [ ] Routines
-		- [ ] Binary Operations
-		- [ ] String Operations
-		- [ ] FFT    (see [numpy.fft](https://docs.scipy.org/doc/numpy/reference/routines.fft.html))
-		- [ ] Financial Functions
-		- [ ] LinAlg (see [numpy.linalg](https://docs.scipy.org/doc/numpy/reference/routines.linalg.html))
-		- [ ] Logic Functions
-		- [ ] Mathematical Functions
-		- [ ] Matlib (see [numpy.matlib](https://docs.scipy.org/doc/numpy/reference/routines.matlib.html))
-		- [ ] Padding Arrays
-		- [ ] Polynomials
-		- [ ] Random (see [numpy.random](https://docs.scipy.org/doc/numpy/reference/routines.random.html))
-		- [ ] Sorting, Searching, and Counting
-		- [ ] Statistics
+	- [ ] Vector Math
+		- [ ] Matrix
+		- [ ] NDArray
+		- [x] Vector
+		- [ ] Line
+		- [ ] Plane
+		- [ ] Polygon
+		- [ ] LinkedList
+	
+
 
 ## Contributing
 
diff --git a/spec/apatite/vector_spec.cr b/spec/apatite/vector_spec.cr
new file mode 100644
index 0000000..6f59cfa
--- /dev/null
+++ b/spec/apatite/vector_spec.cr
@@ -0,0 +1,161 @@
+require "../spec_helper"
+
+describe "Apatite::Vector" do
+  x = Apatite::Vector.new([3, 4])
+
+  describe "#==" do
+    it "is true if elements are the same" do
+      (x == Apatite::Vector.new([3, 4])).should be_true
+    end
+
+    it "is false if elements are different" do
+      (x == Apatite::Vector.new([7, 1])).should be_false
+    end
+  end
+
+  describe "#e" do
+    it "returns the `ith` element in the vector" do
+      x.e(2).should eq(4)
+    end
+
+    it "returns `nil` if `n` is out of range" do
+      x.e(7).should be_nil
+    end
+  end
+
+  describe "#to_unit_vector" do
+    it "successfully converts a vector to a unit vector" do
+      res = [0.6, 0.8]
+      x.to_unit_vector.each_with_index do |e, i|
+        e.should be_close(res[i], 0.1)
+      end
+    end
+  end
+
+  describe "#dimensions" do
+    it "gets vector's dimensions" do
+      x.dimensions.should eq({1, 2})
+    end
+  end
+
+  describe "#rows" do
+    it "gets vector's rows" do
+      x.rows.should eq(1)
+    end
+  end
+
+  describe "#cols" do
+    it "gets vector's columns" do
+      x.cols.should eq(2)
+    end
+  end
+
+  describe "#product" do
+    it "computes the product of a vector" do
+      x.product.should eq(12)
+    end
+  end
+
+  describe "#angle_from" do
+    it "should compute the angle between `y` and `z`" do
+      y = Apatite::Vector.create([1, 1])
+      z = Apatite::Vector.create([1, 0])
+      y.angle_from(z).should be_close(Math::PI / 4, 0.1)
+    end
+  end
+
+  describe "#parallel_to?" do
+    it "correctly determines if a vector is parallel to another" do
+      x.parallel_to?(Apatite::Vector.create([6, 8])).should be_true
+      x.parallel_to?(Apatite::Vector.create([1, 1])).should be_false
+    end
+  end
+
+  describe "#antiparallel_to?" do
+    it "correctly determines if a vector is antiparallel to another" do
+      x.antiparallel_to?(Apatite::Vector.create([-3, -4])).should be_true
+      x.antiparallel_to?(x).should be_false
+    end
+  end
+
+  describe "#perpendicular_to?" do
+    it "correctly determines if a vector is antiparallel to another" do
+      x.perpendicular_to?(Apatite::Vector.create([-4, 3])).should be_true
+      x.perpendicular_to?(x).should be_false
+    end
+  end
+
+  describe "#dot" do
+    it "calculates the dot product of the vector and another vector" do
+      x.dot(Apatite::Vector.create([2, 3])).should eq(18)
+    end
+  end
+
+  describe "#add" do
+    it "adds a number to every item in a vector" do
+      x.add(2).should eq([5, 6])
+    end
+
+    it "adds an enumerable to a vector" do
+      x.add([3, 2]).should eq([6, 6])
+    end
+  end
+
+  describe "#subtract" do
+    it "subtracts a number from every item in a vector" do
+      x.subtract(2).should eq([1, 2])
+    end
+
+    it "subtracts an enumerable from a vector" do
+      x.subtract([3, 2]).should eq([0, 2])
+    end
+  end
+
+  describe "#multiply" do
+    it "multiplies a number with every item in a vector" do
+      x.multiply(2).should eq([6, 8])
+    end
+
+    it "multiplies an enumerable with a vector" do
+      x.multiply([3, 2]).should eq([9, 8])
+    end
+  end
+
+  describe "#sum" do
+    it "sums all items in a vector" do
+      x.sum.should eq(7)
+    end
+  end
+
+  describe "#chomp" do
+    it "returns a new vector with the first `n` items of the old vector" do
+      x.chomp(1).should eq([4])
+    end
+  end
+
+  describe "#top" do
+    it "returns a new vector with the last `n` items of the old vector" do
+      x.top(1).should eq([3])
+    end
+  end
+
+  describe "#augment" do
+    it "creates a new vector with the elements fro vector b appended to those from vector a" do
+      y = x.clone
+      y.augment(Apatite::Vector.create([5])).should eq([3, 4, 5])
+    end
+  end
+
+  describe ".log" do
+    it "should calculate the log of the vector" do
+      pp x
+      x.log.should eq([1.0986122886681098, 1.3862943611198906])
+    end
+  end
+
+  it "should allow for scalar addition" do
+    a = Apatite::Vector.create([2, 3, 4])
+    b = a.add(1)
+    b.should eq([3, 4, 5])
+  end
+end
diff --git a/spec/apatite_spec.cr b/spec/apatite_spec.cr
index 89c0008..6951309 100644
--- a/spec/apatite_spec.cr
+++ b/spec/apatite_spec.cr
@@ -1,9 +1,4 @@
 require "./spec_helper"
 
 describe Apatite do
-  # TODO: Write tests
-
-  it "works" do
-    false.should eq(true)
-  end
 end
diff --git a/src/apatite.cr b/src/apatite.cr
index 9043244..39630c8 100644
--- a/src/apatite.cr
+++ b/src/apatite.cr
@@ -1,5 +1,8 @@
-require "./apatite/core_ext/*"
-require "./apatite/narray"
+require "./apatite/core_ext/array"
+
+require "./apatite/matrix/ndarray"
+require "./apatite/matrix/vector"
+require "./apatite/matrix/matrix"
 
 # Apatite is a fundimental package for scientific computing in Crystal. If that
 # sounds like a modified version of the first line from the NumPy homepage,
@@ -9,16 +12,8 @@ require "./apatite/narray"
 module Apatite
   VERSION = "0.1.0"
 
-  # def self.zeros(shape : Array(Int32))
-  #   curr = shape.shift
-  #   Array.new(curr) do
-  #     if shape.empty?
-  #       Array.new(curr, 0)
-  #     else
-  #       self.zeros(shape)
-  #     end
-  #   end
-  # end
+  class_property precision = 1e-6
+  class_property approx_precision = 1e-5
 
   ## ## ## ## ## ## ## ## ## ## ## ## ##
   # # Array Creation
@@ -133,6 +128,18 @@ module Apatite
   # # Array Manipulation
   ## ## ## ## ## ## ## ## ## ## ## ## ##
 
+  def self.shape(input : Iterable(U) | U) forall U
+    if input.is_a?(Iterable(U))
+      if input.size > 0
+        [input.size].concat(self.shape(input[0]))
+      else
+        [0]
+      end
+    else
+      [] of U
+    end
+  end
+
   def self.copyto
   end
 
@@ -531,10 +538,16 @@ module Apatite
   # # Logic Functions
   ## ## ## ## ## ## ## ## ## ## ## ## ##
 
-  def self.all?(arr : Enumerable)
+  def self.all?(arr : Enumerable, axis : Int32? | Array(Int32)? = nil, keepdims = false)
+    raise "axis of #{axis} out of bounds for array of size #{arr.size}" unless axis.nil? || axis < arr.size
+
+    if axis.nil? && !keepdims
+      arr.flatten.each { |i| return false unless i == true }
+      return true
+    end
   end
 
-  def self.any?
+  def self.any?(arr : Enumerable, axis : Int32? | Array(Int32) = nil, keepdims = false)
   end
 
   def self.finite?
diff --git a/src/apatite/core_ext/array.cr b/src/apatite/core_ext/array.cr
index 4ecb44c..23e2553 100644
--- a/src/apatite/core_ext/array.cr
+++ b/src/apatite/core_ext/array.cr
@@ -20,4 +20,8 @@ class Array(T)
     max = max_by { |i| i.is_a?(Array) ? i.size : i }
     max.is_a?(Array) ? [size, max.size] : [max]
   end
+
+  def to_vec
+    Apatite::Vector.create(self)
+  end
 end
diff --git a/src/apatite/matrix/matrix.cr b/src/apatite/matrix/matrix.cr
new file mode 100644
index 0000000..7831b14
--- /dev/null
+++ b/src/apatite/matrix/matrix.cr
@@ -0,0 +1,269 @@
+require "./vector"
+
+module Apatite
+  class Matrix(T) < Array(T)
+    include Indexable(T)
+
+    @rows : Int32
+
+    @cols : Int32
+
+    # Creates a Matrix with each element initialized as *value*.
+    def self.of(value : T, rows : Int32, cols : Int32)
+      Matrix(T).new(shape) { value }
+    end
+
+    # Creates a Matrix, invoking *initializer* with each pair of indices.
+    def self.build(shape, &initializer : UInt32, UInt32 -> T)
+      Matrix(T).new(shape) do |idx|
+        i = (idx / N).to_u32
+        j = (idx % N).to_u32
+        initializer.call(i, j)
+      end
+    end
+
+    # Creates a single-column `Matrix` from a `Vector`
+    def self.col_vector(vector)
+      col = vector.to_a
+      Matrix(T).from(col, {col.size, 1})
+    end
+
+    # Creates a Matrix from elements contained within a StaticArray.
+    #
+    # The matrix will be filled rows first, such that an array of
+    #
+    #   [1, 2, 3, 4]
+    #
+    # becomes
+    #
+    #   | 1  2 |
+    #   | 3  4 |
+    #
+    def self.from(list : Array(T), shape)
+      raise("Not enough elements to fill matrix") if list.size < shape[0] * shape[1]
+
+      Matrix(T).new(shape) do |idx|
+        list[idx]
+      end
+    end
+
+    # Creates a single-row `Matrix` from a `Vector`
+    def self.row_vector(vector)
+      row = vector.to_a
+      Matrix(T).from(row, {1, row.size})
+    end
+
+    # Build a zero matrix (all elements populated with zero value of the type
+    # isntance).
+    def self.zero(shape)
+      Matrix(T).new(shape) { T.zero }
+    end
+
+    # Build the idenity matrix for the instanced type and dimensions.
+    #
+    # `id` may be used to specify an identity element for the type. If unspecifed
+    # a numeric identity will be assumed.
+    # def self.identity(id = T.zero + 1, shape)
+    #   {{ raise("Matrix dimensions must be square") unless M == N }}
+
+    #   Matrix(T).(shape) do |i, j|
+    #     i == j ? id : T.zero
+    #   end
+    # end
+
+    # Creates Matrix, yielding the linear index for each element to provide an
+    # initial value.
+    def initialize(shape, &block : Int32 -> T)
+      raise("Matrix dimensions must be positive") if shape[0] < 0 || shape[1] < 0
+      @rows = shape[0]
+      @cols = shape[1]
+      @buffer = Pointer(T).malloc(size, &block)
+    end
+
+    # Equality. Returns `true` if each element in `self` is equal to each
+    # corresponding element in *other*.
+    def ==(other : Matrix(U)) forall U
+      {% if other.rows == rows && other.cols == cols %}
+        each_with_index do |e, i|
+          return false unless e == other[i]
+        end
+        true
+      {% else %}
+        false
+      {% end %}
+    end
+
+    # Equality with another object, or differently sized matrix. Always `false`.
+    def ==(other)
+      false
+    end
+
+    # Returns a new Matrix that is the result of performing a matrix addition with
+    # *other*
+    def +(other : Matrix)
+      merge(other) { |a, b| a + b }
+    end
+
+    # Returns a new Matrix that is the result of performing a matrix subtraction
+    # with *other*
+    def -(other : Matrix)
+      merge(other) { |a, b| a - b }
+    end
+
+    # Performs a matrix multiplication with *other*.
+    def *(other : Matrix)
+      raise("Dimension mismatch, cannot multiply a #{rows}x#{cols} by a #{other.rows}x#{other.cols}") unless cols == other.cols
+
+      Matrix(typeof(self[0] * other[0])).build(other.rows, other.cols) do |i, j|
+        pairs = row(i).zip other.col(j)
+        pairs.map(&.product).sum
+      end
+    end
+
+    # Performs a scalar multiplication with *other*.
+    def *(other)
+      map { |x| x * other }
+    end
+
+    # Retrieves the value of the element at *i*,*j*.
+    #
+    # Indicies are zero-based. Negative values may be passed for *i* and *j* to
+    # enable reverse indexing such that `self[-1, -1] == self[M - 1, N - 1]`
+    # (same behaviour as arrays).
+    def [](i : Int, j : Int) : T
+      idx = index i, j
+      to_unsafe[idx]
+    end
+
+    # Sets the value of the element at *i*,*j*.
+    def []=(i : Int, j : Int, value : T)
+      idx = index i, j
+      to_unsafe[idx] = value
+    end
+
+    # Gets the contents of row *i*.
+    def row(i : Int)
+      Vector.new(cols) { |j| self[i, j] }
+    end
+
+    # Gets the contents of column *j*.
+    def col(j : Int)
+      Vector.new(rows) { |i| self[i, j] }
+    end
+
+    # Yields the current element at *i*,*j* and updates the value with the
+    # block's return value.
+    def update(i, j, &block : T -> T)
+      idx = index i, j
+      to_unsafe[idx] = yield to_unsafe[idx]
+    end
+
+    # Apply a morphism to all elements, returning a new Matrix with the result.
+    def map(&block : T -> U) forall U
+      Matrix(U).new({rows, cols}) do |idx|
+        block.call to_unsafe[idx]
+      end
+    end
+
+    # ditto
+    def map_with_indices(&block : T, UInt32, UInt32 -> U) forall U
+      Matrix(U).from({rows, cols}) do |i, j|
+        block.call self[i, j], i, j
+      end
+    end
+
+    # ditto
+    def map_with_index(&block : T, Int32 -> U) forall U
+      Matrix(U).new({rows, cols}) do |idx|
+        block.call to_unsafe[idx], idx
+      end
+    end
+
+    # Apply an endomorphism to `self`, mutating all elements in place.
+    def map!(&block : T -> T)
+      each_with_index do |e, idx|
+        to_unsafe[idx] = yield e
+      end
+      self
+    end
+
+    # Merge with another similarly dimensions matrix, apply the passed block to
+    # each elemenet pair.
+    def merge(other : Matrix(U), &block : T, U -> _) forall U
+      raise("Dimension mismatch") unless other.rows == rows && other.cols == cols
+
+      map_with_index do |e, i|
+        block.call e, other[i]
+      end
+    end
+
+    # Gets an `Vector` of `Vector` representing rows
+    def row_vectors
+      Vector.new(rows) { |i| row(i) }
+    end
+
+    # Gets an `Array` of `Vector` representing columns
+    def col_vectors
+      Vector.new(cols) { |i| col(i) }
+    end
+
+    # Creates a new matrix that is the result of inverting the rows and columns
+    # of `self`.
+    def transpose
+      Matrix(T).build({cols, rows}) do |i, j|
+        self[j, i]
+      end
+    end
+
+    # Returns the dimensions of `self` as a tuple of `{rows, cols}`.
+    def dimensions
+      {rows, cols}
+    end
+
+    # Gets the capacity (total number of elements) of `self`.
+    def size
+      @rows * @cols
+    end
+
+    # Count of rows.
+    def rows
+      @rows
+    end
+
+    # Count of columns.
+    def cols
+      @cols
+    end
+
+    # Returns the element at the given linear index, without doing any bounds
+    # check.
+    #
+    # Used by `Indexable`
+    @[AlwaysInline]
+    def unsafe_fetch(index : Int)
+      to_unsafe[index]
+    end
+
+    # Returns the pointer to the underlying element data.
+    def to_unsafe : Pointer(T)
+      @buffer
+    end
+
+    # Map *i*,*j* coords to an index within the buffer.
+    def index(i : Int, j : Int) : T
+      i += rows if i < 0
+      j += cols if j < 0
+
+      raise IndexError.new if i < 0 || j < 0
+      raise IndexError.new unless i < rows && j < cols
+
+      i * cols + j
+    end
+
+    def inspect
+      ::String.build do |s|
+        s << "Matrix{" << self.join(", ") << "}"
+      end
+    end
+  end
+end
diff --git a/src/apatite/matrix/ndarray.cr b/src/apatite/matrix/ndarray.cr
new file mode 100644
index 0000000..d81ff87
--- /dev/null
+++ b/src/apatite/matrix/ndarray.cr
@@ -0,0 +1,32 @@
+module Apatite
+  class NDArray
+    include Enumerable(Float64)
+    include Indexable(Float64)
+    include Comparable(NDArray)
+
+    getter data : Array(Float64)
+
+    getter shape : Array(Int32)
+
+    delegate :[], to: @data
+    delegate :[]?, to: @data
+    delegate :[]=, to: @data
+    delegate :unsafe_fetch, to: @data
+    delegate :to_unsafe, to: @data
+    delegate :size, to: @data
+
+    def initialize(data : Array(Number), shape : Array(Int32)? = nil)
+      @data = data.is_a?(Array(Float64)) ? data.flatten : data.flatten.map(&.to_f64)
+      @shape = shape || [@data.size]
+    end
+
+    # Returns the absolute value of every item in the array
+    def abs
+      map { |e| e.abs }
+    end
+
+    # Returns the arccosine of each element in the current array.
+    def acos
+    end
+  end
+end
diff --git a/src/apatite/matrix/vector.cr b/src/apatite/matrix/vector.cr
new file mode 100644
index 0000000..b7c4cce
--- /dev/null
+++ b/src/apatite/matrix/vector.cr
@@ -0,0 +1,822 @@
+require "../../apatite"
+
+module Apatite
+  # Represents a mathematical vector, and also constitutes a row or column
+  # of a `Matrix`
+  class Vector
+    include Enumerable(Float64)
+    include Indexable(Float64)
+    include Comparable(Vector)
+
+    I = Vector.create([1.0, 0.0, 0.0])
+    J = Vector.create([0.0, 1.0, 0.0])
+    K = Vector.create([0.0, 0.0, 1.0])
+
+    @buffer : Pointer(Float64)
+    @capacity : Int32
+
+    # Returns the number of elements in the vector
+    getter size : Int32
+
+    # Create a new empty `Vector`.
+    def initialize
+      @size = 0
+      @capacity = 0
+      @buffer = Pointer(Float64).null
+    end
+
+    # Creates a new empty `Vector` backed by a buffer that is initially
+    # `initial_capacity` big.
+    #
+    # The *initial_capacity* is useful to avoid unnecessary reallocations
+    # of the internal buffer in case of growth. If you have an estimate
+    # of the maximum number of elements an vector will hold, the vector should
+    # be initialized with that capacity for improved performance.
+    #
+    # ```
+    # vec = Vector.new(5)
+    # vec.size # => 0
+    # ```
+    def initialize(initial_capacity : Int)
+      if initial_capacity < 0
+        raise ArgumentError.new("Negative array size: #{initial_capacity}")
+      end
+
+      @size = 0
+      @capacity = initial_capacity.to_i
+
+      if initial_capacity == 0
+        @buffer = Pointer(Float64).null
+      else
+        @buffer = Pointer(Float64).malloc(initial_capacity)
+      end
+    end
+
+    # Creates a new `Vector` of the given *size* filled with the same *value* in each position.
+    #
+    # ```
+    # Vector.new(3, 1.0) # => Vector{1.0, 1.0, 1.0}
+    # ```
+    def initialize(size : Int, value : Float64)
+      if size < 0
+        raise ArgumentError.new("Negative vector size: #{size}")
+      end
+
+      @size = size.to_i
+      @capacity = size.to_i
+
+      if size == 0
+        @buffer = Pointer(Float64).null
+      else
+        @buffer = Pointer(Float64).malloc(size, value)
+      end
+    end
+
+    # Creates a new `Vector` of the given *size* and invokes the given block once
+    # for each index of `self`, assigning the block's value in that index.
+    #
+    # ```
+    # Vector.new(3) { |i| (i + 1.0) ** 2.0 } # => Vector{1.0, 4.0, 9.0}
+    #
+    # vec = Vector.new(3) { 5.0 }
+    # vec # => Vector{5.0, 5.0, 5.0}
+    # ```
+    def self.new(size : Int, &block : Int32 -> Float64)
+      Vector.build(size) do |buffer|
+        size.to_i.times do |i|
+          buffer[i] = yield i
+        end
+        size
+      end
+    end
+
+    # Creates a new `Vector`, allocating an internal buffer with the given *capacity*,
+    # and yielding that buffer. The given block must return the desired size of the vector.
+    #
+    # This method is **unsafe**.
+    #
+    # ```
+    # Vector.build(3) do |buffer|
+    #   LibSome.fill_buffer_and_return_number_of_elements_filled(buffer)
+    # end
+    # ```
+    def self.build(capacity : Int) : self
+      vec = Vector.new(capacity)
+      vec.size = (yield vec.to_unsafe).to_i
+      vec
+    end
+
+    # Creates a new `Vector` from the elements of another `Indexable`
+    # collection.
+    #
+    # ```
+    # Vector.create([1, 2, 3]) # Vector{1.0, 2.0, 3.0}
+    # ```
+    def self.create(elements : Indexable(Number))
+      vec = Vector.new
+      elements.each do |el|
+        vec.push(el.to_f64)
+      end
+      vec
+    end
+
+    # Generates a random vector of size `n` with elements
+    # in `range`.
+    def self.random(n, range = Float64::MIN..Float64::MAX)
+      random = Random.new
+      vec = Vector.new(n)
+      n.times do
+        vec.push random.rand(range)
+      end
+      vec
+    end
+
+    # Returns a standard basis n-vector.
+    def self.basis(size, index)
+      raise ArgumentError.new("invalid size (#{size} for 1..)") if size < 1
+      raise ArgumentError.new("invalid index (#{index} for 0...#{size})") unless 0 <= index && index < size
+      vec = Vector.new(size, 0.0)
+      vec[index] = 1.0
+      vec
+    end
+
+    # Create a new vector of size `n` filled with zeros.
+    def self.zeros(n)
+      Vector.new(n, 0.0)
+    end
+
+    # Create a new vector of size `n` filled with ones.
+    def self.ones(n)
+      Vector.new(n, 1.0)
+    end
+
+    def self.[](*array)
+      Vector.create(array)
+    end
+
+    # Returns all elements that are within the given range.
+    #
+    # Negative indices count backward from the end of the vector (-1 is the last
+    # element). Additionally, an empty vector is returned when the starting index
+    # for an element range is at the end of the vector.
+    #
+    # Raises `IndexError` if the starting index is out of range.
+    #
+    # ```
+    # v = Vector{1.0, 2.0, 3.0, 4.0, 5.0}
+    # v[1..3]  # => Vector{2.0, 3.0, 4.0}
+    # v[6..10] # raise IndexError
+    # ```
+    def [](range : Range(Int, Int))
+      self[*Indexable.range_to_index_and_count(range, size)]
+    end
+
+    # Returns count or less (if there aren't enough) elements starting at the
+    # given start index.
+    #
+    # Negative indices count backward from the end of the vector (-1 is the last
+    # element). Additionally, an empty vector is returned when the starting index
+    # for an element range is at the end of the vector.
+    #
+    # Raises `IndexError` if the starting index is out of range.
+    #
+    # ```
+    # v = Vector{1.0, 2.0, 3.0, 4.0, 5.0}
+    # v[1, 3]  # => Vector{2.0, 3.0, 4.0}
+    # v[6, 10] # raise IndexError
+    # ```
+    def [](start : Int, count : Int)
+      raise ArgumentError.new "Negative count: #{count}" if count < 0
+
+      if start == size
+        return Vector.new
+      end
+
+      start += size if start < 0
+      raise IndexError.new unless 0 <= start <= size
+
+      if count == 0
+        return Vector.new
+      end
+
+      count = Math.min(count, size - start)
+
+      Vector.build(count) do |buffer|
+        buffer.copy_from(@buffer + start, count)
+        count
+      end
+    end
+
+    # Sets the given value at the given index.
+    #
+    # Negative indices can be used to start counting from the end of the array.
+    # Raises `IndexError` if trying to set an element outside the array's range.
+    #
+    # ```
+    # vec = Vector{1.0, 2.0, 3.0}
+    # vec[0] = 5.0
+    # p vec # => Vec{5.0, 2.0, 3.0}
+    #
+    # vec[3] = 5.0 # raises IndexError
+    # ```
+    @[AlwaysInline]
+    def []=(index : Int, value)
+      index = check_index_out_of_bounds index
+      @buffer[index] = value.to_f64
+    end
+
+    # Replaces a subrange with a single value. All elements in the range
+    # `index...index+count` are removed and replaced by a single element
+    # *value*.
+    #
+    # If *count* is zero, *value* is inserted at *index*.
+    #
+    # Negative values of *index* count from the end of the vector.
+    #
+    # ```
+    # vec = Vector{1.0, 2.0, 3.0, 4.0, 5.0}
+    # vec[1, 3] = 6.0
+    # vec # => Vector{1.0, 6.0, 5.0}
+    #
+    # vec = Vector{1.0, 2.0, 3.0, 4.0, 5.0}
+    # vec[1, 0] = 6.0
+    # vec # => Vector{1.0, 6.0, 2.0, 3.0, 4.0, 5.0}
+    # ```
+    def []=(index : Int, count : Int, value)
+      raise ArgumentError.new "Negative count: #{count}" if count < 0
+
+      value = value.to_f64
+      index = check_index_out_of_bounds index
+      count = (index + count <= size) ? count : size - index
+
+      case count
+      when 0
+        insert index, value
+      when 1
+        @buffer[index] = value
+      else
+        diff = count - 1
+        (@buffer + index + 1).move_from(@buffer + index + count, size - index - count)
+        (@buffer + @size - diff).clear(diff)
+        @buffer[index] = value
+        @size -= diff
+      end
+
+      value
+    end
+
+    # Replaces a subrange with a single value.
+    #
+    # ```
+    # vec = Vector{1.0, 2.0, 3.0, 4.0, 5.0}
+    # vec[1..3] = 6.0
+    # vec # = Vector{1.0, 6.0, 5.0}
+    #
+    # vec = Vector{1.0, 2.0, 3.0, 4.0, 5.0}
+    # vec[1...1] = 6
+    # vec # = Vector{1.0, 6.0, 2.0, 3.0, 4.0, 5.0}
+    # ```
+    def []=(range : Range(Int, Int), value)
+      self[*Indexable.range_to_index_and_count(range, size)] = value.to_f64
+    end
+
+    # Replaces a subrange with the elements of the given vector.
+    #
+    # ```
+    # vec = Vector{1.0, 2.0, 3.0, 4.0, 5.0}
+    # vec[1, 3] = Vector{6.0, 7.0, 8.0}
+    # vec # = Vector{1.0, 6.0, 7.0, 8.0, 5.0}
+    #
+    # vec = Vector{1.0, 2.0, 3.0, 4.0, 5.0}
+    # vec[1, 3] = Vector{6.0, 7.0}
+    # vec # = Vector{1.0, 6.0, 7.0, 5.0}
+    # ```
+    def []=(index : Int, count : Int, values : Vector)
+      raise ArgumentError.new "Negative count: #{count}" if count < 0
+
+      index = check_index_out_of_bounds index
+      count = (index + count <= size) ? count : size - index
+      diff = values.size - count
+
+      if diff == 0
+        # Replace values directly
+        (@buffer + index).copy_from(values.to_unsafe, values.size)
+      elsif diff < 0
+        # Need to shrink
+        diff = -diff
+        (@buffer + index).copy_from(values.to_unsafe, values.size)
+        (@buffer + index + values.size).move_from(@buffer + index + count, size - index - count)
+        (@buffer + @size - diff).clear(diff)
+        @size -= diff
+      else
+        # Need to grow
+        resize_to_capacity(Math.pw2ceil(@size + diff))
+        (@buffer + index + values.size).move_from(@buffer + index + count, size - index - count)
+        (@buffer + index).copy_from(values.to_unsafe, values.size)
+        @size += diff
+      end
+
+      values
+    end
+
+    # Combined comparison operator. Returns *0* if `self` equals *other*, *1* if
+    # `self` is greater than *other* and *-1* if `self` is smaller than *other*.
+    #
+    # It compares the elements of both vectors in the same position using the
+    # `<=>` operator.  As soon as one of such comparisons returns a non-zero
+    # value, that result is the return value of the comparison.
+    #
+    # If all elements are equal, the comparison is based on the size of the vectors.
+    #
+    # ```
+    # Vector{8.0} <=> Vector(1.0, 2.0, 3.0} # => 1
+    # Vector{2.0} <=> Vector{4.0, 2.0, 3.0} # => -1
+    # Vector{1.0, 2.0} <=> Vector{1.0, 2.0} # => 0
+    # ```
+    def <=>(other)
+      min_size = Math.min(size, other.size)
+      0.upto(min_size - 1) do |i|
+        n = @buffer[i] <=> other.to_unsafe[i]
+        return n if n != 0
+      end
+      size <=> other.size
+    end
+
+    # Alias for `#push`
+    def <<(value : Float64)
+      push(value)
+    end
+
+    def ==(other : Vector)
+      equals?(other) { |x, y| x == y }
+    end
+
+    # :nodoc:
+    def ==(other)
+      false
+    end
+
+    def +(other)
+      add(other)
+    end
+
+    def -(other)
+      subtract(other)
+    end
+
+    def *(other)
+      multiply(other)
+    end
+
+    def clone
+      Vector.create(@elements.clone)
+    end
+
+    # Invokes the given block for each element of `self`.
+    def map(&block)
+      Vector.new(size) { |i| yield @buffer[i] }
+    end
+
+    # Returns true if all vectors are linearly independent.
+    def independent?(*vs)
+      vs.each do |v|
+        raise "Dimension mismatch. Vectors not all the same size." unless v.size == vs.first.size
+      end
+      return false if vs.size > sv.first.size
+      Matrix[*vs].rank.equal?(vs.count)
+    end
+
+
+
+    # Invokes the given block for each element of `self`, replacing the element
+    # with the value returned by the block. Returns `self`.
+    #
+    # ```
+    # vec = Vector{1.0, 2.0, 3.0}
+    # vec.map! { |x| x * x }
+    # a # => Vector{1.0, 4.0, 9.0}
+    # ```
+    def map!
+      @buffer.map!(size) { |e| yield e }
+      self
+    end
+
+    # Optimized version of `Enumerable#map_with_index`.
+    def map_with_index(&block)
+      Vector.new(size) { |i| yield @buffer[i], i }
+    end
+
+    # Like `map_with_index`, but mutates `self` instead of allocating a new object.
+    def map_with_index!(&block)
+      to_unsafe.map_with_index!(size) { |e, i| yield e, i }
+      self
+    end
+
+    # Returns the magnitude/euclidian norm of this vector.
+    #
+    # [https://en.wikipedia.org/wiki/Euclidean_distance](https://en.wikipedia.org/wiki/Euclidean_distance)
+    def magnitude
+      sum = reduce(0.0) { |acc, e| acc += e * e }
+      Math.sqrt(sum)
+    end
+
+    # Returns the `ith` element of the vector. Returns `nil` if `i`
+    # is out of bounds. Indexing starts from 1.
+    def e(i)
+      (i < 1 || i > @size) ? nil : self[i - 1]
+    end
+
+    # Returns a new vector created by normalizing this one
+    # to have a magnitude of `1`. If the vector is a zero
+    # vector, it will not be modified.
+    def to_unit_vector
+      r = magnitude
+      r == 0 ? dup : map { |x| x.to_f64 / r }
+    end
+
+    # Returns the angle between this vector and another in radians.
+    # If the vectors are mirrored across their axes this will return `nil`.
+    def angle_from(vector)
+      v = vector.is_a?(Vector) ? vector : Vector.create(vector)
+
+      unless size == v.size
+        raise "Cannot compute the angle between vectors with different dimensionality"
+      end
+
+      dot = 0_f64
+      mod1 = 0_f64
+      mod2 = 0_f64
+
+      zip(vector).each do |x, v|
+        dot += x * v
+        mod1 += x * x
+        mod2 += v * v
+      end
+
+      mod1 = Math.sqrt(mod1)
+      mod2 = Math.sqrt(mod2)
+
+      if mod2 * mod2 == 0
+        return 0.0
+      end
+
+      theta = (dot / (mod1 * mod2)).clamp(-1, 1)
+      Math.acos(theta)
+    end
+
+    # Returns whether the vectors are parallel to each other.
+    def parallel_to?(vector)
+      angle = angle_from(vector)
+      angle <= Apatite.precision
+    end
+
+    # Returns whether the vectors are antiparallel to each other.
+    def antiparallel_to?(vector)
+      angle = angle_from(vector)
+      (angle - Math::PI).abs <= Apatite.precision
+    end
+
+    # Returns whether the vectors are perpendicular to each other.
+    def perpendicular_to?(vector)
+      (dot(vector)).abs <= Apatite.precision
+    end
+
+    # When the input is a number, this returns the result of adding
+    # it to all vector elements. When it's a vector, the vectors
+    # will be added together.
+    def add(value)
+      run_binary_op(value) { |a, b| a + b }
+    end
+
+    # When the input is a number, this returns the result of subtracting
+    # it to all vector elements. When it's a vector, the vectors
+    # will be subtracted.
+    def subtract(value)
+      run_binary_op(value) { |a, b| a - b }
+    end
+
+    # When the input is a number, this returns the result of multiplying
+    # it to all vector elements. When it's a vector, the vectors
+    # will be element-wise multiplied.
+    def multiply(value)
+      run_binary_op(value) { |a, b| a * b }
+    end
+
+    # Returns the sum of all elements in the vector.
+    def sum
+      reduce(0) { |acc, i| acc + i }
+    end
+
+    # Returns the cross product of this vector with the others.
+    #
+    # ```
+    # v1 = Vector{1.0, 0.0, 0.0}
+    # v2 = Vector{0.0, 1.0, 0.0}
+    # v1.cross(v2) => Vector{0.0, 0.0, 1.0}
+    # ```
+    def cross(*vs)
+      raise "cross product is not defined on vectors of dimension #{size}" unless size >= 2
+      raise ArgumentError.new("wrong number of arguments (#{vs.size} for #{size - 2})") unless vs.size == size - 2
+      vs.each do |v|
+        raise "Dimension mismatch. Vectors not all the same size." unless v.size == size
+      end
+      case size
+      when 2
+        Vector[-@buffer[1], @buffer[0]]
+      when 3
+        v = vs[0]
+        Vector[ v[2] * @buffer[1] - v[1] * @buffer[2],
+          v[0] * @buffer[2] - v[2] * @buffer[0],
+          v[1] * @buffer[0] - v[0] * @buffer[1] ]
+      else
+        # TODO
+        # rows = [self, *vs, Vector.new(size) {|i| Vector.basis(size: size, index: i) }]
+        # Matrix.rows(rows).laplace_expansion(row: size - 1)
+      end
+    end
+
+    # Returns a new vector with the first `n` elements removed from
+    # the beginning.
+    def chomp(n)
+      elements = [] of Float64
+      each_with_index { |e, i| elements << e if i >= n }
+      Vector.create(elements)
+    end
+
+    # Returns a vector containing only the first `n` elements.
+    def top(n)
+      elements = [] of Float64
+      each_with_index { |e, i| elements << e if i < n }
+      Vector.create(elements)
+    end
+
+    # Returns a new vector with the provided `elements` concatenated
+    # on the end.
+    def augment(elements)
+      elements = elements.is_a?(Vector) ? elements : Vector.create(elements)
+      concat(elements)
+    end
+
+    # Return a new `Vector` with the log of every item in `self`.
+    def log
+      map { |x| Math.log(x) }
+    end
+
+    # Get the product of all elements in this vector.
+    def product
+      reduce { |acc, v| acc *= v }
+    end
+
+    # Get the scalar (dot) product of this vector with `vector`.
+    #
+    # [https://en.wikipedia.org/wiki/Scalar_product](https://en.wikipedia.org/wiki/Scalar_product)
+    def dot(other)
+      other = other.is_a?(Vector) ? other : Vector.create(other)
+      unless size == other.size
+        raise "Cannot compute the dot product of vectors with different dimensionality"
+      end
+
+      product = 0
+      (0...size).each do |i|
+        product += self[i] * other[i]
+      end
+
+      product
+    end
+
+    # Returns the (absolute) largest element in this vector.
+    def max
+      reduce { |acc, i| i.abs > acc.abs ? i : acc }
+    end
+
+    # Gets the index of the largest element in this vector.
+    def max_index
+      idx = 0
+      each_with_index { |e, i| idx = e.abs > self[idx].abs ? i : idx }
+      idx
+    end
+
+    # Creates a single-row matrix from this vector.
+    def covector
+      Matrix.row_vector(self)
+    end
+
+    # Returns a diagonal `Matrix` with the vectors elements as its
+    # diagonal elements.
+    def to_diagonal_matrix
+      Matrix::Diagonal.new(elements)
+    end
+
+    # Gets the result of rounding the elements of the vector.
+    def round
+      map { |x| x.round }
+    end
+
+    # Transpose this vector into a 1xn `Matrix`
+    def transpose
+      Matrix.column_vector(self)
+    end
+
+    # Returns a copy of the vector with elements set to `value` if
+    # they differ from it by less than `Apetite.precision`
+    def snap_to(value)
+      map { |y| (y - x).abs <= Apetite.precision ? value : y }
+    end
+
+    # Gets this vector's distance from the argument, when considered
+    # a point in space.
+    def distance_from(obj)
+      if object.is_a?(Plane) || object.is_a?(Line)
+        return object.distance_from(self)
+      end
+
+      v = elements.is_a?(Vector) ? elements.elements : elements
+      unless v.size == @elements.size
+        return nil
+      end
+
+      sum = 0
+      part = 0
+      each_with_index do |x, i|
+        part = x - v[i - 1]
+        sum += part * part
+      end
+      Math.sqrt(sum)
+    end
+
+    # Returns true if the vector is a point on the given line
+    def lies_on(line)
+      line.contains(self)
+    end
+
+    # Returns true if the vector is a point on the given plane.
+    def lies_in(plane)
+      plane.contains(self)
+    end
+
+    def push(value : Float64)
+      check_needs_resize
+      @buffer[@size] = value
+      @size += 1
+      self
+    end
+
+    def push(*values : Float64)
+      new_size = @size + values.size
+      resize_to_capacity(Math.pw2ceil(new_size)) if new_size > @capacity
+      values.each_with_index do |value, i|
+        @buffer[@size + i] = value
+      end
+      @size = new_size
+      self
+    end
+
+    # Rotates the vector about the given `object`. The object should
+    # be a point if the vector is 2D, and a line if it is 3D. Be
+    # careful with line directions!
+    def rotate(t, object)
+      # TODO
+    end
+
+    # Returns the result of reflecting the point in the given `object`
+    # (point, line, or plane).
+    def reflection_in(object)
+      if object.is_a?(Plane) || object.is_a?(Line)
+        # object is a plane or line
+        p = @elements.dup
+        c = object.point_closest_to(p).elements
+        return Vector.create([
+          C[0] + (C[0] - P[0]),
+          C[1] + (C[1] - P[1]),
+          C[2] + (C[2] - (P[2] || 0)),
+        ])
+      end
+
+      # object is a point
+      q = object.is_a?(Vector) ? object.elements : object
+      unless @elements.size == q.size
+        return nil
+      end
+
+      map_with_index { |x, i| q[i - 1] + (q[i - 1] - x) }
+    end
+
+    # Utility to make sure vectors are 3D. If they are 2D, a zero
+    # z-component is added.
+    def to_3d
+      v = dup
+      case v.elements.size
+      when 3
+        break
+      when 2
+        v.elements.push(0)
+      else
+        return nil
+      end
+      return v
+    end
+
+    def set_elements(elements)
+      @elements = elements.is_a?(Vector) ? elements.elements : elements
+      self
+    end
+
+    def concat(other : Vector)
+      other_size = other.size
+      new_size = size + other_size
+      if new_size > @capacity
+        resize_to_capacity(Math.pw2ceil(new_size))
+      end
+
+      (@buffer + @size).copy_from(other.to_unsafe, other_size)
+      @size = new_size
+
+      self
+    end
+
+    def to_a
+      Array(Float64).new(size) { |i| self[i] }
+    end
+
+    def to_s(io)
+      io << "Vector{"
+      join ", ", io, &.inspect(io)
+      io << "}"
+    end
+
+    def pretty_print(pp) : Nil
+      pp.list("Vector{", self, "}")
+    end
+
+    @[AlwaysInline]
+    def unsafe_fetch(index : Int)
+      @buffer[index]
+    end
+
+    # Returns a pointer to the internal buffer where `self`'s elements are stored.
+    #
+    # This method is **unsafe** because it returns a pointer, and the pointed might eventually
+    # not be that of `self` if the array grows and its internal buffer is reallocated.
+    #
+    # ```
+    # vec = Vector{1.0, 2.0, 3.0}
+    # vec.to_unsafe[0] # => 1.0
+    # ```
+    def to_unsafe : Pointer(Float64)
+      @buffer
+    end
+
+    # Removes all elements from self.
+    #
+    # ```
+    # vec = Vector{1.0, 2.0, 3.0}
+    # vec.clear # => Vector{}
+    # ```
+    def clear
+      @buffer.clear(@size)
+      @size = 0
+      self
+    end
+
+    def zero?
+      all?(&.zero?)
+    end
+
+    # :nodoc:
+    protected def size=(size : Int)
+      @size = size.to_i
+    end
+
+    private def check_needs_resize
+      double_capacity if @size == @capacity
+    end
+
+    private def double_capacity
+      resize_to_capacity(@capacity == 0 ? 3 : (@capacity * 2))
+    end
+
+    private def resize_to_capacity(capacity)
+      @capacity = capacity
+      if @buffer
+        @buffer = @buffer.realloc(@capacity)
+      else
+        @buffer = Pointer(Float64).malloc(@capacity)
+      end
+    end
+
+    # Call a block on the value
+    private def run_binary_op(value, &block : (T, T) -> T)
+      if value.is_a?(Number)
+        return map { |v| yield(v, value) }
+      end
+
+      values = value.is_a?(Vector) ? value.elements : value
+
+      unless @elements.size == values.size
+        raise "Cannot perform operations on vectors with different dimensions."
+      end
+
+      map_with_index { |x, i| yield(x, values[i]) }
+    end
+  end
+end
diff --git a/src/apatite/narray.cr b/src/apatite/narray.cr
deleted file mode 100644
index 204a45a..0000000
--- a/src/apatite/narray.cr
+++ /dev/null
@@ -1,255 +0,0 @@
-module Apatite
-  class NArray(T) < Array(T)
-    include Indexable(T)
-    include Comparable(NArray)
-
-    getter size : Int32
-    getter shape : Array(Int32)
-    getter ndim : Int32
-
-    def initialize(@selection : T, @dtype = Int32)
-      @size = @selection.size
-      @shape = @selection.shape
-      @ndim = @shape.size
-      super(1, @selection)
-    end
-
-    # def self.new(shape : Array(Int32), &block : Int32 -> T)
-    # end
-
-    # def ==(other : NArray)
-    # end
-
-    # # :nodoc:
-    # def ==(other)
-    #   false
-    # end
-
-    # def <=>(other : NArray)
-    # end
-
-    # def &(other : NArray(U)) forall U
-    # end
-
-    # def |(other : NArray(U)) forall U
-    # end
-
-    # def +(other : NArray(U)) forall U
-    # end
-
-    # def -(other : NArray(U)) forall U
-    # end
-
-    # def *(times : Int)
-    # end
-
-    # def <<(value : T)
-    # end
-
-    # @[AlwaysInline]
-    # def []=(index : Int, value : T)
-    # end
-
-    # def []=(index : Int, count : Int, value : T)
-    # end
-
-    # def []=(range : Range(Int, Int), value : T)
-    # end
-
-    # def []=(index : Int, count : Int, values : NArray(T))
-    # end
-
-    # def []=(range : Range(Int, Int), values : NArray(T))
-    # end
-
-    # def [](range : Range(Int, Int))
-    # end
-
-    # def [](start : Int, count : Int)
-    # end
-
-    def all?
-      flatten.all?
-    end
-
-    def any?
-      flatten.any?
-    end
-
-    # def arg_max
-    # end
-
-    # def arg_min
-    # end
-
-    # def arg_partition
-    # end
-
-    # def arg_sort
-    # end
-
-    # def as_type
-    # end
-
-    # def byte_swap(inplace = false)
-    # end
-
-    # def choose
-    # end
-
-    # def clip
-    # end
-
-    # def compress
-    # end
-
-    # def conj
-    # end
-
-    # def conjugate
-    # end
-
-    # def copy
-    # end
-
-    # def cum_prod
-    # end
-
-    # def cum_sum
-    # end
-
-    # def diagonal
-    # end
-
-    # def dot
-    # end
-
-    # def dump
-    # end
-
-    # def dumps
-    # end
-
-    # def fill
-    # end
-
-    # def flatten
-    # end
-
-    # def get_field
-    # end
-
-    # def item
-    # end
-
-    # def item_set
-    # end
-
-    # def max
-    # end
-
-    # def mean
-    # end
-
-    # def min
-    # end
-
-    # def new_byte_order
-    # end
-
-    # def non_zero
-    # end
-
-    # def partition
-    # end
-
-    # def prod
-    # end
-
-    # def ptp
-    # end
-
-    # def put
-    # end
-
-    # def ravel
-    # end
-
-    # def repeat
-    # end
-
-    # def repeat
-    # end
-
-    # def reshape
-    # end
-
-    # def resize
-    # end
-
-    # def round
-    # end
-
-    # def search_sorted
-    # end
-
-    # def set_field
-    # end
-
-    # def set_flags
-    # end
-
-    # def sort
-    # end
-
-    # def squeeze
-    # end
-
-    # def std
-    # end
-
-    # def sum
-    # end
-
-    # def swap_axes
-    # end
-
-    # def take
-    # end
-
-    # def to_bytes
-    # end
-
-    # def to_a
-    # end
-
-    # def to_json
-    # end
-
-    # def to_s
-    # end
-
-    # def to_unsafe
-    # end
-
-    # def to_yaml
-    # end
-
-    # def trace
-    # end
-
-    # def transpose
-    # end
-
-    # @[AlwaysInline]
-    # def unsafe_fetch(index : Int)
-    #   @buffer[index]
-    # end
-
-    # def var
-    # end
-
-    # def view
-    # end
-  end
-end