perf: optimize small integer powers in generated targets

codegen/src/render/c.rs

@@ -1,6 +1,6 @@
 use std::fmt;
 
-use crate::semantic::{BinaryOp, Expr, MathFunction, Reference, UnaryOp, Variable};
+use crate::semantic::{BinaryOp, Expr, MathFunction, Number, Reference, UnaryOp, Variable};
 
 #[derive(Clone, Copy)]
 pub(crate) enum Dialect {
@@ -39,12 +39,27 @@ pub(crate) fn requires_math(expression: &Expr) -> bool {
         Expr::Number(_) | Expr::Reference(_) => false,
         Expr::Unary { operand, .. } => requires_math(operand),
         Expr::Binary { op, left, right } => {
-            matches!(op, BinaryOp::Power) || requires_math(left) || requires_math(right)
+            (matches!(op, BinaryOp::Power) && small_integer_exponent(right).is_none())
+                || requires_math(left)
+                || requires_math(right)
         }
         Expr::Call { .. } => true,
     }
 }
 
+pub(crate) fn requires_pow4(expression: &Expr) -> bool {
+    match expression {
+        Expr::Number(_) | Expr::Reference(_) => false,
+        Expr::Unary { operand, .. } => requires_pow4(operand),
+        Expr::Binary { op, left, right } => {
+            (matches!(op, BinaryOp::Power) && small_integer_exponent(right) == Some(4))
+                || requires_pow4(left)
+                || requires_pow4(right)
+        }
+        Expr::Call { args, .. } => args.iter().any(requires_pow4),
+    }
+}
+
 #[derive(Clone, Copy, Eq, Ord, PartialEq, PartialOrd)]
 enum Precedence {
     Sum,
@@ -116,11 +131,19 @@ impl Expression<'_> {
                     self.write_binary(formatter, left, right, Precedence::Product, "/")?
                 }
                 BinaryOp::Power => {
-                    write!(formatter, "{}(", self.math_name("pow"))?;
-                    self.write_expression(formatter, left, None)?;
-                    write!(formatter, ", ")?;
-                    self.write_expression(formatter, right, None)?;
-                    write!(formatter, ")")?;
+                    if small_integer_exponent(right) == Some(4) {
+                        write!(formatter, "ptfkit_pow4(")?;
+                        self.write_expression(formatter, left, None)?;
+                        write!(formatter, ")")?;
+                    } else if let Some(exponent) = small_integer_exponent(right) {
+                        self.write_small_integer_power(formatter, left, exponent)?;
+                    } else {
+                        write!(formatter, "{}(", self.math_name("pow"))?;
+                        self.write_expression(formatter, left, None)?;
+                        write!(formatter, ", ")?;
+                        self.write_expression(formatter, right, None)?;
+                        write!(formatter, ")")?;
+                    }
                 }
             },
             Expr::Call { function, args } => {
@@ -154,6 +177,21 @@ impl Expression<'_> {
         self.write_expression(formatter, right, Some((precedence, true)))
     }
 
+    fn write_small_integer_power(
+        &self,
+        formatter: &mut fmt::Formatter<'_>,
+        base: &Expr,
+        exponent: usize,
+    ) -> fmt::Result {
+        for index in 0..exponent {
+            if index > 0 {
+                write!(formatter, " * ")?;
+            }
+            self.write_expression(formatter, base, Some((Precedence::Product, index > 0)))?;
+        }
+        Ok(())
+    }
+
     fn function_name(&self, function: MathFunction) -> &'static str {
         match function {
             MathFunction::Sqrt => self.math_name("sqrt"),
@@ -202,16 +240,33 @@ fn precedence(expression: &Expr) -> Precedence {
             op: UnaryOp::Plus,
             operand,
         } => precedence(operand),
+        Expr::Binary {
+            op: BinaryOp::Power,
+            right,
+            ..
+        } if small_integer_exponent(right).is_some() => Precedence::Product,
         Expr::Number(_) | Expr::Reference(_) | Expr::Binary { .. } | Expr::Call { .. } => {
             Precedence::Primary
         }
     }
 }
 
+fn small_integer_exponent(expression: &Expr) -> Option<usize> {
+    let Expr::Number(Number { value, .. }) = expression else {
+        return None;
+    };
+    match *value {
+        2.0 => Some(2),
+        3.0 => Some(3),
+        4.0 => Some(4),
+        _ => None,
+    }
+}
+
 #[cfg(test)]
 mod tests {
     use super::*;
-    use crate::semantic::{BinaryOp, Reference};
+    use crate::semantic::{BinaryOp, Number, Reference};
 
     fn inputs() -> Vec<String> {
         vec!["x".into(), "y".into(), "z".into()]
@@ -221,6 +276,13 @@ mod tests {
         Expr::Reference(Reference::Input(index))
     }
 
+    fn number(value: f64) -> Expr {
+        Expr::Number(Number {
+            value,
+            lexeme: value.to_string(),
+        })
+    }
+
     #[test]
     fn preserves_precedence_for_c_and_cpp() {
         let expression = Expr::Binary {
@@ -330,4 +392,42 @@ mod tests {
             "std::fmin(std::pow(x + y, -z), std::sqrt(x * y))"
         );
     }
+
+    #[test]
+    fn renders_small_integer_powers_as_multiplication_for_all_c_dialects() {
+        for (exponent, expected) in [(2.0, "x * x"), (3.0, "x * x * x"), (4.0, "ptfkit_pow4(x)")] {
+            let expression = Expr::Binary {
+                op: BinaryOp::Power,
+                left: Box::new(input(0)),
+                right: Box::new(number(exponent)),
+            };
+            assert_eq!(
+                super::expression(&expression, &inputs(), &[], Dialect::C).to_string(),
+                expected
+            );
+            assert_eq!(
+                super::expression(&expression, &inputs(), &[], Dialect::Cpp).to_string(),
+                expected
+            );
+            assert!(!requires_math(&expression));
+            assert_eq!(requires_pow4(&expression), exponent == 4.0);
+        }
+    }
+
+    #[test]
+    fn parenthesizes_small_integer_powers_when_required_by_division() {
+        let expression = Expr::Binary {
+            op: BinaryOp::Divide,
+            left: Box::new(input(0)),
+            right: Box::new(Expr::Binary {
+                op: BinaryOp::Power,
+                left: Box::new(input(1)),
+                right: Box::new(number(2.0)),
+            }),
+        };
+        assert_eq!(
+            super::expression(&expression, &inputs(), &[], Dialect::C).to_string(),
+            "x / (y * y)"
+        );
+    }
 }

codegen/src/targets/native.rs

@@ -79,8 +79,14 @@ fn c_header(slug: &str, functions: &[&CompiledFunction]) -> Result<String> {
     writer.write(format_args!(
         "{HEADER}\n\n#ifndef {guard}\n#define {guard}\n\n"
     ));
+    if requires_pow4(functions) {
+        writer.write("#include <ptfkit/detail/power.h>\n");
+    }
     if requires_math(functions) {
-        writer.write("#include <math.h>\n\n");
+        writer.write("#include <math.h>\n");
+    }
+    if requires_pow4(functions) || requires_math(functions) {
+        writer.blank_line();
     }
     let first = functions
         .first()
@@ -116,8 +122,15 @@ fn c_header(slug: &str, functions: &[&CompiledFunction]) -> Result<String> {
 fn cpp_module(slug: &str, functions: &[&CompiledFunction]) -> Result<String> {
     let mut writer = Writer::new();
     writer.write(format_args!("{HEADER}\n\n"));
-    if requires_math(functions) {
-        writer.write("module;\n#include <cmath>\n\n");
+    if requires_pow4(functions) || requires_math(functions) {
+        writer.write("module;\n");
+        if requires_pow4(functions) {
+            writer.write("#include <ptfkit/detail/power.h>\n");
+        }
+        if requires_math(functions) {
+            writer.write("#include <cmath>\n");
+        }
+        writer.blank_line();
     }
     writer.write(format_args!("export module ptfkit.{slug};\n\n"));
     let first = functions
@@ -479,6 +492,16 @@ fn requires_math(functions: &[&CompiledFunction]) -> bool {
     })
 }
 
+fn requires_pow4(functions: &[&CompiledFunction]) -> bool {
+    functions.iter().any(|function| {
+        function
+            .ir
+            .variables
+            .iter()
+            .any(|variable| c::requires_pow4(&variable.expression))
+    })
+}
+
 fn c_test(slug: &str, functions: &[&CompiledFunction]) -> Result<String> {
     c_compatibility_test(slug, functions)
 }
@@ -614,6 +637,14 @@ mod tests {
             right: Box::new(value.clone()),
         };
         let power = Expr::Binary {
+            op: BinaryOp::Power,
+            left: Box::new(value.clone()),
+            right: Box::new(Expr::Number(crate::semantic::Number {
+                value: 2.0,
+                lexeme: "2".to_owned(),
+            })),
+        };
+        let non_optimized_power = Expr::Binary {
             op: BinaryOp::Power,
             left: Box::new(value.clone()),
             right: Box::new(value.clone()),
@@ -624,7 +655,8 @@ mod tests {
         };
 
         assert!(!c::requires_math(&arithmetic));
-        assert!(c::requires_math(&power));
+        assert!(!c::requires_math(&power));
+        assert!(c::requires_math(&non_optimized_power));
         assert!(c::requires_math(&logarithm));
     }
 

targets/ptfkit-native/cpp/CMakeLists.txt

@@ -10,6 +10,7 @@ target_compile_definitions(ptfkit_cpp PUBLIC
 target_compile_options(ptfkit_cpp PUBLIC
     $<$<CXX_COMPILER_ID:GNU>:-fmodules-ts>
 )
+target_link_libraries(ptfkit_cpp PRIVATE ptfkit::c)
 file(GLOB_RECURSE PTFKIT_CPP_MODULES CONFIGURE_DEPENDS *.cppm)
 target_sources(ptfkit_cpp PUBLIC
     FILE_SET CXX_MODULES

targets/ptfkit-native/cpp/ferrerjulia2004.cppm

@@ -53,10 +53,9 @@ inline double calc_ptf_ferrerjulia2004_campbell_shiozawa(double sand, double cla
  */
 [[nodiscard]]
 inline double calc_ptf_ferrerjulia2004_saxton(double sand, double clay) {
-    return 10.0 *
-           std::exp(1.01 - 0.0755 * sand +
-                    (-3.895 + 0.03671 * sand - 0.1103 * clay + 0.00087546 * std::pow(clay, 2.0)) /
-                        (0.33 - 0.000751 * sand + 0.176 * std::log10(clay)));
+    return 10.0 * std::exp(1.01 - 0.0755 * sand +
+                           (-3.895 + 0.03671 * sand - 0.1103 * clay + 0.00087546 * (clay * clay)) /
+                               (0.33 - 0.000751 * sand + 0.176 * std::log10(clay)));
 }
 
 /**

targets/ptfkit-native/cpp/hodnett2002.cppm

@@ -87,18 +87,18 @@ inline Hodnett2002PTFResult calc_ptf_hodnett2002(double sand, double silt, doubl
                                                  double cation_exchange_capacity, double ph) {
     const double ln_alpha =
         (-2.294 - 3.526 * silt + 2.440 * organic_carbon - 0.076 * cation_exchange_capacity -
-         11.331 * ph + 0.019 * std::pow(silt, 2.0)) /
+         11.331 * ph + 0.019 * (silt * silt)) /
         100.0;
     const double alpha = std::exp(ln_alpha);
     const double ln_n = (62.986 - 0.833 * clay - 0.529 * organic_carbon + 0.593 * ph +
-                         0.0070 * std::pow(clay, 2.0) - 0.014 * sand * silt) /
+                         0.0070 * (clay * clay) - 0.014 * sand * silt) /
                         100.0;
     const double n = std::exp(ln_n);
     const double theta_s = (81.799 + 0.099 * clay - 31.420 * bulk_density +
                             0.018 * cation_exchange_capacity + 0.451 * ph - 0.0005 * sand * clay) /
                            100.0;
     const double theta_r = (22.733 - 0.164 * sand + 0.235 * cation_exchange_capacity - 0.831 * ph +
-                            0.0018 * std::pow(clay, 2.0) + 0.0026 * sand * clay) /
+                            0.0018 * (clay * clay) + 0.0026 * sand * clay) /
                            100.0;
     return Hodnett2002PTFResult{alpha, n, theta_s, theta_r};
 }

targets/ptfkit-native/cpp/mayr1999.cppm

@@ -71,20 +71,20 @@ struct Mayr1999PTFResult {
 [[nodiscard]]
 inline Mayr1999PTFResult calc_ptf_mayr1999(double sand, double silt, double clay,
                                            double bulk_density, double organic_carbon) {
-    const double log10_a_hc =
-        -4.9840297533 + 0.0509226283 * sand + 0.1575152771 * silt + 0.1240901644 * bulk_density -
-        0.1640033143 * organic_carbon - 0.0021767278 * std::pow(silt, 2.0) +
-        1.438224e-5 * std::pow(silt, 3.0) + 8.040715e-4 * std::pow(clay, 2.0) +
-        0.0044067117 * std::pow(organic_carbon, 2.0);
-    const double log10_inv_b_hc =
-        -0.8466880654 - 0.0046806123 * sand + 0.0092463819 * silt - 0.4542769707 * bulk_density -
-        0.0497915563 * organic_carbon + 3.294687e-4 * std::pow(sand, 2.0) -
-        1.689056e-6 * std::pow(sand, 3.0) + 0.0011225373 * std::pow(organic_carbon, 2.0);
+    const double log10_a_hc = -4.9840297533 + 0.0509226283 * sand + 0.1575152771 * silt +
+                              0.1240901644 * bulk_density - 0.1640033143 * organic_carbon -
+                              0.0021767278 * (silt * silt) + 1.438224e-5 * (silt * silt * silt) +
+                              8.040715e-4 * (clay * clay) +
+                              0.0044067117 * (organic_carbon * organic_carbon);
+    const double log10_inv_b_hc = -0.8466880654 - 0.0046806123 * sand + 0.0092463819 * silt -
+                                  0.4542769707 * bulk_density - 0.0497915563 * organic_carbon +
+                                  3.294687e-4 * (sand * sand) - 1.689056e-6 * (sand * sand * sand) +
+                                  0.0011225373 * (organic_carbon * organic_carbon);
     const double a_hc = std::pow(10.0, log10_a_hc);
     const double b_hc = std::pow(10.0, -log10_inv_b_hc);
     const double theta_s = 0.2345971971 + 0.0046614221 * sand + 0.0088163314 * silt +
                            0.0064338641 * clay - 0.3028160229 * bulk_density +
-                           1.79762e-5 * std::pow(sand, 2.0) - 3.134631e-5 * std::pow(silt, 2.0);
+                           1.79762e-5 * (sand * sand) - 3.134631e-5 * (silt * silt);
     return Mayr1999PTFResult{a_hc, b_hc, theta_s};
 }
 

targets/ptfkit-native/cpp/saxton2006.cppm

@@ -154,7 +154,8 @@ inline Saxton2006PTFResult calc_ptf_saxton2006(double sand, double clay, double
     const double theta_33_preliminary = -0.251 * sand + 0.195 * clay + 0.011 * organic_matter +
                                         0.006 * sand * organic_matter -
                                         0.027 * clay * organic_matter + 0.452 * sand * clay + 0.299;
-    const double theta_33 = theta_33_preliminary + 1.283 * std::pow(theta_33_preliminary, 2.0) -
+    const double theta_33 = theta_33_preliminary +
+                            1.283 * (theta_33_preliminary * theta_33_preliminary) -
                             0.374 * theta_33_preliminary - 0.015;
     const double theta_s_minus_33_preliminary =
         0.278 * sand + 0.034 * clay + 0.022 * organic_matter - 0.018 * sand * organic_matter -
@@ -168,7 +169,7 @@ inline Saxton2006PTFResult calc_ptf_saxton2006(double sand, double clay, double
         71.12 * sand * theta_s_minus_33_preliminary + 8.29 * clay * theta_s_minus_33_preliminary +
         14.05 * sand * clay + 27.16;
     const double air_entry_tension = air_entry_preliminary +
-                                     0.02 * std::pow(air_entry_preliminary, 2.0) -
+                                     0.02 * (air_entry_preliminary * air_entry_preliminary) -
                                      0.113 * air_entry_preliminary - 0.70;
     const double retention_b =
         (std::log(1500.0) - std::log(33.0)) / (std::log(theta_33) - std::log(theta_1500));

targets/ptfkit-native/cpp/varallyai1982.cppm

@@ -1,8 +1,5 @@
 /* @generated by ptfkit-codegen; DO NOT EDIT. */
 
-module;
-#include <cmath>
-
 export module ptfkit.varallyai1982;
 
 /**
@@ -69,8 +66,8 @@ inline Varallyai1982Parameters calc_ptf_varallyai1982_meadow(double bulk_density
                                                              double fine_sand_fraction,
                                                              double fine_fraction) {
     const double theta_0 =
-        -8.78 * std::pow(bulk_density, 2.0) + 14.46 * std::pow(fine_fraction, 2.0) + 62.85;
-    const double m = 0.576 * std::pow(fine_sand_fraction, 2.0) -
+        -8.78 * (bulk_density * bulk_density) + 14.46 * (fine_fraction * fine_fraction) + 62.85;
+    const double m = 0.576 * (fine_sand_fraction * fine_sand_fraction) -
                      1.434 * fine_sand_fraction * fine_fraction + 0.156;
     const double pf_star =
         -1.702 * fine_sand_fraction + 1.103 * bulk_density * fine_fraction + 3.749;
@@ -135,7 +132,7 @@ inline Varallyai1982Parameters calc_ptf_varallyai1982_chernozem_a(double bulk_de
 [[nodiscard]]
 inline Varallyai1982Parameters calc_ptf_varallyai1982_chernozem_b(double bulk_density,
                                                                   double fine_fraction) {
-    const double theta_0 = -62.20 * bulk_density - 49.14 * std::pow(fine_fraction, 2.0) + 140.70;
+    const double theta_0 = -62.20 * bulk_density - 49.14 * (fine_fraction * fine_fraction) + 140.70;
     const double m = 0.635 * bulk_density - 0.482;
     const double pf_star = 4.270 * bulk_density * fine_fraction + 3.509 * bulk_density - 3.075;
     return Varallyai1982Parameters{theta_0, m, pf_star};
@@ -169,7 +166,7 @@ inline Varallyai1982Parameters calc_ptf_varallyai1982_chernozem_c(double bulk_de
     const double theta_0 = -46.80 * bulk_density + 115.39;
     const double m = 0.439 * bulk_density * fine_fraction + 0.625;
     const double pf_star =
-        3.268 * bulk_density * fine_fraction + 0.865 * std::pow(bulk_density, 2.0) + 0.301;
+        3.268 * bulk_density * fine_fraction + 0.865 * (bulk_density * bulk_density) + 0.301;
     return Varallyai1982Parameters{theta_0, m, pf_star};
 }
 

targets/ptfkit-native/cpp/wang2012.cppm

@@ -78,11 +78,12 @@ inline Wang2012PTFResult calc_ptf_wang2012(double sand, double silt, double clay
         46.481 - 4.757 * soil_organic_carbon_g_per_kg - 14.028 * std::log10(clay) -
         13.991 * std::log10(sand) + 42.261 * std::log10(soil_organic_carbon_g_per_kg) -
         11.763 / sand + 19.198 / soil_organic_carbon_g_per_kg -
-        5.448 * std::pow(bulk_density, 2.0) + 0.044 * std::pow(soil_organic_carbon_g_per_kg, 2.0) +
+        5.448 * (bulk_density * bulk_density) +
+        0.044 * (soil_organic_carbon_g_per_kg * soil_organic_carbon_g_per_kg) +
         1.975 * bulk_density * soil_organic_carbon_g_per_kg;
     const double sswc_percent = 98.813 - 21.555 / bulk_density - 39.735 / silt - 2.091 / sand +
                                 3.247 / soil_organic_carbon_g_per_kg -
-                                17.096 * std::pow(bulk_density, 2.0);
+                                17.096 * (bulk_density * bulk_density);
     const double theta_s = sswc_percent / 100.0;
     const double theta_fc = fc_percent / 100.0;
     const double k_sat = k_sat_cm_per_day / 8640000.0;