tor-browser

The Tor Browser
git clone https://git.dasho.dev/tor-browser.git
Log | Files | Refs | README | LICENSE

rational_polynomial_test.cc (8506B)

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