summaryrefslogtreecommitdiffstats
path: root/mathfuncs_log.h
blob: cd71eb3c9651628176f7122e553c11bc8e194e2c (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
// -*-C++-*-

#ifndef MATHFUNCS_LOG_H
#define MATHFUNCS_LOG_H

#include "mathfuncs_base.h"

#include <cmath>



namespace vecmathlib {
  
  template<typename realvec_t>
  realvec_t mathfuncs<realvec_t>::vml_log2(realvec_t x)
  {
    // Algorithm inspired by SLEEF 2.80
    
    // Rescale
    intvec_t ilogb_x = ilogb(x * RV(M_SQRT2));
    x = ldexp(x, -ilogb_x);
    VML_ASSERT(all(x >= RV(M_SQRT1_2) && x <= RV(M_SQRT2)));
    
    realvec_t y = (x - RV(1.0)) / (x + RV(1.0));
    realvec_t y2 = y*y;
    
    realvec_t r;
    switch (sizeof(real_t)) {
    case 4:
      // float, error=7.09807175879142775648452461821e-8
      r = RV(0.59723611417135718739797302426);
      r = mad(r, y2, RV(0.961524413175528426101613434));
      r = mad(r, y2, RV(2.88539097665498228703236701));
      break;
    case 8:
#ifdef VML_HAVE_FP_CONTRACT
      // double, error=1.48294180185938512675770096324e-16
      r = RV(0.243683403415639178527756320773);
      r = mad(r, y2, RV(0.26136626803870009948502658));
      r = mad(r, y2, RV(0.320619429891299265439389));
      r = mad(r, y2, RV(0.4121983452028499242926));
      r = mad(r, y2, RV(0.577078017761894161436));
      r = mad(r, y2, RV(0.96179669392233355927));
      r = mad(r, y2, RV(2.8853900817779295236));
#else
      // double, error=2.1410114030383689267772704676e-14
      r = RV(0.283751646449323373643963474845);
      r = mad(r, y2, RV(0.31983138095551191299118812));
      r = mad(r, y2, RV(0.412211603844146279666022));
      r = mad(r, y2, RV(0.5770779098948940070516));
      r = mad(r, y2, RV(0.961796694295973716912));
      r = mad(r, y2, RV(2.885390081777562819196));
#endif
      break;
    default:
      __builtin_unreachable();
    }
    r *= y;
    
    // Undo rescaling
    r += convert_float(ilogb_x);
    
    return r;
  }
  
  
  
  template<typename realvec_t>
  inline
  realvec_t mathfuncs<realvec_t>::vml_log(realvec_t x)
  {
    return log2(x) * RV(M_LN2);
  }

  template<typename realvec_t>
  inline
  realvec_t mathfuncs<realvec_t>::vml_log10(realvec_t x)
  {
    return log(x) * RV(M_LOG10E);
  }

  template<typename realvec_t>
  inline
  realvec_t mathfuncs<realvec_t>::vml_log1p(realvec_t x)
  {
    // TODO: Check SLEEF 2.80 algorithm
    
    return log(RV(1.0) + x);
#if 0
    // Goldberg, theorem 4
    realvec_t x1 = RV(1.0) + x;
    x1.barrier();
    return ifthen(x1 == x, x, x * log(x1) / (x1 - RV(1.0)));
#endif
  }
  
}; // namespace vecmathlib

#endif  // #ifndef MATHFUNCS_LOG_H
OpenPOWER on IntegriCloud