diff --git a/README.md b/README.md
new file mode 100644
index 0000000..3c11ae6
--- /dev/null
+++ b/README.md
@@ -0,0 +1,49 @@
+# level-index-simulator
+
+# `sli` - A custom precision MATLAB symmetric level-index arithmetic simulator.
+
+## About
+
+`sli` is a MATLAB toolbox that provides sli objects and operations on them.
+
+* [sli.m](sli.m) - the main object file.
+* [tests](./tests) - various testing scripts that test the accuracy of representation and operations.
+* [experiments](./experiments) - scripts to reproduce the experiments with level-index arithmetic in [1].
+
+[![Open in MATLAB Online](https://www.mathworks.com/images/responsive/global/open-in-matlab-online.svg)](https://matlab.mathworks.com/open/github/v1?repo=north-numerical-computing/level-index-simulator)
+
+## Quick start
+
+```
+x = sli(2,12);
+x = x.set_val(pi)
+x = 
+
+  sli with properties:
+
+    level_bits: 2
+    index_bits: 12
+          sign: 0
+    reciprocal: 1
+         level: 2
+         index: 0.135253906250000
+         value: 3.141899100868418
+```
+
+### Installation
+
+Clone the repository and add the root directory to the MATLAB search path.
+
+## Requirements
+
+The code was developed on MATLAB 2023b. Experiments make use of the [CPFloat](https://github.com/north-numerical-computing/cpfloat) package.
+
+## References
+
+[1] Mantas Mikaitis. [*MATLAB Simulator of Level-Index Arithmetic*](https://). In Prep., Feb. 2024.
+
+## Licence
+
+The code is distributed under the terms of the BSD 2-Clause License;
+see [license.txt](license.txt).
+BBB
\ No newline at end of file
diff --git a/experiments/matvec.m b/experiments/matvec.m
new file mode 100644
index 0000000..c593119
--- /dev/null
+++ b/experiments/matvec.m
@@ -0,0 +1,67 @@
+% Compare matrix-vector product accuracy in floating-point and
+% level-index arithmetics.
+%
+% References
+%
+% M. Mikaitis, MATLAB Simulator of Level Index Arithmetic, In Prep., 2024.
+
+clear all;
+
+nlist = round(logspace(1,4,20));
+
+% Set up cpfloat options
+options.format = 'fp16';
+options.round = 1;
+options.explim = 1;
+cpfloat([], options);
+k = 0;
+
+for n = nlist
+    k = k + 1;
+    fprintf('k = %2d, n = %7d\n',k,n);
+    A = rand(10,n)*100;
+    x = rand(n,1);
+
+    % Set up SLI vectors
+    Asli = zeros(10, n, 2, 12, 'sli');
+    xsli = zeros(n, 1, 2, 12, 'sli');
+
+    % Convert vectors to SLI
+    for j = 1:numel(A)
+       Asli(j) = Asli(j).set_val(A(j));
+    end
+    for j = 1:numel(x)
+       xsli(j) = xsli(j).set_val(x(j));
+    end
+
+    % Double reference result
+    ye = A*x;
+    absAx = abs(A)*abs(x);
+
+    % Binary16 result
+    y1 = mtimesv(A,x,options);
+    berr1(k) = max(abs(ye-y1)./absAx);
+
+    % Level index 2.12
+    y2 = Asli*xsli;
+    berr2(k) = max(abs(ye-[y2.value]')./absAx);
+end
+
+filename = sprintf('matvec_binary16_1.dat');
+fid = fopen(filename, 'w');
+for i=1:length(nlist)
+    fprintf(fid, "%d %f %f\n", nlist(i), berr1(i), berr2(i));
+end
+fclose(fid);
+
+loglog(nlist, berr1, '-*', ...
+      nlist, berr2, '-o');
+legend('binary16', 'level index 2.12');
+
+function y = mtimesv(A,x,options)
+    [m,n] = size(A);
+    y = zeros(m,1);
+    for i=1:n
+        y = cpfloat(y + cpfloat(A(:,i)*x(i), options), options);
+    end
+end
\ No newline at end of file
diff --git a/license.txt b/license.txt
new file mode 100644
index 0000000..d9cb936
--- /dev/null
+++ b/license.txt
@@ -0,0 +1,25 @@
+BSD 2-Clause License
+
+Copyright (c) 2024, Mantas Mikaitis
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are met:
+
+1. Redistributions of source code must retain the above copyright notice, this
+   list of conditions and the following disclaimer.
+
+2. Redistributions in binary form must reproduce the above copyright notice,
+   this list of conditions and the following disclaimer in the documentation
+   and/or other materials provided with the distribution.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
+DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
+FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
+SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
+CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
+OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
diff --git a/sli.m b/sli.m
new file mode 100644
index 0000000..72a584c
--- /dev/null
+++ b/sli.m
@@ -0,0 +1,642 @@
+%SLI Symmetric level-index object.
+%
+% SLI object defines a custom precision symmetric level index number and
+% various operations on it.
+%
+% References
+%
+% C. W. Clenshaw and F. W. J. Olver, Beyond floating point,
+% Journal of the ACM, 31 (1984), p. 319–328.
+%
+% M. Mikaitis, MATLAB Simulator of Level Index Arithmetic, In Prep., 2024.
+
+classdef sli
+
+    properties (SetAccess = protected)
+        level_bits = 2
+        index_bits = 12
+        sign
+        reciprocal
+        level
+        index
+        value
+    end
+
+    methods
+        function obj = sli(lbits, ibits)
+            if nargin > 0
+                obj.level_bits = lbits;
+                obj.index_bits = ibits;
+            end
+        end
+
+        function obj = set_val(obj, val)
+            obj.sign = val < 0;
+            temp = abs(val);
+            if (temp == 0)
+                obj.level = 1;
+                obj.index = 0;
+                obj.reciprocal = 0;
+                obj = update_double(obj);
+                return;
+            end
+            if (temp < 1)
+                temp = 1/temp;
+                obj.reciprocal = 0;
+            else
+                obj.reciprocal = 1;
+            end
+            obj.level = 0;
+            while (temp >= 1)
+                obj.level = obj.level + 1;
+                temp = log(temp);
+            end
+            obj.index = temp;
+            obj = update_double(obj);
+        end
+
+        function obj = set_sli(obj, s, r, level, index)
+            obj.sign = s;
+            obj.reciprocal = r;
+            if (level == 0)
+                error('Zero level is not supported.');
+            end
+            if level <= 2^obj.level_bits
+                obj.level = level;
+            else
+                obj.level = 2^obj.level_bits;
+                warning('Specified level is too big.');
+            end
+            if (index < 1)
+                obj.index = floor(2^obj.index_bits * ...
+                    index)/2^(obj.index_bits);
+            end
+            obj = update_double(obj);
+        end
+
+        function r = iszero(obj)
+            indices = find([obj.reciprocal] == 0 & ...
+                [obj.index] == 0 & [obj.level] == 1);
+            if ~isempty(indices)
+                r(indices) = deal(1);
+            else
+                r(1:numel(obj)) = deal(0);
+            end
+        end
+
+        function r = abs(a)
+            [a.sign] = deal(0);
+            r = update_double(a);
+        end
+
+        function r = nextafter(obj)
+            obj.index = obj.index + 2^-obj.index_bits;
+            if obj.index == 1
+                obj.index = 0;
+                obj.level = obj.level + 1;
+            end
+            r = update_double(obj);
+        end
+
+        function r = plus(a, b)
+            [a, b] = check_compatible_sizes(a,b);
+            r = a;
+            for (k = 1:numel(a))
+                if iszero(a(k))
+                    r(k) = b(k); return;
+                elseif iszero(b(k))
+                    r(k) = a(k); return;
+                end
+                if (abs(a(k)) < abs(b(k)))
+                    temp = a(k);
+                    a(k) = b(k);
+                    b(k) = temp;
+                end
+                if (a(k).sign ~= b(k).sign)
+                    r(k) = minus(abs(a(k)),abs(b(k)));
+                    r(k).sign = a(k).sign;
+                    r(k) = update_double(r(k));
+                    continue;
+                end
+                r(k) = main_sli_algorithm(a(k), b(k), 'add');
+            end
+        end
+
+        function r = minus(a, b)
+            [a, b] = check_compatible_sizes(a,b);
+            r = a;
+            for (k = 1:numel(a))
+                if (a(k) == b(k))
+                    r(k) = a(k).set_val(0);
+                    continue;
+                end
+                if iszero(a(k))
+                    r(k) = b(k);
+                    r(k).sign = ~b(k).sign;
+                    r(k) = update_double(r(k));
+                    continue;
+                elseif iszero(b(k))
+                    r(k) = a(k);
+                    continue;
+                end
+                if (a(k).sign ~= b(k).sign)
+                    r(k) = plus(abs(a(k)),abs(b(k)));
+                    r(k).sign = a(k).sign;
+                    r(k) = update_double(r(k));
+                    continue;
+                end
+                if (abs(a(k)) < abs(b(k)))
+                    temp = a(k);
+                    a(k) = b(k);
+                    b(k) = temp;
+                    a(k).sign = ~a(k).sign;
+                end
+                r(k) = main_sli_algorithm(a(k), b(k), 'sub');
+            end
+        end
+
+        function r = uminus(a)
+            a.sign = ~a.sign;
+            r = a;
+        end
+
+        function r = uplus(a)
+            r = a;
+        end
+
+        function r = times(a, b)
+            [a, b] = check_compatible_sizes(a,b);
+            r = a;
+            for (k = 1:numel(a))
+                if iszero(a(k))
+                    r(k) = a(k);
+                    continue;
+                elseif iszero(b(k))
+                    r(k) = b(k);
+                    continue;
+                end
+
+                if (li_le(abs(a(k)), abs(b(k))))
+                    temp = a(k);
+                    a(k) = b(k);
+                    b(k) = temp;
+                end
+                a(k).level = a(k).level - 1;
+                b(k).level = b(k).level - 1;
+                if (a(k).reciprocal == b(k).reciprocal)
+                    r(k) = li_add(abs(a(k)), abs(b(k)));
+                    r(k).reciprocal = a(k).reciprocal;
+                else
+                    r(k) = li_sub(abs(a(k)), abs(b(k)));
+                    r(k).reciprocal = a(k).reciprocal;
+                end
+                r(k).level = r(k).level + 1;
+                if (a(k).sign ~= b(k).sign)
+                    r(k).sign = 1;
+                else
+                    r(k).sign = 0;
+                end
+                if (iszero(r(k)))
+                    r(k).reciprocal = 1;
+                end
+                r(k) = update_double(r(k));
+
+            end
+        end
+
+        function R = mtimes(A, B)
+            lbits = A(1).level_bits;
+            ibits = A(1).index_bits;
+            [m, k1] = size(A);
+            [k2, n] = size(B);
+            R = zeros(m, n, lbits, ibits, 'sli');
+
+            if (k1 ~= k2)
+                error('Incompatible matrix or vector dimensions.');
+            end
+
+            for i = 1:m
+                for j = 1:n
+                    for k = 1:k1
+                        R(i,j) = R(i,j) + A(i, k).*B(k, j);
+                    end
+                end
+            end
+        end
+
+        function r = rdivide(a, b)
+            [a, b] = check_compatible_sizes(a,b);
+            r = a;
+            for (k = 1:numel(a))
+                if iszero(a(k))
+                    r = a(k); continue;
+                elseif iszero(b(k))
+                    error('Division by zero');
+                end
+                invert_recipr = 0;
+                if (li_le(abs(a(k)), abs(b(k))))
+                    temp = a(k);
+                    a(k) = b(k);
+                    b(k) = temp;
+                    invert_recipr = 1;
+                end
+                a(k).level = a(k).level - 1;
+                b(k).level = b(k).level - 1;
+                if (a(k).reciprocal == b(k).reciprocal)
+                    r(k) = li_sub(abs(a(k)), abs(b(k)));
+                    r(k).reciprocal = a(k).reciprocal;
+                else
+                    r(k) = li_add(abs(a(k)), abs(b(k)));
+                    r(k).reciprocal = a(k).reciprocal;
+                end
+                if invert_recipr
+                    r(k).reciprocal = ~r(k).reciprocal;
+                end
+                r(k).level = r(k).level + 1;
+                if (a(k).sign ~= b(k).sign)
+                    r(k).sign = 1;
+                else
+                    r(k).sign = 0;
+                end
+                if (iszero(r(k)))
+                    r(k).reciprocal = 1;
+                end
+                r(k) = update_double(r(k));
+            end
+        end
+
+        function r = ldivide(a, b)
+            r = rdivide(b, a);
+        end
+
+        function r = lt(a, b)
+            for k = 1:numel(a)
+                if (a(k).sign > b(k).sign)
+                    r(k) = true;
+                elseif (a(k).sign < b(k).sign)
+                    r(k) = false;
+                elseif (a(k).reciprocal < b(k).reciprocal)
+                    r(k) = true;
+                elseif (a(k).reciprocal > b(k).reciprocal)
+                    r(k) = false;
+                elseif (a(k).level < b(k).level)
+                    r(k) = a(k).reciprocal;
+                elseif (a(k).level > b(k).level)
+                    r(k) = ~a(k).reciprocal;
+                elseif (a(k).reciprocal)
+                    r(k) = a(k).index < b(k).index;
+                else
+                    r(k) = a(k).index > b(k).index;
+                end
+            end
+        end
+
+        function r = gt(a, b)
+            r = lt(b, a);
+        end
+
+        function r = le(a, b)
+            r = lt(a, b) | eq(a, b);
+        end
+
+        function r = ge(a, b)
+            r = gt(a, b) | eq(a, b);
+        end
+
+        function r = eq(a, b)
+            r = (a.sign == b.sign) & (a.reciprocal == b.reciprocal) & ...
+                (a.level == b.level) & (a.index == b.index);
+        end
+
+        function r = ne(a, b)
+            r = ~eq(a,b);
+        end
+    end
+
+    methods (Static)
+
+        % Create SLI arrays of zeros.
+        % Use:
+        %     zeros(m, lbits, ibits, 'sli')
+        %     zeros(m,k, lbits, ibits, 'sli')
+        function z = zeros(varargin)
+            if (nargin == 3)
+                for i = 1:varargin{1}
+                    z(i) = sli(varargin{2}, varargin{3});
+                    z(i).set_val(0);
+                end
+            elseif (nargin == 4)
+                for i = 1:varargin{1}
+                    for j = 1:varargin{2}
+                        z(i,j) = sli(varargin{3}, varargin{4});
+                        z(i,j) = z(i,j).set_val(0);
+                    end
+                end
+            else
+                error("Please specify 3 or 4 input arguments.");
+            end
+        end
+    end
+
+    methods (Access = private)
+        % Check vector sizes are compatible, as per
+        % https://uk.mathworks.com/help/matlab/matlab_prog/compatible-array-sizes-for-basic-operations.html
+        function [a, b] = check_compatible_sizes(a, b)
+              [m, n] = size(a); ael = numel(a);
+              [k, l] = size(b); bel = numel(b);
+
+              if (m == k && n == l)
+                  return
+              elseif (ael == 1)
+                  temp = a;
+                  a = b;
+                  a(:, :) = temp;
+              elseif (bel == 1)
+                  temp = b;
+                  b = a;
+                  b(:,:) = temp;
+              elseif (m==k && n == 1)
+                  a = repmat(a, 1, l);
+              elseif (m==k && l == 1)
+                  b = repmat(b, 1, n);
+              elseif (n==1 && k==1)
+                  a = repmat(a, 1, l);
+                  b = repmat(b, m, 1);
+              elseif (m==1 && l==1)
+                  a = repmat(a, k, 1);
+                  b = repmat(b, 1, n);
+              else
+                  error('Arrays are not of compatibles sizes.');
+              end
+        end
+
+        function obj = update_double(obj)
+            for s = 1:numel(obj)
+                obj(s) = update_sli(obj(s));
+                if (obj(s).index == 0 && obj(s).level == 1 && obj(s).reciprocal == 0)
+                    obj(s).value = 0;
+                    continue;
+                end
+                obj(s).value = exp([obj(s).index]);
+                for i = 2:obj(s).level
+                    obj.value = exp(obj.value);
+                end
+                if ~obj(s).reciprocal
+                    obj.value = 1/obj.value;
+                end
+                if obj(s).sign
+                    obj.value = -obj.value;
+                end
+            end
+        end
+
+        function r = update_sli(obj)
+            r = obj;
+            if (obj.level > 2^obj.level_bits)
+                r.level = 2^obj.level_bits;
+            end
+            r.index = round(2^obj.index_bits * ...
+                obj.index)/2^(obj.index_bits);
+            % Reciprocal flip if close to 1.
+            if (r.index == 0 && r.level == 1 && obj.index ~= 0)
+                r.reciprocal = 1;
+            end
+        end
+
+        function r = li_le(a, b)
+            if (a.sign > b.sign)
+                r = true;
+            elseif (a.sign < b.sign)
+                r = false;
+            elseif (a.level < b.level)
+                r = true;
+            elseif (a.level > b.level)
+                r = false;
+            else
+                r = a.index <= b.index;
+            end
+        end
+
+        % Following is a non-symmetric level-index addition operation
+        % which is used for symmetric level-index * and / operations.
+        %
+        % Implementation after "Level-Index Arithmetic Operations"
+        % by Clenshaw and Olver, SIAM J. Numer. Anal.24:2, 1987.
+        function r = li_add(a, b)
+            r = sli(a.level_bits, a.index_bits);
+            if (a.level == 0)
+                hzero = a.index + b.index;
+                if (hzero < 1)
+                    r.level = 0;
+                    r.index = hzero;
+                else
+                    r.level = 1;
+                    r.index = log(hzero);
+                end
+            else
+                a_seq(a.level) = exp(-a.index);
+                for (j = a.level-1:-1:1)
+                    a_seq(j) = exp(-1/a_seq(j+1));
+                end
+                if (b.level >= 1)
+                    b_seq(b.level) = a_seq(b.level)*exp(b.index);
+                else
+                    b_seq(1) = a_seq(1)*b.index;
+                end
+                for (j = b.level-1:-1:1)
+                    b_seq(j) = exp(-(1-b_seq(j+1))/a_seq(j+1));
+                end
+                c_seq(1) = 1 + b_seq(1);
+                for (j = 2:a.level)
+                    c_seq(j)=1+a_seq(j)*log(c_seq(j-1));
+                end
+                r.level = a.level;
+                r.index = a.index + log(c_seq(a.level));
+                if (r.index >= 1)
+                    r.level = r.level + 1;
+                    r.index = log(r.index);
+                end
+            end
+        end
+
+        % Following is a non-symmetric level-index subtract operation
+        % which is used for symmetric level-index * and / operations.
+        %
+        % Implementation after "Level-Index Arithmetic Operations"
+        % by Clenshaw and Olver, SIAM J. Numer. Anal.24:2, 1987.
+        function r = li_sub(a, b)
+            r = sli(a.level_bits, a.index_bits);
+            if (a.level == 0)
+                r.level = 0;
+                r.index = a.index - b.index;
+            else
+                a_seq(a.level) = exp(-a.index);
+                for (j = a.level-1:-1:1)
+                    a_seq(j) = exp(-1/a_seq(j+1));
+                end
+                if (b.level >= 1)
+                    b_seq(b.level) = a_seq(b.level)*exp(b.index);
+                else
+                    b_seq(1) = a_seq(1)*b.index;
+                end
+                for (j = b.level-1:-1:1)
+                    b_seq(j) = exp(-(1-b_seq(j+1))/a_seq(j+1));
+                end
+                c_seq(1) = 1 - b_seq(1);
+                if (c_seq(1) < a_seq(1))
+                    r.level = 0;
+                    r.index = c_seq(1)/a_seq(1);
+                    r = update_double(r);
+                    return;
+                end
+                for (j = 2:a.level)
+                    c_seq(j)=1+a_seq(j)*log(c_seq(j-1));
+                    if (c_seq(j) < a_seq(j))
+                        r.level = j-1;
+                        r.index = c_seq(j)/a_seq(j);
+                        r = update_double(r);
+                        return;
+                    end
+                end
+                r.level = a.level;
+                r.index = a.index + log(c_seq(a.level));
+            end
+        end
+
+        % Implementation after "The Symmetrix Level-Index System"
+        % by Clenshaw and Turner, IMA J. Numer. Analysis, 8:4, 1988.
+        function r = main_sli_algorithm(a, b, operation)
+            r = a;
+            % STEP 1
+            a_seq(a.level) = exp(-a.index);
+            for j=a.level-1:-1:1
+                a_seq(j) = exp(-1/a_seq(j+1));
+            end
+            b_seq = [];
+            % STEP 2
+            % Small case
+            if (a.reciprocal == 0) && (b.reciprocal == 0)
+                if (a.level < b.level)
+                    alpha_seq(b.level)=exp(-b.index);
+                    for j=b.level-1:-1:a.level
+                        alpha_seq(j) = exp(-1/alpha_seq(j+1));
+                    end
+                    beta_seq(a.level) = alpha_seq(a.level)/a_seq(a.level);
+                else
+                    beta_seq(a.level) = exp(a.index - b.index);
+                end
+                for j=a.level-1:-1:1
+                    beta_seq(j) = ...
+                        exp((beta_seq(j+1)-1)/(a_seq(j+1)*beta_seq(j+1)));
+                end
+                if strcmp(operation, 'add')
+                    c_seq(1) = 1 + beta_seq(1);
+                elseif strcmp(operation, 'sub')
+                    c_seq(1) = 1 - beta_seq(1);
+                end
+            % Mixed case
+            elseif (a.reciprocal == 1) && (b.reciprocal == 0)
+                alpha_seq(b.level) = exp(-b.index);
+                for j=b.level-1:-1:1
+                    alpha_seq(j) = ...
+                        exp(-1/alpha_seq(j+1));
+                end
+                if strcmp(operation, 'add')
+                    c_seq(1) = 1 + a_seq(1)*alpha_seq(1);
+                elseif strcmp(operation, 'sub')
+                    c_seq(1) = 1 - a_seq(1)*alpha_seq(1);
+                end
+            % Big case
+            else
+                b_seq(b.level) = a_seq(b.level)*exp(b.index);
+                for j=b.level-1:-1:1
+                    b_seq(j) = ...
+                        exp((b_seq(j+1)-1)/a_seq(j+1));
+                end
+                if (operation == 'add')
+                    c_seq(1) = 1 + b_seq(1);
+                elseif (operation == 'sub')
+                    c_seq(1) = 1 - b_seq(1);
+                end
+            end
+            % STEP 3
+            h_index = 1;
+            if (a.reciprocal == 1)
+                if (c_seq(1) < a_seq(1))
+                    h_index = 1;
+                    r.reciprocal = 0;
+                    h_seq(1) = -log(c_seq(1)/a_seq(1));
+                elseif (a.level == 1)
+                    h_index = 1;
+                    h_seq(1) = a.index + log(c_seq(1));
+                else
+                    c_seq(2) = 1 + a_seq(2)*log(c_seq(1));
+                    % STEP 4
+                    for j = 2:a.level-1
+                        if (c_seq(j) < a_seq(j))
+                            r.level = j-1;
+                            r.index = c_seq(j)/a_seq(j);
+                            r = update_double(c);
+                            return;
+                        else
+                            c_seq(j+1) = 1+a_seq(j+1)*log(c_seq(j));
+                        end
+                    end
+                    if (c_seq(a.level) < a_seq(a.level))
+                        r.level = a.level - 1;
+                        r.index = c_seq(a.level)/a_seq(a.level);
+                        r = update_double(r);
+                        return;
+                    else
+                        h_seq(a.level) = a.index + ...
+                            log(c_seq(a.level));
+                        h_index = a.level;
+                    end
+                end
+            else
+                if (a_seq(1) * c_seq(1) > 1)
+                    r.reciprocal = 1;
+                    r.level = 1;
+                    r.index = log(a_seq(1)*c_seq(1));
+                    r = update_double(r);
+                    return;
+                elseif (a.level == 1)
+                    h_seq(1) = a.index - log(c_seq(1));
+                    h_index = 1;
+                else
+                    c_seq(2) = 1 - a_seq(2)*log(c_seq(1));
+                    % STEP 4
+                    for j = 2:a.level-1
+                        if (c_seq(j) < a_seq(j))
+                            r.level = j-1;
+                            r.index = c_seq(j)/a_seq(j);
+                            r = update_double(r);
+                            return;
+                        else
+                            c_seq(j+1) = 1+a_seq(j+1)*log(c_seq(j));
+                        end
+                    end
+                    if (c_seq(a.level) < a_seq(a.level))
+                        r.level = a.level - 1;
+                        r.index = c_seq(a.level)/a_seq(a.level);
+                        r = update_double(r);
+                        return;
+                    else
+                        h_seq(a.level) = a.index + ...
+                            log(c_seq(a.level));
+                        h_index = a.level;
+                    end
+                end
+            end
+            % STEP 5
+            r.level = h_index;
+            h = h_seq(h_index);
+            while (h > 1)
+                h = log(h);
+                r.level = r.level + 1;
+            end
+            r.index = h;
+            r = update_double(r);
+        end
+    end
+end
\ No newline at end of file
diff --git a/tests/addition.m b/tests/addition.m
new file mode 100644
index 0000000..3d8332b
--- /dev/null
+++ b/tests/addition.m
@@ -0,0 +1,42 @@
+clear all;
+
+n = 10000;
+step = 10^-3;
+upper = n/2 * step;
+lower = -upper;
+
+% Set up cpfloat options
+options.format = 'fp16';
+options.round = 1;
+options.explim = 1;
+cpfloat([], options);
+k = 0;
+
+number = lower;
+
+for i = 1:n
+
+    numbers(i) = number;
+    addend = rand(1);
+    addendfp16 = cpfloat(addend, options);
+    addendsli = sli(2, 12);
+    addendsli = addendsli.set_val(addend);
+
+    yfp16 = cpfloat(number, options);
+    yfp16 = cpfloat(yfp16+addend, options);
+
+    ysli = sli(2, 12);
+    ysli = ysli.set_val(number);
+    ysli = ysli + addendsli;
+
+    if (number+addend ~= 0)
+        err1(i) = abs((number+addend - yfp16)/(number+addend));
+        err2(i) = abs((number+addend - ysli.value)/(number+addend));
+    end
+
+    number = number + step;
+end
+
+semilogy(numbers, err2, '-*', ...
+      numbers, err1, '-o');
+legend('level index 2.12', 'binary16');
\ No newline at end of file
diff --git a/tests/division.m b/tests/division.m
new file mode 100644
index 0000000..401faed
--- /dev/null
+++ b/tests/division.m
@@ -0,0 +1,43 @@
+clear all;
+
+n = 10000;
+step = 10^-3;
+upper = n/2 * step;
+lower = -upper;
+
+% Set up cpfloat options
+options.format = 'fp16';
+options.round = 1;
+options.explim = 1;
+cpfloat([], options);
+k = 0;
+
+number = lower;
+
+for i = 1:n
+
+    numbers(i) = number;
+    divisor = rand(1)*10;
+    divisors(i) = divisor;
+    divisorfp16 = cpfloat(divisor, options);
+    divisorsli = sli(2, 12);
+    divisorsli = divisorsli.set_val(divisor);
+
+    yfp16 = cpfloat(number, options);
+    yfp16 = cpfloat(yfp16/divisor, options);
+
+    ysli = sli(2, 12);
+    ysli = ysli.set_val(number);
+    ysli = ysli./divisorsli;
+
+    if (number/divisor ~= 0)
+        err1(i) = abs((number/divisor - yfp16)/(number/divisor));
+        err2(i) = abs((number/divisor - ysli.value)/(number/divisor));
+    end
+
+    number = number + step;
+end
+
+semilogy(numbers, err2, '-*', ...
+      numbers, err1, '-o');
+legend('level index 2.12', 'binary16');
\ No newline at end of file
diff --git a/tests/multiplication.m b/tests/multiplication.m
new file mode 100644
index 0000000..95d0bcc
--- /dev/null
+++ b/tests/multiplication.m
@@ -0,0 +1,45 @@
+clear all;
+
+n = 10000;
+step = 10^-3;
+upper = n/2 * step;
+lower = -upper;
+
+% Set up cpfloat options
+options.format = 'fp16';
+options.round = 1;
+options.explim = 1;
+cpfloat([], options);
+k = 0;
+
+number = lower;
+
+for i = 1:n
+
+    numbers(i) = number;
+    multiplicand = rand(1);
+    multiplicands(i) = multiplicand;
+    multiplicandfp16 = cpfloat(multiplicand, options);
+    multiplicandsli = sli(2, 12);
+    multiplicandsli = multiplicandsli.set_val(multiplicand);
+
+    yfp16 = cpfloat(number, options);
+    yfp16 = cpfloat(yfp16*multiplicand, options);
+
+    ysli = sli(2, 12);
+    ysli = ysli.set_val(number);
+    ysli = ysli * multiplicandsli;
+
+    if (number*multiplicand ~= 0)
+        err1(i) = abs((number*multiplicand - yfp16)/...
+            (number*multiplicand));
+        err2(i) = abs((number*multiplicand - ysli.value)/...
+            (number*multiplicand));
+    end
+
+    number = number + step;
+end
+
+semilogy(numbers, err2, '-*', ...
+      numbers, err1, '-o');
+legend('level index 2.12', 'binary16');
\ No newline at end of file
diff --git a/tests/representation.m b/tests/representation.m
new file mode 100644
index 0000000..7628043
--- /dev/null
+++ b/tests/representation.m
@@ -0,0 +1,43 @@
+clear all;
+
+n = 1000;
+step = 10^-5;
+upper = n/2 * step;
+lower = -upper;
+
+% Set up cpfloat options
+options.format = 'b';
+options.round = 1;
+options.explim = 1;
+cpfloat([], options);
+k = 0;
+
+number = lower;
+
+for i = 1:n
+
+    numbers(i) = number;
+
+    yfp16 = cpfloat(number, options);
+
+    ysli = sli(2, 12);
+    ysli = ysli.set_val(number);
+
+    if (number ~= 0)
+        err1(i) = abs((number - yfp16)/number);
+        err2(i) = abs((number - ysli.value)/number);
+    end
+
+    number = number + step;
+end
+
+filename = sprintf('representation_bfloat16.dat');
+fid = fopen(filename, 'w');
+for i=1:n
+    fprintf(fid, "%d %f %f\n", numbers(i), err1(i), err2(i));
+end
+fclose(fid);
+
+semilogy(numbers, err2, '-*', ...
+      numbers, err1, '-o');
+legend('level index 2.12', 'binary16');
\ No newline at end of file
diff --git a/tests/subtraction.m b/tests/subtraction.m
new file mode 100644
index 0000000..4c8ca1a
--- /dev/null
+++ b/tests/subtraction.m
@@ -0,0 +1,42 @@
+clear all;
+
+n = 10000;
+step = 10^-3;
+upper = n/2 * step;
+lower = -upper;
+
+% Set up cpfloat options
+options.format = 'fp16';
+options.round = 1;
+options.explim = 1;
+cpfloat([], options);
+k = 0;
+
+number = lower;
+
+for i = 1:n
+
+    numbers(i) = number;
+    addend = rand(1);
+    addendfp16 = cpfloat(addend, options);
+    addendsli = sli(2, 12);
+    addendsli = addendsli.set_val(addend);
+
+    yfp16 = cpfloat(number, options);
+    yfp16 = cpfloat(yfp16+addend, options);
+
+    ysli = sli(2, 12);
+    ysli = ysli.set_val(number);
+    ysli = ysli - addendsli;
+
+    if (number-addend ~= 0)
+        err1(i) = abs(((number-addend) - yfp16)/(number-addend));
+        err2(i) = abs(((number-addend) - ysli.value)/(number-addend));
+    end
+
+    number = number + step;
+end
+
+semilogy(numbers, err2, '-*', ...
+      numbers, err1, '-o');
+legend('level index 2.12', 'binary16');
\ No newline at end of file
diff --git a/tests/test_sli.m b/tests/test_sli.m
new file mode 100644
index 0000000..a1e089c
--- /dev/null
+++ b/tests/test_sli.m
@@ -0,0 +1,377 @@
+% Tests for various symmetric level-index arithmetic operations on sli.m
+% objects.
+
+% Addition
+
+x = sli(2, 12);
+y = sli(2, 12);
+
+x = x.set_val(1);
+y = y.set_val(0.2);
+z = x+y;
+assert(round(10*z.value)/10 == 1.2)
+
+x = x.set_val(1);
+y = y.set_val(0.2);
+z = y+x;
+assert(round(10*z.value)/10 == 1.2)
+
+x = x.set_val(1);
+y = y.set_val(-0.2);
+z = x+y;
+assert(round(10*z.value)/10 == 0.8)
+
+x = x.set_val(1);
+y = y.set_val(-0.2);
+z = y+x;
+assert(round(10*z.value)/10 == 0.8)
+
+x = x.set_val(-1);
+y = y.set_val(0.2);
+z = x+y;
+assert(round(10*z.value)/10 == -0.8)
+
+x = x.set_val(-1);
+y = y.set_val(0.2);
+z = y+x;
+assert(round(10*z.value)/10 == -0.8)
+
+x = x.set_val(-1);
+y = y.set_val(-0.2);
+z = x+y;
+assert(round(10*z.value)/10 == -1.2)
+
+x = x.set_val(-1);
+y = y.set_val(-0.2);
+z = y+x;
+assert(round(10*z.value)/10 == -1.2)
+
+% Subtraction
+
+x = x.set_val(1);
+y = y.set_val(0.2);
+z = x-y;
+assert(round(10*z.value)/10 == 0.8)
+
+x = x.set_val(1);
+y = y.set_val(0.2);
+z = y-x;
+assert(round(10*z.value)/10 == -0.8)
+
+x = x.set_val(1);
+y = y.set_val(-0.2);
+z = x-y;
+assert(round(10*z.value)/10 == 1.2)
+
+x = x.set_val(1);
+y = y.set_val(-0.2);
+z = y-x;
+assert(round(10*z.value)/10 == -1.2)
+
+x = x.set_val(-1);
+y = y.set_val(0.2);
+z = x-y;
+assert(round(10*z.value)/10 == -1.2)
+
+x = x.set_val(-1);
+y = y.set_val(0.2);
+z = y-x;
+assert(round(10*z.value)/10 == 1.2)
+
+x = x.set_val(-1);
+y = y.set_val(-0.2);
+z = x-y;
+assert(round(10*z.value)/10 == -0.8)
+
+x = x.set_val(-1);
+y = y.set_val(-0.2);
+z = y-x;
+assert(round(10*z.value)/10 == 0.8)
+
+% Multiplication
+
+x = x.set_val(2);
+y = y.set_val(2.5);
+z = x.*y;
+assert(round(z.value) == 5)
+z = y.*x;
+assert(round(z.value) == 5)
+
+y = y.set_val(-2.5);
+z = x.*y;
+assert(floor(z.value) == -5)
+z = y.*x;
+assert(floor(z.value) == -5)
+
+x = x.set_val(-2);
+z = x.*y;
+assert(ceil(z.value) == 5)
+z = y.*x;
+assert(ceil(z.value) == 5)
+
+x = x.set_val(2);
+y = y.set_val(10);
+z = x.*y;
+assert(round(z.value) == 20)
+z = y.*x;
+assert(round(z.value) == 20)
+
+y = y.set_val(1);
+z = x.*y;
+assert(ceil(z.value) == 2)
+z = y.*x;
+assert(ceil(z.value) == 2)
+
+y = y.set_val(0.5);
+z = x.*y;
+assert(round(z.value) == 1)
+z = y.*x;
+assert(round(z.value) == 1)
+
+x = x.set_val(0.8);
+y = y.set_val(0.5);
+y.value;
+z = x.*y;
+assert(round(z.value,1) == 0.4)
+z = y.*x;
+assert(round(z.value,1) == 0.4)
+
+% Division
+
+x = x.set_val(300);
+y = y.set_val(150);
+z = x./y;
+assert(round(z.value) == 2)
+z = y./x;
+assert(round(10*z.value)/10 == 0.5)
+
+y = y.set_val(-150);
+z = x./y;
+assert(round(z.value) == -2)
+z = y./x;
+assert(round(10*z.value)/10 == -0.5)
+
+x = x.set_val(-300);
+z = x./y;
+assert(round(z.value) == 2)
+z = y./x;
+assert(round(10*z.value)/10 == 0.5)
+
+x = x.set_val(300);
+y = y.set_val(10);
+z = x./y;
+assert(round(z.value) == 30)
+z = y./x;
+assert(floor(1000*z.value)/1000 == 0.033)
+
+y = y.set_val(1);
+z = x./y;
+assert(round(z.value) == 300)
+z = y./x;
+assert(floor(1000*z.value)/1000 == 0.003)
+
+x = x.set_val(2);
+y = y.set_val(0.5);
+z = x./y;
+assert(round(z.value) == 4)
+z = y./x;
+assert(round(z.value, 2) == 0.25)
+
+x = x.set_val(0.8);
+y = y.set_val(0.5);
+z = x./y;
+assert(round(10*z.value)/10 == 1.6)
+z = y./x;
+assert(round(1000*z.value)/1000 == 0.625)
+
+% Special cases
+
+x = x.set_val(0);
+y = y.set_val(2);
+z = x.*y;
+assert(z.value == 0);
+z = y.*z;
+assert(z.value == 0);
+z = x+y;
+assert(round(z.value) == 2);
+z = y+x;
+assert(round(z.value) == 2);
+z = x-y;
+assert(round(z.value) == -2);
+z = y-x;
+assert(round(z.value) == 2);
+z = x./y;
+assert(z.value == 0);
+zero_div_err = 0;
+try z = y./x;
+catch
+    zero_div_err = 1;
+end
+assert(zero_div_err == 1);
+
+% Vectors and matrices
+
+% A = [1 2 3] A' = [1 4]
+%     [4 5 6]      [2 5]
+%                  [3 6]
+%
+A = zeros(2, 3, 2, 12, 'sli');
+A(1,1) = A(1,1).set_val(1);
+A(1,2) = A(1,2).set_val(2);
+A(1,3) = A(1,3).set_val(3);
+A(2,1) = A(2,1).set_val(4);
+A(2,2) = A(2,2).set_val(5);
+A(2,3) = A(2,3).set_val(6);
+
+Y = A*A';
+assert(round(Y(1,1).value) == 14)
+assert(round(Y(1,2).value) == 32)
+assert(round(Y(2,1).value) == 32)
+assert(round(Y(2,2).value) == 77)
+
+% b = [1 2 3]
+b = zeros(3, 2, 12, 'sli');
+b(:) = A(1, :);
+Y = A*b';
+assert(round(Y(1,1).value) == 14)
+assert(round(Y(2,1).value) == 32)
+
+Y = A+A;
+assert(round(Y(1,1).value) == 2)
+assert(round(Y(1,2).value) == 4)
+assert(round(Y(1,3).value) == 6)
+assert(round(Y(2,1).value) == 8)
+assert(round(Y(2,2).value) == 10)
+assert(round(Y(2,3).value) == 12)
+
+Y = A-A;
+assert(round(Y(1,1).value) == 0)
+assert(round(Y(1,2).value) == 0)
+assert(round(Y(1,3).value) == 0)
+assert(round(Y(2,1).value) == 0)
+assert(round(Y(2,2).value) == 0)
+assert(round(Y(2,3).value) == 0)
+
+Y = A.*A;
+assert(round(Y(1,1).value) == 1)
+assert(round(Y(1,2).value) == 4)
+assert(round(Y(1,3).value) == 9)
+assert(round(Y(2,1).value) == 16)
+assert(round(Y(2,2).value) == 25)
+assert(round(Y(2,3).value) == 36)
+
+Y = A./A;
+assert(round(Y(1,1).value) == 1)
+assert(round(Y(1,2).value) == 1)
+assert(round(Y(1,3).value) == 1)
+assert(round(Y(2,1).value) == 1)
+assert(round(Y(2,2).value) == 1)
+assert(round(Y(2,3).value) == 1)
+
+Y = A.\A;
+assert(round(Y(1,1).value) == 1)
+assert(round(Y(1,2).value) == 1)
+assert(round(Y(1,3).value) == 1)
+assert(round(Y(2,1).value) == 1)
+assert(round(Y(2,2).value) == 1)
+assert(round(Y(2,3).value) == 1)
+
+x = x.set_val(1);
+Y = A+x;
+assert(round(Y(1,1).value) == 2)
+assert(round(Y(1,2).value) == 3)
+assert(round(Y(1,3).value) == 4)
+assert(round(Y(2,1).value) == 5)
+assert(round(Y(2,2).value) == 6)
+assert(round(Y(2,3).value) == 7)
+Y = x+A;
+assert(round(Y(1,1).value) == 2)
+assert(round(Y(1,2).value) == 3)
+assert(round(Y(1,3).value) == 4)
+assert(round(Y(2,1).value) == 5)
+assert(round(Y(2,2).value) == 6)
+assert(round(Y(2,3).value) == 7)
+
+ Y = A-x;
+assert(round(Y(1,1).value) == 0)
+assert(round(Y(1,2).value) == 1)
+assert(round(Y(1,3).value) == 2)
+assert(round(Y(2,1).value) == 3)
+assert(round(Y(2,2).value) == 4)
+assert(round(Y(2,3).value) == 5)
+Y = x-A;
+assert(round(Y(1,1).value) == 0)
+assert(round(Y(1,2).value) == -1)
+assert(round(Y(1,3).value) == -2)
+assert(round(Y(2,1).value) == -3)
+assert(round(Y(2,2).value) == -4)
+assert(round(Y(2,3).value) == -5)
+
+x = x.set_val(2);
+Y = A.*x;
+assert(round(Y(1,1).value) == 2)
+assert(round(Y(1,2).value) == 4)
+assert(round(Y(1,3).value) == 6)
+assert(round(Y(2,1).value) == 8)
+assert(round(Y(2,2).value) == 10)
+assert(round(Y(2,3).value) == 12)
+Y = x.*A;
+assert(round(Y(1,1).value) == 2)
+assert(round(Y(1,2).value) == 4)
+assert(round(Y(1,3).value) == 6)
+assert(round(Y(2,1).value) == 8)
+assert(round(Y(2,2).value) == 10)
+assert(round(Y(2,3).value) == 12)
+
+Y = A./x;
+assert(round(10*Y(1,1).value)/10 == 0.5)
+assert(round(Y(1,2).value) == 1)
+assert(round(10*Y(1,3).value)/10 == 1.5)
+assert(round(Y(2,1).value) == 2)
+assert(round(10*Y(2,2).value)/10 == 2.5)
+assert(round(Y(2,3).value) == 3)
+Y = x./A;
+assert(round(Y(1,1).value) == 2)
+assert(round(Y(1,2).value) == 1)
+assert(round(10*Y(1,3).value)/10 == 0.7)
+assert(round(10*Y(2,1).value)/10 == 0.5)
+assert(round(10*Y(2,2).value)/10 == 0.4)
+assert(round(10*Y(2,3).value)/10 == 0.3)
+
+% A = [1 2 3] x = [1]
+%     [4 5 6]     [2]
+%            
+%
+x = zeros(2, 1, 2, 12, 'sli');
+x(1,1) = x(1,1).set_val(1);
+x(2,1) = x(2,1).set_val(2);
+Y = A+x;
+assert(round(Y(1,1).value) == 2)
+assert(round(Y(1,2).value) == 3)
+assert(round(Y(1,3).value) == 4)
+assert(round(Y(2,1).value) == 6)
+assert(round(Y(2,2).value) == 7)
+assert(round(Y(2,3).value) == 8)
+
+% A = [1 2 3] x = [1]
+%                 [2]
+%            
+A = A(1, :);
+Y = A+x;
+assert(round(Y(1,1).value) == 2)
+assert(round(Y(1,2).value) == 3)
+assert(round(Y(1,3).value) == 4)
+assert(round(Y(2,1).value) == 3)
+assert(round(Y(2,2).value) == 4)
+assert(round(Y(2,3).value) == 5)
+
+% A = [1] x = [1, 2]
+%     [2]
+%     [3]       
+Y = A'+x';
+assert(round(Y(1,1).value) == 2)
+assert(round(Y(1,2).value) == 3)
+assert(round(Y(2,1).value) == 3)
+assert(round(Y(2,2).value) == 4)
+assert(round(Y(3,1).value) == 4)
+assert(round(Y(3,2).value) == 5)
\ No newline at end of file