MATHLIB User Guide
MATHLIB_log10.cpp
Go to the documentation of this file.
1 /******************************************************************************
2  * Copyright (C) 2023 Texas Instruments Incorporated - https://www.ti.com/
3  *
4  * Redistribution and use in source and binary forms, with or without
5  * modification, are permitted provided that the following conditions
6  * are met:
7  *
8  * Redistributions of source code must retain the above copyright
9  * notice, this list of conditions and the following disclaimer.
10  *
11  * Redistributions in binary form must reproduce the above copyright
12  * notice, this list of conditions and the following disclaimer in the
13  * documentation and/or other materials provided with the
14  * distribution.
15  *
16  * Neither the name of Texas Instruments Incorporated nor the names of
17  * its contributors may be used to endorse or promote products derived
18  * from this software without specific prior written permission.
19  *
20  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
21  * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
22  * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
23  * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
24  * OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
25  * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
26  * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
27  * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
28  * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
29  * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
30  * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
31  *
32  ******************************************************************************/
33 
34 #define ELEMENT_COUNT(x) c7x::element_count_of<x>::value
35 
36 /******************************************************************************/
37 /* */
38 /* Includes */
39 /* */
40 /******************************************************************************/
41 
42 #include "MATHLIB_ilut.h"
43 #include "MATHLIB_log10_scalar.h"
44 #include "MATHLIB_lut.h"
45 #include "MATHLIB_permute.h"
46 #include "MATHLIB_types.h"
47 #include "MATHLIB_utility.h"
48 
49 /******************************************************************************/
50 /* */
51 /* MATHLIB_log10 */
52 /* */
53 /******************************************************************************/
54 
55 // Division function for double precision
56 static inline c7x::double_vec cmn_DIVDP_opt(c7x::double_vec a, c7x::double_vec b)
57 {
58 
59  c7x::double_vec Two = (c7x::double_vec)(2.0f);
60  c7x::double_vec X;
61  X = __recip(b);
62  X = X * (Two - (b * X));
63  X = X * (Two - (b * X));
64  X = X * (Two - (b * X));
65  X = a * X;
66 
67  return X;
68 }
69 
70 // for length greater than 1 and less than 33 use split - 1 function for better optimization
71 static inline void MATHLIB_log10_vector_split1(size_t length, double *restrict pSrc, double *restrict pDst)
72 {
73 
74  // variables
75  size_t numBlocks = 0; // compute loop's iteration count
76  size_t remNumBlocks = 0; // when numBlocks is not a multiple of SIMD width
77 
78  // derive c7x vector type from template typename
79  typedef typename c7x::make_full_vector<double>::type vec;
80 
81  __SE_TEMPLATE_v1 se0Params = __gen_SE_TEMPLATE_v1();
82  __SA_TEMPLATE_v1 sa0Params = __gen_SA_TEMPLATE_v1();
83 
84  MATHLIB_SE0SA01DSequentialInit(&se0Params, &sa0Params, length, pSrc, pDst);
85 
86  // calculate compute loop's iteration counter
87  numBlocks = length / c7x::element_count_of<vec>::value;
88  remNumBlocks = length % c7x::element_count_of<vec>::value;
89  if (remNumBlocks) {
90  numBlocks++;
91  }
92 
93  // open SE0, SE1, and SA0 for reading and writing operands
94  MATHLIB_SE0SA0Open(&se0Params, &sa0Params, pSrc);
95 
96  /***********************************************************************/
97  /* Create and assign values for constants employed on log10 computation */
98  /***********************************************************************/
99 
100  vec Half, MAXe, srHalf, Half_sq, MINe, a0, a1, a2, b0, b1, b2, c1, c2, c10e, W, X, Y, Z, zn, zd, Rz, Sa, Bd, Cn, Da;
101 
102  Half = (vec) 0.5;
103  Half_sq = (vec) 0.5 * 0.5;
104  MAXe = (vec) 1.7976931348623157e+308;
105  srHalf = (vec) 0.70710678118654752440; /* sqrt(0.5) */
106  MINe = (vec) 2.2250738585072014e-308;
107  a0 = (vec) -0.64124943423745581147e+2;
108  a1 = (vec) 0.16383943563021534222e+2;
109  a2 = (vec) -0.78956112887491257267e+0;
110  b0 = (vec) -0.76949932108494879777e+3;
111  b1 = (vec) 0.31203222091924532844e+3;
112  b2 = (vec) -0.35667977739034646171e+2; /* Note b3 = 1.0 */
113  c1 = (vec) 0.693359375; /* 355/512 */
114  c2 = (vec) -2.121944400546905827679e-4;
115  c10e = (vec) 0.43429448190325182765; /* log (base 10) of e */
116 
117  c7x::long_vec long_zero_vec = (c7x::long_vec) 0;
118  vec double_zero_vec = (vec) 0.0;
119  vec outMAX = (vec)(308.254715974092);
120 
121  // compute loop to perform vector log10
122  for (size_t i = 0; i < numBlocks; i++) {
123 
124  vec a = c7x::strm_eng<0, vec>::get_adv();
125 
126  Y = a;
127  c7x::long_vec exp_ = c7x::as_long_vec((c7x::as_ulong_vec(Y) << 1) >> 53);
128 
129  c7x::ulong_vec upper = c7x::as_ulong_vec(Y) & (0x000FFFFF00000000u);
130  upper = 0x3FE0000000000000u | upper;
131 
132  Z = c7x::as_double_vec((0x00000000FFFFFFFFu & c7x::as_ulong_vec(Y)) | upper);
133 
134  __vpred cmp1 = __cmp_eq_pred(exp_, long_zero_vec);
135  Z = __select(cmp1, double_zero_vec, Z);
136 
137  vec z_minus_half = Z - Half;
138  vec z_mul_half = (Z * Half) + Half;
139  __vpred cmp2 = __cmp_lt_pred(srHalf, Z);
140  zn = __select(cmp2, (z_minus_half - Half), z_minus_half);
141  zd = __select(cmp2, z_mul_half, (z_mul_half - Half_sq));
142 
143  X = cmn_DIVDP_opt(zn, zd);
144  W = X * X;
145  Bd = ((((W + b2) * W) + b1) * W) + b0;
146  Cn = (((W * a2) + a1) * W) + a0;
147  Rz = W * cmn_DIVDP_opt(Cn, Bd);
148  Sa = X + (X * Rz);
149 
150  __vpred tmp = c7x::strm_agen<0, vec>::get_vpred();
151  vec * addr = c7x::strm_agen<0, vec>::get_adv(pDst);
152  __vstore_pred(tmp, addr, Sa);
153  }
154 
155  MATHLIB_SE0SA0Close();
156 
157  MATHLIB_SE0SA0Open(&se0Params, &sa0Params, pSrc);
158  __SE1_OPEN(pDst, se0Params);
159 
160  for (size_t i = 0; i < numBlocks; i++) {
161 
162  vec a = c7x::strm_eng<0, vec>::get_adv();
163  Sa = c7x::strm_eng<1, vec>::get_adv();
164 
165  Y = a;
166  c7x::long_vec exp_ = c7x::as_long_vec((c7x::as_ulong_vec(Y) << 1) >> 53);
167  c7x::long_vec N = exp_ - 1022;
168 
169  c7x::ulong_vec upper = c7x::as_ulong_vec(Y) & (0x000FFFFF00000000u);
170  upper = 0x3FE0000000000000u | upper;
171 
172  Z = c7x::as_double_vec((0x00000000FFFFFFFFu & c7x::as_ulong_vec(Y)) | upper);
173 
174  __vpred cmp1 = __cmp_eq_pred(exp_, long_zero_vec);
175  Z = __select(cmp1, double_zero_vec, Z);
176 
177  __vpred cmp2 = __cmp_lt_pred(srHalf, Z);
178 
179  N = __select(cmp2, N, (N - 1));
180 
181  Cn = __low_int_to_double(c7x::as_int_vec(N));
182  Da = ((Cn * c2) + Sa) + (Cn * c1);
183  Da = c10e * Da;
184 
185  /**********************************************************************/
186  /* Bounds checking */
187  /**********************************************************************/
188 
189  // if (Y < MINe) {
190  // Da = -MAXe;
191  // }
192  __vpred cmp_min = __cmp_lt_pred(Y, MINe);
193  Da = __select(cmp_min, -MAXe, Da);
194 
195  // if (Y > MAXe) {
196  // Da = 308.254715974092;
197  // }
198  __vpred cmp_max = __cmp_lt_pred(MAXe, Y);
199  Da = __select(cmp_max, outMAX, Da);
200 
201  __vpred tmp = c7x::strm_agen<0, vec>::get_vpred();
202  vec * addr = c7x::strm_agen<0, vec>::get_adv(pDst);
203  __vstore_pred(tmp, addr, Da);
204  }
205 
206  __SE1_CLOSE();
207  MATHLIB_SE0SA0Close();
208 }
209 
210 // for length greater than 32 use split - 2 function for better optimization
211 static inline void MATHLIB_log10_vector_split2(size_t length, double *restrict pSrc, double *restrict pDst)
212 {
213  // variables
214  size_t numBlocks = 0; // compute loop's iteration count
215  size_t remNumBlocks = 0; // when numBlocks is not a multiple of SIMD width
216 
217  // derive c7x vector type from template typename
218  typedef typename c7x::make_full_vector<double>::type vec;
219 
220  __SE_TEMPLATE_v1 se0Params = __gen_SE_TEMPLATE_v1();
221  __SA_TEMPLATE_v1 sa0Params = __gen_SA_TEMPLATE_v1();
222 
223  MATHLIB_SE0SA01DSequentialInit(&se0Params, &sa0Params, length, pSrc, pDst);
224 
225  // calculate compute loop's iteration counter
226  numBlocks = length / c7x::element_count_of<vec>::value;
227  remNumBlocks = length % c7x::element_count_of<vec>::value;
228  if (remNumBlocks) {
229  numBlocks++;
230  }
231 
232  // open SE0, SE1, and SA0 for reading and writing operands
233  MATHLIB_SE0SA0Open(&se0Params, &sa0Params, pSrc);
234 
235  /***********************************************************************/
236  /* Create and assign values for constants employed on log10 computation */
237  /***********************************************************************/
238 
239  vec Half, MAXe, srHalf, Half_sq, MINe, a0, a1, a2, b0, b1, b2, c1, c2, c10e, W, X, Y, Z, zn, zd, Rz, Sa, Bd, Cn, Da;
240 
241  Half = (vec) 0.5;
242  Half_sq = (vec) 0.5 * 0.5;
243  MAXe = (vec) 1.7976931348623157e+308;
244  srHalf = (vec) 0.70710678118654752440; /* sqrt(0.5) */
245  MINe = (vec) 2.2250738585072014e-308;
246  a0 = (vec) -0.64124943423745581147e+2;
247  a1 = (vec) 0.16383943563021534222e+2;
248  a2 = (vec) -0.78956112887491257267e+0;
249  b0 = (vec) -0.76949932108494879777e+3;
250  b1 = (vec) 0.31203222091924532844e+3;
251  b2 = (vec) -0.35667977739034646171e+2; /* Note b3 = 1.0 */
252  c1 = (vec) 0.693359375; /* 355/512 */
253  c2 = (vec) -2.121944400546905827679e-4;
254  c10e = (vec) 0.43429448190325182765; /* log (base 10) of e */
255 
256  c7x::long_vec long_zero_vec = (c7x::long_vec) 0;
257  vec double_zero_vec = (vec) 0.0;
258  vec outMAX = (vec)(308.254715974092);
259 
260  // compute loop to perform vector log10
261  for (size_t i = 0; i < numBlocks; i++) {
262 
263  vec a = c7x::strm_eng<0, vec>::get_adv();
264 
265  Y = a;
266  c7x::long_vec exp_ = c7x::as_long_vec((c7x::as_ulong_vec(Y) << 1) >> 53);
267 
268  c7x::ulong_vec upper = c7x::as_ulong_vec(Y) & (0x000FFFFF00000000u);
269  upper = 0x3FE0000000000000u | upper;
270 
271  Z = c7x::as_double_vec((0x00000000FFFFFFFFu & c7x::as_ulong_vec(Y)) | upper);
272 
273  __vpred cmp1 = __cmp_eq_pred(exp_, long_zero_vec);
274  Z = __select(cmp1, double_zero_vec, Z);
275 
276  vec z_minus_half = Z - Half;
277  vec z_mul_half = (Z * Half) + Half;
278  __vpred cmp2 = __cmp_lt_pred(srHalf, Z);
279  zn = __select(cmp2, (z_minus_half - Half), z_minus_half);
280  zd = __select(cmp2, z_mul_half, (z_mul_half - Half_sq));
281 
282  X = cmn_DIVDP_opt(zn, zd);
283 
284  __vpred tmp = c7x::strm_agen<0, vec>::get_vpred();
285  vec * addr = c7x::strm_agen<0, vec>::get_adv(pDst);
286  __vstore_pred(tmp, addr, X);
287  }
288 
289  MATHLIB_SE0SA0Close();
290 
291  MATHLIB_SE0SA0Open(&se0Params, &sa0Params, pDst);
292 
293  for (size_t i = 0; i < numBlocks; i++) {
294 
295  X = c7x::strm_eng<0, vec>::get_adv();
296 
297  W = X * X;
298  Bd = ((((W + b2) * W) + b1) * W) + b0;
299  Cn = (((W * a2) + a1) * W) + a0;
300  Rz = W * cmn_DIVDP_opt(Cn, Bd);
301  Sa = X + (X * Rz);
302 
303  __vpred tmp = c7x::strm_agen<0, vec>::get_vpred();
304  vec * addr = c7x::strm_agen<0, vec>::get_adv(pDst);
305  __vstore_pred(tmp, addr, Sa);
306  }
307 
308  MATHLIB_SE0SA0Close();
309 
310  MATHLIB_SE0SA0Open(&se0Params, &sa0Params, pSrc);
311 
312  __SE1_OPEN(pDst, se0Params);
313 
314  for (size_t i = 0; i < numBlocks; i++) {
315 
316  vec a = c7x::strm_eng<0, vec>::get_adv();
317  Sa = c7x::strm_eng<1, vec>::get_adv();
318 
319  Y = a;
320  c7x::long_vec exp_ = c7x::as_long_vec((c7x::as_ulong_vec(Y) << 1) >> 53);
321  c7x::long_vec N = exp_ - 1022;
322 
323  c7x::ulong_vec upper = c7x::as_ulong_vec(Y) & (0x000FFFFF00000000u);
324  upper = 0x3FE0000000000000u | upper;
325 
326  Z = c7x::as_double_vec((0x00000000FFFFFFFFu & c7x::as_ulong_vec(Y)) | upper);
327 
328  __vpred cmp1 = __cmp_eq_pred(exp_, long_zero_vec);
329  Z = __select(cmp1, double_zero_vec, Z);
330 
331  __vpred cmp2 = __cmp_lt_pred(srHalf, Z);
332 
333  N = __select(cmp2, N, (N - 1));
334 
335  Cn = __low_int_to_double(c7x::as_int_vec(N));
336  Da = ((Cn * c2) + Sa) + (Cn * c1);
337  Da = c10e * Da;
338 
339  /**********************************************************************/
340  /* Bounds checking */
341  /**********************************************************************/
342 
343  // if (Y < MINe) {
344  // Da = -MAXe;
345  // }
346  __vpred cmp_min = __cmp_lt_pred(Y, MINe);
347  Da = __select(cmp_min, -MAXe, Da);
348 
349  // if (Y > MAXe) {
350  // Da = 308.254715974092;
351  // }
352  __vpred cmp_max = __cmp_lt_pred(MAXe, Y);
353  Da = __select(cmp_max, outMAX, Da);
354 
355  __vpred tmp = c7x::strm_agen<0, vec>::get_vpred();
356  vec * addr = c7x::strm_agen<0, vec>::get_adv(pDst);
357  __vstore_pred(tmp, addr, Da);
358  }
359 
360  __SE1_CLOSE();
361 
362  MATHLIB_SE0SA0Close();
363 }
364 
365 
366 template <typename T> static inline void MATHLIB_log10_vector(size_t length, T *pSrc, T *pDst);
367 
368 // this method performs log10 computation of input vector
369 template <> inline void MATHLIB_log10_vector<float>(size_t length, float *restrict pSrc, float *restrict pDst)
370 {
371  // variables
372  size_t numBlocks = 0; // compute loop's iteration count
373  size_t remNumBlocks = 0; // when numBlocks is not a multiple of SIMD width
374 
375  // derive c7x vector type from template typename
376  typedef typename c7x::make_full_vector<float>::type vec;
377 
378  __SE_TEMPLATE_v1 se0Params = __gen_SE_TEMPLATE_v1();
379  __SA_TEMPLATE_v1 sa0Params = __gen_SA_TEMPLATE_v1();
380 
381  MATHLIB_SE0SA01DSequentialInit(&se0Params, &sa0Params, length, pSrc, pDst);
382 
383  // calculate compute loop's iteration counter
384  numBlocks = length / c7x::element_count_of<vec>::value;
385  remNumBlocks = length % c7x::element_count_of<vec>::value;
386  if (remNumBlocks) {
387  numBlocks++;
388  }
389 
390  // open SE0, SE1, and SA0 for reading and writing operands
391  MATHLIB_SE0SA0Open(&se0Params, &sa0Params, pSrc);
392 
393  /***********************************************************************/
394  /* Create and assign values for constants employed on log10 computation */
395  /***********************************************************************/
396 
397  vec C1, C2, C3, C4, C5, eMax, outVecMin, outVecMax;
398  c7x::double_vec ln2, base;
399  c7x::uint_vec zero;
400  zero = (c7x::uint_vec) 0;
401 
402  ln2 = (c7x::double_vec) 0.693147180559945;
403  base = (c7x::double_vec) 0.4342944819033f;
404  C1 = (vec) -0.2302894f;
405  C2 = (vec) 0.1908169f;
406  C3 = (vec) -0.2505905f;
407  C4 = (vec) 0.3333164f;
408  C5 = (vec) -0.5000002f;
409  eMax = (vec) 3.402823466e+38f;
410  outVecMin = (vec) 0xFF800000u;
411  outVecMax = (vec) 308.2547f;
412 
413  // compute loop to perform vector log10
414  for (size_t i = 0; i < numBlocks; i++) {
415  vec inVec = c7x::strm_eng<0, vec>::get_adv();
416 
417  /**********************************************************************/
418  /* Create variables employed on log10 computation */
419  /**********************************************************************/
420  vec pol, r1, r2, r3, r4;
421  c7x::double_vec inVecVals_odd, inVecVals_even, inVecVals_oddReciprocal, inVecVals_evenReciprocal,
422  inVecReciprocalApprox_8_15, inVecReciprocalApprox_0_7, inVecVals_8_15, inVecVals_0_7, rVals_0_7, rVals_8_15,
423  TVals_8_15, TVals_0_7, NVals_odd, NVals_even, NVals_0_7, NVals_8_15, outVec_8_15, outVec_0_7;
424  c7x::uint_vec inVecReciprocal_32_63, inVecReciprocalClr_32_63, inVecReciprocalApprox_32_63, indexT;
425  c7x::int_vec N;
426  vec outVec;
427 
428  /**********************************************************************/
429  /* Calculate Taylor series approximation for log10 */
430  /**********************************************************************/
431 
432  // Split vectors to compute r with double precision
433  inVecVals_odd = __high_float_to_double(inVec);
434  inVecVals_even = __low_float_to_double(inVec);
435  inVecVals_oddReciprocal = __recip(inVecVals_odd);
436  inVecVals_evenReciprocal = __recip(inVecVals_even);
437 
438  // Create floating point reciprocal approximation
439  // Upper 32 bits of all inVec reciprocal values
440  inVecReciprocal_32_63 = c7x::reinterpret<c7x::uint_vec>(
441  __permute_odd_odd_int(MATHLIB_vperm_data_interweave_0_63, c7x::as_uchar_vec(inVecVals_oddReciprocal),
442  c7x::as_uchar_vec(inVecVals_evenReciprocal)));
443 
444  // Clear bits 0-16 inclusive
445  inVecReciprocalClr_32_63 = inVecReciprocal_32_63 & 0xFFFE0000u;
446 
447  // Concatenate cleared bit reciprocal with zero bits
448  inVecReciprocalApprox_8_15 = c7x::reinterpret<c7x::double_vec>(__permute_high_high(
449  MATHLIB_vperm_data_interweave_0_63, c7x::as_uchar_vec(inVecReciprocalClr_32_63), c7x::as_uchar_vec(zero)));
450  inVecReciprocalApprox_0_7 = c7x::reinterpret<c7x::double_vec>(__permute_low_low(
451  MATHLIB_vperm_data_interweave_0_63, c7x::as_uchar_vec(inVecReciprocalClr_32_63), c7x::as_uchar_vec(zero)));
452 
453  // Split inVec into two vectors with double precision
454  inVecVals_0_7 = c7x::reinterpret<c7x::double_vec>(__permute_low_low(
455  MATHLIB_vperm_data_dp_interweave_0_63, c7x::as_uchar_vec(inVecVals_odd), c7x::as_uchar_vec(inVecVals_even)));
456  inVecVals_8_15 = c7x::reinterpret<c7x::double_vec>(__permute_high_high(
457  MATHLIB_vperm_data_dp_interweave_0_63, c7x::as_uchar_vec(inVecVals_odd), c7x::as_uchar_vec(inVecVals_even)));
458 
459  // Calculate r in double precision
460  rVals_0_7 = (inVecReciprocalApprox_0_7 * inVecVals_0_7) - 1.0;
461  rVals_8_15 = (inVecReciprocalApprox_8_15 * inVecVals_8_15) - 1.0;
462 
463  // Convert r to float, compute r to the power of 2, 3, 4
464  r1 = c7x::reinterpret<vec>(__permute_even_even_int(MATHLIB_vperm_data_0_63,
465  c7x::as_uchar_vec(__double_to_float(rVals_8_15)),
466  c7x::as_uchar_vec(__double_to_float(rVals_0_7))));
467  r2 = r1 * r1;
468  r3 = r1 * r2;
469  r4 = r2 * r2;
470 
471  // Compute Taylor series polynomial
472  pol = (C5 * r2) + ((C4 * r3) + ((((C2 * r1) + C3) + (C1 * r2)) * r4));
473  pol = pol * __double_to_float(base);
474 
475  /**********************************************************************/
476  /* Calculate N */
477  /**********************************************************************/
478 
479  // Upper 32 bits of all inVec reciprocal approximation values
480  inVecReciprocalApprox_32_63 = c7x::reinterpret<c7x::uint_vec>(
481  __permute_odd_odd_int(MATHLIB_vperm_data_0_63, c7x::as_uchar_vec(inVecReciprocalApprox_8_15),
482  c7x::as_uchar_vec(inVecReciprocalApprox_0_7)));
483 
484  N = c7x::convert<c7x::int_vec>(((inVecReciprocalApprox_32_63 << 1) >> 21) - 1023);
485 
486  // Covert N to double precision for later calculation with LUT values
487  NVals_odd = __high_int_to_double(N);
488  NVals_even = __low_int_to_double(N);
489  NVals_0_7 = c7x::reinterpret<c7x::double_vec>(__permute_low_low(
490  MATHLIB_vperm_data_dp_interweave_0_63, c7x::as_uchar_vec(NVals_odd), c7x::as_uchar_vec(NVals_even)));
491  NVals_8_15 = c7x::reinterpret<c7x::double_vec>(__permute_high_high(
492  MATHLIB_vperm_data_dp_interweave_0_63, c7x::as_uchar_vec(NVals_odd), c7x::as_uchar_vec(NVals_even)));
493 
494  /**********************************************************************/
495  /* Determine LUT values */
496  /**********************************************************************/
497 
498  // Calculate LUT index
499  indexT = (((inVecReciprocalApprox_32_63 << 12) >> 29) + MATHLIB_LOGTABLE_OFFSET);
500 
501  // Read from LUT or ILUT and reconstruct double values split into two vectors
502  c7x::uint_vec upperBitsIndexT = MATHLIB_LUTReadUpperBits(indexT);
503  c7x::uint_vec lowerBitsIndexT = MATHLIB_LUTReadLowerBits(indexT);
504 
505  // Read from LUT and reconstruct double values split into two vectors
506  TVals_8_15 = c7x::reinterpret<c7x::double_vec>(__permute_high_high(
507  MATHLIB_vperm_data_interweave_0_63, c7x::as_uchar_vec(upperBitsIndexT), c7x::as_uchar_vec(lowerBitsIndexT)));
508  TVals_0_7 = c7x::reinterpret<c7x::double_vec>(__permute_low_low(
509  MATHLIB_vperm_data_interweave_0_63, c7x::as_uchar_vec(upperBitsIndexT), c7x::as_uchar_vec(lowerBitsIndexT)));
510 
511  // Calculate an adjusted T
512  TVals_8_15 = (TVals_8_15 - (ln2 * NVals_8_15)) * base;
513  TVals_0_7 = (TVals_0_7 - (ln2 * NVals_0_7)) * base;
514 
515  /**********************************************************************/
516  /* Calculate output with adjusted LUT and Taylor series values */
517  /**********************************************************************/
518 
519  // Add LUT values and Taylor series values
520  // TODO: Multiple by base
521  outVec_0_7 = (rVals_0_7 * base) + TVals_0_7;
522  outVec_8_15 = (rVals_8_15 * base) + TVals_8_15;
523 
524  // Combine output vector into one floating point result
525  outVec = c7x::reinterpret<vec>(__permute_even_even_int(MATHLIB_vperm_data_0_63,
526  c7x::as_uchar_vec(__double_to_float(outVec_8_15)),
527  c7x::as_uchar_vec(__double_to_float(outVec_0_7))));
528  outVec = outVec + pol;
529 
530  /**********************************************************************/
531  /* Bounds checking */
532  /**********************************************************************/
533 
534  // if (a <= 0.0f) {
535  // res = _itof(0xFF800000u);
536  // }
537  __vpred cmp_min = __cmp_le_pred(inVec, c7x::convert<vec>(zero));
538  outVec = __select(cmp_min, outVecMin, outVec);
539 
540  // if (a > MAXe) {
541  // res = 709.7827f;
542  // }
543  __vpred cmp_max = __cmp_lt_pred(eMax, inVec);
544  outVec = __select(cmp_max, outVecMax, outVec);
545 
546  __vpred tmp = c7x::strm_agen<0, vec>::get_vpred();
547  vec *addr = c7x::strm_agen<0, vec>::get_adv(pDst);
548  __vstore_pred(tmp, addr, outVec);
549  }
550 
551  MATHLIB_SE0SA0Close();
552 }
553 
554 
555 // this method performs log10 computation of input vector
556 template <> inline void MATHLIB_log10_vector<double>(size_t length, double *restrict pSrc, double *restrict pDst)
557 {
558  if (length < 33) {
559  MATHLIB_log10_vector_split1(length, pSrc, pDst);
560  }
561  else {
562  MATHLIB_log10_vector_split2(length, pSrc, pDst);
563  }
564 }
565 
566 // this method performs exponential computation of input vector
567 template <typename T> MATHLIB_STATUS MATHLIB_log10(size_t length, T *restrict pSrc, T *restrict pDst)
568 {
569  MATHLIB_STATUS status = MATHLIB_SUCCESS; // return function status
570 
571  // check for null pointers and non-zero length
572  status = MATHLIB_checkParams(length, pSrc, pDst);
573 
574  // for length equal to 1 use scalar operations othewise use vector operations for better optimization
575  if (status == MATHLIB_SUCCESS) {
576  if (length < 2) {
577  pDst[0] = MATHLIB_log10_scalar_ci<T>(pSrc[0]);
578  }
579  else {
580  MATHLIB_log10_vector<T>(length, pSrc, pDst);
581  }
582  }
583  return status;
584 }
585 
586 /******************************************************************************/
587 /* */
588 /* Explicit templatization for datatypes supported by MATHLIB_log10 */
589 /* */
590 /******************************************************************************/
591 
592 template MATHLIB_STATUS MATHLIB_log10<float>(size_t length, float *pSrc, float *pDst);
593 template MATHLIB_STATUS MATHLIB_log10<double>(size_t length, double *pSrc, double *pDst);
594 
595 /******************************************************************************/
596 /* */
597 /* C-interface wrapper functions */
598 /* */
599 /******************************************************************************/
600 
601 extern "C" {
602 
603 // single-precision wrapper
604 MATHLIB_STATUS MATHLIB_log10_sp(size_t length, float *pSrc, float *pDst)
605 {
606  MATHLIB_STATUS status = MATHLIB_log10<float>(length, pSrc, pDst);
607  return status;
608 }
609 
610 // double-precision wrapper
611 MATHLIB_STATUS MATHLIB_log10_dp(size_t length, double *pSrc, double *pDst)
612 {
613  MATHLIB_STATUS status = MATHLIB_log10<double>(length, pSrc, pDst);
614  return status;
615 }
616 
617 
618 } // extern "C"
MATHLIB_ilut.h
cmn_DIVDP_opt
static c7x::double_vec cmn_DIVDP_opt(c7x::double_vec a, c7x::double_vec b)
Definition: MATHLIB_log10.cpp:56
MATHLIB_log10_vector< double >
void MATHLIB_log10_vector< double >(size_t length, double *restrict pSrc, double *restrict pDst)
Definition: MATHLIB_log10.cpp:556
MATHLIB_log10_vector_split2
static void MATHLIB_log10_vector_split2(size_t length, double *restrict pSrc, double *restrict pDst)
Definition: MATHLIB_log10.cpp:211
MATHLIB_log10_vector< float >
void MATHLIB_log10_vector< float >(size_t length, float *restrict pSrc, float *restrict pDst)
Definition: MATHLIB_log10.cpp:369
MATHLIB_log10
MATHLIB_STATUS MATHLIB_log10(size_t length, T *restrict pSrc, T *restrict pDst)
Definition: MATHLIB_log10.cpp:567
MATHLIB_log10< float >
template MATHLIB_STATUS MATHLIB_log10< float >(size_t length, float *pSrc, float *pDst)
MATHLIB_log10_vector_split1
static void MATHLIB_log10_vector_split1(size_t length, double *restrict pSrc, double *restrict pDst)
Definition: MATHLIB_log10.cpp:71
MATHLIB_log10_vector
static void MATHLIB_log10_vector(size_t length, T *pSrc, T *pDst)
MATHLIB_log10< double >
template MATHLIB_STATUS MATHLIB_log10< double >(size_t length, double *pSrc, double *pDst)
MATHLIB_log10_scalar.h
MATHLIB_lut.h
MATHLIB_permute.h
MATHLIB_types.h
MATHLIB_utility.h
MATHLIB_log10_dp
MATHLIB_STATUS MATHLIB_log10_dp(size_t length, double *pSrc, double *pDst)
This function is the C interface for MATHLIB_log10. Function accepts double pointers.
Definition: MATHLIB_log10.cpp:611
MATHLIB_log10_sp
MATHLIB_STATUS MATHLIB_log10_sp(size_t length, float *pSrc, float *pDst)
This function is the C interface for MATHLIB_log10. Function accepts float pointers.
Definition: MATHLIB_log10.cpp:604
MATHLIB_LOGTABLE_OFFSET
#define MATHLIB_LOGTABLE_OFFSET
Definition: MATHLIB_lut.h:66
MATHLIB_LUTReadLowerBits
static c7x::uint_vec MATHLIB_LUTReadLowerBits(vecType vecOffset)
This method reads bits 31-0 of LUT value at vecOffset.
Definition: MATHLIB_lut.h:111
MATHLIB_LUTReadUpperBits
static c7x::uint_vec MATHLIB_LUTReadUpperBits(vecType vecOffset)
This method reads bits 63-32 of LUT value at vecOffset.
Definition: MATHLIB_lut.h:86
MATHLIB_SE0SA0Close
static void MATHLIB_SE0SA0Close()
This method performs SE0 and SA0 close.
Definition: MATHLIB_utility.h:166
MATHLIB_SE0SA01DSequentialInit
static void MATHLIB_SE0SA01DSequentialInit(__SE_TEMPLATE_v1 *se0Params, __SA_TEMPLATE_v1 *sa0Params, size_t length, T *pSrc, T *pDst)
Definition: MATHLIB_utility.h:66
MATHLIB_checkParams
static MATHLIB_STATUS MATHLIB_checkParams(size_t length, T *pSrc, T *pDst)
This method performs parameter checks for MATHLIB function.
Definition: MATHLIB_utility.h:214
MATHLIB_SE0SA0Open
static void MATHLIB_SE0SA0Open(__SE_TEMPLATE_v1 *se0Params, __SA_TEMPLATE_v1 *sa0Params, T *pSrc)
This method performs SE0 and SA0 open.
Definition: MATHLIB_utility.h:117
MATHLIB_STATUS_NAME
MATHLIB_STATUS_NAME
The enumeration of all status codes.
Definition: MATHLIB_types.h:128
MATHLIB_SUCCESS
@ MATHLIB_SUCCESS
Definition: MATHLIB_types.h:129