Repository: eigency-org/eigency (language: python)

Issue:
Possiblity to use long double as type
Currently there is no possibility to map long double matrices to numpy.
Any idea how to do it?

This would be very beneficial to all users as very often matrix exponentials as computed by Eigen benefit a lot of large precision.



Discussion / hints:
As long as both Eigen and numpy have `long double` support, I don't see a reason why this couldn't be added.
https://numpy.org/doc/stable/reference/arrays.scalars.html#numpy.longdouble

Fix patch:
diff --git a/README.md b/README.md
index 9d37778..9ac6834 100644
--- a/README.md
+++ b/README.md
@@ -497,3 +497,43 @@ def function_filter3(np.ndarray[np.float64_t, ndim=2] array):
 
 In all three cases, the returned array will now be of the same shape
 as the original.
+
+### Long double support
+
+Eigency provides new shorthands for Eigen `long double` and complex `long double` Matrix and Array types.
+Examples:
+```
+Vector4ld
+Matrix3ld
+Vector2cld
+Matrix4cld
+Array3Xld
+ArrayXXcld
+```
+
+These typedefs are available in the `eigency` namespace when including the `eigency` header:
+```c++
+#include "eigency.h"
+
+void receive_long_double_matrix(Eigen::Map<eigency::MatrixXld> &mat) {
+    // use long double eigen matrix
+}
+```
+
+Use Cython (.pyx) to create Python binding to your C++ function:
+```python
+cdef extern from "functions.h":
+     cdef void _receive_long_double_matrix "receive_long_double_matrix"(Map[MatrixXld] &)
+
+def send_long_double_ndarray(np.ndarray[np.longdouble_t, ndim=2] array):
+    return _receive_long_double_matrix(Map[MatrixXld](array))
+```
+
+Invoke in Python:
+```python
+import numpy as np
+import my_module
+
+x = np.array([[1.1, 2.2], [3.3, 4.4]], dtype=np.longdouble)
+my_module.send_long_double_ndarray(x)
+```
diff --git a/eigency/conversions.pxd b/eigency/conversions.pxd
index fcf2d79..b4138b1 100644
--- a/eigency/conversions.pxd
+++ b/eigency/conversions.pxd
@@ -1,6 +1,11 @@
 cimport numpy as np
 
 
+cdef api np.ndarray[long double, ndim=2] ndarray_long_double_C(long double *data, long rows, long cols, long outer_stride, long inner_stride)
+cdef api np.ndarray[long double, ndim=2] ndarray_long_double_F(long double *data, long rows, long cols, long outer_stride, long inner_stride)
+cdef api np.ndarray[long double, ndim=2] ndarray_copy_long_double_C(const long double *data, long rows, long cols, long outer_stride, long inner_stride)
+cdef api np.ndarray[long double, ndim=2] ndarray_copy_long_double_F(const long double *data, long rows, long cols, long outer_stride, long inner_stride)
+
 cdef api np.ndarray[double, ndim=2] ndarray_double_C(double *data, long rows, long cols, long outer_stride, long inner_stride)
 cdef api np.ndarray[double, ndim=2] ndarray_double_F(double *data, long rows, long cols, long outer_stride, long inner_stride)
 cdef api np.ndarray[double, ndim=2] ndarray_copy_double_C(const double *data, long rows, long cols, long outer_stride, long inner_stride)
@@ -51,6 +56,11 @@ cdef api np.ndarray[unsigned char, ndim=2] ndarray_uchar_F(unsigned char *data,
 cdef api np.ndarray[unsigned char, ndim=2] ndarray_copy_uchar_C(const unsigned char *data, long rows, long cols, long outer_stride, long inner_stride)
 cdef api np.ndarray[unsigned char, ndim=2] ndarray_copy_uchar_F(const unsigned char *data, long rows, long cols, long outer_stride, long inner_stride)
 
+cdef api np.ndarray[np.npy_clongdouble, ndim=2] ndarray_complex_long_double_C(np.npy_clongdouble *data, long rows, long cols, long outer_stride, long inner_stride)
+cdef api np.ndarray[np.npy_clongdouble, ndim=2] ndarray_complex_long_double_F(np.npy_clongdouble *data, long rows, long cols, long outer_stride, long inner_stride)
+cdef api np.ndarray[np.npy_clongdouble, ndim=2] ndarray_copy_complex_long_double_C(const np.npy_clongdouble *data, long rows, long cols, long outer_stride, long inner_stride)
+cdef api np.ndarray[np.npy_clongdouble, ndim=2] ndarray_copy_complex_long_double_F(const np.npy_clongdouble *data, long rows, long cols, long outer_stride, long inner_stride)
+
 cdef api np.ndarray[np.complex128_t, ndim=2] ndarray_complex_double_C(np.complex128_t *data, long rows, long cols, long outer_stride, long inner_stride)
 cdef api np.ndarray[np.complex128_t, ndim=2] ndarray_complex_double_F(np.complex128_t *data, long rows, long cols, long outer_stride, long inner_stride)
 cdef api np.ndarray[np.complex128_t, ndim=2] ndarray_copy_complex_double_C(const np.complex128_t *data, long rows, long cols, long outer_stride, long inner_stride)
diff --git a/eigency/conversions.pyx b/eigency/conversions.pyx
index 4093729..8667713 100644
--- a/eigency/conversions.pyx
+++ b/eigency/conversions.pyx
@@ -3,6 +3,41 @@ cimport cython
 import numpy as np
 from numpy.lib.stride_tricks import as_strided
 
+#
+# long double
+#
+
+@cython.boundscheck(False)
+cdef np.ndarray[long double, ndim=2] ndarray_long_double_C(long double *data, long rows, long cols, long row_stride, long col_stride):
+    cdef long double[:,:] mem_view = <long double[:rows,:cols]>data
+    dtype = 'longdouble'
+    cdef int itemsize = np.dtype(dtype).itemsize
+    return as_strided(np.asarray(mem_view, dtype=dtype, order="C"), strides=[row_stride*itemsize, col_stride*itemsize])
+
+@cython.boundscheck(False)
+cdef np.ndarray[long double, ndim=2] ndarray_long_double_F(long double *data, long rows, long cols, long row_stride, long col_stride):
+    cdef long double[::1,:] mem_view = <long double[:rows:1,:cols]>data
+    dtype = 'longdouble'
+    cdef int itemsize = np.dtype(dtype).itemsize
+    return as_strided(np.asarray(mem_view, dtype=dtype, order="F"), strides=[row_stride*itemsize, col_stride*itemsize])
+
+@cython.boundscheck(False)
+cdef np.ndarray[long double, ndim=2] ndarray_copy_long_double_C(const long double *data, long rows, long cols, long row_stride, long col_stride):
+    cdef long double[:,:] mem_view = <long double[:rows,:cols]>data
+    dtype = 'longdouble'
+    cdef int itemsize = np.dtype(dtype).itemsize
+    return np.copy(as_strided(np.asarray(mem_view, dtype=dtype, order="C"), strides=[row_stride*itemsize, col_stride*itemsize]))
+
+@cython.boundscheck(False)
+cdef np.ndarray[long double, ndim=2] ndarray_copy_long_double_F(const long double *data, long rows, long cols, long row_stride, long col_stride):
+    cdef long double[::1,:] mem_view = <long double[:rows:1,:cols]>data
+    dtype = 'longdouble'
+    cdef int itemsize = np.dtype(dtype).itemsize
+    return np.copy(as_strided(np.asarray(mem_view, dtype=dtype, order="F"), strides=[row_stride*itemsize, col_stride*itemsize]))
+
+#
+# double
+#
 
 @cython.boundscheck(False)
 cdef np.ndarray[double, ndim=2] ndarray_double_C(double *data, long rows, long cols, long row_stride, long col_stride):
@@ -10,6 +45,7 @@ cdef np.ndarray[double, ndim=2] ndarray_double_C(double *data, long rows, long c
     dtype = 'double'
     cdef int itemsize = np.dtype(dtype).itemsize
     return as_strided(np.asarray(mem_view, dtype=dtype, order="C"), strides=[row_stride*itemsize, col_stride*itemsize])
+
 @cython.boundscheck(False)
 cdef np.ndarray[double, ndim=2] ndarray_double_F(double *data, long rows, long cols, long row_stride, long col_stride):
     cdef double[::1,:] mem_view = <double[:rows:1,:cols]>data
@@ -23,6 +59,7 @@ cdef np.ndarray[double, ndim=2] ndarray_copy_double_C(const double *data, long r
     dtype = 'double'
     cdef int itemsize = np.dtype(dtype).itemsize
     return np.copy(as_strided(np.asarray(mem_view, dtype=dtype, order="C"), strides=[row_stride*itemsize, col_stride*itemsize]))
+
 @cython.boundscheck(False)
 cdef np.ndarray[double, ndim=2] ndarray_copy_double_F(const double *data, long rows, long cols, long row_stride, long col_stride):
     cdef double[::1,:] mem_view = <double[:rows:1,:cols]>data
@@ -30,6 +67,9 @@ cdef np.ndarray[double, ndim=2] ndarray_copy_double_F(const double *data, long r
     cdef int itemsize = np.dtype(dtype).itemsize
     return np.copy(as_strided(np.asarray(mem_view, dtype=dtype, order="F"), strides=[row_stride*itemsize, col_stride*itemsize]))
 
+#
+# float
+#
 
 @cython.boundscheck(False)
 cdef np.ndarray[float, ndim=2] ndarray_float_C(float *data, long rows, long cols, long row_stride, long col_stride):
@@ -37,6 +77,7 @@ cdef np.ndarray[float, ndim=2] ndarray_float_C(float *data, long rows, long cols
     dtype = 'float'
     cdef int itemsize = np.dtype(dtype).itemsize
     return as_strided(np.asarray(mem_view, dtype=dtype, order="C"), strides=[row_stride*itemsize, col_stride*itemsize])
+
 @cython.boundscheck(False)
 cdef np.ndarray[float, ndim=2] ndarray_float_F(float *data, long rows, long cols, long row_stride, long col_stride):
     cdef float[::1,:] mem_view = <float[:rows:1,:cols]>data
@@ -50,6 +91,7 @@ cdef np.ndarray[float, ndim=2] ndarray_copy_float_C(const float *data, long rows
     dtype = 'float'
     cdef int itemsize = np.dtype(dtype).itemsize
     return np.copy(as_strided(np.asarray(mem_view, dtype=dtype, order="C"), strides=[row_stride*itemsize, col_stride*itemsize]))
+
 @cython.boundscheck(False)
 cdef np.ndarray[float, ndim=2] ndarray_copy_float_F(const float *data, long rows, long cols, long row_stride, long col_stride):
     cdef float[::1,:] mem_view = <float[:rows:1,:cols]>data
@@ -57,6 +99,9 @@ cdef np.ndarray[float, ndim=2] ndarray_copy_float_F(const float *data, long rows
     cdef int itemsize = np.dtype(dtype).itemsize
     return np.copy(as_strided(np.asarray(mem_view, dtype=dtype, order="F"), strides=[row_stride*itemsize, col_stride*itemsize]))
 
+#
+# long
+#
 
 @cython.boundscheck(False)
 cdef np.ndarray[long, ndim=2] ndarray_long_C(long *data, long rows, long cols, long row_stride, long col_stride):
@@ -64,6 +109,7 @@ cdef np.ndarray[long, ndim=2] ndarray_long_C(long *data, long rows, long cols, l
     dtype = 'int_'
     cdef int itemsize = np.dtype(dtype).itemsize
     return as_strided(np.asarray(mem_view, dtype=dtype, order="C"), strides=[row_stride*itemsize, col_stride*itemsize])
+
 @cython.boundscheck(False)
 cdef np.ndarray[long, ndim=2] ndarray_long_F(long *data, long rows, long cols, long row_stride, long col_stride):
     cdef long[::1,:] mem_view = <long[:rows:1,:cols]>data
@@ -77,6 +123,7 @@ cdef np.ndarray[long, ndim=2] ndarray_copy_long_C(const long *data, long rows, l
     dtype = 'int_'
     cdef int itemsize = np.dtype(dtype).itemsize
     return np.copy(as_strided(np.asarray(mem_view, dtype=dtype, order="C"), strides=[row_stride*itemsize, col_stride*itemsize]))
+
 @cython.boundscheck(False)
 cdef np.ndarray[long, ndim=2] ndarray_copy_long_F(const long *data, long rows, long cols, long row_stride, long col_stride):
     cdef long[::1,:] mem_view = <long[:rows:1,:cols]>data
@@ -84,6 +131,9 @@ cdef np.ndarray[long, ndim=2] ndarray_copy_long_F(const long *data, long rows, l
     cdef int itemsize = np.dtype(dtype).itemsize
     return np.copy(as_strided(np.asarray(mem_view, dtype=dtype, order="F"), strides=[row_stride*itemsize, col_stride*itemsize]))
 
+#
+# unsigned long
+#
 
 @cython.boundscheck(False)
 cdef np.ndarray[unsigned long, ndim=2] ndarray_ulong_C(unsigned long *data, long rows, long cols, long row_stride, long col_stride):
@@ -91,6 +141,7 @@ cdef np.ndarray[unsigned long, ndim=2] ndarray_ulong_C(unsigned long *data, long
     dtype = 'uint'
     cdef int itemsize = np.dtype(dtype).itemsize
     return as_strided(np.asarray(mem_view, dtype=dtype, order="C"), strides=[row_stride*itemsize, col_stride*itemsize])
+
 @cython.boundscheck(False)
 cdef np.ndarray[unsigned long, ndim=2] ndarray_ulong_F(unsigned long *data, long rows, long cols, long row_stride, long col_stride):
     cdef unsigned long[::1,:] mem_view = <unsigned long[:rows:1,:cols]>data
@@ -104,6 +155,7 @@ cdef np.ndarray[unsigned long, ndim=2] ndarray_copy_ulong_C(const unsigned long
     dtype = 'uint'
     cdef int itemsize = np.dtype(dtype).itemsize
     return np.copy(as_strided(np.asarray(mem_view, dtype=dtype, order="C"), strides=[row_stride*itemsize, col_stride*itemsize]))
+
 @cython.boundscheck(False)
 cdef np.ndarray[unsigned long, ndim=2] ndarray_copy_ulong_F(const unsigned long *data, long rows, long cols, long row_stride, long col_stride):
     cdef unsigned long[::1,:] mem_view = <unsigned long[:rows:1,:cols]>data
@@ -111,6 +163,9 @@ cdef np.ndarray[unsigned long, ndim=2] ndarray_copy_ulong_F(const unsigned long
     cdef int itemsize = np.dtype(dtype).itemsize
     return np.copy(as_strided(np.asarray(mem_view, dtype=dtype, order="F"), strides=[row_stride*itemsize, col_stride*itemsize]))
 
+#
+# int
+#
 
 @cython.boundscheck(False)
 cdef np.ndarray[int, ndim=2] ndarray_int_C(int *data, long rows, long cols, long row_stride, long col_stride):
@@ -118,6 +173,7 @@ cdef np.ndarray[int, ndim=2] ndarray_int_C(int *data, long rows, long cols, long
     dtype = 'int'
     cdef int itemsize = np.dtype(dtype).itemsize
     return as_strided(np.asarray(mem_view, dtype=dtype, order="C"), strides=[row_stride*itemsize, col_stride*itemsize])
+
 @cython.boundscheck(False)
 cdef np.ndarray[int, ndim=2] ndarray_int_F(int *data, long rows, long cols, long row_stride, long col_stride):
     cdef int[::1,:] mem_view = <int[:rows:1,:cols]>data
@@ -131,6 +187,7 @@ cdef np.ndarray[int, ndim=2] ndarray_copy_int_C(const int *data, long rows, long
     dtype = 'int'
     cdef int itemsize = np.dtype(dtype).itemsize
     return np.copy(as_strided(np.asarray(mem_view, dtype=dtype, order="C"), strides=[row_stride*itemsize, col_stride*itemsize]))
+
 @cython.boundscheck(False)
 cdef np.ndarray[int, ndim=2] ndarray_copy_int_F(const int *data, long rows, long cols, long row_stride, long col_stride):
     cdef int[::1,:] mem_view = <int[:rows:1,:cols]>data
@@ -138,6 +195,9 @@ cdef np.ndarray[int, ndim=2] ndarray_copy_int_F(const int *data, long rows, long
     cdef int itemsize = np.dtype(dtype).itemsize
     return np.copy(as_strided(np.asarray(mem_view, dtype=dtype, order="F"), strides=[row_stride*itemsize, col_stride*itemsize]))
 
+#
+# unsigned int
+#
 
 @cython.boundscheck(False)
 cdef np.ndarray[unsigned int, ndim=2] ndarray_uint_C(unsigned int *data, long rows, long cols, long row_stride, long col_stride):
@@ -145,6 +205,7 @@ cdef np.ndarray[unsigned int, ndim=2] ndarray_uint_C(unsigned int *data, long ro
     dtype = 'uint'
     cdef int itemsize = np.dtype(dtype).itemsize
     return as_strided(np.asarray(mem_view, dtype=dtype, order="C"), strides=[row_stride*itemsize, col_stride*itemsize])
+
 @cython.boundscheck(False)
 cdef np.ndarray[unsigned int, ndim=2] ndarray_uint_F(unsigned int *data, long rows, long cols, long row_stride, long col_stride):
     cdef unsigned int[::1,:] mem_view = <unsigned int[:rows:1,:cols]>data
@@ -158,6 +219,7 @@ cdef np.ndarray[unsigned int, ndim=2] ndarray_copy_uint_C(const unsigned int *da
     dtype = 'uint'
     cdef int itemsize = np.dtype(dtype).itemsize
     return np.copy(as_strided(np.asarray(mem_view, dtype=dtype, order="C"), strides=[row_stride*itemsize, col_stride*itemsize]))
+
 @cython.boundscheck(False)
 cdef np.ndarray[unsigned int, ndim=2] ndarray_copy_uint_F(const unsigned int *data, long rows, long cols, long row_stride, long col_stride):
     cdef unsigned int[::1,:] mem_view = <unsigned int[:rows:1,:cols]>data
@@ -165,6 +227,9 @@ cdef np.ndarray[unsigned int, ndim=2] ndarray_copy_uint_F(const unsigned int *da
     cdef int itemsize = np.dtype(dtype).itemsize
     return np.copy(as_strided(np.asarray(mem_view, dtype=dtype, order="F"), strides=[row_stride*itemsize, col_stride*itemsize]))
 
+#
+# short
+#
 
 @cython.boundscheck(False)
 cdef np.ndarray[short, ndim=2] ndarray_short_C(short *data, long rows, long cols, long row_stride, long col_stride):
@@ -172,6 +237,7 @@ cdef np.ndarray[short, ndim=2] ndarray_short_C(short *data, long rows, long cols
     dtype = 'short'
     cdef int itemsize = np.dtype(dtype).itemsize
     return as_strided(np.asarray(mem_view, dtype=dtype, order="C"), strides=[row_stride*itemsize, col_stride*itemsize])
+
 @cython.boundscheck(False)
 cdef np.ndarray[short, ndim=2] ndarray_short_F(short *data, long rows, long cols, long row_stride, long col_stride):
     cdef short[::1,:] mem_view = <short[:rows:1,:cols]>data
@@ -185,6 +251,7 @@ cdef np.ndarray[short, ndim=2] ndarray_copy_short_C(const short *data, long rows
     dtype = 'short'
     cdef int itemsize = np.dtype(dtype).itemsize
     return np.copy(as_strided(np.asarray(mem_view, dtype=dtype, order="C"), strides=[row_stride*itemsize, col_stride*itemsize]))
+
 @cython.boundscheck(False)
 cdef np.ndarray[short, ndim=2] ndarray_copy_short_F(const short *data, long rows, long cols, long row_stride, long col_stride):
     cdef short[::1,:] mem_view = <short[:rows:1,:cols]>data
@@ -192,6 +259,9 @@ cdef np.ndarray[short, ndim=2] ndarray_copy_short_F(const short *data, long rows
     cdef int itemsize = np.dtype(dtype).itemsize
     return np.copy(as_strided(np.asarray(mem_view, dtype=dtype, order="F"), strides=[row_stride*itemsize, col_stride*itemsize]))
 
+#
+# unsigned short
+#
 
 @cython.boundscheck(False)
 cdef np.ndarray[unsigned short, ndim=2] ndarray_ushort_C(unsigned short *data, long rows, long cols, long row_stride, long col_stride):
@@ -199,6 +269,7 @@ cdef np.ndarray[unsigned short, ndim=2] ndarray_ushort_C(unsigned short *data, l
     dtype = 'ushort'
     cdef int itemsize = np.dtype(dtype).itemsize
     return as_strided(np.asarray(mem_view, dtype=dtype, order="C"), strides=[row_stride*itemsize, col_stride*itemsize])
+
 @cython.boundscheck(False)
 cdef np.ndarray[unsigned short, ndim=2] ndarray_ushort_F(unsigned short *data, long rows, long cols, long row_stride, long col_stride):
     cdef unsigned short[::1,:] mem_view = <unsigned short[:rows:1,:cols]>data
@@ -212,6 +283,7 @@ cdef np.ndarray[unsigned short, ndim=2] ndarray_copy_ushort_C(const unsigned sho
     dtype = 'ushort'
     cdef int itemsize = np.dtype(dtype).itemsize
     return np.copy(as_strided(np.asarray(mem_view, dtype=dtype, order="C"), strides=[row_stride*itemsize, col_stride*itemsize]))
+
 @cython.boundscheck(False)
 cdef np.ndarray[unsigned short, ndim=2] ndarray_copy_ushort_F(const unsigned short *data, long rows, long cols, long row_stride, long col_stride):
     cdef unsigned short[::1,:] mem_view = <unsigned short[:rows:1,:cols]>data
@@ -219,6 +291,9 @@ cdef np.ndarray[unsigned short, ndim=2] ndarray_copy_ushort_F(const unsigned sho
     cdef int itemsize = np.dtype(dtype).itemsize
     return np.copy(as_strided(np.asarray(mem_view, dtype=dtype, order="F"), strides=[row_stride*itemsize, col_stride*itemsize]))
 
+#
+# signed char
+#
 
 @cython.boundscheck(False)
 cdef np.ndarray[signed char, ndim=2] ndarray_schar_C(signed char *data, long rows, long cols, long row_stride, long col_stride):
@@ -226,6 +301,7 @@ cdef np.ndarray[signed char, ndim=2] ndarray_schar_C(signed char *data, long row
     dtype = 'int8'
     cdef int itemsize = np.dtype(dtype).itemsize
     return as_strided(np.asarray(mem_view, dtype=dtype, order="C"), strides=[row_stride*itemsize, col_stride*itemsize])
+
 @cython.boundscheck(False)
 cdef np.ndarray[signed char, ndim=2] ndarray_schar_F(signed char *data, long rows, long cols, long row_stride, long col_stride):
     cdef signed char[::1,:] mem_view = <signed char[:rows:1,:cols]>data
@@ -239,6 +315,7 @@ cdef np.ndarray[signed char, ndim=2] ndarray_copy_schar_C(const signed char *dat
     dtype = 'int8'
     cdef int itemsize = np.dtype(dtype).itemsize
     return np.copy(as_strided(np.asarray(mem_view, dtype=dtype, order="C"), strides=[row_stride*itemsize, col_stride*itemsize]))
+
 @cython.boundscheck(False)
 cdef np.ndarray[signed char, ndim=2] ndarray_copy_schar_F(const signed char *data, long rows, long cols, long row_stride, long col_stride):
     cdef signed char[::1,:] mem_view = <signed char[:rows:1,:cols]>data
@@ -246,6 +323,9 @@ cdef np.ndarray[signed char, ndim=2] ndarray_copy_schar_F(const signed char *dat
     cdef int itemsize = np.dtype(dtype).itemsize
     return np.copy(as_strided(np.asarray(mem_view, dtype=dtype, order="F"), strides=[row_stride*itemsize, col_stride*itemsize]))
 
+#
+# unsigned char
+#
 
 @cython.boundscheck(False)
 cdef np.ndarray[unsigned char, ndim=2] ndarray_uchar_C(unsigned char *data, long rows, long cols, long row_stride, long col_stride):
@@ -253,6 +333,7 @@ cdef np.ndarray[unsigned char, ndim=2] ndarray_uchar_C(unsigned char *data, long
     dtype = 'uint8'
     cdef int itemsize = np.dtype(dtype).itemsize
     return as_strided(np.asarray(mem_view, dtype=dtype, order="C"), strides=[row_stride*itemsize, col_stride*itemsize])
+
 @cython.boundscheck(False)
 cdef np.ndarray[unsigned char, ndim=2] ndarray_uchar_F(unsigned char *data, long rows, long cols, long row_stride, long col_stride):
     cdef unsigned char[::1,:] mem_view = <unsigned char[:rows:1,:cols]>data
@@ -266,6 +347,7 @@ cdef np.ndarray[unsigned char, ndim=2] ndarray_copy_uchar_C(const unsigned char
     dtype = 'uint8'
     cdef int itemsize = np.dtype(dtype).itemsize
     return np.copy(as_strided(np.asarray(mem_view, dtype=dtype, order="C"), strides=[row_stride*itemsize, col_stride*itemsize]))
+
 @cython.boundscheck(False)
 cdef np.ndarray[unsigned char, ndim=2] ndarray_copy_uchar_F(const unsigned char *data, long rows, long cols, long row_stride, long col_stride):
     cdef unsigned char[::1,:] mem_view = <unsigned char[:rows:1,:cols]>data
@@ -273,6 +355,41 @@ cdef np.ndarray[unsigned char, ndim=2] ndarray_copy_uchar_F(const unsigned char
     cdef int itemsize = np.dtype(dtype).itemsize
     return np.copy(as_strided(np.asarray(mem_view, dtype=dtype, order="F"), strides=[row_stride*itemsize, col_stride*itemsize]))
 
+#
+# complex long double
+#
+
+@cython.boundscheck(False)
+cdef np.ndarray[np.npy_clongdouble, ndim=2] ndarray_complex_long_double_C(np.npy_clongdouble *data, long rows, long cols, long row_stride, long col_stride):
+    cdef np.npy_clongdouble[:,:] mem_view = <np.npy_clongdouble[:rows,:cols]>data
+    dtype = 'clongdouble'
+    cdef int itemsize = np.dtype(dtype).itemsize
+    return as_strided(np.asarray(mem_view, dtype=dtype, order="C"), strides=[row_stride*itemsize, col_stride*itemsize])
+
+@cython.boundscheck(False)
+cdef np.ndarray[np.npy_clongdouble, ndim=2] ndarray_complex_long_double_F(np.npy_clongdouble *data, long rows, long cols, long row_stride, long col_stride):
+    cdef np.npy_clongdouble[::1,:] mem_view = <np.npy_clongdouble[:rows:1,:cols]>data
+    dtype = 'clongdouble'
+    cdef int itemsize = np.dtype(dtype).itemsize
+    return as_strided(np.asarray(mem_view, dtype=dtype, order="F"), strides=[row_stride*itemsize, col_stride*itemsize])
+
+@cython.boundscheck(False)
+cdef np.ndarray[np.npy_clongdouble, ndim=2] ndarray_copy_complex_long_double_C(const np.npy_clongdouble *data, long rows, long cols, long row_stride, long col_stride):
+    cdef np.npy_clongdouble[:,:] mem_view = <np.npy_clongdouble[:rows,:cols]>data
+    dtype = 'clongdouble'
+    cdef int itemsize = np.dtype(dtype).itemsize
+    return np.copy(as_strided(np.asarray(mem_view, dtype=dtype, order="C"), strides=[row_stride*itemsize, col_stride*itemsize]))
+
+@cython.boundscheck(False)
+cdef np.ndarray[np.npy_clongdouble, ndim=2] ndarray_copy_complex_long_double_F(const np.npy_clongdouble *data, long rows, long cols, long row_stride, long col_stride):
+    cdef np.npy_clongdouble[::1,:] mem_view = <np.npy_clongdouble[:rows:1,:cols]>data
+    dtype = 'clongdouble'
+    cdef int itemsize = np.dtype(dtype).itemsize
+    return np.copy(as_strided(np.asarray(mem_view, dtype=dtype, order="F"), strides=[row_stride*itemsize, col_stride*itemsize]))
+
+#
+# complex double
+#
 
 @cython.boundscheck(False)
 cdef np.ndarray[np.complex128_t, ndim=2] ndarray_complex_double_C(np.complex128_t *data, long rows, long cols, long row_stride, long col_stride):
@@ -280,6 +397,7 @@ cdef np.ndarray[np.complex128_t, ndim=2] ndarray_complex_double_C(np.complex128_
     dtype = 'complex128'
     cdef int itemsize = np.dtype(dtype).itemsize
     return as_strided(np.asarray(mem_view, dtype=dtype, order="C"), strides=[row_stride*itemsize, col_stride*itemsize])
+
 @cython.boundscheck(False)
 cdef np.ndarray[np.complex128_t, ndim=2] ndarray_complex_double_F(np.complex128_t *data, long rows, long cols, long row_stride, long col_stride):
     cdef np.complex128_t[::1,:] mem_view = <np.complex128_t[:rows:1,:cols]>data
@@ -293,6 +411,7 @@ cdef np.ndarray[np.complex128_t, ndim=2] ndarray_copy_complex_double_C(const np.
     dtype = 'complex128'
     cdef int itemsize = np.dtype(dtype).itemsize
     return np.copy(as_strided(np.asarray(mem_view, dtype=dtype, order="C"), strides=[row_stride*itemsize, col_stride*itemsize]))
+
 @cython.boundscheck(False)
 cdef np.ndarray[np.complex128_t, ndim=2] ndarray_copy_complex_double_F(const np.complex128_t *data, long rows, long cols, long row_stride, long col_stride):
     cdef np.complex128_t[::1,:] mem_view = <np.complex128_t[:rows:1,:cols]>data
@@ -300,6 +419,9 @@ cdef np.ndarray[np.complex128_t, ndim=2] ndarray_copy_complex_double_F(const np.
     cdef int itemsize = np.dtype(dtype).itemsize
     return np.copy(as_strided(np.asarray(mem_view, dtype=dtype, order="F"), strides=[row_stride*itemsize, col_stride*itemsize]))
 
+#
+# complex float
+#
 
 @cython.boundscheck(False)
 cdef np.ndarray[np.complex64_t, ndim=2] ndarray_complex_float_C(np.complex64_t *data, long rows, long cols, long row_stride, long col_stride):
@@ -307,6 +429,7 @@ cdef np.ndarray[np.complex64_t, ndim=2] ndarray_complex_float_C(np.complex64_t *
     dtype = 'complex64'
     cdef int itemsize = np.dtype(dtype).itemsize
     return as_strided(np.asarray(mem_view, dtype=dtype, order="C"), strides=[row_stride*itemsize, col_stride*itemsize])
+
 @cython.boundscheck(False)
 cdef np.ndarray[np.complex64_t, ndim=2] ndarray_complex_float_F(np.complex64_t *data, long rows, long cols, long row_stride, long col_stride):
     cdef np.complex64_t[::1,:] mem_view = <np.complex64_t[:rows:1,:cols]>data
@@ -320,6 +443,7 @@ cdef np.ndarray[np.complex64_t, ndim=2] ndarray_copy_complex_float_C(const np.co
     dtype = 'complex64'
     cdef int itemsize = np.dtype(dtype).itemsize
     return np.copy(as_strided(np.asarray(mem_view, dtype=dtype, order="C"), strides=[row_stride*itemsize, col_stride*itemsize]))
+
 @cython.boundscheck(False)
 cdef np.ndarray[np.complex64_t, ndim=2] ndarray_copy_complex_float_F(const np.complex64_t *data, long rows, long cols, long row_stride, long col_stride):
     cdef np.complex64_t[::1,:] mem_view = <np.complex64_t[:rows:1,:cols]>data
diff --git a/eigency/core.pxd b/eigency/core.pxd
index 8205c20..d54df04 100644
--- a/eigency/core.pxd
+++ b/eigency/core.pxd
@@ -3,6 +3,7 @@ cimport numpy as np
 
 ctypedef signed char schar;
 ctypedef unsigned char uchar;
+ctypedef long double long_double;
 
 ctypedef fused dtype:
     uchar
@@ -12,6 +13,7 @@ ctypedef fused dtype:
     long
     float
     double
+    long_double
 
 ctypedef fused DenseType:
     Matrix
@@ -102,6 +104,18 @@ ctypedef fused DenseTypeShort:
     Matrix3d
     Matrix4d
     MatrixXd
+    Vector2ld
+    Vector3ld
+    Vector4ld
+    VectorXld
+    RowVector2ld
+    RowVector3ld
+    RowVector4ld
+    RowVectorXld
+    Matrix2ld
+    Matrix3ld
+    Matrix4ld
+    MatrixXld
     Vector1cf
     Vector2cf
     Vector3cf
@@ -132,6 +146,18 @@ ctypedef fused DenseTypeShort:
     Matrix3cd
     Matrix4cd
     MatrixXcd
+    Vector2cld
+    Vector3cld
+    Vector4cld
+    VectorXcld
+    RowVector2cld
+    RowVector3cld
+    RowVector4cld
+    RowVectorXcld
+    Matrix2cld
+    Matrix3cld
+    Matrix4cld
+    MatrixXcld
     Array22i
     Array23i
     Array24i
@@ -192,6 +218,26 @@ ctypedef fused DenseTypeShort:
     Array3d
     Array4d
     ArrayXd
+    Array22ld
+    Array23ld
+    Array24ld
+    Array2Xld
+    Array32ld
+    Array33ld
+    Array34ld
+    Array3Xld
+    Array42ld
+    Array43ld
+    Array44ld
+    Array4Xld
+    ArrayX2ld
+    ArrayX3ld
+    ArrayX4ld
+    ArrayXXld
+    Array2ld
+    Array3ld
+    Array4ld
+    ArrayXld
     Array22cf
     Array23cf
     Array24cf
@@ -232,6 +278,26 @@ ctypedef fused DenseTypeShort:
     Array3cd
     Array4cd
     ArrayXcd
+    Array22cld
+    Array23cld
+    Array24cld
+    Array2Xcld
+    Array32cld
+    Array33cld
+    Array34cld
+    Array3Xcld
+    Array42cld
+    Array43cld
+    Array44cld
+    Array4Xcld
+    ArrayX2cld
+    ArrayX3cld
+    ArrayX4cld
+    ArrayXXcld
+    Array2cld
+    Array3cld
+    Array4cld
+    ArrayXcld
 
 ctypedef fused StorageOrder:
     RowMajor
@@ -386,9 +452,6 @@ cdef extern from "eigency_cpp.h" namespace "Eigen":
      cdef cppclass Array(PlainObjectBase):
           pass
 
-     cdef cppclass VectorXd(PlainObjectBase):
-          pass
-
      cdef cppclass Vector1i(PlainObjectBase):
           pass
 
@@ -913,3 +976,197 @@ cdef extern from "eigency_cpp.h" namespace "Eigen":
 
      cdef cppclass ArrayXcd(PlainObjectBase):
           pass
+
+cdef extern from "eigency.h" namespace "eigency":
+
+     cdef cppclass Vector2ld(PlainObjectBase):
+          pass
+
+     cdef cppclass Vector3ld(PlainObjectBase):
+          pass
+
+     cdef cppclass Vector4ld(PlainObjectBase):
+          pass
+
+     cdef cppclass VectorXld(PlainObjectBase):
+          pass
+
+     cdef cppclass RowVector2ld(PlainObjectBase):
+          pass
+
+     cdef cppclass RowVector3ld(PlainObjectBase):
+          pass
+
+     cdef cppclass RowVector4ld(PlainObjectBase):
+          pass
+
+     cdef cppclass RowVectorXld(PlainObjectBase):
+          pass
+
+     cdef cppclass Matrix2ld(PlainObjectBase):
+          pass
+
+     cdef cppclass Matrix3ld(PlainObjectBase):
+          pass
+
+     cdef cppclass Matrix4ld(PlainObjectBase):
+          pass
+
+     cdef cppclass MatrixXld(PlainObjectBase):
+          pass
+
+     cdef cppclass Vector2cld(PlainObjectBase):
+          pass
+
+     cdef cppclass Vector3cld(PlainObjectBase):
+          pass
+
+     cdef cppclass Vector4cld(PlainObjectBase):
+          pass
+
+     cdef cppclass VectorXcld(PlainObjectBase):
+          pass
+
+     cdef cppclass RowVector2cld(PlainObjectBase):
+          pass
+
+     cdef cppclass RowVector3cld(PlainObjectBase):
+          pass
+
+     cdef cppclass RowVector4cld(PlainObjectBase):
+          pass
+
+     cdef cppclass RowVectorXcld(PlainObjectBase):
+          pass
+
+     cdef cppclass Matrix2cld(PlainObjectBase):
+          pass
+
+     cdef cppclass Matrix3cld(PlainObjectBase):
+          pass
+
+     cdef cppclass Matrix4cld(PlainObjectBase):
+          pass
+
+     cdef cppclass MatrixXcld(PlainObjectBase):
+          pass
+
+     cdef cppclass Array22ld(PlainObjectBase):
+          pass
+
+     cdef cppclass Array23ld(PlainObjectBase):
+          pass
+
+     cdef cppclass Array24ld(PlainObjectBase):
+          pass
+
+     cdef cppclass Array2Xld(PlainObjectBase):
+          pass
+
+     cdef cppclass Array32ld(PlainObjectBase):
+          pass
+
+     cdef cppclass Array33ld(PlainObjectBase):
+          pass
+
+     cdef cppclass Array34ld(PlainObjectBase):
+          pass
+
+     cdef cppclass Array3Xld(PlainObjectBase):
+          pass
+
+     cdef cppclass Array42ld(PlainObjectBase):
+          pass
+
+     cdef cppclass Array43ld(PlainObjectBase):
+          pass
+
+     cdef cppclass Array44ld(PlainObjectBase):
+          pass
+
+     cdef cppclass Array4Xld(PlainObjectBase):
+          pass
+
+     cdef cppclass ArrayX2ld(PlainObjectBase):
+          pass
+
+     cdef cppclass ArrayX3ld(PlainObjectBase):
+          pass
+
+     cdef cppclass ArrayX4ld(PlainObjectBase):
+          pass
+
+     cdef cppclass ArrayXXld(PlainObjectBase):
+          pass
+
+     cdef cppclass Array2ld(PlainObjectBase):
+          pass
+
+     cdef cppclass Array3ld(PlainObjectBase):
+          pass
+
+     cdef cppclass Array4ld(PlainObjectBase):
+          pass
+
+     cdef cppclass ArrayXld(PlainObjectBase):
+          pass
+
+     cdef cppclass Array22cld(PlainObjectBase):
+          pass
+
+     cdef cppclass Array23cld(PlainObjectBase):
+          pass
+
+     cdef cppclass Array24cld(PlainObjectBase):
+          pass
+
+     cdef cppclass Array2Xcld(PlainObjectBase):
+          pass
+
+     cdef cppclass Array32cld(PlainObjectBase):
+          pass
+
+     cdef cppclass Array33cld(PlainObjectBase):
+          pass
+
+     cdef cppclass Array34cld(PlainObjectBase):
+          pass
+
+     cdef cppclass Array3Xcld(PlainObjectBase):
+          pass
+
+     cdef cppclass Array42cld(PlainObjectBase):
+          pass
+
+     cdef cppclass Array43cld(PlainObjectBase):
+          pass
+
+     cdef cppclass Array44cld(PlainObjectBase):
+          pass
+
+     cdef cppclass Array4Xcld(PlainObjectBase):
+          pass
+
+     cdef cppclass ArrayX2cld(PlainObjectBase):
+          pass
+
+     cdef cppclass ArrayX3cld(PlainObjectBase):
+          pass
+
+     cdef cppclass ArrayX4cld(PlainObjectBase):
+          pass
+
+     cdef cppclass ArrayXXcld(PlainObjectBase):
+          pass
+
+     cdef cppclass Array2cld(PlainObjectBase):
+          pass
+
+     cdef cppclass Array3cld(PlainObjectBase):
+          pass
+
+     cdef cppclass Array4cld(PlainObjectBase):
+          pass
+
+     cdef cppclass ArrayXcld(PlainObjectBase):
+          pass
diff --git a/eigency/eigency.h b/eigency/eigency.h
new file mode 100644
index 0000000..f93e261
--- /dev/null
+++ b/eigency/eigency.h
@@ -0,0 +1,59 @@
+#ifndef EIGENCY_H
+#define EIGENCY_H
+
+#include <Eigen/Core>
+
+namespace eigency {
+
+//
+// Matrix
+//
+
+#define EIGENCY_MAKE_TYPEDEFS(Type, TypeSuffix, Size, SizeSuffix)           \
+typedef ::Eigen::Matrix<Type, Size, Size> Matrix##SizeSuffix##TypeSuffix;   \
+typedef ::Eigen::Matrix<Type, Size, 1>    Vector##SizeSuffix##TypeSuffix;   \
+typedef ::Eigen::Matrix<Type, 1, Size>    RowVector##SizeSuffix##TypeSuffix;
+
+#define EIGENCY_MAKE_TYPEDEFS_ALL_SIZES(Type, TypeSuffix)   \
+EIGENCY_MAKE_TYPEDEFS(Type, TypeSuffix, 2, 2)               \
+EIGENCY_MAKE_TYPEDEFS(Type, TypeSuffix, 3, 3)               \
+EIGENCY_MAKE_TYPEDEFS(Type, TypeSuffix, 4, 4)               \
+EIGENCY_MAKE_TYPEDEFS(Type, TypeSuffix, ::Eigen::Dynamic, X)
+
+EIGENCY_MAKE_TYPEDEFS_ALL_SIZES(long double,                 ld)
+EIGENCY_MAKE_TYPEDEFS_ALL_SIZES(::std::complex<long double>, cld)
+
+#undef EIGENCY_MAKE_TYPEDEFS_ALL_SIZES
+#undef EIGENCY_MAKE_TYPEDEFS
+
+//
+// Array
+//
+
+#define EIGENCY_MAKE_ARRAY_TYPEDEFS(Type, TypeSuffix, Size, SizeSuffix)             \
+typedef ::Eigen::Array<Type, Size, Size> Array##SizeSuffix##SizeSuffix##TypeSuffix; \
+typedef ::Eigen::Array<Type, Size, 1>    Array##SizeSuffix##TypeSuffix;
+
+#define EIGENCY_MAKE_ARRAY_FIXED_TYPEDEFS(Type, TypeSuffix, Size)                   \
+typedef ::Eigen::Array<Type, Size, ::Eigen::Dynamic> Array##Size##X##TypeSuffix;    \
+typedef ::Eigen::Array<Type, ::Eigen::Dynamic, Size> Array##X##Size##TypeSuffix;
+
+#define EIGENCY_MAKE_ARRAY_TYPEDEFS_ALL_SIZES(Type, TypeSuffix)     \
+EIGENCY_MAKE_ARRAY_TYPEDEFS(Type, TypeSuffix, 2, 2)                 \
+EIGENCY_MAKE_ARRAY_TYPEDEFS(Type, TypeSuffix, 3, 3)                 \
+EIGENCY_MAKE_ARRAY_TYPEDEFS(Type, TypeSuffix, 4, 4)                 \
+EIGENCY_MAKE_ARRAY_TYPEDEFS(Type, TypeSuffix, ::Eigen::Dynamic, X)  \
+EIGENCY_MAKE_ARRAY_FIXED_TYPEDEFS(Type, TypeSuffix, 2)              \
+EIGENCY_MAKE_ARRAY_FIXED_TYPEDEFS(Type, TypeSuffix, 3)              \
+EIGENCY_MAKE_ARRAY_FIXED_TYPEDEFS(Type, TypeSuffix, 4)
+
+EIGENCY_MAKE_ARRAY_TYPEDEFS_ALL_SIZES(long double,                 ld)
+EIGENCY_MAKE_ARRAY_TYPEDEFS_ALL_SIZES(::std::complex<long double>, cld)
+
+#undef EIGENCY_MAKE_ARRAY_TYPEDEFS_ALL_SIZES
+#undef EIGENCY_MAKE_ARRAY_FIXED_TYPEDEFS
+#undef EIGENCY_MAKE_ARRAY_TYPEDEFS
+
+}
+
+#endif
diff --git a/eigency/eigency_cpp.h b/eigency/eigency_cpp.h
index 9224585..b0474e1 100644
--- a/eigency/eigency_cpp.h
+++ b/eigency/eigency_cpp.h
@@ -1,25 +1,26 @@
-#include <Eigen/Dense>
-
 #include <iostream>
 #include <stdexcept>
 #include <complex>
 
-typedef ::std::complex< double > __pyx_t_double_complex;
-typedef ::std::complex< float > __pyx_t_float_complex;
+#ifndef EIGENCY_CPP_H
+#define EIGENCY_CPP_H
 
+#include <Eigen/Core>
 #include <numpy/ndarraytypes.h>
 
-#include "conversions_api.h"
+#include "eigency.h"
 
-#ifndef EIGENCY_CPP
-#define EIGENCY_CPP
+typedef ::std::complex< double > __pyx_t_double_complex;
+typedef ::std::complex< float > __pyx_t_float_complex;
+
+#include "conversions_api.h"
 
 namespace eigency {
 
 template<typename Scalar>
-inline PyArrayObject *_ndarray_view(Scalar *, long rows, long cols, bool is_row_major, long outer_stride=0, long inner_stride=0);
+inline PyArrayObject* _ndarray_view(Scalar *, long rows, long cols, bool is_row_major, long outer_stride=0, long inner_stride=0);
 template<typename Scalar>
-inline PyArrayObject *_ndarray_copy(const Scalar *, long rows, long cols, bool is_row_major, long outer_stride=0, long inner_stride=0);
+inline PyArrayObject* _ndarray_copy(const Scalar *, long rows, long cols, bool is_row_major, long outer_stride=0, long inner_stride=0);
 
 // Strides:
 // Eigen and numpy differ in their way of dealing with strides. Eigen has the concept of outer and
@@ -29,191 +30,64 @@ inline PyArrayObject *_ndarray_copy(const Scalar *, long rows, long cols, bool i
 // In contrast, numpy's stride is simply a measure of how fast each dimension should be incremented.
 // Consequently, a switch in numpy storage order from row-major to column-major involves a switch
 // in strides, while it does not affect the stride in Eigen.
-template<>
-inline PyArrayObject *_ndarray_view<double>(double *data, long rows, long cols, bool is_row_major, long outer_stride, long inner_stride) {
-    if (is_row_major) {
-        // Eigen row-major mode: row_stride=outer_stride, and col_stride=inner_stride
-        // If no stride is given, the row_stride is set to the number of columns.
-        return ndarray_double_C(data, rows, cols, outer_stride>0?outer_stride:cols, inner_stride>0?inner_stride:1);
-    } else {
-        // Eigen column-major mode: row_stride=outer_stride, and col_stride=inner_stride
-        // If no stride is given, the cow_stride is set to the number of rows.
-        return ndarray_double_F(data, rows, cols, inner_stride>0?inner_stride:1, outer_stride>0?outer_stride:rows);
-    }
-}
-template<>
-inline PyArrayObject *_ndarray_copy<double>(const double *data, long rows, long cols, bool is_row_major, long outer_stride, long inner_stride) {
-    if (is_row_major)
-        return ndarray_copy_double_C(data, rows, cols, outer_stride>0?outer_stride:cols, inner_stride>0?inner_stride:1);
-    else
-        return ndarray_copy_double_F(data, rows, cols, inner_stride>0?inner_stride:1, outer_stride>0?outer_stride:rows);
-}
-
-template<>
-inline PyArrayObject *_ndarray_view<float>(float *data, long rows, long cols, bool is_row_major, long outer_stride, long inner_stride) {
-    if (is_row_major)
-        return ndarray_float_C(data, rows, cols, outer_stride>0?outer_stride:cols, inner_stride>0?inner_stride:1);
-    else
-        return ndarray_float_F(data, rows, cols, inner_stride>0?inner_stride:1, outer_stride>0?outer_stride:rows);
-}
-template<>
-inline PyArrayObject *_ndarray_copy<float>(const float *data, long rows, long cols, bool is_row_major, long outer_stride, long inner_stride) {
-    if (is_row_major)
-        return ndarray_copy_float_C(data, rows, cols, outer_stride>0?outer_stride:cols, inner_stride>0?inner_stride:1);
-    else
-        return ndarray_copy_float_F(data, rows, cols, inner_stride>0?inner_stride:1, outer_stride>0?outer_stride:rows);
-}
-
-template<>
-inline PyArrayObject *_ndarray_view<long>(long *data, long rows, long cols, bool is_row_major, long outer_stride, long inner_stride) {
-    if (is_row_major)
-        return ndarray_long_C(data, rows, cols, outer_stride>0?outer_stride:cols, inner_stride>0?inner_stride:1);
-    else
-        return ndarray_long_F(data, rows, cols, inner_stride>0?inner_stride:1, outer_stride>0?outer_stride:rows);
-}
-template<>
-inline PyArrayObject *_ndarray_copy<long>(const long *data, long rows, long cols, bool is_row_major, long outer_stride, long inner_stride) {
-    if (is_row_major)
-        return ndarray_copy_long_C(data, rows, cols, outer_stride>0?outer_stride:cols, inner_stride>0?inner_stride:1);
-    else
-        return ndarray_copy_long_F(data, rows, cols, inner_stride>0?inner_stride:1, outer_stride>0?outer_stride:rows);
-}
-
-template<>
-inline PyArrayObject *_ndarray_view<unsigned long>(unsigned long *data, long rows, long cols, bool is_row_major, long outer_stride, long inner_stride) {
-    if (is_row_major)
-        return ndarray_ulong_C(data, rows, cols, outer_stride>0?outer_stride:cols, inner_stride>0?inner_stride:1);
-    else
-        return ndarray_ulong_F(data, rows, cols, inner_stride>0?inner_stride:1, outer_stride>0?outer_stride:rows);
-}
-template<>
-inline PyArrayObject *_ndarray_copy<unsigned long>(const unsigned long *data, long rows, long cols, bool is_row_major, long outer_stride, long inner_stride) {
-    if (is_row_major)
-        return ndarray_copy_ulong_C(data, rows, cols, outer_stride>0?outer_stride:cols, inner_stride>0?inner_stride:1);
-    else
-        return ndarray_copy_ulong_F(data, rows, cols, inner_stride>0?inner_stride:1, outer_stride>0?outer_stride:rows);
-}
-
-template<>
-inline PyArrayObject *_ndarray_view<int>(int *data, long rows, long cols, bool is_row_major, long outer_stride, long inner_stride) {
-    if (is_row_major)
-        return ndarray_int_C(data, rows, cols, outer_stride>0?outer_stride:cols, inner_stride>0?inner_stride:1);
-    else
-        return ndarray_int_F(data, rows, cols, inner_stride>0?inner_stride:1, outer_stride>0?outer_stride:rows);
-}
-template<>
-inline PyArrayObject *_ndarray_copy<int>(const int *data, long rows, long cols, bool is_row_major, long outer_stride, long inner_stride) {
-    if (is_row_major)
-        return ndarray_copy_int_C(data, rows, cols, outer_stride>0?outer_stride:cols, inner_stride>0?inner_stride:1);
-    else
-        return ndarray_copy_int_F(data, rows, cols, inner_stride>0?inner_stride:1, outer_stride>0?outer_stride:rows);
-}
-
-template<>
-inline PyArrayObject *_ndarray_view<unsigned int>(unsigned int *data, long rows, long cols, bool is_row_major, long outer_stride, long inner_stride) {
-    if (is_row_major)
-        return ndarray_uint_C(data, rows, cols, outer_stride>0?outer_stride:cols, inner_stride>0?inner_stride:1);
-    else
-        return ndarray_uint_F(data, rows, cols, inner_stride>0?inner_stride:1, outer_stride>0?outer_stride:rows);
-}
-template<>
-inline PyArrayObject *_ndarray_copy<unsigned int>(const unsigned int *data, long rows, long cols, bool is_row_major, long outer_stride, long inner_stride) {
-    if (is_row_major)
-        return ndarray_copy_uint_C(data, rows, cols, outer_stride>0?outer_stride:cols, inner_stride>0?inner_stride:1);
-    else
-        return ndarray_copy_uint_F(data, rows, cols, inner_stride>0?inner_stride:1, outer_stride>0?outer_stride:rows);
-}
-
-template<>
-inline PyArrayObject *_ndarray_view<short>(short *data, long rows, long cols, bool is_row_major, long outer_stride, long inner_stride) {
-    if (is_row_major)
-        return ndarray_short_C(data, rows, cols, outer_stride>0?outer_stride:cols, inner_stride>0?inner_stride:1);
-    else
-        return ndarray_short_F(data, rows, cols, inner_stride>0?inner_stride:1, outer_stride>0?outer_stride:rows);
-}
-template<>
-inline PyArrayObject *_ndarray_copy<short>(const short *data, long rows, long cols, bool is_row_major, long outer_stride, long inner_stride) {
-    if (is_row_major)
-        return ndarray_copy_short_C(data, rows, cols, outer_stride>0?outer_stride:cols, inner_stride>0?inner_stride:1);
-    else
-        return ndarray_copy_short_F(data, rows, cols, inner_stride>0?inner_stride:1, outer_stride>0?outer_stride:rows);
-}
-
-template<>
-inline PyArrayObject *_ndarray_view<unsigned short>(unsigned short *data, long rows, long cols, bool is_row_major, long outer_stride, long inner_stride) {
-    if (is_row_major)
-        return ndarray_ushort_C(data, rows, cols, outer_stride>0?outer_stride:cols, inner_stride>0?inner_stride:1);
-    else
-        return ndarray_ushort_F(data, rows, cols, inner_stride>0?inner_stride:1, outer_stride>0?outer_stride:rows);
-}
-template<>
-inline PyArrayObject *_ndarray_copy<unsigned short>(const unsigned short *data, long rows, long cols, bool is_row_major, long outer_stride, long inner_stride) {
-    if (is_row_major)
-        return ndarray_copy_ushort_C(data, rows, cols, outer_stride>0?outer_stride:cols, inner_stride>0?inner_stride:1);
-    else
-        return ndarray_copy_ushort_F(data, rows, cols, inner_stride>0?inner_stride:1, outer_stride>0?outer_stride:rows);
-}
-
-template<>
-inline PyArrayObject *_ndarray_view<signed char>(signed char *data, long rows, long cols, bool is_row_major, long outer_stride, long inner_stride) {
-    if (is_row_major)
-        return ndarray_schar_C(data, rows, cols, outer_stride>0?outer_stride:cols, inner_stride>0?inner_stride:1);
-    else
-        return ndarray_schar_F(data, rows, cols, inner_stride>0?inner_stride:1, outer_stride>0?outer_stride:rows);
-}
-template<>
-inline PyArrayObject *_ndarray_copy<signed char>(const signed char *data, long rows, long cols, bool is_row_major, long outer_stride, long inner_stride) {
-    if (is_row_major)
-        return ndarray_copy_schar_C(data, rows, cols, outer_stride>0?outer_stride:cols, inner_stride>0?inner_stride:1);
-    else
-        return ndarray_copy_schar_F(data, rows, cols, inner_stride>0?inner_stride:1, outer_stride>0?outer_stride:rows);
-}
-
-template<>
-inline PyArrayObject *_ndarray_view<unsigned char>(unsigned char *data, long rows, long cols, bool is_row_major, long outer_stride, long inner_stride) {
-    if (is_row_major)
-        return ndarray_uchar_C(data, rows, cols, outer_stride>0?outer_stride:cols, inner_stride>0?inner_stride:1);
-    else
-        return ndarray_uchar_F(data, rows, cols, inner_stride>0?inner_stride:1, outer_stride>0?outer_stride:rows);
-}
-template<>
-inline PyArrayObject *_ndarray_copy<unsigned char>(const unsigned char *data, long rows, long cols, bool is_row_major, long outer_stride, long inner_stride) {
-    if (is_row_major)
-        return ndarray_copy_uchar_C(data, rows, cols, outer_stride>0?outer_stride:cols, inner_stride>0?inner_stride:1);
-    else
-        return ndarray_copy_uchar_F(data, rows, cols, inner_stride>0?inner_stride:1, outer_stride>0?outer_stride:rows);
-}
-
-template<>
-inline PyArrayObject *_ndarray_view<std::complex<double> >(std::complex<double> *data, long rows, long cols, bool is_row_major, long outer_stride, long inner_stride) {
-    if (is_row_major)
-        return ndarray_complex_double_C(data, rows, cols, outer_stride>0?outer_stride:cols, inner_stride>0?inner_stride:1);
-    else
-        return ndarray_complex_double_F(data, rows, cols, inner_stride>0?inner_stride:1, outer_stride>0?outer_stride:rows);
-}
-template<>
-inline PyArrayObject *_ndarray_copy<std::complex<double> >(const std::complex<double> *data, long rows, long cols, bool is_row_major, long outer_stride, long inner_stride) {
-    if (is_row_major)
-        return ndarray_copy_complex_double_C(data, rows, cols, outer_stride>0?outer_stride:cols, inner_stride>0?inner_stride:1);
-    else
-        return ndarray_copy_complex_double_F(data, rows, cols, inner_stride>0?inner_stride:1, outer_stride>0?outer_stride:rows);
-}
-
-template<>
-inline PyArrayObject *_ndarray_view<std::complex<float> >(std::complex<float> *data, long rows, long cols, bool is_row_major, long outer_stride, long inner_stride) {
-    if (is_row_major)
-        return ndarray_complex_float_C(data, rows, cols, outer_stride>0?outer_stride:cols, inner_stride>0?inner_stride:1);
-    else
-        return ndarray_complex_float_F(data, rows, cols, inner_stride>0?inner_stride:1, outer_stride>0?outer_stride:rows);
-}
-template<>
-inline PyArrayObject *_ndarray_copy<std::complex<float> >(const std::complex<float> *data, long rows, long cols, bool is_row_major, long outer_stride, long inner_stride) {
-    if (is_row_major)
-        return ndarray_copy_complex_float_C(data, rows, cols, outer_stride>0?outer_stride:cols, inner_stride>0?inner_stride:1);
-    else
-        return ndarray_copy_complex_float_F(data, rows, cols, inner_stride>0?inner_stride:1, outer_stride>0?outer_stride:rows);
-}
 
+#define _NDAV(TYPE, FUNC_NAME_C, FUNC_NAME_F) template<>                                                                                    \
+inline PyArrayObject* _ndarray_view< TYPE >(TYPE *data, long rows, long cols, bool is_row_major, long outer_stride, long inner_stride) {    \
+    if (is_row_major) {                                                                                                                     \
+        /* Eigen row-major mode: row_stride=outer_stride, and col_stride=inner_stride */                                                    \
+        /* If no stride is given, the row_stride is set to the number of columns. */                                                        \
+        return FUNC_NAME_C(data, rows, cols, outer_stride>0?outer_stride:cols, inner_stride>0?inner_stride:1);                              \
+    } else {                                                                                                                                \
+        /* Eigen column-major mode: row_stride=outer_stride, and col_stride=inner_stride */                                                 \
+        /* If no stride is given, the cow_stride is set to the number of rows. */                                                           \
+        return FUNC_NAME_F(data, rows, cols, inner_stride>0?inner_stride:1, outer_stride>0?outer_stride:rows);                              \
+    }                                                                                                                                       \
+}
+
+#define _NDAC(TYPE, FUNC_NAME_C, FUNC_NAME_F) template<>                                                                                        \
+inline PyArrayObject* _ndarray_copy< TYPE >(const TYPE *data, long rows, long cols, bool is_row_major, long outer_stride, long inner_stride) {  \
+    if (is_row_major) {                                                                                                                         \
+        return FUNC_NAME_C(data, rows, cols, outer_stride>0?outer_stride:cols, inner_stride>0?inner_stride:1);                                  \
+    } else {                                                                                                                                    \
+        return FUNC_NAME_F(data, rows, cols, inner_stride>0?inner_stride:1, outer_stride>0?outer_stride:rows);                                  \
+    }                                                                                                                                           \
+}
+
+_NDAV(long double, ndarray_long_double_C, ndarray_long_double_F)
+_NDAC(long double, ndarray_copy_long_double_C, ndarray_copy_long_double_F)
+_NDAV(double, ndarray_double_C, ndarray_double_F)
+_NDAC(double, ndarray_copy_double_C, ndarray_copy_double_F)
+_NDAV(float, ndarray_float_C, ndarray_float_F)
+_NDAC(float, ndarray_copy_float_C, ndarray_copy_float_F)
+_NDAV(long, ndarray_long_C, ndarray_long_F)
+_NDAC(long, ndarray_copy_long_C, ndarray_copy_long_F)
+_NDAV(unsigned long, ndarray_ulong_C, ndarray_ulong_F)
+_NDAC(unsigned long, ndarray_copy_ulong_C, ndarray_copy_ulong_F)
+_NDAV(int, ndarray_int_C, ndarray_int_F)
+_NDAC(int, ndarray_copy_int_C, ndarray_copy_int_F)
+_NDAV(unsigned int, ndarray_uint_C, ndarray_uint_F)
+_NDAC(unsigned int, ndarray_copy_uint_C, ndarray_copy_uint_F)
+_NDAV(short, ndarray_short_C, ndarray_short_F)
+_NDAC(short, ndarray_copy_short_C, ndarray_copy_short_F)
+_NDAV(unsigned short, ndarray_ushort_C, ndarray_ushort_F)
+_NDAC(unsigned short, ndarray_copy_ushort_C, ndarray_copy_ushort_F)
+_NDAV(signed char, ndarray_schar_C, ndarray_schar_F)
+_NDAC(signed char, ndarray_copy_schar_C, ndarray_copy_schar_F)
+_NDAV(unsigned char, ndarray_uchar_C, ndarray_uchar_F)
+_NDAC(unsigned char, ndarray_copy_uchar_C, ndarray_copy_uchar_F)
+_NDAV(npy_clongdouble, ndarray_complex_long_double_C, ndarray_complex_long_double_F)
+_NDAC(npy_clongdouble, ndarray_copy_complex_long_double_C, ndarray_copy_complex_long_double_F)
+_NDAV(std::complex<double>, ndarray_complex_double_C, ndarray_complex_double_F)
+_NDAC(std::complex<double>, ndarray_copy_complex_double_C, ndarray_copy_complex_double_F)
+_NDAV(std::complex<float>, ndarray_complex_float_C, ndarray_complex_float_F)
+_NDAC(std::complex<float>, ndarray_copy_complex_float_C, ndarray_copy_complex_float_F)
+
+#undef _NDAV
+#undef _NDAC
+
+//
+// Constructors
+//
 
 template <typename Derived>
 inline PyArrayObject *ndarray(Eigen::PlainObjectBase<Derived> &m) {


Test patch:
diff --git a/tests/eigency_tests/eigency_tests.pyx b/tests/eigency_tests/eigency_tests.pyx
index 90b88c9..d3d6922 100644
--- a/tests/eigency_tests/eigency_tests.pyx
+++ b/tests/eigency_tests/eigency_tests.pyx
@@ -22,8 +22,12 @@ cdef extern from "eigency_tests/eigency_tests_cpp.h":
 
      cdef void _function_w_mat_arg "function_w_mat_arg"(Map[MatrixXd] &)
 
+     cdef void _function_w_ld_mat_arg "function_w_ld_mat_arg"(Map[MatrixXld] &)
+
      cdef void _function_w_complex_mat_arg "function_w_complex_mat_arg"(Map[MatrixXcd] &)
 
+     cdef void _function_w_complex_ld_mat_arg "function_w_complex_ld_mat_arg"(Map[MatrixXcld] &)
+
      cdef void _function_w_fullspec_arg "function_w_fullspec_arg" (FlattenedMap[Array, double, Dynamic, _1] &)
 
      cdef VectorXd _function_w_vec_retval "function_w_vec_retval" ()
@@ -88,10 +92,18 @@ def function_w_vec_arg_no_map2(np.ndarray[np.float64_t] array):
 def function_w_mat_arg(np.ndarray[np.float64_t, ndim=2] array):
     return _function_w_mat_arg(Map[MatrixXd](array))
 
+# Function with long double matrix argument.
+def function_w_ld_mat_arg(np.ndarray[np.longdouble_t, ndim=2] array):
+    return _function_w_ld_mat_arg(Map[MatrixXld](array))
+
 # Function with complex matrix argument.
 def function_w_complex_mat_arg(np.ndarray[np.complex128_t, ndim=2] array):
     return _function_w_complex_mat_arg(Map[MatrixXcd](array))
 
+# Function with complex long double matrix argument.
+def function_w_complex_ld_mat_arg(np.ndarray[np.clongdouble_t, ndim=2] array):
+    return _function_w_complex_ld_mat_arg(Map[MatrixXcld](array))
+
 # Function using a full Map specification, rather than the convenience typedefs
 # Note that since cython does not support nested fused types, the Map has been
 # flattened to include all arguments at once
diff --git a/tests/eigency_tests/eigency_tests_cpp.cpp b/tests/eigency_tests/eigency_tests_cpp.cpp
index 2c5b504..77cf80b 100644
--- a/tests/eigency_tests/eigency_tests_cpp.cpp
+++ b/tests/eigency_tests/eigency_tests_cpp.cpp
@@ -1,5 +1,7 @@
 #include "eigency_tests_cpp.h"
+
 #include <iostream>
+#include <complex>
 
 long function_w_vec_arg(Eigen::Map<Eigen::VectorXd> &vec) {
     vec[0] = 0.;
@@ -25,8 +27,18 @@ void function_w_mat_arg(Eigen::Map<Eigen::MatrixXd> &mat) {
     mat(0,0) = 0.;
 }
 
+void function_w_ld_mat_arg(Eigen::Map<eigency::MatrixXld> &mat) {
+    long double ld = 0.0;
+    mat(0,0) = ld;
+}
+
 void function_w_complex_mat_arg(Eigen::Map<Eigen::MatrixXcd> &mat) {
-    mat(0,0) = 0.+0.j;
+    mat(0,0) = std::complex<double>(0.,0.);
+}
+
+void function_w_complex_ld_mat_arg(Eigen::Map<eigency::MatrixXcld> &mat) {
+    std::complex<long double> cld = std::complex<long double>(0.0, 0.0);
+    mat(0,0) = cld;
 }
 
 void function_w_fullspec_arg(Eigen::Map<Eigen::Array<double, Eigen::Dynamic, 1> > &vec) {
diff --git a/tests/eigency_tests/eigency_tests_cpp.h b/tests/eigency_tests/eigency_tests_cpp.h
index 34734de..b0e3e39 100644
--- a/tests/eigency_tests/eigency_tests_cpp.h
+++ b/tests/eigency_tests/eigency_tests_cpp.h
@@ -1,7 +1,8 @@
 #ifndef TEST1_H
 #define TEST1_H
 
-#include "Eigen/Dense"
+#include "Eigen/Core"
+#include "eigency.h"
 
 long function_w_vec_arg(Eigen::Map<Eigen::VectorXd> &vec);
 
@@ -13,8 +14,12 @@ void function_w_vec_arg_no_map2(const Eigen::VectorXd &vec);
 
 void function_w_mat_arg(Eigen::Map<Eigen::MatrixXd> &mat);
 
+void function_w_ld_mat_arg(Eigen::Map<eigency::MatrixXld> &mat);
+
 void function_w_complex_mat_arg(Eigen::Map<Eigen::MatrixXcd> &mat);
 
+void function_w_complex_ld_mat_arg(Eigen::Map<eigency::MatrixXcld> &mat);
+
 void function_w_fullspec_arg(Eigen::Map<Eigen::Array<double, Eigen::Dynamic, 1> > &vec);
 
 Eigen::VectorXd function_w_vec_retval();
diff --git a/tests/run_tests.py b/tests/run_tests.py
index a8b1384..e7e294f 100644
--- a/tests/run_tests.py
+++ b/tests/run_tests.py
@@ -38,12 +38,24 @@ def test_function_w_mat_arg(self):
         # Shared memory test: Verify that first entry was set to 0 by C++ code.
         self.assertAlmostEqual(x[0], 0.0)
 
+    def test_function_w_ld_mat_arg(self):
+        x = np.array([1.1, 2.2, 3.3, 4.4], dtype=np.longdouble)
+        eigency_tests.function_w_ld_mat_arg(x.reshape([2, 2]))
+        # Shared memory test: Verify that first entry was set to 0 by C++ code.
+        self.assertAlmostEqual(x[0], 0.0)
+
     def test_function_w_complex_mat_arg(self):
         x = np.array([1.0 + 1.0j, 2.0 + 1.0j, 3.0 + 1.0j, 4.0 + 1.0j])
         eigency_tests.function_w_complex_mat_arg(x.reshape([2, 2]))
         # Shared memory test: Verify that first entry was set to 0.+0.j by C++ code.
         self.assertAlmostEqual(x[0], 0.0 + 0.0j)
 
+    def test_function_w_complex_ld_mat_arg(self):
+        x = np.array([1.0 + 1.0j, 2.0 + 1.0j, 3.0 + 1.0j, 4.0 + 1.0j], dtype=np.clongdouble)
+        eigency_tests.function_w_complex_ld_mat_arg(x.reshape([2, 2]))
+        # Shared memory test: Verify that first entry was set to 0.+0.j by C++ code.
+        self.assertAlmostEqual(x[0], 0.0 + 0.0j)
+
     def test_funcion_w_fullspec_arg(self):
         x = np.array([1.0, 2.0, 3.0, 4.0])
         eigency_tests.function_w_fullspec_arg(x)