2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2012,2013,2014,2015,2017 by the GROMACS development team.
5 * Copyright (c) 2018,2019,2020,2021, by the GROMACS development team, led by
6 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
7 * and including many others, as listed in the AUTHORS file in the
8 * top-level source directory and at http://www.gromacs.org.
10 * GROMACS is free software; you can redistribute it and/or
11 * modify it under the terms of the GNU Lesser General Public License
12 * as published by the Free Software Foundation; either version 2.1
13 * of the License, or (at your option) any later version.
15 * GROMACS is distributed in the hope that it will be useful,
16 * but WITHOUT ANY WARRANTY; without even the implied warranty of
17 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18 * Lesser General Public License for more details.
20 * You should have received a copy of the GNU Lesser General Public
21 * License along with GROMACS; if not, see
22 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
23 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
25 * If you want to redistribute modifications to GROMACS, please
26 * consider that scientific software is very special. Version
27 * control is crucial - bugs must be traceable. We will be happy to
28 * consider code for inclusion in the official distribution, but
29 * derived work must not be called official GROMACS. Details are found
30 * in the README & COPYING files - if they are missing, get the
31 * official version at http://www.gromacs.org.
33 * To help us fund GROMACS development, we humbly ask that you cite
34 * the research papers on the package. Check out http://www.gromacs.org.
36 #ifndef GMX_SIMD_SIMD_MATH_H
37 #define GMX_SIMD_SIMD_MATH_H
39 /*! \libinternal \file
41 * \brief Math functions for SIMD datatypes.
43 * \attention This file is generic for all SIMD architectures, so you cannot
44 * assume that any of the optional SIMD features (as defined in simd.h) are
45 * present. In particular, this means you cannot assume support for integers,
46 * logical operations (neither on floating-point nor integer values), shifts,
47 * and the architecture might only have SIMD for either float or double.
48 * Second, to keep this file clean and general, any additions to this file
49 * must work for all possible SIMD architectures in both single and double
50 * precision (if they support it), and you cannot make any assumptions about
53 * \author Erik Lindahl <erik.lindahl@scilifelab.se>
56 * \ingroup module_simd
65 #include "gromacs/math/units.h"
66 #include "gromacs/math/utilities.h"
67 #include "gromacs/simd/simd.h"
68 #include "gromacs/utility/basedefinitions.h"
69 #include "gromacs/utility/real.h"
77 /*! \addtogroup module_simd */
80 /*! \name Implementation accuracy settings
86 # if GMX_SIMD_HAVE_FLOAT
88 /*! \name Single precision SIMD math functions
90 * \note In most cases you should use the real-precision functions instead.
94 /****************************************
95 * SINGLE PRECISION SIMD MATH FUNCTIONS *
96 ****************************************/
98 # if !GMX_SIMD_HAVE_NATIVE_COPYSIGN_FLOAT
99 /*! \brief Composes floating point value with the magnitude of x and the sign of y.
101 * \param x Values to set sign for
102 * \param y Values used to set sign
103 * \return Magnitude of x, sign of y
105 static inline SimdFloat gmx_simdcall copysign(SimdFloat x, SimdFloat y)
107 # if GMX_SIMD_HAVE_LOGICAL
108 return abs(x) | (SimdFloat(GMX_FLOAT_NEGZERO) & y);
110 return blend(abs(x), -abs(x), y < setZero());
115 # if !GMX_SIMD_HAVE_NATIVE_RSQRT_ITER_FLOAT
116 /*! \brief Perform one Newton-Raphson iteration to improve 1/sqrt(x) for SIMD float.
118 * This is a low-level routine that should only be used by SIMD math routine
119 * that evaluates the inverse square root.
121 * \param lu Approximation of 1/sqrt(x), typically obtained from lookup.
122 * \param x The reference (starting) value x for which we want 1/sqrt(x).
123 * \return An improved approximation with roughly twice as many bits of accuracy.
125 static inline SimdFloat gmx_simdcall rsqrtIter(SimdFloat lu, SimdFloat x)
127 SimdFloat tmp1 = x * lu;
128 SimdFloat tmp2 = SimdFloat(-0.5F) * lu;
129 tmp1 = fma(tmp1, lu, SimdFloat(-3.0F));
134 /*! \brief Calculate 1/sqrt(x) for SIMD float.
136 * \param x Argument that must be larger than GMX_FLOAT_MIN and smaller than
137 * GMX_FLOAT_MAX, i.e. within the range of single precision.
138 * For the single precision implementation this is obviously always
139 * true for positive values, but for double precision it adds an
140 * extra restriction since the first lookup step might have to be
141 * performed in single precision on some architectures. Note that the
142 * responsibility for checking falls on you - this routine does not
145 * \return 1/sqrt(x). Result is undefined if your argument was invalid.
147 static inline SimdFloat gmx_simdcall invsqrt(SimdFloat x)
149 SimdFloat lu = rsqrt(x);
150 # if (GMX_SIMD_RSQRT_BITS < GMX_SIMD_ACCURACY_BITS_SINGLE)
151 lu = rsqrtIter(lu, x);
153 # if (GMX_SIMD_RSQRT_BITS * 2 < GMX_SIMD_ACCURACY_BITS_SINGLE)
154 lu = rsqrtIter(lu, x);
156 # if (GMX_SIMD_RSQRT_BITS * 4 < GMX_SIMD_ACCURACY_BITS_SINGLE)
157 lu = rsqrtIter(lu, x);
162 /*! \brief Calculate 1/sqrt(x) for two SIMD floats.
164 * \param x0 First set of arguments, x0 must be in single range (see below).
165 * \param x1 Second set of arguments, x1 must be in single range (see below).
166 * \param[out] out0 Result 1/sqrt(x0)
167 * \param[out] out1 Result 1/sqrt(x1)
169 * In particular for double precision we can sometimes calculate square root
170 * pairs slightly faster by using single precision until the very last step.
172 * \note Both arguments must be larger than GMX_FLOAT_MIN and smaller than
173 * GMX_FLOAT_MAX, i.e. within the range of single precision.
174 * For the single precision implementation this is obviously always
175 * true for positive values, but for double precision it adds an
176 * extra restriction since the first lookup step might have to be
177 * performed in single precision on some architectures. Note that the
178 * responsibility for checking falls on you - this routine does not
181 static inline void gmx_simdcall invsqrtPair(SimdFloat x0, SimdFloat x1, SimdFloat* out0, SimdFloat* out1)
187 # if !GMX_SIMD_HAVE_NATIVE_RCP_ITER_FLOAT
188 /*! \brief Perform one Newton-Raphson iteration to improve 1/x for SIMD float.
190 * This is a low-level routine that should only be used by SIMD math routine
191 * that evaluates the reciprocal.
193 * \param lu Approximation of 1/x, typically obtained from lookup.
194 * \param x The reference (starting) value x for which we want 1/x.
195 * \return An improved approximation with roughly twice as many bits of accuracy.
197 static inline SimdFloat gmx_simdcall rcpIter(SimdFloat lu, SimdFloat x)
199 return lu * fnma(lu, x, SimdFloat(2.0F));
203 /*! \brief Calculate 1/x for SIMD float.
205 * \param x Argument with magnitude larger than GMX_FLOAT_MIN and smaller than
206 * GMX_FLOAT_MAX, i.e. within the range of single precision.
207 * For the single precision implementation this is obviously always
208 * true for positive values, but for double precision it adds an
209 * extra restriction since the first lookup step might have to be
210 * performed in single precision on some architectures. Note that the
211 * responsibility for checking falls on you - this routine does not
214 * \return 1/x. Result is undefined if your argument was invalid.
216 static inline SimdFloat gmx_simdcall inv(SimdFloat x)
218 SimdFloat lu = rcp(x);
219 # if (GMX_SIMD_RCP_BITS < GMX_SIMD_ACCURACY_BITS_SINGLE)
222 # if (GMX_SIMD_RCP_BITS * 2 < GMX_SIMD_ACCURACY_BITS_SINGLE)
225 # if (GMX_SIMD_RCP_BITS * 4 < GMX_SIMD_ACCURACY_BITS_SINGLE)
231 /*! \brief Division for SIMD floats
233 * \param nom Nominator
234 * \param denom Denominator, with magnitude in range (GMX_FLOAT_MIN,GMX_FLOAT_MAX).
235 * For single precision this is equivalent to a nonzero argument,
236 * but in double precision it adds an extra restriction since
237 * the first lookup step might have to be performed in single
238 * precision on some architectures. Note that the responsibility
239 * for checking falls on you - this routine does not check arguments.
243 * \note This function does not use any masking to avoid problems with
244 * zero values in the denominator.
246 static inline SimdFloat gmx_simdcall operator/(SimdFloat nom, SimdFloat denom)
248 return nom * inv(denom);
251 /*! \brief Calculate 1/sqrt(x) for masked entries of SIMD float.
253 * This routine only evaluates 1/sqrt(x) for elements for which mask is true.
254 * Illegal values in the masked-out elements will not lead to
255 * floating-point exceptions.
257 * \param x Argument that must be larger than GMX_FLOAT_MIN and smaller than
258 * GMX_FLOAT_MAX for masked-in entries.
259 * See \ref invsqrt for the discussion about argument restrictions.
261 * \return 1/sqrt(x). Result is undefined if your argument was invalid or
262 * entry was not masked, and 0.0 for masked-out entries.
264 static inline SimdFloat maskzInvsqrt(SimdFloat x, SimdFBool m)
266 SimdFloat lu = maskzRsqrt(x, m);
267 # if (GMX_SIMD_RSQRT_BITS < GMX_SIMD_ACCURACY_BITS_SINGLE)
268 lu = rsqrtIter(lu, x);
270 # if (GMX_SIMD_RSQRT_BITS * 2 < GMX_SIMD_ACCURACY_BITS_SINGLE)
271 lu = rsqrtIter(lu, x);
273 # if (GMX_SIMD_RSQRT_BITS * 4 < GMX_SIMD_ACCURACY_BITS_SINGLE)
274 lu = rsqrtIter(lu, x);
279 /*! \brief Calculate 1/x for SIMD float, masked version.
281 * \param x Argument with magnitude larger than GMX_FLOAT_MIN and smaller than
282 * GMX_FLOAT_MAX for masked-in entries.
283 * See \ref invsqrt for the discussion about argument restrictions.
285 * \return 1/x for elements where m is true, or 0.0 for masked-out entries.
287 static inline SimdFloat gmx_simdcall maskzInv(SimdFloat x, SimdFBool m)
289 SimdFloat lu = maskzRcp(x, m);
290 # if (GMX_SIMD_RCP_BITS < GMX_SIMD_ACCURACY_BITS_SINGLE)
293 # if (GMX_SIMD_RCP_BITS * 2 < GMX_SIMD_ACCURACY_BITS_SINGLE)
296 # if (GMX_SIMD_RCP_BITS * 4 < GMX_SIMD_ACCURACY_BITS_SINGLE)
302 /*! \brief Calculate sqrt(x) for SIMD floats
304 * \tparam opt By default, this function checks if the input value is 0.0
305 * and masks this to return the correct result. If you are certain
306 * your argument will never be zero, and you know you need to save
307 * every single cycle you can, you can alternatively call the
308 * function as sqrt<MathOptimization::Unsafe>(x).
310 * \param x Argument that must be in range 0 <=x <= GMX_FLOAT_MAX, since the
311 * lookup step often has to be implemented in single precision.
312 * Arguments smaller than GMX_FLOAT_MIN will always lead to a zero
313 * result, even in double precision. If you are using the unsafe
314 * math optimization parameter, the argument must be in the range
315 * GMX_FLOAT_MIN <= x <= GMX_FLOAT_MAX.
317 * \return sqrt(x). The result is undefined if the input value does not fall
318 * in the allowed range specified for the argument.
320 template<MathOptimization opt = MathOptimization::Safe>
321 static inline SimdFloat gmx_simdcall sqrt(SimdFloat x)
323 if (opt == MathOptimization::Safe)
325 SimdFloat res = maskzInvsqrt(x, setZero() < x);
330 return x * invsqrt(x);
334 /*! \brief Cube root for SIMD floats
336 * \param x Argument to calculate cube root of. Can be negative or zero,
337 * but NaN or Inf values are not supported. Denormal values will
339 * \return Cube root of x.
341 static inline SimdFloat gmx_simdcall cbrt(SimdFloat x)
343 const SimdFloat signBit(GMX_FLOAT_NEGZERO);
344 const SimdFloat minFloat(std::numeric_limits<float>::min());
345 // Bias is 128-1 = 127, which is not divisible by 3. Since the largest-magnitude
346 // negative exponent from frexp() is -126, we can subtract one more unit to get 126
347 // as offset, which is divisible by 3 (result 42). To avoid clang warnings about fragile integer
348 // division mixed with FP, we let the divided value (42) be the original constant.
349 const std::int32_t offsetDiv3(42);
350 const SimdFloat c2(-0.191502161678719066F);
351 const SimdFloat c1(0.697570460207922770F);
352 const SimdFloat c0(0.492659620528969547F);
353 const SimdFloat one(1.0F);
354 const SimdFloat two(2.0F);
355 const SimdFloat three(3.0F);
356 const SimdFloat oneThird(1.0F / 3.0F);
357 const SimdFloat cbrt2(1.2599210498948731648F);
358 const SimdFloat sqrCbrt2(1.5874010519681994748F);
360 // To calculate cbrt(x) we first take the absolute value of x but save the sign,
361 // since cbrt(-x) = -cbrt(x). Then we only need to consider positive values for
363 // A number x is represented in IEEE754 as fraction*2^e. We rewrite this as
364 // x=fraction*2^(3*n)*2^m, where e=3*n+m, and m is a remainder.
365 // The cube root can the be evaluated by calculating the cube root of the fraction
366 // limited to the mantissa range, multiplied by 2^mod (which is either 1, +/-2^(1/3) or
367 // +/-2^(2/3), and then we load this into a new IEEE754 fp number with the exponent 2^n, where
368 // n is the integer part of the original exponent divided by 3.
370 SimdFloat xSignBit = x & signBit; // create bit mask where the sign bit is 1 for x elements < 0
371 SimdFloat xAbs = andNot(signBit, x); // select everthing but the sign bit => abs(x)
372 SimdFBool xIsNonZero = (minFloat <= xAbs); // treat denormals as 0
375 SimdFloat y = frexp(xAbs, &exponent);
376 // For the mantissa (y) we will use a limited-range approximation of cbrt(y),
377 // by first using a polynomial and then evaluating
378 // Transform y to z = c2*y^2 + c1*y + c0, then w = z^3, and finally
379 // evaluate the quotient q = z * (w + 2 * y) / (2 * w + y).
380 SimdFloat z = fma(fma(y, c2, c1), y, c0);
381 SimdFloat w = z * z * z;
382 SimdFloat nom = z * fma(two, y, w);
383 SimdFloat invDenom = inv(fma(two, w, y));
385 // Handle the exponent. In principle there are beautiful ways to do this with custom 16-bit
386 // division converted to multiplication... but we can't do that since our SIMD layer cannot
387 // assume the presence of integer shift operations!
388 // However, when I first worked with the integer algorithm I still came up with a neat
389 // optimization, so I'll describe the full algorithm here in case we ever want to use it
392 // Our dividend is signed, which is a complication, but let's consider the unsigned case
393 // first: Division by 3 corresponds to multiplication by 1010101... Since we also know
394 // our dividend is less than 16 bits (exponent range) we can accomplish this by
395 // multiplying with 21845 (which is almost 2^16/3 - 21845.333 would be exact) and then
396 // right-shifting by 16 bits to divide out the 2^16 part.
397 // If we add 1 to the dividend to handle the extra 0.333, the integer result will be correct.
398 // To handle the signed exponent one alternative would be to take absolute values, saving
399 // signs, etc - but that gets a bit complicated with 2-complement integers.
400 // Instead, we remember that we don't really want the exact division per se - what we're
401 // really after is only rewriting e = 3*n+m. That will actually be *easier* to handle if
402 // we require that m must be positive (fewer cases to handle) instead of having n as the
404 // To handle this we start by adding 127 to the exponent. This value corresponds to the
405 // exponent bias, minus 1 because frexp() has a different standard for the value it returns,
406 // but then we add 1 back to handle the extra 0.333 in 21845. So, we have offsetExp = e+127
407 // and then multiply by 21845 to get a division result offsetExpDiv3.
408 // A (signed) value for n is then recovered by subtracting 42 (bias-1)/3 from k.
409 // To calculate a strict remainder we should evaluate offsetExp - 3*offsetExpDiv3 - 1, where
410 // the extra 1 corrects for the value we added to the exponent to get correct division.
411 // This remainder would have the value 0,1, or 2, but since we only use it to select
412 // other numbers we can skip the last step and just handle the cases as 1,2 or 3 instead.
414 // OK; end of long detour. Here's how we actually do it in our implementation by using
415 // floating-point for the exponent instead to avoid needing integer shifts:
417 // 1) Convert the exponent (obtained from frexp) to a float
418 // 2) Calculate offsetExp = exp + offset. Note that we should not add the extra 1 here since we
419 // do floating-point division instead of our integer hack, so it's the exponent bias-1, or
420 // the largest exponent minus 2.
421 // 3) Divide the float by 3 by multiplying with 1/3
422 // 4) Truncate it to an integer to get the division result. This is potentially dangerous in
423 // combination with floating-point, because many integers cannot be represented exactly in
424 // floating point, and if we are just epsilon below the result might be truncated to a lower
425 // integer. I have not observed this on x86, but to have a safety margin we can add a small
426 // fraction - since we already know the fraction part should be either 0, 0.333..., or 0.666...
427 // We can even save this extra floating-point addition by adding a small fraction (0.1) when
428 // we introduce the exponent offset - that will correspond to a safety margin of 0.1/3, which is plenty.
429 // 5) Get the remainder part by subtracting the truncated floating-point part.
430 // Here too we will have a plain division, so the remainder is a strict modulus
431 // and will have the values 0, 1 or 2.
433 // Before worrying about the few wasted cycles due to longer fp latency, this has the
434 // additional advantage that we don't use a single integer operation, so the algorithm
435 // will work just A-OK on all SIMD implementations, which avoids diverging code paths.
437 // The 0.1 here is the safety margin due to truncation described in item 4 in the comments above.
438 SimdFloat offsetExp = cvtI2R(exponent) + SimdFloat(static_cast<float>(3 * offsetDiv3) + 0.1);
440 SimdFloat offsetExpDiv3 =
441 trunc(offsetExp * oneThird); // important to truncate here to mimic integer division
443 SimdFInt32 expDiv3 = cvtR2I(offsetExpDiv3 - SimdFloat(static_cast<float>(offsetDiv3)));
445 SimdFloat remainder = offsetExp - offsetExpDiv3 * three;
447 // If remainder is 0 we should just have the factor 1.0,
448 // so first pick 1.0 if it is below 0.5, and 2^(1/3) if it's above 0.5 (i.e., 1 or 2)
449 SimdFloat factor = blend(one, cbrt2, SimdFloat(0.5) < remainder);
450 // Second, we overwrite with 2^(2/3) if rem>1.5 (i.e., 2)
451 factor = blend(factor, sqrCbrt2, SimdFloat(1.5) < remainder);
453 // Assemble the non-signed fraction, and add the sign back by xor
454 SimdFloat fraction = (nom * invDenom * factor) ^ xSignBit;
455 // Load to IEEE754 number, and set result to 0.0 if x was 0.0 or denormal
456 SimdFloat result = selectByMask(ldexp(fraction, expDiv3), xIsNonZero);
461 /*! \brief Inverse cube root for SIMD floats
463 * \param x Argument to calculate cube root of. Can be positive or
464 * negative, but the magnitude cannot be lower than
465 * the smallest normal number.
466 * \return Cube root of x. Undefined for values that don't
467 * fulfill the restriction of abs(x) > minFloat.
469 static inline SimdFloat gmx_simdcall invcbrt(SimdFloat x)
471 const SimdFloat signBit(GMX_FLOAT_NEGZERO);
472 const SimdFloat minFloat(std::numeric_limits<float>::min());
473 // Bias is 128-1 = 127, which is not divisible by 3. Since the largest-magnitude
474 // negative exponent from frexp() is -126, we can subtract one more unit to get 126
475 // as offset, which is divisible by 3 (result 42). To avoid clang warnings about fragile integer
476 // division mixed with FP, we let the divided value (42) be the original constant.
477 const std::int32_t offsetDiv3(42);
478 const SimdFloat c2(-0.191502161678719066F);
479 const SimdFloat c1(0.697570460207922770F);
480 const SimdFloat c0(0.492659620528969547F);
481 const SimdFloat one(1.0F);
482 const SimdFloat two(2.0F);
483 const SimdFloat three(3.0F);
484 const SimdFloat oneThird(1.0F / 3.0F);
485 const SimdFloat invCbrt2(1.0F / 1.2599210498948731648F);
486 const SimdFloat invSqrCbrt2(1.0F / 1.5874010519681994748F);
488 // We use pretty much exactly the same implementation as for cbrt(x),
489 // but to compute the inverse we swap the nominator/denominator
490 // in the quotient, and also swap the sign of the exponent parts.
492 SimdFloat xSignBit = x & signBit; // create bit mask where the sign bit is 1 for x elements < 0
493 SimdFloat xAbs = andNot(signBit, x); // select everthing but the sign bit => abs(x)
496 SimdFloat y = frexp(xAbs, &exponent);
497 // For the mantissa (y) we will use a limited-range approximation of cbrt(y),
498 // by first using a polynomial and then evaluating
499 // Transform y to z = c2*y^2 + c1*y + c0, then w = z^3, and finally
500 // evaluate the quotient q = z * (w + 2 * y) / (2 * w + y).
501 SimdFloat z = fma(fma(y, c2, c1), y, c0);
502 SimdFloat w = z * z * z;
503 SimdFloat nom = fma(two, w, y);
504 SimdFloat invDenom = inv(z * fma(two, y, w));
506 // The 0.1 here is the safety margin due to truncation described in item 4 in the comments above.
507 SimdFloat offsetExp = cvtI2R(exponent) + SimdFloat(static_cast<float>(3 * offsetDiv3) + 0.1);
508 SimdFloat offsetExpDiv3 =
509 trunc(offsetExp * oneThird); // important to truncate here to mimic integer division
511 // We should swap the sign here, so we change order of the terms in the subtraction
512 SimdFInt32 expDiv3 = cvtR2I(SimdFloat(static_cast<float>(offsetDiv3)) - offsetExpDiv3);
514 // Swap sign here too, so remainder is either 0, -1 or -2
515 SimdFloat remainder = offsetExpDiv3 * three - offsetExp;
517 // If remainder is 0 we should just have the factor 1.0,
518 // so first pick 1.0 if it is above -0.5, and 2^(-1/3) if it's below -0.5 (i.e., -1 or -2)
519 SimdFloat factor = blend(one, invCbrt2, remainder < SimdFloat(-0.5));
520 // Second, we overwrite with 2^(-2/3) if rem<-1.5 (i.e., -2)
521 factor = blend(factor, invSqrCbrt2, remainder < SimdFloat(-1.5));
523 // Assemble the non-signed fraction, and add the sign back by xor
524 SimdFloat fraction = (nom * invDenom * factor) ^ xSignBit;
525 // Load to IEEE754 number, and set result to 0.0 if x was 0.0 or denormal
526 SimdFloat result = ldexp(fraction, expDiv3);
531 /*! \brief SIMD float log2(x). This is the base-2 logarithm.
533 * \param x Argument, should be >0.
534 * \result The base-2 logarithm of x. Undefined if argument is invalid.
536 static inline SimdFloat gmx_simdcall log2(SimdFloat x)
538 // This implementation computes log2 by
539 // 1) Extracting the exponent and adding it to...
540 // 2) A 9th-order minimax approximation using only odd
541 // terms of (x-1)/(x+1), where x is the mantissa.
543 # if GMX_SIMD_HAVE_NATIVE_LOG_FLOAT
544 // Just rescale if native log2() is not present, but log() is.
545 return log(x) * SimdFloat(std::log2(std::exp(1.0)));
547 const SimdFloat one(1.0F);
548 const SimdFloat two(2.0F);
549 const SimdFloat invsqrt2(1.0F / std::sqrt(2.0F));
550 const SimdFloat CL9(0.342149508897807708152F);
551 const SimdFloat CL7(0.411570606888219447939F);
552 const SimdFloat CL5(0.577085979152320294183F);
553 const SimdFloat CL3(0.961796550607099898222F);
554 const SimdFloat CL1(2.885390081777926774009F);
555 SimdFloat fExp, x2, p;
559 // For the log2() function, the argument can never be 0, so use the faster version of frexp()
560 x = frexp<MathOptimization::Unsafe>(x, &iExp);
564 // Adjust to non-IEEE format for x<1/sqrt(2): exponent -= 1, mantissa *= 2.0
565 fExp = fExp - selectByMask(one, m);
566 x = x * blend(one, two, m);
568 x = (x - one) * inv(x + one);
571 p = fma(CL9, x2, CL7);
581 # if !GMX_SIMD_HAVE_NATIVE_LOG_FLOAT
582 /*! \brief SIMD float log(x). This is the natural logarithm.
584 * \param x Argument, should be >0.
585 * \result The natural logarithm of x. Undefined if argument is invalid.
587 static inline SimdFloat gmx_simdcall log(SimdFloat x)
589 const SimdFloat one(1.0F);
590 const SimdFloat two(2.0F);
591 const SimdFloat invsqrt2(1.0F / std::sqrt(2.0F));
592 const SimdFloat corr(0.693147180559945286226764F);
593 const SimdFloat CL9(0.2371599674224853515625F);
594 const SimdFloat CL7(0.285279005765914916992188F);
595 const SimdFloat CL5(0.400005519390106201171875F);
596 const SimdFloat CL3(0.666666567325592041015625F);
597 const SimdFloat CL1(2.0F);
598 SimdFloat fExp, x2, p;
602 // For log(), the argument cannot be 0, so use the faster version of frexp()
603 x = frexp<MathOptimization::Unsafe>(x, &iExp);
607 // Adjust to non-IEEE format for x<1/sqrt(2): exponent -= 1, mantissa *= 2.0
608 fExp = fExp - selectByMask(one, m);
609 x = x * blend(one, two, m);
611 x = (x - one) * inv(x + one);
614 p = fma(CL9, x2, CL7);
618 p = fma(p, x, corr * fExp);
624 # if !GMX_SIMD_HAVE_NATIVE_EXP2_FLOAT
625 /*! \brief SIMD float 2^x
627 * \tparam opt If this is changed from the default (safe) into the unsafe
628 * option, input values that would otherwise lead to zero-clamped
629 * results are not allowed and will lead to undefined results.
631 * \param x Argument. For the default (safe) function version this can be
632 * arbitrarily small value, but the routine might clamp the result to
633 * zero for arguments that would produce subnormal IEEE754-2008 results.
634 * This corresponds to inputs below -126 in single or -1022 in double,
635 * and it might overflow for arguments reaching 127 (single) or
636 * 1023 (double). If you enable the unsafe math optimization,
637 * very small arguments will not necessarily be zero-clamped, but
638 * can produce undefined results.
640 * \result 2^x. The result is undefined for very large arguments that cause
641 * internal floating-point overflow. If unsafe optimizations are enabled,
642 * this is also true for very small values.
644 * \note The definition range of this function is just-so-slightly smaller
645 * than the allowed IEEE exponents for many architectures. This is due
646 * to the implementation, which will hopefully improve in the future.
648 * \warning You cannot rely on this implementation returning inf for arguments
649 * that cause overflow. If you have some very large
650 * values and need to rely on getting a valid numerical output,
651 * take the minimum of your variable and the largest valid argument
652 * before calling this routine.
654 template<MathOptimization opt = MathOptimization::Safe>
655 static inline SimdFloat gmx_simdcall exp2(SimdFloat x)
657 const SimdFloat CC6(0.0001534581200287996416911311F);
658 const SimdFloat CC5(0.001339993121934088894618990F);
659 const SimdFloat CC4(0.009618488957115180159497841F);
660 const SimdFloat CC3(0.05550328776964726865751735F);
661 const SimdFloat CC2(0.2402264689063408646490722F);
662 const SimdFloat CC1(0.6931472057372680777553816F);
663 const SimdFloat one(1.0F);
669 // Large negative values are valid arguments to exp2(), so there are two
670 // things we need to account for:
671 // 1. When the exponents reaches -127, the (biased) exponent field will be
672 // zero and we can no longer multiply with it. There are special IEEE
673 // formats to handle this range, but for now we have to accept that
674 // we cannot handle those arguments. If input value becomes even more
675 // negative, it will start to loop and we would end up with invalid
676 // exponents. Thus, we need to limit or mask this.
677 // 2. For VERY large negative values, we will have problems that the
678 // subtraction to get the fractional part loses accuracy, and then we
679 // can end up with overflows in the polynomial.
681 // For now, we handle this by forwarding the math optimization setting to
682 // ldexp, where the routine will return zero for very small arguments.
684 // However, before doing that we need to make sure we do not call cvtR2I
685 // with an argument that is so negative it cannot be converted to an integer.
686 if (opt == MathOptimization::Safe)
688 x = max(x, SimdFloat(std::numeric_limits<std::int32_t>::lowest()));
691 fexppart = ldexp<opt>(one, cvtR2I(x));
695 p = fma(CC6, x, CC5);
706 # if !GMX_SIMD_HAVE_NATIVE_EXP_FLOAT
707 /*! \brief SIMD float exp(x).
709 * In addition to scaling the argument for 2^x this routine correctly does
710 * extended precision arithmetics to improve accuracy.
712 * \tparam opt If this is changed from the default (safe) into the unsafe
713 * option, input values that would otherwise lead to zero-clamped
714 * results are not allowed and will lead to undefined results.
716 * \param x Argument. For the default (safe) function version this can be
717 * arbitrarily small value, but the routine might clamp the result to
718 * zero for arguments that would produce subnormal IEEE754-2008 results.
719 * This corresponds to input arguments reaching
720 * -126*ln(2)=-87.3 in single, or -1022*ln(2)=-708.4 (double).
721 * Similarly, it might overflow for arguments reaching
722 * 127*ln(2)=88.0 (single) or 1023*ln(2)=709.1 (double). If the
723 * unsafe math optimizations are enabled, small input values that would
724 * result in zero-clamped output are not allowed.
726 * \result exp(x). Overflowing arguments are likely to either return 0 or inf,
727 * depending on the underlying implementation. If unsafe optimizations
728 * are enabled, this is also true for very small values.
730 * \note The definition range of this function is just-so-slightly smaller
731 * than the allowed IEEE exponents for many architectures. This is due
732 * to the implementation, which will hopefully improve in the future.
734 * \warning You cannot rely on this implementation returning inf for arguments
735 * that cause overflow. If you have some very large
736 * values and need to rely on getting a valid numerical output,
737 * take the minimum of your variable and the largest valid argument
738 * before calling this routine.
740 template<MathOptimization opt = MathOptimization::Safe>
741 static inline SimdFloat gmx_simdcall exp(SimdFloat x)
743 const SimdFloat argscale(1.44269504088896341F);
744 const SimdFloat invargscale0(-0.693145751953125F);
745 const SimdFloat invargscale1(-1.428606765330187045e-06F);
746 const SimdFloat CC4(0.00136324646882712841033936F);
747 const SimdFloat CC3(0.00836596917361021041870117F);
748 const SimdFloat CC2(0.0416710823774337768554688F);
749 const SimdFloat CC1(0.166665524244308471679688F);
750 const SimdFloat CC0(0.499999850988388061523438F);
751 const SimdFloat one(1.0F);
756 // Large negative values are valid arguments to exp2(), so there are two
757 // things we need to account for:
758 // 1. When the exponents reaches -127, the (biased) exponent field will be
759 // zero and we can no longer multiply with it. There are special IEEE
760 // formats to handle this range, but for now we have to accept that
761 // we cannot handle those arguments. If input value becomes even more
762 // negative, it will start to loop and we would end up with invalid
763 // exponents. Thus, we need to limit or mask this.
764 // 2. For VERY large negative values, we will have problems that the
765 // subtraction to get the fractional part loses accuracy, and then we
766 // can end up with overflows in the polynomial.
768 // For now, we handle this by forwarding the math optimization setting to
769 // ldexp, where the routine will return zero for very small arguments.
771 // However, before doing that we need to make sure we do not call cvtR2I
772 // with an argument that is so negative it cannot be converted to an integer
773 // after being multiplied by argscale.
775 if (opt == MathOptimization::Safe)
777 x = max(x, SimdFloat(std::numeric_limits<std::int32_t>::lowest()) / argscale);
783 fexppart = ldexp<opt>(one, cvtR2I(y));
786 // Extended precision arithmetics
787 x = fma(invargscale0, intpart, x);
788 x = fma(invargscale1, intpart, x);
790 p = fma(CC4, x, CC3);
794 p = fma(x * x, p, x);
795 # if GMX_SIMD_HAVE_FMA
796 x = fma(p, fexppart, fexppart);
798 x = (p + one) * fexppart;
804 /*! \brief SIMD float pow(x,y)
806 * This returns x^y for SIMD values.
808 * \tparam opt If this is changed from the default (safe) into the unsafe
809 * option, there are no guarantees about correct results for x==0.
815 * \result x^y. Overflowing arguments are likely to either return 0 or inf,
816 * depending on the underlying implementation. If unsafe optimizations
817 * are enabled, this is also true for x==0.
819 * \warning You cannot rely on this implementation returning inf for arguments
820 * that cause overflow. If you have some very large
821 * values and need to rely on getting a valid numerical output,
822 * take the minimum of your variable and the largest valid argument
823 * before calling this routine.
825 template<MathOptimization opt = MathOptimization::Safe>
826 static inline SimdFloat gmx_simdcall pow(SimdFloat x, SimdFloat y)
830 if (opt == MathOptimization::Safe)
832 xcorr = max(x, SimdFloat(std::numeric_limits<float>::min()));
839 SimdFloat result = exp2<opt>(y * log2(xcorr));
841 if (opt == MathOptimization::Safe)
843 // if x==0 and y>0 we explicitly set the result to 0.0
844 // For any x with y==0, the result will already be 1.0 since we multiply by y (0.0) and call exp().
845 result = blend(result, setZero(), x == setZero() && setZero() < y);
852 /*! \brief SIMD float erf(x).
854 * \param x The value to calculate erf(x) for.
857 * This routine achieves very close to full precision, but we do not care about
858 * the last bit or the subnormal result range.
860 static inline SimdFloat gmx_simdcall erf(SimdFloat x)
862 // Coefficients for minimax approximation of erf(x)=x*P(x^2) in range [-1,1]
863 const SimdFloat CA6(7.853861353153693e-5F);
864 const SimdFloat CA5(-8.010193625184903e-4F);
865 const SimdFloat CA4(5.188327685732524e-3F);
866 const SimdFloat CA3(-2.685381193529856e-2F);
867 const SimdFloat CA2(1.128358514861418e-1F);
868 const SimdFloat CA1(-3.761262582423300e-1F);
869 const SimdFloat CA0(1.128379165726710F);
870 // Coefficients for minimax approximation of erfc(x)=Exp(-x^2)*P((1/(x-1))^2) in range [0.67,2]
871 const SimdFloat CB9(-0.0018629930017603923F);
872 const SimdFloat CB8(0.003909821287598495F);
873 const SimdFloat CB7(-0.0052094582210355615F);
874 const SimdFloat CB6(0.005685614362160572F);
875 const SimdFloat CB5(-0.0025367682853477272F);
876 const SimdFloat CB4(-0.010199799682318782F);
877 const SimdFloat CB3(0.04369575504816542F);
878 const SimdFloat CB2(-0.11884063474674492F);
879 const SimdFloat CB1(0.2732120154030589F);
880 const SimdFloat CB0(0.42758357702025784F);
881 // Coefficients for minimax approximation of erfc(x)=Exp(-x^2)*(1/x)*P((1/x)^2) in range [2,9.19]
882 const SimdFloat CC10(-0.0445555913112064F);
883 const SimdFloat CC9(0.21376355144663348F);
884 const SimdFloat CC8(-0.3473187200259257F);
885 const SimdFloat CC7(0.016690861551248114F);
886 const SimdFloat CC6(0.7560973182491192F);
887 const SimdFloat CC5(-1.2137903600145787F);
888 const SimdFloat CC4(0.8411872321232948F);
889 const SimdFloat CC3(-0.08670413896296343F);
890 const SimdFloat CC2(-0.27124782687240334F);
891 const SimdFloat CC1(-0.0007502488047806069F);
892 const SimdFloat CC0(0.5642114853803148F);
893 const SimdFloat one(1.0F);
894 const SimdFloat two(2.0F);
897 SimdFloat t, t2, w, w2;
898 SimdFloat pA0, pA1, pB0, pB1, pC0, pC1;
900 SimdFloat res_erf, res_erfc, res;
901 SimdFBool m, maskErf;
907 pA0 = fma(CA6, x4, CA4);
908 pA1 = fma(CA5, x4, CA3);
909 pA0 = fma(pA0, x4, CA2);
910 pA1 = fma(pA1, x4, CA1);
912 pA0 = fma(pA1, x2, pA0);
913 // Constant term must come last for precision reasons
920 maskErf = SimdFloat(0.75F) <= y;
921 t = maskzInv(y, maskErf);
926 // No need for a floating-point sieve here (as in erfc), since erf()
927 // will never return values that are extremely small for large args.
928 expmx2 = exp(-y * y);
930 pB1 = fma(CB9, w2, CB7);
931 pB0 = fma(CB8, w2, CB6);
932 pB1 = fma(pB1, w2, CB5);
933 pB0 = fma(pB0, w2, CB4);
934 pB1 = fma(pB1, w2, CB3);
935 pB0 = fma(pB0, w2, CB2);
936 pB1 = fma(pB1, w2, CB1);
937 pB0 = fma(pB0, w2, CB0);
938 pB0 = fma(pB1, w, pB0);
940 pC0 = fma(CC10, t2, CC8);
941 pC1 = fma(CC9, t2, CC7);
942 pC0 = fma(pC0, t2, CC6);
943 pC1 = fma(pC1, t2, CC5);
944 pC0 = fma(pC0, t2, CC4);
945 pC1 = fma(pC1, t2, CC3);
946 pC0 = fma(pC0, t2, CC2);
947 pC1 = fma(pC1, t2, CC1);
949 pC0 = fma(pC0, t2, CC0);
950 pC0 = fma(pC1, t, pC0);
953 // Select pB0 or pC0 for erfc()
955 res_erfc = blend(pB0, pC0, m);
956 res_erfc = res_erfc * expmx2;
958 // erfc(x<0) = 2-erfc(|x|)
960 res_erfc = blend(res_erfc, two - res_erfc, m);
962 // Select erf() or erfc()
963 res = blend(res_erf, one - res_erfc, maskErf);
968 /*! \brief SIMD float erfc(x).
970 * \param x The value to calculate erfc(x) for.
973 * This routine achieves full precision (bar the last bit) over most of the
974 * input range, but for large arguments where the result is getting close
975 * to the minimum representable numbers we accept slightly larger errors
976 * (think results that are in the ballpark of 10^-30 for single precision)
977 * since that is not relevant for MD.
979 static inline SimdFloat gmx_simdcall erfc(SimdFloat x)
981 // Coefficients for minimax approximation of erf(x)=x*P(x^2) in range [-1,1]
982 const SimdFloat CA6(7.853861353153693e-5F);
983 const SimdFloat CA5(-8.010193625184903e-4F);
984 const SimdFloat CA4(5.188327685732524e-3F);
985 const SimdFloat CA3(-2.685381193529856e-2F);
986 const SimdFloat CA2(1.128358514861418e-1F);
987 const SimdFloat CA1(-3.761262582423300e-1F);
988 const SimdFloat CA0(1.128379165726710F);
989 // Coefficients for minimax approximation of erfc(x)=Exp(-x^2)*P((1/(x-1))^2) in range [0.67,2]
990 const SimdFloat CB9(-0.0018629930017603923F);
991 const SimdFloat CB8(0.003909821287598495F);
992 const SimdFloat CB7(-0.0052094582210355615F);
993 const SimdFloat CB6(0.005685614362160572F);
994 const SimdFloat CB5(-0.0025367682853477272F);
995 const SimdFloat CB4(-0.010199799682318782F);
996 const SimdFloat CB3(0.04369575504816542F);
997 const SimdFloat CB2(-0.11884063474674492F);
998 const SimdFloat CB1(0.2732120154030589F);
999 const SimdFloat CB0(0.42758357702025784F);
1000 // Coefficients for minimax approximation of erfc(x)=Exp(-x^2)*(1/x)*P((1/x)^2) in range [2,9.19]
1001 const SimdFloat CC10(-0.0445555913112064F);
1002 const SimdFloat CC9(0.21376355144663348F);
1003 const SimdFloat CC8(-0.3473187200259257F);
1004 const SimdFloat CC7(0.016690861551248114F);
1005 const SimdFloat CC6(0.7560973182491192F);
1006 const SimdFloat CC5(-1.2137903600145787F);
1007 const SimdFloat CC4(0.8411872321232948F);
1008 const SimdFloat CC3(-0.08670413896296343F);
1009 const SimdFloat CC2(-0.27124782687240334F);
1010 const SimdFloat CC1(-0.0007502488047806069F);
1011 const SimdFloat CC0(0.5642114853803148F);
1012 // Coefficients for expansion of exp(x) in [0,0.1]
1013 // CD0 and CD1 are both 1.0, so no need to declare them separately
1014 const SimdFloat CD2(0.5000066608081202F);
1015 const SimdFloat CD3(0.1664795422874624F);
1016 const SimdFloat CD4(0.04379839977652482F);
1017 const SimdFloat one(1.0F);
1018 const SimdFloat two(2.0F);
1020 /* We need to use a small trick here, since we cannot assume all SIMD
1021 * architectures support integers, and the flag we want (0xfffff000) would
1022 * evaluate to NaN (i.e., it cannot be expressed as a floating-point num).
1023 * Instead, we represent the flags 0xf0f0f000 and 0x0f0f0000 as valid
1024 * fp numbers, and perform a logical or. Since the expression is constant,
1025 * we can at least hope it is evaluated at compile-time.
1027 # if GMX_SIMD_HAVE_LOGICAL
1028 const SimdFloat sieve(SimdFloat(-5.965323564e+29F) | SimdFloat(7.05044434e-30F));
1030 const int isieve = 0xFFFFF000;
1031 alignas(GMX_SIMD_ALIGNMENT) float mem[GMX_SIMD_FLOAT_WIDTH];
1041 SimdFloat x2, x4, y;
1042 SimdFloat q, z, t, t2, w, w2;
1043 SimdFloat pA0, pA1, pB0, pB1, pC0, pC1;
1044 SimdFloat expmx2, corr;
1045 SimdFloat res_erf, res_erfc, res;
1046 SimdFBool m, msk_erf;
1052 pA0 = fma(CA6, x4, CA4);
1053 pA1 = fma(CA5, x4, CA3);
1054 pA0 = fma(pA0, x4, CA2);
1055 pA1 = fma(pA1, x4, CA1);
1057 pA0 = fma(pA0, x4, pA1);
1058 // Constant term must come last for precision reasons
1065 msk_erf = SimdFloat(0.75F) <= y;
1066 t = maskzInv(y, msk_erf);
1071 * We cannot simply calculate exp(-y2) directly in single precision, since
1072 * that will lose a couple of bits of precision due to the multiplication.
1073 * Instead, we introduce y=z+w, where the last 12 bits of precision are in w.
1074 * Then we get exp(-y2) = exp(-z2)*exp((z-y)*(z+y)).
1076 * The only drawback with this is that it requires TWO separate exponential
1077 * evaluations, which would be horrible performance-wise. However, the argument
1078 * for the second exp() call is always small, so there we simply use a
1079 * low-order minimax expansion on [0,0.1].
1081 * However, this neat idea requires support for logical ops (and) on
1082 * FP numbers, which some vendors decided isn't necessary in their SIMD
1083 * instruction sets (Hi, IBM VSX!). In principle we could use some tricks
1084 * in double, but we still need memory as a backup when that is not available,
1085 * and this case is rare enough that we go directly there...
1087 # if GMX_SIMD_HAVE_LOGICAL
1091 for (i = 0; i < GMX_SIMD_FLOAT_WIDTH; i++)
1094 conv.i = conv.i & isieve;
1097 z = load<SimdFloat>(mem);
1099 q = (z - y) * (z + y);
1100 corr = fma(CD4, q, CD3);
1101 corr = fma(corr, q, CD2);
1102 corr = fma(corr, q, one);
1103 corr = fma(corr, q, one);
1105 expmx2 = exp(-z * z);
1106 expmx2 = expmx2 * corr;
1108 pB1 = fma(CB9, w2, CB7);
1109 pB0 = fma(CB8, w2, CB6);
1110 pB1 = fma(pB1, w2, CB5);
1111 pB0 = fma(pB0, w2, CB4);
1112 pB1 = fma(pB1, w2, CB3);
1113 pB0 = fma(pB0, w2, CB2);
1114 pB1 = fma(pB1, w2, CB1);
1115 pB0 = fma(pB0, w2, CB0);
1116 pB0 = fma(pB1, w, pB0);
1118 pC0 = fma(CC10, t2, CC8);
1119 pC1 = fma(CC9, t2, CC7);
1120 pC0 = fma(pC0, t2, CC6);
1121 pC1 = fma(pC1, t2, CC5);
1122 pC0 = fma(pC0, t2, CC4);
1123 pC1 = fma(pC1, t2, CC3);
1124 pC0 = fma(pC0, t2, CC2);
1125 pC1 = fma(pC1, t2, CC1);
1127 pC0 = fma(pC0, t2, CC0);
1128 pC0 = fma(pC1, t, pC0);
1131 // Select pB0 or pC0 for erfc()
1133 res_erfc = blend(pB0, pC0, m);
1134 res_erfc = res_erfc * expmx2;
1136 // erfc(x<0) = 2-erfc(|x|)
1138 res_erfc = blend(res_erfc, two - res_erfc, m);
1140 // Select erf() or erfc()
1141 res = blend(one - res_erf, res_erfc, msk_erf);
1146 /*! \brief SIMD float sin \& cos.
1148 * \param x The argument to evaluate sin/cos for
1149 * \param[out] sinval Sin(x)
1150 * \param[out] cosval Cos(x)
1152 * This version achieves close to machine precision, but for very large
1153 * magnitudes of the argument we inherently begin to lose accuracy due to the
1154 * argument reduction, despite using extended precision arithmetics internally.
1156 static inline void gmx_simdcall sincos(SimdFloat x, SimdFloat* sinval, SimdFloat* cosval)
1158 // Constants to subtract Pi/4*x from y while minimizing precision loss
1159 const SimdFloat argred0(-1.5703125);
1160 const SimdFloat argred1(-4.83751296997070312500e-04F);
1161 const SimdFloat argred2(-7.54953362047672271729e-08F);
1162 const SimdFloat argred3(-2.56334406825708960298e-12F);
1163 const SimdFloat two_over_pi(static_cast<float>(2.0F / M_PI));
1164 const SimdFloat const_sin2(-1.9515295891e-4F);
1165 const SimdFloat const_sin1(8.3321608736e-3F);
1166 const SimdFloat const_sin0(-1.6666654611e-1F);
1167 const SimdFloat const_cos2(2.443315711809948e-5F);
1168 const SimdFloat const_cos1(-1.388731625493765e-3F);
1169 const SimdFloat const_cos0(4.166664568298827e-2F);
1170 const SimdFloat half(0.5F);
1171 const SimdFloat one(1.0F);
1172 SimdFloat ssign, csign;
1173 SimdFloat x2, y, z, psin, pcos, sss, ccc;
1176 # if GMX_SIMD_HAVE_FINT32_ARITHMETICS && GMX_SIMD_HAVE_LOGICAL
1177 const SimdFInt32 ione(1);
1178 const SimdFInt32 itwo(2);
1181 z = x * two_over_pi;
1185 m = cvtIB2B((iy & ione) == SimdFInt32(0));
1186 ssign = selectByMask(SimdFloat(GMX_FLOAT_NEGZERO), cvtIB2B((iy & itwo) == itwo));
1187 csign = selectByMask(SimdFloat(GMX_FLOAT_NEGZERO), cvtIB2B(((iy + ione) & itwo) == itwo));
1189 const SimdFloat quarter(0.25f);
1190 const SimdFloat minusquarter(-0.25f);
1192 SimdFBool m1, m2, m3;
1194 /* The most obvious way to find the arguments quadrant in the unit circle
1195 * to calculate the sign is to use integer arithmetic, but that is not
1196 * present in all SIMD implementations. As an alternative, we have devised a
1197 * pure floating-point algorithm that uses truncation for argument reduction
1198 * so that we get a new value 0<=q<1 over the unit circle, and then
1199 * do floating-point comparisons with fractions. This is likely to be
1200 * slightly slower (~10%) due to the longer latencies of floating-point, so
1201 * we only use it when integer SIMD arithmetic is not present.
1205 // It is critical that half-way cases are rounded down
1206 z = fma(x, two_over_pi, half);
1210 /* z now starts at 0.0 for x=-pi/4 (although neg. values cannot occur), and
1211 * then increased by 1.0 as x increases by 2*Pi, when it resets to 0.0.
1212 * This removes the 2*Pi periodicity without using any integer arithmetic.
1213 * First check if y had the value 2 or 3, set csign if true.
1216 /* If we have logical operations we can work directly on the signbit, which
1217 * saves instructions. Otherwise we need to represent signs as +1.0/-1.0.
1218 * Thus, if you are altering defines to debug alternative code paths, the
1219 * two GMX_SIMD_HAVE_LOGICAL sections in this routine must either both be
1220 * active or inactive - you will get errors if only one is used.
1222 # if GMX_SIMD_HAVE_LOGICAL
1223 ssign = ssign & SimdFloat(GMX_FLOAT_NEGZERO);
1224 csign = andNot(q, SimdFloat(GMX_FLOAT_NEGZERO));
1225 ssign = ssign ^ csign;
1227 ssign = copysign(SimdFloat(1.0f), ssign);
1228 csign = copysign(SimdFloat(1.0f), q);
1230 ssign = ssign * csign; // swap ssign if csign was set.
1232 // Check if y had value 1 or 3 (remember we subtracted 0.5 from q)
1233 m1 = (q < minusquarter);
1234 m2 = (setZero() <= q);
1238 // where mask is FALSE, swap sign.
1239 csign = csign * blend(SimdFloat(-1.0f), one, m);
1241 x = fma(y, argred0, x);
1242 x = fma(y, argred1, x);
1243 x = fma(y, argred2, x);
1244 x = fma(y, argred3, x);
1247 psin = fma(const_sin2, x2, const_sin1);
1248 psin = fma(psin, x2, const_sin0);
1249 psin = fma(psin, x * x2, x);
1250 pcos = fma(const_cos2, x2, const_cos1);
1251 pcos = fma(pcos, x2, const_cos0);
1252 pcos = fms(pcos, x2, half);
1253 pcos = fma(pcos, x2, one);
1255 sss = blend(pcos, psin, m);
1256 ccc = blend(psin, pcos, m);
1257 // See comment for GMX_SIMD_HAVE_LOGICAL section above.
1258 # if GMX_SIMD_HAVE_LOGICAL
1259 *sinval = sss ^ ssign;
1260 *cosval = ccc ^ csign;
1262 *sinval = sss * ssign;
1263 *cosval = ccc * csign;
1267 /*! \brief SIMD float sin(x).
1269 * \param x The argument to evaluate sin for
1272 * \attention Do NOT call both sin & cos if you need both results, since each of them
1273 * will then call \ref sincos and waste a factor 2 in performance.
1275 static inline SimdFloat gmx_simdcall sin(SimdFloat x)
1282 /*! \brief SIMD float cos(x).
1284 * \param x The argument to evaluate cos for
1287 * \attention Do NOT call both sin & cos if you need both results, since each of them
1288 * will then call \ref sincos and waste a factor 2 in performance.
1290 static inline SimdFloat gmx_simdcall cos(SimdFloat x)
1297 /*! \brief SIMD float tan(x).
1299 * \param x The argument to evaluate tan for
1302 static inline SimdFloat gmx_simdcall tan(SimdFloat x)
1304 const SimdFloat argred0(-1.5703125);
1305 const SimdFloat argred1(-4.83751296997070312500e-04F);
1306 const SimdFloat argred2(-7.54953362047672271729e-08F);
1307 const SimdFloat argred3(-2.56334406825708960298e-12F);
1308 const SimdFloat two_over_pi(static_cast<float>(2.0F / M_PI));
1309 const SimdFloat CT6(0.009498288995810566122993911);
1310 const SimdFloat CT5(0.002895755790837379295226923);
1311 const SimdFloat CT4(0.02460087336161924491836265);
1312 const SimdFloat CT3(0.05334912882656359828045988);
1313 const SimdFloat CT2(0.1333989091464957704418495);
1314 const SimdFloat CT1(0.3333307599244198227797507);
1316 SimdFloat x2, p, y, z;
1319 # if GMX_SIMD_HAVE_FINT32_ARITHMETICS && GMX_SIMD_HAVE_LOGICAL
1323 z = x * two_over_pi;
1326 m = cvtIB2B((iy & ione) == ione);
1328 x = fma(y, argred0, x);
1329 x = fma(y, argred1, x);
1330 x = fma(y, argred2, x);
1331 x = fma(y, argred3, x);
1332 x = selectByMask(SimdFloat(GMX_FLOAT_NEGZERO), m) ^ x;
1334 const SimdFloat quarter(0.25f);
1335 const SimdFloat half(0.5f);
1336 const SimdFloat threequarter(0.75f);
1338 SimdFBool m1, m2, m3;
1341 z = fma(w, two_over_pi, half);
1347 m3 = threequarter <= q;
1350 w = fma(y, argred0, w);
1351 w = fma(y, argred1, w);
1352 w = fma(y, argred2, w);
1353 w = fma(y, argred3, w);
1354 w = blend(w, -w, m);
1355 x = w * copysign(SimdFloat(1.0), x);
1358 p = fma(CT6, x2, CT5);
1359 p = fma(p, x2, CT4);
1360 p = fma(p, x2, CT3);
1361 p = fma(p, x2, CT2);
1362 p = fma(p, x2, CT1);
1363 p = fma(x2 * p, x, x);
1365 p = blend(p, maskzInv(p, m), m);
1369 /*! \brief SIMD float asin(x).
1371 * \param x The argument to evaluate asin for
1374 static inline SimdFloat gmx_simdcall asin(SimdFloat x)
1376 const SimdFloat limitlow(1e-4F);
1377 const SimdFloat half(0.5F);
1378 const SimdFloat one(1.0F);
1379 const SimdFloat halfpi(static_cast<float>(M_PI / 2.0F));
1380 const SimdFloat CC5(4.2163199048E-2F);
1381 const SimdFloat CC4(2.4181311049E-2F);
1382 const SimdFloat CC3(4.5470025998E-2F);
1383 const SimdFloat CC2(7.4953002686E-2F);
1384 const SimdFloat CC1(1.6666752422E-1F);
1386 SimdFloat z, z1, z2, q, q1, q2;
1392 z1 = half * (one - xabs);
1394 q1 = z1 * maskzInvsqrt(z1, m2);
1397 z = blend(z2, z1, m);
1398 q = blend(q2, q1, m);
1401 pA = fma(CC5, z2, CC3);
1402 pB = fma(CC4, z2, CC2);
1403 pA = fma(pA, z2, CC1);
1405 z = fma(pB, z2, pA);
1409 z = blend(z, q2, m);
1411 m = limitlow < xabs;
1412 z = blend(xabs, z, m);
1418 /*! \brief SIMD float acos(x).
1420 * \param x The argument to evaluate acos for
1423 static inline SimdFloat gmx_simdcall acos(SimdFloat x)
1425 const SimdFloat one(1.0F);
1426 const SimdFloat half(0.5F);
1427 const SimdFloat pi(static_cast<float>(M_PI));
1428 const SimdFloat halfpi(static_cast<float>(M_PI / 2.0F));
1430 SimdFloat z, z1, z2, z3;
1431 SimdFBool m1, m2, m3;
1437 z = fnma(half, xabs, half);
1439 z = z * maskzInvsqrt(z, m3);
1440 z = blend(x, z, m1);
1446 z = blend(z1, z2, m2);
1447 z = blend(z3, z, m1);
1452 /*! \brief SIMD float asin(x).
1454 * \param x The argument to evaluate atan for
1455 * \result Atan(x), same argument/value range as standard math library.
1457 static inline SimdFloat gmx_simdcall atan(SimdFloat x)
1459 const SimdFloat halfpi(static_cast<float>(M_PI / 2.0F));
1460 const SimdFloat CA17(0.002823638962581753730774F);
1461 const SimdFloat CA15(-0.01595690287649631500244F);
1462 const SimdFloat CA13(0.04250498861074447631836F);
1463 const SimdFloat CA11(-0.07489009201526641845703F);
1464 const SimdFloat CA9(0.1063479334115982055664F);
1465 const SimdFloat CA7(-0.1420273631811141967773F);
1466 const SimdFloat CA5(0.1999269574880599975585F);
1467 const SimdFloat CA3(-0.3333310186862945556640F);
1468 const SimdFloat one(1.0F);
1469 SimdFloat x2, x3, x4, pA, pB;
1475 x = blend(x, maskzInv(x, m2), m2);
1480 pA = fma(CA17, x4, CA13);
1481 pB = fma(CA15, x4, CA11);
1482 pA = fma(pA, x4, CA9);
1483 pB = fma(pB, x4, CA7);
1484 pA = fma(pA, x4, CA5);
1485 pB = fma(pB, x4, CA3);
1486 pA = fma(pA, x2, pB);
1487 pA = fma(pA, x3, x);
1489 pA = blend(pA, halfpi - pA, m2);
1490 pA = blend(pA, -pA, m);
1495 /*! \brief SIMD float atan2(y,x).
1497 * \param y Y component of vector, any quartile
1498 * \param x X component of vector, any quartile
1499 * \result Atan(y,x), same argument/value range as standard math library.
1501 * \note This routine should provide correct results for all finite
1502 * non-zero or positive-zero arguments. However, negative zero arguments will
1503 * be treated as positive zero, which means the return value will deviate from
1504 * the standard math library atan2(y,x) for those cases. That should not be
1505 * of any concern in Gromacs, and in particular it will not affect calculations
1506 * of angles from vectors.
1508 static inline SimdFloat gmx_simdcall atan2(SimdFloat y, SimdFloat x)
1510 const SimdFloat pi(static_cast<float>(M_PI));
1511 const SimdFloat halfpi(static_cast<float>(M_PI / 2.0));
1512 SimdFloat xinv, p, aoffset;
1513 SimdFBool mask_xnz, mask_ynz, mask_xlt0, mask_ylt0;
1515 mask_xnz = x != setZero();
1516 mask_ynz = y != setZero();
1517 mask_xlt0 = x < setZero();
1518 mask_ylt0 = y < setZero();
1520 aoffset = selectByNotMask(halfpi, mask_xnz);
1521 aoffset = selectByMask(aoffset, mask_ynz);
1523 aoffset = blend(aoffset, pi, mask_xlt0);
1524 aoffset = blend(aoffset, -aoffset, mask_ylt0);
1526 xinv = maskzInv(x, mask_xnz);
1534 /*! \brief Calculate the force correction due to PME analytically in SIMD float.
1536 * \param z2 \f$(r \beta)^2\f$ - see below for details.
1537 * \result Correction factor to coulomb force - see below for details.
1539 * This routine is meant to enable analytical evaluation of the
1540 * direct-space PME electrostatic force to avoid tables.
1542 * The direct-space potential should be \f$ \mbox{erfc}(\beta r)/r\f$, but there
1543 * are some problems evaluating that:
1545 * First, the error function is difficult (read: expensive) to
1546 * approxmiate accurately for intermediate to large arguments, and
1547 * this happens already in ranges of \f$(\beta r)\f$ that occur in simulations.
1548 * Second, we now try to avoid calculating potentials in Gromacs but
1549 * use forces directly.
1551 * We can simply things slight by noting that the PME part is really
1552 * a correction to the normal Coulomb force since \f$\mbox{erfc}(z)=1-\mbox{erf}(z)\f$, i.e.
1554 * V = \frac{1}{r} - \frac{\mbox{erf}(\beta r)}{r}
1556 * The first term we already have from the inverse square root, so
1557 * that we can leave out of this routine.
1559 * For pme tolerances of 1e-3 to 1e-8 and cutoffs of 0.5nm to 1.8nm,
1560 * the argument \f$beta r\f$ will be in the range 0.15 to ~4, which is
1561 * the range used for the minimax fit. Use your favorite plotting program
1562 * to realize how well-behaved \f$\frac{\mbox{erf}(z)}{z}\f$ is in this range!
1564 * We approximate \f$f(z)=\mbox{erf}(z)/z\f$ with a rational minimax polynomial.
1565 * However, it turns out it is more efficient to approximate \f$f(z)/z\f$ and
1566 * then only use even powers. This is another minor optimization, since
1567 * we actually \a want \f$f(z)/z\f$, because it is going to be multiplied by
1568 * the vector between the two atoms to get the vectorial force. The
1569 * fastest flops are the ones we can avoid calculating!
1571 * So, here's how it should be used:
1573 * 1. Calculate \f$r^2\f$.
1574 * 2. Multiply by \f$\beta^2\f$, so you get \f$z^2=(\beta r)^2\f$.
1575 * 3. Evaluate this routine with \f$z^2\f$ as the argument.
1576 * 4. The return value is the expression:
1579 * \frac{2 \exp{-z^2}}{\sqrt{\pi} z^2}-\frac{\mbox{erf}(z)}{z^3}
1582 * 5. Multiply the entire expression by \f$\beta^3\f$. This will get you
1585 * \frac{2 \beta^3 \exp(-z^2)}{\sqrt{\pi} z^2} - \frac{\beta^3 \mbox{erf}(z)}{z^3}
1588 * or, switching back to \f$r\f$ (since \f$z=r \beta\f$):
1591 * \frac{2 \beta \exp(-r^2 \beta^2)}{\sqrt{\pi} r^2} - \frac{\mbox{erf}(r \beta)}{r^3}
1594 * With a bit of math exercise you should be able to confirm that
1598 * \frac{\frac{d}{dr}\left( \frac{\mbox{erf}(\beta r)}{r} \right)}{r}
1601 * 6. Add the result to \f$r^{-3}\f$, multiply by the product of the charges,
1602 * and you have your force (divided by \f$r\f$). A final multiplication
1603 * with the vector connecting the two particles and you have your
1604 * vectorial force to add to the particles.
1606 * This approximation achieves an error slightly lower than 1e-6
1607 * in single precision and 1e-11 in double precision
1608 * for arguments smaller than 16 (\f$\beta r \leq 4 \f$);
1609 * when added to \f$1/r\f$ the error will be insignificant.
1610 * For \f$\beta r \geq 7206\f$ the return value can be inf or NaN.
1613 static inline SimdFloat gmx_simdcall pmeForceCorrection(SimdFloat z2)
1615 const SimdFloat FN6(-1.7357322914161492954e-8F);
1616 const SimdFloat FN5(1.4703624142580877519e-6F);
1617 const SimdFloat FN4(-0.000053401640219807709149F);
1618 const SimdFloat FN3(0.0010054721316683106153F);
1619 const SimdFloat FN2(-0.019278317264888380590F);
1620 const SimdFloat FN1(0.069670166153766424023F);
1621 const SimdFloat FN0(-0.75225204789749321333F);
1623 const SimdFloat FD4(0.0011193462567257629232F);
1624 const SimdFloat FD3(0.014866955030185295499F);
1625 const SimdFloat FD2(0.11583842382862377919F);
1626 const SimdFloat FD1(0.50736591960530292870F);
1627 const SimdFloat FD0(1.0F);
1630 SimdFloat polyFN0, polyFN1, polyFD0, polyFD1;
1634 polyFD0 = fma(FD4, z4, FD2);
1635 polyFD1 = fma(FD3, z4, FD1);
1636 polyFD0 = fma(polyFD0, z4, FD0);
1637 polyFD0 = fma(polyFD1, z2, polyFD0);
1639 polyFD0 = inv(polyFD0);
1641 polyFN0 = fma(FN6, z4, FN4);
1642 polyFN1 = fma(FN5, z4, FN3);
1643 polyFN0 = fma(polyFN0, z4, FN2);
1644 polyFN1 = fma(polyFN1, z4, FN1);
1645 polyFN0 = fma(polyFN0, z4, FN0);
1646 polyFN0 = fma(polyFN1, z2, polyFN0);
1648 return polyFN0 * polyFD0;
1652 /*! \brief Calculate the potential correction due to PME analytically in SIMD float.
1654 * \param z2 \f$(r \beta)^2\f$ - see below for details.
1655 * \result Correction factor to coulomb potential - see below for details.
1657 * See \ref pmeForceCorrection for details about the approximation.
1659 * This routine calculates \f$\mbox{erf}(z)/z\f$, although you should provide \f$z^2\f$
1660 * as the input argument.
1662 * Here's how it should be used:
1664 * 1. Calculate \f$r^2\f$.
1665 * 2. Multiply by \f$\beta^2\f$, so you get \f$z^2=\beta^2*r^2\f$.
1666 * 3. Evaluate this routine with z^2 as the argument.
1667 * 4. The return value is the expression:
1670 * \frac{\mbox{erf}(z)}{z}
1673 * 5. Multiply the entire expression by beta and switching back to \f$r\f$ (since \f$z=r \beta\f$):
1676 * \frac{\mbox{erf}(r \beta)}{r}
1679 * 6. Subtract the result from \f$1/r\f$, multiply by the product of the charges,
1680 * and you have your potential.
1682 * This approximation achieves an error slightly lower than 1e-6
1683 * in single precision and 4e-11 in double precision
1684 * for arguments smaller than 16 (\f$ 0.15 \leq \beta r \leq 4 \f$);
1685 * for \f$ \beta r \leq 0.15\f$ the error can be twice as high;
1686 * when added to \f$1/r\f$ the error will be insignificant.
1687 * For \f$\beta r \geq 7142\f$ the return value can be inf or NaN.
1689 static inline SimdFloat gmx_simdcall pmePotentialCorrection(SimdFloat z2)
1691 const SimdFloat VN6(1.9296833005951166339e-8F);
1692 const SimdFloat VN5(-1.4213390571557850962e-6F);
1693 const SimdFloat VN4(0.000041603292906656984871F);
1694 const SimdFloat VN3(-0.00013134036773265025626F);
1695 const SimdFloat VN2(0.038657983986041781264F);
1696 const SimdFloat VN1(0.11285044772717598220F);
1697 const SimdFloat VN0(1.1283802385263030286F);
1699 const SimdFloat VD3(0.0066752224023576045451F);
1700 const SimdFloat VD2(0.078647795836373922256F);
1701 const SimdFloat VD1(0.43336185284710920150F);
1702 const SimdFloat VD0(1.0F);
1705 SimdFloat polyVN0, polyVN1, polyVD0, polyVD1;
1709 polyVD1 = fma(VD3, z4, VD1);
1710 polyVD0 = fma(VD2, z4, VD0);
1711 polyVD0 = fma(polyVD1, z2, polyVD0);
1713 polyVD0 = inv(polyVD0);
1715 polyVN0 = fma(VN6, z4, VN4);
1716 polyVN1 = fma(VN5, z4, VN3);
1717 polyVN0 = fma(polyVN0, z4, VN2);
1718 polyVN1 = fma(polyVN1, z4, VN1);
1719 polyVN0 = fma(polyVN0, z4, VN0);
1720 polyVN0 = fma(polyVN1, z2, polyVN0);
1722 return polyVN0 * polyVD0;
1728 # if GMX_SIMD_HAVE_DOUBLE
1731 /*! \name Double precision SIMD math functions
1733 * \note In most cases you should use the real-precision functions instead.
1737 /****************************************
1738 * DOUBLE PRECISION SIMD MATH FUNCTIONS *
1739 ****************************************/
1741 # if !GMX_SIMD_HAVE_NATIVE_COPYSIGN_DOUBLE
1742 /*! \brief Composes floating point value with the magnitude of x and the sign of y.
1744 * \param x Values to set sign for
1745 * \param y Values used to set sign
1746 * \return Magnitude of x, sign of y
1748 static inline SimdDouble gmx_simdcall copysign(SimdDouble x, SimdDouble y)
1750 # if GMX_SIMD_HAVE_LOGICAL
1751 return abs(x) | (SimdDouble(GMX_DOUBLE_NEGZERO) & y);
1753 return blend(abs(x), -abs(x), (y < setZero()));
1758 # if !GMX_SIMD_HAVE_NATIVE_RSQRT_ITER_DOUBLE
1759 /*! \brief Perform one Newton-Raphson iteration to improve 1/sqrt(x) for SIMD double.
1761 * This is a low-level routine that should only be used by SIMD math routine
1762 * that evaluates the inverse square root.
1764 * \param lu Approximation of 1/sqrt(x), typically obtained from lookup.
1765 * \param x The reference (starting) value x for which we want 1/sqrt(x).
1766 * \return An improved approximation with roughly twice as many bits of accuracy.
1768 static inline SimdDouble gmx_simdcall rsqrtIter(SimdDouble lu, SimdDouble x)
1770 SimdDouble tmp1 = x * lu;
1771 SimdDouble tmp2 = SimdDouble(-0.5) * lu;
1772 tmp1 = fma(tmp1, lu, SimdDouble(-3.0));
1777 /*! \brief Calculate 1/sqrt(x) for SIMD double.
1779 * \param x Argument that must be larger than GMX_FLOAT_MIN and smaller than
1780 * GMX_FLOAT_MAX, i.e. within the range of single precision.
1781 * For the single precision implementation this is obviously always
1782 * true for positive values, but for double precision it adds an
1783 * extra restriction since the first lookup step might have to be
1784 * performed in single precision on some architectures. Note that the
1785 * responsibility for checking falls on you - this routine does not
1788 * \return 1/sqrt(x). Result is undefined if your argument was invalid.
1790 static inline SimdDouble gmx_simdcall invsqrt(SimdDouble x)
1792 SimdDouble lu = rsqrt(x);
1793 # if (GMX_SIMD_RSQRT_BITS < GMX_SIMD_ACCURACY_BITS_DOUBLE)
1794 lu = rsqrtIter(lu, x);
1796 # if (GMX_SIMD_RSQRT_BITS * 2 < GMX_SIMD_ACCURACY_BITS_DOUBLE)
1797 lu = rsqrtIter(lu, x);
1799 # if (GMX_SIMD_RSQRT_BITS * 4 < GMX_SIMD_ACCURACY_BITS_DOUBLE)
1800 lu = rsqrtIter(lu, x);
1802 # if (GMX_SIMD_RSQRT_BITS * 8 < GMX_SIMD_ACCURACY_BITS_DOUBLE)
1803 lu = rsqrtIter(lu, x);
1808 /*! \brief Calculate 1/sqrt(x) for two SIMD doubles.
1810 * \param x0 First set of arguments, x0 must be in single range (see below).
1811 * \param x1 Second set of arguments, x1 must be in single range (see below).
1812 * \param[out] out0 Result 1/sqrt(x0)
1813 * \param[out] out1 Result 1/sqrt(x1)
1815 * In particular for double precision we can sometimes calculate square root
1816 * pairs slightly faster by using single precision until the very last step.
1818 * \note Both arguments must be larger than GMX_FLOAT_MIN and smaller than
1819 * GMX_FLOAT_MAX, i.e. within the range of single precision.
1820 * For the single precision implementation this is obviously always
1821 * true for positive values, but for double precision it adds an
1822 * extra restriction since the first lookup step might have to be
1823 * performed in single precision on some architectures. Note that the
1824 * responsibility for checking falls on you - this routine does not
1827 static inline void gmx_simdcall invsqrtPair(SimdDouble x0, SimdDouble x1, SimdDouble* out0, SimdDouble* out1)
1829 # if GMX_SIMD_HAVE_FLOAT && (GMX_SIMD_FLOAT_WIDTH == 2 * GMX_SIMD_DOUBLE_WIDTH) \
1830 && (GMX_SIMD_RSQRT_BITS < 22)
1831 SimdFloat xf = cvtDD2F(x0, x1);
1832 SimdFloat luf = rsqrt(xf);
1833 SimdDouble lu0, lu1;
1834 // Intermediate target is single - mantissa+1 bits
1835 # if (GMX_SIMD_RSQRT_BITS < GMX_SIMD_ACCURACY_BITS_SINGLE)
1836 luf = rsqrtIter(luf, xf);
1838 # if (GMX_SIMD_RSQRT_BITS * 2 < GMX_SIMD_ACCURACY_BITS_SINGLE)
1839 luf = rsqrtIter(luf, xf);
1841 # if (GMX_SIMD_RSQRT_BITS * 4 < GMX_SIMD_ACCURACY_BITS_SINGLE)
1842 luf = rsqrtIter(luf, xf);
1844 cvtF2DD(luf, &lu0, &lu1);
1845 // Last iteration(s) performed in double - if we had 22 bits, this gets us to 44 (~1e-15)
1846 # if (GMX_SIMD_ACCURACY_BITS_SINGLE < GMX_SIMD_ACCURACY_BITS_DOUBLE)
1847 lu0 = rsqrtIter(lu0, x0);
1848 lu1 = rsqrtIter(lu1, x1);
1850 # if (GMX_SIMD_ACCURACY_BITS_SINGLE * 2 < GMX_SIMD_ACCURACY_BITS_DOUBLE)
1851 lu0 = rsqrtIter(lu0, x0);
1852 lu1 = rsqrtIter(lu1, x1);
1857 *out0 = invsqrt(x0);
1858 *out1 = invsqrt(x1);
1862 # if !GMX_SIMD_HAVE_NATIVE_RCP_ITER_DOUBLE
1863 /*! \brief Perform one Newton-Raphson iteration to improve 1/x for SIMD double.
1865 * This is a low-level routine that should only be used by SIMD math routine
1866 * that evaluates the reciprocal.
1868 * \param lu Approximation of 1/x, typically obtained from lookup.
1869 * \param x The reference (starting) value x for which we want 1/x.
1870 * \return An improved approximation with roughly twice as many bits of accuracy.
1872 static inline SimdDouble gmx_simdcall rcpIter(SimdDouble lu, SimdDouble x)
1874 return lu * fnma(lu, x, SimdDouble(2.0));
1878 /*! \brief Calculate 1/x for SIMD double.
1880 * \param x Argument with magnitude larger than GMX_FLOAT_MIN and smaller than
1881 * GMX_FLOAT_MAX, i.e. within the range of single precision.
1882 * For the single precision implementation this is obviously always
1883 * true for positive values, but for double precision it adds an
1884 * extra restriction since the first lookup step might have to be
1885 * performed in single precision on some architectures. Note that the
1886 * responsibility for checking falls on you - this routine does not
1889 * \return 1/x. Result is undefined if your argument was invalid.
1891 static inline SimdDouble gmx_simdcall inv(SimdDouble x)
1893 SimdDouble lu = rcp(x);
1894 # if (GMX_SIMD_RCP_BITS < GMX_SIMD_ACCURACY_BITS_DOUBLE)
1895 lu = rcpIter(lu, x);
1897 # if (GMX_SIMD_RCP_BITS * 2 < GMX_SIMD_ACCURACY_BITS_DOUBLE)
1898 lu = rcpIter(lu, x);
1900 # if (GMX_SIMD_RCP_BITS * 4 < GMX_SIMD_ACCURACY_BITS_DOUBLE)
1901 lu = rcpIter(lu, x);
1903 # if (GMX_SIMD_RCP_BITS * 8 < GMX_SIMD_ACCURACY_BITS_DOUBLE)
1904 lu = rcpIter(lu, x);
1909 /*! \brief Division for SIMD doubles
1911 * \param nom Nominator
1912 * \param denom Denominator, with magnitude in range (GMX_FLOAT_MIN,GMX_FLOAT_MAX).
1913 * For single precision this is equivalent to a nonzero argument,
1914 * but in double precision it adds an extra restriction since
1915 * the first lookup step might have to be performed in single
1916 * precision on some architectures. Note that the responsibility
1917 * for checking falls on you - this routine does not check arguments.
1921 * \note This function does not use any masking to avoid problems with
1922 * zero values in the denominator.
1924 static inline SimdDouble gmx_simdcall operator/(SimdDouble nom, SimdDouble denom)
1926 return nom * inv(denom);
1930 /*! \brief Calculate 1/sqrt(x) for masked entries of SIMD double.
1932 * This routine only evaluates 1/sqrt(x) for elements for which mask is true.
1933 * Illegal values in the masked-out elements will not lead to
1934 * floating-point exceptions.
1936 * \param x Argument that must be larger than GMX_FLOAT_MIN and smaller than
1937 * GMX_FLOAT_MAX for masked-in entries.
1938 * See \ref invsqrt for the discussion about argument restrictions.
1940 * \return 1/sqrt(x). Result is undefined if your argument was invalid or
1941 * entry was not masked, and 0.0 for masked-out entries.
1943 static inline SimdDouble maskzInvsqrt(SimdDouble x, SimdDBool m)
1945 SimdDouble lu = maskzRsqrt(x, m);
1946 # if (GMX_SIMD_RSQRT_BITS < GMX_SIMD_ACCURACY_BITS_DOUBLE)
1947 lu = rsqrtIter(lu, x);
1949 # if (GMX_SIMD_RSQRT_BITS * 2 < GMX_SIMD_ACCURACY_BITS_DOUBLE)
1950 lu = rsqrtIter(lu, x);
1952 # if (GMX_SIMD_RSQRT_BITS * 4 < GMX_SIMD_ACCURACY_BITS_DOUBLE)
1953 lu = rsqrtIter(lu, x);
1955 # if (GMX_SIMD_RSQRT_BITS * 8 < GMX_SIMD_ACCURACY_BITS_DOUBLE)
1956 lu = rsqrtIter(lu, x);
1961 /*! \brief Calculate 1/x for SIMD double, masked version.
1963 * \param x Argument with magnitude larger than GMX_FLOAT_MIN and smaller than
1964 * GMX_FLOAT_MAX for masked-in entries.
1965 * See \ref invsqrt for the discussion about argument restrictions.
1967 * \return 1/x for elements where m is true, or 0.0 for masked-out entries.
1969 static inline SimdDouble gmx_simdcall maskzInv(SimdDouble x, SimdDBool m)
1971 SimdDouble lu = maskzRcp(x, m);
1972 # if (GMX_SIMD_RCP_BITS < GMX_SIMD_ACCURACY_BITS_DOUBLE)
1973 lu = rcpIter(lu, x);
1975 # if (GMX_SIMD_RCP_BITS * 2 < GMX_SIMD_ACCURACY_BITS_DOUBLE)
1976 lu = rcpIter(lu, x);
1978 # if (GMX_SIMD_RCP_BITS * 4 < GMX_SIMD_ACCURACY_BITS_DOUBLE)
1979 lu = rcpIter(lu, x);
1981 # if (GMX_SIMD_RCP_BITS * 8 < GMX_SIMD_ACCURACY_BITS_DOUBLE)
1982 lu = rcpIter(lu, x);
1988 /*! \brief Calculate sqrt(x) for SIMD doubles.
1990 * \copydetails sqrt(SimdFloat)
1992 template<MathOptimization opt = MathOptimization::Safe>
1993 static inline SimdDouble gmx_simdcall sqrt(SimdDouble x)
1995 if (opt == MathOptimization::Safe)
1997 // As we might use a float version of rsqrt, we mask out small values
1998 SimdDouble res = maskzInvsqrt(x, SimdDouble(GMX_FLOAT_MIN) < x);
2003 return x * invsqrt(x);
2007 /*! \brief Cube root for SIMD doubles
2009 * \param x Argument to calculate cube root of. Can be negative or zero,
2010 * but NaN or Inf values are not supported. Denormal values will
2011 * be treated as 0.0.
2012 * \return Cube root of x.
2014 static inline SimdDouble gmx_simdcall cbrt(SimdDouble x)
2016 const SimdDouble signBit(GMX_DOUBLE_NEGZERO);
2017 const SimdDouble minDouble(std::numeric_limits<double>::min());
2018 // Bias is 1024-1 = 1023, which is divisible by 3, so no need to change it more.
2019 // To avoid clang warnings about fragile integer division mixed with FP, we let
2020 // the divided value (1023/3=341) be the original constant.
2021 const std::int32_t offsetDiv3(341);
2022 const SimdDouble c6(-0.145263899385486377);
2023 const SimdDouble c5(0.784932344976639262);
2024 const SimdDouble c4(-1.83469277483613086);
2025 const SimdDouble c3(2.44693122563534430);
2026 const SimdDouble c2(-2.11499494167371287);
2027 const SimdDouble c1(1.50819193781584896);
2028 const SimdDouble c0(0.354895765043919860);
2029 const SimdDouble one(1.0);
2030 const SimdDouble two(2.0);
2031 const SimdDouble three(3.0);
2032 const SimdDouble oneThird(1.0 / 3.0);
2033 const SimdDouble cbrt2(1.2599210498948731648);
2034 const SimdDouble sqrCbrt2(1.5874010519681994748);
2036 // See the single precision routines for documentation of the algorithm
2038 SimdDouble xSignBit = x & signBit; // create bit mask where the sign bit is 1 for x elements < 0
2039 SimdDouble xAbs = andNot(signBit, x); // select everthing but the sign bit => abs(x)
2040 SimdDBool xIsNonZero = (minDouble <= xAbs); // treat denormals as 0
2042 SimdDInt32 exponent;
2043 SimdDouble y = frexp(xAbs, &exponent);
2044 SimdDouble z = fma(y, c6, c5);
2050 SimdDouble w = z * z * z;
2051 SimdDouble nom = z * fma(two, y, w);
2052 SimdDouble invDenom = inv(fma(two, w, y));
2054 SimdDouble offsetExp = cvtI2R(exponent) + SimdDouble(static_cast<double>(3 * offsetDiv3) + 0.1);
2055 SimdDouble offsetExpDiv3 =
2056 trunc(offsetExp * oneThird); // important to truncate here to mimic integer division
2057 SimdDInt32 expDiv3 = cvtR2I(offsetExpDiv3 - SimdDouble(static_cast<double>(offsetDiv3)));
2058 SimdDouble remainder = offsetExp - offsetExpDiv3 * three;
2059 SimdDouble factor = blend(one, cbrt2, SimdDouble(0.5) < remainder);
2060 factor = blend(factor, sqrCbrt2, SimdDouble(1.5) < remainder);
2061 SimdDouble fraction = (nom * invDenom * factor) ^ xSignBit;
2062 SimdDouble result = selectByMask(ldexp(fraction, expDiv3), xIsNonZero);
2066 /*! \brief Inverse cube root for SIMD doubles.
2068 * \param x Argument to calculate cube root of. Can be positive or
2069 * negative, but the magnitude cannot be lower than
2070 * the smallest normal number.
2071 * \return Cube root of x. Undefined for values that don't
2072 * fulfill the restriction of abs(x) > minDouble.
2074 static inline SimdDouble gmx_simdcall invcbrt(SimdDouble x)
2076 const SimdDouble signBit(GMX_DOUBLE_NEGZERO);
2077 // Bias is 1024-1 = 1023, which is divisible by 3, so no need to change it more.
2078 // To avoid clang warnings about fragile integer division mixed with FP, we let
2079 // the divided value (1023/3=341) be the original constant.
2080 const std::int32_t offsetDiv3(341);
2081 const SimdDouble c6(-0.145263899385486377);
2082 const SimdDouble c5(0.784932344976639262);
2083 const SimdDouble c4(-1.83469277483613086);
2084 const SimdDouble c3(2.44693122563534430);
2085 const SimdDouble c2(-2.11499494167371287);
2086 const SimdDouble c1(1.50819193781584896);
2087 const SimdDouble c0(0.354895765043919860);
2088 const SimdDouble one(1.0);
2089 const SimdDouble two(2.0);
2090 const SimdDouble three(3.0);
2091 const SimdDouble oneThird(1.0 / 3.0);
2092 const SimdDouble invCbrt2(1.0 / 1.2599210498948731648);
2093 const SimdDouble invSqrCbrt2(1.0F / 1.5874010519681994748);
2095 // See the single precision routines for documentation of the algorithm
2097 SimdDouble xSignBit = x & signBit; // create bit mask where the sign bit is 1 for x elements < 0
2098 SimdDouble xAbs = andNot(signBit, x); // select everthing but the sign bit => abs(x)
2100 SimdDInt32 exponent;
2101 SimdDouble y = frexp(xAbs, &exponent);
2102 SimdDouble z = fma(y, c6, c5);
2108 SimdDouble w = z * z * z;
2109 SimdDouble nom = fma(two, w, y);
2110 SimdDouble invDenom = inv(z * fma(two, y, w));
2111 SimdDouble offsetExp = cvtI2R(exponent) + SimdDouble(static_cast<double>(3 * offsetDiv3) + 0.1);
2112 SimdDouble offsetExpDiv3 =
2113 trunc(offsetExp * oneThird); // important to truncate here to mimic integer division
2114 SimdDInt32 expDiv3 = cvtR2I(SimdDouble(static_cast<double>(offsetDiv3)) - offsetExpDiv3);
2115 SimdDouble remainder = offsetExpDiv3 * three - offsetExp;
2116 SimdDouble factor = blend(one, invCbrt2, remainder < SimdDouble(-0.5));
2117 factor = blend(factor, invSqrCbrt2, remainder < SimdDouble(-1.5));
2118 SimdDouble fraction = (nom * invDenom * factor) ^ xSignBit;
2119 SimdDouble result = ldexp(fraction, expDiv3);
2123 /*! \brief SIMD double log2(x). This is the base-2 logarithm.
2125 * \param x Argument, should be >0.
2126 * \result The base-2 logarithm of x. Undefined if argument is invalid.
2128 static inline SimdDouble gmx_simdcall log2(SimdDouble x)
2130 # if GMX_SIMD_HAVE_NATIVE_LOG_DOUBLE
2131 // Just rescale if native log2() is not present, but log is.
2132 return log(x) * SimdDouble(std::log2(std::exp(1.0)));
2134 const SimdDouble one(1.0);
2135 const SimdDouble two(2.0);
2136 const SimdDouble invsqrt2(1.0 / std::sqrt(2.0));
2137 const SimdDouble CL15(0.2138031565795550370534528);
2138 const SimdDouble CL13(0.2208884091496370882801159);
2139 const SimdDouble CL11(0.2623358279761824340958754);
2140 const SimdDouble CL9(0.3205984930182496084327681);
2141 const SimdDouble CL7(0.4121985864521960363227038);
2142 const SimdDouble CL5(0.5770780163410746954610886);
2143 const SimdDouble CL3(0.9617966939260027547931031);
2144 const SimdDouble CL1(2.885390081777926774009302);
2145 SimdDouble fExp, x2, p;
2149 // For log2(), the argument cannot be 0, so use the faster version of frexp()
2150 x = frexp<MathOptimization::Unsafe>(x, &iExp);
2151 fExp = cvtI2R(iExp);
2154 // Adjust to non-IEEE format for x<1/sqrt(2): exponent -= 1, mantissa *= 2.0
2155 fExp = fExp - selectByMask(one, m);
2156 x = x * blend(one, two, m);
2158 x = (x - one) * inv(x + one);
2161 p = fma(CL15, x2, CL13);
2162 p = fma(p, x2, CL11);
2163 p = fma(p, x2, CL9);
2164 p = fma(p, x2, CL7);
2165 p = fma(p, x2, CL5);
2166 p = fma(p, x2, CL3);
2167 p = fma(p, x2, CL1);
2168 p = fma(p, x, fExp);
2174 # if !GMX_SIMD_HAVE_NATIVE_LOG_DOUBLE
2175 /*! \brief SIMD double log(x). This is the natural logarithm.
2177 * \param x Argument, should be >0.
2178 * \result The natural logarithm of x. Undefined if argument is invalid.
2180 static inline SimdDouble gmx_simdcall log(SimdDouble x)
2182 const SimdDouble one(1.0);
2183 const SimdDouble two(2.0);
2184 const SimdDouble invsqrt2(1.0 / std::sqrt(2.0));
2185 const SimdDouble corr(0.693147180559945286226764);
2186 const SimdDouble CL15(0.148197055177935105296783);
2187 const SimdDouble CL13(0.153108178020442575739679);
2188 const SimdDouble CL11(0.181837339521549679055568);
2189 const SimdDouble CL9(0.22222194152736701733275);
2190 const SimdDouble CL7(0.285714288030134544449368);
2191 const SimdDouble CL5(0.399999999989941956712869);
2192 const SimdDouble CL3(0.666666666666685503450651);
2193 const SimdDouble CL1(2.0);
2194 SimdDouble fExp, x2, p;
2198 // For log(), the argument cannot be 0, so use the faster version of frexp()
2199 x = frexp<MathOptimization::Unsafe>(x, &iExp);
2200 fExp = cvtI2R(iExp);
2203 // Adjust to non-IEEE format for x<1/sqrt(2): exponent -= 1, mantissa *= 2.0
2204 fExp = fExp - selectByMask(one, m);
2205 x = x * blend(one, two, m);
2207 x = (x - one) * inv(x + one);
2210 p = fma(CL15, x2, CL13);
2211 p = fma(p, x2, CL11);
2212 p = fma(p, x2, CL9);
2213 p = fma(p, x2, CL7);
2214 p = fma(p, x2, CL5);
2215 p = fma(p, x2, CL3);
2216 p = fma(p, x2, CL1);
2217 p = fma(p, x, corr * fExp);
2223 # if !GMX_SIMD_HAVE_NATIVE_EXP2_DOUBLE
2224 /*! \brief SIMD double 2^x.
2226 * \copydetails exp2(SimdFloat)
2228 template<MathOptimization opt = MathOptimization::Safe>
2229 static inline SimdDouble gmx_simdcall exp2(SimdDouble x)
2231 const SimdDouble CE11(4.435280790452730022081181e-10);
2232 const SimdDouble CE10(7.074105630863314448024247e-09);
2233 const SimdDouble CE9(1.017819803432096698472621e-07);
2234 const SimdDouble CE8(1.321543308956718799557863e-06);
2235 const SimdDouble CE7(0.00001525273348995851746990884);
2236 const SimdDouble CE6(0.0001540353046251466849082632);
2237 const SimdDouble CE5(0.001333355814678995257307880);
2238 const SimdDouble CE4(0.009618129107588335039176502);
2239 const SimdDouble CE3(0.05550410866481992147457793);
2240 const SimdDouble CE2(0.2402265069591015620470894);
2241 const SimdDouble CE1(0.6931471805599453304615075);
2242 const SimdDouble one(1.0);
2245 SimdDouble fexppart;
2248 // Large negative values are valid arguments to exp2(), so there are two
2249 // things we need to account for:
2250 // 1. When the exponents reaches -1023, the (biased) exponent field will be
2251 // zero and we can no longer multiply with it. There are special IEEE
2252 // formats to handle this range, but for now we have to accept that
2253 // we cannot handle those arguments. If input value becomes even more
2254 // negative, it will start to loop and we would end up with invalid
2255 // exponents. Thus, we need to limit or mask this.
2256 // 2. For VERY large negative values, we will have problems that the
2257 // subtraction to get the fractional part loses accuracy, and then we
2258 // can end up with overflows in the polynomial.
2260 // For now, we handle this by forwarding the math optimization setting to
2261 // ldexp, where the routine will return zero for very small arguments.
2263 // However, before doing that we need to make sure we do not call cvtR2I
2264 // with an argument that is so negative it cannot be converted to an integer.
2265 if (opt == MathOptimization::Safe)
2267 x = max(x, SimdDouble(std::numeric_limits<std::int32_t>::lowest()));
2270 fexppart = ldexp<opt>(one, cvtR2I(x));
2274 p = fma(CE11, x, CE10);
2290 # if !GMX_SIMD_HAVE_NATIVE_EXP_DOUBLE
2291 /*! \brief SIMD double exp(x).
2293 * \copydetails exp(SimdFloat)
2295 template<MathOptimization opt = MathOptimization::Safe>
2296 static inline SimdDouble gmx_simdcall exp(SimdDouble x)
2298 const SimdDouble argscale(1.44269504088896340735992468100);
2299 const SimdDouble invargscale0(-0.69314718055966295651160180568695068359375);
2300 const SimdDouble invargscale1(-2.8235290563031577122588448175013436025525412068e-13);
2301 const SimdDouble CE12(2.078375306791423699350304e-09);
2302 const SimdDouble CE11(2.518173854179933105218635e-08);
2303 const SimdDouble CE10(2.755842049600488770111608e-07);
2304 const SimdDouble CE9(2.755691815216689746619849e-06);
2305 const SimdDouble CE8(2.480158383706245033920920e-05);
2306 const SimdDouble CE7(0.0001984127043518048611841321);
2307 const SimdDouble CE6(0.001388888889360258341755930);
2308 const SimdDouble CE5(0.008333333332907368102819109);
2309 const SimdDouble CE4(0.04166666666663836745814631);
2310 const SimdDouble CE3(0.1666666666666796929434570);
2311 const SimdDouble CE2(0.5);
2312 const SimdDouble one(1.0);
2313 SimdDouble fexppart;
2317 // Large negative values are valid arguments to exp2(), so there are two
2318 // things we need to account for:
2319 // 1. When the exponents reaches -1023, the (biased) exponent field will be
2320 // zero and we can no longer multiply with it. There are special IEEE
2321 // formats to handle this range, but for now we have to accept that
2322 // we cannot handle those arguments. If input value becomes even more
2323 // negative, it will start to loop and we would end up with invalid
2324 // exponents. Thus, we need to limit or mask this.
2325 // 2. For VERY large negative values, we will have problems that the
2326 // subtraction to get the fractional part loses accuracy, and then we
2327 // can end up with overflows in the polynomial.
2329 // For now, we handle this by forwarding the math optimization setting to
2330 // ldexp, where the routine will return zero for very small arguments.
2332 // However, before doing that we need to make sure we do not call cvtR2I
2333 // with an argument that is so negative it cannot be converted to an integer
2334 // after being multiplied by argscale.
2336 if (opt == MathOptimization::Safe)
2338 x = max(x, SimdDouble(std::numeric_limits<std::int32_t>::lowest()) / argscale);
2343 fexppart = ldexp<opt>(one, cvtR2I(y));
2346 // Extended precision arithmetics
2347 x = fma(invargscale0, intpart, x);
2348 x = fma(invargscale1, intpart, x);
2350 p = fma(CE12, x, CE11);
2351 p = fma(p, x, CE10);
2360 p = fma(p, x * x, x);
2361 # if GMX_SIMD_HAVE_FMA
2362 x = fma(p, fexppart, fexppart);
2364 x = (p + one) * fexppart;
2371 /*! \brief SIMD double pow(x,y)
2373 * This returns x^y for SIMD values.
2375 * \tparam opt If this is changed from the default (safe) into the unsafe
2376 * option, there are no guarantees about correct results for x==0.
2380 * \param y exponent.
2382 * \result x^y. Overflowing arguments are likely to either return 0 or inf,
2383 * depending on the underlying implementation. If unsafe optimizations
2384 * are enabled, this is also true for x==0.
2386 * \warning You cannot rely on this implementation returning inf for arguments
2387 * that cause overflow. If you have some very large
2388 * values and need to rely on getting a valid numerical output,
2389 * take the minimum of your variable and the largest valid argument
2390 * before calling this routine.
2392 template<MathOptimization opt = MathOptimization::Safe>
2393 static inline SimdDouble gmx_simdcall pow(SimdDouble x, SimdDouble y)
2397 if (opt == MathOptimization::Safe)
2399 xcorr = max(x, SimdDouble(std::numeric_limits<double>::min()));
2406 SimdDouble result = exp2<opt>(y * log2(xcorr));
2408 if (opt == MathOptimization::Safe)
2410 // if x==0 and y>0 we explicitly set the result to 0.0
2411 // For any x with y==0, the result will already be 1.0 since we multiply by y (0.0) and call exp().
2412 result = blend(result, setZero(), x == setZero() && setZero() < y);
2419 /*! \brief SIMD double erf(x).
2421 * \param x The value to calculate erf(x) for.
2424 * This routine achieves very close to full precision, but we do not care about
2425 * the last bit or the subnormal result range.
2427 static inline SimdDouble gmx_simdcall erf(SimdDouble x)
2429 // Coefficients for minimax approximation of erf(x)=x*(CAoffset + P(x^2)/Q(x^2)) in range [-0.75,0.75]
2430 const SimdDouble CAP4(-0.431780540597889301512e-4);
2431 const SimdDouble CAP3(-0.00578562306260059236059);
2432 const SimdDouble CAP2(-0.028593586920219752446);
2433 const SimdDouble CAP1(-0.315924962948621698209);
2434 const SimdDouble CAP0(0.14952975608477029151);
2436 const SimdDouble CAQ5(-0.374089300177174709737e-5);
2437 const SimdDouble CAQ4(0.00015126584532155383535);
2438 const SimdDouble CAQ3(0.00536692680669480725423);
2439 const SimdDouble CAQ2(0.0668686825594046122636);
2440 const SimdDouble CAQ1(0.402604990869284362773);
2442 const SimdDouble CAoffset(0.9788494110107421875);
2444 // Coefficients for minimax approximation of erfc(x)=exp(-x^2)*x*(P(x-1)/Q(x-1)) in range [1.0,4.5]
2445 const SimdDouble CBP6(2.49650423685462752497647637088e-10);
2446 const SimdDouble CBP5(0.00119770193298159629350136085658);
2447 const SimdDouble CBP4(0.0164944422378370965881008942733);
2448 const SimdDouble CBP3(0.0984581468691775932063932439252);
2449 const SimdDouble CBP2(0.317364595806937763843589437418);
2450 const SimdDouble CBP1(0.554167062641455850932670067075);
2451 const SimdDouble CBP0(0.427583576155807163756925301060);
2452 const SimdDouble CBQ7(0.00212288829699830145976198384930);
2453 const SimdDouble CBQ6(0.0334810979522685300554606393425);
2454 const SimdDouble CBQ5(0.2361713785181450957579508850717);
2455 const SimdDouble CBQ4(0.955364736493055670530981883072);
2456 const SimdDouble CBQ3(2.36815675631420037315349279199);
2457 const SimdDouble CBQ2(3.55261649184083035537184223542);
2458 const SimdDouble CBQ1(2.93501136050160872574376997993);
2461 // Coefficients for minimax approximation of erfc(x)=exp(-x^2)/x*(P(1/x)/Q(1/x)) in range [4.5,inf]
2462 const SimdDouble CCP6(-2.8175401114513378771);
2463 const SimdDouble CCP5(-3.22729451764143718517);
2464 const SimdDouble CCP4(-2.5518551727311523996);
2465 const SimdDouble CCP3(-0.687717681153649930619);
2466 const SimdDouble CCP2(-0.212652252872804219852);
2467 const SimdDouble CCP1(0.0175389834052493308818);
2468 const SimdDouble CCP0(0.00628057170626964891937);
2470 const SimdDouble CCQ6(5.48409182238641741584);
2471 const SimdDouble CCQ5(13.5064170191802889145);
2472 const SimdDouble CCQ4(22.9367376522880577224);
2473 const SimdDouble CCQ3(15.930646027911794143);
2474 const SimdDouble CCQ2(11.0567237927800161565);
2475 const SimdDouble CCQ1(2.79257750980575282228);
2477 const SimdDouble CCoffset(0.5579090118408203125);
2479 const SimdDouble one(1.0);
2480 const SimdDouble two(2.0);
2481 const SimdDouble minFloat(std::numeric_limits<float>::min());
2483 SimdDouble xabs, x2, x4, t, t2, w, w2;
2484 SimdDouble PolyAP0, PolyAP1, PolyAQ0, PolyAQ1;
2485 SimdDouble PolyBP0, PolyBP1, PolyBQ0, PolyBQ1;
2486 SimdDouble PolyCP0, PolyCP1, PolyCQ0, PolyCQ1;
2487 SimdDouble res_erf, res_erfcB, res_erfcC, res_erfc, res;
2489 SimdDBool mask, mask_erf, notmask_erf;
2493 mask_erf = (xabs < one);
2494 notmask_erf = (one <= xabs);
2498 PolyAP0 = fma(CAP4, x4, CAP2);
2499 PolyAP1 = fma(CAP3, x4, CAP1);
2500 PolyAP0 = fma(PolyAP0, x4, CAP0);
2501 PolyAP0 = fma(PolyAP1, x2, PolyAP0);
2503 PolyAQ1 = fma(CAQ5, x4, CAQ3);
2504 PolyAQ0 = fma(CAQ4, x4, CAQ2);
2505 PolyAQ1 = fma(PolyAQ1, x4, CAQ1);
2506 PolyAQ0 = fma(PolyAQ0, x4, one);
2507 PolyAQ0 = fma(PolyAQ1, x2, PolyAQ0);
2509 res_erf = PolyAP0 * maskzInv(PolyAQ0, mask_erf && (minFloat <= abs(PolyAQ0)));
2510 res_erf = CAoffset + res_erf;
2511 res_erf = x * res_erf;
2513 // Calculate erfc() in range [1,4.5]
2517 PolyBP0 = fma(CBP6, t2, CBP4);
2518 PolyBP1 = fma(CBP5, t2, CBP3);
2519 PolyBP0 = fma(PolyBP0, t2, CBP2);
2520 PolyBP1 = fma(PolyBP1, t2, CBP1);
2521 PolyBP0 = fma(PolyBP0, t2, CBP0);
2522 PolyBP0 = fma(PolyBP1, t, PolyBP0);
2524 PolyBQ1 = fma(CBQ7, t2, CBQ5);
2525 PolyBQ0 = fma(CBQ6, t2, CBQ4);
2526 PolyBQ1 = fma(PolyBQ1, t2, CBQ3);
2527 PolyBQ0 = fma(PolyBQ0, t2, CBQ2);
2528 PolyBQ1 = fma(PolyBQ1, t2, CBQ1);
2529 PolyBQ0 = fma(PolyBQ0, t2, one);
2530 PolyBQ0 = fma(PolyBQ1, t, PolyBQ0);
2532 // The denominator polynomial can be zero outside the range
2533 res_erfcB = PolyBP0 * maskzInv(PolyBQ0, notmask_erf && (minFloat <= abs(PolyBQ0)));
2535 res_erfcB = res_erfcB * xabs;
2537 // Calculate erfc() in range [4.5,inf]
2538 w = maskzInv(xabs, notmask_erf && (minFloat <= xabs));
2541 PolyCP0 = fma(CCP6, w2, CCP4);
2542 PolyCP1 = fma(CCP5, w2, CCP3);
2543 PolyCP0 = fma(PolyCP0, w2, CCP2);
2544 PolyCP1 = fma(PolyCP1, w2, CCP1);
2545 PolyCP0 = fma(PolyCP0, w2, CCP0);
2546 PolyCP0 = fma(PolyCP1, w, PolyCP0);
2548 PolyCQ0 = fma(CCQ6, w2, CCQ4);
2549 PolyCQ1 = fma(CCQ5, w2, CCQ3);
2550 PolyCQ0 = fma(PolyCQ0, w2, CCQ2);
2551 PolyCQ1 = fma(PolyCQ1, w2, CCQ1);
2552 PolyCQ0 = fma(PolyCQ0, w2, one);
2553 PolyCQ0 = fma(PolyCQ1, w, PolyCQ0);
2557 // The denominator polynomial can be zero outside the range
2558 res_erfcC = PolyCP0 * maskzInv(PolyCQ0, notmask_erf && (minFloat <= abs(PolyCQ0)));
2559 res_erfcC = res_erfcC + CCoffset;
2560 res_erfcC = res_erfcC * w;
2562 mask = (SimdDouble(4.5) < xabs);
2563 res_erfc = blend(res_erfcB, res_erfcC, mask);
2565 res_erfc = res_erfc * expmx2;
2567 // erfc(x<0) = 2-erfc(|x|)
2568 mask = (x < setZero());
2569 res_erfc = blend(res_erfc, two - res_erfc, mask);
2571 // Select erf() or erfc()
2572 res = blend(one - res_erfc, res_erf, mask_erf);
2577 /*! \brief SIMD double erfc(x).
2579 * \param x The value to calculate erfc(x) for.
2582 * This routine achieves full precision (bar the last bit) over most of the
2583 * input range, but for large arguments where the result is getting close
2584 * to the minimum representable numbers we accept slightly larger errors
2585 * (think results that are in the ballpark of 10^-200 for double)
2586 * since that is not relevant for MD.
2588 static inline SimdDouble gmx_simdcall erfc(SimdDouble x)
2590 // Coefficients for minimax approximation of erf(x)=x*(CAoffset + P(x^2)/Q(x^2)) in range [-0.75,0.75]
2591 const SimdDouble CAP4(-0.431780540597889301512e-4);
2592 const SimdDouble CAP3(-0.00578562306260059236059);
2593 const SimdDouble CAP2(-0.028593586920219752446);
2594 const SimdDouble CAP1(-0.315924962948621698209);
2595 const SimdDouble CAP0(0.14952975608477029151);
2597 const SimdDouble CAQ5(-0.374089300177174709737e-5);
2598 const SimdDouble CAQ4(0.00015126584532155383535);
2599 const SimdDouble CAQ3(0.00536692680669480725423);
2600 const SimdDouble CAQ2(0.0668686825594046122636);
2601 const SimdDouble CAQ1(0.402604990869284362773);
2603 const SimdDouble CAoffset(0.9788494110107421875);
2605 // Coefficients for minimax approximation of erfc(x)=exp(-x^2)*x*(P(x-1)/Q(x-1)) in range [1.0,4.5]
2606 const SimdDouble CBP6(2.49650423685462752497647637088e-10);
2607 const SimdDouble CBP5(0.00119770193298159629350136085658);
2608 const SimdDouble CBP4(0.0164944422378370965881008942733);
2609 const SimdDouble CBP3(0.0984581468691775932063932439252);
2610 const SimdDouble CBP2(0.317364595806937763843589437418);
2611 const SimdDouble CBP1(0.554167062641455850932670067075);
2612 const SimdDouble CBP0(0.427583576155807163756925301060);
2613 const SimdDouble CBQ7(0.00212288829699830145976198384930);
2614 const SimdDouble CBQ6(0.0334810979522685300554606393425);
2615 const SimdDouble CBQ5(0.2361713785181450957579508850717);
2616 const SimdDouble CBQ4(0.955364736493055670530981883072);
2617 const SimdDouble CBQ3(2.36815675631420037315349279199);
2618 const SimdDouble CBQ2(3.55261649184083035537184223542);
2619 const SimdDouble CBQ1(2.93501136050160872574376997993);
2622 // Coefficients for minimax approximation of erfc(x)=exp(-x^2)/x*(P(1/x)/Q(1/x)) in range [4.5,inf]
2623 const SimdDouble CCP6(-2.8175401114513378771);
2624 const SimdDouble CCP5(-3.22729451764143718517);
2625 const SimdDouble CCP4(-2.5518551727311523996);
2626 const SimdDouble CCP3(-0.687717681153649930619);
2627 const SimdDouble CCP2(-0.212652252872804219852);
2628 const SimdDouble CCP1(0.0175389834052493308818);
2629 const SimdDouble CCP0(0.00628057170626964891937);
2631 const SimdDouble CCQ6(5.48409182238641741584);
2632 const SimdDouble CCQ5(13.5064170191802889145);
2633 const SimdDouble CCQ4(22.9367376522880577224);
2634 const SimdDouble CCQ3(15.930646027911794143);
2635 const SimdDouble CCQ2(11.0567237927800161565);
2636 const SimdDouble CCQ1(2.79257750980575282228);
2638 const SimdDouble CCoffset(0.5579090118408203125);
2640 const SimdDouble one(1.0);
2641 const SimdDouble two(2.0);
2642 const SimdDouble minFloat(std::numeric_limits<float>::min());
2644 SimdDouble xabs, x2, x4, t, t2, w, w2;
2645 SimdDouble PolyAP0, PolyAP1, PolyAQ0, PolyAQ1;
2646 SimdDouble PolyBP0, PolyBP1, PolyBQ0, PolyBQ1;
2647 SimdDouble PolyCP0, PolyCP1, PolyCQ0, PolyCQ1;
2648 SimdDouble res_erf, res_erfcB, res_erfcC, res_erfc, res;
2650 SimdDBool mask, mask_erf, notmask_erf;
2654 mask_erf = (xabs < one);
2655 notmask_erf = (one <= xabs);
2659 PolyAP0 = fma(CAP4, x4, CAP2);
2660 PolyAP1 = fma(CAP3, x4, CAP1);
2661 PolyAP0 = fma(PolyAP0, x4, CAP0);
2662 PolyAP0 = fma(PolyAP1, x2, PolyAP0);
2663 PolyAQ1 = fma(CAQ5, x4, CAQ3);
2664 PolyAQ0 = fma(CAQ4, x4, CAQ2);
2665 PolyAQ1 = fma(PolyAQ1, x4, CAQ1);
2666 PolyAQ0 = fma(PolyAQ0, x4, one);
2667 PolyAQ0 = fma(PolyAQ1, x2, PolyAQ0);
2669 res_erf = PolyAP0 * maskzInv(PolyAQ0, mask_erf && (minFloat <= abs(PolyAQ0)));
2670 res_erf = CAoffset + res_erf;
2671 res_erf = x * res_erf;
2673 // Calculate erfc() in range [1,4.5]
2677 PolyBP0 = fma(CBP6, t2, CBP4);
2678 PolyBP1 = fma(CBP5, t2, CBP3);
2679 PolyBP0 = fma(PolyBP0, t2, CBP2);
2680 PolyBP1 = fma(PolyBP1, t2, CBP1);
2681 PolyBP0 = fma(PolyBP0, t2, CBP0);
2682 PolyBP0 = fma(PolyBP1, t, PolyBP0);
2684 PolyBQ1 = fma(CBQ7, t2, CBQ5);
2685 PolyBQ0 = fma(CBQ6, t2, CBQ4);
2686 PolyBQ1 = fma(PolyBQ1, t2, CBQ3);
2687 PolyBQ0 = fma(PolyBQ0, t2, CBQ2);
2688 PolyBQ1 = fma(PolyBQ1, t2, CBQ1);
2689 PolyBQ0 = fma(PolyBQ0, t2, one);
2690 PolyBQ0 = fma(PolyBQ1, t, PolyBQ0);
2692 // The denominator polynomial can be zero outside the range
2693 res_erfcB = PolyBP0 * maskzInv(PolyBQ0, notmask_erf && (minFloat <= abs(PolyBQ0)));
2695 res_erfcB = res_erfcB * xabs;
2697 // Calculate erfc() in range [4.5,inf]
2698 // Note that 1/x can only handle single precision!
2699 w = maskzInv(xabs, minFloat <= xabs);
2702 PolyCP0 = fma(CCP6, w2, CCP4);
2703 PolyCP1 = fma(CCP5, w2, CCP3);
2704 PolyCP0 = fma(PolyCP0, w2, CCP2);
2705 PolyCP1 = fma(PolyCP1, w2, CCP1);
2706 PolyCP0 = fma(PolyCP0, w2, CCP0);
2707 PolyCP0 = fma(PolyCP1, w, PolyCP0);
2709 PolyCQ0 = fma(CCQ6, w2, CCQ4);
2710 PolyCQ1 = fma(CCQ5, w2, CCQ3);
2711 PolyCQ0 = fma(PolyCQ0, w2, CCQ2);
2712 PolyCQ1 = fma(PolyCQ1, w2, CCQ1);
2713 PolyCQ0 = fma(PolyCQ0, w2, one);
2714 PolyCQ0 = fma(PolyCQ1, w, PolyCQ0);
2718 // The denominator polynomial can be zero outside the range
2719 res_erfcC = PolyCP0 * maskzInv(PolyCQ0, notmask_erf && (minFloat <= abs(PolyCQ0)));
2720 res_erfcC = res_erfcC + CCoffset;
2721 res_erfcC = res_erfcC * w;
2723 mask = (SimdDouble(4.5) < xabs);
2724 res_erfc = blend(res_erfcB, res_erfcC, mask);
2726 res_erfc = res_erfc * expmx2;
2728 // erfc(x<0) = 2-erfc(|x|)
2729 mask = (x < setZero());
2730 res_erfc = blend(res_erfc, two - res_erfc, mask);
2732 // Select erf() or erfc()
2733 res = blend(res_erfc, one - res_erf, mask_erf);
2738 /*! \brief SIMD double sin \& cos.
2740 * \param x The argument to evaluate sin/cos for
2741 * \param[out] sinval Sin(x)
2742 * \param[out] cosval Cos(x)
2744 * This version achieves close to machine precision, but for very large
2745 * magnitudes of the argument we inherently begin to lose accuracy due to the
2746 * argument reduction, despite using extended precision arithmetics internally.
2748 static inline void gmx_simdcall sincos(SimdDouble x, SimdDouble* sinval, SimdDouble* cosval)
2750 // Constants to subtract Pi/4*x from y while minimizing precision loss
2751 const SimdDouble argred0(-2 * 0.78539816290140151978);
2752 const SimdDouble argred1(-2 * 4.9604678871439933374e-10);
2753 const SimdDouble argred2(-2 * 1.1258708853173288931e-18);
2754 const SimdDouble argred3(-2 * 1.7607799325916000908e-27);
2755 const SimdDouble two_over_pi(2.0 / M_PI);
2756 const SimdDouble const_sin5(1.58938307283228937328511e-10);
2757 const SimdDouble const_sin4(-2.50506943502539773349318e-08);
2758 const SimdDouble const_sin3(2.75573131776846360512547e-06);
2759 const SimdDouble const_sin2(-0.000198412698278911770864914);
2760 const SimdDouble const_sin1(0.0083333333333191845961746);
2761 const SimdDouble const_sin0(-0.166666666666666130709393);
2763 const SimdDouble const_cos7(-1.13615350239097429531523e-11);
2764 const SimdDouble const_cos6(2.08757471207040055479366e-09);
2765 const SimdDouble const_cos5(-2.75573144028847567498567e-07);
2766 const SimdDouble const_cos4(2.48015872890001867311915e-05);
2767 const SimdDouble const_cos3(-0.00138888888888714019282329);
2768 const SimdDouble const_cos2(0.0416666666666665519592062);
2769 const SimdDouble half(0.5);
2770 const SimdDouble one(1.0);
2771 SimdDouble ssign, csign;
2772 SimdDouble x2, y, z, psin, pcos, sss, ccc;
2774 # if GMX_SIMD_HAVE_DINT32_ARITHMETICS && GMX_SIMD_HAVE_LOGICAL
2775 const SimdDInt32 ione(1);
2776 const SimdDInt32 itwo(2);
2779 z = x * two_over_pi;
2783 mask = cvtIB2B((iy & ione) == setZero());
2784 ssign = selectByMask(SimdDouble(GMX_DOUBLE_NEGZERO), cvtIB2B((iy & itwo) == itwo));
2785 csign = selectByMask(SimdDouble(GMX_DOUBLE_NEGZERO), cvtIB2B(((iy + ione) & itwo) == itwo));
2787 const SimdDouble quarter(0.25);
2788 const SimdDouble minusquarter(-0.25);
2790 SimdDBool m1, m2, m3;
2792 /* The most obvious way to find the arguments quadrant in the unit circle
2793 * to calculate the sign is to use integer arithmetic, but that is not
2794 * present in all SIMD implementations. As an alternative, we have devised a
2795 * pure floating-point algorithm that uses truncation for argument reduction
2796 * so that we get a new value 0<=q<1 over the unit circle, and then
2797 * do floating-point comparisons with fractions. This is likely to be
2798 * slightly slower (~10%) due to the longer latencies of floating-point, so
2799 * we only use it when integer SIMD arithmetic is not present.
2803 // It is critical that half-way cases are rounded down
2804 z = fma(x, two_over_pi, half);
2808 /* z now starts at 0.0 for x=-pi/4 (although neg. values cannot occur), and
2809 * then increased by 1.0 as x increases by 2*Pi, when it resets to 0.0.
2810 * This removes the 2*Pi periodicity without using any integer arithmetic.
2811 * First check if y had the value 2 or 3, set csign if true.
2814 /* If we have logical operations we can work directly on the signbit, which
2815 * saves instructions. Otherwise we need to represent signs as +1.0/-1.0.
2816 * Thus, if you are altering defines to debug alternative code paths, the
2817 * two GMX_SIMD_HAVE_LOGICAL sections in this routine must either both be
2818 * active or inactive - you will get errors if only one is used.
2820 # if GMX_SIMD_HAVE_LOGICAL
2821 ssign = ssign & SimdDouble(GMX_DOUBLE_NEGZERO);
2822 csign = andNot(q, SimdDouble(GMX_DOUBLE_NEGZERO));
2823 ssign = ssign ^ csign;
2825 ssign = copysign(SimdDouble(1.0), ssign);
2826 csign = copysign(SimdDouble(1.0), q);
2828 ssign = ssign * csign; // swap ssign if csign was set.
2830 // Check if y had value 1 or 3 (remember we subtracted 0.5 from q)
2831 m1 = (q < minusquarter);
2832 m2 = (setZero() <= q);
2836 // where mask is FALSE, swap sign.
2837 csign = csign * blend(SimdDouble(-1.0), one, mask);
2839 x = fma(y, argred0, x);
2840 x = fma(y, argred1, x);
2841 x = fma(y, argred2, x);
2842 x = fma(y, argred3, x);
2845 psin = fma(const_sin5, x2, const_sin4);
2846 psin = fma(psin, x2, const_sin3);
2847 psin = fma(psin, x2, const_sin2);
2848 psin = fma(psin, x2, const_sin1);
2849 psin = fma(psin, x2, const_sin0);
2850 psin = fma(psin, x2 * x, x);
2852 pcos = fma(const_cos7, x2, const_cos6);
2853 pcos = fma(pcos, x2, const_cos5);
2854 pcos = fma(pcos, x2, const_cos4);
2855 pcos = fma(pcos, x2, const_cos3);
2856 pcos = fma(pcos, x2, const_cos2);
2857 pcos = fms(pcos, x2, half);
2858 pcos = fma(pcos, x2, one);
2860 sss = blend(pcos, psin, mask);
2861 ccc = blend(psin, pcos, mask);
2862 // See comment for GMX_SIMD_HAVE_LOGICAL section above.
2863 # if GMX_SIMD_HAVE_LOGICAL
2864 *sinval = sss ^ ssign;
2865 *cosval = ccc ^ csign;
2867 *sinval = sss * ssign;
2868 *cosval = ccc * csign;
2872 /*! \brief SIMD double sin(x).
2874 * \param x The argument to evaluate sin for
2877 * \attention Do NOT call both sin & cos if you need both results, since each of them
2878 * will then call \ref sincos and waste a factor 2 in performance.
2880 static inline SimdDouble gmx_simdcall sin(SimdDouble x)
2887 /*! \brief SIMD double cos(x).
2889 * \param x The argument to evaluate cos for
2892 * \attention Do NOT call both sin & cos if you need both results, since each of them
2893 * will then call \ref sincos and waste a factor 2 in performance.
2895 static inline SimdDouble gmx_simdcall cos(SimdDouble x)
2902 /*! \brief SIMD double tan(x).
2904 * \param x The argument to evaluate tan for
2907 static inline SimdDouble gmx_simdcall tan(SimdDouble x)
2909 const SimdDouble argred0(-2 * 0.78539816290140151978);
2910 const SimdDouble argred1(-2 * 4.9604678871439933374e-10);
2911 const SimdDouble argred2(-2 * 1.1258708853173288931e-18);
2912 const SimdDouble argred3(-2 * 1.7607799325916000908e-27);
2913 const SimdDouble two_over_pi(2.0 / M_PI);
2914 const SimdDouble CT15(1.01419718511083373224408e-05);
2915 const SimdDouble CT14(-2.59519791585924697698614e-05);
2916 const SimdDouble CT13(5.23388081915899855325186e-05);
2917 const SimdDouble CT12(-3.05033014433946488225616e-05);
2918 const SimdDouble CT11(7.14707504084242744267497e-05);
2919 const SimdDouble CT10(8.09674518280159187045078e-05);
2920 const SimdDouble CT9(0.000244884931879331847054404);
2921 const SimdDouble CT8(0.000588505168743587154904506);
2922 const SimdDouble CT7(0.00145612788922812427978848);
2923 const SimdDouble CT6(0.00359208743836906619142924);
2924 const SimdDouble CT5(0.00886323944362401618113356);
2925 const SimdDouble CT4(0.0218694882853846389592078);
2926 const SimdDouble CT3(0.0539682539781298417636002);
2927 const SimdDouble CT2(0.133333333333125941821962);
2928 const SimdDouble CT1(0.333333333333334980164153);
2929 const SimdDouble minFloat(std::numeric_limits<float>::min());
2931 SimdDouble x2, p, y, z;
2934 # if GMX_SIMD_HAVE_DINT32_ARITHMETICS && GMX_SIMD_HAVE_LOGICAL
2938 z = x * two_over_pi;
2941 m = cvtIB2B((iy & ione) == ione);
2943 x = fma(y, argred0, x);
2944 x = fma(y, argred1, x);
2945 x = fma(y, argred2, x);
2946 x = fma(y, argred3, x);
2947 x = selectByMask(SimdDouble(GMX_DOUBLE_NEGZERO), m) ^ x;
2949 const SimdDouble quarter(0.25);
2950 const SimdDouble half(0.5);
2951 const SimdDouble threequarter(0.75);
2952 const SimdDouble minFloat(std::numeric_limits<float>::min());
2954 SimdDBool m1, m2, m3;
2957 z = fma(w, two_over_pi, half);
2961 m1 = (quarter <= q);
2963 m3 = (threequarter <= q);
2966 w = fma(y, argred0, w);
2967 w = fma(y, argred1, w);
2968 w = fma(y, argred2, w);
2969 w = fma(y, argred3, w);
2971 w = blend(w, -w, m);
2972 x = w * copysign(SimdDouble(1.0), x);
2975 p = fma(CT15, x2, CT14);
2976 p = fma(p, x2, CT13);
2977 p = fma(p, x2, CT12);
2978 p = fma(p, x2, CT11);
2979 p = fma(p, x2, CT10);
2980 p = fma(p, x2, CT9);
2981 p = fma(p, x2, CT8);
2982 p = fma(p, x2, CT7);
2983 p = fma(p, x2, CT6);
2984 p = fma(p, x2, CT5);
2985 p = fma(p, x2, CT4);
2986 p = fma(p, x2, CT3);
2987 p = fma(p, x2, CT2);
2988 p = fma(p, x2, CT1);
2989 p = fma(x2, p * x, x);
2991 p = blend(p, maskzInv(p, m && (minFloat < abs(p))), m);
2995 /*! \brief SIMD double asin(x).
2997 * \param x The argument to evaluate asin for
3000 static inline SimdDouble gmx_simdcall asin(SimdDouble x)
3002 // Same algorithm as cephes library
3003 const SimdDouble limit1(0.625);
3004 const SimdDouble limit2(1e-8);
3005 const SimdDouble one(1.0);
3006 const SimdDouble quarterpi(M_PI / 4.0);
3007 const SimdDouble morebits(6.123233995736765886130e-17);
3009 const SimdDouble P5(4.253011369004428248960e-3);
3010 const SimdDouble P4(-6.019598008014123785661e-1);
3011 const SimdDouble P3(5.444622390564711410273e0);
3012 const SimdDouble P2(-1.626247967210700244449e1);
3013 const SimdDouble P1(1.956261983317594739197e1);
3014 const SimdDouble P0(-8.198089802484824371615e0);
3016 const SimdDouble Q4(-1.474091372988853791896e1);
3017 const SimdDouble Q3(7.049610280856842141659e1);
3018 const SimdDouble Q2(-1.471791292232726029859e2);
3019 const SimdDouble Q1(1.395105614657485689735e2);
3020 const SimdDouble Q0(-4.918853881490881290097e1);
3022 const SimdDouble R4(2.967721961301243206100e-3);
3023 const SimdDouble R3(-5.634242780008963776856e-1);
3024 const SimdDouble R2(6.968710824104713396794e0);
3025 const SimdDouble R1(-2.556901049652824852289e1);
3026 const SimdDouble R0(2.853665548261061424989e1);
3028 const SimdDouble S3(-2.194779531642920639778e1);
3029 const SimdDouble S2(1.470656354026814941758e2);
3030 const SimdDouble S1(-3.838770957603691357202e2);
3031 const SimdDouble S0(3.424398657913078477438e2);
3034 SimdDouble zz, ww, z, q, w, zz2, ww2;
3039 SimdDouble nom, denom;
3040 SimdDBool mask, mask2;
3044 mask = (limit1 < xabs);
3052 RA = fma(R4, zz2, R2);
3053 RB = fma(R3, zz2, R1);
3054 RA = fma(RA, zz2, R0);
3055 RA = fma(RB, zz, RA);
3058 SB = fma(S3, zz2, S1);
3060 SA = fma(SA, zz2, S0);
3061 SA = fma(SB, zz, SA);
3064 PA = fma(P5, ww2, P3);
3065 PB = fma(P4, ww2, P2);
3066 PA = fma(PA, ww2, P1);
3067 PB = fma(PB, ww2, P0);
3068 PA = fma(PA, ww, PB);
3071 QB = fma(Q4, ww2, Q2);
3073 QA = fma(QA, ww2, Q1);
3074 QB = fma(QB, ww2, Q0);
3075 QA = fma(QA, ww, QB);
3080 nom = blend(PA, RA, mask);
3081 denom = blend(QA, SA, mask);
3083 mask2 = (limit2 < xabs);
3084 q = nom * maskzInv(denom, mask2);
3089 zz = fms(zz, q, morebits);
3096 z = blend(w, z, mask);
3098 z = blend(xabs, z, mask2);
3105 /*! \brief SIMD double acos(x).
3107 * \param x The argument to evaluate acos for
3110 static inline SimdDouble gmx_simdcall acos(SimdDouble x)
3112 const SimdDouble one(1.0);
3113 const SimdDouble half(0.5);
3114 const SimdDouble quarterpi0(7.85398163397448309616e-1);
3115 const SimdDouble quarterpi1(6.123233995736765886130e-17);
3118 SimdDouble z, z1, z2;
3121 z1 = half * (one - x);
3123 z = blend(x, z1, mask1);
3129 z2 = quarterpi0 - z;
3130 z2 = z2 + quarterpi1;
3131 z2 = z2 + quarterpi0;
3133 z = blend(z2, z1, mask1);
3138 /*! \brief SIMD double asin(x).
3140 * \param x The argument to evaluate atan for
3141 * \result Atan(x), same argument/value range as standard math library.
3143 static inline SimdDouble gmx_simdcall atan(SimdDouble x)
3145 // Same algorithm as cephes library
3146 const SimdDouble limit1(0.66);
3147 const SimdDouble limit2(2.41421356237309504880);
3148 const SimdDouble quarterpi(M_PI / 4.0);
3149 const SimdDouble halfpi(M_PI / 2.0);
3150 const SimdDouble mone(-1.0);
3151 const SimdDouble morebits1(0.5 * 6.123233995736765886130E-17);
3152 const SimdDouble morebits2(6.123233995736765886130E-17);
3154 const SimdDouble P4(-8.750608600031904122785E-1);
3155 const SimdDouble P3(-1.615753718733365076637E1);
3156 const SimdDouble P2(-7.500855792314704667340E1);
3157 const SimdDouble P1(-1.228866684490136173410E2);
3158 const SimdDouble P0(-6.485021904942025371773E1);
3160 const SimdDouble Q4(2.485846490142306297962E1);
3161 const SimdDouble Q3(1.650270098316988542046E2);
3162 const SimdDouble Q2(4.328810604912902668951E2);
3163 const SimdDouble Q1(4.853903996359136964868E2);
3164 const SimdDouble Q0(1.945506571482613964425E2);
3166 SimdDouble y, xabs, t1, t2;
3168 SimdDouble P_A, P_B, Q_A, Q_B;
3169 SimdDBool mask1, mask2;
3173 mask1 = (limit1 < xabs);
3174 mask2 = (limit2 < xabs);
3176 t1 = (xabs + mone) * maskzInv(xabs - mone, mask1);
3177 t2 = mone * maskzInv(xabs, mask2);
3179 y = selectByMask(quarterpi, mask1);
3180 y = blend(y, halfpi, mask2);
3181 xabs = blend(xabs, t1, mask1);
3182 xabs = blend(xabs, t2, mask2);
3187 P_A = fma(P4, z2, P2);
3188 P_B = fma(P3, z2, P1);
3189 P_A = fma(P_A, z2, P0);
3190 P_A = fma(P_B, z, P_A);
3193 Q_B = fma(Q4, z2, Q2);
3195 Q_A = fma(Q_A, z2, Q1);
3196 Q_B = fma(Q_B, z2, Q0);
3197 Q_A = fma(Q_A, z, Q_B);
3201 z = fma(z, xabs, xabs);
3203 t1 = selectByMask(morebits1, mask1);
3204 t1 = blend(t1, morebits2, mask2);
3214 /*! \brief SIMD double atan2(y,x).
3216 * \param y Y component of vector, any quartile
3217 * \param x X component of vector, any quartile
3218 * \result Atan(y,x), same argument/value range as standard math library.
3220 * \note This routine should provide correct results for all finite
3221 * non-zero or positive-zero arguments. However, negative zero arguments will
3222 * be treated as positive zero, which means the return value will deviate from
3223 * the standard math library atan2(y,x) for those cases. That should not be
3224 * of any concern in Gromacs, and in particular it will not affect calculations
3225 * of angles from vectors.
3227 static inline SimdDouble gmx_simdcall atan2(SimdDouble y, SimdDouble x)
3229 const SimdDouble pi(M_PI);
3230 const SimdDouble halfpi(M_PI / 2.0);
3231 const SimdDouble minFloat(std::numeric_limits<float>::min());
3232 SimdDouble xinv, p, aoffset;
3233 SimdDBool mask_xnz, mask_ynz, mask_xlt0, mask_ylt0;
3235 mask_xnz = x != setZero();
3236 mask_ynz = y != setZero();
3237 mask_xlt0 = (x < setZero());
3238 mask_ylt0 = (y < setZero());
3240 aoffset = selectByNotMask(halfpi, mask_xnz);
3241 aoffset = selectByMask(aoffset, mask_ynz);
3243 aoffset = blend(aoffset, pi, mask_xlt0);
3244 aoffset = blend(aoffset, -aoffset, mask_ylt0);
3246 xinv = maskzInv(x, mask_xnz && (minFloat <= abs(x)));
3255 /*! \brief Calculate the force correction due to PME analytically in SIMD double.
3257 * \param z2 This should be the value \f$(r \beta)^2\f$, where r is your
3258 * interaction distance and beta the ewald splitting parameters.
3259 * \result Correction factor to coulomb force.
3261 * This routine is meant to enable analytical evaluation of the
3262 * direct-space PME electrostatic force to avoid tables. For details, see the
3263 * single precision function.
3265 static inline SimdDouble gmx_simdcall pmeForceCorrection(SimdDouble z2)
3267 const SimdDouble FN10(-8.0072854618360083154e-14);
3268 const SimdDouble FN9(1.1859116242260148027e-11);
3269 const SimdDouble FN8(-8.1490406329798423616e-10);
3270 const SimdDouble FN7(3.4404793543907847655e-8);
3271 const SimdDouble FN6(-9.9471420832602741006e-7);
3272 const SimdDouble FN5(0.000020740315999115847456);
3273 const SimdDouble FN4(-0.00031991745139313364005);
3274 const SimdDouble FN3(0.0035074449373659008203);
3275 const SimdDouble FN2(-0.031750380176100813405);
3276 const SimdDouble FN1(0.13884101728898463426);
3277 const SimdDouble FN0(-0.75225277815249618847);
3279 const SimdDouble FD5(0.000016009278224355026701);
3280 const SimdDouble FD4(0.00051055686934806966046);
3281 const SimdDouble FD3(0.0081803507497974289008);
3282 const SimdDouble FD2(0.077181146026670287235);
3283 const SimdDouble FD1(0.41543303143712535988);
3284 const SimdDouble FD0(1.0);
3287 SimdDouble polyFN0, polyFN1, polyFD0, polyFD1;
3291 polyFD1 = fma(FD5, z4, FD3);
3292 polyFD1 = fma(polyFD1, z4, FD1);
3293 polyFD1 = polyFD1 * z2;
3294 polyFD0 = fma(FD4, z4, FD2);
3295 polyFD0 = fma(polyFD0, z4, FD0);
3296 polyFD0 = polyFD0 + polyFD1;
3298 polyFD0 = inv(polyFD0);
3300 polyFN0 = fma(FN10, z4, FN8);
3301 polyFN0 = fma(polyFN0, z4, FN6);
3302 polyFN0 = fma(polyFN0, z4, FN4);
3303 polyFN0 = fma(polyFN0, z4, FN2);
3304 polyFN0 = fma(polyFN0, z4, FN0);
3305 polyFN1 = fma(FN9, z4, FN7);
3306 polyFN1 = fma(polyFN1, z4, FN5);
3307 polyFN1 = fma(polyFN1, z4, FN3);
3308 polyFN1 = fma(polyFN1, z4, FN1);
3309 polyFN0 = fma(polyFN1, z2, polyFN0);
3312 return polyFN0 * polyFD0;
3316 /*! \brief Calculate the potential correction due to PME analytically in SIMD double.
3318 * \param z2 This should be the value \f$(r \beta)^2\f$, where r is your
3319 * interaction distance and beta the ewald splitting parameters.
3320 * \result Correction factor to coulomb force.
3322 * This routine is meant to enable analytical evaluation of the
3323 * direct-space PME electrostatic potential to avoid tables. For details, see the
3324 * single precision function.
3326 static inline SimdDouble gmx_simdcall pmePotentialCorrection(SimdDouble z2)
3328 const SimdDouble VN9(-9.3723776169321855475e-13);
3329 const SimdDouble VN8(1.2280156762674215741e-10);
3330 const SimdDouble VN7(-7.3562157912251309487e-9);
3331 const SimdDouble VN6(2.6215886208032517509e-7);
3332 const SimdDouble VN5(-4.9532491651265819499e-6);
3333 const SimdDouble VN4(0.00025907400778966060389);
3334 const SimdDouble VN3(0.0010585044856156469792);
3335 const SimdDouble VN2(0.045247661136833092885);
3336 const SimdDouble VN1(0.11643931522926034421);
3337 const SimdDouble VN0(1.1283791671726767970);
3339 const SimdDouble VD5(0.000021784709867336150342);
3340 const SimdDouble VD4(0.00064293662010911388448);
3341 const SimdDouble VD3(0.0096311444822588683504);
3342 const SimdDouble VD2(0.085608012351550627051);
3343 const SimdDouble VD1(0.43652499166614811084);
3344 const SimdDouble VD0(1.0);
3347 SimdDouble polyVN0, polyVN1, polyVD0, polyVD1;
3351 polyVD1 = fma(VD5, z4, VD3);
3352 polyVD0 = fma(VD4, z4, VD2);
3353 polyVD1 = fma(polyVD1, z4, VD1);
3354 polyVD0 = fma(polyVD0, z4, VD0);
3355 polyVD0 = fma(polyVD1, z2, polyVD0);
3357 polyVD0 = inv(polyVD0);
3359 polyVN1 = fma(VN9, z4, VN7);
3360 polyVN0 = fma(VN8, z4, VN6);
3361 polyVN1 = fma(polyVN1, z4, VN5);
3362 polyVN0 = fma(polyVN0, z4, VN4);
3363 polyVN1 = fma(polyVN1, z4, VN3);
3364 polyVN0 = fma(polyVN0, z4, VN2);
3365 polyVN1 = fma(polyVN1, z4, VN1);
3366 polyVN0 = fma(polyVN0, z4, VN0);
3367 polyVN0 = fma(polyVN1, z2, polyVN0);
3369 return polyVN0 * polyVD0;
3375 /*! \name SIMD math functions for double prec. data, single prec. accuracy
3377 * \note In some cases we do not need full double accuracy of individual
3378 * SIMD math functions, although the data is stored in double precision
3379 * SIMD registers. This might be the case for special algorithms, or
3380 * if the architecture does not support single precision.
3381 * Since the full double precision evaluation of math functions
3382 * typically require much more expensive polynomial approximations
3383 * these functions implement the algorithms used in the single precision
3384 * SIMD math functions, but they operate on double precision
3390 /*********************************************************************
3391 * SIMD MATH FUNCTIONS WITH DOUBLE PREC. DATA, SINGLE PREC. ACCURACY *
3392 *********************************************************************/
3394 /*! \brief Calculate 1/sqrt(x) for SIMD double, but in single accuracy.
3396 * \param x Argument that must be larger than GMX_FLOAT_MIN and smaller than
3397 * GMX_FLOAT_MAX, i.e. within the range of single precision.
3398 * For the single precision implementation this is obviously always
3399 * true for positive values, but for double precision it adds an
3400 * extra restriction since the first lookup step might have to be
3401 * performed in single precision on some architectures. Note that the
3402 * responsibility for checking falls on you - this routine does not
3405 * \return 1/sqrt(x). Result is undefined if your argument was invalid.
3407 static inline SimdDouble gmx_simdcall invsqrtSingleAccuracy(SimdDouble x)
3409 SimdDouble lu = rsqrt(x);
3410 # if (GMX_SIMD_RSQRT_BITS < GMX_SIMD_ACCURACY_BITS_SINGLE)
3411 lu = rsqrtIter(lu, x);
3413 # if (GMX_SIMD_RSQRT_BITS * 2 < GMX_SIMD_ACCURACY_BITS_SINGLE)
3414 lu = rsqrtIter(lu, x);
3416 # if (GMX_SIMD_RSQRT_BITS * 4 < GMX_SIMD_ACCURACY_BITS_SINGLE)
3417 lu = rsqrtIter(lu, x);
3422 /*! \brief 1/sqrt(x) for masked-in entries of SIMD double, but in single accuracy.
3424 * This routine only evaluates 1/sqrt(x) for elements for which mask is true.
3425 * Illegal values in the masked-out elements will not lead to
3426 * floating-point exceptions.
3428 * \param x Argument that must be larger than GMX_FLOAT_MIN and smaller than
3429 * GMX_FLOAT_MAX, i.e. within the range of single precision.
3430 * For the single precision implementation this is obviously always
3431 * true for positive values, but for double precision it adds an
3432 * extra restriction since the first lookup step might have to be
3433 * performed in single precision on some architectures. Note that the
3434 * responsibility for checking falls on you - this routine does not
3438 * \return 1/sqrt(x). Result is undefined if your argument was invalid or
3439 * entry was not masked, and 0.0 for masked-out entries.
3441 static inline SimdDouble maskzInvsqrtSingleAccuracy(SimdDouble x, SimdDBool m)
3443 SimdDouble lu = maskzRsqrt(x, m);
3444 # if (GMX_SIMD_RSQRT_BITS < GMX_SIMD_ACCURACY_BITS_SINGLE)
3445 lu = rsqrtIter(lu, x);
3447 # if (GMX_SIMD_RSQRT_BITS * 2 < GMX_SIMD_ACCURACY_BITS_SINGLE)
3448 lu = rsqrtIter(lu, x);
3450 # if (GMX_SIMD_RSQRT_BITS * 4 < GMX_SIMD_ACCURACY_BITS_SINGLE)
3451 lu = rsqrtIter(lu, x);
3456 /*! \brief Calculate 1/sqrt(x) for two SIMD doubles, but single accuracy.
3458 * \param x0 First set of arguments, x0 must be in single range (see below).
3459 * \param x1 Second set of arguments, x1 must be in single range (see below).
3460 * \param[out] out0 Result 1/sqrt(x0)
3461 * \param[out] out1 Result 1/sqrt(x1)
3463 * In particular for double precision we can sometimes calculate square root
3464 * pairs slightly faster by using single precision until the very last step.
3466 * \note Both arguments must be larger than GMX_FLOAT_MIN and smaller than
3467 * GMX_FLOAT_MAX, i.e. within the range of single precision.
3468 * For the single precision implementation this is obviously always
3469 * true for positive values, but for double precision it adds an
3470 * extra restriction since the first lookup step might have to be
3471 * performed in single precision on some architectures. Note that the
3472 * responsibility for checking falls on you - this routine does not
3475 static inline void gmx_simdcall invsqrtPairSingleAccuracy(SimdDouble x0,
3480 # if GMX_SIMD_HAVE_FLOAT && (GMX_SIMD_FLOAT_WIDTH == 2 * GMX_SIMD_DOUBLE_WIDTH) \
3481 && (GMX_SIMD_RSQRT_BITS < 22)
3482 SimdFloat xf = cvtDD2F(x0, x1);
3483 SimdFloat luf = rsqrt(xf);
3484 SimdDouble lu0, lu1;
3485 // Intermediate target is single - mantissa+1 bits
3486 # if (GMX_SIMD_RSQRT_BITS < GMX_SIMD_ACCURACY_BITS_SINGLE)
3487 luf = rsqrtIter(luf, xf);
3489 # if (GMX_SIMD_RSQRT_BITS * 2 < GMX_SIMD_ACCURACY_BITS_SINGLE)
3490 luf = rsqrtIter(luf, xf);
3492 # if (GMX_SIMD_RSQRT_BITS * 4 < GMX_SIMD_ACCURACY_BITS_SINGLE)
3493 luf = rsqrtIter(luf, xf);
3495 cvtF2DD(luf, &lu0, &lu1);
3496 // We now have single-precision accuracy values in lu0/lu1
3500 *out0 = invsqrtSingleAccuracy(x0);
3501 *out1 = invsqrtSingleAccuracy(x1);
3505 /*! \brief Calculate 1/x for SIMD double, but in single accuracy.
3507 * \param x Argument with magnitude larger than GMX_FLOAT_MIN and smaller than
3508 * GMX_FLOAT_MAX, i.e. within the range of single precision.
3509 * For the single precision implementation this is obviously always
3510 * true for positive values, but for double precision it adds an
3511 * extra restriction since the first lookup step might have to be
3512 * performed in single precision on some architectures. Note that the
3513 * responsibility for checking falls on you - this routine does not
3516 * \return 1/x. Result is undefined if your argument was invalid.
3518 static inline SimdDouble gmx_simdcall invSingleAccuracy(SimdDouble x)
3520 SimdDouble lu = rcp(x);
3521 # if (GMX_SIMD_RCP_BITS < GMX_SIMD_ACCURACY_BITS_SINGLE)
3522 lu = rcpIter(lu, x);
3524 # if (GMX_SIMD_RCP_BITS * 2 < GMX_SIMD_ACCURACY_BITS_SINGLE)
3525 lu = rcpIter(lu, x);
3527 # if (GMX_SIMD_RCP_BITS * 4 < GMX_SIMD_ACCURACY_BITS_SINGLE)
3528 lu = rcpIter(lu, x);
3533 /*! \brief 1/x for masked entries of SIMD double, single accuracy.
3535 * \param x Argument with magnitude larger than GMX_FLOAT_MIN and smaller than
3536 * GMX_FLOAT_MAX, i.e. within the range of single precision.
3537 * For the single precision implementation this is obviously always
3538 * true for positive values, but for double precision it adds an
3539 * extra restriction since the first lookup step might have to be
3540 * performed in single precision on some architectures. Note that the
3541 * responsibility for checking falls on you - this routine does not
3545 * \return 1/x for elements where m is true, or 0.0 for masked-out entries.
3547 static inline SimdDouble gmx_simdcall maskzInvSingleAccuracy(SimdDouble x, SimdDBool m)
3549 SimdDouble lu = maskzRcp(x, m);
3550 # if (GMX_SIMD_RCP_BITS < GMX_SIMD_ACCURACY_BITS_SINGLE)
3551 lu = rcpIter(lu, x);
3553 # if (GMX_SIMD_RCP_BITS * 2 < GMX_SIMD_ACCURACY_BITS_SINGLE)
3554 lu = rcpIter(lu, x);
3556 # if (GMX_SIMD_RCP_BITS * 4 < GMX_SIMD_ACCURACY_BITS_SINGLE)
3557 lu = rcpIter(lu, x);
3563 /*! \brief Calculate sqrt(x) (correct for 0.0) for SIMD double, with single accuracy.
3565 * \copydetails sqrt(SimdFloat)
3567 template<MathOptimization opt = MathOptimization::Safe>
3568 static inline SimdDouble gmx_simdcall sqrtSingleAccuracy(SimdDouble x)
3570 if (opt == MathOptimization::Safe)
3572 SimdDouble res = maskzInvsqrt(x, SimdDouble(GMX_FLOAT_MIN) < x);
3577 return x * invsqrtSingleAccuracy(x);
3581 /*! \brief Cube root for SIMD doubles, single accuracy.
3583 * \param x Argument to calculate cube root of. Can be negative or zero,
3584 * but NaN or Inf values are not supported. Denormal values will
3585 * be treated as 0.0.
3586 * \return Cube root of x.
3588 static inline SimdDouble gmx_simdcall cbrtSingleAccuracy(SimdDouble x)
3590 const SimdDouble signBit(GMX_DOUBLE_NEGZERO);
3591 const SimdDouble minDouble(std::numeric_limits<double>::min());
3592 // Bias is 1024-1 = 1023, which is divisible by 3, so no need to change it more.
3593 // Use the divided value as original constant to avoid division warnings.
3594 const std::int32_t offsetDiv3(341);
3595 const SimdDouble c2(-0.191502161678719066);
3596 const SimdDouble c1(0.697570460207922770);
3597 const SimdDouble c0(0.492659620528969547);
3598 const SimdDouble one(1.0);
3599 const SimdDouble two(2.0);
3600 const SimdDouble three(3.0);
3601 const SimdDouble oneThird(1.0 / 3.0);
3602 const SimdDouble cbrt2(1.2599210498948731648);
3603 const SimdDouble sqrCbrt2(1.5874010519681994748);
3605 // See the single precision routines for documentation of the algorithm
3607 SimdDouble xSignBit = x & signBit; // create bit mask where the sign bit is 1 for x elements < 0
3608 SimdDouble xAbs = andNot(signBit, x); // select everthing but the sign bit => abs(x)
3609 SimdDBool xIsNonZero = (minDouble <= xAbs); // treat denormals as 0
3611 SimdDInt32 exponent;
3612 SimdDouble y = frexp(xAbs, &exponent);
3613 SimdDouble z = fma(fma(y, c2, c1), y, c0);
3614 SimdDouble w = z * z * z;
3615 SimdDouble nom = z * fma(two, y, w);
3616 SimdDouble invDenom = inv(fma(two, w, y));
3618 SimdDouble offsetExp = cvtI2R(exponent) + SimdDouble(static_cast<double>(3 * offsetDiv3) + 0.1);
3619 SimdDouble offsetExpDiv3 =
3620 trunc(offsetExp * oneThird); // important to truncate here to mimic integer division
3621 SimdDInt32 expDiv3 = cvtR2I(offsetExpDiv3 - SimdDouble(static_cast<double>(offsetDiv3)));
3622 SimdDouble remainder = offsetExp - offsetExpDiv3 * three;
3623 SimdDouble factor = blend(one, cbrt2, SimdDouble(0.5) < remainder);
3624 factor = blend(factor, sqrCbrt2, SimdDouble(1.5) < remainder);
3625 SimdDouble fraction = (nom * invDenom * factor) ^ xSignBit;
3626 SimdDouble result = selectByMask(ldexp(fraction, expDiv3), xIsNonZero);
3630 /*! \brief Inverse cube root for SIMD doubles, single accuracy.
3632 * \param x Argument to calculate cube root of. Can be positive or
3633 * negative, but the magnitude cannot be lower than
3634 * the smallest normal number.
3635 * \return Cube root of x. Undefined for values that don't
3636 * fulfill the restriction of abs(x) > minDouble.
3638 static inline SimdDouble gmx_simdcall invcbrtSingleAccuracy(SimdDouble x)
3640 const SimdDouble signBit(GMX_DOUBLE_NEGZERO);
3641 // Bias is 1024-1 = 1023, which is divisible by 3, so no need to change it more.
3642 // Use the divided value as original constant to avoid division warnings.
3643 const std::int32_t offsetDiv3(341);
3644 const SimdDouble c2(-0.191502161678719066);
3645 const SimdDouble c1(0.697570460207922770);
3646 const SimdDouble c0(0.492659620528969547);
3647 const SimdDouble one(1.0);
3648 const SimdDouble two(2.0);
3649 const SimdDouble three(3.0);
3650 const SimdDouble oneThird(1.0 / 3.0);
3651 const SimdDouble invCbrt2(1.0 / 1.2599210498948731648);
3652 const SimdDouble invSqrCbrt2(1.0F / 1.5874010519681994748);
3654 // See the single precision routines for documentation of the algorithm
3656 SimdDouble xSignBit = x & signBit; // create bit mask where the sign bit is 1 for x elements < 0
3657 SimdDouble xAbs = andNot(signBit, x); // select everthing but the sign bit => abs(x)
3659 SimdDInt32 exponent;
3660 SimdDouble y = frexp(xAbs, &exponent);
3661 SimdDouble z = fma(fma(y, c2, c1), y, c0);
3662 SimdDouble w = z * z * z;
3663 SimdDouble nom = fma(two, w, y);
3664 SimdDouble invDenom = inv(z * fma(two, y, w));
3665 SimdDouble offsetExp = cvtI2R(exponent) + SimdDouble(static_cast<double>(3 * offsetDiv3) + 0.1);
3666 SimdDouble offsetExpDiv3 =
3667 trunc(offsetExp * oneThird); // important to truncate here to mimic integer division
3668 SimdDInt32 expDiv3 = cvtR2I(SimdDouble(static_cast<double>(offsetDiv3)) - offsetExpDiv3);
3669 SimdDouble remainder = offsetExpDiv3 * three - offsetExp;
3670 SimdDouble factor = blend(one, invCbrt2, remainder < SimdDouble(-0.5));
3671 factor = blend(factor, invSqrCbrt2, remainder < SimdDouble(-1.5));
3672 SimdDouble fraction = (nom * invDenom * factor) ^ xSignBit;
3673 SimdDouble result = ldexp(fraction, expDiv3);
3677 /*! \brief SIMD log2(x). Double precision SIMD data, single accuracy.
3679 * \param x Argument, should be >0.
3680 * \result The base 2 logarithm of x. Undefined if argument is invalid.
3682 static inline SimdDouble gmx_simdcall log2SingleAccuracy(SimdDouble x)
3684 # if GMX_SIMD_HAVE_NATIVE_LOG_DOUBLE
3685 return log(x) * SimdDouble(std::log2(std::exp(1.0)));
3687 const SimdDouble one(1.0);
3688 const SimdDouble two(2.0);
3689 const SimdDouble sqrt2(std::sqrt(2.0));
3690 const SimdDouble CL9(0.342149508897807708152F);
3691 const SimdDouble CL7(0.411570606888219447939F);
3692 const SimdDouble CL5(0.577085979152320294183F);
3693 const SimdDouble CL3(0.961796550607099898222F);
3694 const SimdDouble CL1(2.885390081777926774009F);
3695 SimdDouble fexp, x2, p;
3699 // For log2(), the argument cannot be 0, so use the faster version of frexp
3700 x = frexp<MathOptimization::Unsafe>(x, &iexp);
3701 fexp = cvtI2R(iexp);
3704 // Adjust to non-IEEE format for x<sqrt(2): exponent -= 1, mantissa *= 2.0
3705 fexp = fexp - selectByMask(one, mask);
3706 x = x * blend(one, two, mask);
3708 x = (x - one) * invSingleAccuracy(x + one);
3711 p = fma(CL9, x2, CL7);
3712 p = fma(p, x2, CL5);
3713 p = fma(p, x2, CL3);
3714 p = fma(p, x2, CL1);
3715 p = fma(p, x, fexp);
3721 /*! \brief SIMD log(x). Double precision SIMD data, single accuracy.
3723 * \param x Argument, should be >0.
3724 * \result The natural logarithm of x. Undefined if argument is invalid.
3726 static inline SimdDouble gmx_simdcall logSingleAccuracy(SimdDouble x)
3728 # if GMX_SIMD_HAVE_NATIVE_LOG_DOUBLE
3731 const SimdDouble one(1.0);
3732 const SimdDouble two(2.0);
3733 const SimdDouble invsqrt2(1.0 / std::sqrt(2.0));
3734 const SimdDouble corr(0.693147180559945286226764);
3735 const SimdDouble CL9(0.2371599674224853515625);
3736 const SimdDouble CL7(0.285279005765914916992188);
3737 const SimdDouble CL5(0.400005519390106201171875);
3738 const SimdDouble CL3(0.666666567325592041015625);
3739 const SimdDouble CL1(2.0);
3740 SimdDouble fexp, x2, p;
3744 // For log(), the argument cannot be 0, so use the faster version of frexp
3745 x = frexp<MathOptimization::Unsafe>(x, &iexp);
3746 fexp = cvtI2R(iexp);
3748 mask = x < invsqrt2;
3749 // Adjust to non-IEEE format for x<1/sqrt(2): exponent -= 1, mantissa *= 2.0
3750 fexp = fexp - selectByMask(one, mask);
3751 x = x * blend(one, two, mask);
3753 x = (x - one) * invSingleAccuracy(x + one);
3756 p = fma(CL9, x2, CL7);
3757 p = fma(p, x2, CL5);
3758 p = fma(p, x2, CL3);
3759 p = fma(p, x2, CL1);
3760 p = fma(p, x, corr * fexp);
3766 /*! \brief SIMD 2^x. Double precision SIMD, single accuracy.
3768 * \copydetails exp2(SimdFloat)
3770 template<MathOptimization opt = MathOptimization::Safe>
3771 static inline SimdDouble gmx_simdcall exp2SingleAccuracy(SimdDouble x)
3773 # if GMX_SIMD_HAVE_NATIVE_EXP2_DOUBLE
3776 const SimdDouble CC6(0.0001534581200287996416911311);
3777 const SimdDouble CC5(0.001339993121934088894618990);
3778 const SimdDouble CC4(0.009618488957115180159497841);
3779 const SimdDouble CC3(0.05550328776964726865751735);
3780 const SimdDouble CC2(0.2402264689063408646490722);
3781 const SimdDouble CC1(0.6931472057372680777553816);
3782 const SimdDouble one(1.0);
3788 // Large negative values are valid arguments to exp2(), so there are two
3789 // things we need to account for:
3790 // 1. When the exponents reaches -1023, the (biased) exponent field will be
3791 // zero and we can no longer multiply with it. There are special IEEE
3792 // formats to handle this range, but for now we have to accept that
3793 // we cannot handle those arguments. If input value becomes even more
3794 // negative, it will start to loop and we would end up with invalid
3795 // exponents. Thus, we need to limit or mask this.
3796 // 2. For VERY large negative values, we will have problems that the
3797 // subtraction to get the fractional part loses accuracy, and then we
3798 // can end up with overflows in the polynomial.
3800 // For now, we handle this by forwarding the math optimization setting to
3801 // ldexp, where the routine will return zero for very small arguments.
3803 // However, before doing that we need to make sure we do not call cvtR2I
3804 // with an argument that is so negative it cannot be converted to an integer.
3805 if (opt == MathOptimization::Safe)
3807 x = max(x, SimdDouble(std::numeric_limits<std::int32_t>::lowest()));
3814 p = fma(CC6, x, CC5);
3820 x = ldexp<opt>(p, ix);
3827 /*! \brief SIMD exp(x). Double precision SIMD, single accuracy.
3829 * \copydetails exp(SimdFloat)
3831 template<MathOptimization opt = MathOptimization::Safe>
3832 static inline SimdDouble gmx_simdcall expSingleAccuracy(SimdDouble x)
3834 # if GMX_SIMD_HAVE_NATIVE_EXP_DOUBLE
3837 const SimdDouble argscale(1.44269504088896341);
3838 // Lower bound: Clamp args that would lead to an IEEE fp exponent below -1023.
3839 const SimdDouble smallArgLimit(-709.0895657128);
3840 const SimdDouble invargscale(-0.69314718055994528623);
3841 const SimdDouble CC4(0.00136324646882712841033936);
3842 const SimdDouble CC3(0.00836596917361021041870117);
3843 const SimdDouble CC2(0.0416710823774337768554688);
3844 const SimdDouble CC1(0.166665524244308471679688);
3845 const SimdDouble CC0(0.499999850988388061523438);
3846 const SimdDouble one(1.0);
3851 // Large negative values are valid arguments to exp2(), so there are two
3852 // things we need to account for:
3853 // 1. When the exponents reaches -1023, the (biased) exponent field will be
3854 // zero and we can no longer multiply with it. There are special IEEE
3855 // formats to handle this range, but for now we have to accept that
3856 // we cannot handle those arguments. If input value becomes even more
3857 // negative, it will start to loop and we would end up with invalid
3858 // exponents. Thus, we need to limit or mask this.
3859 // 2. For VERY large negative values, we will have problems that the
3860 // subtraction to get the fractional part loses accuracy, and then we
3861 // can end up with overflows in the polynomial.
3863 // For now, we handle this by forwarding the math optimization setting to
3864 // ldexp, where the routine will return zero for very small arguments.
3866 // However, before doing that we need to make sure we do not call cvtR2I
3867 // with an argument that is so negative it cannot be converted to an integer
3868 // after being multiplied by argscale.
3870 if (opt == MathOptimization::Safe)
3872 x = max(x, SimdDouble(std::numeric_limits<std::int32_t>::lowest()) / argscale);
3878 intpart = round(y); // use same rounding algorithm here
3880 // Extended precision arithmetics not needed since
3881 // we have double precision and only need single accuracy.
3882 x = fma(invargscale, intpart, x);
3884 p = fma(CC4, x, CC3);
3888 p = fma(x * x, p, x);
3890 x = ldexp<opt>(p, iy);
3895 /*! \brief SIMD pow(x,y). Double precision SIMD data, single accuracy.
3897 * This returns x^y for SIMD values.
3899 * \tparam opt If this is changed from the default (safe) into the unsafe
3900 * option, there are no guarantees about correct results for x==0.
3904 * \param y exponent.
3906 * \result x^y. Overflowing arguments are likely to either return 0 or inf,
3907 * depending on the underlying implementation. If unsafe optimizations
3908 * are enabled, this is also true for x==0.
3910 * \warning You cannot rely on this implementation returning inf for arguments
3911 * that cause overflow. If you have some very large
3912 * values and need to rely on getting a valid numerical output,
3913 * take the minimum of your variable and the largest valid argument
3914 * before calling this routine.
3916 template<MathOptimization opt = MathOptimization::Safe>
3917 static inline SimdDouble gmx_simdcall powSingleAccuracy(SimdDouble x, SimdDouble y)
3921 if (opt == MathOptimization::Safe)
3923 xcorr = max(x, SimdDouble(std::numeric_limits<double>::min()));
3930 SimdDouble result = exp2SingleAccuracy<opt>(y * log2SingleAccuracy(xcorr));
3932 if (opt == MathOptimization::Safe)
3934 // if x==0 and y>0 we explicitly set the result to 0.0
3935 // For any x with y==0, the result will already be 1.0 since we multiply by y (0.0) and call exp().
3936 result = blend(result, setZero(), x == setZero() && setZero() < y);
3942 /*! \brief SIMD erf(x). Double precision SIMD data, single accuracy.
3944 * \param x The value to calculate erf(x) for.
3947 * This routine achieves very close to single precision, but we do not care about
3948 * the last bit or the subnormal result range.
3950 static inline SimdDouble gmx_simdcall erfSingleAccuracy(SimdDouble x)
3952 // Coefficients for minimax approximation of erf(x)=x*P(x^2) in range [-1,1]
3953 const SimdDouble CA6(7.853861353153693e-5);
3954 const SimdDouble CA5(-8.010193625184903e-4);
3955 const SimdDouble CA4(5.188327685732524e-3);
3956 const SimdDouble CA3(-2.685381193529856e-2);
3957 const SimdDouble CA2(1.128358514861418e-1);
3958 const SimdDouble CA1(-3.761262582423300e-1);
3959 const SimdDouble CA0(1.128379165726710);
3960 // Coefficients for minimax approximation of erfc(x)=Exp(-x^2)*P((1/(x-1))^2) in range [0.67,2]
3961 const SimdDouble CB9(-0.0018629930017603923);
3962 const SimdDouble CB8(0.003909821287598495);
3963 const SimdDouble CB7(-0.0052094582210355615);
3964 const SimdDouble CB6(0.005685614362160572);
3965 const SimdDouble CB5(-0.0025367682853477272);
3966 const SimdDouble CB4(-0.010199799682318782);
3967 const SimdDouble CB3(0.04369575504816542);
3968 const SimdDouble CB2(-0.11884063474674492);
3969 const SimdDouble CB1(0.2732120154030589);
3970 const SimdDouble CB0(0.42758357702025784);
3971 // Coefficients for minimax approximation of erfc(x)=Exp(-x^2)*(1/x)*P((1/x)^2) in range [2,9.19]
3972 const SimdDouble CC10(-0.0445555913112064);
3973 const SimdDouble CC9(0.21376355144663348);
3974 const SimdDouble CC8(-0.3473187200259257);
3975 const SimdDouble CC7(0.016690861551248114);
3976 const SimdDouble CC6(0.7560973182491192);
3977 const SimdDouble CC5(-1.2137903600145787);
3978 const SimdDouble CC4(0.8411872321232948);
3979 const SimdDouble CC3(-0.08670413896296343);
3980 const SimdDouble CC2(-0.27124782687240334);
3981 const SimdDouble CC1(-0.0007502488047806069);
3982 const SimdDouble CC0(0.5642114853803148);
3983 const SimdDouble one(1.0);
3984 const SimdDouble two(2.0);
3986 SimdDouble x2, x4, y;
3987 SimdDouble t, t2, w, w2;
3988 SimdDouble pA0, pA1, pB0, pB1, pC0, pC1;
3990 SimdDouble res_erf, res_erfc, res;
3991 SimdDBool mask, msk_erf;
3997 pA0 = fma(CA6, x4, CA4);
3998 pA1 = fma(CA5, x4, CA3);
3999 pA0 = fma(pA0, x4, CA2);
4000 pA1 = fma(pA1, x4, CA1);
4002 pA0 = fma(pA1, x2, pA0);
4003 // Constant term must come last for precision reasons
4010 msk_erf = (SimdDouble(0.75) <= y);
4011 t = maskzInvSingleAccuracy(y, msk_erf);
4016 expmx2 = expSingleAccuracy(-y * y);
4018 pB1 = fma(CB9, w2, CB7);
4019 pB0 = fma(CB8, w2, CB6);
4020 pB1 = fma(pB1, w2, CB5);
4021 pB0 = fma(pB0, w2, CB4);
4022 pB1 = fma(pB1, w2, CB3);
4023 pB0 = fma(pB0, w2, CB2);
4024 pB1 = fma(pB1, w2, CB1);
4025 pB0 = fma(pB0, w2, CB0);
4026 pB0 = fma(pB1, w, pB0);
4028 pC0 = fma(CC10, t2, CC8);
4029 pC1 = fma(CC9, t2, CC7);
4030 pC0 = fma(pC0, t2, CC6);
4031 pC1 = fma(pC1, t2, CC5);
4032 pC0 = fma(pC0, t2, CC4);
4033 pC1 = fma(pC1, t2, CC3);
4034 pC0 = fma(pC0, t2, CC2);
4035 pC1 = fma(pC1, t2, CC1);
4037 pC0 = fma(pC0, t2, CC0);
4038 pC0 = fma(pC1, t, pC0);
4041 // Select pB0 or pC0 for erfc()
4043 res_erfc = blend(pB0, pC0, mask);
4044 res_erfc = res_erfc * expmx2;
4046 // erfc(x<0) = 2-erfc(|x|)
4047 mask = (x < setZero());
4048 res_erfc = blend(res_erfc, two - res_erfc, mask);
4050 // Select erf() or erfc()
4051 mask = (y < SimdDouble(0.75));
4052 res = blend(one - res_erfc, res_erf, mask);
4057 /*! \brief SIMD erfc(x). Double precision SIMD data, single accuracy.
4059 * \param x The value to calculate erfc(x) for.
4062 * This routine achieves singleprecision (bar the last bit) over most of the
4063 * input range, but for large arguments where the result is getting close
4064 * to the minimum representable numbers we accept slightly larger errors
4065 * (think results that are in the ballpark of 10^-30) since that is not
4068 static inline SimdDouble gmx_simdcall erfcSingleAccuracy(SimdDouble x)
4070 // Coefficients for minimax approximation of erf(x)=x*P(x^2) in range [-1,1]
4071 const SimdDouble CA6(7.853861353153693e-5);
4072 const SimdDouble CA5(-8.010193625184903e-4);
4073 const SimdDouble CA4(5.188327685732524e-3);
4074 const SimdDouble CA3(-2.685381193529856e-2);
4075 const SimdDouble CA2(1.128358514861418e-1);
4076 const SimdDouble CA1(-3.761262582423300e-1);
4077 const SimdDouble CA0(1.128379165726710);
4078 // Coefficients for minimax approximation of erfc(x)=Exp(-x^2)*P((1/(x-1))^2) in range [0.67,2]
4079 const SimdDouble CB9(-0.0018629930017603923);
4080 const SimdDouble CB8(0.003909821287598495);
4081 const SimdDouble CB7(-0.0052094582210355615);
4082 const SimdDouble CB6(0.005685614362160572);
4083 const SimdDouble CB5(-0.0025367682853477272);
4084 const SimdDouble CB4(-0.010199799682318782);
4085 const SimdDouble CB3(0.04369575504816542);
4086 const SimdDouble CB2(-0.11884063474674492);
4087 const SimdDouble CB1(0.2732120154030589);
4088 const SimdDouble CB0(0.42758357702025784);
4089 // Coefficients for minimax approximation of erfc(x)=Exp(-x^2)*(1/x)*P((1/x)^2) in range [2,9.19]
4090 const SimdDouble CC10(-0.0445555913112064);
4091 const SimdDouble CC9(0.21376355144663348);
4092 const SimdDouble CC8(-0.3473187200259257);
4093 const SimdDouble CC7(0.016690861551248114);
4094 const SimdDouble CC6(0.7560973182491192);
4095 const SimdDouble CC5(-1.2137903600145787);
4096 const SimdDouble CC4(0.8411872321232948);
4097 const SimdDouble CC3(-0.08670413896296343);
4098 const SimdDouble CC2(-0.27124782687240334);
4099 const SimdDouble CC1(-0.0007502488047806069);
4100 const SimdDouble CC0(0.5642114853803148);
4101 const SimdDouble one(1.0);
4102 const SimdDouble two(2.0);
4104 SimdDouble x2, x4, y;
4105 SimdDouble t, t2, w, w2;
4106 SimdDouble pA0, pA1, pB0, pB1, pC0, pC1;
4108 SimdDouble res_erf, res_erfc, res;
4109 SimdDBool mask, msk_erf;
4115 pA0 = fma(CA6, x4, CA4);
4116 pA1 = fma(CA5, x4, CA3);
4117 pA0 = fma(pA0, x4, CA2);
4118 pA1 = fma(pA1, x4, CA1);
4120 pA0 = fma(pA0, x4, pA1);
4121 // Constant term must come last for precision reasons
4128 msk_erf = (SimdDouble(0.75) <= y);
4129 t = maskzInvSingleAccuracy(y, msk_erf);
4134 expmx2 = expSingleAccuracy(-y * y);
4136 pB1 = fma(CB9, w2, CB7);
4137 pB0 = fma(CB8, w2, CB6);
4138 pB1 = fma(pB1, w2, CB5);
4139 pB0 = fma(pB0, w2, CB4);
4140 pB1 = fma(pB1, w2, CB3);
4141 pB0 = fma(pB0, w2, CB2);
4142 pB1 = fma(pB1, w2, CB1);
4143 pB0 = fma(pB0, w2, CB0);
4144 pB0 = fma(pB1, w, pB0);
4146 pC0 = fma(CC10, t2, CC8);
4147 pC1 = fma(CC9, t2, CC7);
4148 pC0 = fma(pC0, t2, CC6);
4149 pC1 = fma(pC1, t2, CC5);
4150 pC0 = fma(pC0, t2, CC4);
4151 pC1 = fma(pC1, t2, CC3);
4152 pC0 = fma(pC0, t2, CC2);
4153 pC1 = fma(pC1, t2, CC1);
4155 pC0 = fma(pC0, t2, CC0);
4156 pC0 = fma(pC1, t, pC0);
4159 // Select pB0 or pC0 for erfc()
4161 res_erfc = blend(pB0, pC0, mask);
4162 res_erfc = res_erfc * expmx2;
4164 // erfc(x<0) = 2-erfc(|x|)
4165 mask = (x < setZero());
4166 res_erfc = blend(res_erfc, two - res_erfc, mask);
4168 // Select erf() or erfc()
4169 mask = (y < SimdDouble(0.75));
4170 res = blend(res_erfc, one - res_erf, mask);
4175 /*! \brief SIMD sin \& cos. Double precision SIMD data, single accuracy.
4177 * \param x The argument to evaluate sin/cos for
4178 * \param[out] sinval Sin(x)
4179 * \param[out] cosval Cos(x)
4181 static inline void gmx_simdcall sinCosSingleAccuracy(SimdDouble x, SimdDouble* sinval, SimdDouble* cosval)
4183 // Constants to subtract Pi/4*x from y while minimizing precision loss
4184 const SimdDouble argred0(2 * 0.78539816290140151978);
4185 const SimdDouble argred1(2 * 4.9604678871439933374e-10);
4186 const SimdDouble argred2(2 * 1.1258708853173288931e-18);
4187 const SimdDouble two_over_pi(2.0 / M_PI);
4188 const SimdDouble const_sin2(-1.9515295891e-4);
4189 const SimdDouble const_sin1(8.3321608736e-3);
4190 const SimdDouble const_sin0(-1.6666654611e-1);
4191 const SimdDouble const_cos2(2.443315711809948e-5);
4192 const SimdDouble const_cos1(-1.388731625493765e-3);
4193 const SimdDouble const_cos0(4.166664568298827e-2);
4195 const SimdDouble half(0.5);
4196 const SimdDouble one(1.0);
4197 SimdDouble ssign, csign;
4198 SimdDouble x2, y, z, psin, pcos, sss, ccc;
4201 # if GMX_SIMD_HAVE_DINT32_ARITHMETICS && GMX_SIMD_HAVE_LOGICAL
4202 const SimdDInt32 ione(1);
4203 const SimdDInt32 itwo(2);
4206 z = x * two_over_pi;
4210 mask = cvtIB2B((iy & ione) == setZero());
4211 ssign = selectByMask(SimdDouble(GMX_DOUBLE_NEGZERO), cvtIB2B((iy & itwo) == itwo));
4212 csign = selectByMask(SimdDouble(GMX_DOUBLE_NEGZERO), cvtIB2B(((iy + ione) & itwo) == itwo));
4214 const SimdDouble quarter(0.25);
4215 const SimdDouble minusquarter(-0.25);
4217 SimdDBool m1, m2, m3;
4219 /* The most obvious way to find the arguments quadrant in the unit circle
4220 * to calculate the sign is to use integer arithmetic, but that is not
4221 * present in all SIMD implementations. As an alternative, we have devised a
4222 * pure floating-point algorithm that uses truncation for argument reduction
4223 * so that we get a new value 0<=q<1 over the unit circle, and then
4224 * do floating-point comparisons with fractions. This is likely to be
4225 * slightly slower (~10%) due to the longer latencies of floating-point, so
4226 * we only use it when integer SIMD arithmetic is not present.
4230 // It is critical that half-way cases are rounded down
4231 z = fma(x, two_over_pi, half);
4235 /* z now starts at 0.0 for x=-pi/4 (although neg. values cannot occur), and
4236 * then increased by 1.0 as x increases by 2*Pi, when it resets to 0.0.
4237 * This removes the 2*Pi periodicity without using any integer arithmetic.
4238 * First check if y had the value 2 or 3, set csign if true.
4241 /* If we have logical operations we can work directly on the signbit, which
4242 * saves instructions. Otherwise we need to represent signs as +1.0/-1.0.
4243 * Thus, if you are altering defines to debug alternative code paths, the
4244 * two GMX_SIMD_HAVE_LOGICAL sections in this routine must either both be
4245 * active or inactive - you will get errors if only one is used.
4247 # if GMX_SIMD_HAVE_LOGICAL
4248 ssign = ssign & SimdDouble(GMX_DOUBLE_NEGZERO);
4249 csign = andNot(q, SimdDouble(GMX_DOUBLE_NEGZERO));
4250 ssign = ssign ^ csign;
4252 ssign = copysign(SimdDouble(1.0), ssign);
4253 csign = copysign(SimdDouble(1.0), q);
4255 ssign = ssign * csign; // swap ssign if csign was set.
4257 // Check if y had value 1 or 3 (remember we subtracted 0.5 from q)
4258 m1 = (q < minusquarter);
4259 m2 = (setZero() <= q);
4263 // where mask is FALSE, swap sign.
4264 csign = csign * blend(SimdDouble(-1.0), one, mask);
4266 x = fnma(y, argred0, x);
4267 x = fnma(y, argred1, x);
4268 x = fnma(y, argred2, x);
4271 psin = fma(const_sin2, x2, const_sin1);
4272 psin = fma(psin, x2, const_sin0);
4273 psin = fma(psin, x * x2, x);
4274 pcos = fma(const_cos2, x2, const_cos1);
4275 pcos = fma(pcos, x2, const_cos0);
4276 pcos = fms(pcos, x2, half);
4277 pcos = fma(pcos, x2, one);
4279 sss = blend(pcos, psin, mask);
4280 ccc = blend(psin, pcos, mask);
4281 // See comment for GMX_SIMD_HAVE_LOGICAL section above.
4282 # if GMX_SIMD_HAVE_LOGICAL
4283 *sinval = sss ^ ssign;
4284 *cosval = ccc ^ csign;
4286 *sinval = sss * ssign;
4287 *cosval = ccc * csign;
4291 /*! \brief SIMD sin(x). Double precision SIMD data, single accuracy.
4293 * \param x The argument to evaluate sin for
4296 * \attention Do NOT call both sin & cos if you need both results, since each of them
4297 * will then call \ref sincos and waste a factor 2 in performance.
4299 static inline SimdDouble gmx_simdcall sinSingleAccuracy(SimdDouble x)
4302 sinCosSingleAccuracy(x, &s, &c);
4306 /*! \brief SIMD cos(x). Double precision SIMD data, single accuracy.
4308 * \param x The argument to evaluate cos for
4311 * \attention Do NOT call both sin & cos if you need both results, since each of them
4312 * will then call \ref sincos and waste a factor 2 in performance.
4314 static inline SimdDouble gmx_simdcall cosSingleAccuracy(SimdDouble x)
4317 sinCosSingleAccuracy(x, &s, &c);
4321 /*! \brief SIMD tan(x). Double precision SIMD data, single accuracy.
4323 * \param x The argument to evaluate tan for
4326 static inline SimdDouble gmx_simdcall tanSingleAccuracy(SimdDouble x)
4328 const SimdDouble argred0(2 * 0.78539816290140151978);
4329 const SimdDouble argred1(2 * 4.9604678871439933374e-10);
4330 const SimdDouble argred2(2 * 1.1258708853173288931e-18);
4331 const SimdDouble two_over_pi(2.0 / M_PI);
4332 const SimdDouble CT6(0.009498288995810566122993911);
4333 const SimdDouble CT5(0.002895755790837379295226923);
4334 const SimdDouble CT4(0.02460087336161924491836265);
4335 const SimdDouble CT3(0.05334912882656359828045988);
4336 const SimdDouble CT2(0.1333989091464957704418495);
4337 const SimdDouble CT1(0.3333307599244198227797507);
4339 SimdDouble x2, p, y, z;
4342 # if GMX_SIMD_HAVE_DINT32_ARITHMETICS && GMX_SIMD_HAVE_LOGICAL
4346 z = x * two_over_pi;
4349 mask = cvtIB2B((iy & ione) == ione);
4351 x = fnma(y, argred0, x);
4352 x = fnma(y, argred1, x);
4353 x = fnma(y, argred2, x);
4354 x = selectByMask(SimdDouble(GMX_DOUBLE_NEGZERO), mask) ^ x;
4356 const SimdDouble quarter(0.25);
4357 const SimdDouble half(0.5);
4358 const SimdDouble threequarter(0.75);
4360 SimdDBool m1, m2, m3;
4363 z = fma(w, two_over_pi, half);
4367 m1 = (quarter <= q);
4369 m3 = (threequarter <= q);
4372 w = fnma(y, argred0, w);
4373 w = fnma(y, argred1, w);
4374 w = fnma(y, argred2, w);
4376 w = blend(w, -w, mask);
4377 x = w * copysign(SimdDouble(1.0), x);
4380 p = fma(CT6, x2, CT5);
4381 p = fma(p, x2, CT4);
4382 p = fma(p, x2, CT3);
4383 p = fma(p, x2, CT2);
4384 p = fma(p, x2, CT1);
4385 p = fma(x2, p * x, x);
4387 p = blend(p, maskzInvSingleAccuracy(p, mask), mask);
4391 /*! \brief SIMD asin(x). Double precision SIMD data, single accuracy.
4393 * \param x The argument to evaluate asin for
4396 static inline SimdDouble gmx_simdcall asinSingleAccuracy(SimdDouble x)
4398 const SimdDouble limitlow(1e-4);
4399 const SimdDouble half(0.5);
4400 const SimdDouble one(1.0);
4401 const SimdDouble halfpi(M_PI / 2.0);
4402 const SimdDouble CC5(4.2163199048E-2);
4403 const SimdDouble CC4(2.4181311049E-2);
4404 const SimdDouble CC3(4.5470025998E-2);
4405 const SimdDouble CC2(7.4953002686E-2);
4406 const SimdDouble CC1(1.6666752422E-1);
4408 SimdDouble z, z1, z2, q, q1, q2;
4410 SimdDBool mask, mask2;
4413 mask = (half < xabs);
4414 z1 = half * (one - xabs);
4415 mask2 = (xabs < one);
4416 q1 = z1 * maskzInvsqrtSingleAccuracy(z1, mask2);
4419 z = blend(z2, z1, mask);
4420 q = blend(q2, q1, mask);
4423 pA = fma(CC5, z2, CC3);
4424 pB = fma(CC4, z2, CC2);
4425 pA = fma(pA, z2, CC1);
4427 z = fma(pB, z2, pA);
4431 z = blend(z, q2, mask);
4433 mask = (limitlow < xabs);
4434 z = blend(xabs, z, mask);
4440 /*! \brief SIMD acos(x). Double precision SIMD data, single accuracy.
4442 * \param x The argument to evaluate acos for
4445 static inline SimdDouble gmx_simdcall acosSingleAccuracy(SimdDouble x)
4447 const SimdDouble one(1.0);
4448 const SimdDouble half(0.5);
4449 const SimdDouble pi(M_PI);
4450 const SimdDouble halfpi(M_PI / 2.0);
4452 SimdDouble z, z1, z2, z3;
4453 SimdDBool mask1, mask2, mask3;
4456 mask1 = (half < xabs);
4457 mask2 = (setZero() < x);
4459 z = half * (one - xabs);
4460 mask3 = (xabs < one);
4461 z = z * maskzInvsqrtSingleAccuracy(z, mask3);
4462 z = blend(x, z, mask1);
4463 z = asinSingleAccuracy(z);
4468 z = blend(z1, z2, mask2);
4469 z = blend(z3, z, mask1);
4474 /*! \brief SIMD asin(x). Double precision SIMD data, single accuracy.
4476 * \param x The argument to evaluate atan for
4477 * \result Atan(x), same argument/value range as standard math library.
4479 static inline SimdDouble gmx_simdcall atanSingleAccuracy(SimdDouble x)
4481 const SimdDouble halfpi(M_PI / 2);
4482 const SimdDouble CA17(0.002823638962581753730774);
4483 const SimdDouble CA15(-0.01595690287649631500244);
4484 const SimdDouble CA13(0.04250498861074447631836);
4485 const SimdDouble CA11(-0.07489009201526641845703);
4486 const SimdDouble CA9(0.1063479334115982055664);
4487 const SimdDouble CA7(-0.1420273631811141967773);
4488 const SimdDouble CA5(0.1999269574880599975585);
4489 const SimdDouble CA3(-0.3333310186862945556640);
4490 SimdDouble x2, x3, x4, pA, pB;
4491 SimdDBool mask, mask2;
4493 mask = (x < setZero());
4495 mask2 = (SimdDouble(1.0) < x);
4496 x = blend(x, maskzInvSingleAccuracy(x, mask2), mask2);
4501 pA = fma(CA17, x4, CA13);
4502 pB = fma(CA15, x4, CA11);
4503 pA = fma(pA, x4, CA9);
4504 pB = fma(pB, x4, CA7);
4505 pA = fma(pA, x4, CA5);
4506 pB = fma(pB, x4, CA3);
4507 pA = fma(pA, x2, pB);
4508 pA = fma(pA, x3, x);
4510 pA = blend(pA, halfpi - pA, mask2);
4511 pA = blend(pA, -pA, mask);
4516 /*! \brief SIMD atan2(y,x). Double precision SIMD data, single accuracy.
4518 * \param y Y component of vector, any quartile
4519 * \param x X component of vector, any quartile
4520 * \result Atan(y,x), same argument/value range as standard math library.
4522 * \note This routine should provide correct results for all finite
4523 * non-zero or positive-zero arguments. However, negative zero arguments will
4524 * be treated as positive zero, which means the return value will deviate from
4525 * the standard math library atan2(y,x) for those cases. That should not be
4526 * of any concern in Gromacs, and in particular it will not affect calculations
4527 * of angles from vectors.
4529 static inline SimdDouble gmx_simdcall atan2SingleAccuracy(SimdDouble y, SimdDouble x)
4531 const SimdDouble pi(M_PI);
4532 const SimdDouble halfpi(M_PI / 2.0);
4533 SimdDouble xinv, p, aoffset;
4534 SimdDBool mask_xnz, mask_ynz, mask_xlt0, mask_ylt0;
4536 mask_xnz = x != setZero();
4537 mask_ynz = y != setZero();
4538 mask_xlt0 = (x < setZero());
4539 mask_ylt0 = (y < setZero());
4541 aoffset = selectByNotMask(halfpi, mask_xnz);
4542 aoffset = selectByMask(aoffset, mask_ynz);
4544 aoffset = blend(aoffset, pi, mask_xlt0);
4545 aoffset = blend(aoffset, -aoffset, mask_ylt0);
4547 xinv = maskzInvSingleAccuracy(x, mask_xnz);
4549 p = atanSingleAccuracy(p);
4555 /*! \brief Analytical PME force correction, double SIMD data, single accuracy.
4557 * \param z2 \f$(r \beta)^2\f$ - see below for details.
4558 * \result Correction factor to coulomb force - see below for details.
4560 * This routine is meant to enable analytical evaluation of the
4561 * direct-space PME electrostatic force to avoid tables.
4563 * The direct-space potential should be \f$ \mbox{erfc}(\beta r)/r\f$, but there
4564 * are some problems evaluating that:
4566 * First, the error function is difficult (read: expensive) to
4567 * approxmiate accurately for intermediate to large arguments, and
4568 * this happens already in ranges of \f$(\beta r)\f$ that occur in simulations.
4569 * Second, we now try to avoid calculating potentials in Gromacs but
4570 * use forces directly.
4572 * We can simply things slight by noting that the PME part is really
4573 * a correction to the normal Coulomb force since \f$\mbox{erfc}(z)=1-\mbox{erf}(z)\f$, i.e.
4575 * V = \frac{1}{r} - \frac{\mbox{erf}(\beta r)}{r}
4577 * The first term we already have from the inverse square root, so
4578 * that we can leave out of this routine.
4580 * For pme tolerances of 1e-3 to 1e-8 and cutoffs of 0.5nm to 1.8nm,
4581 * the argument \f$beta r\f$ will be in the range 0.15 to ~4. Use your
4582 * favorite plotting program to realize how well-behaved \f$\frac{\mbox{erf}(z)}{z}\f$ is
4585 * We approximate \f$f(z)=\mbox{erf}(z)/z\f$ with a rational minimax polynomial.
4586 * However, it turns out it is more efficient to approximate \f$f(z)/z\f$ and
4587 * then only use even powers. This is another minor optimization, since
4588 * we actually \a want \f$f(z)/z\f$, because it is going to be multiplied by
4589 * the vector between the two atoms to get the vectorial force. The
4590 * fastest flops are the ones we can avoid calculating!
4592 * So, here's how it should be used:
4594 * 1. Calculate \f$r^2\f$.
4595 * 2. Multiply by \f$\beta^2\f$, so you get \f$z^2=(\beta r)^2\f$.
4596 * 3. Evaluate this routine with \f$z^2\f$ as the argument.
4597 * 4. The return value is the expression:
4600 * \frac{2 \exp{-z^2}}{\sqrt{\pi} z^2}-\frac{\mbox{erf}(z)}{z^3}
4603 * 5. Multiply the entire expression by \f$\beta^3\f$. This will get you
4606 * \frac{2 \beta^3 \exp(-z^2)}{\sqrt{\pi} z^2} - \frac{\beta^3 \mbox{erf}(z)}{z^3}
4609 * or, switching back to \f$r\f$ (since \f$z=r \beta\f$):
4612 * \frac{2 \beta \exp(-r^2 \beta^2)}{\sqrt{\pi} r^2} - \frac{\mbox{erf}(r \beta)}{r^3}
4615 * With a bit of math exercise you should be able to confirm that
4619 * \frac{\frac{d}{dr}\left( \frac{\mbox{erf}(\beta r)}{r} \right)}{r}
4622 * 6. Add the result to \f$r^{-3}\f$, multiply by the product of the charges,
4623 * and you have your force (divided by \f$r\f$). A final multiplication
4624 * with the vector connecting the two particles and you have your
4625 * vectorial force to add to the particles.
4627 * This approximation achieves an accuracy slightly lower than 1e-6; when
4628 * added to \f$1/r\f$ the error will be insignificant.
4631 static inline SimdDouble gmx_simdcall pmeForceCorrectionSingleAccuracy(SimdDouble z2)
4633 const SimdDouble FN6(-1.7357322914161492954e-8);
4634 const SimdDouble FN5(1.4703624142580877519e-6);
4635 const SimdDouble FN4(-0.000053401640219807709149);
4636 const SimdDouble FN3(0.0010054721316683106153);
4637 const SimdDouble FN2(-0.019278317264888380590);
4638 const SimdDouble FN1(0.069670166153766424023);
4639 const SimdDouble FN0(-0.75225204789749321333);
4641 const SimdDouble FD4(0.0011193462567257629232);
4642 const SimdDouble FD3(0.014866955030185295499);
4643 const SimdDouble FD2(0.11583842382862377919);
4644 const SimdDouble FD1(0.50736591960530292870);
4645 const SimdDouble FD0(1.0);
4648 SimdDouble polyFN0, polyFN1, polyFD0, polyFD1;
4652 polyFD0 = fma(FD4, z4, FD2);
4653 polyFD1 = fma(FD3, z4, FD1);
4654 polyFD0 = fma(polyFD0, z4, FD0);
4655 polyFD0 = fma(polyFD1, z2, polyFD0);
4657 polyFD0 = invSingleAccuracy(polyFD0);
4659 polyFN0 = fma(FN6, z4, FN4);
4660 polyFN1 = fma(FN5, z4, FN3);
4661 polyFN0 = fma(polyFN0, z4, FN2);
4662 polyFN1 = fma(polyFN1, z4, FN1);
4663 polyFN0 = fma(polyFN0, z4, FN0);
4664 polyFN0 = fma(polyFN1, z2, polyFN0);
4666 return polyFN0 * polyFD0;
4670 /*! \brief Analytical PME potential correction, double SIMD data, single accuracy.
4672 * \param z2 \f$(r \beta)^2\f$ - see below for details.
4673 * \result Correction factor to coulomb potential - see below for details.
4675 * This routine calculates \f$\mbox{erf}(z)/z\f$, although you should provide \f$z^2\f$
4676 * as the input argument.
4678 * Here's how it should be used:
4680 * 1. Calculate \f$r^2\f$.
4681 * 2. Multiply by \f$\beta^2\f$, so you get \f$z^2=\beta^2*r^2\f$.
4682 * 3. Evaluate this routine with z^2 as the argument.
4683 * 4. The return value is the expression:
4686 * \frac{\mbox{erf}(z)}{z}
4689 * 5. Multiply the entire expression by beta and switching back to \f$r\f$ (since \f$z=r \beta\f$):
4692 * \frac{\mbox{erf}(r \beta)}{r}
4695 * 6. Subtract the result from \f$1/r\f$, multiply by the product of the charges,
4696 * and you have your potential.
4698 * This approximation achieves an accuracy slightly lower than 1e-6; when
4699 * added to \f$1/r\f$ the error will be insignificant.
4701 static inline SimdDouble gmx_simdcall pmePotentialCorrectionSingleAccuracy(SimdDouble z2)
4703 const SimdDouble VN6(1.9296833005951166339e-8);
4704 const SimdDouble VN5(-1.4213390571557850962e-6);
4705 const SimdDouble VN4(0.000041603292906656984871);
4706 const SimdDouble VN3(-0.00013134036773265025626);
4707 const SimdDouble VN2(0.038657983986041781264);
4708 const SimdDouble VN1(0.11285044772717598220);
4709 const SimdDouble VN0(1.1283802385263030286);
4711 const SimdDouble VD3(0.0066752224023576045451);
4712 const SimdDouble VD2(0.078647795836373922256);
4713 const SimdDouble VD1(0.43336185284710920150);
4714 const SimdDouble VD0(1.0);
4717 SimdDouble polyVN0, polyVN1, polyVD0, polyVD1;
4721 polyVD1 = fma(VD3, z4, VD1);
4722 polyVD0 = fma(VD2, z4, VD0);
4723 polyVD0 = fma(polyVD1, z2, polyVD0);
4725 polyVD0 = invSingleAccuracy(polyVD0);
4727 polyVN0 = fma(VN6, z4, VN4);
4728 polyVN1 = fma(VN5, z4, VN3);
4729 polyVN0 = fma(polyVN0, z4, VN2);
4730 polyVN1 = fma(polyVN1, z4, VN1);
4731 polyVN0 = fma(polyVN0, z4, VN0);
4732 polyVN0 = fma(polyVN1, z2, polyVN0);
4734 return polyVN0 * polyVD0;
4740 /*! \name SIMD4 math functions
4742 * \note Only a subset of the math functions are implemented for SIMD4.
4747 # if GMX_SIMD4_HAVE_FLOAT
4749 /*************************************************************************
4750 * SINGLE PRECISION SIMD4 MATH FUNCTIONS - JUST A SMALL SUBSET SUPPORTED *
4751 *************************************************************************/
4753 /*! \brief Perform one Newton-Raphson iteration to improve 1/sqrt(x) for SIMD4 float.
4755 * This is a low-level routine that should only be used by SIMD math routine
4756 * that evaluates the inverse square root.
4758 * \param lu Approximation of 1/sqrt(x), typically obtained from lookup.
4759 * \param x The reference (starting) value x for which we want 1/sqrt(x).
4760 * \return An improved approximation with roughly twice as many bits of accuracy.
4762 static inline Simd4Float gmx_simdcall rsqrtIter(Simd4Float lu, Simd4Float x)
4764 Simd4Float tmp1 = x * lu;
4765 Simd4Float tmp2 = Simd4Float(-0.5F) * lu;
4766 tmp1 = fma(tmp1, lu, Simd4Float(-3.0F));
4770 /*! \brief Calculate 1/sqrt(x) for SIMD4 float.
4772 * \param x Argument that must be larger than GMX_FLOAT_MIN and smaller than
4773 * GMX_FLOAT_MAX, i.e. within the range of single precision.
4774 * For the single precision implementation this is obviously always
4775 * true for positive values, but for double precision it adds an
4776 * extra restriction since the first lookup step might have to be
4777 * performed in single precision on some architectures. Note that the
4778 * responsibility for checking falls on you - this routine does not
4780 * \return 1/sqrt(x). Result is undefined if your argument was invalid.
4782 static inline Simd4Float gmx_simdcall invsqrt(Simd4Float x)
4784 Simd4Float lu = rsqrt(x);
4785 # if (GMX_SIMD_RSQRT_BITS < GMX_SIMD_ACCURACY_BITS_SINGLE)
4786 lu = rsqrtIter(lu, x);
4788 # if (GMX_SIMD_RSQRT_BITS * 2 < GMX_SIMD_ACCURACY_BITS_SINGLE)
4789 lu = rsqrtIter(lu, x);
4791 # if (GMX_SIMD_RSQRT_BITS * 4 < GMX_SIMD_ACCURACY_BITS_SINGLE)
4792 lu = rsqrtIter(lu, x);
4798 # endif // GMX_SIMD4_HAVE_FLOAT
4801 # if GMX_SIMD4_HAVE_DOUBLE
4802 /*************************************************************************
4803 * DOUBLE PRECISION SIMD4 MATH FUNCTIONS - JUST A SMALL SUBSET SUPPORTED *
4804 *************************************************************************/
4806 /*! \brief Perform one Newton-Raphson iteration to improve 1/sqrt(x) for SIMD4 double.
4808 * This is a low-level routine that should only be used by SIMD math routine
4809 * that evaluates the inverse square root.
4811 * \param lu Approximation of 1/sqrt(x), typically obtained from lookup.
4812 * \param x The reference (starting) value x for which we want 1/sqrt(x).
4813 * \return An improved approximation with roughly twice as many bits of accuracy.
4815 static inline Simd4Double gmx_simdcall rsqrtIter(Simd4Double lu, Simd4Double x)
4817 Simd4Double tmp1 = x * lu;
4818 Simd4Double tmp2 = Simd4Double(-0.5F) * lu;
4819 tmp1 = fma(tmp1, lu, Simd4Double(-3.0F));
4823 /*! \brief Calculate 1/sqrt(x) for SIMD4 double.
4825 * \param x Argument that must be larger than GMX_FLOAT_MIN and smaller than
4826 * GMX_FLOAT_MAX, i.e. within the range of single precision.
4827 * For the single precision implementation this is obviously always
4828 * true for positive values, but for double precision it adds an
4829 * extra restriction since the first lookup step might have to be
4830 * performed in single precision on some architectures. Note that the
4831 * responsibility for checking falls on you - this routine does not
4833 * \return 1/sqrt(x). Result is undefined if your argument was invalid.
4835 static inline Simd4Double gmx_simdcall invsqrt(Simd4Double x)
4837 Simd4Double lu = rsqrt(x);
4838 # if (GMX_SIMD_RSQRT_BITS < GMX_SIMD_ACCURACY_BITS_DOUBLE)
4839 lu = rsqrtIter(lu, x);
4841 # if (GMX_SIMD_RSQRT_BITS * 2 < GMX_SIMD_ACCURACY_BITS_DOUBLE)
4842 lu = rsqrtIter(lu, x);
4844 # if (GMX_SIMD_RSQRT_BITS * 4 < GMX_SIMD_ACCURACY_BITS_DOUBLE)
4845 lu = rsqrtIter(lu, x);
4847 # if (GMX_SIMD_RSQRT_BITS * 8 < GMX_SIMD_ACCURACY_BITS_DOUBLE)
4848 lu = rsqrtIter(lu, x);
4854 /**********************************************************************
4855 * SIMD4 MATH FUNCTIONS WITH DOUBLE PREC. DATA, SINGLE PREC. ACCURACY *
4856 **********************************************************************/
4858 /*! \brief Calculate 1/sqrt(x) for SIMD4 double, but in single accuracy.
4860 * \param x Argument that must be larger than GMX_FLOAT_MIN and smaller than
4861 * GMX_FLOAT_MAX, i.e. within the range of single precision.
4862 * For the single precision implementation this is obviously always
4863 * true for positive values, but for double precision it adds an
4864 * extra restriction since the first lookup step might have to be
4865 * performed in single precision on some architectures. Note that the
4866 * responsibility for checking falls on you - this routine does not
4868 * \return 1/sqrt(x). Result is undefined if your argument was invalid.
4870 static inline Simd4Double gmx_simdcall invsqrtSingleAccuracy(Simd4Double x)
4872 Simd4Double lu = rsqrt(x);
4873 # if (GMX_SIMD_RSQRT_BITS < GMX_SIMD_ACCURACY_BITS_SINGLE)
4874 lu = rsqrtIter(lu, x);
4876 # if (GMX_SIMD_RSQRT_BITS * 2 < GMX_SIMD_ACCURACY_BITS_SINGLE)
4877 lu = rsqrtIter(lu, x);
4879 # if (GMX_SIMD_RSQRT_BITS * 4 < GMX_SIMD_ACCURACY_BITS_SINGLE)
4880 lu = rsqrtIter(lu, x);
4886 # endif // GMX_SIMD4_HAVE_DOUBLE
4890 # if GMX_SIMD_HAVE_FLOAT
4891 /*! \brief Calculate 1/sqrt(x) for SIMD float, only targeting single accuracy.
4893 * \param x Argument that must be larger than GMX_FLOAT_MIN and smaller than
4894 * GMX_FLOAT_MAX, i.e. within the range of single precision.
4895 * For the single precision implementation this is obviously always
4896 * true for positive values, but for double precision it adds an
4897 * extra restriction since the first lookup step might have to be
4898 * performed in single precision on some architectures. Note that the
4899 * responsibility for checking falls on you - this routine does not
4901 * \return 1/sqrt(x). Result is undefined if your argument was invalid.
4903 static inline SimdFloat gmx_simdcall invsqrtSingleAccuracy(SimdFloat x)
4908 /*! \brief Calculate 1/sqrt(x) for masked SIMD floats, only targeting single accuracy.
4910 * This routine only evaluates 1/sqrt(x) for elements for which mask is true.
4911 * Illegal values in the masked-out elements will not lead to
4912 * floating-point exceptions.
4914 * \param x Argument that must be larger than GMX_FLOAT_MIN and smaller than
4915 * GMX_FLOAT_MAX, i.e. within the range of single precision.
4916 * For the single precision implementation this is obviously always
4917 * true for positive values, but for double precision it adds an
4918 * extra restriction since the first lookup step might have to be
4919 * performed in single precision on some architectures. Note that the
4920 * responsibility for checking falls on you - this routine does not
4923 * \return 1/sqrt(x). Result is undefined if your argument was invalid or
4924 * entry was not masked, and 0.0 for masked-out entries.
4926 static inline SimdFloat maskzInvsqrtSingleAccuracy(SimdFloat x, SimdFBool m)
4928 return maskzInvsqrt(x, m);
4931 /*! \brief Calculate 1/sqrt(x) for two SIMD floats, only targeting single accuracy.
4933 * \param x0 First set of arguments, x0 must be in single range (see below).
4934 * \param x1 Second set of arguments, x1 must be in single range (see below).
4935 * \param[out] out0 Result 1/sqrt(x0)
4936 * \param[out] out1 Result 1/sqrt(x1)
4938 * In particular for double precision we can sometimes calculate square root
4939 * pairs slightly faster by using single precision until the very last step.
4941 * \note Both arguments must be larger than GMX_FLOAT_MIN and smaller than
4942 * GMX_FLOAT_MAX, i.e. within the range of single precision.
4943 * For the single precision implementation this is obviously always
4944 * true for positive values, but for double precision it adds an
4945 * extra restriction since the first lookup step might have to be
4946 * performed in single precision on some architectures. Note that the
4947 * responsibility for checking falls on you - this routine does not
4950 static inline void gmx_simdcall invsqrtPairSingleAccuracy(SimdFloat x0, SimdFloat x1, SimdFloat* out0, SimdFloat* out1)
4952 return invsqrtPair(x0, x1, out0, out1);
4955 /*! \brief Calculate 1/x for SIMD float, only targeting single accuracy.
4957 * \param x Argument with magnitude larger than GMX_FLOAT_MIN and smaller than
4958 * GMX_FLOAT_MAX, i.e. within the range of single precision.
4959 * For the single precision implementation this is obviously always
4960 * true for positive values, but for double precision it adds an
4961 * extra restriction since the first lookup step might have to be
4962 * performed in single precision on some architectures. Note that the
4963 * responsibility for checking falls on you - this routine does not
4965 * \return 1/x. Result is undefined if your argument was invalid.
4967 static inline SimdFloat gmx_simdcall invSingleAccuracy(SimdFloat x)
4973 /*! \brief Calculate 1/x for masked SIMD floats, only targeting single accuracy.
4975 * \param x Argument with magnitude larger than GMX_FLOAT_MIN and smaller than
4976 * GMX_FLOAT_MAX, i.e. within the range of single precision.
4977 * For the single precision implementation this is obviously always
4978 * true for positive values, but for double precision it adds an
4979 * extra restriction since the first lookup step might have to be
4980 * performed in single precision on some architectures. Note that the
4981 * responsibility for checking falls on you - this routine does not
4984 * \return 1/x for elements where m is true, or 0.0 for masked-out entries.
4986 static inline SimdFloat maskzInvSingleAccuracy(SimdFloat x, SimdFBool m)
4988 return maskzInv(x, m);
4991 /*! \brief Calculate sqrt(x) for SIMD float, always targeting single accuracy.
4993 * \copydetails sqrt(SimdFloat)
4995 template<MathOptimization opt = MathOptimization::Safe>
4996 static inline SimdFloat gmx_simdcall sqrtSingleAccuracy(SimdFloat x)
4998 return sqrt<opt>(x);
5001 /*! \brief Calculate cbrt(x) for SIMD float, always targeting single accuracy.
5003 * \copydetails cbrt(SimdFloat)
5005 static inline SimdFloat gmx_simdcall cbrtSingleAccuracy(SimdFloat x)
5010 /*! \brief Calculate 1/cbrt(x) for SIMD float, always targeting single accuracy.
5012 * \copydetails cbrt(SimdFloat)
5014 static inline SimdFloat gmx_simdcall invcbrtSingleAccuracy(SimdFloat x)
5019 /*! \brief SIMD float log2(x), only targeting single accuracy. This is the base-2 logarithm.
5021 * \param x Argument, should be >0.
5022 * \result The base-2 logarithm of x. Undefined if argument is invalid.
5024 static inline SimdFloat gmx_simdcall log2SingleAccuracy(SimdFloat x)
5029 /*! \brief SIMD float log(x), only targeting single accuracy. This is the natural logarithm.
5031 * \param x Argument, should be >0.
5032 * \result The natural logarithm of x. Undefined if argument is invalid.
5034 static inline SimdFloat gmx_simdcall logSingleAccuracy(SimdFloat x)
5039 /*! \brief SIMD float 2^x, only targeting single accuracy.
5041 * \copydetails exp2(SimdFloat)
5043 template<MathOptimization opt = MathOptimization::Safe>
5044 static inline SimdFloat gmx_simdcall exp2SingleAccuracy(SimdFloat x)
5046 return exp2<opt>(x);
5049 /*! \brief SIMD float e^x, only targeting single accuracy.
5051 * \copydetails exp(SimdFloat)
5053 template<MathOptimization opt = MathOptimization::Safe>
5054 static inline SimdFloat gmx_simdcall expSingleAccuracy(SimdFloat x)
5059 /*! \brief SIMD pow(x,y), only targeting single accuracy.
5061 * \copydetails pow(SimdFloat,SimdFloat)
5063 template<MathOptimization opt = MathOptimization::Safe>
5064 static inline SimdFloat gmx_simdcall powSingleAccuracy(SimdFloat x, SimdFloat y)
5066 return pow<opt>(x, y);
5069 /*! \brief SIMD float erf(x), only targeting single accuracy.
5071 * \param x The value to calculate erf(x) for.
5074 * This routine achieves very close to single precision, but we do not care about
5075 * the last bit or the subnormal result range.
5077 static inline SimdFloat gmx_simdcall erfSingleAccuracy(SimdFloat x)
5082 /*! \brief SIMD float erfc(x), only targeting single accuracy.
5084 * \param x The value to calculate erfc(x) for.
5087 * This routine achieves singleprecision (bar the last bit) over most of the
5088 * input range, but for large arguments where the result is getting close
5089 * to the minimum representable numbers we accept slightly larger errors
5090 * (think results that are in the ballpark of 10^-30) since that is not
5093 static inline SimdFloat gmx_simdcall erfcSingleAccuracy(SimdFloat x)
5098 /*! \brief SIMD float sin \& cos, only targeting single accuracy.
5100 * \param x The argument to evaluate sin/cos for
5101 * \param[out] sinval Sin(x)
5102 * \param[out] cosval Cos(x)
5104 static inline void gmx_simdcall sinCosSingleAccuracy(SimdFloat x, SimdFloat* sinval, SimdFloat* cosval)
5106 sincos(x, sinval, cosval);
5109 /*! \brief SIMD float sin(x), only targeting single accuracy.
5111 * \param x The argument to evaluate sin for
5114 * \attention Do NOT call both sin & cos if you need both results, since each of them
5115 * will then call \ref sincos and waste a factor 2 in performance.
5117 static inline SimdFloat gmx_simdcall sinSingleAccuracy(SimdFloat x)
5122 /*! \brief SIMD float cos(x), only targeting single accuracy.
5124 * \param x The argument to evaluate cos for
5127 * \attention Do NOT call both sin & cos if you need both results, since each of them
5128 * will then call \ref sincos and waste a factor 2 in performance.
5130 static inline SimdFloat gmx_simdcall cosSingleAccuracy(SimdFloat x)
5135 /*! \brief SIMD float tan(x), only targeting single accuracy.
5137 * \param x The argument to evaluate tan for
5140 static inline SimdFloat gmx_simdcall tanSingleAccuracy(SimdFloat x)
5145 /*! \brief SIMD float asin(x), only targeting single accuracy.
5147 * \param x The argument to evaluate asin for
5150 static inline SimdFloat gmx_simdcall asinSingleAccuracy(SimdFloat x)
5155 /*! \brief SIMD float acos(x), only targeting single accuracy.
5157 * \param x The argument to evaluate acos for
5160 static inline SimdFloat gmx_simdcall acosSingleAccuracy(SimdFloat x)
5165 /*! \brief SIMD float atan(x), only targeting single accuracy.
5167 * \param x The argument to evaluate atan for
5168 * \result Atan(x), same argument/value range as standard math library.
5170 static inline SimdFloat gmx_simdcall atanSingleAccuracy(SimdFloat x)
5175 /*! \brief SIMD float atan2(y,x), only targeting single accuracy.
5177 * \param y Y component of vector, any quartile
5178 * \param x X component of vector, any quartile
5179 * \result Atan(y,x), same argument/value range as standard math library.
5181 * \note This routine should provide correct results for all finite
5182 * non-zero or positive-zero arguments. However, negative zero arguments will
5183 * be treated as positive zero, which means the return value will deviate from
5184 * the standard math library atan2(y,x) for those cases. That should not be
5185 * of any concern in Gromacs, and in particular it will not affect calculations
5186 * of angles from vectors.
5188 static inline SimdFloat gmx_simdcall atan2SingleAccuracy(SimdFloat y, SimdFloat x)
5193 /*! \brief SIMD Analytic PME force correction, only targeting single accuracy.
5195 * \param z2 \f$(r \beta)^2\f$ - see default single precision version for details.
5196 * \result Correction factor to coulomb force.
5198 static inline SimdFloat gmx_simdcall pmeForceCorrectionSingleAccuracy(SimdFloat z2)
5200 return pmeForceCorrection(z2);
5203 /*! \brief SIMD Analytic PME potential correction, only targeting single accuracy.
5205 * \param z2 \f$(r \beta)^2\f$ - see default single precision version for details.
5206 * \result Correction factor to coulomb force.
5208 static inline SimdFloat gmx_simdcall pmePotentialCorrectionSingleAccuracy(SimdFloat z2)
5210 return pmePotentialCorrection(z2);
5212 # endif // GMX_SIMD_HAVE_FLOAT
5214 # if GMX_SIMD4_HAVE_FLOAT
5215 /*! \brief Calculate 1/sqrt(x) for SIMD4 float, only targeting single accuracy.
5217 * \param x Argument that must be larger than GMX_FLOAT_MIN and smaller than
5218 * GMX_FLOAT_MAX, i.e. within the range of single precision.
5219 * For the single precision implementation this is obviously always
5220 * true for positive values, but for double precision it adds an
5221 * extra restriction since the first lookup step might have to be
5222 * performed in single precision on some architectures. Note that the
5223 * responsibility for checking falls on you - this routine does not
5225 * \return 1/sqrt(x). Result is undefined if your argument was invalid.
5227 static inline Simd4Float gmx_simdcall invsqrtSingleAccuracy(Simd4Float x)
5231 # endif // GMX_SIMD4_HAVE_FLOAT
5233 /*! \} end of addtogroup module_simd */
5234 /*! \endcond end of condition libabl */
5240 #endif // GMX_SIMD_SIMD_MATH_H