<SYSTEM_TASK:>
Solve the following problem using Python, implementing the functions described below, one line at a time
<END_TASK>
<USER_TASK:>
Description:
 def orth_ttr( order, dist, normed=False, sort="GR", retall=False, cross_truncation=1., **kws):
 """ Create orthogonal polynomial expansion from three terms recursion formula. Args: order (int):
 Order of polynomial expansion. dist (Dist):
 Distribution space where polynomials are orthogonal If dist.ttr exists, it will be used. Must be stochastically independent. normed (bool):
 If True orthonormal polynomials will be used. sort (str):
 Polynomial sorting. Same as in basis. retall (bool):
 If true return numerical stabilized norms as well. Roughly the same as ``cp.E(orth**2, dist)``. cross_truncation (float):
 Use hyperbolic cross truncation scheme to reduce the number of terms in expansion. only include terms where the exponents ``K`` satisfied the equation ``order >= sum(K**(1/cross_truncation))**cross_truncation``. Returns: (Poly, numpy.ndarray):
 Orthogonal polynomial expansion and norms of the orthogonal expansion on the form ``E(orth**2, dist)``. Calculated using recurrence coefficients for stability. Examples: [1.0, q0, q0^2-1.0, q0^3-3.0q0, q0^4-6.0q0^2+3.0] """

    polynomials, norms, _, _ = chaospy.quad.generate_stieltjes(
        dist=dist, order=numpy.max(order), retall=True, **kws)

    if normed:
        for idx, poly in enumerate(polynomials):
            polynomials[idx] = poly / numpy.sqrt(norms[:, idx])
        norms = norms**0

    dim = len(dist)
    if dim > 1:
        mv_polynomials = []
        mv_norms = []
        indices = chaospy.bertran.bindex(
            start=0, stop=order, dim=dim, sort=sort,
            cross_truncation=cross_truncation,
        )

        for index in indices:
            poly = polynomials[index[0]][0]
            for idx in range(1, dim):
                poly = poly * polynomials[index[idx]][idx]
            mv_polynomials.append(poly)

        if retall:
            for index in indices:
                mv_norms.append(
                    numpy.prod([norms[idx, index[idx]] for idx in range(dim)]))

    else:
        mv_norms = norms[0]
        mv_polynomials = polynomials

    polynomials = chaospy.poly.flatten(chaospy.poly.Poly(mv_polynomials))

    if retall:
        return polynomials, numpy.array(mv_norms)
    return polynomials