Bugfix for Fortran+C++ code generation

Author: satfra

modules/COEN/CppOptimize.m

@@ -50,6 +50,44 @@ These are global memory reads (~400 cycle latency) and must be computed once.
         <|"Expr" -> newExpr, "Definitions" -> Table[{names[[i]], uniqueCalls[[i]]}, {i, 1, Length[uniqueCalls]}], "Count" -> Length[uniqueCalls]|>
     ];
 
+(**********************************************************************************
+    Pass 1b: Composite-Denominator Hoisting
+    Extract every distinct negative-integer power of a composite (Plus/Times) base
+    -- the regulated-propagator denominators 1/(M^2 + q*(...)^2) etc. These are NOT
+    interpolator leaves, so without global hoisting they fall to the per-sub-kernel
+    CSE (Pass 3) and are recomputed once per chunk when the term-sum is split.
+    Hoisting them to shared "_den" defs (like interpolators) computes each once.
+
+    Bases are sorted by LeafCount ascending so a denominator nested inside another
+    gets a lower index; each def body then substitutes the already-assigned inner
+    "_denK" (rules without its own), making the def list topologically ordered.
+    Runs AFTER interpolator hoisting, so bases already carry "_interp" placeholders.
+**********************************************************************************)
+
+hoistDivisions[expr_] :=
+    Module[{cands, uniqueDens, names, rules, newExpr, defs},
+        If[Not @ TrueQ[$codeHoistDivisions],
+            Return[<|"Expr" -> expr, "Definitions" -> {}, "Count" -> 0|>]
+        ];
+        (* composite reciprocals: Power[base, n<0] with base a Plus or Times (a shared,
+           non-trivial denominator), not a bare symbol/number/placeholder *)
+        cands = Cases[expr,
+            Power[b_, n_Integer] /; n < 0 && (Head[b] === Plus || Head[b] === Times),
+            Infinity];
+        (* inner (smaller) denominators first -> def list is dependency-ordered *)
+        uniqueDens = SortBy[DeleteDuplicates[cands], LeafCount];
+        If[Length[uniqueDens] === 0,
+            Return[<|"Expr" -> expr, "Definitions" -> {}, "Count" -> 0|>]
+        ];
+        names = Table["_den" <> ToString[i], {i, 1, Length[uniqueDens]}];
+        rules = Thread[uniqueDens -> names];
+        newExpr = expr //. rules;
+        (* def body for _denI: replace NESTED denominators (all rules except its own,
+           which would otherwise collapse the whole body to its own name) *)
+        defs = Table[{names[[i]], uniqueDens[[i]] //. Drop[rules, {i}]}, {i, 1, Length[uniqueDens]}];
+        <|"Expr" -> newExpr, "Definitions" -> defs, "Count" -> Length[uniqueDens]|>
+    ];
+
 (**********************************************************************************
     Power Basis Normalization
     After CSE, rewrite Power[base, m] in terms of an already-hoisted Power[base, n]
@@ -407,11 +445,35 @@ Partitions interps into shared (referenced by 2+ kernels) and local (1 kernel).
         interpNames = #[[1]]& /@ interpDefs;
         interpsByName = Association @ Table[interpDefs[[i, 1]] -> interpDefs[[i]], {i, Length[interpDefs]}];
         Module[
-            {chunkRefs, useCounts, sharedNames, sharedDefs}
+            {chunkRefs, useCounts, sharedNames, sharedDefs, defsByName}
             ,
-            (* Find which interp names each chunk references *)
-            chunkRefs = Map[Intersection[interpNames, DeleteDuplicates @ Cases[#, _String, Infinity]]&, chunks];
-            (* Count how many chunks reference each interp *)
+            (* name -> definition, for transitive reference expansion *)
+            defsByName = Association @ Table[interpDefs[[i, 1]] -> interpDefs[[i]], {i, Length[interpDefs]}];
+            (* Find which def names each chunk references -- TRANSITIVELY: a chunk that
+               references _den5 also "uses" every _interp/_den that _den5's body reads,
+               so a denominator shared across chunks (and its inputs) is correctly counted
+               as shared rather than dropped or recomputed per chunk. *)
+            chunkRefs =
+                Map[
+                    Function[{chunk},
+                        Module[{refs, prevLen = -1},
+                            (* level {0, Infinity}: a chunk (or def body) that IS a bare
+                               placeholder string -- e.g. a term that is a single hoisted
+                               reciprocal _denK with no other factors -- sits at level 0 and
+                               would be missed by the default Infinity (= {1, Infinity}). *)
+                            refs = Intersection[interpNames, DeleteDuplicates @ Cases[chunk, _String, {0, Infinity}]];
+                            While[Length[refs] =!= prevLen,
+                                prevLen = Length[refs];
+                                refs = DeleteDuplicates @ Join[refs,
+                                    Intersection[interpNames,
+                                        Flatten @ Map[Cases[defsByName[#][[2]], _String, {0, Infinity}]&, refs]]];
+                            ];
+                            refs
+                        ]
+                    ],
+                    chunks
+                ];
+            (* Count how many chunks reference each interp/den *)
             useCounts = Counts[Flatten[chunkRefs]];
             sharedNames = Keys @ Select[useCounts, # > 1&];
             sharedDefs = Select[interpDefs, MemberQ[sharedNames, #[[1]]]&];
@@ -583,7 +645,7 @@ Partitions interps into shared (referenced by 2+ kernels) and local (1 kernel).
 **********************************************************************************)
 
 optimizeExpression[equation_] :=
-    Module[{expr, interpResult, interpCount, splitResult, sharedDefs, subKernels, result, allDefs},
+    Module[{expr, interpResult, interpCount, divResult, divCount, globalDefs, splitResult, sharedDefs, subKernels, result, allDefs},
         FunKitDebug[1, "Starting optimization pipeline (optimize = ", $codeOptimize, ")"];
         (* If optimization is disabled, return raw expression with no passes *)
         If[!TrueQ[$codeOptimize],
@@ -596,8 +658,15 @@ Partitions interps into shared (referenced by 2+ kernels) and local (1 kernel).
         expr = interpResult["Expr"];
         interpCount = interpResult["Count"];
         FunKitDebug[2, "Hoisted ", interpCount, " interpolator calls"];
+        (* Pass 1b: Composite-denominator hoisting — share reciprocals across chunks *)
+        divResult = cgTimed[$ProfileCgHoist, hoistDivisions[expr]];
+        expr = divResult["Expr"];
+        divCount = divResult["Count"];
+        FunKitDebug[2, "Hoisted ", divCount, " composite denominators"];
+        (* global defs: interps FIRST, then denominators (which reference interps) *)
+        globalDefs = Join[interpResult["Definitions"], divResult["Definitions"]];
         (* Pass 2: Early split decision *)
-        splitResult = cgTimed[$ProfileCgSplit, earlySplit[interpResult["Definitions"], expr]];
+        splitResult = cgTimed[$ProfileCgSplit, earlySplit[globalDefs, expr]];
         FunKitDebug[2, "Early split: ", splitResult["Split"]];
         (* === PER-KERNEL PASSES === *)
         If[TrueQ[splitResult["Split"]],
@@ -623,8 +692,8 @@ Partitions interps into shared (referenced by 2+ kernels) and local (1 kernel).
             <|"UseSubKernels" -> True, "SharedDefinitions" -> sharedDefs, "SubKernels" -> subKernels|>
             ,
             (* Single-kernel path: optimize the whole expression *)
-            result = optimizeSubKernel[expr, interpCount];
-            allDefs = Join[interpResult["Definitions"], result["Definitions"]];
+            result = optimizeSubKernel[expr, interpCount + divCount];
+            allDefs = Join[globalDefs, result["Definitions"]];
             FunKitDebug[2, "Single-kernel optimization complete: ", Length[allDefs], " total defs"];
             (* Try splitting for registers if expression is large *)
             splitResult = cgTimed[$ProfileCgSplit, splitIntoSubKernels[allDefs, result["Expr"]]];

modules/COEN/Fortran.m

@@ -31,20 +31,20 @@
     ];
 
 (* Fortran identifiers must start with a letter. The optimizer produces names
-   like _interp1, _cse1, _tran1, _result1 which are invalid. The sub-kernel
+   like _interp1, _cse1, _tran1, _result1, _den1 which are invalid. The sub-kernel
    accumulator _acc and the ReturnTransform placeholder postAcc<n> are also
    coerced. This function renames them by stripping any leading underscore and
    adding an "fk" prefix. *)
 fortranFixNames[code_String] :=
     StringReplace[code, {
-        "_" ~~ prefix:("interp"|"cse"|"tran"|"result") ~~ num:DigitCharacter.. :> "fk" <> prefix <> num,
+        "_" ~~ prefix:("interp"|"cse"|"tran"|"result"|"den") ~~ num:DigitCharacter.. :> "fk" <> prefix <> num,
         WordBoundary ~~ "_acc" ~~ WordBoundary :> "fkacc"
     }];
 
 (* Generate double precision declarations for all optimizer variables in the code *)
 fortranVarDeclarations[code_String] :=
     Module[{numbered, accVar, vars},
-        numbered = Union @ StringCases[code, "fk" ~~ ("interp"|"cse"|"tran"|"result") ~~ DigitCharacter..];
+        numbered = Union @ StringCases[code, "fk" ~~ ("interp"|"cse"|"tran"|"result"|"den") ~~ DigitCharacter..];
         accVar = If[StringContainsQ[code, WordBoundary ~~ "fkacc" ~~ WordBoundary], {"fkacc"}, {}];
         vars = Join[numbered, accVar];
         If[Length[vars] === 0, "",

modules/COEN/Tools.m

@@ -188,6 +188,12 @@ Codegen profiling instrumentation (zero overhead unless enabled)
 
 $codeOptimizeInterps = {a_Symbol[__] /; Not @ MatchQ[a, Times | Plus | Power | Rational | Complex | Real | Integer]};
 
+(* Composite-denominator (negative integer power of a Plus/Times base) global hoisting.
+   These reciprocals are NOT interpolator leaves, so without this they fall to the
+   PER-SUB-KERNEL CSE and get recomputed once per chunk when the term-sum is split.
+   Hoisting them to shared defs (like interpolators) computes each once across chunks. *)
+$codeHoistDivisions = True;
+
 $availableRegisters = 32;
 
 $codeOptimize = True;

tests/COEN/CppTests.m

@@ -268,6 +268,112 @@ int main () {
 
 AppendTo[tests, VerificationTest[StringContainsQ[tranCode, "_tran"], True, TestID -> "Verify transcendental hoisting creates _tran variables"]];
 
+(**********************************************************************************
+    Composite-denominator hoisting (Pass 1b)
+    Negative integer powers of a Plus/Times base — the regulated-propagator
+    denominators 1/(M^2 + q*(...)^2) — are hoisted to shared "_den" defs so a
+    denominator shared across split sub-kernels is computed once, not recomputed
+    per chunk by the per-sub-kernel CSE. See CppOptimize.m hoistDivisions and
+    Tools.m $codeHoistDivisions.
+**********************************************************************************)
+
+ClearAll[a, mm, b, c];
+
+(* A composite denominator is hoisted to a _den variable. *)
+
+Block[{FunKit`Private`$codeOptimize = True, FunKit`Private`$codeHoistDivisions = True},
+    denHoistCode = CppCode[Sin[a] / (mm + a^2) + Cos[a] / (mm + a^2)];
+];
+
+AppendTo[tests, VerificationTest[StringContainsQ[denHoistCode, "_den"], True, TestID -> "Composite denominator hoisting creates _den variables"]];
+
+(* With hoisting disabled the same denominator falls to the per-kernel CSE, so no
+   _den appears — guards the $codeHoistDivisions toggle. *)
+
+Block[{FunKit`Private`$codeOptimize = True, FunKit`Private`$codeHoistDivisions = False},
+    denNoHoistCode = CppCode[Sin[a] / (mm + a^2) + Cos[a] / (mm + a^2)];
+];
+
+AppendTo[tests, VerificationTest[StringFreeQ[denNoHoistCode, "_den"], True, TestID -> "No _den variables when $codeHoistDivisions is False"]];
+
+(* A denominator shared across many terms that split into sub-kernels must be
+   hoisted to exactly ONE shared def (the efficiency invariant: not recomputed
+   once per chunk). Exercises the transitive reference tracking in earlySplit. *)
+
+ClearAll[a];
+denSharedExpr = Sum[Sin[a + i] Cos[a - i] / (3 + a^2), {i, 1, 600}];
+
+Block[{FunKit`Private`$codeOptimize = True, FunKit`Private`$codeHoistDivisions = True, FunKit`Private`$codeMaxKernelTerms = 200},
+    denSharedCode = CppCode[denSharedExpr];
+];
+
+AppendTo[tests, VerificationTest[StringContainsQ[denSharedCode, "// subkernel 1"], True, TestID -> "Shared denominator: expression does split into sub-kernels"]];
+
+AppendTo[tests, VerificationTest[Length @ StringCases[denSharedCode, "const auto _den" ~~ DigitCharacter.. ~~ " ="], 1, TestID -> "Shared denominator hoisted to a single def across sub-kernels"]];
+
+(* ...and the split kernel still computes the right value. *)
+
+Block[{FunKit`Private`$codeOptimize = True, FunKit`Private`$codeHoistDivisions = True, FunKit`Private`$codeMaxKernelTerms = 200},
+    funBodyDen = MakeCppFunction[denSharedExpr, "Name" -> "funDen", "Body" -> "using namespace std; const auto a = in;", "Parameters" -> {"in"}];
+];
+
+execDen = CreateExecutable["
+#include <iostream>
+#include <iomanip>
+#include <cmath>
+using NumberType = double;
+
+" <> fmaCode <> "
+" <> powrCode <> "
+" <> funBodyDen <> "
+
+int main () {
+  std::cout << std::setprecision (10) << funDen (1.5) << std::endl;
+}
+", "FunKitCppTestDen", "CompilerName" -> CppCompiler, "SystemCompileOptions" -> "-std=c++20"];
+
+outputDen = Import["!" <> QuoteFile[execDen], "Text"];
+
+expectedDen = ToString[NumberForm[denSharedExpr /. a -> 1.5, 10]];
+
+AppendTo[tests, VerificationTest[execDen =!= $Failed, True, TestID -> "Shared denominator hoisting compiles with sub-kernel splitting"]];
+
+AppendTo[tests, VerificationTest[outputDen, expectedDen, TestID -> "Shared denominator hoisting preserves numerical value across sub-kernels"]];
+
+(* Nested denominators: inner is hoisted first and the outer def references it
+   (topological ordering by LeafCount + nested substitution). Compile & run to
+   confirm the dependency-ordered defs produce the correct value. *)
+
+ClearAll[a, b, c];
+denNestedExpr = 1 / (b + 1 / (c + a^2)) + a / (c + a^2);
+
+Block[{FunKit`Private`$codeOptimize = True, FunKit`Private`$codeHoistDivisions = True},
+    funBodyNested = MakeCppFunction[denNestedExpr /. {b -> 0.7, c -> 1.3}, "Name" -> "funNested", "Body" -> "using namespace std; const auto a = in;", "Parameters" -> {"in"}];
+];
+
+execNested = CreateExecutable["
+#include <iostream>
+#include <iomanip>
+#include <cmath>
+using NumberType = double;
+
+" <> fmaCode <> "
+" <> powrCode <> "
+" <> funBodyNested <> "
+
+int main () {
+  std::cout << std::setprecision (10) << funNested (1.5) << std::endl;
+}
+", "FunKitCppTestNested", "CompilerName" -> CppCompiler, "SystemCompileOptions" -> "-std=c++20"];
+
+outputNested = Import["!" <> QuoteFile[execNested], "Text"];
+
+expectedNested = ToString[NumberForm[denNestedExpr /. {b -> 0.7, c -> 1.3, a -> 1.5}, 10]];
+
+AppendTo[tests, VerificationTest[execNested =!= $Failed, True, TestID -> "Nested denominator hoisting compiles"]];
+
+AppendTo[tests, VerificationTest[outputNested, expectedNested, TestID -> "Nested denominator hoisting preserves numerical value"]];
+
 (**********************************************************************************
     ReturnTransform option (post-processing injection)
 **********************************************************************************)

tests/COEN/FortranTests.m

@@ -121,6 +121,73 @@ Basic Fortran Function Test (MakeFortranFunction)
 
 AppendTo[tests, VerificationTest[Abs[fortranOptVal - expectedOptVal] < 1*^-8, True, TestID -> "Optimized Fortran function returns correct value"]];
 
+(**********************************************************************************
+    Composite-denominator hoisting: _den names must be coerced to valid Fortran
+    identifiers (fkden), since Fortran identifiers cannot start with "_". See
+    Fortran.m fortranFixNames / fortranVarDeclarations and CppOptimize.m
+    hoistDivisions.
+**********************************************************************************)
+
+ClearAll[a, mm]
+
+Block[{FunKit`Private`$codeOptimize = True, FunKit`Private`$codeHoistDivisions = True},
+    fortranDenCode = FortranCode[Sin[a] / (mm + a^2) + Cos[a] / (mm + a^2)];
+];
+
+AppendTo[tests, VerificationTest[
+    StringContainsQ[fortranDenCode, "fkden"],
+    True,
+    TestID -> "Composite denominator hoisting produces fkden variables in Fortran"
+]];
+
+(* No leading-underscore _den survives the rename (would be invalid Fortran). *)
+
+AppendTo[tests, VerificationTest[
+    StringFreeQ[fortranDenCode, "_den"],
+    True,
+    TestID -> "No invalid _den identifiers remain in Fortran output"
+]];
+
+(* ...and the result is numerically correct end-to-end. *)
+
+funBodyDen = MakeFortranFunction[Sin[a] / (mm + a^2) + Cos[a] / (mm + a^2), "Name" -> "funden",
+    "Body" -> "double precision :: a, mm\na = in1\nmm = in2",
+    "Parameters" -> {"in1", "in2"}];
+
+(* The hoisted denominator is declared with double precision (required: the
+   generated function uses implicit none). *)
+
+AppendTo[tests, VerificationTest[
+    StringContainsQ[funBodyDen, "double precision" ~~ Shortest[__] ~~ "fkden"],
+    True,
+    TestID -> "Hoisted denominator gets a double precision declaration in Fortran"
+]];
+
+codeDen = funBodyDen <> "
+
+program main
+  implicit none
+  double precision :: res, funden
+  res = funden(1.2d0, 0.7d0)
+  write(*,'(F25.15)') res
+end program main
+";
+
+execFileDen = $TemporaryDirectory <> "/FunKitFortranTestDen.f90";
+execPathDen = $TemporaryDirectory <> "/FunKitFortranTestDen";
+Export[execFileDen, codeDen, "Text"];
+
+compileDen = RunProcess[{"gfortran", "-ffree-form", "-o", execPathDen, execFileDen}];
+
+AppendTo[tests, VerificationTest[compileDen["ExitCode"], 0, TestID -> "Verify compilation of Fortran function with hoisted denominator"]];
+
+outputDen = If[compileDen["ExitCode"] === 0, RunProcess[{execPathDen}], <|"StandardOutput" -> ""|>];
+
+expectedDenVal = N[(Sin[a] / (mm + a^2) + Cos[a] / (mm + a^2)) /. {a -> 1.2, mm -> 0.7}, 15];
+fortranDenVal = ToExpression[StringTrim[outputDen["StandardOutput"]]];
+
+AppendTo[tests, VerificationTest[Abs[fortranDenVal - expectedDenVal] < 1*^-8, True, TestID -> "Fortran function with hoisted denominator returns correct value"]];
+
 (**********************************************************************************
     Simple expression: plain return, no CSE
 **********************************************************************************)