From c47785e1590bf432af9c18f65f26a408075b5949 Mon Sep 17 00:00:00 2001 From: Petr Tsymbarovich Date: Mon, 24 Aug 2026 13:46:50 +0300 Subject: [PATCH] perf: optimize small integer powers in generated targets --- codegen/src/render/c.rs | 116 ++++++++++++++++-- codegen/src/targets/native.rs | 40 +++++- targets/ptfkit-native/cpp/CMakeLists.txt | 1 + .../ptfkit-native/cpp/ferrerjulia2004.cppm | 7 +- targets/ptfkit-native/cpp/hodnett2002.cppm | 6 +- targets/ptfkit-native/cpp/mayr1999.cppm | 20 +-- targets/ptfkit-native/cpp/saxton2006.cppm | 5 +- targets/ptfkit-native/cpp/varallyai1982.cppm | 11 +- targets/ptfkit-native/cpp/wang2012.cppm | 5 +- .../include/ptfkit/detail/power.h | 9 ++ .../include/ptfkit/ferrerjulia2004.h | 2 +- .../include/ptfkit/hodnett2002.h | 6 +- .../ptfkit-native/include/ptfkit/mayr1999.h | 12 +- .../ptfkit-native/include/ptfkit/saxton2006.h | 5 +- .../include/ptfkit/varallyai1982.h | 13 +- .../ptfkit-native/include/ptfkit/wang2012.h | 14 +-- .../ptfkit-py/src/ptfkit/ferrerjulia2004.c | 2 +- targets/ptfkit-py/src/ptfkit/hodnett2002.c | 6 +- targets/ptfkit-py/src/ptfkit/mayr1999.c | 13 +- targets/ptfkit-py/src/ptfkit/saxton2006.c | 5 +- targets/ptfkit-py/src/ptfkit/ufunc.h | 5 + targets/ptfkit-py/src/ptfkit/varallyai1982.c | 9 +- targets/ptfkit-py/src/ptfkit/wang2012.c | 6 +- 23 files changed, 233 insertions(+), 85 deletions(-) create mode 100644 targets/ptfkit-native/include/ptfkit/detail/power.h diff --git a/codegen/src/render/c.rs b/codegen/src/render/c.rs index 1939645..ad29a9d 100644 --- a/codegen/src/render/c.rs +++ b/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 { + 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 { 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)" + ); + } } diff --git a/codegen/src/targets/native.rs b/codegen/src/targets/native.rs index d25856c..ab2c24f 100644 --- a/codegen/src/targets/native.rs +++ b/codegen/src/targets/native.rs @@ -79,8 +79,14 @@ fn c_header(slug: &str, functions: &[&CompiledFunction]) -> Result { writer.write(format_args!( "{HEADER}\n\n#ifndef {guard}\n#define {guard}\n\n" )); + if requires_pow4(functions) { + writer.write("#include \n"); + } if requires_math(functions) { - writer.write("#include \n\n"); + writer.write("#include \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 { fn cpp_module(slug: &str, functions: &[&CompiledFunction]) -> Result { let mut writer = Writer::new(); writer.write(format_args!("{HEADER}\n\n")); - if requires_math(functions) { - writer.write("module;\n#include \n\n"); + if requires_pow4(functions) || requires_math(functions) { + writer.write("module;\n"); + if requires_pow4(functions) { + writer.write("#include \n"); + } + if requires_math(functions) { + writer.write("#include \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 { 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)); } diff --git a/targets/ptfkit-native/cpp/CMakeLists.txt b/targets/ptfkit-native/cpp/CMakeLists.txt index c91cbdd..9ce1908 100644 --- a/targets/ptfkit-native/cpp/CMakeLists.txt +++ b/targets/ptfkit-native/cpp/CMakeLists.txt @@ -10,6 +10,7 @@ target_compile_definitions(ptfkit_cpp PUBLIC target_compile_options(ptfkit_cpp PUBLIC $<$:-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 diff --git a/targets/ptfkit-native/cpp/ferrerjulia2004.cppm b/targets/ptfkit-native/cpp/ferrerjulia2004.cppm index 02669e7..0d5436a 100644 --- a/targets/ptfkit-native/cpp/ferrerjulia2004.cppm +++ b/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))); } /** diff --git a/targets/ptfkit-native/cpp/hodnett2002.cppm b/targets/ptfkit-native/cpp/hodnett2002.cppm index 20fbf2a..b28060c 100644 --- a/targets/ptfkit-native/cpp/hodnett2002.cppm +++ b/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}; } diff --git a/targets/ptfkit-native/cpp/mayr1999.cppm b/targets/ptfkit-native/cpp/mayr1999.cppm index 10615fe..21208bf 100644 --- a/targets/ptfkit-native/cpp/mayr1999.cppm +++ b/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}; } diff --git a/targets/ptfkit-native/cpp/saxton2006.cppm b/targets/ptfkit-native/cpp/saxton2006.cppm index e5b4e9c..d9d87d3 100644 --- a/targets/ptfkit-native/cpp/saxton2006.cppm +++ b/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)); diff --git a/targets/ptfkit-native/cpp/varallyai1982.cppm b/targets/ptfkit-native/cpp/varallyai1982.cppm index 657754b..1cee58b 100644 --- a/targets/ptfkit-native/cpp/varallyai1982.cppm +++ b/targets/ptfkit-native/cpp/varallyai1982.cppm @@ -1,8 +1,5 @@ /* @generated by ptfkit-codegen; DO NOT EDIT. */ -module; -#include - 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}; } diff --git a/targets/ptfkit-native/cpp/wang2012.cppm b/targets/ptfkit-native/cpp/wang2012.cppm index d9884a9..1198947 100644 --- a/targets/ptfkit-native/cpp/wang2012.cppm +++ b/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; diff --git a/targets/ptfkit-native/include/ptfkit/detail/power.h b/targets/ptfkit-native/include/ptfkit/detail/power.h new file mode 100644 index 0000000..4e9e5cf --- /dev/null +++ b/targets/ptfkit-native/include/ptfkit/detail/power.h @@ -0,0 +1,9 @@ +#ifndef PTFKIT_DETAIL_POWER_H +#define PTFKIT_DETAIL_POWER_H + +static inline double ptfkit_pow4(double value) { + const double square = value * value; + return square * square; +} + +#endif diff --git a/targets/ptfkit-native/include/ptfkit/ferrerjulia2004.h b/targets/ptfkit-native/include/ptfkit/ferrerjulia2004.h index c9f8dae..e51be27 100644 --- a/targets/ptfkit-native/include/ptfkit/ferrerjulia2004.h +++ b/targets/ptfkit-native/include/ptfkit/ferrerjulia2004.h @@ -50,7 +50,7 @@ static inline double calc_ptf_ferrerjulia2004_campbell_shiozawa(double sand, dou */ static inline double calc_ptf_ferrerjulia2004_saxton(double sand, double clay) { return 10.0 * exp(1.01 - 0.0755 * sand + - (-3.895 + 0.03671 * sand - 0.1103 * clay + 0.00087546 * pow(clay, 2.0)) / + (-3.895 + 0.03671 * sand - 0.1103 * clay + 0.00087546 * (clay * clay)) / (0.33 - 0.000751 * sand + 0.176 * log10(clay))); } diff --git a/targets/ptfkit-native/include/ptfkit/hodnett2002.h b/targets/ptfkit-native/include/ptfkit/hodnett2002.h index cc8a353..54a3ceb 100644 --- a/targets/ptfkit-native/include/ptfkit/hodnett2002.h +++ b/targets/ptfkit-native/include/ptfkit/hodnett2002.h @@ -84,18 +84,18 @@ calc_ptf_hodnett2002(double sand, double silt, double clay, double organic_carbo double bulk_density, 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 * pow(silt, 2.0)) / + 11.331 * ph + 0.019 * (silt * silt)) / 100.0; const double alpha = exp(ln_alpha); const double ln_n = (62.986 - 0.833 * clay - 0.529 * organic_carbon + 0.593 * ph + - 0.0070 * pow(clay, 2.0) - 0.014 * sand * silt) / + 0.0070 * (clay * clay) - 0.014 * sand * silt) / 100.0; const double n = 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 * pow(clay, 2.0) + 0.0026 * sand * clay) / + 0.0018 * (clay * clay) + 0.0026 * sand * clay) / 100.0; #ifdef __cplusplus return hodnett2002_ptf_result{alpha, n, theta_s, theta_r}; diff --git a/targets/ptfkit-native/include/ptfkit/mayr1999.h b/targets/ptfkit-native/include/ptfkit/mayr1999.h index 61681c8..7077341 100644 --- a/targets/ptfkit-native/include/ptfkit/mayr1999.h +++ b/targets/ptfkit-native/include/ptfkit/mayr1999.h @@ -70,18 +70,18 @@ static inline mayr1999_ptf_result calc_ptf_mayr1999(double sand, double silt, do 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 * pow(silt, 2.0) + 1.438224e-5 * pow(silt, 3.0) + - 8.040715e-4 * pow(clay, 2.0) + - 0.0044067117 * pow(organic_carbon, 2.0); + 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 * pow(sand, 2.0) - 1.689056e-6 * pow(sand, 3.0) + - 0.0011225373 * pow(organic_carbon, 2.0); + 3.294687e-4 * (sand * sand) - 1.689056e-6 * (sand * sand * sand) + + 0.0011225373 * (organic_carbon * organic_carbon); const double a_hc = pow(10.0, log10_a_hc); const double b_hc = 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 * pow(sand, 2.0) - 3.134631e-5 * pow(silt, 2.0); + 1.79762e-5 * (sand * sand) - 3.134631e-5 * (silt * silt); #ifdef __cplusplus return mayr1999_ptf_result{a_hc, b_hc, theta_s}; #else diff --git a/targets/ptfkit-native/include/ptfkit/saxton2006.h b/targets/ptfkit-native/include/ptfkit/saxton2006.h index e767cd9..44b1d1e 100644 --- a/targets/ptfkit-native/include/ptfkit/saxton2006.h +++ b/targets/ptfkit-native/include/ptfkit/saxton2006.h @@ -152,7 +152,8 @@ static inline saxton2006_ptf_result calc_ptf_saxton2006(double sand, double clay 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 * 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 - @@ -166,7 +167,7 @@ static inline saxton2006_ptf_result calc_ptf_saxton2006(double sand, double clay 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 * 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 = (log(1500.0) - log(33.0)) / (log(theta_33) - log(theta_1500)); const double retention_a = exp(log(33.0) + retention_b * log(theta_33)); diff --git a/targets/ptfkit-native/include/ptfkit/varallyai1982.h b/targets/ptfkit-native/include/ptfkit/varallyai1982.h index c78f593..5a7545e 100644 --- a/targets/ptfkit-native/include/ptfkit/varallyai1982.h +++ b/targets/ptfkit-native/include/ptfkit/varallyai1982.h @@ -3,8 +3,6 @@ #ifndef PTFKIT_VARALLYAI1982_H #define PTFKIT_VARALLYAI1982_H -#include - /** * @brief Varallyai et al. (1982), water-retention parameter regressions for Hungarian soils. * @@ -65,9 +63,10 @@ typedef struct { static inline varallyai1982_parameters calc_ptf_varallyai1982_meadow(double bulk_density, double fine_sand_fraction, double fine_fraction) { - const double theta_0 = -8.78 * pow(bulk_density, 2.0) + 14.46 * pow(fine_fraction, 2.0) + 62.85; - const double m = - 0.576 * pow(fine_sand_fraction, 2.0) - 1.434 * fine_sand_fraction * fine_fraction + 0.156; + const double theta_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; #ifdef __cplusplus @@ -145,7 +144,7 @@ static inline varallyai1982_parameters calc_ptf_varallyai1982_chernozem_a(double */ static inline varallyai1982_parameters calc_ptf_varallyai1982_chernozem_b(double bulk_density, double fine_fraction) { - const double theta_0 = -62.20 * bulk_density - 49.14 * 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; #ifdef __cplusplus @@ -186,7 +185,7 @@ static inline varallyai1982_parameters calc_ptf_varallyai1982_chernozem_c(double 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 * pow(bulk_density, 2.0) + 0.301; + 3.268 * bulk_density * fine_fraction + 0.865 * (bulk_density * bulk_density) + 0.301; #ifdef __cplusplus return varallyai1982_parameters{theta_0, m, pf_star}; #else diff --git a/targets/ptfkit-native/include/ptfkit/wang2012.h b/targets/ptfkit-native/include/ptfkit/wang2012.h index 9997eda..cdabfe4 100644 --- a/targets/ptfkit-native/include/ptfkit/wang2012.h +++ b/targets/ptfkit-native/include/ptfkit/wang2012.h @@ -71,15 +71,15 @@ static inline wang2012_ptf_result calc_ptf_wang2012(double sand, double silt, do 0.865 / sand - 0.030 * bulk_density * silt - 0.00000995 * soil_organic_carbon_g_per_kg * altitude; const double k_sat_cm_per_day = pow(10.0, log10_k_sat_cm_per_day); - const double fc_percent = 46.481 - 4.757 * soil_organic_carbon_g_per_kg - 14.028 * log10(clay) - - 13.991 * log10(sand) + 42.261 * log10(soil_organic_carbon_g_per_kg) - - 11.763 / sand + 19.198 / soil_organic_carbon_g_per_kg - - 5.448 * pow(bulk_density, 2.0) + - 0.044 * pow(soil_organic_carbon_g_per_kg, 2.0) + - 1.975 * bulk_density * soil_organic_carbon_g_per_kg; + const double fc_percent = + 46.481 - 4.757 * soil_organic_carbon_g_per_kg - 14.028 * log10(clay) - + 13.991 * log10(sand) + 42.261 * log10(soil_organic_carbon_g_per_kg) - 11.763 / sand + + 19.198 / soil_organic_carbon_g_per_kg - 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 * 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; diff --git a/targets/ptfkit-py/src/ptfkit/ferrerjulia2004.c b/targets/ptfkit-py/src/ptfkit/ferrerjulia2004.c index 3c8e67c..af32d66 100644 --- a/targets/ptfkit-py/src/ptfkit/ferrerjulia2004.c +++ b/targets/ptfkit-py/src/ptfkit/ferrerjulia2004.c @@ -26,7 +26,7 @@ static void calc_ptf_ferrerjulia2004_saxton_loop(char **args, const npy_intp *di const double clay = *(const double *)args[1]; const double k_sat = 10.0 * exp(1.01 - 0.0755 * sand + - (-3.895 + 0.03671 * sand - 0.1103 * clay + 0.00087546 * pow(clay, 2.0)) / + (-3.895 + 0.03671 * sand - 0.1103 * clay + 0.00087546 * (clay * clay)) / (0.33 - 0.000751 * sand + 0.176 * log10(clay))); *(double *)args[2] = k_sat; for (int arg = 0; arg < 3; arg++) diff --git a/targets/ptfkit-py/src/ptfkit/hodnett2002.c b/targets/ptfkit-py/src/ptfkit/hodnett2002.c index 6226e1f..be5f9e8 100644 --- a/targets/ptfkit-py/src/ptfkit/hodnett2002.c +++ b/targets/ptfkit-py/src/ptfkit/hodnett2002.c @@ -14,11 +14,11 @@ static void calc_ptf_hodnett2002_loop(char **args, const npy_intp *dimensions, const double ph = *(const double *)args[6]; const double ln_alpha = (-2.294 - 3.526 * silt + 2.440 * organic_carbon - 0.076 * cation_exchange_capacity - - 11.331 * ph + 0.019 * pow(silt, 2.0)) / + 11.331 * ph + 0.019 * (silt * silt)) / 100.0; const double alpha = exp(ln_alpha); const double ln_n = (62.986 - 0.833 * clay - 0.529 * organic_carbon + 0.593 * ph + - 0.0070 * pow(clay, 2.0) - 0.014 * sand * silt) / + 0.0070 * (clay * clay) - 0.014 * sand * silt) / 100.0; const double n = exp(ln_n); const double theta_s = @@ -26,7 +26,7 @@ static void calc_ptf_hodnett2002_loop(char **args, const npy_intp *dimensions, 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 * pow(clay, 2.0) + 0.0026 * sand * clay) / + 0.831 * ph + 0.0018 * (clay * clay) + 0.0026 * sand * clay) / 100.0; *(double *)args[7] = alpha; *(double *)args[8] = n; diff --git a/targets/ptfkit-py/src/ptfkit/mayr1999.c b/targets/ptfkit-py/src/ptfkit/mayr1999.c index 073373a..21ae5c3 100644 --- a/targets/ptfkit-py/src/ptfkit/mayr1999.c +++ b/targets/ptfkit-py/src/ptfkit/mayr1999.c @@ -12,18 +12,19 @@ static void calc_ptf_mayr1999_loop(char **args, const npy_intp *dimensions, cons const double organic_carbon = *(const double *)args[4]; const double log10_a_hc = -4.9840297533 + 0.0509226283 * sand + 0.1575152771 * silt + 0.1240901644 * bulk_density - 0.1640033143 * organic_carbon - - 0.0021767278 * pow(silt, 2.0) + 1.438224e-5 * pow(silt, 3.0) + - 8.040715e-4 * pow(clay, 2.0) + - 0.0044067117 * pow(organic_carbon, 2.0); + 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 * pow(sand, 2.0) - 1.689056e-6 * pow(sand, 3.0) + - 0.0011225373 * pow(organic_carbon, 2.0); + 3.294687e-4 * (sand * sand) - + 1.689056e-6 * (sand * sand * sand) + + 0.0011225373 * (organic_carbon * organic_carbon); const double a_hc = pow(10.0, log10_a_hc); const double b_hc = 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 * pow(sand, 2.0) - 3.134631e-5 * pow(silt, 2.0); + 1.79762e-5 * (sand * sand) - 3.134631e-5 * (silt * silt); *(double *)args[5] = a_hc; *(double *)args[6] = b_hc; *(double *)args[7] = theta_s; diff --git a/targets/ptfkit-py/src/ptfkit/saxton2006.c b/targets/ptfkit-py/src/ptfkit/saxton2006.c index 06ce898..83fcadc 100644 --- a/targets/ptfkit-py/src/ptfkit/saxton2006.c +++ b/targets/ptfkit-py/src/ptfkit/saxton2006.c @@ -15,7 +15,8 @@ static void calc_ptf_saxton2006_loop(char **args, const npy_intp *dimensions, co 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 * 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 - @@ -29,7 +30,7 @@ static void calc_ptf_saxton2006_loop(char **args, const npy_intp *dimensions, co 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 * 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 = (log(1500.0) - log(33.0)) / (log(theta_33) - log(theta_1500)); const double retention_a = exp(log(33.0) + retention_b * log(theta_33)); diff --git a/targets/ptfkit-py/src/ptfkit/ufunc.h b/targets/ptfkit-py/src/ptfkit/ufunc.h index fb7e115..e03b839 100644 --- a/targets/ptfkit-py/src/ptfkit/ufunc.h +++ b/targets/ptfkit-py/src/ptfkit/ufunc.h @@ -2,6 +2,11 @@ #include +static inline double ptfkit_pow4(double value) { + const double square = value * value; + return square * square; +} + static inline int ptfkit_add_ufunc(PyObject *module, const char *name, PyUFuncGenericFunction *functions, char *types, int nin, int nout) { diff --git a/targets/ptfkit-py/src/ptfkit/varallyai1982.c b/targets/ptfkit-py/src/ptfkit/varallyai1982.c index 3b28727..b7f656f 100644 --- a/targets/ptfkit-py/src/ptfkit/varallyai1982.c +++ b/targets/ptfkit-py/src/ptfkit/varallyai1982.c @@ -9,8 +9,8 @@ static void calc_ptf_varallyai1982_meadow_loop(char **args, const npy_intp *dime const double fine_sand_fraction = *(const double *)args[1]; const double fine_fraction = *(const double *)args[2]; const double theta_0 = - -8.78 * pow(bulk_density, 2.0) + 14.46 * pow(fine_fraction, 2.0) + 62.85; - const double m = 0.576 * 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; @@ -53,7 +53,8 @@ static void calc_ptf_varallyai1982_chernozem_b_loop(char **args, const npy_intp for (index = 0; index < dimensions[0]; index++) { const double bulk_density = *(const double *)args[0]; const double fine_fraction = *(const double *)args[1]; - const double theta_0 = -62.20 * bulk_density - 49.14 * 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; *(double *)args[2] = theta_0; @@ -77,7 +78,7 @@ static void calc_ptf_varallyai1982_chernozem_c_loop(char **args, const npy_intp 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 * pow(bulk_density, 2.0) + 0.301; + 3.268 * bulk_density * fine_fraction + 0.865 * (bulk_density * bulk_density) + 0.301; *(double *)args[2] = theta_0; *(double *)args[3] = m; *(double *)args[4] = pf_star; diff --git a/targets/ptfkit-py/src/ptfkit/wang2012.c b/targets/ptfkit-py/src/ptfkit/wang2012.c index 3f253fb..240164c 100644 --- a/targets/ptfkit-py/src/ptfkit/wang2012.c +++ b/targets/ptfkit-py/src/ptfkit/wang2012.c @@ -19,12 +19,12 @@ static void calc_ptf_wang2012_loop(char **args, const npy_intp *dimensions, cons const double fc_percent = 46.481 - 4.757 * soil_organic_carbon_g_per_kg - 14.028 * log10(clay) - 13.991 * log10(sand) + 42.261 * log10(soil_organic_carbon_g_per_kg) - 11.763 / sand + - 19.198 / soil_organic_carbon_g_per_kg - 5.448 * pow(bulk_density, 2.0) + - 0.044 * pow(soil_organic_carbon_g_per_kg, 2.0) + + 19.198 / soil_organic_carbon_g_per_kg - 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 * 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;