libjxl

FORK: libjxl patches used on blog
git clone https://git.neptards.moe/blog/libjxl.git
Log | Files | Refs | Submodules | README | LICENSE

rational_polynomial_test.cc (8445B)

      1 // Copyright (c) the JPEG XL Project Authors. All rights reserved.
      2 //
      3 // Use of this source code is governed by a BSD-style
      4 // license that can be found in the LICENSE file.
      5 
      6 #include <cmath>
      7 #include <string>
      8 
      9 #undef HWY_TARGET_INCLUDE
     10 #define HWY_TARGET_INCLUDE "lib/jxl/rational_polynomial_test.cc"
     11 #include <hwy/foreach_target.h>
     12 #include <hwy/highway.h>
     13 #include <hwy/tests/hwy_gtest.h>
     14 
     15 #include "lib/jxl/base/common.h"
     16 #include "lib/jxl/base/rational_polynomial-inl.h"
     17 #include "lib/jxl/base/status.h"
     18 #include "lib/jxl/testing.h"
     19 
     20 HWY_BEFORE_NAMESPACE();
     21 namespace jxl {
     22 namespace HWY_NAMESPACE {
     23 
     24 using T = float;  // required by EvalLog2
     25 using D = HWY_FULL(T);
     26 
     27 // These templates are not found via ADL.
     28 using hwy::HWY_NAMESPACE::Add;
     29 using hwy::HWY_NAMESPACE::GetLane;
     30 using hwy::HWY_NAMESPACE::ShiftLeft;
     31 using hwy::HWY_NAMESPACE::ShiftRight;
     32 using hwy::HWY_NAMESPACE::Sub;
     33 
     34 // Generic: only computes polynomial
     35 struct EvalPoly {
     36   template <size_t NP, size_t NQ>
     37   T operator()(T x, const T (&p)[NP], const T (&q)[NQ]) const {
     38     const HWY_FULL(T) d;
     39     const auto vx = Set(d, x);
     40     const auto approx = EvalRationalPolynomial(d, vx, p, q);
     41     return GetLane(approx);
     42   }
     43 };
     44 
     45 // Range reduction for log2
     46 struct EvalLog2 {
     47   template <size_t NP, size_t NQ>
     48   T operator()(T x, const T (&p)[NP], const T (&q)[NQ]) const {
     49     const HWY_FULL(T) d;
     50     auto vx = Set(d, x);
     51 
     52     const HWY_FULL(int32_t) di;
     53     const auto x_bits = BitCast(di, vx);
     54     // Cannot handle negative numbers / NaN.
     55     JXL_DASSERT(AllTrue(di, Eq(Abs(x_bits), x_bits)));
     56 
     57     // Range reduction to [-1/3, 1/3] - 3 integer, 2 float ops
     58     const auto exp_bits = Sub(x_bits, Set(di, 0x3f2aaaab));  // = 2/3
     59     // Shifted exponent = log2; also used to clear mantissa.
     60     const auto exp_shifted = ShiftRight<23>(exp_bits);
     61     const auto mantissa = BitCast(d, Sub(x_bits, ShiftLeft<23>(exp_shifted)));
     62     const auto exp_val = ConvertTo(d, exp_shifted);
     63     vx = Sub(mantissa, Set(d, 1.0f));
     64 
     65     const auto approx = Add(EvalRationalPolynomial(d, vx, p, q), exp_val);
     66     return GetLane(approx);
     67   }
     68 };
     69 
     70 // Functions to approximate:
     71 
     72 T LinearToSrgb8Direct(T val) {
     73   if (val < 0.0) return 0.0;
     74   if (val >= 255.0) return 255.0;
     75   if (val <= 10.0 / 12.92) return val * 12.92;
     76   return 255.0 * (std::pow(val / 255.0, 1.0 / 2.4) * 1.055 - 0.055);
     77 }
     78 
     79 T SimpleGamma(T v) {
     80   static const T kGamma = 0.387494322593;
     81   static const T limit = 43.01745241042018;
     82   T bright = v - limit;
     83   if (bright >= 0) {
     84     static const T mul = 0.0383723643799;
     85     v -= bright * mul;
     86   }
     87   static const T limit2 = 94.68634353321337;
     88   T bright2 = v - limit2;
     89   if (bright2 >= 0) {
     90     static const T mul = 0.22885405968;
     91     v -= bright2 * mul;
     92   }
     93   static const T offset = 0.156775786057;
     94   static const T scale = 8.898059160493739;
     95   T retval = scale * (offset + pow(v, kGamma));
     96   return retval;
     97 }
     98 
     99 // Runs CaratheodoryFejer and verifies the polynomial using a lot of samples to
    100 // return the biggest error.
    101 template <size_t NP, size_t NQ, class Eval>
    102 T RunApproximation(T x0, T x1, const T (&p)[NP], const T (&q)[NQ],
    103                    const Eval& eval, T func_to_approx(T)) {
    104   float maxerr = 0;
    105   T lastPrint = 0;
    106   // NOLINTNEXTLINE(clang-analyzer-security.FloatLoopCounter)
    107   for (T x = x0; x <= x1; x += (x1 - x0) / 10000.0) {
    108     const T f = func_to_approx(x);
    109     const T g = eval(x, p, q);
    110     maxerr = std::max(fabsf(g - f), maxerr);
    111     if (x == x0 || x - lastPrint > (x1 - x0) / 20.0) {
    112       printf("x: %11.6f, f: %11.6f, g: %11.6f, e: %11.6f\n", x, f, g,
    113              fabs(g - f));
    114       lastPrint = x;
    115     }
    116   }
    117   return maxerr;
    118 }
    119 
    120 void TestSimpleGamma() {
    121   const T p[4 * (6 + 1)] = {
    122       HWY_REP4(-5.0646949363741811E-05), HWY_REP4(6.7369380528439771E-05),
    123       HWY_REP4(8.9376652530412794E-05),  HWY_REP4(2.1153513301520462E-06),
    124       HWY_REP4(-6.9130322970386449E-08), HWY_REP4(3.9424752749293728E-10),
    125       HWY_REP4(1.2360288207619576E-13)};
    126 
    127   const T q[4 * (6 + 1)] = {
    128       HWY_REP4(-6.6389733798591366E-06), HWY_REP4(1.3299859726565908E-05),
    129       HWY_REP4(3.8538748358398873E-06),  HWY_REP4(-2.8707687262928236E-08),
    130       HWY_REP4(-6.6897385800005434E-10), HWY_REP4(6.1428748869186003E-12),
    131       HWY_REP4(-2.5475738169252870E-15)};
    132 
    133   const T err = RunApproximation(0.77, 274.579999999999984, p, q, EvalPoly(),
    134                                  SimpleGamma);
    135   EXPECT_LT(err, 0.05);
    136 }
    137 
    138 void TestLinearToSrgb8Direct() {
    139   const T p[4 * (5 + 1)] = {
    140       HWY_REP4(-9.5357499040105154E-05), HWY_REP4(4.6761186249798248E-04),
    141       HWY_REP4(2.5708174333943594E-04),  HWY_REP4(1.5250087770436082E-05),
    142       HWY_REP4(1.1946768008931187E-07),  HWY_REP4(5.9916446295972850E-11)};
    143 
    144   const T q[4 * (4 + 1)] = {
    145       HWY_REP4(1.8932479758079768E-05), HWY_REP4(2.7312342474687321E-05),
    146       HWY_REP4(4.3901204783327006E-06), HWY_REP4(1.0417787306920273E-07),
    147       HWY_REP4(3.0084206762140419E-10)};
    148 
    149   const T err =
    150       RunApproximation(0.77, 255, p, q, EvalPoly(), LinearToSrgb8Direct);
    151   EXPECT_LT(err, 0.05);
    152 }
    153 
    154 void TestExp() {
    155   const T p[4 * (2 + 1)] = {HWY_REP4(9.6266879665530902E-01),
    156                             HWY_REP4(4.8961265681586763E-01),
    157                             HWY_REP4(8.2619259189548433E-02)};
    158   const T q[4 * (2 + 1)] = {HWY_REP4(9.6259895571622622E-01),
    159                             HWY_REP4(-4.7272457588933831E-01),
    160                             HWY_REP4(7.4802088567547664E-02)};
    161   const T err = RunApproximation(-1, 1, p, q, EvalPoly(),
    162                                  [](T x) { return static_cast<T>(exp(x)); });
    163   EXPECT_LT(err, 1E-4);
    164 }
    165 
    166 void TestNegExp() {
    167   // 4,3 is the min required for monotonicity; max error in 0,10: 751 ppm
    168   // no benefit for k>50.
    169   const T p[4 * (4 + 1)] = {
    170       HWY_REP4(5.9580258551150123E-02), HWY_REP4(-2.5073728806886408E-02),
    171       HWY_REP4(4.1561830213689248E-03), HWY_REP4(-3.1815408488900372E-04),
    172       HWY_REP4(9.3866690094906802E-06)};
    173   const T q[4 * (3 + 1)] = {
    174       HWY_REP4(5.9579108238812878E-02), HWY_REP4(3.4542074345478582E-02),
    175       HWY_REP4(8.7263562483501714E-03), HWY_REP4(1.4095109143061216E-03)};
    176 
    177   const T err = RunApproximation(0, 10, p, q, EvalPoly(),
    178                                  [](T x) { return static_cast<T>(exp(-x)); });
    179   EXPECT_LT(err, sizeof(T) == 8 ? 2E-5 : 3E-5);
    180 }
    181 
    182 void TestSin() {
    183   const T p[4 * (6 + 1)] = {
    184       HWY_REP4(1.5518122109203780E-05),  HWY_REP4(2.3388958643675966E+00),
    185       HWY_REP4(-8.6705520940849157E-01), HWY_REP4(-1.9702294764873535E-01),
    186       HWY_REP4(1.2193404314472320E-01),  HWY_REP4(-1.7373966109788839E-02),
    187       HWY_REP4(7.8829435883034796E-04)};
    188   const T q[4 * (5 + 1)] = {
    189       HWY_REP4(2.3394371422557279E+00), HWY_REP4(-8.7028221081288615E-01),
    190       HWY_REP4(2.0052872219658430E-01), HWY_REP4(-3.2460335995264836E-02),
    191       HWY_REP4(3.1546157932479282E-03), HWY_REP4(-1.6692542019380155E-04)};
    192 
    193   const T err = RunApproximation(0, Pi<T>(1) * 2, p, q, EvalPoly(),
    194                                  [](T x) { return static_cast<T>(sin(x)); });
    195   EXPECT_LT(err, sizeof(T) == 8 ? 5E-4 : 7E-4);
    196 }
    197 
    198 void TestLog() {
    199   HWY_ALIGN const T p[4 * (2 + 1)] = {HWY_REP4(-1.8503833400518310E-06),
    200                                       HWY_REP4(1.4287160470083755E+00),
    201                                       HWY_REP4(7.4245873327820566E-01)};
    202   HWY_ALIGN const T q[4 * (2 + 1)] = {HWY_REP4(9.9032814277590719E-01),
    203                                       HWY_REP4(1.0096718572241148E+00),
    204                                       HWY_REP4(1.7409343003366853E-01)};
    205   const T err = RunApproximation(1E-6, 1000, p, q, EvalLog2(), std::log2);
    206   printf("%E\n", err);
    207 }
    208 
    209 HWY_NOINLINE void TestRationalPolynomial() {
    210   TestSimpleGamma();
    211   TestLinearToSrgb8Direct();
    212   TestExp();
    213   TestNegExp();
    214   TestSin();
    215   TestLog();
    216 }
    217 
    218 // NOLINTNEXTLINE(google-readability-namespace-comments)
    219 }  // namespace HWY_NAMESPACE
    220 }  // namespace jxl
    221 HWY_AFTER_NAMESPACE();
    222 
    223 #if HWY_ONCE
    224 namespace jxl {
    225 
    226 class RationalPolynomialTest : public hwy::TestWithParamTarget {};
    227 HWY_TARGET_INSTANTIATE_TEST_SUITE_P(RationalPolynomialTest);
    228 
    229 HWY_EXPORT_AND_TEST_P(RationalPolynomialTest, TestSimpleGamma);
    230 HWY_EXPORT_AND_TEST_P(RationalPolynomialTest, TestLinearToSrgb8Direct);
    231 HWY_EXPORT_AND_TEST_P(RationalPolynomialTest, TestExp);
    232 HWY_EXPORT_AND_TEST_P(RationalPolynomialTest, TestNegExp);
    233 HWY_EXPORT_AND_TEST_P(RationalPolynomialTest, TestSin);
    234 HWY_EXPORT_AND_TEST_P(RationalPolynomialTest, TestLog);
    235 
    236 }  // namespace jxl
    237 #endif  // HWY_ONCE