diff --git a/src/recast/fortran/constants.py b/src/recast/fortran/constants.py index 8c0803b..b8d2883 100644 --- a/src/recast/fortran/constants.py +++ b/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 diff --git a/src/recast/fortran/expr.py b/src/recast/fortran/expr.py index e08fd06..6403575 100644 --- a/src/recast/fortran/expr.py +++ b/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,22 +102,85 @@ 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) @@ -101,7 +188,11 @@ def render( 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}") diff --git a/src/recast/fortran/tree.py b/src/recast/fortran/tree.py index e05173e..1068cc9 100644 --- a/src/recast/fortran/tree.py +++ b/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,6 +157,10 @@ 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: @@ -162,16 +168,17 @@ def _evaluate( 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()) diff --git a/src/recast/fortran/use.py b/src/recast/fortran/use.py index ecb0afa..42f64c5 100644 --- a/src/recast/fortran/use.py +++ b/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( diff --git a/src/recast/transform/numpy/constants.py b/src/recast/transform/numpy/constants.py index a63a852..7ffe6c5 100644 --- a/src/recast/transform/numpy/constants.py +++ b/src/recast/transform/numpy/constants.py @@ -23,7 +23,7 @@ from pathlib import PurePath, PurePosixPath from typing import Any -from recast.fortran.expr import Expr, render +from recast.fortran.expr import Expr, python_call, render, with_integer_division __all__ = [ "constants_module", @@ -299,20 +299,12 @@ def use_constants_module(resolved: list[dict[str, Any]], module_name: str) -> st def _python(expr: Expr) -> str: - text = render( - expr, + # Fortran divides two integers to an integer; ``with_integer_division`` + # spells those quotients ``//`` from the tree's own types. + return render( + with_integer_division(expr), real=lambda text: f"np.float64('{text}')", integer=lambda text: text, name=lambda text: text.upper(), + call=python_call, ) - return _integer_division(text) - - -def _integer_division(text: str) -> str: - """Fortran divides two integers to an integer. An expression with no - real literal in it is integer arithmetic throughout (a name in it is an - integer parameter, or it would have been a real one), and ``/`` has to - be ``//`` -- ``nrk = runge_kutta_type / 10`` is 4, not 4.1.""" - if "np.float64(" in text or "float(" in text: - return text - return text.replace("//", "/").replace("/", "//") diff --git a/tests/test_fortran_analysis.py b/tests/test_fortran_analysis.py index a5c048c..e86b611 100644 --- a/tests/test_fortran_analysis.py +++ b/tests/test_fortran_analysis.py @@ -502,12 +502,14 @@ def test_a_missing_constant_fails_rather_than_vanishing(tmp_path: Path) -> None: def test_an_initializer_too_rich_to_model_refuses(tmp_path: Path) -> None: - """These are sums, products and powers over literals. A function call is - not approximated, it is declined.""" + """These are sums, products and powers over literals, and the intrinsics + both sides fold to the same bits. A transcendental is not approximated, + it is declined: gfortran folds ``exp`` with MPFR and libm need not agree + in the last bit.""" src = """\ module rich_mod implicit none - real, parameter :: weird = sqrt(2.0) + real, parameter :: weird = exp(2.0) end module rich_mod """ _write(tmp_path, "rich.f90", src) @@ -515,6 +517,65 @@ def test_an_initializer_too_rich_to_model_refuses(tmp_path: Path) -> None: use.resolve(["weird"], [tmp_path / "rich.f90"]) +CONSTS_INTRINSIC = """\ +module tol_mod + implicit none + integer, parameter :: r8 = 8 + real(r8), parameter :: pi = 3.14159265358979_r8 + real(r8), parameter :: eps = max( 1.0e-10_r8, epsilon(pi) ) + real(r8), parameter :: tol = max( 1.e-10_r8, epsilon(tol) ) + real(r8), parameter :: three = real( 3, kind = r8 ) + real(r8), parameter :: root = sqrt( 2.0_r8 ) + integer, parameter :: half = int( 7.9_r8 ) / 2 +end module tol_mod +""" + + +def test_intrinsic_calls_in_initializers_fold_the_same_on_both_sides(tmp_path: Path) -> None: + """CLUBB's ``eps = max( 1.0e-10, epsilon(pi) )`` (constants_clubb) and + ``bicgstab_tol = max( 1.e-10, epsilon(bicgstab_tol) )``: a whitelist of + intrinsics the compiler and NumPy fold to the same value, including the + one legal self-reference, a kind inquiry on the constant being declared.""" + import numpy as np + + from recast.transform.numpy.constants import use_constants_module + + _write(tmp_path, "tol.f90", CONSTS_INTRINSIC) + resolved = use.resolve(["eps", "tol", "three", "root", "half"], [tmp_path / "tol.f90"]) + got = {r["name"]: r["expr"] for r in resolved} + assert got["eps"].kind == "call" and got["eps"].text == "max" + # The self-reference is gone from the tree, so it is not a dependency. + assert "tol" not in expr.names_used(got["tol"]) + scope: dict[str, object] = {} + exec(use_constants_module(resolved, "tol_mod"), scope) # generated text under test + assert scope["EPS"] == np.float64(1.0e-10) + assert scope["TOL"] == np.float64(1.0e-10) + assert scope["THREE"] == np.float64(3.0) and type(scope["THREE"]) is np.float64 + assert scope["ROOT"] == np.sqrt(np.float64(2.0)) + assert scope["HALF"] == 3 + + +def test_a_renderer_without_a_call_spelling_refuses_a_call(tmp_path: Path) -> None: + _write(tmp_path, "tol.f90", CONSTS_INTRINSIC) + resolved = use.resolve(["eps"], [tmp_path / "tol.f90"]) + tree = {r["name"]: r["expr"] for r in resolved}["eps"] + with pytest.raises(expr.UnsupportedExpression): + expr.render(tree, real=str, integer=str, name=str) + + +def test_a_conversion_to_a_kind_the_fold_cannot_honour_refuses(tmp_path: Path) -> None: + src = """\ +module sp_mod + implicit none + integer, parameter :: sp = 4 + real(sp), parameter :: x = real( 3, kind = sp ) +end module sp_mod +""" + _write(tmp_path, "sp.f90", src) + with pytest.raises(expr.UnsupportedExpression): + use.resolve(["x"], [tmp_path / "sp.f90"]) + + # --- block boundaries as shared vocabulary ----------------------------------- @@ -1716,6 +1777,29 @@ def test_character_parameters_fold_and_fit_or_stay_a_skip() -> None: assert char_length("CHARACTER", "c*6") == 6 +def test_a_kind_inquiry_on_the_constant_itself_is_kept(tmp_path: Path) -> None: + """``tol = max( 1.e-10_core_rknd, epsilon(tol) )`` (CLUBB's + penta_bicgstab_solver): the only legal self-reference, asking the + constant's own kind. It cannot be a ``ref`` to a name not yet defined, and + the target renders the kind's epsilon with no argument.""" + src = """\ +module self_mod + implicit none + integer, parameter :: core_rknd = 8 + real( kind = core_rknd ), parameter :: tol = max( 1.e-10_core_rknd, epsilon(tol) ) +end module self_mod +""" + _write(tmp_path, "self.f90", src) + payload = _param(constants.extract(tmp_path / "self.f90"), "tol")["payload"] + assert payload == [ + { + "t": "call", + "v": "max", + "args": [[{"t": "real", "v": "1.e-10"}], [{"t": "call", "v": "epsilon", "args": []}]], + } + ] + + CALLBACK = """\ module callback_mod implicit none