diff --git a/hnswlib/space_ip.h b/hnswlib/space_ip.h
index c0029bde..ac50f93c 100644
--- a/hnswlib/space_ip.h
+++ b/hnswlib/space_ip.h
@@ -2,16 +2,44 @@
 #include "hnswlib.h"
 
 namespace hnswlib {
-
-    static float
+    /**
+     * For a given loop unrolling factor K, distance type dist_t, and data type data_t,
+     * calculate the inner product distance between two vectors.
+     * The compiler should automatically do the loop unrolling for us here and vectorize as appropriate.
+     */
+    template<typename dist_t, typename data_t = dist_t, int K = 1>
+    static dist_t
     InnerProduct(const void *pVect1, const void *pVect2, const void *qty_ptr) {
         size_t qty = *((size_t *) qty_ptr);
-        float res = 0;
-        for (unsigned i = 0; i < qty; i++) {
-            res += ((float *) pVect1)[i] * ((float *) pVect2)[i];
+        dist_t res = 0;
+        data_t *a = (data_t *) pVect1;
+        data_t *b = (data_t *) pVect2;
+
+        qty = qty / K;
+
+        for (size_t i = 0; i < qty; i++) {
+            for (size_t j = 0; j < K; j++) {
+                const size_t index = (i * K) + j;
+                const dist_t _a = a[index];
+                const dist_t _b = b[index];
+                res += _a * _b;
+            }
         }
-        return (1.0f - res);
 
+        return (static_cast<dist_t>(1.0f) - res);
+    }
+
+    template<typename dist_t, typename data_t = dist_t, int K>
+    static dist_t
+    InnerProductAtLeast(const void *__restrict pVect1, const void *__restrict pVect2, const void *__restrict qty_ptr) {
+        size_t k = K;
+        size_t remainder = *((size_t *) qty_ptr) - K;
+
+        data_t *a = (data_t *) pVect1;
+        data_t *b = (data_t *) pVect2;
+
+        return InnerProduct<dist_t, data_t, K>(a, b, &k)
+             + InnerProduct<dist_t, data_t, 1>(a + K, b + K, &remainder);
     }
 
 #if defined(USE_AVX)
@@ -254,7 +282,7 @@ namespace hnswlib {
         float *pVect2 = (float *) pVect2v + qty16;
 
         size_t qty_left = qty - qty16;
-        float res_tail = InnerProduct(pVect1, pVect2, &qty_left);
+        float res_tail = InnerProduct<float, float>(pVect1, pVect2, &qty_left);
         return res + res_tail - 1.0f;
     }
 
@@ -268,48 +296,75 @@ namespace hnswlib {
 
         float *pVect1 = (float *) pVect1v + qty4;
         float *pVect2 = (float *) pVect2v + qty4;
-        float res_tail = InnerProduct(pVect1, pVect2, &qty_left);
+        float res_tail = InnerProduct<float, float>(pVect1, pVect2, &qty_left);
 
         return res + res_tail - 1.0f;
     }
 #endif
 
-    class InnerProductSpace : public SpaceInterface<float> {
+    template<typename dist_t, typename data_t = dist_t>
+    class InnerProductSpace : public SpaceInterface<dist_t> {
 
-        DISTFUNC<float> fstdistfunc_;
+        DISTFUNC<dist_t> fstdistfunc_;
         size_t data_size_;
         size_t dim_;
     public:
-        InnerProductSpace(size_t dim) {
-            fstdistfunc_ = InnerProduct;
-    #if defined(USE_AVX) || defined(USE_SSE) || defined(USE_AVX512)
-            if (dim % 16 == 0)
-                fstdistfunc_ = InnerProductSIMD16Ext;
+        InnerProductSpace(size_t dim) : dim_(dim), data_size_(dim * sizeof(data_t)) {
+            if (dim % 128 == 0)
+                fstdistfunc_ = InnerProduct<dist_t, data_t, 128>;
+            else if (dim % 64 == 0)
+                fstdistfunc_ = InnerProduct<dist_t, data_t, 64>;
+            else if (dim % 32 == 0)
+                fstdistfunc_ = InnerProduct<dist_t, data_t, 32>;
+            else if (dim % 16 == 0)
+                fstdistfunc_ = InnerProduct<dist_t, data_t, 16>;
+            else if (dim % 8 == 0)
+                fstdistfunc_ = InnerProduct<dist_t, data_t, 8>;
             else if (dim % 4 == 0)
-                fstdistfunc_ = InnerProductSIMD4Ext;
+                fstdistfunc_ = InnerProduct<dist_t, data_t, 4>;
+
+            else if (dim > 128)
+                fstdistfunc_ = InnerProductAtLeast<dist_t, data_t, 128>;            
+            else if (dim > 64)
+                fstdistfunc_ = InnerProductAtLeast<dist_t, data_t, 64>;
+            else if (dim > 32)
+                fstdistfunc_ = InnerProductAtLeast<dist_t, data_t, 32>;
             else if (dim > 16)
-                fstdistfunc_ = InnerProductSIMD16ExtResiduals;
+                fstdistfunc_ = InnerProductAtLeast<dist_t, data_t, 16>;
+            else if (dim > 8)
+                fstdistfunc_ = InnerProductAtLeast<dist_t, data_t, 8>;
             else if (dim > 4)
-                fstdistfunc_ = InnerProductSIMD4ExtResiduals;
-    #endif
-            dim_ = dim;
-            data_size_ = dim * sizeof(float);
+                fstdistfunc_ = InnerProductAtLeast<dist_t, data_t, 4>;
+            else
+                fstdistfunc_ = InnerProduct<dist_t, data_t>;
         }
 
         size_t get_data_size() {
             return data_size_;
         }
 
-        DISTFUNC<float> get_dist_func() {
+        DISTFUNC<dist_t> get_dist_func() {
             return fstdistfunc_;
         }
 
         void *get_dist_func_param() {
             return &dim_;
         }
-
-    ~InnerProductSpace() {}
+        ~InnerProductSpace() {}
     };
 
+    template<> InnerProductSpace<float, float>::InnerProductSpace(size_t dim) : dim_(dim), data_size_(dim * sizeof(float)) {
+        fstdistfunc_ = InnerProduct<float, float>;
+    #if defined(USE_SSE) || defined(USE_AVX) || defined(USE_AVX512)
+        if (dim % 16 == 0)
+            fstdistfunc_ = InnerProductSIMD16Ext;
+        else if (dim % 4 == 0)
+            fstdistfunc_ = InnerProductSIMD4Ext;
+        else if (dim > 16)
+            fstdistfunc_ = InnerProductSIMD16ExtResiduals;
+        else if (dim > 4)
+            fstdistfunc_ = InnerProductSIMD4ExtResiduals;
+    #endif
+    }
 
 }
diff --git a/hnswlib/space_l2.h b/hnswlib/space_l2.h
index 3b6a49ef..7d9415c1 100644
--- a/hnswlib/space_l2.h
+++ b/hnswlib/space_l2.h
@@ -2,23 +2,45 @@
 #include "hnswlib.h"
 
 namespace hnswlib {
-
-    static float
-    L2Sqr(const void *pVect1v, const void *pVect2v, const void *qty_ptr) {
-        float *pVect1 = (float *) pVect1v;
-        float *pVect2 = (float *) pVect2v;
+    /**
+     * For a given loop unrolling factor K, distance type dist_t, and data type data_t,
+     * calculate the L2 squared distance between two vectors.
+     * The compiler should automatically do the loop unrolling for us here and vectorize as appropriate.
+     */
+    template<typename dist_t, typename data_t = dist_t, int K = 1>
+    static dist_t
+    L2Sqr(const void *__restrict pVect1, const void *__restrict pVect2, const void *__restrict qty_ptr) {
         size_t qty = *((size_t *) qty_ptr);
+        dist_t res = 0;
+        data_t *a = (data_t *) pVect1;
+        data_t *b = (data_t *) pVect2;
+
+        qty = qty / K;
 
-        float res = 0;
         for (size_t i = 0; i < qty; i++) {
-            float t = *pVect1 - *pVect2;
-            pVect1++;
-            pVect2++;
-            res += t * t;
+            for (size_t j = 0; j < K; j++) {
+                const size_t index = (i * K) + j;
+                const dist_t _a = a[index];
+                const dist_t _b = b[index];
+                res += (_a - _b) * (_a - _b);
+            }
         }
         return (res);
     }
 
+    template<typename dist_t, typename data_t = dist_t, int K>
+    static dist_t
+    L2SqrAtLeast(const void *__restrict pVect1, const void *__restrict pVect2, const void *__restrict qty_ptr) {
+        size_t k = K;
+        size_t remainder = *((size_t *) qty_ptr) - K;
+
+        data_t *a = (data_t *) pVect1;
+        data_t *b = (data_t *) pVect2;
+
+        return L2Sqr<dist_t, data_t, K>(a, b, &k)
+             + L2Sqr<dist_t, data_t, 1>(a + K, b + K, &remainder);
+    }
+
 #if defined(USE_AVX512)
 
     // Favor using AVX512 if available.
@@ -150,7 +172,7 @@ namespace hnswlib {
         float *pVect2 = (float *) pVect2v + qty16;
 
         size_t qty_left = qty - qty16;
-        float res_tail = L2Sqr(pVect1, pVect2, &qty_left);
+        float res_tail = L2Sqr<float, float>(pVect1, pVect2, &qty_left);
         return (res + res_tail);
     }
 #endif
@@ -194,39 +216,54 @@ namespace hnswlib {
 
         float *pVect1 = (float *) pVect1v + qty4;
         float *pVect2 = (float *) pVect2v + qty4;
-        float res_tail = L2Sqr(pVect1, pVect2, &qty_left);
+        float res_tail = L2Sqr<float, float>(pVect1, pVect2, &qty_left);
 
         return (res + res_tail);
     }
 #endif
 
-    class L2Space : public SpaceInterface<float> {
+    template <typename dist_t, typename data_t = dist_t>
+    class L2Space : public SpaceInterface<dist_t> {
 
-        DISTFUNC<float> fstdistfunc_;
+        DISTFUNC<dist_t> fstdistfunc_;
         size_t data_size_;
         size_t dim_;
     public:
-        L2Space(size_t dim) {
-            fstdistfunc_ = L2Sqr;
-        #if defined(USE_SSE) || defined(USE_AVX) || defined(USE_AVX512)
-            if (dim % 16 == 0)
-                fstdistfunc_ = L2SqrSIMD16Ext;
+        L2Space(size_t dim) : dim_(dim), data_size_(dim * sizeof(data_t)) {
+            if (dim % 128 == 0)
+                fstdistfunc_ = L2Sqr<dist_t, data_t, 128>;
+            else if (dim % 64 == 0)
+                fstdistfunc_ = L2Sqr<dist_t, data_t, 64>;
+            else if (dim % 32 == 0)
+                fstdistfunc_ = L2Sqr<dist_t, data_t, 32>;
+            else if (dim % 16 == 0)
+                fstdistfunc_ = L2Sqr<dist_t, data_t, 16>;
+            else if (dim % 8 == 0)
+                fstdistfunc_ = L2Sqr<dist_t, data_t, 8>;
             else if (dim % 4 == 0)
-                fstdistfunc_ = L2SqrSIMD4Ext;
+                fstdistfunc_ = L2Sqr<dist_t, data_t, 4>;
+
+            else if (dim > 128)
+                fstdistfunc_ = L2SqrAtLeast<dist_t, data_t, 128>;            
+            else if (dim > 64)
+                fstdistfunc_ = L2SqrAtLeast<dist_t, data_t, 64>;
+            else if (dim > 32)
+                fstdistfunc_ = L2SqrAtLeast<dist_t, data_t, 32>;
             else if (dim > 16)
-                fstdistfunc_ = L2SqrSIMD16ExtResiduals;
+                fstdistfunc_ = L2SqrAtLeast<dist_t, data_t, 16>;
+            else if (dim > 8)
+                fstdistfunc_ = L2SqrAtLeast<dist_t, data_t, 8>;
             else if (dim > 4)
-                fstdistfunc_ = L2SqrSIMD4ExtResiduals;
-        #endif
-            dim_ = dim;
-            data_size_ = dim * sizeof(float);
+                fstdistfunc_ = L2SqrAtLeast<dist_t, data_t, 4>;
+            else
+                fstdistfunc_ = L2Sqr<dist_t, data_t>;
         }
 
         size_t get_data_size() {
             return data_size_;
         }
 
-        DISTFUNC<float> get_dist_func() {
+        DISTFUNC<dist_t> get_dist_func() {
             return fstdistfunc_;
         }
 
@@ -237,79 +274,17 @@ namespace hnswlib {
         ~L2Space() {}
     };
 
-    static int
-    L2SqrI4x(const void *__restrict pVect1, const void *__restrict pVect2, const void *__restrict qty_ptr) {
-
-        size_t qty = *((size_t *) qty_ptr);
-        int res = 0;
-        unsigned char *a = (unsigned char *) pVect1;
-        unsigned char *b = (unsigned char *) pVect2;
-
-        qty = qty >> 2;
-        for (size_t i = 0; i < qty; i++) {
-
-            res += ((*a) - (*b)) * ((*a) - (*b));
-            a++;
-            b++;
-            res += ((*a) - (*b)) * ((*a) - (*b));
-            a++;
-            b++;
-            res += ((*a) - (*b)) * ((*a) - (*b));
-            a++;
-            b++;
-            res += ((*a) - (*b)) * ((*a) - (*b));
-            a++;
-            b++;
-        }
-        return (res);
+    template<> L2Space<float, float>::L2Space(size_t dim) : dim_(dim), data_size_(dim * sizeof(float)) {
+        fstdistfunc_ = L2Sqr<float, float>;
+    #if defined(USE_SSE) || defined(USE_AVX) || defined(USE_AVX512)
+        if (dim % 16 == 0)
+            fstdistfunc_ = L2SqrSIMD16Ext;
+        else if (dim % 4 == 0)
+            fstdistfunc_ = L2SqrSIMD4Ext;
+        else if (dim > 16)
+            fstdistfunc_ = L2SqrSIMD16ExtResiduals;
+        else if (dim > 4)
+            fstdistfunc_ = L2SqrSIMD4ExtResiduals;
+    #endif
     }
-
-    static int L2SqrI(const void* __restrict pVect1, const void* __restrict pVect2, const void* __restrict qty_ptr) {
-        size_t qty = *((size_t*)qty_ptr);
-        int res = 0;
-        unsigned char* a = (unsigned char*)pVect1;
-        unsigned char* b = (unsigned char*)pVect2;
-
-        for(size_t i = 0; i < qty; i++)
-        {
-            res += ((*a) - (*b)) * ((*a) - (*b));
-            a++;
-            b++;
-        }
-        return (res);
-    }
-
-    class L2SpaceI : public SpaceInterface<int> {
-
-        DISTFUNC<int> fstdistfunc_;
-        size_t data_size_;
-        size_t dim_;
-    public:
-        L2SpaceI(size_t dim) {
-            if(dim % 4 == 0) {
-                fstdistfunc_ = L2SqrI4x;
-            }
-            else {
-                fstdistfunc_ = L2SqrI;
-            }
-            dim_ = dim;
-            data_size_ = dim * sizeof(unsigned char);
-        }
-
-        size_t get_data_size() {
-            return data_size_;
-        }
-
-        DISTFUNC<int> get_dist_func() {
-            return fstdistfunc_;
-        }
-
-        void *get_dist_func_param() {
-            return &dim_;
-        }
-
-        ~L2SpaceI() {}
-    };
-
-
 }
diff --git a/python_bindings/bindings.cpp b/python_bindings/bindings.cpp
index 90d26161..3ad911c7 100644
--- a/python_bindings/bindings.cpp
+++ b/python_bindings/bindings.cpp
@@ -80,20 +80,20 @@ inline void assert_true(bool expr, const std::string & msg) {
 }
 
 
-template<typename dist_t, typename data_t=float>
+template<typename dist_t, typename data_t=dist_t>
 class Index {
 public:
   Index(const std::string &space_name, const int dim) :
   space_name(space_name), dim(dim) {
     normalize=false;
     if(space_name=="l2") {
-      l2space = new hnswlib::L2Space(dim);
+      l2space = new hnswlib::L2Space<dist_t, data_t>(dim);
     }
     else if(space_name=="ip") {
-      l2space = new hnswlib::InnerProductSpace(dim);
+      l2space = new hnswlib::InnerProductSpace<dist_t, data_t>(dim);
     }
     else if(space_name=="cosine") {
-      l2space = new hnswlib::InnerProductSpace(dim);
+      l2space = new hnswlib::InnerProductSpace<dist_t, data_t>(dim);
       normalize=true;
     } else {
       throw new std::runtime_error("Space name must be one of l2, ip, or cosine.");
@@ -119,7 +119,7 @@ class Index {
   int num_threads_default;
   hnswlib::labeltype cur_l;
   hnswlib::HierarchicalNSW<dist_t> *appr_alg;
-  hnswlib::SpaceInterface<float> *l2space;
+  hnswlib::SpaceInterface<dist_t> *l2space;
 
   ~Index() {
       delete l2space;
@@ -163,8 +163,8 @@ class Index {
       index_inited = true;
     }
 
-    void normalize_vector(float *data, float *norm_array){
-        float norm=0.0f;
+    void normalize_vector(data_t *data, data_t *norm_array){
+        data_t norm=0.0f;
         for(int i=0;i<dim;i++)
             norm+=data[i]*data[i];
         norm= 1.0f / (sqrtf(norm) + 1e-30f);
@@ -173,7 +173,7 @@ class Index {
     }
 
     void addItems(py::object input, py::object ids_ = py::none(), int num_threads = -1) {
-        py::array_t < dist_t, py::array::c_style | py::array::forcecast > items(input);
+        py::array_t < data_t, py::array::c_style | py::array::forcecast > items(input);
         auto buffer = items.request();
         if (num_threads <= 0)
             num_threads = num_threads_default;
@@ -223,8 +223,8 @@ class Index {
           int start = 0;
           if (!ep_added) {
             size_t id = ids.size() ? ids.at(0) : (cur_l);
-            float *vector_data = (float *) items.data(0);
-            std::vector<float> norm_array(dim);
+            data_t *vector_data = (data_t *) items.data(0);
+            std::vector<data_t> norm_array(dim);
             if(normalize){
               normalize_vector(vector_data, norm_array.data());
               vector_data = norm_array.data();
@@ -241,11 +241,11 @@ class Index {
                     appr_alg->addPoint((void *) items.data(row), (size_t) id);
                 });
             } else{
-                std::vector<float> norm_array(num_threads * dim);
+                std::vector<data_t> norm_array(num_threads * dim);
                 ParallelFor(start, rows, num_threads, [&](size_t row, size_t threadId) {
                     // normalize vector:
                     size_t start_idx = threadId * dim;
-                    normalize_vector((float *) items.data(row), (norm_array.data()+start_idx));
+                    normalize_vector((data_t *) items.data(row), (norm_array.data()+start_idx));
 
                     size_t id = ids.size() ? ids.at(row) : (cur_l+row);
                     appr_alg->addPoint((void *) (norm_array.data()+start_idx), (size_t) id);
@@ -255,7 +255,7 @@ class Index {
         }
     }
 
-    std::vector<std::vector<data_t>> getDataReturnList(py::object ids_ = py::none()) {
+    py::array_t<data_t> getDataReturnNumpy(py::object ids_ = py::none()) {
         std::vector<size_t> ids;
         if (!ids_.is_none()) {
             py::array_t < size_t, py::array::c_style | py::array::forcecast > items(ids_);
@@ -267,11 +267,21 @@ class Index {
             ids.swap(ids1);
         }
 
-        std::vector<std::vector<data_t>> data;
-        for (auto id : ids) {
-            data.push_back(appr_alg->template getDataByLabel<data_t>(id));
+        data_t *data_numpy = (data_t *)malloc(sizeof(data_t) * dim * ids.size());
+        for (int i = 0; i < ids.size(); i++) {
+            std::vector<data_t> vector = appr_alg->template getDataByLabel<data_t>(ids[i]);
+            std::copy(vector.begin(), vector.end(), data_numpy + (i * dim));
         }
-        return data;
+
+        py::capsule free_when_done(data_numpy, [](void *f) {
+            delete[] f;
+        });
+
+        return py::array_t<data_t>(
+            {(size_t) ids.size(), (size_t) dim}, // shape
+            {dim * sizeof(data_t), sizeof(data_t)}, // strides
+            data_numpy, // the data pointer
+            free_when_done);
     }
 
     std::vector<hnswlib::labeltype> getIdsList() {
@@ -400,7 +410,7 @@ class Index {
 
     py::dict getIndexParams() const { /* WARNING: Index::getAnnData is not thread-safe with Index::addItems */
         auto params = py::dict(
-                            "ser_version"_a=py::int_(Index<float>::ser_version), //serialization version
+                            "ser_version"_a=py::int_(Index<dist_t, data_t>::ser_version), //serialization version
                             "space"_a=space_name,
                             "dim"_a=dim,
                             "index_inited"_a=index_inited,
@@ -418,16 +428,15 @@ class Index {
         return py::dict(**params, **ann_params);
     }
 
-
-    static Index<float> * createFromParams(const py::dict d) {
+    static Index<dist_t, data_t> * createFromParams(const py::dict d) {
       // check serialization version
-      assert_true(((int)py::int_(Index<float>::ser_version)) >= d["ser_version"].cast<int>(), "Invalid serialization version!");
+      assert_true(((int)py::int_(Index<dist_t, data_t>::ser_version)) >= d["ser_version"].cast<int>(), "Invalid serialization version!");
 
       auto space_name_=d["space"].cast<std::string>();
       auto dim_=d["dim"].cast<int>();
       auto index_inited_=d["index_inited"].cast<bool>();
 
-      Index<float> *new_index = new Index<float>(space_name_, dim_);
+      Index<dist_t, data_t> *new_index = new Index<dist_t, data_t>(space_name_, dim_);
 
       /*  TODO: deserialize state of random generators into new_index->level_generator_ and new_index->update_probability_generator_  */
       /*        for full reproducibility / state of generators is serialized inside Index::getIndexParams                      */
@@ -449,7 +458,7 @@ class Index {
       return new_index;
     }
 
-    static Index<float> * createFromIndex(const Index<float> & index) {
+    static Index<dist_t, data_t> * createFromIndex(const Index<dist_t, data_t> & index) {
         return createFromParams(index.getIndexParams());
     }
 
@@ -534,7 +543,7 @@ class Index {
 }
 
     py::object knnQuery_return_numpy(py::object input, size_t k = 1, int num_threads = -1) {
-        py::array_t < dist_t, py::array::c_style | py::array::forcecast > items(input);
+        py::array_t < data_t, py::array::c_style | py::array::forcecast > items(input);
         auto buffer = items.request();
         hnswlib::labeltype *data_numpy_l;
         dist_t *data_numpy_d;
@@ -582,12 +591,12 @@ class Index {
                 );
             }
             else{
-                std::vector<float> norm_array(num_threads*features);
+                std::vector<data_t> norm_array(num_threads*features);
                 ParallelFor(0, rows, num_threads, [&](size_t row, size_t threadId) {
-                                float *data= (float *) items.data(row);
+                                data_t *data= (data_t *) items.data(row);
 
                                 size_t start_idx = threadId * dim;
-                                normalize_vector((float *) items.data(row), (norm_array.data()+start_idx));
+                                normalize_vector((data_t *) items.data(row), (norm_array.data()+start_idx));
 
                                 std::priority_queue<std::pair<dist_t, hnswlib::labeltype >> result = appr_alg->searchKnn(
                                         (void *) (norm_array.data()+start_idx), k);
@@ -615,12 +624,12 @@ class Index {
                 py::array_t<hnswlib::labeltype>(
                         {rows, k}, // shape
                         {k * sizeof(hnswlib::labeltype),
-                         sizeof(hnswlib::labeltype)}, // C-style contiguous strides for double
+                         sizeof(hnswlib::labeltype)}, // C-style contiguous strides for labeltype
                         data_numpy_l, // the data pointer
                         free_when_done_l),
                 py::array_t<dist_t>(
                         {rows, k}, // shape
-                        {k * sizeof(dist_t), sizeof(dist_t)}, // C-style contiguous strides for double
+                        {k * sizeof(dist_t), sizeof(dist_t)}, // C-style contiguous strides for dist_t
                         data_numpy_d, // the data pointer
                         free_when_done_d));
 
@@ -650,17 +659,17 @@ class Index {
 template<typename dist_t, typename data_t=float>
 class BFIndex {
 public:
-    BFIndex(const std::string &space_name, const int dim) :
+    BFIndex(const std::string &space_name, const size_t dim) :
             space_name(space_name), dim(dim) {
         normalize=false;
         if(space_name=="l2") {
-            space = new hnswlib::L2Space(dim);
+            space = new hnswlib::L2Space<dist_t, data_t>(dim);
         }
         else if(space_name=="ip") {
-            space = new hnswlib::InnerProductSpace(dim);
+            space = new hnswlib::InnerProductSpace<dist_t, data_t>(dim);
         }
         else if(space_name=="cosine") {
-            space = new hnswlib::InnerProductSpace(dim);
+            space = new hnswlib::InnerProductSpace<dist_t, data_t>(dim);
             normalize=true;
         } else {
             throw new std::runtime_error("Space name must be one of l2, ip, or cosine.");
@@ -678,7 +687,7 @@ class BFIndex {
 
     hnswlib::labeltype cur_l;
     hnswlib::BruteforceSearch<dist_t> *alg;
-    hnswlib::SpaceInterface<float> *space;
+    hnswlib::SpaceInterface<dist_t> *space;
 
     ~BFIndex() {
         delete space;
@@ -695,8 +704,8 @@ class BFIndex {
         index_inited = true;
     }
 
-    void normalize_vector(float *data, float *norm_array){
-        float norm=0.0f;
+    void normalize_vector(data_t *data, data_t *norm_array){
+        data_t norm=0.0f;
         for(int i=0;i<dim;i++)
             norm+=data[i]*data[i];
         norm= 1.0f / (sqrtf(norm) + 1e-30f);
@@ -744,8 +753,8 @@ class BFIndex {
                 if (!normalize) {
                     alg->addPoint((void *) items.data(row), (size_t) id);
                 } else {
-                    std::vector<float> normalized_vector(dim);
-                    normalize_vector((float *)items.data(row), normalized_vector.data());
+                    std::vector<data_t> normalized_vector(dim);
+                    normalize_vector((data_t *)items.data(row), normalized_vector.data());
                     alg->addPoint((void *) normalized_vector.data(), (size_t) id);
                 }
             }
@@ -830,68 +839,103 @@ class BFIndex {
 
 };
 
-PYBIND11_PLUGIN(hnswlib) {
-        py::module m("hnswlib");
 
-        py::class_<Index<float>>(m, "Index")
-        .def(py::init(&Index<float>::createFromParams), py::arg("params"))
+template<typename T>
+inline const std::string typeName();
+
+template<> const std::string typeName<char>() { return "int8"; }
+template<> const std::string typeName<unsigned char>() { return "uint8"; }
+template<> const std::string typeName<short>() { return "int16"; }
+template<> const std::string typeName<unsigned short>() { return "uint16"; }
+template<> const std::string typeName<int>() { return "int32"; }
+template<> const std::string typeName<unsigned int>() { return "uint32"; }
+template<> const std::string typeName<float>() { return "float32"; }
+template<> const std::string typeName<long>() { return "int64"; }
+template<> const std::string typeName<unsigned long>() { return "uint64"; }
+template<> const std::string typeName<double>() { return "float64"; }
+
+template<typename dist_t, typename data_t>
+inline void register_index_class(py::module &m, std::string className) {
+    py::class_<Index<dist_t, data_t>>(m, className.c_str())
+        .def(py::init(&Index<dist_t, data_t>::createFromParams), py::arg("params"))
            /* WARNING: Index::createFromIndex is not thread-safe with Index::addItems */
-        .def(py::init(&Index<float>::createFromIndex), py::arg("index"))
+        .def(py::init(&Index<dist_t, data_t>::createFromIndex), py::arg("index"))
         .def(py::init<const std::string &, const int>(), py::arg("space"), py::arg("dim"))
-        .def("init_index", &Index<float>::init_new_index, py::arg("max_elements"), py::arg("M")=16, py::arg("ef_construction")=200, py::arg("random_seed")=100)
-        .def("knn_query", &Index<float>::knnQuery_return_numpy, py::arg("data"), py::arg("k")=1, py::arg("num_threads")=-1)
-        .def("add_items", &Index<float>::addItems, py::arg("data"), py::arg("ids") = py::none(), py::arg("num_threads")=-1)
-        .def("get_items", &Index<float, float>::getDataReturnList, py::arg("ids") = py::none())
-        .def("get_ids_list", &Index<float>::getIdsList)
-        .def("set_ef", &Index<float>::set_ef, py::arg("ef"))
-        .def("set_num_threads", &Index<float>::set_num_threads, py::arg("num_threads"))
-        .def("save_index", &Index<float>::saveIndex, py::arg("path_to_index"))
-        .def("load_index", &Index<float>::loadIndex, py::arg("path_to_index"), py::arg("max_elements")=0)
-        .def("mark_deleted", &Index<float>::markDeleted, py::arg("label"))
-        .def("unmark_deleted", &Index<float>::unmarkDeleted, py::arg("label"))
-        .def("resize_index", &Index<float>::resizeIndex, py::arg("new_size"))
-        .def("get_max_elements", &Index<float>::getMaxElements)
-        .def("get_current_count", &Index<float>::getCurrentCount)
-        .def_readonly("space", &Index<float>::space_name)
-        .def_readonly("dim", &Index<float>::dim)
-        .def_readwrite("num_threads", &Index<float>::num_threads_default)
+        .def("init_index", &Index<dist_t, data_t>::init_new_index, py::arg("max_elements"), py::arg("M")=16, py::arg("ef_construction")=200, py::arg("random_seed")=100)
+        .def("knn_query", &Index<dist_t, data_t>::knnQuery_return_numpy, py::arg("data"), py::arg("k")=1, py::arg("num_threads")=-1)
+        .def("add_items", &Index<dist_t, data_t>::addItems, py::arg("data"), py::arg("ids") = py::none(), py::arg("num_threads")=-1)
+        .def("get_items", &Index<dist_t, data_t>::getDataReturnNumpy, py::arg("ids") = py::none())
+        .def("get_ids_list", &Index<dist_t, data_t>::getIdsList)
+        .def("set_ef", &Index<dist_t, data_t>::set_ef, py::arg("ef"))
+        .def("set_num_threads", &Index<dist_t, data_t>::set_num_threads, py::arg("num_threads"))
+        .def("save_index", &Index<dist_t, data_t>::saveIndex, py::arg("path_to_index"))
+        .def("load_index", &Index<dist_t, data_t>::loadIndex, py::arg("path_to_index"), py::arg("max_elements")=0)
+        .def("mark_deleted", &Index<dist_t, data_t>::markDeleted, py::arg("label"))
+        .def("unmark_deleted", &Index<dist_t, data_t>::unmarkDeleted, py::arg("label"))
+        .def("resize_index", &Index<dist_t, data_t>::resizeIndex, py::arg("new_size"))
+        .def("get_max_elements", &Index<dist_t, data_t>::getMaxElements)
+        .def("get_current_count", &Index<dist_t, data_t>::getCurrentCount)
+        .def_readonly("space", &Index<dist_t, data_t>::space_name)
+        .def_readonly("dim", &Index<dist_t, data_t>::dim)
+        .def_readwrite("num_threads", &Index<dist_t, data_t>::num_threads_default)
         .def_property("ef",
-          [](const Index<float> & index) {
+          [](const Index<dist_t, data_t> & index) {
             return index.index_inited ? index.appr_alg->ef_ : index.default_ef;
           },
-          [](Index<float> & index, const size_t ef_) {
+          [](Index<dist_t, data_t> & index, const size_t ef_) {
             index.default_ef=ef_;
             if (index.appr_alg)
               index.appr_alg->ef_ = ef_;
         })
-        .def_property_readonly("max_elements", [](const Index<float> & index) {
+        .def_property_readonly("max_elements", [](const Index<dist_t, data_t> & index) {
             return index.index_inited ? index.appr_alg->max_elements_ : 0;
         })
-        .def_property_readonly("element_count", [](const Index<float> & index) {
+        .def_property_readonly("element_count", [](const Index<dist_t, data_t> & index) {
             return index.index_inited ? index.appr_alg->cur_element_count : 0;
         })
-        .def_property_readonly("ef_construction", [](const Index<float> & index) {
+        .def_property_readonly("ef_construction", [](const Index<dist_t, data_t> & index) {
           return index.index_inited ? index.appr_alg->ef_construction_ : 0;
         })
-        .def_property_readonly("M",  [](const Index<float> & index) {
+        .def_property_readonly("M",  [](const Index<dist_t, data_t> & index) {
           return index.index_inited ? index.appr_alg->M_ : 0;
         })
+        .def_property_readonly_static("distance_type",  [](py::object /* cls */) {
+          return typeName<dist_t>();
+        })
+        .def_property_readonly_static("data_type",  [](py::object /* cls */) {
+          return typeName<data_t>();
+        })
 
         .def(py::pickle(
-            [](const Index<float> &ind) { // __getstate__
+            [](const Index<dist_t, data_t> &ind) { // __getstate__
                 return py::make_tuple(ind.getIndexParams()); /* Return dict (wrapped in a tuple) that fully encodes state of the Index object */
             },
             [](py::tuple t) { // __setstate__
                 if (t.size() != 1)
                     throw std::runtime_error("Invalid state!");
 
-                return Index<float>::createFromParams(t[0].cast<py::dict>());
+                return Index<dist_t, data_t>::createFromParams(t[0].cast<py::dict>());
             }
         ))
 
-        .def("__repr__", [](const Index<float> &a) {
-            return "<hnswlib.Index(space='" + a.space_name + "', dim="+std::to_string(a.dim)+")>";
+        .def("__repr__", [className](const Index<dist_t, data_t> &a) {
+            return "<hnswlib." + className + "(space='" + a.space_name + "', dim=" + std::to_string(a.dim) + ") distance type " + typeName<dist_t>() + ", data type " + typeName<data_t>() + ">";
         });
+};
+
+PYBIND11_PLUGIN(hnswlib) {
+        py::module m("hnswlib");
+
+        register_index_class<float, float>(m, "Index");
+        register_index_class<double, double>(m, "DoubleIndex");
+        
+        // Using an unsigned int as dist_t with char as data_t allows up to 65,536 dimensions before we overflow.
+        register_index_class<unsigned int, char>(m, "Int8Index");
+        register_index_class<unsigned int, unsigned char>(m, "UInt8Index");
+
+        // Using an unsigned int as dist_t with short as data_t allows up to 2^32 dimensions before we overflow.
+        register_index_class<unsigned long, short>(m, "Int16Index");
+        register_index_class<unsigned long, unsigned short>(m, "UInt16Index");
 
         py::class_<BFIndex<float>>(m, "BFIndex")
         .def(py::init<const std::string &, const int>(), py::arg("space"), py::arg("dim"))
diff --git a/python_bindings/tests/bindings_test_getdata.py b/python_bindings/tests/bindings_test_getdata.py
index 2985c1dd..5290da7b 100644
--- a/python_bindings/tests/bindings_test_getdata.py
+++ b/python_bindings/tests/bindings_test_getdata.py
@@ -43,4 +43,4 @@ def testGettingItems(self):
 
         # After adding them, all labels should be retrievable
         returned_items = p.get_items(labels)
-        self.assertSequenceEqual(data.tolist(), returned_items)
+        np.testing.assert_equal(data, returned_items)
diff --git a/python_bindings/tests/bindings_test_labels.py b/python_bindings/tests/bindings_test_labels.py
index 2b091371..54c8d560 100644
--- a/python_bindings/tests/bindings_test_labels.py
+++ b/python_bindings/tests/bindings_test_labels.py
@@ -8,127 +8,144 @@
 
 class RandomSelfTestCase(unittest.TestCase):
     def testRandomSelf(self):
-        for idx in range(2):
-            print("\n**** Index save-load test ****\n")
-
-            np.random.seed(idx)
-            dim = 16
-            num_elements = 10000
-
-            # Generating sample data
-            data = np.float32(np.random.random((num_elements, dim)))
-
-            # Declaring index
-            p = hnswlib.Index(space='l2', dim=dim)  # possible options are l2, cosine or ip
-
-            # Initiating index
-            # max_elements - the maximum number of elements, should be known beforehand
-            #     (probably will be made optional in the future)
-            #
-            # ef_construction - controls index search speed/build speed tradeoff
-            # M - is tightly connected with internal dimensionality of the data
-            #     strongly affects the memory consumption
-
-            p.init_index(max_elements=num_elements, ef_construction=100, M=16)
-
-            # Controlling the recall by setting ef:
-            # higher ef leads to better accuracy, but slower search
-            p.set_ef(100)
-
-            p.set_num_threads(4)  # by default using all available cores
-
-            # We split the data in two batches:
-            data1 = data[:num_elements // 2]
-            data2 = data[num_elements // 2:]
-
-            print("Adding first batch of %d elements" % (len(data1)))
-            p.add_items(data1)
-
-            # Query the elements for themselves and measure recall:
-            labels, distances = p.knn_query(data1, k=1)
-
-            items = p.get_items(labels)
-
-            # Check the recall:
-            self.assertAlmostEqual(np.mean(labels.reshape(-1) == np.arange(len(data1))), 1.0, 3)
-
-            # Check that the returned element data is correct:
-            diff_with_gt_labels=np.mean(np.abs(data1-items))
-            self.assertAlmostEqual(diff_with_gt_labels, 0, delta=1e-4)
-
-            # Serializing and deleting the index.
-            # We need the part to check that serialization is working properly.
-
-            index_path = 'first_half.bin'
-            print("Saving index to '%s'" % index_path)
-            p.save_index(index_path)
-            print("Saved. Deleting...")
-            del p
-            print("Deleted")
-
-            print("\n**** Mark delete test ****\n")
-            # Re-initiating, loading the index
-            print("Re-initiating")
-            p = hnswlib.Index(space='l2', dim=dim)
-
-            print("\nLoading index from '%s'\n" % index_path)
-            p.load_index(index_path)
-            p.set_ef(100)
-
-            print("Adding the second batch of %d elements" % (len(data2)))
-            p.add_items(data2)
-
-            # Query the elements for themselves and measure recall:
-            labels, distances = p.knn_query(data, k=1)
-            items = p.get_items(labels)
-
-            # Check the recall:
-            self.assertAlmostEqual(np.mean(labels.reshape(-1) == np.arange(len(data))), 1.0, 3)
-
-            # Check that the returned element data is correct:
-            diff_with_gt_labels = np.mean(np.abs(data-items))
-            self.assertAlmostEqual(diff_with_gt_labels, 0, delta=1e-4) # deleting index.
-
-            # Checking that all labels are returned correctly:
-            sorted_labels = sorted(p.get_ids_list())
-            self.assertEqual(np.sum(~np.asarray(sorted_labels) == np.asarray(range(num_elements))), 0)
-
-            # Delete data1
-            labels1_deleted, _ = p.knn_query(data1, k=1)
-
-            for l in labels1_deleted:
-                p.mark_deleted(l[0])
-            labels2, _ = p.knn_query(data2, k=1)
-            items = p.get_items(labels2)
-            diff_with_gt_labels = np.mean(np.abs(data2-items))
-            self.assertAlmostEqual(diff_with_gt_labels, 0, delta=1e-3) # console
-
-            labels1_after, _ = p.knn_query(data1, k=1)
-            for la in labels1_after:
-                for lb in labels1_deleted:
-                    if la[0] == lb[0]:
-                        self.assertTrue(False)
-            print("All the data in data1 are removed")
-
-            # Checking saving/loading index with elements marked as deleted
-            del_index_path = "with_deleted.bin"
-            p.save_index(del_index_path)
-            p = hnswlib.Index(space='l2', dim=dim)
-            p.load_index(del_index_path)
-            p.set_ef(100)
-
-            labels1_after, _ = p.knn_query(data1, k=1)
-            for la in labels1_after:
-                for lb in labels1_deleted:
-                    if la[0] == lb[0]:
-                        self.assertTrue(False)
-
-            # Unmark deleted data
-            for l in labels1_deleted:
-                p.unmark_deleted(l[0])
-            labels_restored, _ = p.knn_query(data1, k=1)
-            self.assertAlmostEqual(np.mean(labels_restored.reshape(-1) == np.arange(len(data1))), 1.0, 3)
-            print("All the data in data1 are restored")
-
-        os.remove(index_path)
-        os.remove(del_index_path)
+        index_classes = (
+            hnswlib.Index,
+            hnswlib.DoubleIndex,
+            hnswlib.Int8Index,
+            hnswlib.Int16Index,
+            hnswlib.UInt8Index,
+            hnswlib.UInt16Index,
+        )
+        for index_class in index_classes:
+            for idx in range(2):
+                print("\n**** %s save-load test ****\n" % index_class.__name__)
+
+                np.random.seed(idx)
+                dim = 16
+                num_elements = 10000
+
+                # Generating sample data
+                data = np.random.random((num_elements, dim))
+                if np.issubdtype(index_class.data_type, np.integer):
+                    # Scale the data to fit the integer bounds:
+                    info = np.iinfo(index_class.data_type)
+                    data = (data * (info.max - info.min)) + info.min
+                data = data.astype(index_class.data_type)
+
+                # Declaring index
+                p = index_class(space='l2', dim=dim)  # possible options are l2, cosine or ip
+
+                # Initiating index
+                # max_elements - the maximum number of elements, should be known beforehand
+                #     (probably will be made optional in the future)
+                #
+                # ef_construction - controls index search speed/build speed tradeoff
+                # M - is tightly connected with internal dimensionality of the data
+                #     strongly affects the memory consumption
+
+                p.init_index(max_elements=num_elements, ef_construction=100, M=16)
+
+                # Controlling the recall by setting ef:
+                # higher ef leads to better accuracy, but slower search
+                p.set_ef(100)
+
+                p.set_num_threads(4)  # by default using all available cores
+
+                # We split the data in two batches:
+                data1 = data[:num_elements // 2]
+                data2 = data[num_elements // 2:]
+
+                print("Adding first batch of %d elements" % (len(data1)))
+                p.add_items(data1)
+
+                # Query the elements for themselves and measure recall:
+                labels, distances = p.knn_query(data1)
+                distances = distances.reshape(-1)
+
+                # Check that the distances are indeed monotonically increasing:
+                is_sorted = np.all(np.diff(distances) >= 0)
+                assert is_sorted, ("Expected distances to be sorted in ascending order, but got: " + repr(distances))
+                
+                items = p.get_items(labels)
+
+                # Check the recall:
+                recall = np.mean(labels.reshape(-1) == np.arange(len(data1)))
+                self.assertAlmostEqual(recall, 1.0, 3)
+
+                # Check that the returned element data is correct:
+                np.testing.assert_allclose(data1, items)
+
+                # Serializing and deleting the index.
+                # We need the part to check that serialization is working properly.
+
+                index_path = 'first_half.bin'
+                print("Saving index to '%s'" % index_path)
+                p.save_index(index_path)
+                print("Saved. Deleting...")
+                del p
+                print("Deleted")
+
+                print("\n**** Mark delete test ****\n")
+                # Re-initiating, loading the index
+                print("Re-initiating")
+                p = index_class(space='l2', dim=dim)
+
+                print("\nLoading index from '%s'\n" % index_path)
+                p.load_index(index_path)
+                p.set_ef(100)
+
+                print("Adding the second batch of %d elements" % (len(data2)))
+                p.add_items(data2)
+
+                # Query the elements for themselves and measure recall:
+                labels, distances = p.knn_query(data, k=1)
+                items = p.get_items(labels)
+
+                # Check the recall:
+                self.assertAlmostEqual(np.mean(labels.reshape(-1) == np.arange(len(data))), 1.0, 3)
+
+                # Check that the returned element data is correct:
+                np.testing.assert_allclose(data, items) # deleting index.
+
+                # Checking that all labels are returned correctly:
+                sorted_labels = sorted(p.get_ids_list())
+                self.assertEqual(np.sum(~np.asarray(sorted_labels) == np.asarray(range(num_elements))), 0)
+
+                # Delete data1
+                labels1_deleted, _ = p.knn_query(data1, k=1)
+
+                for l in labels1_deleted:
+                    p.mark_deleted(l[0])
+                labels2, _ = p.knn_query(data2, k=1)
+                items = p.get_items(labels2)
+                np.testing.assert_allclose(data2, items)
+
+                labels1_after, _ = p.knn_query(data1, k=1)
+                for la in labels1_after:
+                    for lb in labels1_deleted:
+                        if la[0] == lb[0]:
+                            self.assertTrue(False)
+                print("All the data in data1 are removed")
+
+                # Checking saving/loading index with elements marked as deleted
+                del_index_path = "with_deleted.bin"
+                p.save_index(del_index_path)
+                p = index_class(space='l2', dim=dim)
+                p.load_index(del_index_path)
+                p.set_ef(100)
+
+                labels1_after, _ = p.knn_query(data1, k=1)
+                for la in labels1_after:
+                    for lb in labels1_deleted:
+                        if la[0] == lb[0]:
+                            self.assertTrue(False)
+
+                # Unmark deleted data
+                for l in labels1_deleted:
+                    p.unmark_deleted(l[0])
+                labels_restored, _ = p.knn_query(data1, k=1)
+                self.assertAlmostEqual(np.mean(labels_restored.reshape(-1) == np.arange(len(data1))), 1.0, 3)
+                print("All the data in data1 are restored")
+
+            os.remove(index_path)
+            os.remove(del_index_path)
diff --git a/sift_1b.cpp b/sift_1b.cpp
index 2739490c..adabbc63 100644
--- a/sift_1b.cpp
+++ b/sift_1b.cpp
@@ -146,7 +146,7 @@ static size_t getCurrentRSS() {
 
 
 static void
-get_gt(unsigned int *massQA, unsigned char *massQ, unsigned char *mass, size_t vecsize, size_t qsize, L2SpaceI &l2space,
+get_gt(unsigned int *massQA, unsigned char *massQ, unsigned char *mass, size_t vecsize, size_t qsize, L2Space<int> &l2space,
        size_t vecdim, vector<std::priority_queue<std::pair<int, labeltype >>> &answers, size_t k) {
 
 
@@ -288,7 +288,7 @@ void sift_test1B() {
     unsigned char *mass = new unsigned char[vecdim];
     ifstream input(path_data, ios::binary);
     int in = 0;
-    L2SpaceI l2space(vecdim);
+    L2Space<int> l2space(vecdim);
 
     HierarchicalNSW<int> *appr_alg;
     if (exists_test(path_index)) {