/*!\file cairy.c
   \brief compute complex airy function

   [1]  "Handbook of Mathematical Functions with
        Formulas, Graphs, and Mathematical Tables"
        Milton Abramowitz and Irene A. Stegun
        public domain (work of US government)
        online http://www.math.sfu.ca/~cbm/aands/
	
   [2]  "Computation of complex Airy functions and their zeros 
        using asymptotics and the differential equation"
	Fabijonas, B. R., Lozier, D. W., and Olver, F. W.   
        ACM Trans. Math. Softw. 30, 4 (Dec. 2004), 471-490. 
	http://dx.doi.org/10.1145/1039813.1039818

   [3]  "Computing complex Airy functions by numerical quadrature"
        Gil A., J. Segura, N.M. Temme. 
        Numerical Algorithms. 
	Volume 30, Number 1, May 2002 , pp. 11-23(13)
	http://dx.doi.org/10.1023/A:1015636825525
	http://citeseer.ist.psu.edu/gil01computing.html  
	[online 20070922]

   [4]  "Airy Functions" 
        Weisstein, Eric W. 
	From MathWorld -- A Wolfram Web Resource. 
	http://mathworld.wolfram.com/AiryFunctions.html	
	
   [5]  "Computing Complex Airy Functions by Numerical Quadrature"
        Amparo Gil1, Javier Segura, and Nico M. Temme
	Numerical Algorithms (issn 1017-1398)
	volume 30, Number 1, may 2002, pages 11-23
	htpp://dx.doi.org/10.1023/A:1015636825525
	http://www.cwi.nl/ftp/CWIreports/MAS/MAS-R0117.pdf
	[online 20070922]

   [6] "Bessel, Airy, Struve Functions, AiryAiPrime[z], Series representations,
        Generalized power series, Expansions at z==0"
	wolfram research
	http://functions.wolfram.com/03.07.06.0002.01
	[online 20070922]

*/

/*
** Copyright (C) 2007 bastien ROUCARIES 
**
** This library is free software; you can redistribute it and/or
** modify it under the terms of the GNU Lesser General Public
** License as published by the Free Software Foundation; either
** version 2.1 of the License, or (at your option) any later version.
**
** This library 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
** Lesser General Public License for more details.

** You should have received a copy of the GNU Lesser General Public
** License along with this library; if not, write to the Free Software
** Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
**
*/

/* Please send comments to roucaries.bastien+mathsoft at gmail.com
 *
 * ABOUT SPAM and SPAM BOT
 * If you want to send me spam (you don't and by UE law forbid you to do so)
 * please send me a mail to my kill list address 
 * address@hidden
 * END SPAM and SPAM BOT section
 */



#include "cairy.h"

#include <complex.h>
#include <float.h>
#include <math.h>
#include <limits.h>

#include <fenv.h>

#include <stdbool.h>

#include <errno.h>

#include <assert.h>
#include <stdlib.h>
#include <stdio.h>

#define fmtc "%+.8e%+.8e*I"
#define efmtc "%+.16e%+.16e*I"
#define cprint(z) creal(z), cimag(z)

//#define DEBUG printf
#define DEBUG(x,...) 

#define sqr(z) ((z)*(z))

#define norm2(z) (sqr(creal(z))+sqr(cimag(z)))
#define ARRAY_SIZE(x) (sizeof(x)/sizeof(x[0]))

#define max(x,y) (x) > (y) ? (x) : (y)

/*!\brief Force a compilation error if condition is true */
#define BUILD_BUG_ON(condition) ((void)sizeof(char[1 - 2*!!(condition)]))

/*!\brief ASSERT but a build time */
#define BUILD_ASSERT(condition) BUILD_BUG_ON(!(condition))

/*!\brief DBL_MAX^(2/3) */
#define DBL_MAX_POW_2_3 3.18525133652251465091489984275398832987840324458967e205

/*!\brief square root of 3 (\f$\sqrt(3)\f$) */
#define M_SQRT3 1.7320508075688772935274463415058723669428

/*!\brief square root of 3 divided by two (\f$\frac{\sqrt(3)}{2}\f$) */
#define M_SQRT3_2 0.86602540378443864676372317075293618347140262690519

/*!\brief exponential of \f$i\pi/3\f$ (\f$e^{i\frac{\pi}{3}}\f$) */
#define M_EXP_I_PI_3 (0.5 + I * M_SQRT3_2)

/*!\brief exponential of \f$2i\pi/3\f$ (\f$e^{i\frac{2\pi}{3}}\f$) */
#define M_EXP_I_2PI_3 (-0.5 + I * M_SQRT3_2)

/*!\brief exponential of \f$-2i\pi/3\f$ (\f$-e^{i\frac{2\pi}{3}}\f$) */
#define M_EXP_MI_2PI_3 (-0.5 - I * M_SQRT3_2)  

/*!\brief exponential of \f$i\pi/6\f$ (\f$e^{i\frac{\pi}{6}}\f$) */
#define M_EXP_I_PI_6 (M_SQRT3_2 + I * 0.5)

/*!\brief exponential of \f$-i\pi/6\f$ (\f$e^{-i\frac{\pi}{6}}\f$) */
#define M_EXP_MI_PI_6 (M_SQRT3_2 - I * 0.5)

/*!\brief The max of \f$|\Im\text{m}\;z|\f$ used for Mc Laurin Serie */
#define MAX_IMAG_SMALL 2.0

/*!\brief The min of \f$\Re\text{e}\;z\f$ used for Mc Laurin Serie */
#define MIN_REAL_SMALL -2.6

/*!\brief The max of \f$\Re\text{e}\;z\f$ used for Mc Laurin Serie */
#define MAX_REAL_SMALL 1.3

/*!\brief Minimum of \f$|z|\f$ in order to use asymptotic expansion 

  Use value given in [2, tab I pp 480] ie:
  - 5.394496440887451 for ieee single precision
  - 8.818819021035965 for ieee double precision
  - 14.4578           for ieee quadruple precision
*/
#define AIRY_LARGE 8.818819021035965

/*!\brief Tolerance for safety check in percent */
#define TOL_ABORT 0.1

#ifndef NDEBUG
#define ABORTNDEBUG() do { } while(0)
#else 
#define ABORTNDEBUG() abort()
#endif

/*!\brief Safety check for small argument function 
   \note for use in assert
*/
inline static bool
issmall (const complex z)
{
  double x, y;

  ABORTNDEBUG();

  x = creal (z);
  y = cimag (z);

  if (fabs (y) <= MAX_IMAG_SMALL * (1 + TOL_ABORT) &&
      x >= MIN_REAL_SMALL * (1 + TOL_ABORT) &&
      x <= MAX_REAL_SMALL * (1 + TOL_ABORT))
    return true;

  fprintf (stderr, "!issmall(" fmtc ")\n", creal (z), cimag (z));
  return false;
}

/*!\brief Safety check for medium argument function 
   \note for use in assert
*/
inline static bool
ismedium (const complex z)
{
  double x, y;

  x = creal (z);
  y = cimag (z);

  if (fabs (y) <= MAX_IMAG_SMALL * (1 - TOL_ABORT) &&
      x >= MIN_REAL_SMALL * (1 - TOL_ABORT) &&
      x <= MAX_REAL_SMALL * (1 - TOL_ABORT))
    return false;

  if(cabs(z) > AIRY_LARGE *( 1 + TOL_ABORT))
    return false;

  return true;
}

/*!\brief Safety check for large argument function 
   \note for use in assert
*/
inline static bool
islarge (const complex z)
{
  ABORTNDEBUG();
  return cabs(z) > AIRY_LARGE *( 1 - TOL_ABORT);
}

/*!\brief Test if the argument of \f$\arg z\f$ is in region I ie 
  than \f$|\arg z|< 2\pi/3\f$ 
  \note Should only be used in assert
*/
static inline bool isregionI(const complex z)
{
  if(cimag(z) > 0 && creal(z) < 0 && cimag(z) <=  -M_SQRT3 * creal(z))
    return false;
  if(cimag(z) < 0 && creal(z) < 0 && cimag(z) >=  M_SQRT3 * creal(z))
    return false;
  if(cimag(z) == 0 && creal(z) <= 0)
    return false;
  return true;
}

/*!\brief Test if the argument of \f$\arg z\f$ is in region II ie 
  than \f$|\arg z|> 2\pi/3\f$ 
  \note Should only be used in assert
*/
static inline bool isregionII(const complex z)
{
  ABORTNDEBUG();
  return !isregionI(z);
}


/*!\brief Test if the argument of \f$\arg z\f$ is in region A ie 
  than \f$|\arg z |< \pi/3\f$ 
  \note Should only be used in assert
*/
static inline bool isregionA(const complex z)
{
  ABORTNDEBUG();
  if(cimag(z) > 0 && creal(z) > 0 && cimag(z) <  M_SQRT3 * creal(z))
    return true;
  if(cimag(z) < 0 && creal(z) > 0 && cimag(z) >  -M_SQRT3 * creal(z))
    return true;
  if(cimag(z) == 0 && creal(z) >= 0)
    return true;
  return false;
}

/*!\brief Test if the argument of \f$\arg z\f$ is in region B ie 
  than \f$\pi/3 \le |\arg z| \le 2\pi/3\f$ 
  \note Should only be used in assert
*/
static inline bool isregionB(const complex z)
{
  return isregionI(z) && !isregionA(z);
}

/*!\brief Test if the argument of \f$\arg z\f$ is in region C ie 
  than \f$2\pi/3 \le |\arg z| \le \pi\f$ 
  \note Should only be used in assert
*/
static inline bool isregionC(const complex z)
{
  ABORTNDEBUG();
  return !isregionI(z);
}


/*!\brief Test is the argument of \f$z\f$ is larger than \f$\pi/3\f$ 

  Use than fact that if \f$\arg z>\pi/3\f$ and \f$\Im\text{m}\;z >0\f$
  then \f$R\text{e}\;z < 0\f$ or 
  \f$R\text{e}\;z > 0 \wedge \Im\text{m}\;z \ge \sqrt{3} \Re\text{e}\;z\f$
  \note use simplified formulation (negated)
  \attention assert \f$ \Im\text{m}\;z > 0 \f$

*/ 
static inline bool arglargerpi_3(const complex z)
{
  assert(cimag(z) >= 0);  

  if(creal(z) > 0 && cimag(z) <=  M_SQRT3 * creal(z))
    return false;
  return true;
}


/*!\brief Test if the argument of \f$\arg z\f$ is larger than \f$2\pi/3\f$ 

  Use than fact that if \f$\arg z>2\pi/3\f$ and \f$\Im\text{m}\;z >0\f$
  then  \f$R\text{e}\;z < 0 \wedge \Im\text{m}\;z < -\sqrt{3} \Re\text{e}\;z\f$
  \attention assert \f$ \Im\text{m}\;z > 0 \f$

*/ 
static inline bool arglargerpi2_3(const complex z)
{
  assert(cimag(z) >= 0);
  
  if(creal(z) < 0 && cimag(z) <  -M_SQRT3 * creal(z))
    return true;
  return false;
}


/*!\brief Compute \f$z^{3/2}\f$ assume no overflow
  \attention assert no overflow
*/
static complex cpow32fast(const complex z)
{
  assert(cabs(z) <= DBL_MAX_POW_2_3);
  return z * csqrt(z);
}

/*!\brief Compute \f$z^{3/2}\f$ 
  \note raise FE_OVERFLOW and set errno to ERANGE on 
  overflow
*/
static complex cpow32(const complex z)
{
  if(cabs(z) <= DBL_MAX_POW_2_3)
    return cpow32fast(z);

  /* overflow */
  errno = ERANGE;
  (void) feraiseexcept (FE_OVERFLOW);
  return cpow(z, 3.0/2.0);
}

/*!\brief Compute zeta assuming no overflow

  zeta is defined by:
  \[
  \zeta = \frac{2}{3} z^{3/2}
  \]
*/
static complex czetafast(const complex z)
{
  return (2.0/3.0) * cpow32fast(z);
}

/*!\brief Compute zeta

  zeta is defined by:
  \[
  \zeta = \frac{2}{3} z^{3/2}
  \]
  \note raise FE_OVERFLOW and set errno to ERANGE on 
  overflow
*/
static complex czeta(const complex z)
{
  if(cabs(z) <= DBL_MAX_POW_2_3)
    return czetafast(z);

  /* overflow */
  errno = ERANGE;
  (void) feraiseexcept (FE_OVERFLOW);
  return (2.0/3.0) * cpow(z, 3.0/2.0);
}

/*!\brief Compute \f$\exp \left[\frac{4}{3} z^{3/2}\right]\f$ */
static complex cexp_43z32(const complex z)
{
  return cexp((4.0/3.0) * cpow32(z));
}

/*!\brief Compute \f$z^{1/4}\f$ */
static complex cz1_4(const complex z)
{
  return csqrt(csqrt(z));
}

/*!\brief Computed coefficient of f auxilary function 

  According to (10.4.3) [1-p446]  auxiliary function \f$f\f$ 
  is:
  \f[
  f(z)=\sum_{k=0}^\infty 3^k \left(\frac{1}{3}\right)_k \frac{z^{3k}}{(3k)!} 
  \f]
  Using \f$3^k(\alpha+\frac{1}{3})_k = (3\alpha + 1)(3\alpha + 4)\dotsm(3\alpha + 3k-2)\f$:
  \f[
  f(z)=1+\frac{1}{3!} z^3 + \frac{1\times4}{6!} z^6 + \frac{1\times4\times7}{9!} z^9 + \dotsb
  \f]
  Therefore compute \f$1/3!\f$, \f$\frac{1\times4}{6!}\f$, etc...
  \note computed using a arbitrary precision program
*/
static const double fcoeff[] = 
{
	/*(1)/3! */
	0.1666666666666666666666666666666666666667e0,
	/*(1*4)/6! */
	0.5555555555555555555555555555555555555556e-2,
	/*(1*4*7)/9! */
	0.7716049382716049382716049382716049382716e-4,
	/*(1*4*7*10)/12! */
	0.5845491956603067714178825289936401047512e-6,
	/*(1*4*7*10*13)/15! */
	0.2783567598382413197228012042826857641672e-8,
	/*(1*4*7*10*13*16)/18! */
	0.9096626138504618291594810597473390985858e-11,
	/*(1*4*7*10*13*16*19)/21! */
	0.2165863366310623402760669189874616901395e-13,
	/*(1*4*7*10*13*16*19*22)/24! */
	0.3923665518678665584711357228033726270643e-16,
	/*(1*4*7*10*13*16*19*22*25)/27! */
	0.558926712062487975030107867241271548523e-19,
	/*(1*4*7*10*13*16*19*22*25*28)/30! */
	0.6424444966235493965863308818865190212926e-22,
	/*(1*4*7*10*13*16*19*22*25*28*31)/33! */
	0.6083754702874520801006921229985975580423e-25,
	/*(1*4*7*10*13*16*19*22*25*28*31*34)/36! */
	0.4828376748313111746830889865068234587637e-28,
	/*(1*4*7*10*13*16*19*22*25*28*31*34*37)/39! */
	0.3258014000211276482342030948089227117164e-31,
	/*(1*4*7*10*13*16*19*22*25*28*31*34*37*40)/42! */
	0.189199419292176334630779962142231539905e-34,
	/*(1*4*7*10*13*16*19*22*25*28*31*34*37*...*43)/45! */
	0.9555526226877592658120200108193512116407e-38,
	/*(1*4*7*10*13*16*19*22*25*28*31*34*37*...*46)/48! */
	0.4235605597020209511578102884837549696992e-41,
	/*(1*4*7*10*13*16*19*22*25*28*31*34*37*...*49)/51! */
	0.1661021802753023337873765837191195959605e-44,
	/*(1*4*7*10*13*16*19*22*25*28*31*34*37*...*52)/54! */
	0.5803710002631108797602256593959454785481e-48,
	/*(1*4*7*10*13*16*19*22*25*28*31*34*37*...*55)/57! */
	0.1818204888042327317544566602117623679662e-51,
	/*(1*4*7*10*13*16*19*22*25*28*31*34*37*...*58)/60! */
	0.5136172000119568693628719215021535818254e-55,
	/*(1*4*7*10*13*16*19*22*25*28*31*34*37*...*61)/63! */
	0.1314944188458670940509144704306588791156e-58,
	/*(1*4*7*10*13*16*19*22*25*28*31*34*37*...*64)/66! */
	0.3065137968435130397457213762952421424605e-62,
	/*(1*4*7*10*13*16*19*22*25*28*31*34*37*...*67)/69! */
	0.6532689617295674333881529759063131765996e-66,
	/*(1*4*7*10*13*16*19*22*25*28*31*34*37*...*70)/72! */
	0.1277912679439685902558984694652412317292e-69,
	/*(1*4*7*10*13*16*19*22*25*28*31*34*37*...*73)/75! */
	0.2302545368359794419025197648022364535661e-73,
	/*(1*4*7*10*13*16*19*22*25*28*31*34*37*...*76)/78! */
	0.3833741872060929768606722690679927631803e-77,
	/*(1*4*7*10*13*16*19*22*25*28*31*34*37*..*79)/81! */
	0.5916268321081681741677041189320875975005e-81,
	/*(1*4*7*10*13*16*19*22*25*28*31*34*37*...*82)/84! */
	0.8485754906887093720133449783879627043897e-85,
	/*(1*4*7*10*13*16*19*22*25*28*31*34*37*...*85)/87! */
	0.1134155961893490205845155009874315295896e-88,
	/*(1*4*7*10*13*16*19*22*25*28*31*34*37*...*88)/90! */
	0.1415925046059288646498320861266311230832e-92,
	/*(1*4*7*10*13*16*19*22*25*28*31*34*37*...*91)/93! */
	0.1654891358180561765425807458235520372641e-96,
	/*(1*4*7*10*13*16*19*22*25*28*31*34*37*...*94)/96! */
	0.181457385765412474279145554631087760158e-100,
	/*(1*4*7*10*13*16*19*22*25*28*31*34*37*...*97)/99! */
	0.1870309067876855022460786998877424862482e-104,
	/*(1*4*7*10*13*16*19*22*25*28*31*34*37*..*100)/102! */
	0.1815481525797762592177040379418971910777e-108,
	/*(1*4*7*10*13*16*19*22*25*28*31*34*37*..*103)/105! */
	0.1662528869778170871956996684449607976902e-112
};

/*!\brief Auxiliary function \f$f\f$ 
  
  According to (10.4.3) [1-p446]  auxiliary function \f$f\f$ 
  is:
  \f[
  f(z)=\sum_{k=0}^\infty 3^k \left(\frac{1}{3}\right)_k \frac{z^{3k}}{(3k)!} 
  \f]
  Using \f$3^k(\alpha+\frac{1}{3})_k = (3\alpha + 1)(3\alpha + 4)\dotsm(3\alpha + 3k-2)\f$:
  \f[
  f(z)=1+\frac{1}{3!} z^3 + \frac{1\times4}{6!} z^6 + \frac{1\times4\times7}{9!} z^9 + \dotsb
  \f]

  \param[in] z complex number
  \return f(z)
*/
static complex
f (const complex z)
{
  unsigned int k;

  complex sum;
  complex fk;

  complex z3;
  complex z3n = 1;

  z3 = z * z * z;

  /* f_0 = 1 */
  fk = 1;
  sum = 1;

  for (k = 0; k < ARRAY_SIZE (fcoeff); k++)
    {
      z3n *= z3;
      fk = z3n * fcoeff[k];

      if (norm2 (fk) < norm2 (sum) * sqr (DBL_EPSILON)
	  || norm2 (fk) < DBL_MIN)
	return sum;

      sum += fk;
      DEBUG ("fcoeff=%.8e fk=" fmtc " sum=" fmtc "\n",
	     fcoeff[k], creal (fk), cimag (fk), creal (sum), cimag (sum));
    }

  /* Do not converge */
  assert (0 && sizeof ("f do not converge!"));
  errno = ERANGE;
  (void) feraiseexcept (FE_INVALID);
  return sum;
}

/*!\brief Computed coefficient of g auxilary function 

  According to (10.4.3) [1-p446]  auxiliary function \f$f\f$ 
  is:
  \f[
  g(z)=\sum_{k=0}^\infty 3^k \left(\frac{2}{3}\right)_k \frac{z^{3k+1}}{(3k+1)!} 
  \f]
  Using \f$3^k(\alpha+\frac{1}{3})_k = (3\alpha + 1)(3\alpha + 4)\dotsm(3\alpha + 3k-2)\f$:
  \f[
  g(z)=z+\frac{2}{4!} z^4 + \frac{2\times5}{7!} z^7 + \frac{2\times5\times8}{9!} z^{10} + \dotsb
  \f]
  Therefore compute \f$2/4!\f$, \f$\frac{2\times5}{7!}\f$, etc...
  \note computed using a arbitrary precision program
*/
static const double gcoeff[] = 
{
	/*(2)/4!*/
	0.8333333333333333333333333333333333333333e-1,
	/*(2*5)/7!*/
	0.1984126984126984126984126984126984126984e-2,
	/*(2*5*8)/10!*/
	0.2204585537918871252204585537918871252205e-4,
	/*(2*5*8*11)/13!*/
	0.141319585764030208474652919097363541808e-6,
	/*(2*5*8*11*14)/16!*/
	0.5888316073501258686443871629056814241999e-9,
	/*(2*5*8*11*14*17)/19!*/
	0.1721729846052999615919260710250530480117e-11,
	/*(2*5*8*11*14*17*20)/22!*/
	0.3726687978469696138353378160715433939647e-14,
	/*(2*5*8*11*14*17*20*23)/25!*/
	0.6211146630782826897255630267859056566078e-17,
	/*(2*5*8*11*14*17*20*23*26)/28!*/
	0.8215802421670405948750833687644254717034e-20,
	/*(2*5*8*11*14*17*20*23*26*29)/31!*/
	0.8834196152333769837366487836176617975306e-23,
	/*(2*5*8*11*14*17*20*23*26*29*32)/34!*/
	0.7873615109031880425460327839729606038597e-26,
	/*(2*5*8*11*14*17*20*23*26*29*32*35)/37!*/
	0.5911122454228138457552798678475680209157e-29,
	/*(2*5*8*11*14*17*20*23*26*29*32*35*38)/40!*/
	0.3789181060402652857405640178510051416126e-32,
	/*(2*5*8*11*14*17*20*23*26*29*32*35*38*41)/43!*/
	0.2098106899447759057256722136495045080911e-35,
	/*(2*5*8*11*14*17*20*23*26*29*32*35*38*41*44)/46!*/
	0.1013578212293603409302764317147364773387e-38,
	/*(2*5*8*11*14*17*20*23*26*29*32*35*38*41*..*47)/49!*/
	0.4309431174717701570164814273585734580728e-42,
	/*(2*5*8*11*14*17*20*23*26*29*32*35*38*41*...*50)/52!*/
	0.1624974047781938751947516694413927066639e-45,
	/*(2*5*8*11*14*17*20*23*26*29*32*35*38*41*...*53)/55!*/
	0.5471293090174877952685241395333087766463e-49,
	/*(2*5*8*11*14*17*20*23*26*29*32*35*38*41*...*56)/58!*/
	0.1654958587469715049209086931437715597841e-52,
	/*(2*5*8*11*14*17*20*23*26*29*32*35*38*41*...*59)/61!*/
	0.4521744774507418167237942435622173764593e-56,
	/*(2*5*8*11*14*17*20*23*26*29*32*35*38*41*...*62)/64!*/
	0.112146447780441918830306111994597563606e-59,
	/*(2*5*8*11*14*17*20*23*26*29*32*35*38*41*...*65)/67!*/
	0.2536102392140251443471418181696010031795e-63,
	/*(2*5*8*11*14*17*20*23*26*29*32*35*38*41*...*68)/70!*/
	0.525072959035248746060335027266254664968e-67,
	/*(2*5*8*11*14*17*20*23*26*29*32*35*38*41*...*71)/73!*/
	0.9989972584384489080295567489845027872297e-71,
	/*(2*5*8*11*14*17*20*23*26*29*32*35*38*41*...*74)/76!*/
	0.1752626769190261242157117103481583837245e-74,
	/*(2*5*8*11*14*17*20*23*26*29*32*35*38*41*...*77)/79!*/
	0.2844249868857937750985259823890918268817e-78,
	/*(2*5*8*11*14*17*20*23*26*29*32*35*38*41*...*80)/82!*/
	0.4282219013637364876521017500588555056936e-82,
	/*(2*5*8*11*14*17*20*23*26*29*32*35*38*41*...*83)/85!*/
	0.5997505621340847165995822829955959463496e-86,
	/*(2*5*8*11*14*17*20*23*26*29*32*35*38*41*...*86)/88!*/
	0.7833732525262339558510740373505694179071e-90,
	/*(2*5*8*11*14*17*20*23*26*29*32*35*38*41*...*89)/91!*/
	0.956499697834229494323655723260768520033e-94,
	/*(2*5*8*11*14*17*20*23*26*29*32*35*38*41*...*92)/94!*/
	0.1094142871006897156627380145573974513879e-97,
	/*(2*5*8*11*14*17*20*23*26*29*32*35*38*41*...*95)/97!*/
	0.1174981605462733200845554280040780191021e-101,
	/*(2*5*8*11*14*17*20*23*26*29*32*35*38*41*...*98)/100!*/
	0.1186850106528013334187428565697757768708e-105,
	/*(2*5*8*11*14*17*20*23*26*29*32*35*38*41*...*101)/103!*/
	0.1129687898846386192830219460972546895782e-109,
	/*(2*5*8*11*14*17*20*23*26*29*32*35*38*41*...*104)/106!*/
	0.1014993619808073848005588015249368280127e-113
};


/*!\brief Auxiliary function \f$g\f$ 
  
  According to (10.4.3) [1-p446]  auxiliary function \f$f\f$ 
  is:
  \f[
  g(z)=\sum_{k=0}^\infty 3^k \left(\frac{2}{3}\right)_k \frac{z^{3k+1}}{(3k+1)!} 
  \f]
  Using \f$3^k(\alpha+\frac{1}{3})_k = (3\alpha + 1)(3\alpha + 4)\dotsm(3\alpha + 3k-2)\f$:
  \f[
  g(z)=z+\frac{2}{4!} z^4 + \frac{2\times5}{7!} z^7 + \frac{2\times5\times8}{9!} z^{10} + \dotsb
  \f]

  \param[in] z complex number
  \return f(z)
*/
static complex
g (const complex z)
{
  unsigned int k;

  complex sum;
  complex gk;

  complex z3;
  complex z3np1;

  z3 = z * z * z;

  sum = z;
  z3np1 = z;

  for (k = 0; k < ARRAY_SIZE (gcoeff); k++)
    {
      z3np1 *= z3;
      gk = z3np1 * gcoeff[k];

      if (norm2 (gk) < norm2 (sum) * sqr (DBL_EPSILON)
	  || norm2 (gk) < DBL_MIN)
	return sum;

      sum += gk;
      DEBUG ("gcoeff=%.8e gk=" fmtc " sum=" fmtc "\n",
	     gcoeff[k], creal (gk), cimag (gk), creal (sum), cimag (sum));
    }

  /* Do not converge */
  assert (0 && sizeof ("g do not converge!"));
  errno = ERANGE;
  (void) feraiseexcept (FE_INVALID);
  return sum;
}


/*!\brief \f$Ai(0)\f$

  Constant  used in computation of \f$Ai\f$

  According to (10.4.4) [1-p446]:
  \f[
   Ai(0) = \frac{1}{3^{2/3} \Gamma(2/3)}
  \f]
*/
static const double ai0 =
  0.35502805388781723926006318600418317639797917419918;

/*!\brief \f$A'i(0)\f$
  
  Constant used in computation of \f$Ai\f$ by definition:
  \f[
  Ai'(0) = - \frac{1}{3^{1/3} \Gamma(\frac{1}{3})}
  \f]
*/
static const double aip0 =
  -0.25881940379280679840518356018920396347909113835493;


/*!\brief Compute airy function of first kind  \f$Ai\f$ 
          for small argument 

	  Compute airy function of first kind  \f$Ai\f$ 
          for small argument using the relation 
	  (10.4.2) [1-p446]:
	  \f[
	  Ai=Ai(0) f(z) + A'i(0) g(z)
	  \f]
	  \param[in] z: Where to evaluate \f$Ai\f$
	  \return \f$Ai(z)\f$
*/
static complex
cairyAi_smallarg (const complex z)
{
  assert (issmall (z));

  /* underflow */
  if (norm2 (z) < sqr (DBL_MIN))
    {
      DEBUG ("Underflow\n");
      errno = ERANGE;
      (void) feraiseexcept (FE_UNDERFLOW);
      return ai0;
    }

  return ai0 * f (z) + aip0 * g (z);
}


/*!\brief Compute scaled airy function of first kind  \f$Aie\f$ 
          for small argument

  Compute the scalled airy function of first kind \f$Aie(z)=e^{\zeta}Ai\f$ where 
  \f$\zeta\f$ is \f$\frac{2}{3}z^\frac{3}{2}\f$	  

  \note Use non scaled version
*/
static complex
cairyAie_smallarg (const complex z)
{
  complex zeta;
  complex expzeta, nonscaled;

  zeta = czetafast(z);

  nonscaled = cairyAi_smallarg(z);
  expzeta = cexp(zeta);
  DEBUG("zeta=" fmtc  " non scaled=" fmtc " exp(zeta)=" fmtc "\n",
	cprint(zeta), cprint(nonscaled), cprint(expmzeta));
  return nonscaled * expzeta;
}

/*!\brief Constant used in computation of \f$Bi\f$

  According to (10.4.4) [1-p446]:
  \f[
  Bi(0) = \sqrt{3} \left[3^{-2/3} \Gamma(2/3)\right]
  \f]
*/
static const double bi0 =
  0.61492662744600073515092236909361355359472818864860;

/*!\brief Constant used in computation of \f$Bi\f$

  According to (10.4.5) [1-p446]:
  \f[
  Bi'(0) = - \sqrt{3} \left[ 3^{-2/3} \Gamma(2/3) \right]
  \f]
*/
static const double bip0 =
  0.44828835735382635791482371039882839086622679921226;

/*!\brief Compute airy function of second kind  \f$Bi\f$ 
          for small argument 

	  Compute airy function of first kind  \f$Bi\f$ 
          for small argument using the relation 
	  (10.4.2) [1-p446]:
	  \f[
	  Bi=Bi(0) f(z) + Bi'(0) f(z)
	  \f]
*/
static complex
cairyBi_smallarg (const complex z)
{
  assert (issmall (z));

  /* underflow */
  if (norm2 (z) < sqr (DBL_MIN))
    {
      DEBUG ("Underflow\n");
      errno = ERANGE;
      (void) feraiseexcept (FE_UNDERFLOW);
      return bi0;
    }

  return  bi0 * f (z) + bip0 * g (z);
}

/*!\brief Compute scaled airy function of second kind  \f$Bie\f$ 
          for small argument

  Compute the scalled airy function of second kind \f$Bie(z)\f$ is defined as:
  \f[
  Bie(z) = \begin{cases}
  e^{-\zeta} Bi(z) & \text{$|\arg z| \le \frac{\pi}{3}$} \\
  e^{\zeta} Bi(z) & \text{else}
  \end{cases}
  \f]
  where \f$\zeta\f$ is \f$\frac{2}{3}z^\frac{3}{2}\f$	  

  \note Use non scaled version
*/
static complex
cairyBie_smallarg (const complex z)
{
  complex zeta;
  complex exppmzeta, nonscaled;

  zeta = czetafast(z);

  nonscaled = cairyAi_smallarg(z);
  if(arglargerpi_3(z))
    exppmzeta = cexp(-zeta);
  else
    exppmzeta = cexp(zeta);

  DEBUG("zeta=" fmtc  " non scaled=" fmtc " exp(zeta)=" fmtc "\n",
	cprint(zeta), cprint(nonscaled), cprint(exppmzeta));
  return nonscaled * exppmzeta;
}


/*\brief wiegh

\todo regenerate these data usign http://www.alglib.net/integral/gq/glaguerre.php
*/
static double x[25] = { .0283891417994567679, .170985378860034935,
			    .43587167834177046, .823518257913030858, 1.33452543254227372,
			    1.96968293206435071, 2.72998134002859938, 3.61662161916100897,
			    4.63102611052654146, 5.77485171830547694, 7.05000568630218682,
			    8.45866437513237792, 10.0032955242749393, 11.6866845947722423,
			    13.5119659344693551, 15.482659695937714, 17.6027156808069112,
			    19.8765656022785451, 22.309185677396278, 24.9061720212974207,
			    27.673832073949719, 30.6192963295084111, 33.7506560850239946,
			    37.0771349708391198, 40.6093049694341322
    };
static double w[25] = { .143720408803313866, .230407559241880881,
			    .242253045521327626, .203636639103440807, .14376063062292141,
			    .086912883470607812, .0454175001832915883, .0206118031206069497,
			    .00814278821268606972, .00280266075663377634,
			    8.40337441621719716e-4, 2.1930373290776502e-4,
			    4.9740165900925776e-5, 9.78508095920717661e-6,
			    1.66542824603725563e-6, 2.44502736801316287e-7,
			    3.08537034236207072e-8, 3.33296072940112245e-9,
			    3.06781892316295828e-10, 2.39331309885375719e-11,
			    1.57294707710054952e-12, 8.64936011664392267e-14,
			    3.94819815638647111e-15, 1.48271173082850884e-16,
			    4.53390377327054458e-18
};


/*!\brief coefficient \f$a(z)\f$ 

  Defined in [5, pp 10] eq (6.3) as :
  \f[
  azcoeff = \frac{1}{\sqrt{\pi}48^{1/6}\Gamma(\frac{5}{6})}
  \f]
*/
static const double azcoeff = 0.2621839970883229496786247788550868016857;


/*!\brief Compute scaled airy function (\f$Aie\f$) for medium argument 
  and in first quadrant (\f$0 < \arg z < \pi/2\f$)

   Use integral representation of Airy function [3, pp3] eq (3.3) or 
   [5, pp 8] eq (6.3)
   \f{align*}
   Ai(z) &= a(z) \int_0^\infty (2+\frac{t}{\zeta})^{-1/6} t^{-1/6} e^{-t}\;\text{d}t &
   a(z) &= \frac{1}{\sqrt{\pi}48^{1/6}\Gamma(\frac{5}{6})} e^{-\zeta} \zeta^{-1/6}
   \f}
   Therefore:
   \f{align*}
   AiE(z) &= a(z) \int_0^\infty (2+\frac{t}{\zeta})^{-1/6} t^{-1/6} e^{-t}\;\text{d}t &
   a(z) &= \frac{1}{\sqrt{\pi}48^{1/6}\Gamma(\frac{5}{6})} \zeta^{-1/6}
   \f}
   Compute using a Gauss Laguerre quadrature
   \attention assert \f$0 < \arg z < \pi/2\f$
*/     
static complex cairyAie_mediumarg_Q1(const complex z)
{
  complex zeta, zetam16;
  complex az;
  double t;
  complex integrand;
  complex sum;

  unsigned int k;
  
  assert(creal(z) >= 0);
  assert(cimag(z) >= 0);
  assert(ismedium(z));
  BUILD_ASSERT(ARRAY_SIZE(w) == ARRAY_SIZE(x));

  zeta = czetafast(z);
  zetam16 = cpow(zeta, - 1.0 / 6.0);
  /* here az is az scaled ie without \f$e^{-\zeta}\f$ */
  az = azcoeff * zetam16;

  sum = 0;
  for(k = 0; k < ARRAY_SIZE(x); k++)
  {
    t = x[k];
    integrand = cpow( 2 + t/zeta, -1.0 /6.0);
    sum += integrand * w[k];
  }

  return sum * az;
}                       

/*!\brief Compute scaled airy function (\f$Aie\f$) for medium argument 
  and in second quadrant (\f$\pi/3 \ge \arg z > \pi/2\f$)

   Use integral representation of Airy function [3, pp3] eq (3.3) or 
   [5, pp 8] eq (6.3)
   \f{align*}
   Ai(z) &= a(z) \int_0^\infty (2+\frac{t}{\zeta})^{-1/6} t^{-1/6} e^{-t}\;\text{d}t &
   a(z) &= \frac{1}{\sqrt{\pi}48^{1/6}\Gamma(\frac{5}{6})} e^{-\zeta} \zeta^{-1/6}
   \f}
   But this representation is not numerically stable for \f$\pi/3 \ge \arg z > \pi/2\f$
   The reason (see [5, pp9] is the presence of singularities in the integrand at 
   \f$t=-2\zeta\f$. This problem could be solved by turning the integration path over a 
   given angle \f$\tau=\frac{3}{2}(\arg z -\frac{\pi}{2})\f$, that can be done using 
   the substituation \f$t\rightarrow (t+i\tan \tau)\f$. Therefore:
   \f{align*}
   Ai(z) &= a(z) \left(\frac{e^{i\tau}}{\cos \tau}\right)^{5/6}
         \int_0^\infty (2+\frac{t}{\tilde{\zeta}})^{-1/6} 
	       e^{-it\tan\tau} t^{-1/6} e^{-t}\;\text{d}t &
   a(z) &= \frac{1}{\sqrt{\pi}48^{1/6}\Gamma(\frac{5}{6})} e^{-\zeta} \zeta^{-1/6}
   \f}
   Where \f$\tilde\zeta=\cos \tau e^{-i\tau} \zeta\f$
   Therefore:
   \f{align*}
   Aie(z) &= a(z) \left(\frac{e^{i\tau}}{\cos \tau}\right)^{5/6}
         \int_0^\infty (2+\frac{t}{\zeta})^{-1/6} 
	       e^{-it\tan\tau} t^{-1/6} e^{-t}\;\text{d}t &
   a(z) &= \frac{1}{\sqrt{\pi}48^{1/6}\Gamma(\frac{5}{6})} \zeta^{-1/6}

   \f}

   Compute using a Gauss Laguerre quadrature
   \attention assert \f$0 < \arg z < \pi/2\f$
*/     
static complex cairyAie_mediumarg_Q2(const complex z)
{
  complex zeta, zetam16 ,zetatilde;
  double tau;
  complex eitaucostau16;
  complex az;
  double t;
  complex integrand;
  complex sum;

  unsigned int k;
  
  assert(cimag(z) >= 0);
  assert(ismedium(z));
  assert(creal(z) <= 0 && !arglargerpi2_3(z));
  BUILD_ASSERT(ARRAY_SIZE(w) == ARRAY_SIZE(x));

  tau = 3.0 / 2.0 * (carg(z) - M_PI_2);
  eitaucostau16 = cpow(cexp(I * tau)/cos(tau), 5.0 / 6.0);


  zeta = czetafast(z);
  zetam16 = cpow(zeta, - 1.0 / 6.0);
  zetatilde = cos(tau) * cexp(-I * tau) * zeta;
  /* here az is az scaled ie without \f$e^{-\zeta}\f$ */
  az = azcoeff * zetam16;

  sum = 0;
  for(k = 0; k < ARRAY_SIZE(x); k++)
  {
    t = x[k];
    integrand = cpow( 2 + t/zetatilde, -1.0 /6.0);
    integrand *= cexp(-I * t * tan(tau));
    sum += integrand * w[k];
  }

  return sum * az * eitaucostau16;
}                       
/*!\brief Compute scaled Airy function of first kind in region I ie 
    \f$|\arg z|< 2\pi/3\f$ 
*/
static complex cairyAie_mediumarg_regionI(const complex z)
{
  assert(cimag(z) >= 0);
  assert(isregionI(z));
  assert(ismedium(z));

  if(creal(z) > 0)
    return cairyAie_mediumarg_Q1(z);
  return cairyAie_mediumarg_Q2(z);
}

/*!\brief Compute scaled airy Ai function of first kind (\f$Aie\f$)
          in upper region II upper complex plane (\f$\frac{2\pi}{3} \le \arg z \le \pi\f$)

   This function use connexion formula given by
   [3] eq (7.1) for \f$\frac{2\pi}{3} \le \arg z \le \pi \f$

   The connexion formula is (see [3, pp 10] eq (7.1)) for 
   \f$\frac{2\pi}{3} \le \arg z \le \pi \f$:
   \f[
   Aie(z)=-e^{-\frac{2\pi i }{3}} e^{\frac{4}{3}z^\frac{3}{2}} 
                                   Aie\left(e^{-\frac{2\pi i }{3}}z\right)
		        -e^{\frac{2\pi i }{3}}
			        Aie\left(e^{\frac{2\pi i }{3}}z\right)
				   
   \f]

   \attention assert \f$\Im\text{m}\; z \ge 0\f$
   \note Because we assert \f$\Im\text{m}\; z \ge 0\f$ then
   \f$-\frac{\pi}{3} \le \arg e^{\frac{2\pi i }{3}} \le 0\f$
   we could therefore use the mirror property of scaled Airy function. 
   \note Because we assert \f$\Im\text{m}\; z \ge 0\f$  then
   \f$0 \le \arg e^{-\frac{2\pi i }{3}} \le \frac{\pi}{3}\f$
*/
static complex cairyAie_mediumarg_regionII(const complex z)
{
  complex exp_43z32;
  complex ze2pi_3, zem2pi_3;
  complex Aizem2pi_3, Aize2pi_3;
  complex first, second;

  assert(cimag(z) >= 0);
  assert(ismedium(z));
  assert(isregionII(z));
  
  exp_43z32 = cexp_43z32(z);
  
  ze2pi_3 = z * M_EXP_I_2PI_3;
  zem2pi_3 = z * M_EXP_MI_2PI_3;

  /* quick test */
  assert(cimag(ze2pi_3) <= 0);
  assert(cimag(zem2pi_3) >= 0);
  
  /* Use reflexion property */
  Aize2pi_3 = conj(cairyAie_mediumarg_regionI(conj(ze2pi_3)));
  Aizem2pi_3 = cairyAie_mediumarg_regionI(zem2pi_3);

  first = - exp_43z32 * M_EXP_MI_2PI_3 * Aizem2pi_3;
  second =  M_EXP_I_2PI_3 * Aize2pi_3;
  
  return first - second;  
}


/*!\brief Compute scaled airy Ai function of first kind (\f$Aie\f$)
          in upper complex plane (\f$\Im\text{m}\; z > 0\f$)

   This function use integral representation in 
   upper region I (\f$0 \le \arg z < \frac{2\pi}{3}\f$) and connexion formula given by
   [3] eq (7.1) for \f$\frac{2\pi}{3} \le \arg z \le \pi \f$

   \attention assert \f$\Im\text{m}\; z \ge 0\f$
*/
static complex cairyAie_mediumarg(const complex z)
{
  assert(cimag(z) >= 0);
  assert(ismedium(z));

  DEBUG("z= " fmtc "arg=%e (%e°) norm2(z)=%e cabs(z)=%e\n", 
	cprint(z), carg(z), carg(z) * 180 / M_PI, norm2(z), cabs(z));

  /* region I */
  if(!arglargerpi2_3(z))
  {
    DEBUG("z is in region I\n");
    return cairyAie_mediumarg_regionI(z);
  }

  /* Region II */
  DEBUG("z is in region II\n");
  return cairyAie_mediumarg_regionII(z);
}


/*!\brief Compute airy function of first kind (\f$Ai\f$)
          in upper complex plane (\f$\Im\text{m}\; z > 0\f$)

   \attention assert \f$\Im\text{m}\; z \ge 0\f$
*/
static complex cairyAi_mediumarg(const complex z)
{
  complex zeta;
  complex expmzeta, scaled;

  assert(ismedium(z));

  zeta = czetafast(z);

  scaled = cairyAie_mediumarg(z);
  expmzeta = cexp(-zeta);
  return scaled * expmzeta;
}


/*!\brief Computed \f$(-1)^k c_k\f$ coefficient of large argument expansion 

  According to (10.4.58) [1-p448] \f$c_k\f$ is defined by:
  \f{align*}
  c_1 &= 1 \\
  c_k &=  \frac{\Gamma(3k+\frac{1}{2})}{54^k k! \Gamma(k + \frac{1}{2})}
       =  \frac{(2k + 1)(2k + 3)\dotsm(6k-1)}{216^k k!}
  \f}

  \note Use a table of 37 coefficient as recommended by [2]
*/
static const double m1kck[] = 
  {
    /*-(3*5)/(216^1*1!)*/
    -0.6944444444444444444444444444444444444444e-1,
    /*+(5*7*9*11)/(216^2*2!)*/
    +0.3713348765432098765432098765432098765432e-1,
    /*-(7*9*11*13*15*17)/(216^3*3!)*/
    -0.3799305912780064014631915866483767718336e-1,
    /*+(9*11*13*15*17*19*21*23)/(216^4*4!)*/
    +0.576491904126697213331301122796321698928e-1,
    /*-(11*13*15*17*19*21*23*25*27*29)/(216^5*5!)*/
    -0.1160990640255154110181092538964814532563e0,
    /*+(13*15*17*19*21*23*25*27*29*...*35)/(216^6*6!)*/
    +0.2915913992307505114690938436983388351461e0,
    /*-(15*17*19*21*23*25*27*29*31*...*41)/(216^7*7!)*/
    -0.877666969510016916465506668433293676422e0,
    /*+(17*19*21*23*25*27*29*31*33*...*47)/(216^8*8!)*/
    +0.3079453030173166993362480862680011319529e1,
    /*-(19*21*23*25*27*29*31*33*35*...*53)/(216^9*9!)*/
    -0.1234157333234523870642339938330245277287e2,
    /*+(21*23*25*27*29*31*33*35*37*...*59)/(216^(10)*10!)*/
    +0.5562278536591708278103323749835619339993e2,
    /*-(23*25*27*29*31*33*35*37*39*...*65)/(216^(11)*11!)*/
    -0.2784650807776025672055514983345736197358e3,
    /*+(25*27*29*31*33*35*37*39*41*...*71)/(216^(12)*12!)*/
    +0.1533169432012795615968528330529591098476e4,
    /*-(27*29*31*33*35*37*39*41*43*...*77)/(216^(13)*13!)*/
    -0.9207206599726414698033224087507298680056e4,
    /*+(29*31*33*35*37*39*41*43*45*...*83)/(216^(14)*14!)*/
    +0.5989251356587906862599588327558071175112e5,
    /*-(31*33*35*37*39*41*43*45*47*...*89)/(216^(15)*15!)*/
    -0.419524875116551068662647089796081559627e6,
    /*+(33*35*37*39*41*43*45*47*49*...*95)/(216^(16)*16!)*/
    +0.3148257417866826378983145912575629412305e7,
    /*-(35*37*39*41*43*45*47*49*51*...*101)/(216^(17)*17!)*/
    -0.2519891987160236767557016381168581809833e8,
    /*+(37*39*41*43*45*47*49*51*53*...*107)/(216^(18)*18!)*/
    +0.2142880369636803195620823884016893528254e9,
    /*-(39*41*43*45*47*49*51*53*55*...*113)/(216^(19)*19!)*/
    -0.1929375549182493052665328054052344852887e10,
    /*+(41*43*45*47*49*51*53*55*57*...*119)/(216^(20)*20!)*/
    +0.1833576693789056765675348223590718007761e11,
    /*-(*43*45*47*49*51*53*55*57*59...*123)/(216^(21)*21!)*/
    -0.1834183035288325633653415468373651446256e12,
    /*+(*45*47*49*51*53*55*57*59*61*...*129)/(216^(22)*22!)*/
    +0.1926471158970446563579032395664926711576e13,
    /*-(*47*49*51*53*55*57*59*61*63*...*135)/(216^(23)*23!)*/
    -0.2119699938864764905493571816510303720509e14,
    /*+(*49*51*53*55*57*59*61*63*65*...*141)/(216^(24)*24!)*/
    +0.2438268268797160418199984201202274713689e15,
    /*-(*51*53*55*57*59*61*63*65*67*...*147)/(216^(25)*25!)*/
    -0.2926599219297925046400592148165285843848e16,
    /*+(*53*55*57*59*61*63*65*67*69*...*153)/(216^(26)*26!)*/
    +0.3659030701264312805075099317724813844832e17,
    /*-(*55*57*59*61*63*65*67*69*71*...*159)/(216^(27)*27!)*/
    -0.4757681020363067632401403572743318907192e18,
    /*+(*57*59*61*63*65*67*69*71*73*...*165)/(216^(28)*28!)*/
    +0.6424049357901937699484057869724498212931e19,
    /*-(*59*61*63*65*67*69*71*73*75*...*171)/(216^(29)*29!)*/
    -0.8995207427058378952098438694307239188288e20,
    /*+(*61*63*65*67*69*71*73*75*77*...*177)/(216^(30)*30!)*/
    +0.130451329931760981793742403749617716469e22,
    /*-(*63*65*67*69*71*73*75*77*79*...*183)/(216^(31)*31!)*/
    -0.1957062178658161503299043185284923492816e23,
    /*+(*65*67*69*71*73*75*77*79*81*...*189)/(216^(32)*32!)*/
    +0.3033871086594338299189753708716215815665e24,
    /*-(*67*69*71*73*75*77*79*81*83*...*195)/(216^(33)*33!)*/
    -0.4854832179436167359995522969659058986843e25,
    /*+(*69*71*73*75*77*79*81*83*85*...*201)/(216^(34)*34!)*/
    +0.8011464687609593661835749240413276384457e26,
    /*-(*71*73*75*77*79*81*83*85*87*...*207)/(216^(35)*35!)*/
    -0.1362107954526321589052986810339312804334e28,
    /*+(*73*75*77*79*81*83*85*87*89*...*213)/(216^(36)*36!)*/
    +0.2383951672727105666951726336845791873788e29,
    /*-(*75*77*79*81*83*85*87*89*91*...*219)/(216^(37)*37!)*/
    -0.4291560449285803546171319066670932465888e30,
  };



/*!\brief Constant \f$\frac{1}{2} \pi^{-\frac{1}{2}}\f$ */
static const double halfinversesquarerootpi = 
  0.282094791773878143474039725780386292922e0;

/*!\brief Compute scalled airy function of first kind (\f$e^{\zeta}Ai\f$) 
  in upper region I (\f$0 < \arg z < \frac{2\pi}{3}\f$)

  Compute the scalled airy function of first kind (\f$e^{\zeta}Ai\f$) where 
  \f$\zeta\f$ is \f$\frac{2}{3}z^\frac{3}{2}\f$

  For large argument according to (10.4.59) [1-p448] 
  \f$|z|\f$ \f$A_i\f$ could be approximated by (\f$|\arg z| < \pi\f$):
  \f[
  Ai = \frac{1}{2} \pi^\frac{-1}{2} z^{-1}{4} e^{-\zeta} \sum_0^\infty (-1)^k c_k \zeta^{-k}
  \f]
  and c_k defined by (10.4.58) [1-p448]:
  \f{align*}
  c_1 &= 1 \\
  c_k &=  \frac{\Gamma(3k+\frac{1}{2})}{54^k k! \Gamma(k + \frac{1}{2})}
       =  \frac{(2k + 1)(2k + 3)\dotsm(6k-1)}{216^k k!}
  \f}
  Therefore scaled \f$Ai\f$ is defined by:
  \f[
  Ai = \frac{1}{2} \pi^\frac{-1}{2} z^{-1}{4} \sum_0^\infty (-1)^k c_k \zeta^{-k}
  \f]

  According to [2, pp 475-476] this expansion is valid in region I, ie 
  (\f$|\arg z| < \frac{2\pi}{3}\f$) and \f$|z| > \rho\f$ where 
  \f$\rho \f$ is defined in [2, tab I pp 480].

  \attention assert \f$\Im\text{m}\; z \ge 0\f$
  \attention Computation are only valid for region I ie \f$|\arg z| < \frac{2\pi}{3}\f$
             (asserted)
*/
static complex cairyAie_largearg_regionI(const complex z)
{
  complex z1_4;
  complex z32;
  complex coeff;
  complex invzeta;
  complex zmk;

  complex sum;
  complex aik;

  unsigned int k;

  assert(islarge(z));
  assert(cimag(z) >= 0);
  assert(isregionI(z));

  /* z1_4 = z^(1/4) */
  z1_4 = cz1_4(z);
  /* coeff = 1/2 * sqrt(pi) * z^(-1/4) */
  coeff = halfinversesquarerootpi / z1_4;

  DEBUG("constant=%e coeff=" fmtc " z^1/4=" fmtc "\n", 
	halfinversesquarerootpi, cprint(coeff), cprint(z1_4));

  /* underflow */
  if(cabs(z) > DBL_MAX_POW_2_3)
    goto overflow;

  /* overflow is avoided because |z| < DBL_MAX^(2/3) */
  z32 = cpow32fast(z);
  invzeta = (3 / 2.0) / z32;

  sum = 1;
  zmk = 1;
  
  for ( k = 0 ; k < ARRAY_SIZE(m1kck) ; k++)
    {
      zmk *= invzeta;
      aik = zmk * m1kck[k];
      
      if (norm2 (aik) < norm2 (sum) * sqr (DBL_EPSILON)
	  || norm2 (aik) < DBL_MIN)
	return coeff * sum;
      
      sum += aik;
      DEBUG("zmk=" fmtc " aik=" fmtc " sum=" fmtc "\n",
	    cprint(zmk), cprint(aik), cprint(sum));
    }
  /* Do not converge */
  assert (0 && sizeof ("large expansion do not converge!"));
  errno = ERANGE;
  (void) feraiseexcept (FE_INVALID);
  return coeff * sum;

 overflow:
  return coeff;
}


/*!\brief Compute scaled airy Ai function of first kind (\f$e^{\zeta}Ai\f$)
          in upper region II upper complex plane (\f$\frac{2\pi}{3} \le \arg z \le \pi\f$)

   This function use connexion formula given by
   [3] eq (7.1) for \f$\frac{2\pi}{3} \le \arg z \le \pi \f$

   The connexion formula is (see [3, pp 10] eq (7.1)) for 
   \f$\frac{2\pi}{3} \le \arg z \le \pi \f$:
   \f[
   Ai_\text{scalled}(z)=-e^{-\frac{2\pi i }{3}} e^{\frac{4}{3}z^\frac{3}{2}} 
                                   Ai_\text{scalled}\left(e^{-\frac{2\pi i }{3}}z\right)
		        -e^{\frac{2\pi i }{3}}
			        Ai_\text{scalled}\left(e^{\frac{2\pi i }{3}}z\right)
				   
   \f]

   \attention assert \f$\Im\text{m}\; z \ge 0\f$
   \note Because we assert \f$\Im\text{m}\; z \ge 0\f$ then
   \f$-\frac{\pi}{3} \le \arg e^{\frac{2\pi i }{3}} \le 0\f$
   we could therefore use the mirror property of scaled Airy function. 
   \note Because we assert \f$\Im\text{m}\; z \ge 0\f$  then
   \f$0 \le \arg e^{-\frac{2\pi i }{3}} \le \frac{\pi}{3}\f$
*/
static complex cairyAie_largearg_regionII(const complex z)
{
  complex exp_43z32;
  complex ze2pi_3, zem2pi_3;
  complex Aizem2pi_3, Aize2pi_3;
  complex first, second;

  assert(cimag(z) >= 0);
  assert(islarge(z));
  assert(isregionII(z));

  DEBUG("z= " efmtc " arg=%e (%e°) norm2(z)=%e cabs(z)=%e\n", 
	cprint(z), carg(z), carg(z) * 180 / M_PI, norm2(z), cabs(z));


  exp_43z32 = cexp_43z32(z);

  DEBUG("exp(3/4 z^(3/2))=" fmtc "\n", cprint(exp_z32));
  
  ze2pi_3 = z * M_EXP_I_2PI_3;
  zem2pi_3 = z * M_EXP_MI_2PI_3;

  DEBUG("z e^(2pi/3)= " fmtc " z e^(-2pi/3)=" fmtc "\n",
	cprint(ze2pi_3), cprint(zem2pi_3));

  /* quick test */
  assert(cimag(ze2pi_3) <= 0);
  assert(cimag(zem2pi_3) >= 0);

  /* Use reflexion property */
  Aize2pi_3 = conj(cairyAie_largearg_regionI(conj(ze2pi_3)));
  Aizem2pi_3 = cairyAie_largearg_regionI(zem2pi_3);

  DEBUG("Ai(z e^(2pi/3))= " fmtc " Ai(z e^(-2pi/3))=" fmtc "\n",
	cprint(Aize2pi_3), cprint(Aizem2pi_3));
  
  first = - exp_43z32 * M_EXP_MI_2PI_3 * Aizem2pi_3;
  second =  M_EXP_I_2PI_3 * Aize2pi_3;
  
  DEBUG("first=" fmtc " second=" fmtc " total=" fmtc "\n", 
	cprint(first), cprint(second), cprint(first - second));
  
  return first - second;  
}

/*!\brief Compute scaled airy Ai function of first kind (\f$e^{\zeta}Ai\f$)
          in upper complex plane (\f$\Im\text{m}\; z > 0\f$)

   This function use the asymtotic expansion in 
   upper region I (\f$0 \le \arg z < \frac{2\pi}{3}\f$) and connexion formula given by
   [3] eq (7.1) for \f$\frac{2\pi}{3} \le \arg z \le \pi \f$

   \attention assert \f$\Im\text{m}\; z \ge 0\f$
*/
static complex cairyAie_largearg(const complex z)
{
  assert(cimag(z) >= 0);
  assert(islarge(z));

  DEBUG("z= " fmtc "arg=%e (%e°) norm2(z)=%e cabs(z)=%e\n", 
	cprint(z), carg(z), carg(z) * 180 / M_PI, norm2(z), cabs(z));

  /* region I */
  if(!arglargerpi2_3(z))
  {
    DEBUG("z is in region I\n");
    return cairyAie_largearg_regionI(z);
  }

  /* Region II */
  DEBUG("z is in region II\n");
  return cairyAie_largearg_regionII(z);
}

/*!\brief Compute airy function of first kind (\f$Ai\f$)
          in upper complex plane (\f$\Im\text{m}\; z > 0\f$)

   \attention assert \f$\Im\text{m}\; z \ge 0\f$
*/
static complex cairyAi_largearg(const complex z)
{
  complex zeta;
  complex expmzeta, scaled;

  //  assert(islarge(z));

  zeta = czetafast(z);

  scaled = cairyAie_largearg(z);
  expmzeta = cexp(-zeta);
  return scaled * expmzeta;
}



/*!\brief Compute scaled airy function of second kind (\f$Bie(z)\f$) 
  in upper region A (\f$0 \le \arg z < \frac{\pi}{3}\f$)

  Use the connection formula given in [5, pp3] eq (2.5), in region A:
  \f[
  Bie(z)=i e^{-2\zeta} Aie(z) 
           + 2 e^{-i\frac{\pi}{6}} Aie(e^{-\frac{2\pi}{3}}z)
  \f]
  Where \f$\zeta\f$ is \f$\frac{2}{3}z^\frac{3}{2}\f$. The previous equation could 
  be simplified using reflexion properties:
  \f[
  Bie(z)=i \overline{e^{\frac{4}{3} z^\frac{3}{2}}} Aie(z) 
              + 2 e^{-i\frac{\pi}{6}} \overline{Aie(\overline{e^{-2\frac{\pi}{3}}z})}
  \f]
  \attention assert z in upper region A
 
*/
static complex cairyBie_regionA(const complex z)
{
  complex exp_m43z32;
  complex Aiz, Aiem2pi3z;
  complex first, second;

  assert (cimag(z) >= 0);
  assert (isregionA(z));

  exp_m43z32 = cexp(-(4.0/3.0) * cpow32(z));

  Aiz = cairyAie(z);
  Aiem2pi3z = conj(cairyAie(conj(M_EXP_MI_2PI_3*z)));

  first = I * exp_m43z32 * Aiz;
  second = 2 * M_EXP_MI_PI_6 * Aiem2pi3z;

  return first + second;
}


/*!\brief Compute scalled airy function of second kind (\f$Bie(z)\f$) 
  in upper region B (\f$\frac{\pi}{3} \le \arg z < \frac{2\pi}{3}\f$)

  Use the connection formula given in [5, pp3] eq (2.5), in region B:
  \f[
  Bie(z)=i Aie(z) + 2 e^{-i\frac{\pi}{6}} 
                      e^{2\zeta} Aie(e^{-\frac{2\pi}{3}}z)
  \f]
  Where \f$\zeta\f$ is \f$\frac{2}{3}z^\frac{3}{2}\f$. The previous equation could 
  be simplified using reflexion properties:
  \f[
  Bie(z)=i Aie(z) + 2 e^{-i\frac{\pi}{6}} 
                      e^{\frac{4}{3} z^\frac{3}{2}} 
		      \overline{Aie(\overline{e^{-2\frac{\pi}{3}}z})}
  \f]
  \attention assert z in upper region B
  \todo Use Aie Bie always
*/
static complex cairyBie_regionB(const complex z)
{
  complex exp_43z32;
  complex Aiz, Aiem2pi3z;
  complex first, second;


  assert (cimag(z) >= 0);
  assert (isregionB(z));

  exp_43z32 = cexp_43z32(z);

  Aiz = cairyAie(z);
  Aiem2pi3z = conj(cairyAie(conj(M_EXP_MI_2PI_3*z)));

  first = I * Aiz;
  second = 2 * M_EXP_MI_PI_6 * exp_43z32 * Aiem2pi3z;

  return first + second;
}


/*!\brief Compute scalled airy function of second kind (\f$Bie(z)\f$) 
  in upper region C (\f$\frac{2\pi}{3} \le \arg z \le \pi\f$)

  Use the connection formula given in [5, pp3] eq (2.5), in region C:
  \f[
  Bie(z)=e^{i\frac{\pi}{6}} Aie(z e^{i\frac{2\pi}{3}}) + 2 e^{-i\frac{\pi}{6}} 
                      e^{2\zeta} Aie(e^{-\frac{2\pi}{3}}z)
  \f]
  Where \f$\zeta\f$ is \f$\frac{2}{3}z^\frac{3}{2}\f$. The previous equation could 
  be simplified using reflexion properties:
  \f[
  Bie(z)=e^{i\frac{\pi}{6}} \overline{Aie(\overline{z e^{i\frac{2\pi}{3}}})} + e^{-i\frac{\pi}{6}} 
                      e^{\frac{4}{3} z^\frac{3}{2}} Aie(e^{-2\frac{\pi}{3}}z)
  \f]
  \attention assert z in upper region C
*/
static complex cairyBie_regionC(const complex z)
{
  complex exp_43z32;
  complex Aie2pi3z, Aiem2pi3z;
  complex z2pi3, zm2pi3;
  complex first, second;

  assert (cimag(z) >= 0);
  assert (isregionC(z));

  exp_43z32 = cexp_43z32(z);

  z2pi3 = M_EXP_I_2PI_3*z;
  zm2pi3 = M_EXP_MI_2PI_3*z;
  
  /* quick test */
  assert(cimag(z2pi3) <= 0);
  assert(cimag(zm2pi3) >= 0);

  Aie2pi3z = conj(cairyAie(conj(z2pi3)));
  Aiem2pi3z = cairyAie(zm2pi3);

  first =  M_EXP_I_PI_6 * Aie2pi3z;
  second = M_EXP_MI_PI_6 * exp_43z32 * Aiem2pi3z;

  return first + second;
}

/*!\brief Compute scalled airy function of second kind (\f$Bie(z)\f$)
  for large argument or medium argument

  \f$Bie(z)\f$ is defined as ([5, pp3] eq (2.4)):
  \f[
  Bie(z) = \begin{cases}
  e^{-\zeta} Bi(z) & \text{$|\arg z| \le \frac{\pi}{3}$} \\
  e^{\zeta} Bi(z) & \text{else}
  \end{cases}
  \f]
  \note assert z in upper plane
*/
static complex cairyBie_largemediumarg(const complex z)
{
  assert(cimag(z) > 0);
  assert(islarge(z) || ismedium(z));

  if(!arglargerpi_3(z))
    return cairyBie_regionA(z);
  if(!arglargerpi2_3(z))
    return cairyBie_regionB(z);
  return cairyBie_regionC(z);
}


/*!\brief Compute airy function of second kind (\f$e^{\zeta}Bi\f$) 
  for large argument using scaled version
  
*/
static complex cairyBi_largemediumarg(const complex z)
{
  complex zeta;
  complex exppmzeta;
  complex scalled;

  assert(islarge(z) || ismedium(z));
  assert (cimag(z) >= 0);

  zeta = czeta(z);
  
  scalled = cairyBie_largemediumarg(z);
  if(!arglargerpi_3(z))
    exppmzeta = cexp(zeta);
  else
    exppmzeta = cexp(-zeta);
  
  DEBUG("zeta=" fmtc  " scalled=" fmtc " exp(-zeta)=" fmtc "\n",
	cprint(zeta), cprint(scalled), cprint(exppmzeta));
  return scalled * exppmzeta;
}

/*!\brief Test in order to use serie representation */
static inline bool
use_small (const complex z)
{
  double x, y;

  x = creal (z);
  y = cimag (z);

  return fabs (y) <= MAX_IMAG_SMALL &&
    x >= MIN_REAL_SMALL && x <= MAX_REAL_SMALL;
}


/*!\brief Compute the Airy function of first kind (\f$Ai\f$)
  
  Airy function of first kind [4] is one of the two 
  independant solution (with \f$Bi\f$) of the following 
  differential equation:
  \f[
  y''-yz=0.
  \f]
  \f$Ai\f$ could be written:
  \f[
  Ai(z)=\frac{1}{3^{2/3}\Gamma(\frac{2}{3})}\leftidx{_0}{F}{_1}(;\frac{2}{3};\frac{1}{9}z^3)
         -\frac{z}{3^(1/3)\Gamma(\frac{1}{3})}\leftidx{_0}{F}{_1}(;\frac{4}{3};\frac{1}{9} z^3)
  \f]
  where \f$\leftidx{_0}{F_1}(;a;z)\f$ is a confluent hypergeometric limit function.

  Another equivalent definition is:
  \f[
  Ai(z) = \frac{1}{2\pi i} \int_{\infty e^{-i\frac{\pi}{3}}}^{\infty e^{i\frac{\pi}{3}}}
              \exp(\frac{1}{3} t^3-zt) \;\text{d}t 
  \f]
*/
complex
cairyAi (const complex z)
{
  /* mirror symetry */
  if(cimag(z) < 0.0)
    return conj(cairyAi(conj(z)));

  /* after this point Im z > 0 */

  /* small */
  if (use_small (z))
    return cairyAi_smallarg (z);
  
  /* large */
  if (norm2(z) > sqr(AIRY_LARGE))
    return cairyAi_largearg(z);

  /* medium */
  return cairyAi_mediumarg(z);
}

/*!\brief Compute the scales Airy function of first kind (\f$e^\zeta Ai(z)\f$)
  
  Compute the scalled airy function of first kind (\f$e^{\zeta}Ai\f$) where 
  \f$\zeta\f$ is \f$\frac{2}{3}z^\frac{3}{2}\f$	  

  Airy function of first kind [4] is one of the two 
  independant solution (with \f$Bi\f$) of the following 
  differential equation:
  \f[
  y''-yz=0.
  \f]
  \f$Ai\f$ could be written:
  \f[
  Ai(z)=\frac{1}{3^{2/3}\Gamma(\frac{2}{3})}\leftidx{_0}{F}{_1}(;\frac{2}{3};\frac{1}{9}z^3)
         -\frac{z}{3^(1/3)\Gamma(\frac{1}{3})}\leftidx{_0}{F}{_1}(;\frac{4}{3};\frac{1}{9} z^3)
  \f]
  where \f$\leftidx{_0}{F_1}(;a;z)\f$ is a confluent hypergeometric limit function.
  
  Another equivalent definition is:
  \f[
  Ai(z) = \frac{1}{2\pi i} \int_{\infty e^{-i\frac{\pi}{3}}}^{\infty e^{i\frac{\pi}{3}}}
              \exp(\frac{1}{3}t^3-zt) \;\text{d}t 
  \f]
*/
complex cairyAie(const complex z)
{
 /* mirror symetry */
  if(cimag(z) < 0.0)
    return conj(cairyAie(conj(z)));

  /* after this point Im z > 0 */
  /* small */
  if (use_small (z))
    return cairyAie_smallarg(z);
  
  /* large */
  if (norm2(z) > sqr(AIRY_LARGE))
    return cairyAie_largearg(z);

  /* medium */
  return cairyAie_mediumarg(z);
}

/*!\brief Compute the Airy function of second kind (\f$Bi\f$)
  
  Airy function of second kind [4] is one of the two 
  independant solution (with \f$Ai\f$) of the following 
  differential equation:
  \f[
  y''-yz=0.
  \f]
  \f$Bi\f$ could be written:
  \f[
  Bi(z)=\frac{1}{3^{1/6}\Gamma(\frac{2}{3})}\leftidx{_0}{F}{_1}(;\frac{2}{3};\frac{1}{9}z^3)
         + \frac{3^{1/6} z}{\Gamma(\frac{1}{3})}\leftidx{_0}{F}{_1}(;\frac{4}{3};\frac{1}{9} z^3)
  \f]
  where \f$\leftidx{_0}{F_1}(;a;z)\f$ is a confluent hypergeometric limit function.

  Another equivalent definition is:
  \f[
  Bi(z) = -\frac{1}{2\pi i} \int_{\infty e^{-i\frac{\pi}{3}}}^{\infty e^{i\frac{\pi}{3}}}
              t\exp(\frac{1}{3} t^3-zt) \;\text{d}t 
  \f]
*/
complex
cairyBi (const complex z)
{
  /* mirror symetry */
  if(cimag(z) < 0.0)
    return conj(cairyBi(conj(z)));

  /* after this point Im z > 0 */

  /* small */
  if (use_small (z))
    return cairyBi_smallarg (z);
  
  return cairyBi_largemediumarg(z);
}

/*!\brief Compute the scaled Airy function of second kind (\f$Bie(z)\f$)
  
  \f$Bie(z)\f$ is defined as ([5, pp3] eq (2.4)):
  \f[
  Bie(z) = \begin{cases}
  e^{-\zeta} Bi(z) & \text{$|\arg z| \le \frac{\pi}{3}$} \\
  e^{\zeta} Bi(z) & \text{else}
  \end{cases}
  \f]
  

  Airy function of second kind [4] is one of the two 
  independant solution (with \f$Ai\f$) of the following 
  differential equation:
  \f[
  y''-yz=0.
  \f]
  \f$Bi\f$ could be written:
  \f[
  Bi(z)=\frac{1}{3^{1/6}\Gamma(\frac{2}{3})}\leftidx{_0}{F}{_1}(;\frac{2}{3};\frac{1}{9}z^3)
         + \frac{3^{1/6} z}{\Gamma(\frac{1}{3})}\leftidx{_0}{F}{_1}(;\frac{4}{3};\frac{1}{9} z^3)
  \f]
  where \f$\leftidx{_0}{F_1}(;a;z)\f$ is a confluent hypergeometric limit function.

  Another equivalent definition is:
  \f[
  Bi(z) = -\frac{1}{2\pi i} \int_{\infty e^{-i\frac{\pi}{3}}}^{\infty e^{i\frac{\pi}{3}}}
              t\exp(\frac{1}{3} t^3-zt) \;\text{d}t 
  \f]
*/
complex cairyBie(const complex z)
{
 /* mirror symetry */
  if(cimag(z) < 0.0)
    return conj(cairyBie(conj(z)));

  /* after this point Im z > 0 */
  /* small */
  if (use_small (z))
    return cairyBie_smallarg(z);
  
  return cairyBie_largemediumarg(z);
}

/*!\brief Coefficient used in auxiliary function \f$f'\f$ 
  
  Auxiliary function could be computed from is defined as using [6]:
  is:
  \f[
  f'(z)=1+\frac{2}{3!} z^3 + \frac{2\times5}{6!} z^6 + \frac{2\times4\times7}{9!} z^9 + \dotsb
  \f]
  Therefore coefficient are \f$\frac{2}{3!}\f$, \f$\frac{2\times5}{6!}\f$, etc
*/
static const double fpcoeff[] = 
{
	/*(2)/3! */
	0.3333333333333333333333333333333333333333e0,
	/*(2*5)/6! */
	0.1388888888888888888888888888888888888889e-1,
	/*(2*5*8)/9! */
	0.2204585537918871252204585537918871252205e-3,
	/*(2*5*8*11)/12! */
	0.1837154614932392710170487948265726043504e-5,
	/*(2*5*8*11*14)/15! */
	0.9421305717602013898310194606490902787199e-8,
	/*(2*5*8*11*14*17)/18! */
	0.3271286707500699270246595349476007912222e-10,
	/*(2*5*8*11*14*17*20)/21! */
	0.8198713552633331504377431953573954667223e-13,
	/*(2*5*8*11*14*17*20*23)/24! */
	0.1552786657695706724313907566964764141519e-15,
	/*(2*5*8*11*14*17*20*23*26)/27! */
	0.2300424678067713665650233432540391320769e-18,
	/*(2*5*8*11*14*17*20*23*26*29)/30! */
	0.2738600807223468649583611229214751572345e-21,
	/*(2*5*8*11*14*17*20*23*26*29*32)/33! */
	0.2677029137070839344656511465508066053123e-24,
	/*(2*5*8*11*14*17*20*23*26*29*32*35)/36! */
	0.2187115308064411229294535511036001677388e-27,
	/*(2*5*8*11*14*17*20*23*26*29*32*35*38)/39! */
	0.151567242416106114296225607140402056645e-30,
	/*(2*5*8*11*14*17*20*23*26*29*32*...*41)/42! */
	0.9021859667625363946203905186928693847919e-34,
	/*(2*5*8*11*14*17*20*23*26*29*32*...*44)/45! */
	0.4662459776550575682792715858877877957581e-37,
	/*(2*5*8*11*14*17*20*23*26*29*32*...*47)/48! */
	0.2111621275611673769380758994057009944557e-40,
	/*(2*5*8*11*14*17*20*23*26*29*32*...*50)/51! */
	0.8449865048466081510127086810952420746525e-44,
	/*(2*5*8*11*14*17*20*23*26*29*32*...*53)/54! */
	0.3009211199596182873976882767433198271555e-47,
	/*(2*5*8*11*14*17*20*23*26*29*32*...*56)/57! */
	0.9598759807324347285412704202338750467479e-51,
	/*(2*5*8*11*14*17*20*23*26*29*32*...*59)/60! */
	0.2758264312449525082015144885729525996402e-54,
	/*(2*5*8*11*14*17*20*23*26*29*32*...*62)/63! */
	0.7177372657948282805139591167654244070783e-58,
	/*(2*5*8*11*14*17*20*23*26*29*32*...*65)/66! */
	0.1699188602733968467125850181736326721303e-61,
	/*(2*5*8*11*14*17*20*23*26*29*32*...*68)/69! */
	0.3675510713246741222422345190863782654776e-65,
	/*(2*5*8*11*14*17*20*23*26*29*32*...*71)/72! */
	0.7292679986600677028615764267586870346777e-69,
	/*(2*5*8*11*14*17*20*23*26*29*32*...*74)/75! */
	0.1331996344584598544039408998646003716306e-72,
	/*(2*5*8*11*14*17*20*23*26*29*32*...*77)/78! */
	0.2246957396397770823278355260873825432365e-76,
	/*(2*5*8*11*14*17*20*23*26*29*32*...*80)/81! */
	0.3511419591182639198747234350482615146688e-80,
	/*(2*5*8*11*14*17*20*23*26*29*32*...*83)/84! */
	0.5097879778139720091096449405462565543972e-84,
	/*(2*5*8*11*14*17*20*23*26*29*32*...*86)/87! */
	0.6893684622230858811489451528685010877582e-88,
	/*(2*5*8*11*14*17*20*23*26*29*32*...*89)/90! */
	0.8704147250291488398345267081672993532301e-92,
	/*(2*5*8*11*14*17*20*23*26*29*32*...*92)/93! */
	0.1028494298746483327229737336839536043046e-95,
	/*(2*5*8*11*14*17*20*23*26*29*32*...*95)/96! */
	0.1139732157298851204820187651639556785291e-99,
	/*(2*5*8*11*14*17*20*23*26*29*32*...*98)/99! */
	0.1186850106528013334187428565697757768708e-103,
	/*(2*5*8*11*14*17*20*23*26*29*32*...*101)/102! */
	0.1163578535811777778615126044801723302655e-107,
	/*(2*5*8*11*14*17*20*23*26*29*32*...*104)/105! */
	0.1075893236996558278885923296164330376935e-111
};

/*!\brief Auxiliary function \f$f'\f$ 
  
  Auxiliary \f$f'\f$ function could be computed from is defined as using [6]:
  is:
  \f[
  f'(z)=1+\frac{2}{3!} z^3 + \frac{2\times5}{6!} z^6 + \frac{2\times4\times7}{9!} z^9 + \dotsb
  \f]

  \param[in] z complex number
  \return f(z)
*/
static complex
fp (const complex z)
{
  unsigned int k;

  complex sum;
  complex fpk;

  complex z3;
  complex z3n = 1;

  z3 = z * z * z;

  /* f_0 = 1 */
  fpk = 1;
  sum = 1;

  for (k = 0; k < ARRAY_SIZE (fpcoeff); k++)
    {
      z3n *= z3;
      fpk = z3n * fpcoeff[k];

      if (norm2 (fpk) < norm2 (sum) * sqr (DBL_EPSILON)
	  || norm2 (fpk) < DBL_MIN)
	return sum;

      sum += fpk;
      DEBUG ("fpcoeff=%.8e fpk=" fmtc " sum=" fmtc "\n",
	     fpcoeff[k], creal (fpk), cimag (fpk), creal (sum), cimag (sum));
    }

  /* Do not converge */
  assert (0 && sizeof ("f do not converge!"));
  errno = ERANGE;
  (void) feraiseexcept (FE_INVALID);
  return sum;
}

/*!\brief Coefficient used in auxiliary function \f$f'\f$ 
  
  Auxiliary function could be computed from is defined as using [6]:
  is:
  \f[
  g'(z)=\frac{1}{2!}z^2 +\frac{1\times4}{5!} z^5 + 
       \frac{1\times4\times7}{8!} z^8 + \dotsb
  \f]
  Therefore coefficient are \f$\frac{1\times4}{5!}\f$, \f$\frac{1\times4\times7}{8!}\f$, etc
*/
static const double gpcoeff[] = 
{
	/*(1*4)/5! */
	0.3333333333333333333333333333333333333333e-1,
	/*(1*4*7)/8! */
	0.6944444444444444444444444444444444444444e-3,
	/*(1*4*7*10)/11! */
	0.7014590347923681257014590347923681257015e-5,
	/*(1*4*7*10*13)/14! */
	0.4175351397573619795842018064240286462509e-7,
	/*(1*4*7*10*13*16)/17! */
	0.1637392704930831292487065907545210377454e-9,
	/*(1*4*7*10*13*16*19)/20! */
	0.4548313069252309145797405298736695492929e-12,
	/*(1*4*7*10*13*16*19*22)/23! */
	0.9416797244828797403307257347280943049542e-15,
	/*(1*4*7*10*13*16*19*22*25)/26! */
	0.1509102122568717532581291241551433181016e-17,
	/*(1*4*7*10*13*16*19*22*25*28)/29! */
	0.1927333489870648189758992645659557063878e-20,
	/*(1*4*7*10*13*16*19*22*25*28*31)/32! */
	0.2007639051948591864332284005895371941539e-23,
	/*(1*4*7*10*13*16*19*22*25*28*31*34)/35! */
	0.1738215629392720228859120351424564451549e-26,
	/*(1*4*7*10*13*16*19*22*25*28*31*34*37)/38! */
	0.1270625460082397828113392069754798575694e-29,
	/*(1*4*7*10*13*16*19*22*25*28*31*...*40)/41! */
	0.794637561027140605449275840997372467601e-33,
	/*(1*4*7*10*13*16*19*22*25*28*31*...*43)/44! */
	0.4299986802094916696154090048687080452386e-36,
	/*(1*4*7*10*13*16*19*22*25*28*31*...*46)/47! */
	0.2033090686569700565557489384722023854556e-39,
	/*(1*4*7*10*13*16*19*22*25*28*31*...*49)/50! */
	0.8471211194040419023156205769675099393984e-43,
	/*(1*4*7*10*13*16*19*22*25*28*31*...*52)/53! */
	0.313400340142079875070521856073810558416e-46,
	/*(1*4*7*10*13*16*19*22*25*28*31*...*55)/56! */
	0.1036376786184126571000402963207045497407e-49,
	/*(1*4*7*10*13*16*19*22*25*28*31*...*58)/59! */
	0.3081703200071741216177231529012921490952e-53,
	/*(1*4*7*10*13*16*19*22*25*28*31*...*61)/62! */
	0.8284148387289626925207611637131509384278e-57,
	/*(1*4*7*10*13*16*19*22*25*28*31*...*64)/65! */
	0.202299105916718606232176108354859814024e-60,
	/*(1*4*7*10*13*16*19*22*25*28*31*...*67)/68! */
	0.4507555835934015290378255533753560918537e-64,
	/*(1*4*7*10*13*16*19*22*25*28*31*...*70)/71! */
	0.9200971291965738498424689801497368684501e-68,
	/*(1*4*7*10*13*16*19*22*25*28*31*...*73)/74! */
	0.1726909026269845814268898236016773401746e-71,
	/*(1*4*7*10*13*16*19*22*25*28*31*...*76)/77! */
	0.2990318660207525219513243698730343552806e-75,
	/*(1*4*7*10*13*16*19*22*25*28*31*...*79)/80! */
	0.4792177340076162210758403363349909539754e-79,
	/*(1*4*7*10*13*16*19*22*25*28*31*...*82)/83! */
	0.7128034121785158724912097818458886716873e-83,
	/*(1*4*7*10*13*16*19*22*25*28*31*...*85)/86! */
	0.9867156868473364790852848585906543074299e-87,
	/*(1*4*7*10*13*16*19*22*25*28*31*...*88)/89! */
	0.1274332541453359781848488775139680107749e-90,
	/*(1*4*7*10*13*16*19*22*25*28*31*...*91)/92! */
	0.1539048963107922441846000936159033946556e-94,
	/*(1*4*7*10*13*16*19*22*25*28*31*...*94)/95! */
	0.1741990903347959753079797324458442497517e-98,
	/*(1*4*7*10*13*16*19*22*25*28*31*...*97)/98! */
	0.1851605977198086472236179128888650613858e-102,
	/*(1*4*7*10*13*16*19*22*25*28*31*...*100)/101! */
	0.1851791156313717844020581187007351348993e-106,
	/*(1*4*7*10*13*16*19*22*25*28*31*...*103)/104! */
	0.1745655313267079415554846518672088375747e-110
};

/*!\brief Auxiliary function \f$g'\f$ 
  
  Auxiliary \f$f'\f$ function could be computed from is defined as using [6]:
  is:
  \f[
  f'(z)=\frac{1}{2!}z^2 +\frac{1\times4}{5!} z^5 + 
       \frac{1\times4\times7}{8!} z^8 + \dotsb
  \f]

  \param[in] z complex number
  \return f(z)
*/
static complex
gp (const complex z)
{
  unsigned int k;

  complex sum;
  complex gpk;

  complex z3;
  complex z3np1;

  z3 = z * z * z;

  sum = 0.5 * z * z;
  z3np1 = z * z;

  for (k = 0; k < ARRAY_SIZE (gcoeff); k++)
    {
      z3np1 *= z3;
      gpk = z3np1 * gpcoeff[k];

      if (norm2 (gpk) < norm2 (sum) * sqr (DBL_EPSILON)
	  || norm2 (gpk) < DBL_MIN)
	return sum;

      sum += gpk;
      DEBUG ("gpcoeff=%.8e gpk=" fmtc " sum=" fmtc "\n",
	     gpcoeff[k], creal (gpk), cimag (gpk), creal (sum), cimag (sum));
    }

  /* Do not converge */
  assert (0 && sizeof ("g do not converge!"));
  errno = ERANGE;
  (void) feraiseexcept (FE_INVALID);
  return sum;
}

/*!\brief Compute derivative of airy function of first kind  \f$Ai'\f$ 
          for small argument 

	  Compute derivative of airy function of first kind  \f$Ai'\f$ 
          for small argument using the relation [6]
	  \f[
	  Ai'=Ai'(0) f'(z) + Ai(0) g'(z)
	  \f]
	  \param[in] z: Where to evaluate \f$Ai'\f$
	  \return \f$Ai(z)\f$
*/
static complex
cairyAip_smallarg (const complex z)
{
  assert (issmall (z));

  /* underflow */
  if (norm2 (z) < sqr (DBL_MIN))
    {
      DEBUG ("Underflow\n");
      errno = ERANGE;
      (void) feraiseexcept (FE_UNDERFLOW);
      return aip0;
    }

  return aip0 * fp (z) + ai0 * gp (z);
}


/*!\brief Compute scaled derivative airy function of first kind  \f$Aie'\f$ 
          for small argument

  Compute the scalled airy function of first kind \f$Aie'(z)=e^{\zeta}Ai'(z)\f$ where 
  \f$\zeta\f$ is \f$\frac{2}{3}z^\frac{3}{2}\f$	  

  \note Use non scaled version
*/
static complex
cairyAiep_smallarg (const complex z)
{
  complex zeta;
  complex expzeta, nonscaled;

  zeta = czetafast(z);

  nonscaled = cairyAip_smallarg(z);
  expzeta = cexp(zeta);
  DEBUG("zeta=" fmtc  " non scaled=" fmtc " exp(zeta)=" fmtc "\n",
	cprint(zeta), cprint(nonscaled), cprint(expmzeta));
  return nonscaled * expzeta;
}


/*!\brief Compute derivative of airy function of first kind  \f$Ai'\f$ 
          for small argument 

	  Compute derivative of airy function of first kind  \f$Ai'\f$ 
          for small argument using the relation [6]
	  \f[
	  Bi'=Bi'(0) f'(z) + Bi(0) g'(z)
	  \f]
	  \param[in] z: Where to evaluate \f$Ai'\f$
	  \return \f$Ai(z)\f$
	  \todo Check why not working
	  \bug not working
*/
static complex
cairyBip_smallarg (const complex z)
{
  assert (issmall (z));

  /* underflow */
  if (norm2 (z) < sqr (DBL_MIN))
    {
      DEBUG ("Underflow\n");
      errno = ERANGE;
      (void) feraiseexcept (FE_UNDERFLOW);
      return bip0;
    }

  return bip0 * fp (z) + bi0 * gp (z);
}

/*!\brief Compute scaled derivative airy function of first kind  \f$Bie'\f$ 
          for small argument

  \note Use non scaled version
*/
static complex
cairyBiep_smallarg (const complex z)
{
  complex zeta;
  complex exppmzeta, nonscaled;

  zeta = czetafast(z);

  nonscaled = cairyBip_smallarg(z);
  if(arglargerpi_3(z))
    exppmzeta = cexp(-zeta);
  else
    exppmzeta = cexp(zeta);

  DEBUG("zeta=" fmtc  " non scaled=" fmtc " exp(zeta)=" fmtc "\n",
	cprint(zeta), cprint(nonscaled), cprint(exppmzeta));
  return nonscaled * exppmzeta;
}

static double xp[25] = { .0435079659953445, .205779160144678,
      .489916161318751, .896390483211727, 1.4258249673758,
      2.07903190767599, 2.85702335104978, 3.76102058198275,
      4.79246521225895, 5.95303247470003, 7.24464710774066,
      8.66950223642504, 10.2300817341775, 11.9291866622602,
      13.7699665302828, 15.7559563095946, 17.8911203751898,
      20.1799048700978, 22.6273004064466, 25.2389175786164,
      28.0210785229929, 30.9809287996116, 34.1265753192057,
      37.4672580871163, 41.0135664833476
    };
static double wp[25] = { .0576354557898966, .139560003272262,
      .187792315011311, .187446935256946, .150716717316301,
      .10106990445338, .0575274105486025, .0280625783448681,
      .0117972164134041, .00428701743297432, .00134857915232883,
      3.67337337105948e-4, 8.65882267841931e-5, 1.76391622890609e-5,
      3.09929190938078e-6, 4.68479653648208e-7, 6.07273267228907e-8,
      6.72514812555074e-9, 6.33469931761606e-10, 5.04938861248542e-11,
      3.38602527895834e-12, 1.89738532450555e-13, 8.81618802142698e-15,
      3.36676636121976e-16, 1.04594827170761e-17
    };

/*!\brief coefficient \f$b(z)\f$ 

  Defined in [5, pp 10] eq (6.3) as :
  \f[
  bzcoeff = - \frac{1}{\sqrt{3\pi}2^{2/3}\Gamma(\frac{7}{6})}
  \f]
*/
static const double bzcoeff = - 0.2211877980269592949839044600746471027181;                       


/*!\brief Compute derivative of scaled airy function (\f$Aie'\f$) for medium argument 
  and in first quadrant (\f$0 < \arg z < \pi/2\f$)

  Use integral representation of Airy function [3, pp3] eq (3.3) or 
  [5, pp 8] eq (6.3)
  \f{align*}
   Ai'(z) &= b(z) \int_0^\infty (2+\frac{t}{\zeta})^{1/6} t^{1/6} e^{-t}\;\text{d}t &
   b(z) &= \frac{z}{\sqrt{3\pi}2^{2/3}\Gamma(\frac{7}{6})} \frac{e^{-\zeta}}{\sqrt{\zeta}}
   \f}
   Therefore:
   \f{align*}
   Aie(z) &= b(z) \int_0^\infty (2+\frac{t}{\zeta})^{1/6} t^{1/6} e^{-t}\;\text{d}t &
   b(z) &= \frac{z}{\sqrt{3\pi}2^{2/3}\Gamma(\frac{7}{6})} \frac{1}{\sqrt{\zeta}}
   \f}
   Compute using a Gauss Laguerre quadrature
   \attention assert \f$0 < \arg z < \pi/2\f$
*/     
static complex cairyAiep_mediumarg_Q1(const complex z)
{
  complex zeta, zetam12;
  complex bz;
  double t;
  complex integrand;
  complex sum;

  unsigned int k;
  
  assert(creal(z) >= 0);
  assert(cimag(z) >= 0);
  assert(ismedium(z));
  BUILD_ASSERT(ARRAY_SIZE(wp) == ARRAY_SIZE(xp));

  zeta = czetafast(z);
  zetam12 = 1 / csqrt(zeta);
  /* here bz is bz scaled ie without \f$e^{-\zeta}\f$ */
  bz = bzcoeff * z * zetam12;

  sum = 0;
  for(k = 0; k < ARRAY_SIZE(xp); k++)
  {
    t = xp[k];
    integrand = cpow( 2 + t/zeta, 1.0 /6.0);
    sum += integrand * wp[k];
  }

  return sum * bz;
}                       

/*!\brief Compute derivative of scaled airy function (\f$Aie\f$) for medium argument 
  and in second quadrant (\f$\pi/3 \ge \arg z > \pi/2\f$)

   Use integral representation of Airy function [3, pp3] eq (3.3) or 
   [5, pp 8] eq (6.3)
   \f{align*}
   Ai'(z) &= b(z) \int_0^\infty (2+\frac{t}{\zeta})^{1/6} t^{1/6} e^{-t}\;\text{d}t &
   b(z) &= \frac{z}{\sqrt{3\pi}2^{2/3}\Gamma(\frac{7}{6})} \frac{e^{-\zeta}}{\sqrt{\zeta}}
   \f}
   But this representation is not numerically stable for \f$\pi/3 \ge \arg z > \pi/2\f$
   The reason (see [5, pp9] is the presence of singularities in the integrand at 
   \f$t=-2\zeta\f$. This problem could be solved by turning the integration path over a 
   given angle \f$\tau=\frac{3}{2}(\arg z -\frac{\pi}{2})\f$, that can be done using 
   the substituation \f$t\rightarrow (t+i\tan \tau)\f$. Therefore:
   \f{align*}
   Ai'(z) &= b(z) \left(\frac{e^{i\tau}}{\cos \tau}\right)^{7/6}
         \int_0^\infty (2+\frac{t}{\tilde{\zeta}})^{1/6} 
	       e^{-it\tan\tau} t^{1/6} e^{-t}\;\text{d}t &
   b(z) &= \frac{z}{\sqrt{3\pi}2^{2/3}\Gamma(\frac{7}{6})} \frac{e^{-\zeta}}{\sqrt{\zeta}}
   \f}
   Where \f$\tilde\zeta=\cos \tau e^{-i\tau} \zeta\f$
   Therefore:
   \f{align*}
   Aie(z) &=  b(z) \left(\frac{e^{i\tau}}{\cos \tau}\right)^{7/6}
         \int_0^\infty (2+\frac{t}{\tilde{\zeta}})^{1/6} 
	       e^{-it\tan\tau} t^{1/6} e^{-t}\;\text{d}t &
   b(z) &= \frac{z}{\sqrt{3\pi}2^{2/3}\Gamma(\frac{7}{6})} \frac{1}{\sqrt{\zeta}}
   \f}

   Compute using a Gauss Laguerre quadrature
   \attention assert \f$0 < \arg z < \pi/2\f$
*/     
static complex cairyAiep_mediumarg_Q2(const complex z)
{
  complex zeta, zetam12 ,zetatilde;
  double tau;
  complex eitaucostau16;
  complex bz;
  double t;
  complex integrand;
  complex sum;

  unsigned int k;
  
  assert(cimag(z) >= 0);
  assert(ismedium(z));
  assert(creal(z) <= 0 && !arglargerpi2_3(z));
  BUILD_ASSERT(ARRAY_SIZE(wp) == ARRAY_SIZE(xp));

  tau = 3.0 / 2.0 * (carg(z) - M_PI_2);
  eitaucostau16 = cpow(cexp(I * tau)/cos(tau), 7.0 / 6.0);


  zeta = czetafast(z);
  zetam12 = 1 / csqrt(zeta);
  zetatilde = cos(tau) * cexp(-I * tau) * zeta;
  /* here az is az scaled ie without \f$e^{-\zeta}\f$ */
  bz = bzcoeff * z * zetam12;

  sum = 0;
  for(k = 0; k < ARRAY_SIZE(xp); k++)
  {
    t = xp[k];
    integrand = cpow( 2 + t/zetatilde, 1.0 /6.0);
    integrand *= cexp(-I * t * tan(tau));
    sum += integrand * wp[k];
  }

  return sum * bz * eitaucostau16;
}                       

/*!\brief Compute derivative of scaled Airy function of first kind in region I ie 
    \f$|\arg z|< 2\pi/3\f$ 
*/
static complex cairyAiep_mediumarg_regionI(const complex z)
{
  assert(cimag(z) >= 0);
  assert(isregionI(z));
  assert(ismedium(z));

  if(creal(z) > 0)
    return cairyAiep_mediumarg_Q1(z);
  return cairyAiep_mediumarg_Q2(z);
}

/*!\brief Compute scaled derivative of airy Ai function of first kind (\f$Aie'\f$)
          in upper region II upper complex plane (\f$\frac{2\pi}{3} \le \arg z \le \pi\f$)

   This function use connexion formula given by
   [3] eq (7.1) (see comment just after)
   for \f$\frac{2\pi}{3} \le \arg z \le \pi \f$

   The connexion formula is derivate using derivation of [3, pp 10] eq (7.1) for 
   \f$\frac{2\pi}{3} \le \arg z \le \pi \f$:
   \f[
   Aie'(z)=-e^{-\frac{4\pi i }{3}} e^{\frac{4}{3}z^\frac{3}{2}} 
                                   Aie'\left(e^{-\frac{2\pi i }{3}}z\right)
		        -e^{\frac{4\pi i }{3}}
			        Aie'\left(e^{\frac{2\pi i }{3}}z\right)
				   
   \f]

   \attention assert \f$\Im\text{m}\; z \ge 0\f$
   \note Because we assert \f$\Im\text{m}\; z \ge 0\f$ then
   \f$-\frac{\pi}{3} \le \arg e^{\frac{2\pi i }{3}} \le 0\f$
   we could therefore use the mirror property of derivative of scaled Airy function. 
   \note Because we assert \f$\Im\text{m}\; z \ge 0\f$  then
   \f$0 \le \arg e^{-\frac{2\pi i }{3}} \le \frac{\pi}{3}\f$
*/
static complex cairyAiep_mediumarg_regionII(const complex z)
{
  complex exp_43z32;
  complex ze2pi_3, zem2pi_3;
  complex Aizem2pi_3, Aize2pi_3;
  complex first, second;

  assert(cimag(z) >= 0);
  assert(ismedium(z));
  assert(isregionII(z));
  
  exp_43z32 = cexp_43z32(z);
  
  ze2pi_3 = z * M_EXP_I_2PI_3;
  zem2pi_3 = z * M_EXP_MI_2PI_3;

  /* quick test */
  assert(cimag(ze2pi_3) <= 0);
  assert(cimag(zem2pi_3) >= 0);
  
  /* Use reflexion property */
  Aize2pi_3 = conj(cairyAiep_mediumarg_regionI(conj(ze2pi_3)));
  Aizem2pi_3 = cairyAiep_mediumarg_regionI(zem2pi_3);

  first = - exp_43z32 * M_EXP_I_2PI_3 * Aizem2pi_3;
  second =  M_EXP_MI_2PI_3 * Aize2pi_3;
  
  return first - second;  
}

/*!\brief Compute derivative of scaled airy Ai function of first kind (\f$Aie'\f$)
          in upper complex plane (\f$\Im\text{m}\; z > 0\f$)

   This function use integral representation in 
   upper region I (\f$0 \le \arg z < \frac{2\pi}{3}\f$) and connexion formula given by
   [3] eq (7.1) for \f$\frac{2\pi}{3} \le \arg z \le \pi \f$

   \attention assert \f$\Im\text{m}\; z \ge 0\f$
*/
static complex cairyAiep_mediumarg(const complex z)
{
  assert(cimag(z) >= 0);
  assert(ismedium(z));

  DEBUG("z= " fmtc "arg=%e (%e°) norm2(z)=%e cabs(z)=%e\n", 
	cprint(z), carg(z), carg(z) * 180 / M_PI, norm2(z), cabs(z));

  /* region I */
  if(!arglargerpi2_3(z))
  {
    DEBUG("z is in region I\n");
    return cairyAiep_mediumarg_regionI(z);
  }

  /* Region II */
  DEBUG("z is in region II\n");
  return cairyAiep_mediumarg_regionII(z);
}


/*!\brief Compute derivative airy function of first kind (\f$Ai'\f$)
          in upper complex plane (\f$\Im\text{m}\; z > 0\f$)

   \attention assert \f$\Im\text{m}\; z \ge 0\f$
*/
static complex cairyAip_mediumarg(const complex z)
{
  complex zeta;
  complex expmzeta, scaled;

  assert(ismedium(z));

  zeta = czetafast(z);

  scaled = cairyAiep_mediumarg(z);
  expmzeta = cexp(-zeta);
  return scaled * expmzeta;
}


/*!\brief Compute \f$(-1)^k d_k\f$ coefficient of large argument expansion 

  According to (10.4.58) [1-p448] \f$d_k\f$ is defined by:
  \f{align*}
  c_1 &= 1 \\
  c_k &=  \frac{\Gamma(3k+\frac{1}{2})}{54^k k! \Gamma(k + \frac{1}{2})}
      &=  \frac{(2k + 1)(2k + 3)\dotsm(6k-1)}{216^k k!}
  \f}
  And
  \f{align*}
  d_0 &= 1 \\
  d_k &= - \frac{6k+1}{6k-1}c_k
  \f}
*/
static const double m1kdk[] = 
{
	/* 1 */
	1.0,
	/*-(6*1+1)/(6*1-1)*(3*5)/(216^1*1!)*/
	-0.9722222222222222222222222222222222222222e-1,
	/*-(6*2+1)/(6*2-1)*(5*7*9*11)/(216^2*2!)*/
	0.438850308641975308641975308641975308642e-1,
	/*-(6*3+1)/(6*3-1)*(7*9*11*13*15*17)/(216^3*3!)*/
	-0.4246283078989483310470964791952446273434e-1,
	/*-(6*4+1)/(6*4-1)*(9*11*13*15*17*19*21*23)/(216^4*4!)*/
	0.6266216349203230579688055682568714118783e-1,
	/*-(6*5+1)/(6*5-1)*(11*13*15*17*19*21*23*25*27*29)/(216^5*5!)*/
	-0.1241058960272750945365995472686525879637e0,
	/*-(6*6+1)/(6*6-1)*(13*15*17*19*21*23*25*27*29*31*33*35)/(216^6*6!)*/
	0.3082537649010791121244706347668153400116e0,
	/*-(6*7+1)/(6*7-1)*(15*17*19*21*23*25*27*...*39*41)/(216^7*7!)*/
	-0.9204799924129445709272387010397958069791e0,
	/*-(6*8+1)/(6*8-1)*(17*19*21*23*25*27*29*...*45*47)/(216^8*8!)*/
	0.3210493584648620907973650261091926694828e1,
	/*-(6*9+1)/(6*9-1)*(19*21*23*25*27*29*31*...*51*53)/(216^9*9!)*/
	-0.1280729308073562507270352766191763966996e2,
	/*-(6*10+1)/(6*10-1)*(21*23*25*27*29*31*...*57*59)/(216^(10)*10!)*/
	0.5750830351391427202784792351524962368467e2,
	/*-(6*11+1)/(6*11-1)*(23*25*27*29*31*33*...*63*65)/(216^(11)*11!)*/
	-0.2870332371092211077349530828987143464969e3,
	/*-(6*12+1)/(6*12-1)*(25*27*29*31*33*35*...*69*71)/(216^(12)*12!)*/
	0.1576357303337099717826796734206480988574e4,
	/*-(6*13+1)/(6*13-1)*(27*29*31*33*35*37*...*75*77)/(216^(13)*13!)*/
	-0.9446354823095931962917203933936059684732e4,
	/*-(6*14+1)/(6*14-1)*(29*31*33*35*37*39*...*81*83)/(216^(14)*14!)*/
	0.6133570666385205823144156720993205420295e5,
	/*-(6*15+1)/(6*15-1)*(31*33*35*37*39*41*...*87*89)/(216^(15)*15!)*/
	-0.4289524004000690702056279232746451901804e6,
	/*-(6*16+1)/(6*16-1)*(33*35*37*39*41*43*...*93*95)/(216^(16)*16!)*/
	0.3214536521400864829067001615998274242038e7,
	/*-(6*17+1)/(6*17-1)*(35*37*39*41*43*45*...*99*101)/(216^(17)*17!)*/
	-0.2569790838391132545132402844162019073394e8,
	/*-(6*18+1)/(6*18-1)*(37*39*41*43*45*47*...*105*107)/(216^(18)*18!)*/
	0.2182934208321603255352054236989171911959e9,
	/*-(6*19+1)/(6*19-1)*(39*41*43*45*47*49*...*111*113)/(216^(19)*19!)*/
	-0.1963523788991032752712502001911678390107e10,
	/*-(6*20+1)/(6*20-1)*(41*43*45*47*49*51*...*117*119)/(216^(20)*20!)*/
	0.1864393108810721585266530546676276293606e11,
	/*-(6*21+1)/(6*21-1)*(43*45*47*49*51*53*...*123*125)/(216^(21)*21!)*/
	-0.1863529963852938843791870115867629869396e12,
	/*-(6*22+1)/(6*22-1)*(45*47*49*51*53*55*...*129*131)/(216^(22)*22!)*/
	0.1955882932389842694320697012392635516333e13,
	/*-(6*23+1)/(6*23-1)*(47*49*51*53*55*57*...*135*137)/(216^(23)*23!)*/
	-0.2150644463519724977106616660546950490151e14,
	/*-(6*24+1)/(6*24-1)*(49*51*53*55*57*59*...*141*143)/(216^(24)*24!)*/
	0.2472369922906211612860123840379928905489e15,
	/*-(6*25+1)/(6*25-1)*(51*53*55*57*59*61*...*147*149)/(216^(25)*25!)*/
	-0.2965882430295212630916036338073544714235e16,
	/*-(6*26+1)/(6*26-1)*(53*55*57*59*61*63*...*153*155)/(216^(26)*26!)*/
	0.3706244000635465228366390921824488862185e17,
	/*-(6*27+1)/(6*27-1)*(55*57*59*61*63*65*...*159*161)/(216^(27)*27!)*/
	-0.4816782647945217540878439641969943986785e18,
	/*-(6*28+1)/(6*28-1)*(57*59*61*63*65*67*...*165*167)/(216^(28)*28!)*/
	0.6500984080751062701873088502894851484942e19,
	/*-(6*29+1)/(6*29-1)*(59*61*63*65*67*69*...*171*173)/(216^(29)*29!)*/
	-0.90991982643654122347816576387500974448e20,
	/*-(6*30+1)/(6*30-1)*(61*63*65*67*69*71*...*177*179)/(216^(30)*30!)*/
	0.1319088866907750709757953915010100931894e22,
	/*-(6*31+1)/(6*31-1)*(63*65*67*69*71*73*...*183*185)/(216^(31)*31!)*/
	-0.1978219607616628114145519327828544287333e23,
	/*-(6*32+1)/(6*32-1)*(65*67*69*71*73*75*...*189*191)/(216^(32)*32!)*/
	0.3065639370223598386092264218755129070279e24,
	/*-(6*33+1)/(6*33-1)*(67*69*71*73*75*77*...*195*197)/(216^(33)*33!)*/
	-0.4904119815775620835731518126711435220212e25,
	/*-(6*34+1)/(6*34-1)*(69*71*73*75*77*79*...*201*203)/(216^(34)*34!)*/
	0.809039537418702808214940194228926925524e26,
	/*-(6*35+1)/(6*35-1)*(71*73*75*77*79*81*...*207*209)/(216^(35)*35!)*/
	-0.1375142480406956245407560846801889960356e28,
	/*-(6*36+1)/(6*36-1)*(73*75*77*79*81*83*...*213*215)/(216^(36)*36!)*/
	0.2406127967357125254551277279514124821451e29,
	/*-(6*37+1)/(6*37-1)*(75*77*79*81*83*85*...*219*221)/(216^(37)*37!)*/
	-0.4330398100410561949304091184921348144313e30
};


/*!\brief Compute scalled airy function of first kind (\f$e^{\zeta}Ai\f$) 
  in upper region I (\f$0 < \arg z < \frac{2\pi}{3}\f$)

  Compute the scalled airy function of first kind (\f$e^{\zeta}Ai\f$) where 
  \f$\zeta\f$ is \f$\frac{2}{3}z^\frac{3}{2}\f$

  For large argument according to (10.4.61) [1-p448] 
  \f$|z|\f$ \f$A_i\f$ could be approximated by (\f$|\arg z| < \pi\f$):
  \f[
  Ai' = - \frac{1}{2} \pi^\frac{-1}{2} z^{1}{4} e^{-\zeta} \sum_0^\infty (-1)^k d_k \zeta^{-k}
  \f]
  and d_k defined by (10.4.58) [1-p448]:
  \f{align*}
  c_1 &= 1 \\
  c_k &=  \frac{\Gamma(3k+\frac{1}{2})}{54^k k! \Gamma(k + \frac{1}{2})}
       =  \frac{(2k + 1)(2k + 3)\dotsm(6k-1)}{216^k k!}
  \f}
  And 
  \f{align*}
  d_0 &= 1 \\
  d_k &= - \frac{6k+1}{6k-1}c_k
  \f}
  Therefore scaled \f$Ai'\f$ is defined by:
  \f[
  Ai' = -\frac{1}{2} \pi^\frac{-1}{2} z^{1}{4} \sum_0^\infty (-1)^k d_k \zeta^{-k}
  \f]

  According to [2, pp 475-476] this expansion is valid in region I, ie 
  (\f$|\arg z| < \frac{2\pi}{3}\f$) and \f$|z| > \rho\f$ where 
  \f$\rho \f$ is defined in [2, tab I pp 480].

  \attention assert \f$\Im\text{m}\; z \ge 0\f$
  \attention Computation are only valid for region I ie \f$|\arg z| < \frac{2\pi}{3}\f$
             (asserted)
*/
static complex cairyAiep_largearg_regionI(const complex z)
{
  complex z1_4;
  complex z32;
  complex coeff;
  complex invzeta;
  complex zmk;

  complex sum;
  complex aik;

  unsigned int k;

  assert(islarge(z));
  assert(cimag(z) >= 0);
  assert(isregionI(z));

  /* z1_4 = z^(1/4) */
  z1_4 = cz1_4(z);
  /* coeff = 1/2 * sqrt(pi) * z^(1/4) */
  coeff = - halfinversesquarerootpi * z1_4;

  DEBUG("constant=%e coeff=" fmtc " z^1/4=" fmtc "\n", 
	halfinversesquarerootpi, cprint(coeff), cprint(z1_4));

  /* underflow */
  if(cabs(z) > DBL_MAX_POW_2_3)
    goto overflow;

  /* overflow is avoided because |z| < DBL_MAX^(2/3) */
  z32 = cpow32fast(z);
  invzeta = (3 / 2.0) / z32;

  sum = 1;
  zmk = 1;
  
  for ( k = 0 ; k < ARRAY_SIZE(m1kdk) ; k++)
    {
      zmk *= invzeta;
      aik = zmk * m1kdk[k];
      
      if (norm2 (aik) < norm2 (sum) * sqr (DBL_EPSILON)
	  || norm2 (aik) < DBL_MIN)
	return coeff * sum;
      
      sum += aik;
      DEBUG("zmk=" fmtc " aik=" fmtc " sum=" fmtc "\n",
	    cprint(zmk), cprint(aik), cprint(sum));
    }
  /* Do not converge */
  assert (0 && sizeof ("large expansion do not converge!"));
  errno = ERANGE;
  (void) feraiseexcept (FE_INVALID);
  return coeff * sum;

 overflow:
  return coeff;
}






/*!\brief Derivative of Airy function of first kind (\$fAi'\f$)
   \todo Document
*/
complex cairyAip(const complex z)
{
  /* mirror symetry */
  if(cimag(z) < 0.0)
    return conj(cairyAip(conj(z)));

  if(use_small (z))
    return cairyAip_smallarg(z);

  return cairyAip_mediumarg(z);
}

/*!\brief Scaled derivative of Airy function of first kind (\f$Aie'\f$)
   \todo Document
*/
complex cairyAiep(const complex z)
{
  /* mirror symetry */
  if(cimag(z) < 0.0)
    return conj(cairyAiep(conj(z)));

  if(use_small (z))
    return cairyAiep_smallarg(z);

  return cairyAiep_mediumarg(z);
}

/*!\brief Derivative of Airy function of second kind (\f$Bi'\f$)
   \todo Document
*/
complex cairyBip(const complex z)
{
  /* mirror symetry */
  if(cimag(z) < 0.0)
    return conj(cairyBip(conj(z)));

  if(use_small (z))
    return cairyBip_smallarg(z);

}

/*!\brief Scaled derivative of Airy function of second kind (\f$Bie'\f$)
   \todo Document
*/
complex cairyBiep(const complex z)
{
  /* mirror symetry */
  if(cimag(z) < 0.0)
    return conj(cairyBiep(conj(z)));

  if(use_small (z))
    return cairyBiep_smallarg(z);

}



#ifdef TEST_CAIRY
#include "test_cairy.c"
#endif