1  /*
  2   * Copyright (C) 2019, Advanced Micro Devices. All rights reserved.
  3   *
  4   * Redistribution and use in source and binary forms, with or without modification,
  5   * are permitted provided that the following conditions are met:
  6   * 1. Redistributions of source code must retain the above copyright notice,
  7   *    this list of conditions and the following disclaimer.
  8   * 2. Redistributions in binary form must reproduce the above copyright notice,
  9   *    this list of conditions and the following disclaimer in the documentation
 10   *    and/or other materials provided with the distribution.
 11   * 3. Neither the name of the copyright holder nor the names of its contributors
 12   *    may be used to endorse or promote products derived from this software without
 13   *    specific prior written permission.
 14   *
 15   * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
 16   * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
 17   * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
 18   * IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
 19   * INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
 20   * BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA,
 21   * OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,
 22   * WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
 23   * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
 24   * POSSIBILITY OF SUCH DAMAGE.
 25   *
 26   */
 27  #include <stdint.h>
 28  #include <libm_util_amd.h>
 29  #include <libm_special.h>
 30  
 31  #include <libm_macros.h>
 32  #include <libm/amd_funcs_internal.h>
 33  #include <libm/types.h>
 34  #include <libm/typehelper.h>
 35  #include <libm/compiler.h>
 36  #include <libm/poly.h>
 37  
 38  static struct {
 39      float poly[9];
 40      float ln2;
 41  } logf_data = {
 42      .ln2 = 0x1.62e42fefp-1f,
 43      // Polynomial coefficients obtained using fpminimax algorithm from Sollya
 44      .poly = {
 45          0x1.0000000f500p0f,
 46          -0x1.ffffd272p-2f,
 47          0x1.5554e6dap-2f,
 48          -0x1.00150b01p-2f,
 49          0x1.99d8391ap-3f,
 50          -0x1.50224b33p-3f,
 51          0x1.1e779278p-3f,
 52          -0x1.3b751efcp-3f,
 53          0x1.1fde870dp-3f,
 54      },
 55  };
 56  
 57  #define ln2 logf_data.ln2
 58  #define C1 logf_data.poly[0]
 59  #define C2 logf_data.poly[1]
 60  #define C3 logf_data.poly[2]
 61  #define C4 logf_data.poly[3]
 62  #define C5 logf_data.poly[4]
 63  #define C6 logf_data.poly[5]
 64  #define C7 logf_data.poly[6]
 65  #define C8 logf_data.poly[7]
 66  #define C9 logf_data.poly[8]
 67  
 68  
 69  /*
 70   * Implementation Notes:
 71   *     reduce 'x' into the form:
 72   *        x = (-1)^s*2^n*m
 73   * 's' will be always zero, as log is defined for posititve numbers
 74   * 'n' is an integer known as the exponent
 75   * 'm' is mantissa
 76   *
 77   * x is reduced such that the mantissa, m lies in [2/3,4/3]
 78   *      x = 2^n*m where m is in [2/3,4/3]
 79   *      log(x) = log(2^n*m)         We have log(a*b) = log(a)+log(b)
 80   *             = log(2^n) + log(m)   We have log(a^n) = n*log(a)
 81   *             = n*log(2) + log(m)
 82   *             = n*log(2) + log(1+(m-1))
 83   *             = n*log(2) + log(1+f)        Where f = m-1
 84   *             = n*log(2) + log1p(f)        f lies in [-1/3,+1/3]
 85   *
 86   * Thus we have :
 87   * log(x) = n*log(2) + log1p(f)
 88   *            log1p(F) is approximated by using a polynomial
 89   */
 90  
 91  float
 92  FN_PROTOTYPE_FAST(logf)(float x)
 93  {
 94      float m, r, n, f;
 95  
 96      // 0 and FLT_MAX
 97      if ((x > 0.0f) && (x <= 0x1.fffffep+127)) {
 98          // Get the exponent
 99          int32_t ux = (asuint32(x) - 0x3f2aaaab) & 0xff800000;
100          n = (float)(ux >> EXPSHIFTBITS_SP32) ;
101          // Get the mantissa, m is in [2/3, 4/3]
102          m = asfloat(asuint32(x) - ux);
103  
104          f = m - 1.0f;// f is in [-1/3,+1/3]
105  
106          // Compute log1p(f) using Polynomial approximation
107          r = POLY_EVAL_9(f, C1, C2, C3, C4, C5, C6, C7, C8, C9);
108  
109          r = n*ln2 + r;
110  
111      } else {
112          r = x + x;  // Output NaNs for the inputs NaN's
113          if (x  < 0.0f) r =  0.0f / 0.0f; //  NaN
114          if (x == 0.0f) r = -1.0f / 0.0f; // -Inf
115      }
116  
117      return r;
118  
119  }
120  
121  strong_alias (__logf_finite, FN_PROTOTYPE_FAST(logf))
122  weak_alias (logf, FN_PROTOTYPE_FAST(logf))
123  strong_alias (__ieee754_logf, FN_PROTOTYPE_FAST(logf))