bessel_i0.c

/*
 * Copyright (c) 2002, 2017 Jens Keiner, Stefan Kunis, Daniel Potts
 *
 * 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 2 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.
 */

#include "infft.h"

#if defined(NFFT_LDOUBLE)
  #if LDBL_MANT_DIG > 64
    /* long double 128 bit wide */
    static const R P1[] =
    {
        K(0.9999999999999999999999999999999999962027889818871465625924752326583154315551642695323965510694779151206544639605419513240736857),
        K(0.25000000000000000000000000000000012196714800029501487764875563265263401524033425368818847175854607800541392646833606522273821246),
        K(0.027777777777777777777777777777777144587988290342214362953013414854374623029362784740516243952857824585136440104808521394795213943),
        K(0.0017361111111111111111111111111123869249041395763984631620599418268304349163419724426016945727467811299651660803976432220176426406),
        K(0.000069444444444444444444444444443107647328522467237563132750155356213911633630315806134293192377451924142086027094615621024770286162),
        K(1.9290123456790123456790123465252951528964587234675843733257549721655185850604449629037749798050593867815808367493607761074167203e-6),
        K(3.9367598891408415217939027108074887847673564509601917363278330920730719763184952569786096533271821267912464336083679709264700463e-8),
        K(6.1511873267825648778029740886343589006492244222098376295896480615672663845071908228080782331139431967897922330634666102008303163e-10),
        K(7.5940584281266233059295735665953083409170414099680216865037645091739362117435423650805312126318743025058577203204754258976678014e-12),
        K(7.5940584281266233059299735717776013802752093873942144943735559285381124771290874275774972580263377448832380260417572227398010273e-14),
        K(6.2760813455591928148129664638468741352313877508511272487605818260638787040214741939853005199347727650099145873785771162609470873e-16),
        K(4.3583898233049950103405756560043488004237104566659409254314063662323845080352835082013472245469229428130314480254874472462460891e-18),
        K(2.5789288895295828459033545114345754110126610917070986842981263284383138335900724842812739092696462831628856077304156728937367615e-20),
        K(1.3157800456783585979726388980387458854840810376002467594059415767303181634829326622143751699233015215909660383669969973896437826e-22),
        K(5.8479113141260380335975041071277913723234662896226038462639284334574178411617861304919195972489113065590753369218560772342930726e-25),
        K(2.2843403570804844741559425131287256347316408180969139215424927280400161659234508601992572475246984021932920101195981533775809011e-27),
        K(7.9042918930120215557655530936488590704671411306518762754706457865191035559604608171265261834240025123226121413351951869227469481e-30),
        K(2.439596263275438744691212721375550679559038036180973780451914524892069264667965537779114429543772487937962814842584430315822309e-32),
        K(6.7578843858008565460738956555556721994562125360817163818492848944114465748014992064357014275747373903974186259850477923287772269e-35),
        K(1.6894710964596677914300599135928934608666852785592175435009841929849562545281838798552390215907162912732614665877854602723350365e-37),
        K(3.8310002185198913939699183167260082132386992335359252310560205887083199288222921373551164604122527855575096799590794332372889982e-40),
        K(7.9152897117050638868250808160461112764395266893341033364145136187242187903457065607810464366513129569605469553757332577108931927e-43),
        K(1.4962740408949741994371775249428502243459110838196377209217623910923455903816887608569441907818721649385471112081388662207931364e-45),
        K(2.5976980842076886545162534185856248239556535520445635413487233605795098016339073754464291911162237045775316398305220029774563812e-48),
        K(4.1563157246970089177790104414625258959749615665778227563200792024682955941012716142727278222955329547858593509607233532321977384e-51),
        K(6.1484091972930398448317003607725995381550723813563853679777776077080177504939862858361387931698777740352242553154617386041598546e-54),
        K(8.4339079896649387396805165150345921345758626392364604080732341511582436140021851014611248449264294960348295127030417414763783526e-57),
        K(1.0758542181237264915749101505104765885049322002226174095616335030835886896108338436649327399785636341219913669804985122152824265e-59),
        K(1.2785524522454856200134470274184700289665086209215001086254221172250431156408812081801869424331341843840158057362920094547656938e-62),
        K(1.4248179766163819508194806049285708797690012879420281636927945272874304904429404793261110930103216776338690899423424418356698994e-65),
        K(1.4613184522449821802557940850011777432583484917091332887731152191646500305736827143107335993192474581511762327392034830200299768e-68),
        K(1.517919158523286579991864653430437909257242986598610110778281992820044480031975647280367518640551811313326095597412581231791542e-71),
        K(1.0730639360074958234700168111189609812309927776575154251790558216554191498594054258151950239282462897703667491995027177113192561e-74),
        K(1.849733759307828545234442543657764632352967627228530812064253357302934215881683992754041627944673908746817007111464947383730042e-77),
        K(-5.8729229858051647671731463714537532835823360515481187134639456051237156797256723676348818396781843804444226329124520413341278697e-81),
        K(3.1842186240879590225313774846938043014439524114338193689962549518448769929496194164836586422709285896399036894999980974150749013e-83),
        K(-2.1881719673472495157660923989130852804625696106789978549832744474725166266000647439983140838502086580706730083730364920779307716e-86),
        K(2.0535558773071560722118718580083555594558665120203432096590057810260358624306996796538168725229696967831406758396551788353456152e-89),
    };
    static const R Q1[] =
    {
        K(1),
    };
    static const R P2[] =
    {
        K(0.39894228040143267793994605993438189141018035442705729869804407689568381581734064686884602636060801485954621762740457550454120697),
        K(0.049867785050179084742493257490426145214512545591609659268248950196534280305961469882423461068790635445162491599891668241427694247),
        K(0.028050629090725735167652470998240899244071027169401361478749442415263988321141666521103085221687123192693334717091473761538817597),
        K(0.029219405302839307466250377876541919632987444632087277804127165783408826159066976653503074036487303121478723913535606767065816568),
        K(0.044742214369972689672757506995683668911838579935127513656738903414493578300728048844259081530992294260985826760124616165540729805),
        K(0.090602984099194545717108245690464709907435886237578028801057472781356468016933973445902282489417082489842418887592185592440521463),
        K(0.22839502241685343551568202900938534408323646867712940217756161458845324824124440104346632880269277474229622722857526929839330408),
        K(0.68926354970839797649853996399096479758550548316271823261307182611573812144139459011679425469995979044875363036072019650427109275),
        K(2.423192207724660758719355443491407013949377273831782675793140405839150024194700004683477388942832249119599378305563068242785029),
        K(9.7264091894977322855634221908603466759090961796678318811908721617317858636550474385646395183519977007947556596807399846944242623),
        K(43.894849749018396087574371484153328058059670508806822372081594668827438357942108803133491029755284933660373288867373132032368835),
        K(218.92224061880098667718426029789638820474517627422355720085391021586878854687601367885067299893531247137789283401227975375634307),
        K(1411.0714378841016321410168079012058915863951925541781772306639164464058932052840516833410793638465388787557712463940136004684829),
        K(-24669.417143852155941925922310645066937535848529101153480534808775058200117108762353716556845844109846670042694143992197668096542),
        K(4.3456187215737626225180053878439968236045760017182272244210822162169787026396322773165192567887238575971635413521445809876980004e6),
        K(-4.8794244314757594930125129712472496055896630924509326926769483634245544755261389240465702901747116533723557473550656594631532836e8),
        K(4.7135622802198833573595004986499079081107263879749738524409185091193774357674279922423541892450668395500169866908527469556459572e10),
        K(-3.8845767816711551934241716931447885535253872729659980739269139137736111723828410537289978592243029663663867611966517996124536087e12),
        K(2.7431597072135274612373797235764178284960675893163510470039767771950309546569261187500692221566731027281621740913265115102310496e14),
        K(-1.6636598713596741944604248189852920881501346018142892635481495198437840156774926185383215014490693563834394957950230048940030269e16),
        K(8.6761557906756247348638118281136603405500414721831534457612925481123322024593371759387146931591096846424375193986594474897877506e17),
        K(-3.8913176955569588730031997156307239438827348326480575220510931691527300070325174335141194550785118636626722266446637917200705305e19),
        K(1.4994387859623891999511364539836016709264090275867767459706869909005029610415881414010795232967045454552995180823348444416265002e21),
        K(-4.9526316726474597668971927656830715565011030190318549395742059758921820863775022631233335917983607537011812481817492614325105297e22),
        K(1.397064884652902423392749086230445982786614531772675459975976690233879103001760917728872951577136897092365891622780984709402973e24),
        K(-3.3476761876177350335523912452608780785827233458182277101687250287742574272022970793769978972224508627617386323901820464201569753e25),
        K(6.7638961792609295046405395451813410863531398189017078811550463864945380021979484321848645213471320221939039055909450745691397385e26),
        K(-1.14079573340166931747184530752445311008079545441036388481953511060965509254548158980345036081716115934507061629230171580151953e28),
        K(1.5843626737066263794105803441416923917274696837109592895480698389316466299714783209728565162325083868829306539854083641550288868e29),
        K(-1.7782859979199278216778586806675418350823245944594339589914373884037721073371915283109990177929377444043444423433544695887093251e30),
        K(1.5708853802374339215889671052352242220779637148891880270680676413885516142759102997242314905431899904918009952268613041942043472e31),
        K(-1.0500350399986598062815891870047338843859846801171623906896291522399576189691430301012749176131933567692272239383081900006983104e32),
        K(4.9857172763962740855781912335959709090990833630483482545791462951647915758552039690361109105750366582497609235282417166004295666e32),
        K(-1.4964974779697884725668300877703598873484331649647945507167852672912478283006052834187237212505395652815341504129608210970635177e33),
        K(2.1318910685882714113074231521717786024327977165170135019860427349218908241618243742167588675566067235117191319364522120898974407e33),
    };
    static const R Q2[] =
    {
        K(1.0),
    };
  #elif LDBL_MANT_DIG == 64
    /* long double 96 bit wide */
    static const R P1[] =
    {
        K(1.00696388290874250231638626673686646317801154370159972703168538),
        K(0.243352848727738955738908687369450214577920342918851509272408866),
        K(0.006964401160721188186398281247079919082283450941469460164590432),
        K(0.000083047334117897959145500056901191736030823531931380263686302),
        K(5.18256420384764810882467760619532575731801821889985626099e-7),
        K(1.90790611016475818883461118145629029943434891680660527e-9),
        K(4.44170587990105074420754325358582895345307949815573e-12),
        K(6.805150196466153819995090798791966304827558189423e-15),
        K(6.985104315031938858779570788468047860794936128e-18),
        K(4.785507068734939741097928056648844894386614e-21),
        K(2.117077490896605677726199140622837572025e-24),
        K(5.52919580174986488729896702518475621e-28),
        K(6.5666588969169003434516942087381e-32),
    };
    static const R Q1[] =
    {
        K(1.000011863675914860400478598182318948452642744176732473923183795358),
        K(-0.006896324225185339751945118908659032102601115341265224641280055928),
        K(0.00001186368725103095674191039208189536570056102297912907572890034),
        K(-1.3496696876875206170114218872940978330152830864048922600366e-8),
        K(1.1336174449932022007831556161183407392648142591512923821e-11),
        K(-7.422841655569707018644701823047396523966712871897812e-15),
        K(3.893669273036094904159100761627937282479044304922e-18),
        K(-1.651720545895290413869725701665259282787265642e-21),
        K(5.6323805635535562808481781647661895955232e-25),
        K(-1.50728066570923164846664348266414336373e-28),
        K(3.006044492319661074666639642233229e-32),
        K(-4.010148023149017379419706572721e-36),
        K(2.70282874465984817539266054e-40),
    };
    static const R P2[] =
    {
      K( 1.30090423521760256476093919023146864017751590623897710895862681),
      K(-1.981041925270972574120174940817336830170017871902975653312750388),
      K( 0.956892580228917795561363651877698243164566364537052353014543669),
      K(-0.295476285312266394050596510402082979039773201845265239542019439),
      K( 0.056978837924988815165935230495950981635872574537538384147033652),
      K(-0.006299149197554616295736173514236214970859775932020376086036399),
      K( 0.000353716966863384475462973243411450895641022763240331882363443),
      K(-8.707624424632528381900923003415938761710942641810978203625e-6),
      K( 7.283705999222063845686558855093093825421931939071318202e-8),
      K(-9.7967727386492889920273780071218382357131320542055799e-11)
    };
    static const R Q2[] =
    {
      K( 3.257608431020108786259398271424889402309379351594793640349460063),
      K(-4.96363276525502538609792324882976732173260916421950408693842731),
      K( 2.400495835659089927333294199555080092801133193497330702140754591),
      K(-0.742868968166381852162379299256973953894545292197384361382965225),
      K( 0.143801810439830068463911726822151703498931831272162081681231077),
      K(-0.016019224718850575023820322478614758671031103220377245057110662),
      K( 0.000914623505897601721718970098041677534130323750396506936364887),
      K(-0.000023411644633126949191317085153966622167096685843127825287574),
      K( 2.17705048674331703171406080664526952334380771487046428339e-7),
      K(-4.47580289731041130181939560179689655281441839562189718e-10)
    };
  #else
    #error Unsupported size of long double
  #endif
#elif defined(NFFT_SINGLE)
  /* float */
  static const R P1[] =
  {
      K(1.006634511033311726164163027592274220828216885723379609007274761),
      K(0.240606487720090757394176928596156553834296465200311569457994763),
      K(0.006634921274522227156198202198389031672287220144321235665461021),
      K(0.000073749622820821337100502174723273851941734199062726870961819),
      K(4.10243517822171814488230564074819973544765129449450710122e-7),
      K(1.262110026222369902633819303536802438120823461060572684e-9),
      K(2.218532296437410634454463125960648541194468552527652e-12),
      K(2.141504045536019682125761418851096299425878119158e-15),
      K(9.19584570350722374435337612379408707845677156e-19),
  };
  static const R Q1[] =
  {
      K(1.000022624782705275228334312456728477812835742762369533496905023937),
      K(-0.009614857078745003693609489751018087358244444264456521971379273084),
      K(0.000022624818652773047747424411495054891627754515915461183178099877),
      K(-3.4080521639954323706277061786236961377055349443081338572762e-8),
      K(3.5947512112800645225066705862453058797853924958888263259e-11),
      K(-2.7149805873212658218594464017972758572144265290831215e-14),
      K(1.4293388301569282795540255590126107486209476445158e-17),
      K(-4.771887851505849942903948600229238419570937509e-21),
      K(7.68298982666756594543081799488936861257839e-25),
  };
  static const R P2[] =
  {
    K( 0.400758393969643840397216812932361963736749407866811083462461),
    K(-0.0312216150704950438088565774064329777860642477326179964345542),
    K( 0.0001215451718646727844117193541329442989170354233955281424116)
  };
  static const R Q2[] =
  {
    K( 1.00043733569136882353241680221279480297575523819814430369272934),
    K(-0.0822433017391967535749382764476705160129315137731445852657631),
    K( 0.00043733569136882353241680221279480297575523819814430369272934)
  };
#else
  /* double */
  static const R P1[] =
  {
      K(1.006897990143384859657820271920512961153421109156614230747188622),
      K(0.242805341483041870658834102275462978674549112393424086979586278),
      K(0.006898486035482686938510112687043665965094733332210445239567379),
      K(0.000081165067173822070066416843139523709162208390998449005642346),
      K(4.95896034564955471201271060753697747487292805350402943964e-7),
      K(1.769262324717844587819564151110983803173733141412266849e-9),
      K(3.936742942676484111899247866083681245613312522754135e-12),
      K(5.65030097981781148787580946077568408874044779529e-15),
      K(5.267856044117588097078633338366456262960465052e-18),
      K(3.111192981528832405775039015470693622536939e-21),
      K(1.071238669051606108411504195862449904664e-24),
      K(1.66685455020362122704904175079692613e-28),
  };
  static const R Q1[] =
  {
      K(1.000013770640886533569435896302721489503868900260448440877422679934),
      K(-0.007438195256024963574139196893944950727405523418354136393367554385),
      K(0.000013770655915064256304772604385297068669909609091264440116789601),
      K(-1.6794623118559896441239590667288215019925076196457659206142e-8),
      K(1.50285363491992136130760477001818578470292828225498818e-11),
      K(-1.0383232801211938342796582949062551517465351830706356e-14),
      K(5.66233115275307483428203764087829782195312564006e-18),
      K(-2.44062252162491829675666639093292109472275754e-21),
      K(8.15441695513966815222186223740016719597617e-25),
      K(-2.01117218503954384746303760121365911698e-28),
      K(3.2919820158429806312377323449729691e-32),
      K(-2.70343047912331415988664032397e-36),
  };
  static const R P2[] =
  {
    K( 0.4305671332839579065931339658100499864903788418438938270811),
    K(-0.2897224581554843285637983312103876003389911968369470222427),
    K( 0.0299419330186508349765969995362253891383950029259740306077),
    K(-0.0010756807437990349677633120240742396555192749710627626584),
    K( 0.0000116485185631252780743187413946316104574410146692335443),
    K(-1.89995137955806752293614125586568854200245376235433e-08)
  };
  static const R Q2[] =
  {
    K(1.0762291019783101702628805159947862543863829764738274558421),
    K(-0.7279167074883770739509279847502106137135422309409220238564),
    K(0.0762629142282649564822465976300194596092279190843683614797),
    K(-0.0028345107938479082322784040228834113914746923069059932628),
    K(0.0000338122499547862193660816352332052228449426105409056376),
    K(-8.28850093512263912295888947693700479250899073022595e-08)
  };
#endif

static const INT N1 = sizeof(P1)/sizeof(P1[0]);
static const INT M1 = sizeof(Q1)/sizeof(Q1[0]);
static const INT N2 = sizeof(P2)/sizeof(P2[0]);
static const INT M2 = sizeof(Q2)/sizeof(Q2[0]);

static inline R evaluate_chebyshev(const INT n, const R *c, const R x)
{
  R a = c[n-2], b = c[n-1], t;
  INT j;
  
  A(n >= 2);
  
  for (j = n - 2; j > 0; j--)
  {
    t = c[j-1] - b;
    b = a + K(2.0) * x * b;
    a = t;
  }
  return a + x * b;
}

static inline R evaluate_polynomial(const INT n, const R *c, const R x)
{
  R r = c[n-1];
  INT j;

  A(n >= 2);

  for (j = n - 2; j >= 0; j--)
  {
    r = r * x + c[j];
  }

  return r;
}

R Y(bessel_i0)(R x)
{
  if (x < 0)
  {
    /* even function */
    x = -x;
  }

  if (x == K(0.0))
    return K(1.0);

#if defined(NFFT_LDOUBLE) && LDBL_MANT_DIG > 64
  if (x <= K(25.0))
  {
    R y = (x / K(2.0));
    y = y * y;

    return (K(1.0) + y * evaluate_polynomial(N1, P1, y));
  }
  else
  {
    return (EXP(x) / SQRT(x)) * evaluate_polynomial(N2, P2, 1 / x);
  }
#else
  if (x <= K(15.0))
  {
    /* x in (0, 15] */
    const R y = x * x;
    return evaluate_chebyshev(N1, P1, y) / evaluate_chebyshev(M1, Q1, y);
  }
  else
  {
    /* x in (15, \infty) */
    const R y = (K(30.0) - x) / x;
    return (EXP(x) / SQRT(x)) * (evaluate_chebyshev(N2, P2, y) /
      evaluate_chebyshev(M2, Q2, y));
  }
#endif
}