ENH: [f2py] Symbolic solver for dimension specifications.

1 files changed, 1258 insertions, 0 deletions

diff --git a/numpy/f2py/symbolic.py b/numpy/f2py/symbolic.py
new file mode 100644
index 000000000..d303da0a0
--- /dev/null
+++ b/numpy/f2py/symbolic.py
@@ -0,0 +1,1258 @@
+"""Fortran/C symbolic expressions
+
+References:
+- J3/21-007: Draft Fortran 202x. https://j3-fortran.org/doc/year/21/21-007.pdf
+"""
+
+# To analyze Fortran expressions to solve dimensions specifications,
+# for instances, we implement a minimal symbolic engine for parsing
+# expressions into a tree of expression instances. As a first
+# instance, we care only about arithmetic expressions involving
+# integers and operations like addition (+), subtraction (-),
+# multiplication (*), division (Fortran / is Python //, Fortran // is
+# concatenate), and exponentiation (**).  In addition, .pyf files may
+# contain C expressions that support here is implemented as well.
+#
+# TODO: support logical constants (Op.BOOLEAN)
+# TODO: support logical operators (.AND., ...)
+# TODO: support relational operators (<, >, ..., .LT., ...)
+# TODO: support defined operators (.MYOP., ...)
+#
+__all__ = ['Expr']
+
+
+import re
+import warnings
+from enum import Enum
+from math import gcd
+
+
+class Language(Enum):
+    """
+    Used as Expr.tostring language argument.
+    """
+    Python = 0
+    Fortran = 1
+    C = 2
+
+
+class Op(Enum):
+    """
+    Used as Expr op attribute.
+    """
+    INTEGER = 10
+    REAL = 12
+    COMPLEX = 15
+    STRING = 20
+    ARRAY = 30
+    SYMBOL = 40
+    APPLY = 200
+    INDEXING = 210
+    CONCAT = 220
+    TERMS = 1000
+    FACTORS = 2000
+
+
+class ArithOp(Enum):
+    """
+    Used in Op.APPLY expression to specify the function part.
+    """
+    POS = 1
+    NEG = 2
+    ADD = 3
+    SUB = 4
+    MUL = 5
+    DIV = 6
+    POW = 7
+
+
+class OpError(Exception):
+    pass
+
+
+class Precedence(Enum):
+    """
+    Used as Expr.tostring precedence argument.
+    """
+    NONE = 0
+    TUPLE = 1
+    SUM = 2
+    PRODUCT = 3
+    POWER = 4
+    ATOM = 5
+
+
+integer_types = (int,)
+number_types = (int, float)
+
+
+def _pairs_add(d, k, v):
+    # Internal utility method for updating terms and factors data.
+    c = d.get(k)
+    if c is None:
+        d[k] = v
+    else:
+        c = c + v
+        if c:
+            d[k] = c
+        else:
+            del d[k]
+
+
+class Expr:
+    """Represents a Fortran expression as a op-data pair.
+
+    Expr instances are hashable and sortable.
+    """
+
+    @staticmethod
+    def parse(s):
+        """Parse a Fortran expression to a Expr.
+        """
+        return fromstring(s)
+
+    def __init__(self, op, data):
+        assert isinstance(op, Op)
+
+        # sanity checks
+        if op is Op.INTEGER:
+            # data is a 2-tuple of numeric object and a kind value
+            # (default is 4)
+            assert isinstance(data, tuple) and len(data) == 2
+            assert isinstance(data[0], int)
+            assert isinstance(data[1], (int, str)), data
+        elif op is Op.REAL:
+            # data is a 2-tuple of numeric object and a kind value
+            # (default is 4)
+            assert isinstance(data, tuple) and len(data) == 2
+            assert isinstance(data[0], float)
+            assert isinstance(data[1], (int, str)), data
+        elif op is Op.COMPLEX:
+            # data is a 2-tuple of constant expressions
+            assert isinstance(data, tuple) and len(data) == 2
+        elif op is Op.STRING:
+            # data is a 2-tuple of quoted string and a kind value
+            # (default is 1)
+            assert isinstance(data, tuple) and len(data) == 2
+            assert (isinstance(data[0], str)
+                    and data[0][::len(data[0])-1] in ('""', "''", '@@'))
+            assert isinstance(data[1], (int, str)), data
+        elif op is Op.SYMBOL:
+            # data is any hashable object
+            assert hash(data) is not None
+        elif op in (Op.ARRAY, Op.CONCAT):
+            # data is a tuple of expressions
+            assert isinstance(data, tuple)
+            assert all(isinstance(item, Expr) for item in data), data
+        elif op in (Op.TERMS, Op.FACTORS):
+            # data is {<term|base>:<coeff|exponent>} where dict values
+            # are nonzero Python integers
+            assert isinstance(data, dict)
+        elif op is Op.APPLY:
+            # data is (<function>, <operands>, <kwoperands>) where
+            # operands are Expr instances
+            assert isinstance(data, tuple) and len(data) == 3
+            # function is any hashable object
+            assert hash(data[0]) is not None
+            assert isinstance(data[1], tuple)
+            assert isinstance(data[2], dict)
+        elif op is Op.INDEXING:
+            # data is (<object>, <indices>)
+            assert isinstance(data, tuple) and len(data) == 2
+            # function is any hashable object
+            assert hash(data[0]) is not None
+        else:
+            raise NotImplementedError(
+                f'unknown op or missing sanity check: {op}')
+
+        self.op = op
+        self.data = data
+
+    def __eq__(self, other):
+        return (isinstance(other, Expr)
+                and self.op is other.op
+                and self.data == other.data)
+
+    def __hash__(self):
+        if self.op in (Op.TERMS, Op.FACTORS):
+            data = tuple(sorted(self.data.items()))
+        elif self.op is Op.APPLY:
+            data = self.data[:2] + tuple(sorted(self.data[2].items()))
+        else:
+            data = self.data
+        return hash((self.op, data))
+
+    def __lt__(self, other):
+        if isinstance(other, Expr):
+            if self.op is not other.op:
+                return self.op.value < other.op.value
+            if self.op in (Op.TERMS, Op.FACTORS):
+                return (tuple(sorted(self.data.items()))
+                        < tuple(sorted(other.data.items())))
+            if self.op is Op.APPLY:
+                if self.data[:2] != other.data[:2]:
+                    return self.data[:2] < other.data[:2]
+                return tuple(sorted(self.data[2].items())) < tuple(
+                    sorted(other.data[2].items()))
+            return self.data < other.data
+        return NotImplemented
+
+    def __repr__(self):
+        return f'{type(self).__name__}({self.op}, {self.data!r})'
+
+    def __str__(self):
+        return self.tostring()
+
+    def tostring(self, parent_precedence=Precedence.NONE,
+                 language=Language.Fortran):
+        """Return a string representation of Expr.
+        """
+        if self.op in (Op.INTEGER, Op.REAL):
+            precedence = (Precedence.SUM if self.data[0] < 0
+                          else Precedence.ATOM)
+            r = str(self.data[0]) + (f'_{self.data[1]}'
+                                     if self.data[1] != 4 else '')
+        elif self.op is Op.COMPLEX:
+            r = ', '.join(item.tostring(Precedence.TUPLE, language=language)
+                          for item in self.data)
+            r = '(' + r + ')'
+            precedence = Precedence.ATOM
+        elif self.op is Op.SYMBOL:
+            precedence = Precedence.ATOM
+            r = str(self.data)
+        elif self.op is Op.STRING:
+            r = self.data[0]
+            if self.data[1] != 1:
+                r = self.data[1] + '_' + r
+            precedence = Precedence.ATOM
+        elif self.op is Op.ARRAY:
+            r = ', '.join(item.tostring(Precedence.TUPLE, language=language)
+                          for item in self.data)
+            r = '[' + r + ']'
+            precedence = Precedence.ATOM
+        elif self.op is Op.TERMS:
+            terms = []
+            for term, coeff in sorted(self.data.items()):
+                if coeff < 0:
+                    op = ' - '
+                    coeff = -coeff
+                else:
+                    op = ' + '
+                if coeff == 1:
+                    term = term.tostring(Precedence.SUM, language=language)
+                else:
+                    if term == as_number(1):
+                        term = str(coeff)
+                    else:
+                        term = f'{coeff} * ' + term.tostring(
+                            Precedence.PRODUCT, language=language)
+                if terms:
+                    terms.append(op)
+                elif op == ' - ':
+                    terms.append('-')
+                terms.append(term)
+            r = ''.join(terms) or '0'
+            precedence = Precedence.SUM if terms else Precedence.ATOM
+        elif self.op is Op.FACTORS:
+            factors = []
+            tail = []
+            for base, exp in sorted(self.data.items()):
+                op = ' * '
+                if exp == 1:
+                    factor = base.tostring(Precedence.PRODUCT,
+                                           language=language)
+                elif language is Language.C:
+                    if exp in range(2, 10):
+                        factor = base.tostring(Precedence.PRODUCT,
+                                               language=language)
+                        factor = ' * '.join([factor] * exp)
+                    elif exp in range(-10, 0):
+                        factor = base.tostring(Precedence.PRODUCT,
+                                               language=language)
+                        tail += [factor] * -exp
+                        continue
+                    else:
+                        factor = base.tostring(Precedence.TUPLE,
+                                               language=language)
+                        factor = f'pow({factor}, {exp})'
+                else:
+                    factor = base.tostring(Precedence.POWER,
+                                           language=language) + f' ** {exp}'
+                if factors:
+                    factors.append(op)
+                factors.append(factor)
+            if tail:
+                if not factors:
+                    factors += ['1']
+                factors += ['/', '(', ' * '.join(tail), ')']
+            r = ''.join(factors) or '1'
+            precedence = Precedence.PRODUCT if factors else Precedence.ATOM
+        elif self.op is Op.APPLY:
+            name, args, kwargs = self.data
+            if name is ArithOp.DIV and language is Language.C:
+                numer, denom = [arg.tostring(Precedence.PRODUCT,
+                                             language=language)
+                                for arg in args]
+                r = f'{numer} / {denom}'
+                precedence = Precedence.PRODUCT
+            else:
+                args = [arg.tostring(Precedence.TUPLE, language=language)
+                        for arg in args]
+                args += [k + '=' + v.tostring(Precedence.NONE)
+                         for k, v in kwargs.items()]
+                r = f'{name}({", ".join(args)})'
+                precedence = Precedence.ATOM
+        elif self.op is Op.INDEXING:
+            name = self.data[0]
+            args = [arg.tostring(Precedence.TUPLE, language=language)
+                    for arg in self.data[1:]]
+            r = f'{name}[{", ".join(args)}]'
+            precedence = Precedence.ATOM
+        elif self.op is Op.CONCAT:
+            args = [arg.tostring(Precedence.PRODUCT, language=language)
+                    for arg in self.data]
+            r = " // ".join(args)
+            precedence = Precedence.PRODUCT
+        else:
+            raise NotImplementedError(f'tostring for op {self.op}')
+        if parent_precedence.value > precedence.value:
+            # If parent precedence is higher than operand precedence,
+            # operand will be enclosed in parenthesis.
+            return '(' + r + ')'
+        return r
+
+    def __pos__(self):
+        return self
+
+    def __neg__(self):
+        return self * -1
+
+    def __add__(self, other):
+        other = as_expr(other)
+        if isinstance(other, Expr):
+            if self.op is other.op:
+                if self.op in (Op.INTEGER, Op.REAL):
+                    return as_number(
+                        self.data[0] + other.data[0],
+                        max(self.data[1], other.data[1]))
+                if self.op is Op.COMPLEX:
+                    r1, i1 = self.data
+                    r2, i2 = other.data
+                    return as_complex(r1 + r2, i1 + i2)
+                if self.op is Op.TERMS:
+                    r = Expr(self.op, dict(self.data))
+                    for k, v in other.data.items():
+                        _pairs_add(r.data, k, v)
+                    return normalize(r)
+            if self.op is Op.COMPLEX and other.op in (Op.INTEGER, Op.REAL):
+                return self + as_complex(other)
+            elif self.op is (Op.INTEGER, Op.REAL) and other.op is Op.COMPLEX:
+                return as_complex(self) + other
+            elif self.op is Op.REAL and other.op is Op.INTEGER:
+                return self + as_real(other, kind=self.data[1])
+            elif self.op is Op.INTEGER and other.op is Op.REAL:
+                return as_real(self, kind=other.data[1]) + other
+            return as_terms(self) + as_terms(other)
+        return NotImplemented
+
+    def __radd__(self, other):
+        if isinstance(other, number_types):
+            return as_number(other) + self
+        return NotImplemented
+
+    def __sub__(self, other):
+        return self + (-other)
+
+    def __rsub__(self, other):
+        if isinstance(other, number_types):
+            return as_number(other) - self
+        return NotImplemented
+
+    def __mul__(self, other):
+        other = as_expr(other)
+        if isinstance(other, Expr):
+            if self.op is other.op:
+                if self.op in (Op.INTEGER, Op.REAL):
+                    return as_number(self.data[0] * other.data[0],
+                                     max(self.data[1], other.data[1]))
+                elif self.op is Op.COMPLEX:
+                    r1, i1 = self.data
+                    r2, i2 = other.data
+                    return as_complex(r1 * r2 - i1 * i2, r1 * i2 + r2 * i1)
+
+                if self.op is Op.FACTORS:
+                    r = Expr(self.op, dict(self.data))
+                    for k, v in other.data.items():
+                        _pairs_add(r.data, k, v)
+                    return normalize(r)
+                elif self.op is Op.TERMS:
+                    r = Expr(self.op, {})
+                    for t1, c1 in self.data.items():
+                        for t2, c2 in other.data.items():
+                            _pairs_add(r.data, t1 * t2, c1 * c2)
+                    return normalize(r)
+
+            if self.op is Op.COMPLEX and other.op in (Op.INTEGER, Op.REAL):
+                return self * as_complex(other)
+            elif other.op is Op.COMPLEX and self.op in (Op.INTEGER, Op.REAL):
+                return as_complex(self) * other
+            elif self.op is Op.REAL and other.op is Op.INTEGER:
+                return self * as_real(other, kind=self.data[1])
+            elif self.op is Op.INTEGER and other.op is Op.REAL:
+                return as_real(self, kind=other.data[1]) * other
+
+            if self.op is Op.TERMS:
+                return self * as_terms(other)
+            elif other.op is Op.TERMS:
+                return as_terms(self) * other
+
+            return as_factors(self) * as_factors(other)
+        return NotImplemented
+
+    def __rmul__(self, other):
+        if isinstance(other, number_types):
+            return as_number(other) * self
+        return NotImplemented
+
+    def __pow__(self, other):
+        other = as_expr(other)
+        if isinstance(other, Expr):
+            if other.op is Op.INTEGER:
+                exponent = other.data[0]
+                # TODO: other kind not used
+                if exponent == 0:
+                    return as_number(1)
+                if exponent == 1:
+                    return self
+                if exponent > 0:
+                    if self.op is Op.FACTORS:
+                        r = Expr(self.op, {})
+                        for k, v in self.data.items():
+                            r.data[k] = v * exponent
+                        return normalize(r)
+                    return self * (self ** (exponent - 1))
+                elif exponent != -1:
+                    return (self ** (-exponent)) ** -1
+                return Expr(Op.FACTORS, {self: exponent})
+            return as_apply(ArithOp.POW, self, other)
+        return NotImplemented
+
+    def __truediv__(self, other):
+        other = as_expr(other)
+        if isinstance(other, Expr):
+            # Fortran / is different from Python /:
+            # - `/` is a truncate operation for integer operands
+            return normalize(as_apply(ArithOp.DIV, self, other))
+        return NotImplemented
+
+    def __rtruediv__(self, other):
+        other = as_expr(other)
+        if isinstance(other, Expr):
+            return other / self
+        return NotImplemented
+
+    def __floordiv__(self, other):
+        other = as_expr(other)
+        if isinstance(other, Expr):
+            # Fortran // is different from Python //:
+            # - `//` is a concatenate operation for string operands
+            return normalize(Expr(Op.CONCAT, (self, other)))
+        return NotImplemented
+
+    def __rfloordiv__(self, other):
+        other = as_expr(other)
+        if isinstance(other, Expr):
+            return other // self
+        return NotImplemented
+
+    def __call__(self, *args, **kwargs):
+        # In Fortran, parenthesis () are use for both function call as
+        # well as indexing operations.
+        #
+        # TODO: implement a method for deciding when __call__ should
+        # return an INDEXING expression.
+        return as_apply(self, *map(as_expr, args),
+                        **dict((k, as_expr(v)) for k, v in kwargs.items()))
+
+    def __getitem__(self, index):
+        # Provided to support C indexing operations that .pyf files
+        # may contain.
+        index = as_expr(index)
+        if not isinstance(index, tuple):
+            index = index,
+        if len(index) > 1:
+            warnings.warn(
+                f'C-index should be a single expression but got `{index}`')
+        return Expr(Op.INDEXING, (self,) + index)
+
+    def substitute(self, symbols_map):
+        """Recursively substitute symbols with values in symbols map.
+
+        Symbols map is a dictionary of symbol-expression pairs.
+        """
+        if self.op is Op.SYMBOL:
+            value = symbols_map.get(self)
+            if value is None:
+                return self
+            m = re.match(r'\A(@__f2py_PARENTHESIS_(\w+)_\d+@)\Z', self.data)
+            if m:
+                # complement to fromstring method
+                items, paren = m.groups()
+                if paren in ['ROUNDDIV', 'SQUARE']:
+                    return as_array(value)
+                assert paren == 'ROUND', (paren, value)
+            return value
+        if self.op in (Op.INTEGER, Op.REAL, Op.STRING):
+            return self
+        if self.op in (Op.ARRAY, Op.COMPLEX):
+            return Expr(self.op, tuple(item.substitute(symbols_map)
+                                       for item in self.data))
+        if self.op is Op.CONCAT:
+            return normalize(Expr(self.op, tuple(item.substitute(symbols_map)
+                                                 for item in self.data)))
+        if self.op is Op.TERMS:
+            r = None
+            for term, coeff in self.data.items():
+                if r is None:
+                    r = term.substitute(symbols_map) * coeff
+                else:
+                    r += term.substitute(symbols_map) * coeff
+            if r is None:
+                warnings.warn('substitute: empty TERMS expression interpreted'
+                              ' as int-literal 0')
+                return as_number(0)
+            return r
+        if self.op is Op.FACTORS:
+            r = None
+            for base, exponent in self.data.items():
+                if r is None:
+                    r = base.substitute(symbols_map) ** exponent
+                else:
+                    r *= base.substitute(symbols_map) ** exponent
+            if r is None:
+                warnings.warn(
+                    'substitute: empty FACTORS expression interpreted'
+                    ' as int-literal 1')
+                return as_number(1)
+            return r
+        if self.op is Op.APPLY:
+            target, args, kwargs = self.data
+            if isinstance(target, Expr):
+                target = target.substitute(symbols_map)
+            args = tuple(a.substitute(symbols_map) for a in args)
+            kwargs = dict((k, v.substitute(symbols_map))
+                          for k, v in kwargs.items())
+            return normalize(Expr(self.op, (target, args, kwargs)))
+        if self.op is Op.INDEXING:
+            func = self.data[0]
+            if isinstance(func, Expr):
+                func = func.substitute(symbols_map)
+            args = tuple(a.substitute(symbols_map) for a in self.data[1:])
+            return normalize(Expr(self.op, (func,) + args))
+        raise NotImplementedError(f'substitute method for {self.op}: {self!r}')
+
+    def traverse(self, visit, *args, **kwargs):
+        """Traverse expression tree with visit function.
+
+        The visit function is applied to an expression with given args
+        and kwargs.
+
+        Traverse call returns an expression returned by visit when not
+        None, otherwise return a new normalized expression with
+        traverse-visit sub-expressions.
+        """
+        result = visit(self, *args, **kwargs)
+        if result is not None:
+            return result
+
+        if self.op in (Op.INTEGER, Op.REAL, Op.STRING, Op.SYMBOL):
+            return self
+        elif self.op in (Op.COMPLEX, Op.ARRAY, Op.CONCAT):
+            return normalize(Expr(self.op, tuple(
+                item.traverse(visit, *args, **kwargs)
+                for item in self.data)))
+        elif self.op in (Op.TERMS, Op.FACTORS):
+            data = {}
+            for k, v in self.data.items():
+                k = k.traverse(visit, *args, **kwargs)
+                v = (v.traverse(visit, *args, **kwargs)
+                     if isinstance(v, Expr) else v)
+                if k in data:
+                    v = data[k] + v
+                data[k] = v
+            return normalize(Expr(self.op, data))
+        elif self.op is Op.APPLY:
+            obj = self.data[0]
+            func = (obj.traverse(visit, *args, **kwargs)
+                    if isinstance(obj, Expr) else obj)
+            operands = tuple(operand.traverse(visit, *args, **kwargs)
+                             for operand in self.data[1])
+            kwoperands = dict((k, v.traverse(visit, *args, **kwargs))
+                              for k, v in self.data[2].items())
+            return normalize(Expr(self.op, (func, operands, kwoperands)))
+        elif self.op is Op.INDEXING:
+            obj = self.data[0]
+            obj = (obj.traverse(visit, *args, **kwargs)
+                   if isinstance(obj, Expr) else obj)
+            indices = tuple(index.traverse(visit, *args, **kwargs)
+                            for index in self.data[1:])
+            return normalize(Expr(self.op, (obj,) + indices))
+        raise NotImplementedError(f'traverse method for {self.op}')
+
+    def contains(self, other):
+        """Check if self contains other.
+        """
+        found = []
+
+        def visit(expr, found=found):
+            if found:
+                return expr
+            elif expr == other:
+                found.append(1)
+                return expr
+
+        self.traverse(visit)
+
+        return len(found) != 0
+
+    def symbols(self):
+        """Return a set of symbols contained in self.
+        """
+        found = set()
+
+        def visit(expr, found=found):
+            if expr.op is Op.SYMBOL:
+                found.add(expr)
+
+        self.traverse(visit)
+
+        return found
+
+    def polynomial_atoms(self):
+        """Return a set of expressions used as atoms in polynomial self.
+        """
+        found = set()
+
+        def visit(expr, found=found):
+            if expr.op is Op.FACTORS:
+                for b in expr.data:
+                    b.traverse(visit)
+                return expr
+            if expr.op in (Op.TERMS, Op.COMPLEX):
+                return
+            if expr.op is Op.APPLY and isinstance(expr.data[0], ArithOp):
+                if expr.data[0] is ArithOp.POW:
+                    expr.data[1][0].traverse(visit)
+                    return expr
+                return
+            if expr.op in (Op.INTEGER, Op.REAL):
+                return expr
+
+            found.add(expr)
+
+            if expr.op in (Op.INDEXING, Op.APPLY):
+                return expr
+
+        self.traverse(visit)
+
+        return found
+
+    def linear_solve(self, symbol):
+        """Return a, b such that a * symbol + b == self.
+
+        If self is not linear with respect to symbol, raise RuntimeError.
+        """
+        b = self.substitute({symbol: as_number(0)})
+        ax = self - b
+        a = ax.substitute({symbol: as_number(1)})
+
+        zero, _ = as_numer_denom(a * symbol - ax)
+
+        if zero != as_number(0):
+            raise RuntimeError(f'not a {symbol}-linear equation:'
+                               f' {a} * {symbol} + {b} == {self}')
+        return a, b
+
+
+def normalize(obj):
+    """Normalize Expr and apply basic evaluation methods.
+    """
+    if not isinstance(obj, Expr):
+        return obj
+
+    if obj.op is Op.TERMS:
+        d = {}
+        for t, c in obj.data.items():
+            if c == 0:
+                continue
+            if t.op is Op.COMPLEX and c != 1:
+                t = t * c
+                c = 1
+            if t.op is Op.TERMS:
+                for t1, c1 in t.data.items():
+                    _pairs_add(d, t1, c1 * c)
+            else:
+                _pairs_add(d, t, c)
+        if len(d) == 0:
+            # TODO: deterimine correct kind
+            return as_number(0)
+        elif len(d) == 1:
+            (t, c), = d.items()
+            if c == 1:
+                return t
+        return Expr(Op.TERMS, d)
+
+    if obj.op is Op.FACTORS:
+        coeff = 1
+        d = {}
+        for b, e in obj.data.items():
+            if e == 0:
+                continue
+            if b.op is Op.TERMS and isinstance(e, integer_types) and e > 1:
+                # expand integer powers of sums
+                b = b * (b ** (e - 1))
+                e = 1
+
+            if b.op in (Op.INTEGER, Op.REAL):
+                if e == 1:
+                    coeff *= b.data[0]
+                elif e > 0:
+                    coeff *= b.data[0] ** e
+                else:
+                    _pairs_add(d, b, e)
+            elif b.op is Op.FACTORS:
+                if e > 0 and isinstance(e, integer_types):
+                    for b1, e1 in b.data.items():
+                        _pairs_add(d, b1, e1 * e)
+                else:
+                    _pairs_add(d, b, e)
+            else:
+                _pairs_add(d, b, e)
+        if len(d) == 0 or coeff == 0:
+            # TODO: deterimine correct kind
+            assert isinstance(coeff, number_types)
+            return as_number(coeff)
+        elif len(d) == 1:
+            (b, e), = d.items()
+            if e == 1:
+                t = b
+            else:
+                t = Expr(Op.FACTORS, d)
+            if coeff == 1:
+                return t
+            return Expr(Op.TERMS, {t: coeff})
+        elif coeff == 1:
+            return Expr(Op.FACTORS, d)
+        else:
+            return Expr(Op.TERMS, {Expr(Op.FACTORS, d): coeff})
+
+    if obj.op is Op.APPLY and obj.data[0] is ArithOp.DIV:
+        dividend, divisor = obj.data[1]
+        t1, c1 = as_term_coeff(dividend)
+        t2, c2 = as_term_coeff(divisor)
+        if isinstance(c1, integer_types) and isinstance(c2, integer_types):
+            g = gcd(c1, c2)
+            c1, c2 = c1//g, c2//g
+        else:
+            c1, c2 = c1/c2, 1
+
+        if t1.op is Op.APPLY and t1.data[0] is ArithOp.DIV:
+            numer = t1.data[1][0] * c1
+            denom = t1.data[1][1] * t2 * c2
+            return as_apply(ArithOp.DIV, numer, denom)
+
+        if t2.op is Op.APPLY and t2.data[0] is ArithOp.DIV:
+            numer = t2.data[1][1] * t1 * c1
+            denom = t2.data[1][0] * c2
+            return as_apply(ArithOp.DIV, numer, denom)
+
+        d = dict(as_factors(t1).data)
+        for b, e in as_factors(t2).data.items():
+            _pairs_add(d, b, -e)
+        numer, denom = {}, {}
+        for b, e in d.items():
+            if e > 0:
+                numer[b] = e
+            else:
+                denom[b] = -e
+        numer = normalize(Expr(Op.FACTORS, numer)) * c1
+        denom = normalize(Expr(Op.FACTORS, denom)) * c2
+
+        if denom.op in (Op.INTEGER, Op.REAL) and denom.data[0] == 1:
+            # TODO: denom kind not used
+            return numer
+        return as_apply(ArithOp.DIV, numer, denom)
+
+    if obj.op is Op.CONCAT:
+        lst = [obj.data[0]]
+        for s in obj.data[1:]:
+            last = lst[-1]
+            if (
+                    last.op is Op.STRING
+                    and s.op is Op.STRING
+                    and last.data[0][0] in '"\''
+                    and s.data[0][0] == last.data[0][-1]
+            ):
+                new_last = as_string(last.data[0][:-1] + s.data[0][1:],
+                                     max(last.data[1], s.data[1]))
+                lst[-1] = new_last
+            else:
+                lst.append(s)
+        if len(lst) == 1:
+            return lst[0]
+        return Expr(Op.CONCAT, tuple(lst))
+    return obj
+
+
+def as_expr(obj):
+    """Convert non-Expr objects to Expr objects.
+    """
+    if isinstance(obj, complex):
+        return as_complex(obj.real, obj.imag)
+    if isinstance(obj, number_types):
+        return as_number(obj)
+    if isinstance(obj, str):
+        # STRING expression holds string with boundary quotes, hence
+        # applying repr:
+        return as_string(repr(obj))
+    if isinstance(obj, tuple):
+        return tuple(map(as_expr, obj))
+    return obj
+
+
+def as_symbol(obj):
+    """Return object as SYMBOL expression (variable or unparsed expression).
+    """
+    return Expr(Op.SYMBOL, obj)
+
+
+def as_number(obj, kind=4):
+    """Return object as INTEGER or REAL constant.
+    """
+    if isinstance(obj, int):
+        return Expr(Op.INTEGER, (obj, kind))
+    if isinstance(obj, float):
+        return Expr(Op.REAL, (obj, kind))
+    if isinstance(obj, Expr):
+        if obj.op in (Op.INTEGER, Op.REAL):
+            return obj
+    raise OpError(f'cannot convert {obj} to INTEGER or REAL constant')
+
+
+def as_integer(obj, kind=4):
+    """Return object as INTEGER constant.
+    """
+    if isinstance(obj, int):
+        return Expr(Op.INTEGER, (obj, kind))
+    if isinstance(obj, Expr):
+        if obj.op is Op.INTEGER:
+            return obj
+    raise OpError(f'cannot convert {obj} to INTEGER constant')
+
+
+def as_real(obj, kind=4):
+    """Return object as REAL constant.
+    """
+    if isinstance(obj, int):
+        return Expr(Op.REAL, (float(obj), kind))
+    if isinstance(obj, float):
+        return Expr(Op.REAL, (obj, kind))
+    if isinstance(obj, Expr):
+        if obj.op is Op.REAL:
+            return obj
+        elif obj.op is Op.INTEGER:
+            return Expr(Op.REAL, (float(obj.data[0]), kind))
+    raise OpError(f'cannot convert {obj} to REAL constant')
+
+
+def as_string(obj, kind=1):
+    """Return object as STRING expression (string literal constant).
+    """
+    return Expr(Op.STRING, (obj, kind))
+
+
+def as_array(obj):
+    """Return object as ARRAY expression (array constant).
+    """
+    if isinstance(obj, Expr):
+        obj = obj,
+    return Expr(Op.ARRAY, obj)
+
+
+def as_complex(real, imag=0):
+    """Return object as COMPLEX expression (complex literal constant).
+    """
+    return Expr(Op.COMPLEX, (as_expr(real), as_expr(imag)))
+
+
+def as_apply(func, *args, **kwargs):
+    """Return object as APPLY expression (function call, constructor, etc.)
+    """
+    return Expr(Op.APPLY,
+                (func, tuple(map(as_expr, args)),
+                 dict((k, as_expr(v)) for k, v in kwargs.items())))
+
+
+def as_terms(obj):
+    """Return expression as TERMS expression.
+    """
+    if isinstance(obj, Expr):
+        obj = normalize(obj)
+        if obj.op is Op.TERMS:
+            return obj
+        if obj.op is Op.INTEGER:
+            return Expr(Op.TERMS, {as_integer(1, obj.data[1]): obj.data[0]})
+        if obj.op is Op.REAL:
+            return Expr(Op.TERMS, {as_real(1, obj.data[1]): obj.data[0]})
+        return Expr(Op.TERMS, {obj: 1})
+    raise OpError(f'cannot convert {type(obj)} to terms Expr')
+
+
+def as_factors(obj):
+    """Return expression as FACTORS expression.
+    """
+    if isinstance(obj, Expr):
+        obj = normalize(obj)
+        if obj.op is Op.FACTORS:
+            return obj
+        if obj.op is Op.TERMS:
+            if len(obj.data) == 1:
+                (term, coeff), = obj.data.items()
+                if coeff == 1:
+                    return Expr(Op.FACTORS, {term: 1})
+                return Expr(Op.FACTORS, {term: 1, Expr.number(coeff): 1})
+        if ((obj.op is Op.APPLY
+             and obj.data[0] is ArithOp.DIV
+             and not obj.data[2])):
+            return Expr(Op.FACTORS, {obj.data[1][0]: 1, obj.data[1][1]: -1})
+        return Expr(Op.FACTORS, {obj: 1})
+    raise OpError(f'cannot convert {type(obj)} to terms Expr')
+
+
+def as_term_coeff(obj):
+    """Return expression as term-coefficient pair.
+    """
+    if isinstance(obj, Expr):
+        obj = normalize(obj)
+        if obj.op is Op.INTEGER:
+            return as_integer(1, obj.data[1]), obj.data[0]
+        if obj.op is Op.REAL:
+            return as_real(1, obj.data[1]), obj.data[0]
+        if obj.op is Op.TERMS:
+            if len(obj.data) == 1:
+                (term, coeff), = obj.data.items()
+                return term, coeff
+            # TODO: find common divisior of coefficients
+        if obj.op is Op.APPLY and obj.data[0] is ArithOp.DIV:
+            t, c = as_term_coeff(obj.data[1][0])
+            return as_apply(ArithOp.DIV, t, obj.data[1][1]), c
+        return obj, 1
+    raise OpError(f'cannot convert {type(obj)} to term and coeff')
+
+
+def as_numer_denom(obj):
+    """Return expression as numer-denom pair.
+    """
+    if isinstance(obj, Expr):
+        obj = normalize(obj)
+        if obj.op in (Op.INTEGER, Op.REAL, Op.COMPLEX, Op.SYMBOL, Op.INDEXING):
+            return obj, as_number(1)
+        elif obj.op is Op.APPLY:
+            if obj.data[0] is ArithOp.DIV and not obj.data[2]:
+                numers, denoms = map(as_numer_denom, obj.data[1])
+                return numers[0] * denoms[1], numers[1] * denoms[0]
+            return obj, as_number(1)
+        elif obj.op is Op.TERMS:
+            numers, denoms = [], []
+            for term, coeff in obj.data.items():
+                n, d = as_numer_denom(term)
+                n = n * coeff
+                numers.append(n)
+                denoms.append(d)
+            numer, denom = as_number(0), as_number(1)
+            for i in range(len(numers)):
+                n = numers[i]
+                for j in range(len(numers)):
+                    if i != j:
+                        n *= denoms[j]
+                numer += n
+                denom *= denoms[i]
+            if denom.op in (Op.INTEGER, Op.REAL) and denom.data[0] < 0:
+                numer, denom = -numer, -denom
+            return numer, denom
+        elif obj.op is Op.FACTORS:
+            numer, denom = as_number(1), as_number(1)
+            for b, e in obj.data.items():
+                bnumer, bdenom = as_numer_denom(b)
+                if e > 0:
+                    numer *= bnumer ** e
+                    denom *= bdenom ** e
+                elif e < 0:
+                    numer *= bdenom ** (-e)
+                    denom *= bnumer ** (-e)
+            return numer, denom
+    raise OpError(f'cannot convert {type(obj)} to numer and denom')
+
+
+def _counter():
+    # Used internally to generate unique dummy symbols
+    counter = 0
+    while True:
+        counter += 1
+        yield counter
+
+
+COUNTER = _counter()
+
+
+def replace_quotes(s):
+    """Replace quoted substrings of input.
+
+    Return a new string and a mapping of replacements.
+    """
+    d = {}
+
+    def repl(m):
+        kind, value = m.groups()[:2]
+        if kind:
+            # remove trailing underscore
+            kind = kind[:-1]
+        p = {"'": "SINGLE", '"': "DOUBLE"}[value[0]]
+        k = f'@__f2py_QUOTES_{p}_{COUNTER.__next__()}@'
+        d[as_symbol(k)] = as_string(value, kind or 1)
+        return k
+
+    return (re.sub(r'({kind}_|)({single_quoted}|{double_quoted})'.format(
+        kind=r'\w[\w\d_]*',
+        single_quoted=r"('([^'\\]|(\\.))*')",
+        double_quoted=r'("([^"\\]|(\\.))*")'),
+        repl, s), d)
+
+
+def replace_parenthesis(s):
+    """Replace substrings of input that are enclosed in parenthesis.
+
+    Return a new string and a mapping of replacements.
+    """
+    # Find a parenthesis pair that appears first.
+
+    # Fortran deliminator are `(`, `)`, `[`, `]`, `(/', '/)`, `/`.
+    # We don't handle `/` deliminator because it is not a part of an
+    # expression.
+    left, right = None, None
+    mn_i = len(s)
+    for left_, right_ in (('(/', '/)'),
+                          '()',
+                          '{}',  # to support C literal structs
+                          '[]'):
+        i = s.find(left_)
+        if i == -1:
+            continue
+        if i < mn_i:
+            mn_i = i
+            left, right = left_, right_
+
+    if left is None:
+        return s, {}
+
+    i = mn_i
+    j = s.find(right, i)
+
+    while s.count(left, i + 1, j) != s.count(right, i + 1, j):
+        j = s.find(right, j + 1)
+        if j == -1:
+            raise ValueError(f'Mismatch of {left+right} parenthesis in {s!r}')
+
+    p = {'(': 'ROUND', '[': 'SQUARE', '{': 'CURLY', '(/': 'ROUNDDIV'}[left]
+
+    k = f'@__f2py_PARENTHESIS_{p}_{COUNTER.__next__()}@'
+    v = s[i+len(left):j]
+    r, d = replace_parenthesis(s[j+len(right):])
+    d[k] = v
+    return s[:i] + k + r, d
+
+
+def fromstring(s):
+    """Create an expression from string.
+
+    This is a "lazy" parser, that is, only arithmetic operations are
+    resolved, non-arithmetic operations are treated as symbols.
+    """
+    # We replace string literal constants that may include parenthesis
+    # that may confuse the arithmetic expression parser.
+    s, quotes_map = replace_quotes(s)
+
+    # Apply arithmetic expression parser.
+    r = _fromstring_worker(s)
+
+    # Restore string constants.
+    return r.substitute(quotes_map)
+
+
+class _Pair:
+    # Internal class to represent a pair of expressions
+
+    def __init__(self, left, right):
+        self.left = left
+        self.right = right
+
+    def substitute(self, symbols_map):
+        left, right = self.left, self.right
+        if isinstance(left, Expr):
+            left = left.substitute(symbols_map)
+        if isinstance(right, Expr):
+            right = right.substitute(symbols_map)
+        return _Pair(left, right)
+
+    def __repr__(self):
+        return f'{type(self).__name__}({self.left}, {self.right})'
+
+
+def _fromstring_worker(s, dummy=None):
+    # Internal method used by fromstring to convert a string of
+    # arithmetic expressions to an expression tree.
+
+    # Replace subexpression enclosed in parenthesis with dummy symbols
+    # and recursively parse subexpression into symbol-expression map:
+    r, raw_symbols_map = replace_parenthesis(s)
+    r = r.strip()
+
+    symbols_map = dict([(as_symbol(k), _fromstring_worker(v))
+                        for k, v in raw_symbols_map.items()])
+
+    if ',' in r:
+        # comma-separate tuple of expressions
+        return tuple(_fromstring_worker(r_).substitute(symbols_map)
+                     for r_ in r.split(','))
+
+    m = re.match(r'\A(\w[\w\d_]*)\s*[=](.*)\Z', r)
+    if m:
+        # keyword argument
+        keyname, value = m.groups()
+        return _Pair(keyname,
+                     _fromstring_worker(value).substitute(symbols_map))
+
+    operands = re.split(r'((?<!\d[edED])[+-])', r)
+    if len(operands) > 1:
+        # Expression is an arithmetic sum
+        result = _fromstring_worker(operands[0] or '0')
+        for op, operand in zip(operands[1::2], operands[2::2]):
+            op = op.strip()
+            if op == '+':
+                result += _fromstring_worker(operand)
+            else:
+                assert op == '-'
+                result -= _fromstring_worker(operand)
+        return result.substitute(symbols_map)
+
+    operands = re.split(r'//', r)
+    if len(operands) > 1:
+        # Expression is string concatenate operation
+        return Expr(Op.CONCAT,
+                    tuple(_fromstring_worker(operand).substitute(symbols_map)
+                          for operand in operands))
+
+    operands = re.split(r'([*]|/)', r.replace('**', '@__f2py_DOUBLE_STAR@'))
+    operands = [operand.replace('@__f2py_DOUBLE_STAR@', '**')
+                for operand in operands]
+    if len(operands) > 1:
+        # Expression is an arithmetic product
+        result = _fromstring_worker(operands[0])
+        for op, operand in zip(operands[1::2], operands[2::2]):
+            op = op.strip()
+            if op == '*':
+                result *= _fromstring_worker(operand)
+            else:
+                assert op == '/'
+                result /= _fromstring_worker(operand)
+        return result.substitute(symbols_map)
+
+    operands = list(reversed(re.split(r'[*][*]', r)))
+    if len(operands) > 1:
+        # Expression is an arithmetic exponentiation
+        result = _fromstring_worker(operands[0])
+        for operand in operands[1:]:
+            result = _fromstring_worker(operand) ** result
+        return result
+
+    m = re.match(r'\A({digit_string})({kind}|)\Z'.format(
+        digit_string=r'\d+',
+        kind=r'_(\d+|\w[\w\d_]*)'), r)
+    if m:
+        # Expression is int-literal-constant
+        value, _, kind = m.groups()
+        if kind and kind.isdigit():
+            kind = int(kind)
+        return as_integer(int(value), kind or 4)
+
+    m = re.match(r'\A({significant}({exponent}|)|\d+{exponent})({kind}|)\Z'
+                 .format(
+                     significant=r'[.]\d+|\d+[.]\d*',
+                     exponent=r'[edED][+-]?\d+',
+                     kind=r'_(\d+|\w[\w\d_]*)'), r)
+    if m:
+        # Expression is real-literal-constant
+        value, _, _, kind = m.groups()
+        if kind and kind.isdigit():
+            kind = int(kind)
+        value = value.lower()
+        if 'd' in value:
+            return as_real(float(value.replace('d', 'e')), kind or 8)
+        return as_real(float(value), kind or 4)
+
+    m = re.match(r'\A(\w[\w\d_]*)_(@__f2py_QUOTES_(\w+)_\d+@)\Z', r)
+    if m:
+        # string-literal-constant with kind parameter specification
+        kind, name, _ = m.groups()
+        return as_string(name, kind)
+
+    m = re.match(r'\A(\w[\w\d+]*)\s*(@__f2py_PARENTHESIS_(\w+)_\d+@)\Z', r)
+    if m:
+        # apply or indexing expression
+        name, args, paren = m.groups()
+        target = as_symbol(name).substitute(symbols_map)
+        if args in raw_symbols_map:
+            args = symbols_map[as_symbol(args)]
+            if not isinstance(args, tuple):
+                args = args,
+            if paren == 'ROUND':
+                # Expression is a function call or a Fortran indexing
+                # operation
+                kwargs = dict((a.left, a.right) for a in args
+                              if isinstance(a, _Pair))
+                args = tuple(a for a in args if not isinstance(a, _Pair))
+                return as_apply(target, *args, **kwargs)
+            if paren == 'SQUARE':
+                # Expression is a C indexing operation (e.g. used in
+                # .pyf files)
+                return target[args]
+        # A non-Fortran/C expression?
+        assert 0, (s, r)
+
+    m = re.match(r'\A(@__f2py_PARENTHESIS_(\w+)_\d+@)\Z', r)
+    if m:
+        # array constructor or literal complex constant
+        items, paren = m.groups()
+        if items in raw_symbols_map:
+            items = symbols_map[as_symbol(items)]
+            if (paren == 'ROUND' and isinstance(items,
+                                                tuple) and len(items) == 2):
+                # Expression is a literal complex constant
+                return as_complex(items[0], items[1])
+            if paren in ['ROUNDDIV', 'SQUARE']:
+                # Expression is a array constructor
+                return as_array(items)
+            assert 0, (r, items, paren)
+
+    m = re.match(r'\A\w[\w\d_]*\Z', r)
+    if m:
+        # Fortran standard conforming name
+        return as_symbol(r)
+
+    m = re.match(r'\A@\w[\w\d_]*@\Z', r)
+    if m:
+        # f2py special dummy name
+        return as_symbol(r)
+
+    warnings.warn(f'fromstring: treating {r!r} as symbol')
+    return as_symbol(r)