From 9a5def616b0339be1a6ff326801039aded032bf3 Mon Sep 17 00:00:00 2001
From: Roland Conybeare <rconybeare@gmail.com>
Date: Sat, 4 May 2024 14:50:58 -0400
Subject: [PATCH] xo-unit: bugfix: need handling for power=-1/2

---
 include/xo/unit/natural_unit.hpp | 23 ++++++++++++++++-------
 utest/natural_unit.test.cpp      |  9 +++++----
 2 files changed, 21 insertions(+), 11 deletions(-)

diff --git a/include/xo/unit/natural_unit.hpp b/include/xo/unit/natural_unit.hpp
index 47a25bd4..c06faa4c 100644
--- a/include/xo/unit/natural_unit.hpp
+++ b/include/xo/unit/natural_unit.hpp
@@ -215,20 +215,29 @@ namespace xo {
             bpu2_rescale(const bpu<Int> & orig,
                          const scalefactor_ratio_type & new_scalefactor)
             {
+                /* we have orig,  representing qty (b.u)^p,
+                 * with b=orig.scalefactor, u=native dimension, p=orig.power
+                 */
+
                 ratio::ratio<Int> mult = (orig.scalefactor() / new_scalefactor);
 
-                /* inv: p_frac in [0, 1) */
+                /* inv: p_frac in (-1, 1) */
                 auto p_frac = orig.power().frac();
 
                 /* asof c++26: replace mult_sq with ::pow(mult, p_frac) */
                 double mult_sq = std::numeric_limits<double>::quiet_NaN();
 
-                if (p_frac.den() == 1) {
-                    mult_sq = 1.0;
-                } else if(p_frac.den() == 2) {
-                    mult_sq = mult.template convert_to<double>();
-                } else {
-                    // remaining possibilities not supported until c++26
+                /* pre-c++26 workaround */
+                {
+                    if (p_frac.den() == 1) {
+                        mult_sq = 1.0;
+                    } else if(p_frac.num() == 1 && p_frac.den() == 2) {
+                        mult_sq = mult.template convert_to<double>();
+                    } else if(p_frac.num() == -1 && p_frac.den() == 2) {
+                        mult_sq = 1.0 / mult.template convert_to<double>();
+                    } else {
+                        // remaining possibilities not supported until c++26
+                    }
                 }
 
                 ratio::ratio<Int> mult_p = mult.power(orig.power().floor());
diff --git a/utest/natural_unit.test.cpp b/utest/natural_unit.test.cpp
index b4a23d68..195e7221 100644
--- a/utest/natural_unit.test.cpp
+++ b/utest/natural_unit.test.cpp
@@ -232,12 +232,13 @@ namespace xo {
                 constexpr auto p = power_ratio_type(-1, 2);
 
                 constexpr auto orig_bpu = bpu<int64_t>(dim::time,
-                                                        scalefactor_ratio_type(360*24*3600, 1),
-                                                        power_ratio_type(-1, 2));
+                                                       scalefactor_ratio_type(360*24*3600, 1),
+                                                       power_ratio_type(-1, 2));
                 static_assert(orig_bpu.native_dim() == dim::time);
 
                 constexpr auto new_scalefactor = scalefactor_ratio_type(30*24*3600, 1);
 
+                /* orig ~ 360d volatility,  new = 30d volatility */
                 constexpr auto mult = orig_bpu.scalefactor() / new_scalefactor;
                 log && log(xtag("mult", mult));
                 static_assert(mult.num() == 12);
@@ -263,7 +264,7 @@ namespace xo {
 
                 static_assert(rr.bpu_rescaled_.power() == power_ratio_type(-1,2));
                 static_assert(rr.outer_scale_factor_ == outer_sf_exact);
-                static_assert(rr.outer_scale_sq_ == 12.0);
+                static_assert(rr.outer_scale_sq_ == 1 / 12.0);
             }
 
             /* keep spelled-out test.  Will generalize to other fractional powers when c++26 available */
@@ -302,7 +303,7 @@ namespace xo {
 
                 static_assert(rr.bpu_rescaled_.power() == power_ratio_type(-3,2));
                 static_assert(rr.outer_scale_factor_ == outer_sf_exact);
-                static_assert(rr.outer_scale_sq_ == 12.0);
+                static_assert(rr.outer_scale_sq_ == 1 / 12.0);
             }
         } /*TEST_CASE(bpu_rescale)*/