fracmem

Fractional-order derivatives as a fixed-size recursive filter. Fit it once in Python, then run it on a microcontroller with constant compute and constant memory per sample, for as long as the device stays on.

Why it exists

A fractional derivative of order α has memory. The value at time t depends on the whole past of the signal, with weights that fall off slowly, like j^(-α-1). Computed exactly, every new sample costs a little more than the last, and memory grows without limit. That rules it out on a flight controller, a motor driver, or anything else that runs for hours on a few kilobytes of RAM.

fracmem splits the history in two:

• The latest L samples go through an exact FIR filter.
• Everything older goes through p exponential modes ("leaky buckets"). Each mode costs one multiply-add per sample.

The decay rates of the modes come from a Gamma-function identity, not from data. Only the p readout weights are fitted, by cross-validated ridge regression on a few example signals. The result is O(L+p) compute and O(L+p) memory per sample, however long the signal runs.

Explore the docs

At a glance

import numpy as np
from fracmem import CompressedFractionalFilter

train = [np.random.randn(3000).cumsum() * 0.01 for _ in range(8)]

f = CompressedFractionalFilter(alpha=0.5, h=0.01, L=32, p=16)
f.fit(train, j_max=10_000)       # once, offline
y = f.predict(signal)            # batch

f.reset_stream()
y_k = f.step(x_k)                # or one sample at a time

#include "fracmemfilter.h"

filtSetup();                        /* once at startup */
float y = fracmemStep(&filt, x_k);  /* every sample, forever */

Who it is for

• Control engineers running fractional-order PID (PI^λD^μ) or CRONE controllers on real hardware.
• Anyone modelling viscoelastic, diffusive or battery dynamics who needs a fractional derivative inside a real-time loop.
• Researchers who want a fast, checked reference for the Grünwald–Letnikov, Riemann–Liouville or Caputo derivative on long signals.

Installation

fracmem is a pure-Python package on PyPI. The C runtime ships inside it, so a single install covers both fitting and deployment.

Requirements

Install the package

pip install fracmem

python -m venv fracmem-env
source fracmem-env/bin/activate      # Windows: fracmem-env\Scripts\activate
pip install fracmem

git clone https://github.com/edith-lang/fracmem.git
cd fracmem
pip install -e ".[dev]"              # editable install plus pytest

Check the install

python -c "import fracmem; print(fracmem.__version__)"

0.2.0

To check that the bundled C runtime is present, print its folder:

python -c "from fracmem.embedded.export_c import _C_DIR; print(_C_DIR)"

That folder holds fracmemfilter.h and fracmemfilter.c. You copy these two files into your firmware project once. See Embedded C.

Install a C compiler (optional)

You only need a desktop compiler to run the C round-trip example and the C tests. Your firmware toolchain does the real build.

sudo apt install build-essential

xcode-select --install

Development install and tests

git clone https://github.com/edith-lang/fracmem.git
cd fracmem
python -m venv .venv && source .venv/bin/activate
pip install -e ".[dev]"
pytest tests/ -v

......................                                    [100%]
22 passed in 1.40s

Quickstart

Fit a half-derivative filter on eight training signals, run it on a signal 15 times longer, and compare it with the exact computation.

1. Install
   bash

   pip install fracmem

2. Collect training signals

   Use a handful of signals that look like what the filter will see in service. Eight to 32 is usually enough. They must all be the same length. Here they are random walks.

   python

   import numpy as np
   from fracmem import CompressedFractionalFilter
   
   rng = np.random.default_rng(0)
   train_signals = [np.cumsum(rng.standard_normal(3000)) * 0.01 for _ in range(8)]

3. Fit the filter

   alpha is the order of the derivative, h is the sample period in seconds, L is the exact window and p is the number of tail modes. j_max is how far back, in samples, the tail is designed to reach.

   python

   filt = CompressedFractionalFilter(alpha=0.5, h=0.01, L=32, p=16)
   filt.fit(train_signals, j_max=10_000)
   print(f"fitted: L={filt.L}, p={filt.p}, ridge reg={filt.reg_used:.1e}")

4. Run it and compare with the exact answer
   python

   from fracmem import full_gl_derivative, gl_weights
   
   test_signal = np.cumsum(rng.standard_normal(50_000)) * 0.01
   y_hat = filt.predict(test_signal)
   
   w = gl_weights(0.5, len(test_signal) + 1)
   y_exact = full_gl_derivative(test_signal, 0.5, 0.01, w)
   
   rmse = np.sqrt(np.mean((y_hat - y_exact) ** 2))
   scale = np.sqrt(np.mean(y_exact ** 2))
   print(f"RMSE vs exact: {rmse:.6f} ({100 * rmse / scale:.3f}% relative)")

   This is examples/basic-usage.py. On our machine it prints:

   output

   fitted: L=32, p=16, ridge reg=1.0e-10
   
   signal length:      50,000 samples
   compressed predict:  83.81 ms  (O(L+p) per sample, fixed)
   exact reference:     157.09 ms  (O(n) per sample, grows forever)
   RMSE vs exact:       0.002226  (signal RMS: 0.193502, 1.150% relative)

5. Stream one sample at a time

   In a real-time loop, call step() instead of predict(). It gives the same numbers.

   python

   filt.reset_stream()
   for x_k in sensor_samples():
       d_half = filt.step(x_k)

Next steps

• Export the filter to C and run it on a board.
• Pick L, p and j_max for your accuracy and RAM budget.
• Build a filter with no training data when you have none yet.

Fitting a filter

How fit() works, what data it needs, and which parts of the filter the data can and cannot change.

What gets fitted

A fitted filter has three sets of numbers:

Because lam is fixed before any data is seen, fitting c is a linear least-squares problem. It is convex and has one solution, so there is no random initialization and no local minimum to land in.

Choosing training signals

• Look like deployment. Use logs from the real sensor, or a simulation of it. Random walks are a good default for slowly drifting signals such as positions and integrated rates.
• Same length. fit() stacks the signals into one array, so every signal must have the same number of samples.
• Enough of them. Cross-validation uses up to 4 folds over the signals. Use at least 4 signals, and 8 to 32 in practice.
• Longer than L. The tail only starts after the first L samples. Signals of a few thousand samples work well.

The j_max argument

j_max sets the slowest decay rate. The slowest mode is designed to still matter about j_max samples back. If you leave it out, it defaults to the training signal length.

Regularization

fit() picks the ridge strength by k-fold cross-validation across your training signals. It tries each value in reg_candidates, keeps the one with the lowest held-out error, then refits on all signals. The chosen value is stored in filt.reg_used.

filt.fit(train_signals, j_max=10_000,
         reg_candidates=(1e-10, 1e-9, 1e-8, 1e-7))
print(filt.reg_used)

The default grid is 1e-11 to 1e-5. The penalty is scaled by each mode's own energy, so one value works across modes of very different size.

What the fit minimizes

For each training signal, fit() computes the exact derivative with the full convolution, subtracts the exact local term, and fits c so the modes explain what is left:

This minimizes end-to-end error on the kind of signal you train on. That is why a fitted filter often beats the purely analytic one on drifting signals at small L and p.

Streaming, one sample at a time

predict() processes a whole array. step() processes one sample and keeps its state between calls, which is what a control loop needs.

The streaming API

filt.reset_stream()          # clear the window and the modes
for x_k in samples:
    y_k = filt.step(x_k)     # one estimate per sample

step() calls reset_stream() for you the first time. Call it yourself when you start a new run, for example after re-arming a vehicle.

Same numbers as predict()

The tests check that streaming every sample of a signal gives the same output as predict() on the whole signal, for both the "gl" and "caputo" definitions. The C runtime is checked against the same reference.

What one step does

1. For Caputo only: remember the first sample x0, and subtract it from every sample.
2. Push x_k into the length-L window. The sample that falls out, x_{k-L}, feeds the tail.
3. Local term: dot product of the window with w, times h^-α.
4. Tail: for each mode, m_i = lam_i * m_i + x_{k-L}. Then sum c_i * m_i.
5. Return local + tail.

That is L + p multiply-adds, and nothing is allocated.

Example: a fixed-rate loop

import time
from fracmem import CompressedFractionalFilter

filt = CompressedFractionalFilter.load("half_derivative.npz")
filt.reset_stream()

period = filt.h
next_t = time.monotonic()
while running:
    x_k = read_sensor()
    u = kp * x_k + kd * filt.step(x_k)   # fractional PD term
    write_actuator(u)
    next_t += period
    time.sleep(max(0.0, next_t - time.monotonic()))

GL, RL and Caputo definitions

Pick the definition with the definition argument. All three use the same filter; Caputo only changes the input.

Why Caputo matters in practice

If a sensor reads 1.0 at power-on, the GL/RL derivative sees a step from 0 to 1.0 at t = 0. That step leaves a t^-α transient that decays slowly. Caputo subtracts the first sample, so a constant input gives zero output, as you would expect from a derivative.

Usage

filt = CompressedFractionalFilter(alpha=0.5, h=0.01, L=32, p=16,
                                  definition="caputo")
filt.fit(train_signals, j_max=10_000)

The definition is stored by save() and baked into the C export. On the device, caputo = 1 makes fracmemStep record the first sample after each fracmemReset.

Filters without training data

from_analytic() builds a ready-to-deploy filter from math alone. Use it to flash a device before you have any logs.

Usage

from fracmem import CompressedFractionalFilter

filt = CompressedFractionalFilter.from_analytic(alpha=0.5, h=0.01, L=32, p=16,
                                                j_max=20_000)
y_hat = filt.predict(signal)

Both lam and c come from the sum-of-exponentials construction. c is fitted to the exact kernel weights, not to any signal. The result works with predict(), step(), save() and export_c() like a fitted filter.

The two design knobs

tail_fit_points int = 400

How many log-spaced points along the tail the kernel fit checks. More points give a more careful fit at a one-time cost.

tail_fit_reg float = 0.1

How strongly the kernel fit is kept away from extreme weights. Lower trusts the fit more.

j_max int = max(2000, 100·L)

How far back the design reaches. Set it to cover the longest history you need.

How it compares with fit()

Relative RMS error against the exact derivative on a 20,000-sample random walk, α = 0.5, h = 0.01. The fitted filters used eight 3,000-sample training walks.

Two things stand out. At small sizes the analytic filter fails badly and fit() is far better. At L=32, p=16 and above, with j_max set to the signal length, the analytic filter is excellent. The default j_max (3,200 for L=32) is too short for a 20,000-sample signal, and the error shows it.

Straight to C

from fracmem.embedded import export_c

filt = CompressedFractionalFilter.from_analytic(
    alpha=0.5, h=0.01, L=16, p=8,
    tail_fit_points=800,   # more careful than the default 400
    tail_fit_reg=0.05,     # trust the fit a bit more than the default 0.1
)
export_c(filt, "device_filter.c", instance_name="deployed")

This is examples/embedded-export-c.py.

Choosing L, p and j_max

Three numbers set accuracy, CPU cost and RAM. All figures on this page were measured with fracmem 0.2.0.

Cost per sample

The C struct reserves room for the maximum size (L=64, p=32) whatever you use. sizeof(FracmemFilter) is 924 bytes on x86-64 with gcc. Lower fracmemMaxL and fracmemMaxP in the header to shrink it.

Kernel accuracy

Worst-case relative error of the tail weights against the exact Grünwald–Letnikov weights, from soe_tail_error with j_max = 10,000:

Doubling both L and p cuts the kernel error by roughly 10 to 10,000 times. Lower orders (small α) have longer memory and are harder.

from fracmem import gl_weights, soe_tail_kernel, soe_tail_error

alpha, L, p, j_max = 0.5, 32, 16, 10_000
w = gl_weights(alpha, j_max + 1)
lam, c = soe_tail_kernel(alpha, L, p, w, j_max)
print(soe_tail_error(alpha, L, w, lam, c, j_max))   # 0.00674...

Rules of thumb

• Start at L=32, p=16. It is the default, it costs 48 multiply-adds, and it fits well on drifting signals.
• Keep p about L/2. All measured configurations use this ratio.
• Set j_max to the memory you need. For 100 Hz and 100 s of relevant history, use j_max=10_000.
• Stay within L ≤ 64 and p ≤ 32 for the C export, or raise the limits in fracmemfilter.h.
• Always check end-to-end. Run predict() on a realistic signal and compare with full_gl_derivative.

The sample period h

h only scales the output by h^-α. The filter shape depends on α, L, p and j_max. You still need a filter per sample rate, because j_max is counted in samples and the scale factor is baked into the C export.

Saving and loading

Save a fitted filter to a small .npz file, then load it in another process without refitting.

filt.save("half_derivative.npz")

from fracmem import CompressedFractionalFilter
filt = CompressedFractionalFilter.load("half_derivative.npz")
y = filt.predict(signal)

What the file holds

The local weights w are not stored; they are recomputed from α. reg_used and streaming state are not saved. A loaded filter starts with a fresh stream.

Deploy to embedded C

Fitting needs NumPy and SciPy. Running the filter needs about 60 lines of plain C, no heap and no libraries.

1. Export the constants
   python

   from fracmem.embedded import export_c
   
   export_c(filt, "device_filter.c")              # instance_name="filt" by default

2. Copy three files into your project

   File	Where from	How often
   fracmemfilter.h	src/fracmem/embedded/c/	Once
   fracmemfilter.c	src/fracmem/embedded/c/	Once
   device_filter.c	Generated by export_c	Each time you refit

3. Call it from your firmware
   c

   #include "fracmemfilter.h"
   
   extern FracmemFilter filt;
   void filtSetup(void);
   
   void app_init(void)
   {
       filtSetup();                       /* copies constants, clears state */
   }
   
   void control_tick(void)                /* called every h seconds */
   {
       float x = read_sensor();
       float d_half = fracmemStep(&filt, x);
       /* use d_half in your controller */
   }

4. Build
   bash

   gcc -O2 -o app main.c device_filter.c fracmemfilter.c -I path/to/fracmem/c

   With a firmware toolchain, add the two .c files to your sources and the header folder to your include path.

What export_c generates

Real output for from_analytic(alpha=0.5, h=0.01, L=16, p=8), with arrays shortened here:

/* Auto-generated by fracmem.embedded.export_c -- do not edit by hand.
   alpha=0.5  h=0.01  L=16  p=8  definition='gl'
   Compile together with fracmemfilter.c (see embedded/c/). */
#include "fracmemfilter.h"

static const float filtW[16] = {1.00000000f, -0.50000000f, -0.12500000f, ...};
static const float filtLam[8] = {...};
static const float filtC[8] = {...};

FracmemFilter filt;

void filtSetup(void)
{
    fracmemInit(&filt, 16, 8, 10.00000000f, 0,
                filtW, filtLam, filtC);
}

10.0 is h^-α = 0.01^-0.5. The 0 after it means Caputo is off.

Several filters in one firmware

Give each export its own instance_name. Compile fracmemfilter.c once.

export_c(roll_filter,  "roll_filter.c",  instance_name="rollD")
export_c(pitch_filter, "pitch_filter.c", instance_name="pitchD")

rollDSetup();
pitchDSetup();
float dr = fracmemStep(&rollD, roll);
float dp = fracmemStep(&pitchD, pitch);

Limits and memory

• fracmemMaxL = 64 and fracmemMaxP = 32 are set in the header. export_c raises ValueError if your filter is larger.
• Each FracmemFilter uses 924 bytes at the default limits. Lower the limits to save RAM.
• The runtime uses float (32-bit). On cores with a single-precision FPU, such as Cortex-M4F and ESP32, it runs in hardware.
• The source is C89: no VLAs, no malloc, no <math.h>.

Platform notes

extern "C" {
#include "fracmemfilter.h"
extern FracmemFilter filt;
void filtSetup(void);
}

void setup() { filtSetup(); }
void loop()  { float y = fracmemStep(&filt, analogRead(A0) * 0.001f); delay(10); }

idf_component_register(SRCS "main.c" "device_filter.c" "fracmemfilter.c"
                       INCLUDE_DIRS ".")

Verifying the C build

Before flashing, compile the exported file on your desktop and check that it gives the same numbers as Python.

Run the round-trip example

python examples/embedded-export-c.py

built with no training data: L=16 p=8 tail_fit_points=800 tail_fit_reg=0.05
exported 1044 bytes to device_filter.c

compiled and ran the exported C file successfully
max |desktop - compiled C| over 200 samples: 3.55e-07  (float32 rounding only)

The script exports a filter, writes a small main.c that reads samples from stdin, compiles it with cc or gcc, and compares the output with predict().

Check your own filter

#include <stdio.h>
#include "fracmemfilter.h"

extern FracmemFilter filt;
void filtSetup(void);

int main(void)
{
    float xk;
    filtSetup();
    while (scanf("%f", &xk) == 1) {
        printf("%.8f\n", fracmemStep(&filt, xk));
    }
    return 0;
}

import subprocess, numpy as np

subprocess.run(["gcc", "-O2", "-o", "check", "main.c", "device_filter.c",
                "fracmemfilter.c", "-I", "path/to/fracmem/c"], check=True)

x = np.cumsum(np.random.default_rng(7).standard_normal(200)) * 0.01
out = subprocess.run(["./check"], input="\n".join(map(repr, x)),
                     capture_output=True, text=True, check=True).stdout
y_c = np.array(out.split(), dtype=float)
print(np.max(np.abs(y_c - filt.predict(x))))

How it works

An exact fractional derivative is a convolution with a slowly decaying kernel. fracmem keeps the head of that kernel exact and replaces the tail with a short sum of exponentials.

The target

The Grünwald–Letnikov derivative of order α with sample period h:

The sum grows with k. The weights decay as a power law, so no fixed window is long enough to truncate it safely.

Split the history

Sum of exponentials, from an identity

The tail rests on an exact integral for a power law, with β = α + 1:

Discretize the integral at p log-spaced nodes s_i between 1/j_max and 8. Each node becomes one exponential mode with decay rate λ_i = e^−s_i:

An exponential needs no history. Its running sum obeys a one-line recursion:

Data only sets the readout

The decay rates λ are fixed by the quadrature and never fitted. The weights c then have two possible sources:

• Analytic (from_analytic): least squares against the exact kernel weights at log-spaced lags, with relative weighting and a small ridge term.
• Data (fit): ridge regression against the exact derivative of your training signals, with the penalty chosen by cross-validation.

With λ fixed, both are linear problems with a single solution. Refitting λ from data could only add noise to numbers that already come from an exact identity, so the library never does it.

Background

This builds on the sum-of-exponentials method of Jiang, Zhang, Zhang and Zhang for the Caputo derivative, and related work by Lubich and Schädle and by Baffet and Hesthaven on fast convolution for fractional operators. For the GL/RL relationship, see Podlubny, Fractional Differential Equations, chapters 2 and 3.

Accuracy and limits

What is guaranteed, what is measured, and where the filter can go wrong.

Exact parts

• The local weights w_0 … w_{L-1} are exact, from the stable recursion.
• The first L outputs of a fresh filter are exact, because the tail has not received any input yet.
• step() and predict() agree, and the C runtime agrees with both up to float32 rounding.

Kernel bound vs end-to-end error

soe_tail_error reports the worst relative mismatch between the modes and the exact tail weights. That is a statement about the kernel. End-to-end error on a signal also depends on the signal: a long random walk adds up small kernel errors across thousands of samples.

Measured numbers

Known limits

• Fixed sample rate. h is built into the filter. Jittery or variable-rate sampling is not modelled.
• Orders tested. The test suite runs α = 0.5. Other orders in (0, 1) work in the measurements above; check your own α.
• Memory horizon. History much older than j_max is under-weighted. Pick j_max for your use.
• float32 on device. Very long runs with large, slowly growing inputs can lose precision in the slowest modes. Feed the filter signals near zero mean when you can, or use Caputo.
• Python streaming speed. step() is a reference implementation in Python, fine at hundreds of hertz. Use the C runtime for fast loops.

CompressedFractionalFilter

The main class: build, fit, run, stream, save and load. from fracmem import CompressedFractionalFilter

Constructor

alpha float

Order of the derivative, typically 0 < α < 1.

h float

Sample period in seconds.

L int = 32

Length of the exact local window, in samples.

p int = 16

Number of exponential tail modes.

definition str = "gl"

"gl", "rl" or "caputo". Anything else raises ValueError.

fit method

train_signals list of 1-D arrays

Same-length example signals.

j_max int | None

Memory horizon in samples. Defaults to the training signal length.

reg_candidates tuple of float

Ridge strengths to cross-validate. Default 1e-11 … 1e-5.

Sets lam, c and reg_used. Returns the filter, so you can chain .fit(...).predict(...).

from_analytic classmethod

Builds a ready filter with no training data. j_max defaults to max(2000, 100*L). See Filters without training data.

predict method

Runs the filter over a whole 1-D signal and returns an array of the same length. w overrides the local weights; leave it as None. Raises RuntimeError before fitting.

reset_stream, step method

reset_stream clears the window, the modes and the stored Caputo start value. step feeds one sample and returns one estimate. Both raise RuntimeError before fitting.

save, load

Writes or reads an .npz file. See Saving and loading.

Attributes

Kernel and SOE functions

Lower-level building blocks, all importable from fracmem. Use them for reference computations, checks and experiments.

Kernel (fracmem.kernel)

Exact GL weights w_0 … w_{n-1} from the stable recursion.

The exact derivative by full convolution. Needs len(w) >= len(x). Cost grows with signal length; use it as a reference only.

Named versions of the exact derivative. rl_derivative equals full_gl_derivative. caputo_derivative subtracts x[0] first, and also accepts a 2-D batch.

The exact length-L FIR part, and a zero-padded shift by L samples (1-D or batch).

Sum of exponentials (fracmem.soe)

Raw quadrature nodes and weights, before any refit.

Least-squares refit of c against the exact tail weights.

Quadrature plus refit: the full classical construction with no data.

Worst relative error of c · lam^m against w[L+m], probed at log-spaced lags.

Filter helpers (fracmem.filter)

The deployed tail recursion, written as a plain loop.

The pieces fit() uses: batched mode traces, energy-scaled ridge solve, and k-fold selection of the ridge strength.

fracmem.embedded

Python-side export of a filter's constants to C.

export_c

Writes the output of generate_source_c to path.

generate_source_c

filt CompressedFractionalFilter

A fitted or analytic filter. Raises RuntimeError if neither.

w array | None

Local weights to embed. Defaults to gl_weights(alpha, L).

instance_name str = "filt"

Prefix for every generated symbol: <name>W, <name>Lam, <name>C, the FracmemFilter <name> object, and <name>Setup().

Raises ValueError if L > 64 or p > 32. Constants are printed with 8 decimals as float literals.

C runtime API

Everything in fracmemfilter.h. Three functions and one struct.

Limits

Compile-time maxima for every filter. Change them and rebuild to trade RAM for size.

FracmemFilter

fracmemInit

Copies the constants into f and resets its state. The generated <name>Setup() calls this for you.

fracmemReset

Zeros the window and modes and forgets the Caputo start value. Constants are kept.

fracmemStep

Feeds one sample and returns one derivative estimate. Exactly L + p multiply-adds. Not reentrant on the same f; use one filter per signal and per thread.

Examples

Both scripts live in examples/ and run in CI on every push.

git clone https://github.com/edith-lang/fracmem.git
cd fracmem
pip install -e .
python examples/basic-usage.py
python examples/embedded-export-c.py

Their output is shown on the Quickstart and Verifying the C build pages.

Testing and contributing

fracmem is small and checked by tests. Changes should keep it that way.

The test suite

CI runs the tests and both examples on Python 3.9, 3.10, 3.11 and 3.12 for every push and pull request to main.

Design rules

• Never fit λ from data. The decay rates come from an exact identity. If you believe a case needs it, open an issue with data first.
• No heavy dependencies. NumPy and SciPy on the desktop, plain C on the device.
• Every fix comes with a test that fails without it.

Pull requests

1. Fork, branch, and make the change with tests.
2. Run pytest tests/ -v.
3. Explain why in the PR, especially for anything in the fitting logic.

Reports from real hardware are especially welcome. Include the board, α, h, L, p and what you measured.

FAQ and troubleshooting

Common errors and questions.

RuntimeError: call .fit(...) before .predict(...)

The filter has no lam and c yet. Call fit(), build it with from_analytic(), or load() a saved one. The message mentions fit in every case.

fit() fails with a shape error

Your training signals have different lengths. Trim them to a common length.

ValueError: L=… p=… exceed fracmemMaxL=64 / fracmemMaxP=32

Use a smaller filter, or raise the limits in fracmemfilter.h and use that header in your firmware.

pytest crashes on import of an unrelated plugin

Other software on your machine, such as ROS, can put pytest plugins on PYTHONPATH. Run the tests in a clean environment:

env -u PYTHONPATH PYTEST_DISABLE_PLUGIN_AUTOLOAD=1 python -m pytest tests/ -q

The output drifts on a long run

• Check that j_max covers the history you care about.
• If the signal starts away from zero, use definition="caputo".
• If you used from_analytic at a small size, fit on real data instead.

Can α be larger than 1?

The method and tests target 0 < α < 1. For orders above 1, take integer derivatives separately and apply fracmem to the fractional part.

Can it compute fractional integrals?

Not in 0.2.0. The kernel and the C runtime are written for derivatives.

License

MIT. Author: Adithyan Lalu. Source and issues on GitHub.