def ortho_basis_at_mp(self, p, q, r):
        r = r if r != 1 else r + mp.eps

        a = 2*p/(1 - r)
        b = 2*q/(1 - r)
        c = r

        sk = [mp.mpf(2)**(-k - 0.25)*mp.sqrt(k + 0.5)
              for k in xrange(self.order)]
        pa = [s*jp for s, jp in zip(sk, jacobi(self.order - 1, 0, 0, a))]
        pb = [s*jp for s, jp in zip(sk, jacobi(self.order - 1, 0, 0, b))]

        ob = []
        for i, pi in enumerate(pa):
            for j, pj in enumerate(pb):
                cij = (1 - c)**(i + j)
                pij = pi*pj

                pc = jacobi(self.order - max(i, j) - 1, 2*(i + j + 1), 0, c)
                for k, pk in enumerate(pc):
                    ck = mp.sqrt(2*(k + j + i) + 3)

                    ob.append(cij*ck*pij*pk)

        return ob

    def jac_ortho_basis_at_mp(self, p, q, r):
        a = 2 * p / (1 - r) if r != 1 else 0
        b = 2 * q / (1 - r) if r != 1 else 0
        c = r

        sk = [mp.mpf(2) ** (-k - 0.25) * mp.sqrt(k + 0.5) for k in range(self.order)]
        fc = [s * jp for s, jp in zip(sk, jacobi(self.order - 1, 0, 0, a))]
        gc = [s * jp for s, jp in zip(sk, jacobi(self.order - 1, 0, 0, b))]

        dfc = [s * jp for s, jp in zip(sk, jacobi_diff(self.order - 1, 0, 0, a))]
        dgc = [s * jp for s, jp in zip(sk, jacobi_diff(self.order - 1, 0, 0, b))]

        ob = []
        for i, (fi, dfi) in enumerate(zip(fc, dfc)):
            for j, (gj, dgj) in enumerate(zip(gc, dgc)):
                h = jacobi(self.order - max(i, j) - 1, 2 * (i + j + 1), 0, c)
                dh = jacobi_diff(self.order - max(i, j) - 1, 2 * (i + j + 1), 0, c)

                for k, (hk, dhk) in enumerate(zip(h, dh)):
                    ck = mp.sqrt(2 * (k + j + i) + 3)

                    tmp = (1 - c) ** (i + j - 1) if i + j > 0 else 1

                    pijk = 2 * tmp * dfi * gj * hk
                    qijk = 2 * tmp * fi * dgj * hk
                    rijk = (
                        tmp * (a * dfi * gj + b * fi * dgj - (i + j) * fi * gj) * hk
                        + (1 - c) ** (i + j) * fi * gj * dhk
                    )

                    ob.append([ck * pijk, ck * qijk, ck * rijk])

        return ob

def test_svd_test_case():
    # a test case from Golub and Reinsch
    #  (see wilkinson/reinsch: handbook for auto. comp., vol ii-linear algebra, 134-151(1971).)

    eps = mp.exp(0.8 * mp.log(mp.eps))

    a = [[22, 10,  2,   3,  7],
         [14,  7, 10,   0,  8],
         [-1, 13, -1, -11,  3],
         [-3, -2, 13,  -2,  4],
         [ 9,  8,  1,  -2,  4],
         [ 9,  1, -7,   5, -1],
         [ 2, -6,  6,   5,  1],
         [ 4,  5,  0,  -2,  2]]

    a = mp.matrix(a)
    b = mp.matrix([mp.sqrt(1248), 20, mp.sqrt(384), 0, 0])

    S = mp.svd_r(a, compute_uv = False)
    S -= b
    assert mp.mnorm(S) < eps

    S = mp.svd_c(a, compute_uv = False)
    S -= b
    assert mp.mnorm(S) < eps

def test_gauss_quadrature_dynamic(verbose = False):
    n = 5

    A = mp.randmatrix(2 * n, 1)

    def F(x):
        r = 0
        for i in xrange(len(A) - 1, -1, -1):
            r = r * x + A[i]
        return r

    def run(qtype, FW, R, alpha = 0, beta = 0):
        X, W = mp.gauss_quadrature(n, qtype, alpha = alpha, beta = beta)

        a = 0
        for i in xrange(len(X)):
            a += W[i] * F(X[i])

        b = mp.quad(lambda x: FW(x) * F(x), R)

        c = mp.fabs(a - b)

        if verbose:
            print(qtype, c, a, b)

        assert c < 1e-5

    run("legendre", lambda x: 1, [-1, 1])
    run("legendre01", lambda x: 1, [0, 1])
    run("hermite", lambda x: mp.exp(-x*x), [-mp.inf, mp.inf])
    run("laguerre", lambda x: mp.exp(-x), [0, mp.inf])
    run("glaguerre", lambda x: mp.sqrt(x)*mp.exp(-x), [0, mp.inf], alpha = 1 / mp.mpf(2))
    run("chebyshev1", lambda x: 1/mp.sqrt(1-x*x), [-1, 1])
    run("chebyshev2", lambda x: mp.sqrt(1-x*x), [-1, 1])
    run("jacobi", lambda x: (1-x)**(1/mp.mpf(3)) * (1+x)**(1/mp.mpf(5)), [-1, 1], alpha = 1 / mp.mpf(3), beta = 1 / mp.mpf(5) )

def _CalculateLeastUpperBoundInoperativeInterval(x0, x1, v0, v1, vm, am):
    # All input must already be of mp.mpf type
    d = x1 - x0
    temp1 = Prod([number('2'), Neg(Sqr(am)), Sub(Prod([number('2'), am, d]), Add(Sqr(v0), Sqr(v1)))])
    if temp1 < zero:
        T0 = number('-1')
        T1 = number('-1')
    else:
        term1 = mp.fdiv(Add(v0, v1), am)
        term2 = mp.fdiv(mp.sqrt(temp1), Sqr(am))
        T0 = Add(term1, term2)
        T1 = Sub(term1, term2)

    temp2 = Prod([number('2'), Sqr(am), Add(Prod([number('2'), am, d]), Add(Sqr(v0), Sqr(v1)))])
    if temp2 < zero:
        T2 = number('-1')
        T3 = number('-1')
    else:
        term1 = Neg(mp.fdiv(Add(v0, v1), am))
        term2 = mp.fdiv(mp.sqrt(temp2), Sqr(am))
        T2 = Add(term1, term2)
        T3 = Sub(term1, term2)

    newDuration = max(max(T0, T1), max(T2, T3))
    if newDuration > zero:
        dStraight = Prod([pointfive, Add(v0, v1), newDuration])
        if Sub(d, dStraight) > 0:
            amNew = am
            vmNew = vm
        else:
            amNew = -am
            vmNew = -vm

        # import IPython; IPython.embed()

        vp = Mul(pointfive, Sum([Mul(newDuration, amNew), v0, v1])) # the peak velocity
        if (Abs(vp) > vm):
            dExcess = mp.fdiv(Sqr(Sub(vp, vmNew)), am)
            assert(dExcess > 0)
            deltaTime = mp.fdiv(dExcess, vm)
            newDuration = Add(newDuration, deltaTime)

        log.debug('Calculation successful: T0 = {0}; T1 = {1}; T2 = {2}; T3 = {3}'.format(mp.nstr(T0, n=_prec), mp.nstr(T1, n=_prec), mp.nstr(T2, n=_prec), mp.nstr(T3, n=_prec)))
        
        newDuration = Mul(newDuration, number('1.01')) # add 1% safety bound
        return newDuration
    else:
        if (FuzzyEquals(x0, x1, epsilon) and FuzzyZero(v0, epsilon) and FuzzyZero(v1, epsilon)):
            # t = 0 is actually a correct solution
            newDuration = 0
            return newDuration
        log.debug('Unable to calculate the least upper bound: T0 = {0}; T1 = {1}; T2 = {2}; T3 = {3}'.\
                  format(mp.nstr(T0, n=_prec), mp.nstr(T1, n=_prec), mp.nstr(T2, n=_prec), mp.nstr(T3, n=_prec)))
        return number('-1')

    def __init__(self, npts):
        if not mp.isint(mp.sqrt(npts)):
            raise ValueError('Invalid number of points for quad rule')

        rulecls = subclass_where(BaseLineQuadRule, name=self.name)
        rule = rulecls(int(mp.sqrt(npts)))

        self.points = [(i, j) for j in rule.points for i in rule.points]

        if hasattr(rule, 'weights'):
            self.weights = [i*j for j in rule.weights for i in rule.weights]

def z_x123_frm_m(N, m):
    """Function to get x1, x2 and x3 (eq 3, 5 and 6, [McNamara93]_)."""
    M = -ellipk(m) / N
    snMM = ellipfun('sn', u= -M, m=m)
    snM = ellipfun('sn', u=M, m=m)
    cnM = ellipfun('cn', u=M, m=m)
    dnM = ellipfun('dn', u=M, m=m)
    znM = z_zn(M, m)
    x3 = snMM
    x1 = x3 * mp.sqrt(1 - m) / dnM
    x2 = x3 * mp.sqrt(1 - (cnM * znM) / (snM * dnM))  
    return x1, x2, x3

def _hex_orthob_at(order, p, q, r):
    sk = [mp.sqrt(k + 0.5) for k in xrange(order)]
    pa = [c*jp for c, jp in zip(sk, jacobi(order - 1, 0, 0, p))]
    pb = [c*jp for c, jp in zip(sk, jacobi(order - 1, 0, 0, q))]
    pc = [c*jp for c, jp in zip(sk, jacobi(order - 1, 0, 0, r))]

    return [pi*pj*pk for pi in pa for pj in pb for pk in pc]

    def ortho_basis_at_mp(self, p, q, r):
        sk = [mp.sqrt(k + 0.5) for k in range(self.order)]
        pa = [c * jp for c, jp in zip(sk, jacobi(self.order - 1, 0, 0, p))]
        pb = [c * jp for c, jp in zip(sk, jacobi(self.order - 1, 0, 0, q))]
        pc = [c * jp for c, jp in zip(sk, jacobi(self.order - 1, 0, 0, r))]

        return [pi * pj * pk for pi in pa for pj in pb for pk in pc]

def _pri_orthob_at(order, p, q, r):
    a = 2*(1 + p)/(1 - q) - 1 if q != 1 else 0
    b = q
    c = r

    pab = []
    for i, pi in enumerate(jacobi(order - 1, 0, 0, a)):
        ci = (1 - b)**i / 2**(i + 1)

        for j, pj in enumerate(jacobi(order - i - 1, 2*i + 1, 0, b)):
            cij = mp.sqrt((2*i + 1)*(2*i + 2*j + 2))*ci

            pab.append(cij*pi*pj)

    sk = [mp.sqrt(k + 0.5) for k in xrange(order)]
    pc = [s*jp for s, jp in zip(sk, jacobi(order - 1, 0, 0, c))]

    return [pij*pk for pij in pab for pk in pc]

def z_x123_frm_m(N, m):
    r"""
    Function to get x1, x2 and x3 (eq 3, 5 and 6, [McNamara93]_).
    
    :param N:    Order of the Zolotarev polynomial
    :param m:    m is the elliptic parameter (not the modulus k and not the nome q)
                  
    :rtype:      Returns [x1,x2,x3] ... where x1, x2, and x3 are defined in Fig. 1, [McNamara93]_
    """
    M = -ellipk(m) / N
    snMM = ellipfun('sn', u= -M, m=m)
    snM = ellipfun('sn', u=M, m=m)
    cnM = ellipfun('cn', u=M, m=m)
    dnM = ellipfun('dn', u=M, m=m)
    znM = z_zn(M, m)
    x3 = snMM
    x1 = x3 * mp.sqrt(1 - m) / dnM
    x2 = x3 * mp.sqrt(1 - (cnM * znM) / (snM * dnM))  
    return x1, x2, x3

    def ortho_basis_at_mp(self, p, q, r):
        q = q if q != 1 else q + mp.eps

        a = 2*(1 + p)/(1 - q) - 1
        b = q
        c = r

        pab = []
        for i, pi in enumerate(jacobi(self.order - 1, 0, 0, a)):
            ci = (1 - b)**i / 2**(i + 1)

            for j, pj in enumerate(jacobi(self.order - i - 1, 2*i + 1, 0, b)):
                cij = mp.sqrt((2*i + 1)*(2*i + 2*j + 2))*ci

                pab.append(cij*pi*pj)

        sk = [mp.sqrt(k + 0.5) for k in xrange(self.order)]
        pc = [s*jp for s, jp in zip(sk, jacobi(self.order - 1, 0, 0, c))]

        return [pij*pk for pij in pab for pk in pc]

def _tet_orthob_at(order, p, q, r):
    a = -2*(1 + p)/(q + r) - 1 if q + r != 0 else 0
    b = 2*(1 + q)/(1 - r) - 1 if r != 1 else 0
    c = r

    ob = []
    for i, pi in enumerate(jacobi(order - 1, 0, 0, a)):
        ci = mp.mpf(2)**(-2*i - 1.5)*mp.sqrt(4*i + 2)*(1 - b)**i

        for j, pj in enumerate(jacobi(order - i - 1, 2*i + 1, 0, b)):
            cj = mp.sqrt(i + j + 1)*2**-j*(1 - c)**(i + j)
            cij = ci*cj
            pij = pi*pj

            jp = jacobi(order - i - j - 1, 2*(i + j + 1), 0, c)
            for k, pk in enumerate(jp):
                ck = mp.sqrt(2*(k + j + i) + 3)

                ob.append(cij*ck*pij*pk)

    return ob

def _Interpolate1DNoVelocityLimit(x0, x1, v0, v1, am):
    # Check types
    if type(x0) is not mp.mpf:
        x0 = mp.mpf("{:.15e}".format(x0))
    if type(x1) is not mp.mpf:
        x1 = mp.mpf("{:.15e}".format(x1))
    if type(v0) is not mp.mpf:
        v0 = mp.mpf("{:.15e}".format(v0))
    if type(v1) is not mp.mpf:
        v1 = mp.mpf("{:.15e}".format(v1))
    if type(am) is not mp.mpf:
        am = mp.mpf("{:.15e}".format(am))

    # Check inputs
    assert(am > zero)

    # Check for an appropriate acceleration direction of the first ramp
    d = Sub(x1, x0)
    dv = Sub(v1, v0)
    difVSqr = Sub(v1**2, v0**2)
    
    if Abs(dv) < epsilon:
        if Abs(d) < epsilon:
            # Stationary ramp
            ramp0 = Ramp(zero, zero, zero, x0)
            return ParabolicCurve([ramp0])

        else:
            dStraight = zero
    else:    
        dStraight = mp.fdiv(difVSqr, Prod([2, mp.sign(dv), am]))
    
    if IsEqual(d, dStraight):
        # With the given distance, v0 and v1 can be directly connected using max/min
        # acceleration. Here the resulting profile has only one ramp.
        a0 = mp.sign(dv) * am
        ramp0 = Ramp(v0, a0, mp.fdiv(dv, a0), x0)
        return ParabolicCurve([ramp0])

    sumVSqr = Add(v0**2, v1**2)
    sigma = mp.sign(Sub(d, dStraight))
    a0 = sigma * am # acceleration of the first ramp
    vp = sigma * mp.sqrt(Add(Mul(pointfive, sumVSqr), Mul(a0, d)))
    t0 = mp.fdiv(Sub(vp, v0), a0)
    t1 = mp.fdiv(Sub(vp, v1), a0)
    ramp0 = Ramp(v0, a0, t0, x0)    
    assert(IsEqual(ramp0.v1, vp)) # check soundness
    ramp1 = Ramp(vp, Neg(a0), t1)

    curve = ParabolicCurve([ramp0, ramp1])
    assert(IsEqual(curve.d, d)) # check soundness
    return curve

def _tri_orthob_at(order, p, q):
    a = 2*(1 + p)/(1 - q) - 1 if q != 1 else 0
    b = q

    ob = []
    for i, pi in enumerate(jacobi(order - 1, 0, 0, a)):
        pa = pi*(1 - b)**i

        for j, pj in enumerate(jacobi(order - i - 1, 2*i + 1, 0, b)):
            cij = mp.sqrt((2*i + 1)*(2*i + 2*j + 2)) / 2**(i + 1)

            ob.append(cij*pa*pj)

    return ob

def z_Zolotarev(N, x, m):
    """Function to evaluate the Zolotarev polynomial (eq 1, [McNamara93]_)."""
    M = -ellipk(m) / N
    x3 = ellipfun('sn', u= -M, m=m)  
    xbar = x3 * mp.sqrt((x ** 2 - 1) / (x ** 2 - x3 ** 2)) # rearranged eq 21, [Levy70]_
    u = ellipf(mp.asin(xbar), m) # rearranged eq 20, [Levy70]_, asn(x) = F(asin(x)|m)     
    f = mp.cosh((N / 2) * mp.log(z_eta(M + u, m) / z_eta(M - u, m)))
    if (f.imag / f.real > 1e-10):
        print "imaginary part of the Zolotarev function is not negligible!"
        print "f_imaginary = ", f.imag
    else:
        if (x > 0): # no idea why I am doing this ... anyhow, it seems working
            f = -f.real  
        else:
            f = f.real        
    return f

def SolveQuartic(a, b, c, d, e):
    """
    SolveQuartic solves a quartic (fouth order) equation of the form
            ax^4 + bx^3 + cx^2 + dx + e = 0.
    For the detail of formulae presented here, see https://en.wikipedia.org/wiki/Quartic_function
    """
    # Check types
    if type(a) is not mp.mpf:
        a = mp.mpf("{:.15e}".format(a))
    if type(b) is not mp.mpf:
        b = mp.mpf("{:.15e}".format(b))
    if type(c) is not mp.mpf:
        c = mp.mpf("{:.15e}".format(c))
    if type(d) is not mp.mpf:
        d = mp.mpf("{:.15e}".format(d))
    if type(e) is not mp.mpf:
        e = mp.mpf("{:.15e}".format(e))

    """
    # Working code (more readable but probably less precise)
    p = (8*a*c - 3*b*b)/(8*a*a)
    q = (b**3 - 4*a*b*c + 8*a*a*d)/(8*a*a*a)
    delta0 = c*c - 3*b*d + 12*a*e
    delta1 = 2*(c**3) - 9*b*c*d + 27*b*b*e + 27*a*d*d - 72*a*c*e
    Q = mp.nthroot(pointfive*(delta1 + mp.sqrt(delta1*delta1 - 4*mp.power(delta0, 3))), 3)
    S = pointfive*mp.sqrt(-mp.fdiv(mp.mpf('2'), mp.mpf('3'))*p + (one/(3*a))*(Q + delta0/Q))

    x1 = -b/(4*a) - S + pointfive*mp.sqrt(-4*S*S - 2*p + q/S)
    x2 = -b/(4*a) - S - pointfive*mp.sqrt(-4*S*S - 2*p + q/S)
    x3 = -b/(4*a) + S + pointfive*mp.sqrt(-4*S*S - 2*p - q/S)
    x4 = -b/(4*a) + S - pointfive*mp.sqrt(-4*S*S - 2*p - q/S)
    """
    p = mp.fdiv(Sub(Prod([number('8'), a, c]), Mul(number('3'), mp.power(b, 2))), Mul(number('8'), mp.power(a, 2)))
    q = mp.fdiv(Sum([mp.power(b, 3), Prod([number('-4'), a, b, c]), Prod([number('8'), mp.power(a, 2), d])]), Mul(8, mp.power(a, 3)))
    delta0 = Sum([mp.power(c, 2), Prod([number('-3'), b, d]), Prod([number('12'), a, e])])
    delta1 = Sum([Mul(2, mp.power(c, 3)), Prod([number('-9'), b, c, d]), Prod([number('27'), mp.power(b, 2), e]), Prod([number('27'), a, mp.power(d, 2)]), Prod([number('-72'), a, c, e])])
    Q = mp.nthroot(Mul(pointfive, Add(delta1, mp.sqrt(Add(mp.power(delta1, 2), Mul(number('-4'), mp.power(delta0, 3)))))), 3)
    S = Mul(pointfive, mp.sqrt(Mul(mp.fdiv(mp.mpf('-2'), mp.mpf('3')), p) + Mul(mp.fdiv(one, Mul(number('3'), a)), Add(Q, mp.fdiv(delta0, Q)))))

    # log.debug("p = {0}".format(mp.nstr(p, n=_prec)))
    # log.debug("q = {0}".format(mp.nstr(q, n=_prec)))
    # log.debug("delta0 = {0}".format(mp.nstr(delta0, n=_prec)))
    # log.debug("delta1 = {0}".format(mp.nstr(delta1, n=_prec)))
    # log.debug("Q = {0}".format(mp.nstr(Q, n=_prec)))
    # log.debug("S = {0}".format(mp.nstr(S, n=_prec)))

    x1 = Sum([mp.fdiv(b, Mul(number('-4'), a)), Neg(S), Mul(pointfive, mp.sqrt(Sum([Mul(number('-4'), mp.power(S, 2)), Mul(number('-2'), p), mp.fdiv(q, S)])))])
    x2 = Sum([mp.fdiv(b, Mul(number('-4'), a)), Neg(S), Neg(Mul(pointfive, mp.sqrt(Sum([Mul(number('-4'), mp.power(S, 2)), Mul(number('-2'), p), mp.fdiv(q, S)]))))])
    x3 = Sum([mp.fdiv(b, Mul(number('-4'), a)), S, Mul(pointfive, mp.sqrt(Sum([Mul(number('-4'), mp.power(S, 2)), Mul(number('-2'), p), Neg(mp.fdiv(q, S))])))])
    x4 = Sum([mp.fdiv(b, Mul(number('-4'), a)), S, Neg(Mul(pointfive, mp.sqrt(Sum([Mul(number('-4'), mp.power(S, 2)), Mul(number('-2'), p), Neg(mp.fdiv(q, S))]))))])
    
    return [x1, x2, x3, x4]

def test_levin_3():
    mp.dps = 17
    z=mp.mpf(2)
    eps = mp.mpf(mp.eps)
    with mp.extraprec(7*mp.prec):  # we need copious amount of precision to sum this highly divergent series
        L = mp.levin(method = "levin", variant = "t")
        n, s = 0, 0
        while 1:
            s += (-z)**n * mp.fac(4 * n) / (mp.fac(n) * mp.fac(2 * n) * (4 ** n))
            n += 1
            v, e = L.step_psum(s)
            if e < eps:
                break
            if n > 1000: raise RuntimeError("iteration limit exceeded")
    eps = mp.exp(0.8 * mp.log(eps))
    exact = mp.quad(lambda x: mp.exp( -x * x / 2 - z * x ** 4), [0,mp.inf]) * 2 / mp.sqrt(2 * mp.pi)
    # there is also a symbolic expression for the integral:
    #   exact = mp.exp(mp.one / (32 * z)) * mp.besselk(mp.one / 4, mp.one / (32 * z)) / (4 * mp.sqrt(z * mp.pi))
    err = abs(v - exact)
    assert err < eps
    w = mp.nsum(lambda n: (-z)**n * mp.fac(4 * n) / (mp.fac(n) * mp.fac(2 * n) * (4 ** n)), [0, mp.inf], method = "levin", levin_variant = "t", workprec = 8*mp.prec, steps = [2] + [1 for x in xrange(1000)])
    err = abs(v - w)
    assert err < eps

def norm_q(prob):
    r"""
    A multi-precision calculation of the
    standard normal quantile function:

    .. math::

       \int_{-\infty}^{q(p)} \frac{e^{-z^2/2}}{\sqrt{2\pi}} \; dz = p

    where $p$ is `prob`.

    Parameters
    ----------

    prob : float

    Returns
    -------

    quantile : float

    """
    return np.array(mp.erfinv(2*prob-1)*mp.sqrt(2))

def z_Zolotarev(N, x, m):
    r"""
    Function to evaluate the Zolotarev polynomial (eq 1, [McNamara93]_).
    
    :param N:    Order of the Zolotarev polynomial
    :param x:    The argument at which one would like to evaluate the Zolotarev polynomial
    :param m:    m is the elliptic parameter (not the modulus k and not the nome q)
                  
    :rtype:      Returns a float, the value of Zolotarev polynomial at x
    """
    M = -ellipk(m) / N
    x3 = ellipfun('sn', u= -M, m=m)  
    xbar = x3 * mp.sqrt((x ** 2 - 1) / (x ** 2 - x3 ** 2)) # rearranged eq 21, [Levy70]_
    u = ellipf(mp.asin(xbar), m) # rearranged eq 20, [Levy70]_, asn(x) = F(asin(x)|m)     
    f = mp.cosh((N / 2) * mp.log(z_eta(M + u, m) / z_eta(M - u, m)))
    if (f.imag / f.real > 1e-10):
        print "imaginary part of the Zolotarev function is not negligible!"
        print "f_imaginary = ", f.imag
    else:
        if (x > 0): # no idea why I am doing this ... anyhow, it seems working
            f = -f.real  
        else:
            f = f.real        
    return f