kopia lustrzana https://github.com/Dsplib/libdspl-2.0
188 wiersze
6.6 KiB
C
188 wiersze
6.6 KiB
C
/*
|
||
* Copyright (c) 2015-2024 Sergey Bakhurin
|
||
* Digital Signal Processing Library [http://dsplib.org]
|
||
*
|
||
* This file is part of DSPL.
|
||
*
|
||
* 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.
|
||
*
|
||
* DSPL 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 Foobar. If not, see <http://www.gnu.org/licenses/>.
|
||
*/
|
||
|
||
|
||
#include <stdio.h>
|
||
#include <stdlib.h>
|
||
#include "dspl.h"
|
||
|
||
|
||
|
||
#ifdef DOXYGEN_ENGLISH
|
||
/*! ****************************************************************************
|
||
\ingroup SPEC_MATH_TRANSCEND
|
||
|
||
\brief Modified Bessel Function of the First Kind – \f$ I_0(x)\f$ [1].
|
||
|
||
\param[in] x
|
||
Pointer to the function argument vector \f$ x \f$. \n
|
||
Vector size is `[n x 1]`. \n
|
||
Input vector must contain nonnegative values. \n \n
|
||
|
||
\param[in] n
|
||
Input vector size `x`. \n \n
|
||
|
||
\param[out] y
|
||
Pointer to \f$ I_0(x)\f$ function vector. \n
|
||
Vector size is `[n x 1]`. \n
|
||
Memory must be allocated. \n \n
|
||
|
||
\return
|
||
`RES_OK` if function calculated successfully. \n
|
||
Else \ref ERROR_CODE_GROUP "code error". \n
|
||
|
||
\note
|
||
[1] Rational Approximations for the Modified Bessel Function
|
||
of the First Kind – I0(x) for Computations with Double Precision
|
||
by PAVEL HOLOBORODKO on NOVEMBER 11, 2015
|
||
|
||
Example:
|
||
|
||
\include bessel_i0.c
|
||
|
||
Program calcultes \f$ I_0(x)\f$ function for `x`
|
||
in \f$[0 \ 3]\f$ interval.
|
||
Data saved if `dat/dat0.txt` file and shows on the plot
|
||
|
||
\image html bessel_i0.png
|
||
|
||
\author Sergey Bakhurin www.dsplib.org
|
||
***************************************************************************** */
|
||
#endif
|
||
#ifdef DOXYGEN_RUSSIAN
|
||
/*! ****************************************************************************
|
||
\ingroup SPEC_MATH_TRANSCEND
|
||
|
||
\brief
|
||
Модифицированная функция Бесселя первого рода \f$ I_0(x)\f$.
|
||
|
||
Функция рассчитывает значения функции для вещественного вектора `x`,
|
||
который должен принимать неотрицательные значения. \n
|
||
|
||
\param[in] x
|
||
Указатель на вектор переменной \f$ x \f$. \n
|
||
Размер вектора `[n x 1]`. \n
|
||
Память должна быть выделена. \n \n
|
||
|
||
\param[in] n
|
||
Размер входного вектора `x`. \n \n
|
||
|
||
\param[out] y
|
||
Указатель на вектор значений функции \f$ I_0(x)\f$. \n
|
||
Размер вектора `[n x 1]`. \n
|
||
Память должна быть выделена. \n \n
|
||
|
||
\return
|
||
`RES_OK` --- расчёт произведен успешно. \n
|
||
В противном случае \ref ERROR_CODE_GROUP "код ошибки". \n
|
||
|
||
\note
|
||
Используемый алгоритм описан в статье:
|
||
Rational Approximations for the Modified Bessel Function
|
||
of the First Kind – I0(x) for Computations with Double Precision
|
||
by PAVEL HOLOBORODKO on NOVEMBER 11, 2015
|
||
|
||
Пример использования функции `bessel_i0`:
|
||
|
||
\include bessel_i0.c
|
||
|
||
Данная программа рассчитывает значения функции \f$ I_0(x)\f$ переменной `x`
|
||
в интервале \f$[0 \ 3]\f$.
|
||
Рассчитанные данные сохраняются в текстовый файл `dat/dat0.txt`
|
||
и выводятся на график `img/bessel_i0.png`
|
||
|
||
\image html bessel_i0.png
|
||
|
||
\author Бахурин Сергей www.dsplib.org
|
||
***************************************************************************** */
|
||
#endif
|
||
int DSPL_API bessel_i0(double* x, int n, double* y)
|
||
{
|
||
double P16[17] = { 1.0000000000000000000000801e+00,
|
||
2.4999999999999999999629693e-01,
|
||
2.7777777777777777805664954e-02,
|
||
1.7361111111111110294015271e-03,
|
||
6.9444444444444568581891535e-05,
|
||
1.9290123456788994104574754e-06,
|
||
3.9367598891475388547279760e-08,
|
||
6.1511873265092916275099070e-10,
|
||
7.5940584360755226536109511e-12,
|
||
7.5940582595094190098755663e-14,
|
||
6.2760839879536225394314453e-16,
|
||
4.3583591008893599099577755e-18,
|
||
2.5791926805873898803749321e-20,
|
||
1.3141332422663039834197910e-22,
|
||
5.9203280572170548134753422e-25,
|
||
2.0732014503197852176921968e-27,
|
||
1.1497640034400735733456400e-29};
|
||
|
||
double P22[23] = { 3.9894228040143265335649948e-01,
|
||
4.9867785050353992900698488e-02,
|
||
2.8050628884163787533196746e-02,
|
||
2.9219501690198775910219311e-02,
|
||
4.4718622769244715693031735e-02,
|
||
9.4085204199017869159183831e-02,
|
||
-1.0699095472110916094973951e-01,
|
||
2.2725199603010833194037016e+01,
|
||
-1.0026890180180668595066918e+03,
|
||
3.1275740782277570164423916e+04,
|
||
-5.9355022509673600842060002e+05,
|
||
2.6092888649549172879282592e+06,
|
||
2.3518420447411254516178388e+08,
|
||
-8.9270060370015930749184222e+09,
|
||
1.8592340458074104721496236e+11,
|
||
-2.6632742974569782078420204e+12,
|
||
2.7752144774934763122129261e+13,
|
||
-2.1323049786724612220362154e+14,
|
||
1.1989242681178569338129044e+15,
|
||
-4.8049082153027457378879746e+15,
|
||
1.3012646806421079076251950e+16,
|
||
-2.1363029690365351606041265e+16,
|
||
1.6069467093441596329340754e+16};
|
||
|
||
double x2;
|
||
int k;
|
||
|
||
if(!x || !y)
|
||
return ERROR_PTR;
|
||
if(n < 1)
|
||
return ERROR_SIZE;
|
||
|
||
for(k =0; k < n; k++)
|
||
{
|
||
if(x[k] < 0.0)
|
||
return ERROR_NEGATIVE;
|
||
|
||
if(x[k] < 7.75)
|
||
{
|
||
x2 = x[k] * x[k] * 0.25;
|
||
polyval(P16, 16, &x2, 1, y+k);
|
||
y[k] = x2 * y[k] + 1.0;
|
||
}
|
||
else
|
||
{
|
||
x2 = 1.0 / x[k];
|
||
polyval(P22, 22, &x2, 1, y+k);
|
||
y[k] *= exp(x[k]) / sqrt(x[k]);
|
||
}
|
||
}
|
||
return RES_OK;
|
||
}
|