initial commit

LICENSE

@@ -0,0 +1,21 @@
+MIT License
+
+Copyright (c) 2021 Nicolas Tessore
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.

MANIFEST.in

@@ -0,0 +1,3 @@
+include flt.pyx
+include dctdlt.c
+exclude flt.c

README.md

@@ -0,0 +1,20 @@
+
+flt
+===
+
+**fast Legendre transform**
+
+This is a minimal Python package for fast discrete Legendre transforms (DLTs).
+The implementation uses a recursive version of the matrix relations by Alpert &
+Rokhlin (1991) to compute the DLT via a discrete cosine transform (DCT).
+
+The package can be installed using pip:
+
+    pip install flt
+
+For more information, please see the [documentation].
+
+Current functionality covers the absolutely minimal use case.  Please open an
+issue on GitHub if you would like to see anything added.
+
+[documentation]: https://cltools.readthedocs.io/flt/

dctdlt.c

@@ -0,0 +1,107 @@
+// dctdlt.c
+// ========
+// discrete Legendre transform via DCT
+//
+// author: Nicolas Tessore <[email protected]>
+// license: MIT
+//
+// Synopsis
+// --------
+// The `dctdlt` and `dltdct` functions convert the coefficients of a discrete
+// cosine transform (DCT) to the coefficients of a discrete Legendre transform
+// (DLT) and vice versa [1].
+//
+// References
+// ----------
+// [1] Alpert, B. K., & Rokhlin, V. (1991). A fast algorithm for the evaluation
+//     of Legendre expansions. SIAM Journal on Scientific and Statistical
+//     Computing, 12(1), 158-179.
+//
+
+#define DCTDLT_VERSION 20210318L
+
+
+// dctdlt
+// ======
+// convert DCT coefficients to DLT coefficients
+//
+// Parameters
+// ----------
+// n : unsigned int
+//     Length of the input array.
+// dct : (n,) array of double
+//     Input DCT coefficients.
+// dlt : (n,) array of double, output
+//     Output DLT coefficients.
+//
+void dctdlt(unsigned int n, const double* dct, double* dlt)
+{
+    double a, b;
+    unsigned int k, l;
+
+    // first row
+    a = 1.;
+    b = a;
+    dlt[0] = 0.5*b*dct[0];
+    for(k = 2; k < n; k += 2)
+    {
+        b *= (k-3.)/(k+1.);
+        dlt[0] += b*dct[k];
+    }
+
+    // remaining rows
+    for(l = 1; l < n; ++l)
+    {
+        a /= (1. - 0.5/l);
+        b = a;
+        dlt[l] = b*dct[l];
+        for(k = l+2; k < n; k += 2)
+        {
+            b *= (k*(k+l-2.)*(k-l-3.))/((k-2.)*(k+l+1.)*(k-l));
+            dlt[l] += b*dct[k];
+        }
+    }
+}
+
+
+// dltdct
+// ======
+// convert DLT coefficients to DCT coefficients
+//
+// Parameters
+// ----------
+// n : unsigned int
+//     Length of the input array.
+// dlt : (n,) array of double
+//     Input DLT coefficients.
+// dct : (n,) array of double, output
+//     Output DCT coefficients.
+//
+void dltdct(unsigned int n, const double* dlt, double* dct)
+{
+    double a, b;
+    unsigned int k, l;
+
+    // first row
+    a = 1.;
+    b = a;
+    dct[0] = b*dlt[0];
+    for(l = 2; l < n; l += 2)
+    {
+        b *= ((l-1.)*(l-1.))/(l*l);
+        dct[0] += b*dlt[l];
+    }
+
+    // remaining rows
+    for(k = 1; k < n; ++k)
+    {
+        a *= (1. - 0.5/k);
+        b = a;
+        dct[k] = b*dlt[k];
+        for(l = k+2; l < n; l += 2)
+        {
+            b *= ((l-k-1.)*(l+k-1.))/((l-k)*(l+k));
+            dct[k] += b*dlt[l];
+        }
+    }
+}