Constant initializers may call the intrinsics both sides fold identically

src/recast/fortran/constants.py

@@ -191,16 +191,25 @@ def _split_arguments(tokens: list[str]) -> list[list[str]]:
 
 def _argument_tokens(
     tokens: list[str], known_names: set[str], aliases: dict[str, str] | None = None
-) -> list[dict[str, str]] | None:
+) -> list[dict[str, Any]] | None:
     """One argument as tokens of the same vocabulary; ``None`` if it names
     something no earlier constant defines."""
-    spelled: list[dict[str, str]] = []
-    for piece in tokens:
+    spelled: list[dict[str, Any]] = []
+    at = 0
+    while at < len(tokens):
+        piece = tokens[at]
         if re.match(r"[A-Za-z_]", piece):
+            if piece.lower() in INTRINSICS and at + 1 < len(tokens) and tokens[at + 1] == "(":
+                # A call inside an argument: ``max( 1.e-10, epsilon(tol) )``.
+                call, at = _intrinsic_call(tokens, at, known_names, aliases)
+                if call is None:
+                    return None
+                spelled.append(call)
+                continue
             if piece.lower() in known_names:
                 spelled.append({"t": "ref", "v": _spelled_ref(piece, aliases)})
             elif _KIND_ARGUMENT.match(piece):
-                continue
+                pass
             else:
                 return None
         elif re.match(r"\d", piece):
@@ -211,9 +220,13 @@ def _argument_tokens(
                 spelled.append({"t": "int", "v": base})
         else:
             spelled.append({"t": "op", "v": piece})
+        at += 1
     return spelled
 
 
+_KIND_INQUIRIES = frozenset({"epsilon", "huge", "tiny"})
+
+
 def _spelled_ref(name: str, aliases: dict[str, str] | None) -> str:
     """The name a reference is recorded under.
 
@@ -233,7 +246,12 @@ def _intrinsic_call(
 
     A trailing kind argument -- ``real(x, r8)`` -- is dropped: it says what
     precision the compiler evaluated in, which the target's own float64 is,
-    and it is not a value to pass on.
+    and it is not a value to pass on. The argument of a kind inquiry --
+    ``epsilon(pi)``, ``huge(x)`` -- contributes its kind and no value, and
+    may legally be the constant being declared
+    (``tol = max( 1.e-10_core_rknd, epsilon(tol) )``, CLUBB); when it names
+    nothing yet defined the call is kept with no argument, which is what the
+    target renders anyway.
     """
     name = tokens[at].lower()
     depth = 0
@@ -256,6 +274,10 @@ def _intrinsic_call(
         kept.append(argument)
     spelled = [_argument_tokens(argument, known_names, aliases) for argument in kept]
     if any(text is None for text in spelled):
+        if name in _KIND_INQUIRIES:
+            # The argument names nothing yet defined -- the constant itself,
+            # legally -- and only its kind was ever asked for.
+            return {"t": "call", "v": name, "args": []}, end + 1
         return None, end + 1
     return {"t": "call", "v": name, "args": spelled}, end + 1
 

src/recast/fortran/expr.py

@@ -13,8 +13,9 @@
 to prevent.
 
 Deliberately small. These are physical-constant initializers -- sums, products
-and powers over literals and earlier constants. Anything richer raises
-``UnsupportedExpression`` rather than being approximated.
+and powers over literals and earlier constants, and a short list of intrinsic
+calls both languages fold to the same bits (``INTRINSICS``). Anything richer
+raises ``UnsupportedExpression`` rather than being approximated.
 """
 
 from __future__ import annotations
@@ -29,6 +30,28 @@
 BINARY_OPS = ("+", "-", "*", "/", "**")
 UNARY_OPS = ("+", "-")
 
+INTRINSICS = frozenset(
+    {"max", "min", "abs", "sqrt", "epsilon", "huge", "tiny", "real", "dble", "int"}
+)
+"""The intrinsic calls an initializer may make.
+
+Each is one the compiler folds to the same value NumPy computes: ``max`` /
+``min`` / ``abs`` select or negate, ``sqrt`` is correctly rounded on both
+sides (IEEE 754 requires it), the three kind inquiries are the kind's own
+constants, and ``real`` / ``dble`` / ``int`` are conversions. ``exp``,
+``log`` and the trigonometric functions are not here: gfortran folds them
+with MPFR and libm need not agree in the last bit, so a constant made of one
+is declined rather than approximated.
+"""
+
+KIND_INQUIRIES = frozenset({"epsilon", "huge", "tiny"})
+"""Calls whose argument contributes its kind, not its value."""
+
+KINDS_64 = frozenset({"8", "dp", "r8", "core_rknd", "shr_kind_r8", "real64"})
+"""Spellings of a 64-bit real kind a ``kind=`` argument may name. Any other
+spelling refuses: the fold renders every real as 64-bit and must not claim a
+conversion it cannot honour."""
+
 
 class UnsupportedExpression(RecastError):
     """An initializer this frontend will not claim to understand.
@@ -43,8 +66,9 @@ class UnsupportedExpression(RecastError):
 class Expr:
     """One node of a constant initializer.
 
-    ``kind`` is ``real``, ``int``, ``name``, ``paren``, ``unary`` or ``binary``.
-    ``text`` carries the literal text, the identifier, or the operator.
+    ``kind`` is ``real``, ``int``, ``name``, ``paren``, ``unary``, ``binary``
+    or ``call``. ``text`` carries the literal text, the identifier, the
+    operator, or the intrinsic's lower-case name.
     """
 
     kind: str
@@ -78,30 +102,97 @@ def build(node: Any) -> Expr:
     if children and len(children) == 3 and isinstance(children[1], str):
         if children[1] in BINARY_OPS:
             return Expr("binary", children[1], (build(children[0]), build(children[2])))
+    if isinstance(node, f03.Intrinsic_Function_Reference):
+        return _call(node)
     raise UnsupportedExpression(f"unsupported initializer node {type(node).__name__}: {node}")
 
 
+def _call(node: Any) -> Expr:
+    """An intrinsic call from ``INTRINSICS``, its arguments built in order.
+
+    A ``kind=`` argument on a conversion is checked and dropped: the fold
+    renders every real as 64-bit, so a 64-bit kind is the identity and any
+    other kind is refused. Every other keyword argument refuses -- the
+    intrinsics listed take positional arguments in every initializer seen.
+    """
+    fname = str(node.children[0]).lower()
+    if fname not in INTRINSICS:
+        raise UnsupportedExpression(f"unsupported intrinsic in initializer: {node}")
+    spec = node.children[1]
+    args: list[Expr] = []
+    for item in spec.items if spec is not None else ():
+        if isinstance(item, f03.Actual_Arg_Spec):
+            keyword, value = (str(c).lower() for c in item.children)
+            if keyword == "kind" and fname in {"real", "dble", "int"} and value in KINDS_64:
+                continue
+            raise UnsupportedExpression(f"unsupported keyword argument in initializer: {node}")
+        args.append(build(item))
+    if fname in {"real", "dble"} and len(args) == 2:
+        # ``real(x, r8)``: the positional form of the same kind argument.
+        if args[1].kind == "name" and args[1].text in KINDS_64:
+            args = args[:1]
+        elif args[1].kind == "int" and args[1].text in KINDS_64:
+            args = args[:1]
+        else:
+            raise UnsupportedExpression(f"unsupported kind in initializer: {node}")
+    arity = {
+        "abs": 1,
+        "sqrt": 1,
+        "epsilon": 1,
+        "huge": 1,
+        "tiny": 1,
+        "real": 1,
+        "dble": 1,
+        "int": 1,
+    }
+    if len(args) != arity.get(fname, len(args)) or (fname in {"max", "min"} and len(args) < 2):
+        raise UnsupportedExpression(f"unsupported argument count in initializer: {node}")
+    return Expr("call", fname, tuple(args))
+
+
+def substitute(expr: Expr, name: str, replacement: Expr) -> Expr:
+    """``expr`` with every reference to ``name`` replaced.
+
+    For the one legal self-reference in an initializer, a kind inquiry on the
+    constant being declared (``tol = max( 1.e-10_r8, epsilon(tol) )``): the
+    reference carries the constant's kind and nothing else, and the fold
+    renders reals as 64-bit, so a 64-bit literal stands in for it.
+    """
+    if expr.kind == "name" and expr.text == name:
+        return replacement
+    if not expr.args:
+        return expr
+    return Expr(expr.kind, expr.text, tuple(substitute(a, name, replacement) for a in expr.args))
+
+
 def render(
     expr: Expr,
     *,
     real: Callable[[str], str],
     integer: Callable[[str], str],
     name: Callable[[str], str],
+    call: Callable[[str, list[str]], str] | None = None,
 ) -> str:
-    """Fold an ``Expr`` to text, given how to spell its three kinds of atom.
+    """Fold an ``Expr`` to text, given how to spell its three kinds of atom
+    and, optionally, an intrinsic call over already-rendered arguments.
 
     Grouping and spacing are fixed here so that every target language brackets
     the arithmetic identically. That is the whole point: two renderings of one
     tree can differ in how a literal is spelled and not in what is multiplied
-    by what.
+    by what. A renderer given no ``call`` refuses a tree with one in it
+    rather than guessing a spelling.
     """
     if expr.kind == "real":
         return real(expr.text)
     if expr.kind == "int":
         return integer(expr.text)
     if expr.kind == "name":
         return name(expr.text)
-    sub = [render(a, real=real, integer=integer, name=name) for a in expr.args]
+    sub = [render(a, real=real, integer=integer, name=name, call=call) for a in expr.args]
+    if expr.kind == "call":
+        if call is None:
+            raise UnsupportedExpression(f"no rendering for intrinsic {expr.text!r} in this target")
+        return call(expr.text, sub)
     if expr.kind == "paren":
         return f"({sub[0]})"
     if expr.kind == "unary":
@@ -111,6 +202,60 @@ def render(
     raise UnsupportedExpression(f"unknown Expr kind {expr.kind!r}")
 
 
+REAL_CALLS = frozenset({"real", "dble", "sqrt"} | KIND_INQUIRIES)
+
+
+def typed(expr: Expr) -> str | None:
+    """``"real"``, ``"int"``, or ``None`` when a bare name leaves it open.
+
+    Type inference the fold needs for exactly one decision: whether a ``/``
+    is Fortran's integer division. A real literal or a real-valued call
+    anywhere in an operand makes the quotient real; ``int(...)`` and integer
+    literals make it integer; a name is whatever its initializer was, which
+    this tree does not carry.
+    """
+    if expr.kind == "real":
+        return "real"
+    if expr.kind == "int":
+        return "int"
+    if expr.kind == "name":
+        return None
+    if expr.kind == "call":
+        if expr.text == "int":
+            return "int"
+        if expr.text in REAL_CALLS:
+            return "real"
+    kinds = {typed(a) for a in expr.args}
+    if "real" in kinds:
+        return "real"
+    if kinds == {"int"}:
+        return "int"
+    return None
+
+
+def with_integer_division(expr: Expr, *, default_integer: bool | None = None) -> Expr:
+    """The tree with every integer ``/`` spelled ``//``.
+
+    Fortran divides two integers to an integer: ``nrk = runge_kutta_type / 10``
+    is 4, not 4.1. A quotient whose operands are both known integers is
+    marked; one with a name in it falls back to ``default_integer``, which
+    the caller sets from the whole initializer -- an expression with no real
+    literal and no real-valued call in it is integer arithmetic throughout,
+    because a name in it is an integer parameter or it would have had one.
+    """
+    if default_integer is None:
+        default_integer = typed(expr) != "real"
+    if not expr.args:
+        return expr
+    args = tuple(with_integer_division(a, default_integer=default_integer) for a in expr.args)
+    text = expr.text
+    if expr.kind == "binary" and expr.text == "/":
+        kinds = {typed(a) for a in args}
+        if kinds == {"int"} or ("real" not in kinds and default_integer):
+            text = "//"
+    return Expr(expr.kind, text, args)
+
+
 def names_used(expr: Expr) -> list[str]:
     """Every identifier the expression depends on, in traversal order."""
     if expr.kind == "name":
@@ -119,3 +264,28 @@ def names_used(expr: Expr) -> list[str]:
     for a in expr.args:
         out.extend(names_used(a))
     return out
+
+
+def python_call(fname: str, args: list[str], *, real64: str = "np.float64") -> str:
+    """Spell one whitelisted intrinsic in Python over rendered arguments.
+
+    ``real64`` names the 64-bit real constructor the caller renders literals
+    with (``np.float64`` in an emitted module, ``float`` in an evaluator that
+    imports nothing). A kind inquiry is spelled over the argument's own value
+    where NumPy is available -- ``np.finfo(PI).eps`` is the epsilon of PI's
+    kind -- and over the 64-bit kind otherwise.
+    """
+    numpy = real64 == "np.float64"
+    if fname in {"max", "min", "abs", "int"}:
+        return f"{fname}({', '.join(args)})"
+    if fname == "sqrt":
+        return f"np.sqrt({args[0]})" if numpy else f"math.sqrt({args[0]})"
+    if fname in {"real", "dble"}:
+        return f"{real64}({args[0]})"
+    if fname in KIND_INQUIRIES:
+        attr = {"epsilon": "eps", "huge": "max", "tiny": "tiny"}[fname]
+        if numpy:
+            return f"np.finfo({args[0]}).{attr}"
+        field = {"epsilon": "epsilon", "huge": "max", "tiny": "min"}[fname]
+        return f"sys.float_info.{field}"
+    raise UnsupportedExpression(f"no Python spelling for intrinsic {fname!r}")

src/recast/fortran/tree.py

@@ -11,7 +11,9 @@
 
 from __future__ import annotations
 
+import math
 import re
+import sys
 from pathlib import Path
 from typing import Any
 
@@ -155,23 +157,28 @@ def _evaluate(
 ) -> Any:
     """The value of one named constant, or ``None`` where the tree does not
     initialize it with something a parameter can be folded from."""
+    # Lazy, like ``render`` above: ``expr`` parses, and this module is imported
+    # by paths that must stay importable without the ``fortran`` extra.
+    from recast.fortran.expr import python_call, with_integer_division
+
     try:
         records = resolve([name], files)
     except unresolved:
         return None
     env: dict[str, Any] = {}
     try:
         for entry in records:
+            # Integer arithmetic where Fortran's ``/`` truncates.
             text = render(
-                entry["expr"],
+                with_integer_division(entry["expr"]),
                 real=lambda t: f"float('{t}')",
                 integer=lambda t: t,
                 name=lambda t: t.upper(),
+                call=lambda f, a: python_call(f, a, real64="float"),
             )
-            if "float(" not in text:
-                # Integer arithmetic throughout: Fortran's ``/`` truncates.
-                text = text.replace("/", "//")
-            env[entry["name"].upper()] = eval(text, {"__builtins__": {}}, dict(env))  # noqa: S307
+            scope = {"__builtins__": {}, "max": max, "min": min, "abs": abs, "int": int}
+            scope.update({"float": float, "math": math, "sys": sys})
+            env[entry["name"].upper()] = eval(text, scope, dict(env))  # noqa: S307
     except Exception:  # an initializer shape the renderer has no rule for
         return None
     return env.get(name.upper())

src/recast/fortran/use.py

@@ -18,7 +18,7 @@
 
 from recast.errors import RecastError
 from recast.fortran._parse import f03, parse, walk
-from recast.fortran.expr import Expr, build, names_used
+from recast.fortran.expr import Expr, build, names_used, substitute
 
 
 class UnresolvedConstant(RecastError):
@@ -81,7 +81,9 @@ def need(name: str) -> None:
             raise UnresolvedConstant(f"no initializer for {name!r} in {[str(s) for s in sources]}")
         seen.add(name)
         node, line = table[name]
-        expr: Expr = build(node)
+        # The one legal self-reference, a kind inquiry on the constant being
+        # declared, stands for its kind alone; see ``expr.substitute``.
+        expr: Expr = substitute(build(node), name, Expr("real", "1.0"))
         for dep in names_used(expr):
             need(dep)
         ordered.append(