Transverse_velocity_gradient_moments.ipynb

Kernel: Python 3 (CoCalc)

from __future__ import annotations

import random
from functools import lru_cache
from typing import List, Sequence, Tuple, Optional

import sympy as sp
from sympy import ZZ
from sympy.polys.matrices import DM


# ======================================================================================
#  Exact orientational moment for transverse gradient:
#    A_perp = u_i A_{ij} v_j
#    u ~ uniform on S^2, v ~ uniform on unit circle in plane \perp u
#
#  Compute exactly (SymPy Rational):
#    < (u^T A v)^n >_orien   (odd n = 0)
#
#  Then decompose into paper invariants:
#    I1 = tr(A)
#    I2 = tr(A^2)
#    I3 = tr(A A^T)
#    I4 = tr(A^3)
#    I5 = tr(A A A^T)
#    I6 = tr(A A A^T A^T)
#
#  Return exact polynomial in I1..I6 and format as:
#    sum coeff * <monomial(I)>_space
#
#  Memory-friendly exact solve:
#    - choose K independent equations mod p (incremental)
#    - solve integer system using DomainMatrix.solve_den (fraction-free)
# ======================================================================================


# ============================================================
# Utilities: double factorial, sphere monomial moments, circle constant
# ============================================================

@lru_cache(None)
def double_factorial(k: int) -> int:
    """k!! for integer k >= -1. Convention: (-1)!!=1, 0!!=1."""
    if k < -1:
        raise ValueError("double factorial needs k >= -1")
    if k in (0, 1):
        return 1
    out = 1
    x = k
    while x >= 2:
        out *= x
        x -= 2
    return out


@lru_cache(None)
def sphere_monomial_moment(exps: Tuple[int, int, int]) -> sp.Rational:
    """
    Exact moment for u ~ Uniform(S^2):
      E[u1^a u2^b u3^c]
    If any exponent is odd -> 0.
    If all even: a=2A,b=2B,c=2C:
      E = ( (2A-1)!! (2B-1)!! (2C-1)!! ) / ( (2(A+B+C)+1)!! )
    """
    a, b, c = exps
    if (a | b | c) & 1:
        return sp.Rational(0, 1)
    A = a // 2
    B = b // 2
    C = c // 2
    num = (
        double_factorial(2 * A - 1)
        * double_factorial(2 * B - 1)
        * double_factorial(2 * C - 1)
    )
    den = double_factorial(2 * (A + B + C) + 1)
    return sp.Rational(num, den)


@lru_cache(None)
def K_circle_even_moment(m: int) -> sp.Rational:
    """
    For v uniform on unit circle in 2D plane and fixed planar vector a:
      E[(a·v)^{2m}] = K_m * ||a||^{2m}
    where:
      K_m = (2m)! / (2^{2m} (m!)^2)
    """
    return sp.Rational(sp.factorial(2 * m), 2 ** (2 * m) * sp.factorial(m) ** 2)


def common_denominator_for_moment(n: int) -> int:
    """
    Safe common denominator D(n) so that D(n) * <(u^T A v)^n>_orien is integer
    for all integer A, when n is even.

    n=2m:
      <...> = K_m * E_u[(q(u))^m]
      denom(K_m) divides 2^{2m}(m!)^2
      denom(E_u[poly degree 4m]) divides (4m+1)!!
    so take:
      D = 2^{2m} (m!)^2 (4m+1)!!
    """
    if n % 2 == 1:
        return 1
    m = n // 2
    return (2 ** (2 * m)) * int(sp.factorial(m) ** 2) * double_factorial(4 * m + 1)


# ============================================================
# Exact orientational average: <(u^T A v)^n>_orien
# ============================================================

_u1, _u2, _u3 = sp.symbols("u1 u2 u3")
_uvec = sp.Matrix([_u1, _u2, _u3])


def moment_uAv_power_exact(A_int: Sequence[Sequence[int]], n: int) -> sp.Rational:
    """
    Exact orientational average:
        <(u^T A v)^n>_orien
    with u~S^2 uniform, v~unit circle in plane ⟂ u uniform.
    Odd n -> 0.
    """
    if n < 0:
        raise ValueError("n must be nonnegative")
    if n % 2 == 1:
        return sp.Rational(0, 1)
    if n == 0:
        return sp.Rational(1, 1)

    m = n // 2
    A = sp.Matrix([[int(A_int[i][j]) for j in range(3)] for i in range(3)])

    # q(u) = u^T (A A^T) u - (u^T A^T u)^2
    Q = A * A.T
    term1 = (_uvec.T * Q * _uvec)[0]
    term2 = (_uvec.T * A.T * _uvec)[0]
    q = sp.expand(term1 - term2 ** 2)

    poly_q = sp.Poly(q, _u1, _u2, _u3, domain="ZZ")
    poly_qm = poly_q ** m

    Eu = sp.Rational(0, 1)
    for (a, b, c), coeff in poly_qm.terms():
        Eu += sp.Rational(int(coeff), 1) * sphere_monomial_moment((a, b, c))

    return sp.simplify(K_circle_even_moment(m) * Eu)


# ============================================================
# Paper invariants I1..I6
# ============================================================

I1, I2, I3, I4, I5, I6 = sp.symbols("I1 I2 I3 I4 I5 I6")
InvariantSyms = (I1, I2, I3, I4, I5, I6)
InvariantDegrees = (1, 2, 2, 3, 3, 4)


def invariant_blocks(A: sp.Matrix) -> Tuple[sp.Expr, ...]:
    """
    Paper (Eq. 3.1):
      I1 = tr(A)
      I2 = tr(A^2)
      I3 = tr(A A^T)
      I4 = tr(A^3)
      I5 = tr(A A A^T)
      I6 = tr(A A A^T A^T)
    """
    AT = A.T
    return (
        sp.trace(A),
        sp.trace(A * A),
        sp.trace(A * AT),
        sp.trace(A * A * A),
        sp.trace(A * A * AT),
        sp.trace(A * A * AT * AT),
    )


def degree_ntuples(weight: int) -> List[Tuple[int, int, int, int, int, int]]:
    """
    All 6-tuples alpha=(a1..a6) with weighted degree:
      a1 + 2a2 + 2a3 + 3a4 + 3a5 + 4a6 = weight
    """
    if weight < 0:
        raise ValueError("weight must be nonnegative")
    d1, d2, d3, d4, d5, d6 = InvariantDegrees
    out: List[Tuple[int, int, int, int, int, int]] = []
    for a6 in range(weight // d6 + 1):
        for a5 in range((weight - d6 * a6) // d5 + 1):
            for a4 in range((weight - d6 * a6 - d5 * a5) // d4 + 1):
                for a3 in range((weight - d6 * a6 - d5 * a5 - d4 * a4) // d3 + 1):
                    for a2 in range((weight - d6 * a6 - d5 * a5 - d4 * a4 - d3 * a3) // d2 + 1):
                        rem = weight - d6 * a6 - d5 * a5 - d4 * a4 - d3 * a3 - d2 * a2
                        if rem >= 0 and rem % d1 == 0:
                            a1 = rem // d1
                            out.append((a1, a2, a3, a4, a5, a6))
    return out


def build_basis_tuples(weight: int, incompressible: bool = False) -> List[Tuple[int, int, int, int, int, int]]:
    t = degree_ntuples(weight)
    if incompressible:
        t = [a for a in t if a[0] == 0]
    return t


def monomial_symb(alpha: Sequence[int]) -> sp.Expr:
    expr = sp.Integer(1)
    for sym, p in zip(InvariantSyms, alpha):
        p = int(p)
        if p:
            expr *= sym ** p
    return expr


# ============================================================
# Random integer matrices
# ============================================================

def random_matrix(d: int = 3, low: int = -2, high: int = 2) -> List[List[int]]:
    return [[random.randint(low, high) for _ in range(d)] for _ in range(d)]


def random_trace_free_matrix(d: int = 3, low: int = -2, high: int = 2) -> List[List[int]]:
    M = [[random.randint(low, high) for _ in range(d)] for _ in range(d)]
    tr_first = sum(M[i][i] for i in range(d - 1))
    M[d - 1][d - 1] = -tr_first
    return M


# ============================================================
# Mod-p incremental independence selector (memory-light)
# ============================================================

class ModRankSelector:
    """Incremental row-independence test over GF(p), storing echelon basis."""
    def __init__(self, ncols: int, p: int):
        self.ncols = ncols
        self.p = p
        self.pivots: List[int] = []
        self.basis: List[List[int]] = []

    def try_add(self, row: Sequence[int]) -> bool:
        p = self.p
        v = [int(x) % p for x in row]

        for piv, br in zip(self.pivots, self.basis):
            if v[piv] != 0:
                f = v[piv]
                for j in range(piv, self.ncols):
                    v[j] = (v[j] - f * br[j]) % p

        for piv in range(self.ncols):
            if v[piv] != 0:
                inv = pow(v[piv], -1, p)
                for j in range(piv, self.ncols):
                    v[j] = (v[j] * inv) % p
                self.pivots.append(piv)
                self.basis.append(v)
                return True

        return False


# ============================================================
# Exact verification
# ============================================================

def verify_identity(
    expr_in_I: sp.Expr,
    n: int,
    *,
    tests: int,
    incompressible: bool,
    seed: int = 0,
    low: int = -4,
    high: int = 4,
) -> None:
    random.seed(seed)
    for _ in range(tests):
        A_int = random_trace_free_matrix(3, low=low, high=high) if incompressible else random_matrix(3, low=low, high=high)
        A = sp.Matrix(A_int)
        Ivals = invariant_blocks(A)
        rhs = sp.expand(expr_in_I.subs({s: v for s, v in zip(InvariantSyms, Ivals)}))
        lhs = moment_uAv_power_exact(A_int, n)
        if sp.simplify(lhs - rhs) != 0:
            raise RuntimeError(f"Verification failed.\nA={A_int}\nLHS={lhs}\nRHS={rhs}")


# ============================================================
# Decomposition: exact, memory-friendly
#   Solve M * (D*coeffs) = (D*y) using DomainMatrix.solve_den
# ============================================================

def decompose_to_invariant_basis_exact(
    n: int,
    *,
    incompressible: bool = False,
    seed: int = 1234,
    low: int = -2,
    high: int = 2,
    max_tries: int = 30,
    row_cap_factor: int = 14,
    verify_tests: int = 3,
) -> sp.Expr:
    if n % 2 == 1:
        return sp.Integer(0)

    basis_tuples = build_basis_tuples(n, incompressible=incompressible)
    K = len(basis_tuples)
    if K == 0:
        return sp.Integer(0)

    # evaluation plan: monomial = Π I[idx]^pow
    monom_plan: List[List[Tuple[int, int]]] = []
    for alpha in basis_tuples:
        plan = [(idx, int(p)) for idx, p in enumerate(alpha) if p]
        monom_plan.append(plan)

    D = common_denominator_for_moment(n)

    # large prime for independence test
    p = 1_000_000_007

    random.seed(seed)

    for attempt in range(1, max_tries + 1):
        selector = ModRankSelector(ncols=K, p=p)
        rows_int: List[List[int]] = []
        rhs_int: List[int] = []

        cap = row_cap_factor * (K + 5)
        tries = 0

        while len(rows_int) < K and tries < cap:
            tries += 1
            A_int = random_trace_free_matrix(3, low=low, high=high) if incompressible else random_matrix(3, low=low, high=high)
            A = sp.Matrix(A_int)

            # invariants are integers for integer A
            Ivals = invariant_blocks(A)
            Iint = [int(v) for v in Ivals]

            row: List[int] = []
            for plan in monom_plan:
                val = 1
                for idx, pow_ in plan:
                    val *= pow(Iint[idx], pow_)
                row.append(val)

            if not selector.try_add(row):
                continue

            y = moment_uAv_power_exact(A_int, n)  # Rational
            b_num = int(y.p) * int(D)
            b_den = int(y.q)
            if b_num % b_den != 0:
                # should not happen with this D, but keep safe
                continue
            b = b_num // b_den

            rows_int.append(row)
            rhs_int.append(b)

        if len(rows_int) < K:
            continue

        try:
            # Solve integer system M * x = b in ZZ, returning x = xnum/xden exactly.
            A_dm = DM([[ZZ(v) for v in row] for row in rows_int], ZZ)
            b_dm = DM([[ZZ(v)] for v in rhs_int], ZZ)

            xnum_dm, xden = A_dm.solve_den(b_dm)  # x = xnum / xden
            xnum_mat = xnum_dm.to_Matrix()
            xden_int = int(xden)

            x_scaled = [sp.Rational(int(xnum_mat[i, 0]), xden_int) for i in range(K)]  # x = D*coeffs
            coeffs = [sp.simplify(x / sp.Integer(D)) for x in x_scaled]

            expr = sp.expand(sum(c * monomial_symb(a) for c, a in zip(coeffs, basis_tuples)))
            if incompressible:
                expr = sp.simplify(expr.subs(I1, 0))

            if verify_tests > 0:
                verify_identity(
                    expr, n,
                    tests=verify_tests,
                    incompressible=incompressible,
                    seed=seed + 10_000 + attempt,
                    low=low, high=high
                )

            return sp.simplify(expr)

        except Exception:
            # try another batch
            continue

    raise RuntimeError(
        "Failed to find a full-rank sample set / stable solve. "
        "Try changing seed, increasing max_tries/row_cap_factor, or widening low/high."
    )


# ============================================================
# Formatting as ensemble/space averages: coeff * <monomial>
# ============================================================

class AngleAverage(sp.Function):
    nargs = 1

    @classmethod
    def eval(cls, arg):
        return None


def _angleaverage_str(self, printer):
    return "<" + printer.doprint(self.args[0]) + ">"


AngleAverage._sympystr = _angleaverage_str  # type: ignore


def format_invariant_moment(expr: sp.Expr) -> sp.Expr:
    expr = sp.expand(expr)
    out = sp.Integer(0)
    for term in expr.as_ordered_terms():
        coeff, rest = term.as_coeff_Mul()
        out += sp.simplify(coeff) * AngleAverage(sp.simplify(rest))
    return sp.simplify(out)


def enforce_betchov_second_order(expr: sp.Expr, incompressible: bool = False) -> sp.Expr:
    """
    Optional: apply Betchov/homogeneity constraint ONLY for n=2:
      incompressible: <I2> = 0
      compressible:   <I2> = <I1^2>
    (Here I2 = tr(A^2), I1 = tr(A).)
    """
    if incompressible:
        return sp.simplify(expr.subs(I2, 0))
    return sp.simplify(expr.subs(I2, I1 ** 2))


def A_perp_moment(
    n: int,
    *,
    incompressible: bool = False,
    seed: int = 1234,
    low: int = -2,
    high: int = 2,
    max_tries: int = 30,
    row_cap_factor: int = 14,
    verify_tests: int = 3,
    use_betchov: bool = False,
) -> sp.Expr:
    """
    Return analytical expression for:
      < (u^T A v)^n >_space = << (u^T A v)^n >_orien >_space
    in terms of invariant averages <monomial(I1..I6)>.

    Odd n -> 0.
    """
    if n % 2 == 1:
        return sp.Integer(0)

    poly = decompose_to_invariant_basis_exact(
        n,
        incompressible=incompressible,
        seed=seed,
        low=low,
        high=high,
        max_tries=max_tries,
        row_cap_factor=row_cap_factor,
        verify_tests=verify_tests,
    )

    if n == 2 and use_betchov:
        poly = enforce_betchov_second_order(poly, incompressible=incompressible)

    return format_invariant_moment(poly)


# ============================================================
# Demo
# ============================================================

if __name__ == "__main__":
    # Paper’s even moments <A_perp^{2n_paper}> correspond to calling here with n = 2*n_paper.
    # Example: paper n_paper=5 -> here n=10.

    for n in [2, 8]:
        print(f"\n=== n={n} incompressible (tr A = 0) ===")
        expr_i = A_perp_moment(
            n,
            incompressible=True,
            low=-2, high=2,
            verify_tests=3,     
            max_tries=40,
            row_cap_factor=16,
            use_betchov=(n == 2),
        )
        print(expr_i)
        
        print(f"\n=== n={n} compressible ===") # upto 16
        expr_c = A_perp_moment(
            n,
            incompressible=False,
            low=-2, high=2,
            verify_tests=3,    
            max_tries=40,
            row_cap_factor=16,
            use_betchov=(n == 2),
        )
        print(expr_c)

=== n=2 incompressible (tr A = 0) ===
2*<I3>/15

=== n=2 compressible ===
-<I1**2>/15 + 2*<I3>/15

=== n=8 incompressible (tr A = 0) ===
751*<I2**4>/350064 + 490*<I3**4>/21879 + 5*<I6**2>/143 - 245*<I2*I3**3>/43758 - 8*<I2*I4**2>/7293 + 10*<I2*I5**2>/2431 + 55*<I2**2*I3**2>/2652 - 5*<I2**2*I6>/286 - 235*<I2**3*I3>/87516 + 20*<I3*I4**2>/2431 - 140*<I3*I5**2>/7293 - 35*<I3**2*I6>/429 + 5*<I2*I3*I6>/429 + 20*<I2*I4*I5>/7293 + 20*<I3*I4*I5>/7293

=== n=8 compressible ===
703*<I1**8>/350064 + 751*<I2**4>/350064 + 490*<I3**4>/21879 + 5*<I6**2>/143 - 122*<I1**2*I2**3>/21879 - 245*<I1**2*I3**3>/3366 - 4*<I1**2*I4**2>/1683 + 490*<I1**2*I5**2>/7293 + 853*<I1**4*I2**2>/175032 + 5365*<I1**4*I3**2>/87516 - 5*<I1**4*I6>/286 - 32*<I1**5*I4>/21879 + 50*<I1**5*I5>/1989 - 217*<I1**6*I2>/43758 - 1675*<I1**6*I3>/87516 - 245*<I2*I3**3>/43758 - 8*<I2*I4**2>/7293 + 10*<I2*I5**2>/2431 + 55*<I2**2*I3**2>/2652 - 5*<I2**2*I6>/286 - 235*<I2**3*I3>/87516 + 20*<I3*I4**2>/2431 - 140*<I3*I5**2>/7293 - 35*<I3**2*I6>/429 + 28*<I1*I2**2*I4>/21879 + 10*<I1*I2**2*I5>/561 - 40*<I1*I3**2*I4>/21879 + 980*<I1*I3**2*I5>/7293 - 40*<I1*I5*I6>/429 - 190*<I1**2*I2*I3**2>/7293 + 10*<I1**2*I2*I6>/429 - 5*<I1**2*I2**2*I3>/884 + 35*<I1**2*I3*I6>/429 - 40*<I1**2*I4*I5>/21879 + 40*<I1**3*I2*I4>/7293 - 20*<I1**3*I2*I5>/663 + 40*<I1**3*I3*I4>/7293 - 300*<I1**3*I3*I5>/2431 + 5*<I1**4*I2*I3>/204 + 5*<I2*I3*I6>/429 + 20*<I2*I4*I5>/7293 + 20*<I3*I4*I5>/7293 - 140*<I1*I2*I3*I4>/7293 - 80*<I1*I2*I3*I5>/7293

#Formulas in terms of invariants based on S and W
I1p, I2p, I3p, I4p, I5p, I6p = sp.symbols("I1' I2' I3' I4' I5' I6'")

def subs_I_to_Ip():
    return {
        I1: I1p,
        I2: I2p + I3p,
        I3: I2p - I3p,
        I4: I4p + 3*I5p,
        I5: I4p - I5p,
        I6: (sp.Rational(1,6)*I1p**4
             - I1p**2*I2p
             + sp.Rational(1,2)*I2p**2
             + sp.Rational(4,3)*I1p*I4p
             + sp.Rational(1,2)*I3p**2
             + I3p*(I1p**2 - I2p)
             + 4*I6p
             - 4*I1p*I5p),
    }

def poly_I_to_Ip(poly_in_I: sp.Expr) -> sp.Expr:
    return sp.expand(poly_in_I.subs(subs_I_to_Ip()))


# ============================================================
# Betchov constraints (prime basis) -- ONLY for second order (n=2)
#   From identities:
#     I1 = I1'
#     I2 = I2' + I3'
#   Paper constraints used in your I-basis code at n=2:
#     incompressible: <I2' + I3'> = 0  -> eliminate I3': I3' -> -I2'
#     compressible:   <I2' + I3' - I1'^2>=0 -> eliminate I3': I3' -> I1'^2 - I2'
#   This elimination is valid as a symbolic substitution only because n=2 result is linear in invariants.
# ============================================================
def enforce_betchov_second_order_prime(poly_Ip: sp.Expr, *, incompressible: bool) -> sp.Expr:
    if incompressible:
        # tr(A)=0 -> I1'=0, and <I2'+I3'>=0 -> eliminate I3'
        poly_Ip = sp.expand(poly_Ip.subs({I1p: 0}))
        return sp.simplify(sp.expand(poly_Ip.subs({I3p: -I2p})))
    # compressible: <I2'+I3' - I1'^2>=0 -> eliminate I3'
    return sp.simplify(sp.expand(poly_Ip.subs({I3p: I1p**2 - I2p})))


def A_perp_moment_prime_via_substitution(
    n: int,
    *,
    incompressible: bool = False,
    seed: int = 1234,
    low: int = -2,
    high: int = 2,
    max_tries: int = 30,
    row_cap_factor: int = 14,
    verify_tests: int = 3,
    use_betchov: bool = True,
) -> sp.Expr:
    # 1) Use previous I-basis decomposition
    poly_I = decompose_to_invariant_basis_exact(
        n,
        incompressible=incompressible,
        seed=seed,
        low=low, high=high,
        max_tries=max_tries,
        row_cap_factor=row_cap_factor,
        verify_tests=verify_tests,
    )

    # 2) Keep original second-order Betchov in I-basis
    #    This enforces <I2>=0 or <I2>=<I1^2>, but after converting to prime,
    #    you'll still have I3' floating around. So we ALSO do prime elimination below.
    if n == 2:
        poly_I = enforce_betchov_second_order(poly_I, incompressible=incompressible)

    # 3) Convert polynomial from I to I'
    poly_Ip = poly_I_to_Ip(poly_I)

    # 4) incompressible: I1' = 0 at the polynomial level
    if incompressible:
        poly_Ip = sp.simplify(poly_Ip.subs(I1p, 0))

    # 5)For n=2 only: apply prime-basis Betchov to eliminate I3'
    if n == 2 and use_betchov:
        poly_Ip = enforce_betchov_second_order_prime(poly_Ip, incompressible=incompressible)

    # 6) Wrap each monomial as <...>
    return format_invariant_moment(poly_Ip)


# ============================================================
# Demo / output
# ============================================================
if __name__ == "__main__":
    for n in [2, 4, 6]:
        # second order: show Betchov-simplified prime output
        print(f"\n=== n={n} incompressible ===")
        print(A_perp_moment_prime_via_substitution(n, incompressible=True))
  
        print(f"\n=== n={n} compressible ===")
        print(A_perp_moment_prime_via_substitution(n, incompressible=False))

=== n=2 incompressible (tr A = 0) ===
4*<I2'>/15

=== n=2 compressible (tr A = 0) ===
-<I1'**2>/5 + 4*<I2'>/15

=== n=4 incompressible (tr A = 0) ===
3*<I2'**2>/140 + <I3'**2>/20 - 12*<I6'>/35 + <I2'*I3'>/70

=== n=4 compressible (tr A = 0) ===
<I1'**4>/420 + 3*<I2'**2>/140 + <I3'**2>/20 - 12*<I6'>/35 + 8*<I1'*I5'>/35 - <I1'**2*I2'>/70 - 3*<I1'**2*I3'>/70 + <I2'*I3'>/70

=== n=6 incompressible (tr A = 0) ===
7*<I2'**3>/1144 - <I3'**3>/56 - 3*<I4'**2>/1001 + <I5'**2>/21 - 3*<I2'*I3'**2>/56 - 24*<I2'*I6'>/77 + 31*<I2'**2*I3'>/616 + 8*<I3'*I6'>/21 + 6*<I4'*I5'>/77

=== n=6 compressible (tr A = 0) ===
-3*<I1'**6>/8008 + 7*<I2'**3>/1144 - <I3'**3>/56 - 3*<I4'**2>/1001 + <I5'**2>/21 - 73*<I1'**2*I2'**2>/8008 + 11*<I1'**2*I3'**2>/168 + 8*<I1'**2*I6'>/77 - 4*<I1'**3*I4'>/3003 - 4*<I1'**3*I5'>/77 + 27*<I1'**4*I2'>/8008 + <I1'**4*I3'>/88 - 3*<I2'*I3'**2>/56 - 24*<I2'*I6'>/77 + 31*<I2'**2*I3'>/616 + 8*<I3'*I6'>/21 + 6*<I4'*I5'>/77 + 6*<I1'*I2'*I4'>/1001 + 10*<I1'*I2'*I5'>/77 - 2*<I1'*I3'*I4'>/77 - 2*<I1'*I3'*I5'>/7 - 13*<I1'**2*I2'*I3'>/308