<commit_before>class InverseLaplaceTransform(object):
    r"""
    Inverse Laplace transform methods are implemented using this
    class, in order to simplify the code and provide a common
    infrastructure.

    Implement a custom inverse Laplace transform algorithm by
    subclassing :class:`InverseLaplaceTransform` and implementing the
    appropriate methods. The subclass can then be used by
    :func:`~mpmath.invertlaplace` by passing it as the *method*
    argument.
    """

    def __init__(self,ctx):
        self.ctx = ctx

    def calc_laplace_parameter(self,t,**kwargs):
        """
        Determine the vector of Laplace parameter values needed for an
        algorithm, this will depend on the choice of algorithm (de
        Hoog is default), the algorithm-specific parameters passed (or
        default ones), and desired time.
        """
        raise NotImplementedError

    def calc_time_domain_solution(self,fp):
        """
        Compute the time domain solution, after computing the
        Laplace-space function evaluations at the abscissa required
        for the algorithm. Abscissa computed for one algorithm are
        typically not useful for another algorithm.
        """
        raise NotImplementedError

class FixedTalbot(InverseLaplaceTransform):

    def calc_laplace_parameter(self,t,**kwargs):
        r"""
        The "fixed" Talbot method deforms the Bromwich contour towards
        `-\infty` in the shape of a parabola. Traditionally the Talbot
        algorithm has adjustable parameters, but the "fixed" version
        does not. The `r` parameter could be passed in as a parameter,
        if you want to override the default given by (Abate & Valko,
        2004).

        The Laplace parameter is sampled along a parabola opening
        along the negative imaginary axis, with the bottom of the
        parabola at `p=\frac{r}{t_\mathrm{max}}+0j`

        **Optional arguments**

        :class:`~mpmath.calculus.inverselaplace.FixedTalbot.calc_laplace_parameter`
        recognizes the following keywords

        *tmax*
            maximum time associated with vector of times
            (typically just the time requested)
        *degree*
            integer order of approximation (M = number of terms)
        *r*
            abscissa for `p_0` (otherwise computed using rule
            of thumb `2M/5`)

        The working precision will be increased according to a rule of
        thumb, in relation to the degree; dps_goal`=1.7M`.  dps_goal
        can be overridden with a user-specified value.

        `p_0=\frac{r}{t}`

        `p_i=\frac{i r \pi}{Mt_\mathrm{max}}\left[\cot\left(
        \frac{i\pi}{M}\right) + j \right] \qquad 1\le i <M`

        where `j=\sqrt{-1}`, `r=2M/5`, and `t_\mathrm{max}` is the
        maximum specified time.
        """

        # required
        # ------------------------------
        # time of desired approximation
        self.t = self.ctx.convert(t)

        # optional
        # ------------------------------
        # maximum time desired (used for scaling)
        # default is requested time.
        self.tmax = self.ctx.convert(kwargs.get('tmax',self.t))

        # empirical relationships used here based on a linear fit of
        # requested and delivered dps for exponentially decaying time
        # functions for requested dps up to 512.

        if 'degree' in kwargs:
            self.degree = kwargs['degree']
            self.dps_goal = self.degree
        else:
            self.dps_goal = int(1.72*self.ctx.dps)
            self.degree = max(12,int(1.38*self.dps_goal))

        M = self.degree

        # this is adjusting the dps of the calling context hopefully
        # the caller doesn't monkey around with it between calling
        # this routine and calc_time_domain_solution()
        self.dps_orig = self.ctx.dps
        self.ctx.dps = self.dps_goal

        # Abate & Valko rule of thumb for r parameter
        self.r = kwargs.get('r',self.ctx.fraction(2,5)*M)

        self.theta = self.ctx.linspace(0.0, self.ctx.pi, M+1)

        self.cot_theta = self.ctx.matrix(M,1)
        self.cot_theta[0] = 0 # not used

        # all but time-dependent part of p
        self.delta = self.ctx.matrix(M,1)
        self.delta[0] = self.r

        for i in range(1,M):
            self.cot_theta[i] = self.ctx.cot(self.theta[i])
            self.delta[i] = self.r*self.theta[i]*(self.cot_theta[i] + 1j)

        self.p = self.ctx.matrix(M,1)
        self.p = self.delta/self.tmax

        # NB: p is complex (mpc)

    def calc_time_domain_solution(self,fp,t,manual_prec=False):
        r"""
        The fixed Talbot time-domain solution is computed from the Laplace-space
        function evaluations using

        `f(t,M)=\frac{2}{5t}\sum_{k=0}^{M-1}\mathrm{Re}
        \left[ \gamma_k \bar{f}(p_k)\right]`

        where

        `\gamma_0 = \frac{1}{2}e^{r}\bar{f}(p_0)`

        `\gamma_k = e^{tp_k}\left\lbrace 1 + \frac{jk\pi}{M}\left[1 +
        \cot \left( \frac{k \pi}{M} \right)^2 \right] -
        j\cot\left( \frac{k \pi}{M}\right)\right \rbrace \qquad 1\le k<M`.

        Again, `j=\sqrt{-1}`.

        Before calling this function, call
        :class:`~mpmath.calculus.inverselaplace.FixedTalbot.calc_laplace_parameter`
        to set the parameters and compute the required coefficients.
        """

        # required
        # ------------------------------
        self.t = self.ctx.convert(t)

        # assume fp was computed from p matrix returned from
        # calc_laplace_parameter(), so is already
        # a list or matrix of mpmath 'mpc' types

        # these were computed in previous call to calc_laplace_parameter()
        theta = self.theta
        delta = self.delta
        M = self.degree
        p = self.p
        r = self.r
            
        ans = self.ctx.matrix(M,1)
        ans[0] = self.ctx.exp(delta[0])*fp[0]/2

        for i in range(1,M):
            ans[i] = self.ctx.exp(delta[i])*fp[i]*(
                1 + 1j*theta[i]*(1 + self.cot_theta[i]**2) -
                1j*self.cot_theta[i])

        result = self.ctx.fraction(2,5)*self.ctx.fsum(ans)/self.t

        # setting dps back to value when calc_laplace_parameter was called
        if not manual_prec:
            self.ctx.dps = self.dps_orig

        return result.real

# ****************************************

class Stehfest(InverseLaplaceTransform):

    def calc_laplace_parameter(self,t,**kwargs):
        r"""
        The Gaver-Stehfest method is a discrete approximation of the
        Widder-Post inversion algorithm, rather than a direct
        approximation of the Bromwich contour integral.

        The method only uses values of `p` along the real axis, and
        therefore has issues inverting oscillatory functions (which
        have poles in pairs away from the real axis).

        The Stehfest coefficients only depend on the order of
        approximation and not the desired time, so they could be
        computed once and re-used for a given degree and working
        precision.
        """

        # required
        # ------------------------------
        # time of desired approximation
        self.t = self.ctx.convert(t)

        # optional
        # ------------------------------

        # empirical relationships used here based on a linear fit of
        # requested and delivered dps for exponentially decaying time
        # functions for requested dps up to 512.

        if 'degree' in kwargs:
            self.degree = kwargs['degree']
            self.dps_goal = int(1.38*self.degree)
        else:
            self.dps_goal = int(2.93*self.ctx.dps)
            self.degree = max(16,self.dps_goal)

        # _coeff routine requires even degree
        if self.degree%2 > 0:
            self.degree += 1

        M = self.degree

        # this is adjusting the dps of the calling context
        # hopefully the caller doesn't monkey around with it
        # between calling this routine and calc_time_domain_solution()
        self.dps_orig = self.ctx.dps
        self.ctx.dps = self.dps_goal

        self.V = self._coeff()
        self.p = self.ctx.matrix(self.ctx.arange(1,M+1))*self.ctx.ln2/self.t

        # NB: p is real (mpf)

    def _coeff(self):
        """These Salzer summation weights (aka, "Stehfest coefficients")
        only depend on the approximation order (M) and the precision"""

        M = self.degree
        M2 = M/2

        V = self.ctx.matrix(M,1)

        # Salzer summation weights
        # get very large in magnitude and oscillate in sign,
        # if the precision is not high enough, there will be
        # catastrophic cancellation
        for k in range(1,M+1):
            z = self.ctx.matrix(int(min(k,M2)+1),1)
            for j in range(int(self.ctx.floor((k+1)/2.0)),min(k,M2)+1):
                z[j] = (self.ctx.power(j,M2)*self.ctx.fac(2*j)/
                        (self.ctx.fac(M2-j)*self.ctx.fac(j)*
                         self.ctx.fac(j-1)*self.ctx.fac(k-j)*
                         self.ctx.fac(2*j-k)))
            V[k-1] = self.ctx.power(-1,k+M2)*self.ctx.fsum(z)

        return V

    def calc_time_domain_solution(self,fp,t,manual_prec=False):
        """Compute time-domain Stehfest algorithm solution using f(p)
        and coefficients"""

        # required
        self.t = self.ctx.convert(t)

        # assume fp was computed from p matrix returned from
        # calc_laplace_parameter(), so is already
        # a list or matrix of mpmath 'mpf' types

        result = self.ctx.fdot(self.V,fp)*self.ctx.ln2/self.t

        # setting dps back to value when calc_laplace_parameter was called
        if not manual_prec:
            self.ctx.dps = self.dps_orig

        # ignore any small imaginary part
        return result.real

# ****************************************

class deHoog(InverseLaplaceTransform):

    def calc_laplace_parameter(self,t,**kwargs):
        r"""the de Hoog, Knight & Stokes algorithm is an
        accelerated form of the Fourier series numerical
        inverse Laplace transform algorithms, like Crump (YEAR)
        or Dubner & Abate (YEAR).
        """

        self.t = self.ctx.convert(t)
        self.tmax = self.t

        # empirical relationships used here based on a linear fit of
        # requested and delivered dps for exponentially decaying time
        # functions for requested dps up to 512.

        if 'degree' in kwargs:
            self.degree = kwargs['degree']
            self.dps_goal = int(1.38*self.degree)
        else:
            self.dps_goal = int(self.ctx.dps*1.36)
            self.degree = max(10,self.dps_goal)

        # 2*M+1 terms in approximation
        M = self.degree

        # adjust alpha component of abscissa of convergence for higher precision
        tmp = self.ctx.power(10.0,-self.dps_goal)
        self.alpha = self.ctx.convert(kwargs.get('alpha',tmp))

        # desired tolerance (here simply related to alpha)
        self.tol = self.ctx.convert(kwargs.get('tol',self.alpha*10.0))
        self.np = 2*self.degree+1 # number of terms in approximation

        # this is adjusting the dps of the calling context
        # hopefully the caller doesn't monkey around with it
        # between calling this routine and calc_time_domain_solution()
        self.dps_orig = self.ctx.dps
        self.ctx.dps = self.dps_goal

        # scaling factor (likely tun-able, but 2 is typical)
        self.scale = kwargs.get('scale',2)
        self.T = self.ctx.convert(kwargs.get('T',self.scale*self.tmax))

        self.p = self.ctx.matrix(2*M+1,1)
        self.gamma = self.alpha - self.ctx.log(self.tol)/(self.scale*self.T)
        self.p = (self.gamma + self.ctx.pi*
                  self.ctx.matrix(self.ctx.arange(self.np))/self.T*1j)

        # NB: p is complex (mpc)

    def calc_time_domain_solution(self,fp,t,manual_prec=False):
        r"""Calculate time-domain solution for
        de Hoog, Knight & Stokes algorithm.
        """
        
        M = self.degree
        np = self.np
        T = self.T

        self.t = self.ctx.convert(t)

        # would it be useful to try re-using
        # space between e&q and A&B?
        e = self.ctx.zeros(np,M+1)
        q = self.ctx.matrix(np,M)
        d = self.ctx.matrix(np,1)
        A = self.ctx.zeros(np+2,1)
        B = self.ctx.ones(np+2,1)

        # initialize Q-D table
        # e[0:2*M,0] = 0.0 + 0.0j
        q[0,0] = fp[1]/(fp[0]/2)
        for i in range(1,2*M):
            q[i,0] = fp[i+1]/fp[i]

        # rhombus rule for filling triangular Q-D table (e & q)
        for r in range(1,M+1):
            # start with e, column 1, 0:2*M-2
            mr = 2*(M-r)<commit_msg>adding some formulae to documentation of individual methods<commit_after>class InverseLaplaceTransform(object):
    r"""
    Inverse Laplace transform methods are implemented using this
    class, in order to simplify the code and provide a common
    infrastructure.

    Implement a custom inverse Laplace transform algorithm by
    subclassing :class:`InverseLaplaceTransform` and implementing the
    appropriate methods. The subclass can then be used by
    :func:`~mpmath.invertlaplace` by passing it as the *method*
    argument.
    """

    def __init__(self,ctx):
        self.ctx = ctx

    def calc_laplace_parameter(self,t,**kwargs):
        """
        Determine the vector of Laplace parameter values needed for an
        algorithm, this will depend on the choice of algorithm (de
        Hoog is default), the algorithm-specific parameters passed (or
        default ones), and desired time.
        """
        raise NotImplementedError

    def calc_time_domain_solution(self,fp):
        """
        Compute the time domain solution, after computing the
        Laplace-space function evaluations at the abscissa required
        for the algorithm. Abscissa computed for one algorithm are
        typically not useful for another algorithm.
        """
        raise NotImplementedError

class FixedTalbot(InverseLaplaceTransform):

    def calc_laplace_parameter(self,t,**kwargs):
        r"""
        The "fixed" Talbot method deforms the Bromwich contour towards
        `-\infty` in the shape of a parabola. Traditionally the Talbot
        algorithm has adjustable parameters, but the "fixed" version
        does not. The `r` parameter could be passed in as a parameter,
        if you want to override the default given by (Abate & Valko,
        2004).

        The Laplace parameter is sampled along a parabola opening
        along the negative imaginary axis, with the base of the
        parabola along the real axis at
        `p=\frac{r}{t_\mathrm{max}}`. As the number of terms used in
        the approximation (degree) grows, the abscissa required for
        function evaluation tend towards `-\infty`, requiring high
        precision to prevent overflow.

        **Optional arguments**

        :class:`~mpmath.calculus.inverselaplace.FixedTalbot.calc_laplace_parameter`
        recognizes the following keywords

        *tmax*
            maximum time associated with vector of times
            (typically just the time requested)
        *degree*
            integer order of approximation (M = number of terms)
        *r*
            abscissa for `p_0` (otherwise computed using rule
            of thumb `2M/5`)

        The working precision will be increased according to a rule of
        thumb. If 'degree' is not specified, the working precision and
        degree are chosen to hopefully achieve the dps of the calling
        context. If 'degree' is specified, the working precision is
        chosen to achieve maximum resulting precision for the
        specified degree.

        `p_0=\frac{r}{t}`

        `p_i=\frac{i r \pi}{Mt_\mathrm{max}}\left[\cot\left(
        \frac{i\pi}{M}\right) + j \right] \qquad 1\le i <M`

        where `j=\sqrt{-1}`, `r=2M/5`, and `t_\mathrm{max}` is the
        maximum specified time.
        """

        # required
        # ------------------------------
        # time of desired approximation
        self.t = self.ctx.convert(t)

        # optional
        # ------------------------------
        # maximum time desired (used for scaling) default is requested
        # time.
        self.tmax = self.ctx.convert(kwargs.get('tmax',self.t))

        # empirical relationships used here based on a linear fit of
        # requested and delivered dps for exponentially decaying time
        # functions for requested dps up to 512.

        if 'degree' in kwargs:
            self.degree = kwargs['degree']
            self.dps_goal = self.degree
        else:
            self.dps_goal = int(1.72*self.ctx.dps)
            self.degree = max(12,int(1.38*self.dps_goal))

        M = self.degree

        # this is adjusting the dps of the calling context hopefully
        # the caller doesn't monkey around with it between calling
        # this routine and calc_time_domain_solution()
        self.dps_orig = self.ctx.dps
        self.ctx.dps = self.dps_goal

        # Abate & Valko rule of thumb for r parameter
        self.r = kwargs.get('r',self.ctx.fraction(2,5)*M)

        self.theta = self.ctx.linspace(0.0, self.ctx.pi, M+1)

        self.cot_theta = self.ctx.matrix(M,1)
        self.cot_theta[0] = 0 # not used

        # all but time-dependent part of p
        self.delta = self.ctx.matrix(M,1)
        self.delta[0] = self.r

        for i in range(1,M):
            self.cot_theta[i] = self.ctx.cot(self.theta[i])
            self.delta[i] = self.r*self.theta[i]*(self.cot_theta[i] + 1j)

        self.p = self.ctx.matrix(M,1)
        self.p = self.delta/self.tmax

        # NB: p is complex (mpc)

    def calc_time_domain_solution(self,fp,t,manual_prec=False):
        r"""
        The fixed Talbot time-domain solution is computed from the
        Laplace-space function evaluations using

        `f(t,M)=\frac{2}{5t}\sum_{k=0}^{M-1}\mathrm{Re} \left[
        \gamma_k \bar{f}(p_k)\right]`

        where

        `\gamma_0 = \frac{1}{2}e^{r}\bar{f}(p_0)`

        `\gamma_k = e^{tp_k}\left\lbrace 1 + \frac{jk\pi}{M}\left[1 +
        \cot \left( \frac{k \pi}{M} \right)^2 \right] - j\cot\left(
        \frac{k \pi}{M}\right)\right \rbrace \qquad 1\le k<M`.

        Again, `j=\sqrt{-1}`.

        Before calling this function, call
        :class:`~mpmath.calculus.inverselaplace.FixedTalbot.calc_laplace_parameter`
        to set the parameters and compute the required coefficients.
        """

        # required
        # ------------------------------
        self.t = self.ctx.convert(t)

        # assume fp was computed from p matrix returned from
        # calc_laplace_parameter(), so is already a list or matrix of
        # mpmath 'mpc' types

        # these were computed in previous call to
        # calc_laplace_parameter()
        theta = self.theta
        delta = self.delta
        M = self.degree
        p = self.p
        r = self.r
            
        ans = self.ctx.matrix(M,1)
        ans[0] = self.ctx.exp(delta[0])*fp[0]/2

        for i in range(1,M):
            ans[i] = self.ctx.exp(delta[i])*fp[i]*(
                1 + 1j*theta[i]*(1 + self.cot_theta[i]**2) -
                1j*self.cot_theta[i])

        result = self.ctx.fraction(2,5)*self.ctx.fsum(ans)/self.t

        # setting dps back to value when calc_laplace_parameter was
        # called, unless flag is set.
        if not manual_prec:
            self.ctx.dps = self.dps_orig

        return result.real

# ****************************************

class Stehfest(InverseLaplaceTransform):

    def calc_laplace_parameter(self,t,**kwargs):
        r"""
        The Gaver-Stehfest method is a discrete approximation of the
        Widder-Post inversion algorithm, rather than a direct
        approximation of the Bromwich contour integral.

        The method abscissa along the real axis, and therefore has
        issues inverting oscillatory functions (which have poles in
        pairs away from the real axis).

        The working precision will be increased according to a rule of
        thumb. If 'degree' is not specified, the working precision and
        degree are chosen to hopefully achieve the dps of the calling
        context. If 'degree' is specified, the working precision is
        chosen to achieve maximum resulting precision for the
        specified degree.

        `p_j = \frac{j \log 2}{t} \qquad 1 \le j \le M`
        """

        # required
        # ------------------------------
        # time of desired approximation
        self.t = self.ctx.convert(t)

        # optional
        # ------------------------------

        # empirical relationships used here based on a linear fit of
        # requested and delivered dps for exponentially decaying time
        # functions for requested dps up to 512.

        if 'degree' in kwargs:
            self.degree = kwargs['degree']
            self.dps_goal = int(1.38*self.degree)
        else:
            self.dps_goal = int(2.93*self.ctx.dps)
            self.degree = max(16,self.dps_goal)

        # _coeff routine requires even degree
        if self.degree%2 > 0:
            self.degree += 1

        M = self.degree

        # this is adjusting the dps of the calling context
        # hopefully the caller doesn't monkey around with it
        # between calling this routine and calc_time_domain_solution()
        self.dps_orig = self.ctx.dps
        self.ctx.dps = self.dps_goal

        self.V = self._coeff()
        self.p = self.ctx.matrix(self.ctx.arange(1,M+1))*self.ctx.ln2/self.t

        # NB: p is real (mpf)

    def _coeff(self):
        r"""Salzer summation weights (aka, "Stehfest coefficients")
        only depend on the approximation order (M) and the precision"""

        M = self.degree
        M2 = M/2

        V = self.ctx.matrix(M,1)

        # Salzer summation weights
        # get very large in magnitude and oscillate in sign,
        # if the precision is not high enough, there will be
        # catastrophic cancellation
        for k in range(1,M+1):
            z = self.ctx.matrix(int(min(k,M2)+1),1)
            for j in range(int(self.ctx.floor((k+1)/2.0)),min(k,M2)+1):
                z[j] = (self.ctx.power(j,M2)*self.ctx.fac(2*j)/
                        (self.ctx.fac(M2-j)*self.ctx.fac(j)*
                         self.ctx.fac(j-1)*self.ctx.fac(k-j)*
                         self.ctx.fac(2*j-k)))
            V[k-1] = self.ctx.power(-1,k+M2)*self.ctx.fsum(z)

        return V

    def calc_time_domain_solution(self,fp,t,manual_prec=False):
        r"""
        Compute time-domain Stehfest algorithm solution.

        `f(t,M) = \frac{\log 2}{t} \sum_{k=1}^{M} V_k \bar{f}\left(
        p_k \right)`

        where

        `V_k=(-1)^{k + N/2}\sum^{\min(k,N/2)}_{j=\lfloor(k+1)/2
        \rfloor}\frac{j^{\frac{N}{2}}(2j)!}{(\frac{N}{2}-j)!\,j!\,
        (j-1)!\,(k-j)!\,(2j-k)!}.
        """

        # required
        self.t = self.ctx.convert(t)

        # assume fp was computed from p matrix returned from
        # calc_laplace_parameter(), so is already
        # a list or matrix of mpmath 'mpf' types

        result = self.ctx.fdot(self.V,fp)*self.ctx.ln2/self.t

        # setting dps back to value when calc_laplace_parameter was called
        if not manual_prec:
            self.ctx.dps = self.dps_orig

        # ignore any small imaginary part
        return result.real

# ****************************************

class deHoog(InverseLaplaceTransform):

    def calc_laplace_parameter(self,t,**kwargs):
        r"""the de Hoog, Knight & Stokes algorithm is an
        accelerated form of the Fourier series numerical
        inverse Laplace transform algorithms, like Crump (YEAR)
        or Dubner & Abate (YEAR).

        `p = `
        """

        self.t = self.ctx.convert(t)
        self.tmax = self.t

        # empirical relationships used here based on a linear fit of
        # requested and delivered dps for exponentially decaying time
        # functions for requested dps up to 512.

        if 'degree' in kwargs:
            self.degree = kwargs['degree']
            self.dps_goal = int(1.38*self.degree)
        else:
            self.dps_goal = int(self.ctx.dps*1.36)
            self.degree = max(10,self.dps_goal)

        # 2*M+1 terms in approximation
        M = self.degree

        # adjust alpha component of abscissa of convergence for higher
        # precision
        tmp = self.ctx.power(10.0,-self.dps_goal)
        self.alpha = self.ctx.convert(kwargs.get('alpha',tmp))

        # desired tolerance (here simply related to alpha)
        self.tol = self.ctx.convert(kwargs.get('tol',self.alpha*10.0))
        self.np = 2*self.degree+1 # number of terms in approximation

        # this is adjusting the dps of the calling context
        # hopefully the caller doesn't monkey around with it
        # between calling this routine and calc_time_domain_solution()
        self.dps_orig = self.ctx.dps
        self.ctx.dps = self.dps_goal

        # scaling factor (likely tun-able, but 2 is typical)
        self.scale = kwargs.get('scale',2)
        self.T = self.ctx.convert(kwargs.get('T',self.scale*self.tmax))

        self.p = self.ctx.matrix(2*M+1,1)
        self.gamma = self.alpha - self.ctx.log(self.tol)/(self.scale*self.T)
        self.p = (self.gamma + self.ctx.pi*
                  self.ctx.matrix(self.ctx.arange(self.np))/self.T*1j)

        # NB: p is complex (mpc)

    def calc_time_domain_solution(self,fp,t,manual_prec=False):
        r"""Calculate time-domain solution for
        de Hoog, Knight & Stokes algorithm.

        The un-accelerated Fourier series approach:

        `f(t,2M+1) = \frac{e^{\gamma t}}{T} \sum_{k=0}^{2M}{}^{'}
        \Re\left[\bar{f}\left(\gamma_0 +\frac{i\pi k}{T}\right)
        \exp\left(\frac{i\pi t}{T}\right)\right],` 
        """
        
        M = self.degree
        np = self.np
        T = self.T

        self.t = self.ctx.convert(t)

        # would it be useful to try re-using
        # space between e&q and A&B?
        e = self.ctx.zeros(np,M+1)
        q = self.ctx.matrix(np,M)
        d = self.ctx.matrix(np,1)
        A = self.ctx.zeros(np+2,1)
        B = self.ctx.ones(np+2,1)

        # initialize Q-D table
        # e[0:2*M,0] = 0.0 + 0.0j
        q[0,0] = fp[1]/(fp[0]/2)
        for i in range(1,2*M):
            q[i,0] = fp[i+1]/fp[i]

        # rhombus rule for filling triangular Q-D table (e & q)
        for r in range(1,M+1):
            # start with e, column 1, 0:2*M-2
            mr = 2*(M-r)
{'fraction_non_alphanumeric': 0.12131689717378231, 'fraction_numerical': 0.012627781118460614, 'mean_word_length': 2.026501364877161, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '5304223', 'n_tokens_mistral': 8336, 'n_tokens_neox': 7843, 'n_words': 2318}