text: #pragma once
  
  /*
   * From
   * https://github.com/scipy/scipy/blob/5caee5d4ad564cfae4596f8dfa8b45997767035b/scipy/special/cephes/i1.c
   */
  
  /*                                                     i1.c
   *
   *     Modified Bessel function of order one
   *
   *
   *
   * SYNOPSIS:
   *
   * double x, y, i1();
   *
   * y = i1( x );
   *
   *
   *
   * DESCRIPTION:
   *
   * Returns modified Bessel function of order one of the
   * argument.
   *
   * The function is defined as i1(x) = -i j1( ix ).
   *
   * The range is partitioned into the two intervals [0,8] and
   * (8, infinity).  Chebyshev polynomial expansions are employed
   * in each interval.
   *
   *
   *
   * ACCURACY:
   *
   *                      Relative error:
   * arithmetic   domain     # trials      peak         rms
   *    IEEE      0, 30       30000       1.9e-15     2.1e-16
   *
   *
   */
  /*							i1e.c
   *
   *	Modified Bessel function of order one,
   *	exponentially scaled
   *
   *
   *
   * SYNOPSIS:
   *
   * double x, y, i1e();
   *
   * y = i1e( x );
   *
   *
   *
   * DESCRIPTION:
   *
   * Returns exponentially scaled modified Bessel function
   * of order one of the argument.
   *
   * The function is defined as i1(x) = -i exp(-|x|) j1( ix ).
   *
   *
   *
   * ACCURACY:
   *
   *                      Relative error:
   * arithmetic   domain     # trials      peak         rms
   *    IEEE      0, 30       30000       2.0e-15     2.0e-16
   * See i1().
   *
   */
  
  /*                                                     i1.c 2          */
  
  
  /*
   * Cephes Math Library Release 2.8:  June, 2000
   * Copyright 1985, 1987, 2000 by Stephen L. Moshier
   */
  
  #include "chbevl.h"
  #include <cmath>
  
  
  namespace cephes { namespace cpu {
  
  namespace {
  
  using std::abs;
  using std::exp;
  using std::sqrt;
  
  
  template<typename scalar_t, bool e>
  static inline scalar_t calc_i1(scalar_t x) {
      /* Chebyshev coefficients for exp(-x) I1(x) / x
      * in the interval [0,8].
      *
      * lim(x->0){ exp(-x) I1(x) / x } = 1/2.
      */
  
      static const scalar_t A[] = {
          2.77791411276104639959E-18,
          -2.11142121435816608115E-17,
          1.55363195773620046921E-16,
          -1.10559694773538630805E-15,
          7.60068429473540693410E-15,
          -5.04218550472791168711E-14,
          3.22379336594557470981E-13,
          -1.98397439776494371520E-12,
          1.17361862988909016308E-11,
          -6.66348972350202774223E-11,
          3.62559028155211703701E-10,
          -1.88724975172282928790E-9,
          9.38153738649577178388E-9,
          -4.44505912879632808065E-8,
          2.00329475355213526229E-7,
          -8.56872026469545474066E-7,
          3.47025130813767847674E-6,
          -1.32731636560394358279E-5,
          4.78156510755005422638E-5,
          -1.61760815825896745588E-4,
          5.12285956168575772895E-4,
          -1.51357245063125314899E-3,
          4.15642294431288815669E-3,
          -1.05640848946261981558E-2,
          2.47264490306265168283E-2,
          -5.29459812080949914269E-2,
          1.02643658689847095384E-1,
          -1.76416518357834055153E-1,
          2.52587186443633654823E-1
      };
  
      /* Chebyshev coefficients for exp(-x) sqrt(x) I1(x)
      * in the inverted interval [8,infinity].
      *
      * lim(x->inf){ exp(-x) sqrt(x) I1(x) } = 1/sqrt(2pi).
      */
      static const scalar_t B[] = {
          7.51729631084210481353E-18,
          4.41434832307170791151E-18,
          -4.65030536848935832153E-17,
          -3.20952592199342395980E-17,
          2.96262899764595013876E-16,
          3.30820231092092828324E-16,
          -1.88035477551078244854E-15,
          -3.81440307243700780478E-15,
          1.04202769841288027642E-14,
          4.27244001671195135429E-14,
          -2.10154184277266431302E-14,
          -4.08355111109219731823E-13,
          -7.19855177624590851209E-13,
          2.03562854414708950722E-12,
          1.41258074366137813316E-11,
          3.25260358301548823856E-11,
          -1.89749581235054123450E-11,
          -5.58974346219658380687E-10,
          -3.83538038596423702205E-9,
          -2.63146884688951950684E-8,
          -2.51223623787020892529E-7,
          -3.88256480887769039346E-6,
          -1.10588938762623716291E-4,
          -9.76109749136146840777E-3,
          7.78576235018280120474E-1
      };
  
  
      scalar_t y, z;
  
      z = abs(x);
      if (z <= 8.0) {
          y = (z / 2.0) - 2.0;
          if (e) {
              z = chbevl<scalar_t>(y, A, 29) * z;
          } else {
              z = chbevl<scalar_t>(y, A, 29) * z * exp(z);
          }
      }
      else {
          if (e) {
  	        z = chbevl<scalar_t>(32.0 / z - 2.0, B, 25) / sqrt(z);
          } else {
  	        z = exp(z) * chbevl<scalar_t>(32.0 / z - 2.0, B, 25) / sqrt(z);
          }
      }
      if (x < 0.0) {
  	    z = -z;
      }
      return (z);
  
  }
  
  }
  
  template<typename scalar_t>
  static inline scalar_t i1(scalar_t x) {
      return calc_i1<scalar_t, false>(x);
  }
  
  template<typename scalar_t>
  static inline scalar_t i1e(scalar_t x) {
      return calc_i1<scalar_t, true>(x);
  }
  
  } // namespace cpu
  } // namespace cephes
  