/*
 * Copyright (C) Bernard Gingold, 2020-2026 
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 3 of the License, or (at
 * your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
*/

#ifndef __GMS_BESSEL_I0_SSE_H__
#define __GMS_BESSEL_I0_SSE_H__ 050820260642

#include <cstdint>
#include <immintrin.h>
#include "GMS_config.h"


namespace file_info 
{

     static const unsigned int GMS_BESSEL_I0_SSE_MAJOR = 1;
     static const unsigned int GMS_BESSEL_I0_SSE_MINOR = 1;
     static const unsigned int GMS_BESSEL_I0_SSE_MICRO = 0;
     static const unsigned int GMS_BESSEL_I0_SSE_FULLVER =
       1000U*GMS_BESSEL_I0_SSE_MAJOR+100U*GMS_BESSEL_I0_SSE_MINOR+
       10U*GMS_BESSEL_I0_SSE_MICRO;
     static const char GMS_BESSEL_I0_SSE_CREATION_DATE[] = "05-08-2026 06:43AM +00200 (WED 05 AUG 2026 GMT+2)";
     static const char GMS_BESSEL_I0_SSE_BUILD_DATE[]    = __DATE__; 
     static const char GMS_BESSEL_I0_SSE_BUILD_TIME[]    = __TIME__;
     static const char GMS_BESSEL_I0_SSE_SYNOPSIS[]      = "Manual vectorization of the Bessel-I0 scalar function. Authors: \
                                                            Original FORTRAN77 version by William Cody, Laura Stoltz. FORTRAN90 version by John Burkardt. \
                                                            The SSE version created by Bernard Gingold";
                                                           
}

namespace gms
{

namespace math
{

/*
 !*****************************************************************************80
!
!! BESSEL_I0 evaluates the modified Bessel function I0(X).
!
!  Discussion:
!
!    The main computation evaluates slightly modified forms of
!    minimax approximations generated by Blair and Edwards, Chalk
!    River (Atomic Energy of Canada Limited) Report AECL-4928,
!    October, 1974.  This transportable program is patterned after
!    the machine dependent FUNPACK packet NATSI0, but cannot match
!    that version for efficiency or accuracy.  This version uses
!    rational functions that theoretically approximate I-SUB-0(X)
!    to at least 18 significant decimal digits.
!
!  Machine dependent constants:
!
!    beta   = Radix for the floating-point system
!    maxexp = Smallest power of beta that overflows
!    XMAX =   Largest argument acceptable to BESI0;  Solution to
!             equation:
!               W(X) * (1+1/(8*X)+9/(128*X^2) = beta^maxexp
!             where  W(X) = EXP(X)/sqrt(2*PI*X)
!
!    Approximate values for some important machines are:
!
!                             beta       maxexp       XMAX
!
!    CRAY-1        (S.P.)       2         8191       5682.810
!    Cyber 180/855
!      under NOS   (S.P.)       2         1070        745.893
!    IEEE (IBM/XT,
!      SUN, etc.)  (S.P.)       2          128         91.900
!    IEEE (IBM/XT,
!      SUN, etc.)  (D.P.)       2         1024        713.986
!    IBM 3033      (D.P.)      16           63        178.182
!    VAX           (S.P.)       2          127         91.203
!    VAX D-Format  (D.P.)       2          127         91.203
!    VAX G-Format  (D.P.)       2         1023        713.293
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license.
!
!  Modified:
!
!    17 October 2008
!
!  Author:
!
!    Original FORTRAN77 version by William Cody, Laura Stoltz.
!    FORTRAN90 version by John Burkardt.
!
!  Parameters:
!
!    Input, real ( kind = 8 ) ARG, the argument.
!
!    Output, real ( kind = 8 ) BESSEL_I0, the value of the modified
!    Bessel function of the first kind.
!
*/
/*
   The original implementation is splitted into two parts:
   1) ABS(ARG) <= 15.0
   2) 15.0 <= ABS(ARG)
*/
#if defined(__INTEL_COMPILER) || defined(__ICC)
#pragma intel optimization_level 3 
#pragma intel optimization_parameter target_arch=SSE
#elif defined (__GNUC__) && (!defined (__INTEL_COMPILER) || !defined(__ICC))
#pragma GCC optimize("O3")
#pragma GCC target("sse")
#endif
__ATTR_ALWAYS_INLINE__
static inline 
__m128d
bessel_i0_le15_sse_pd(const __m128d arg)
{
  __ATTR_ALIGN__(16) const __m128d p[15] = 
  {
    _mm_set1_pd(-5.2487866627945699800E-18),
    _mm_set1_pd(-1.5982226675653184646E-14), 
    _mm_set1_pd(-2.6843448573468483278E-11), 
    _mm_set1_pd(-3.0517226450451067446E-08), 
    _mm_set1_pd(-2.5172644670688975051E-05), 
    _mm_set1_pd(-1.5453977791786851041E-02), 
    _mm_set1_pd(-7.0935347449210549190E+00), 
    _mm_set1_pd(-2.4125195876041896775E+03), 
    _mm_set1_pd(-5.9545626019847898221E+05), 
    _mm_set1_pd(-1.0313066708737980747E+08), 
    _mm_set1_pd(-1.1912746104985237192E+10), 
    _mm_set1_pd(-8.4925101247114157499E+11), 
    _mm_set1_pd(-3.2940087627407749166E+13), 
    _mm_set1_pd(-5.5050369673018427753E+14), 
    _mm_set1_pd(-2.2335582639474375249E+15)
  };
  __ATTR_ALIGN__(16) const __m128d q[5] = 
  {
    _mm_set1_pd(-3.7277560179962773046E+03), 
    _mm_set1_pd(6.5158506418655165707E+06), 
    _mm_set1_pd(-6.5626560740833869295E+09), 
    _mm_set1_pd(3.7604188704092954661E+12), 
    _mm_set1_pd(-9.7087946179594019126E+14)
  };
  const __m128d C225{_mm_set1_pd(225.0)};
  __m128d sump;
  __m128d sumq;
  __m128d result;
  //__m128d x  = _mm_abs_pd(arg);
  __m128d xx = _mm_mul_pd(arg,arg);
  sump = p[0];
  sump = _mm_fmadd_pd(sump,xx,p[1]);
  sump = _mm_fmadd_pd(sump,xx,p[2]);
  sump = _mm_fmadd_pd(sump,xx,p[3]);
  sump = _mm_fmadd_pd(sump,xx,p[4]);
  sump = _mm_fmadd_pd(sump,xx,p[5]);
  sump = _mm_fmadd_pd(sump,xx,p[6]);
  sump = _mm_fmadd_pd(sump,xx,p[7]);
  sump = _mm_fmadd_pd(sump,xx,p[8]);
  sump = _mm_fmadd_pd(sump,xx,p[9]);
  sump = _mm_fmadd_pd(sump,xx,p[10]);
  sump = _mm_fmadd_pd(sump,xx,p[11]);
  sump = _mm_fmadd_pd(sump,xx,p[12]);
  sump = _mm_fmadd_pd(sump,xx,p[13]);
  sump = _mm_fmadd_pd(sump,xx,p[14]);
  xx   = _mm_sub_pd(xx,C225);
  /*
     double sumq = ((((xx + q[0]) * xx + q[1]) * xx + q[2]) * xx + q[3]) * xx + q[4];
  */	
  sumq   = _mm_fmadd_pd(_mm_fmadd_pd(_mm_fmadd_pd(_mm_fmadd_pd(_mm_add_pd(xx,q[0]),xx,q[1]),xx,q[2]),xx,q[3]),xx,q[4]);    
  result = _mm_div_pd(sump,sumq);
  return (result);
}

#if defined(__INTEL_COMPILER) || defined(__ICC)
#pragma intel optimization_level 3 
#pragma intel optimization_parameter target_arch=SSE
#elif defined (__GNUC__) && (!defined (__INTEL_COMPILER) || !defined(__ICC))
#pragma GCC optimize("O3")
#pragma GCC target("sse")
#endif
__ATTR_ALWAYS_INLINE__
static inline 
__m128d
bessel_i0_ge15_sse_pd(const __m128d arg)
{
  __ATTR_ALIGN__(16) const __m128d pp[8] = 
  {
    _mm_set1_pd(-3.9843750000000000000E-01), 
    _mm_set1_pd(2.9205384596336793945E+00), 
    _mm_set1_pd(-2.4708469169133954315E+00), 
    _mm_set1_pd(4.7914889422856814203E-01), 
    _mm_set1_pd(-3.7384991926068969150E-03), 
    _mm_set1_pd(-2.6801520353328635310E-03), 
    _mm_set1_pd(9.9168777670983678974E-05), 
    _mm_set1_pd(-2.1877128189032726730E-06)
  };
  
  __ATTR_ALIGN__(16) const __m128d qq[7] =
  {
    _mm_set1_pd(-3.1446690275135491500E+01), 
    _mm_set1_pd(8.5539563258012929600E+01), 
    _mm_set1_pd(-6.0228002066743340583E+01), 
    _mm_set1_pd(1.3982595353892851542E+01), 
    _mm_set1_pd(-1.1151759188741312645E+00), 
    _mm_set1_pd(3.2547697594819615062E-02), 
    _mm_set1_pd(-5.5194330231005480228E-04)
  };
  __m128d sump;
  __m128d sumq;
  __m128d a;
  __m128d b;
  __m128d result;
  const __m128d rec15 = _mm_set1_pd(6.6666666666666666666E-02);
  const __m128d one   = _mm_set1_pd(1.0e+00);
  const __m128d C15   = _mm_set1_pd(15.0e+00);
  const __m128d exp40 = _mm_set1_pd(2.353852668370199854E+17);
  const __m128d C40   = _mm_set1_pd(40.0e+00);
  const __m128d xmax  = _mm_set1_pd(91.9E+00);
  __m128d       xx    = _mm_sub_pd(_mm_div_pd(one,arg),rec15);
  const __m128d t0    = _mm_fmadd_pd(pp[0],xx,pp[1]);
  const __m128d c0    = _mm_fmadd_pd(_mm_add_pd(xx,qq[0]),xx,qq[1]);
  const __m128d t1    = _mm_fmadd_pd(t0,xx,pp[2]);
  const __m128d c1    = _mm_fmadd_pd(c0,xx,qq[2]);
  const __m128d t2    = _mm_fmadd_pd(t1,xx,pp[3]);
  const __m128d c2    = _mm_fmadd_pd(c1,xx,qq[3]);
  const __m128d t3    = _mm_fmadd_pd(t2,xx,pp[4]);
  const __m128d c3    = _mm_fmadd_pd(c2,xx,qq[4]);
  const __m128d t4    = _mm_fmadd_pd(t3,xx,pp[5]);
  const __m128d c4    = _mm_fmadd_pd(c3,xx,qq[5]);
  const __m128d t5    = _mm_fmadd_pd(t4,xx,pp[6]);
  const __m128d c5    = _mm_fmadd_pd(c4,xx,qq[6]);
  const __m128d t6    = _mm_fmadd_pd(t5,xx,pp[7]);
  sump                = t6;
  sumq                = c5;
  result              = _mm_div_pd(sump,sumq);
  const __mmask8 m    = _mm_cmp_pd_mask(arg,_mm_sub_pd(xmax,C15),_CMP_LE_OQ);
  a                   = _mm_mask_blend_pd(m,_mm_exp_pd(_mm_sub_pd(arg,C40)),
					                                             _mm_exp_pd(arg));
  b                   = _mm_mask_blend_pd(m,exp40,one);
  const __m128d tmp   = _mm_fmsub_pd(result,a,_mm_mul_pd(pp[0],a));
  result              = _mm_mul_pd(_mm_div_pd(tmp,_mm_sqrt_pd(arg)),b);
  return (result);
}


} // math

} // gms

#endif /*__GMS_BESSEL_I0_SSE_H__*/