AuZex.cin

// AuZex.cin is numerical implementation in C++ of function AuZex, with is inverse function of SuZex and the Abel function of function LambertW; the AuZex approximations are used in the algorithm.

// #include "tania.cin"
// #include "SuZex.cin"
// z_type zex(z_type z) { return z*exp(z) ; }
// LambertW is suppozed to be included separately

z_type AuZexAsy(z_type z){ int m=15,n; z_type s;
DB AuZexAsyCo[66]={
//     1.1259817765745026,
       1.1259817765744977,
-0.16666666666666667,0.0625,-0.035185185185185185,   0.020833333333333333,-0.0097619047619047619,
  0.00035686728395061728,0.0057788485764676241,    -0.0054935515873015873,-0.0025850528358244408,
  0.012198640046296296,-0.0064941110501851844,    -0.026451479667987191,0.047851552440450232,
  0.053758729874794383,-0.27073626193208126,    -0.0065521186641040204,1.6278812626136699,
  -1.6076876900940989,-10.838187174665133,    24.785092910583443,78.1860907900493,
  -340.82033939567854,-562.09669036948494,    4908.2972195975399,2830.3348202280227,
  -77032.682625139301,31753.861334668919,    1.3308075281671887e6,-1.9538825394585991e6,
  -2.5276072666444939e7,6.8487621133842322e7,    5.2387590408716527e8,-2.2452803031477711e9,
  -1.1695088667110365e10,7.5450114071773863e10,    2.7518770516521089e11,-2.679843770347832e12,
  -6.5513897892809159e12,1.0183199942874624e14,    1.4245100909222102e14,-4.1590749536758109e15,
  -1.7035860994968719e15,1.8277639374581802e17,    -1.0679787375232343e17,-8.6360007004109277e18,
  1.4065535508735497e19,4.3787928905903592e20,    -1.2109904910255305e21,-2.3763215966256672e22,
  9.6017875634223775e22,1.3758361324208134e24,    -7.5814699303356225e24,-8.4659910118430712e25,
  6.1398187270988087e26,5.5112074423601013e27,    -5.1693736994996277e28,-3.773953811181416e29,
  4.5556026027743353e30,2.6980978199281418e31,    -4.2167875664888744e32,-1.9922345654966198e33,
  4.1064414540582216e34,1.4932212480440397e35,    -4.2098201366546206e36 };
 s=AuZexAsyCo[m]*z;   for(n=m-1;n>0;n--){s+=AuZexAsyCo[n]; s*=z;}
 return s + AuZexAsyCo[0] + .5*log(z)-1./z;}

z_type AuZex01(z_type z){ int n, m=39; 
 double AuZexCo[46]={0,
  1.4011764331478447,-1.2313176379841106,   1.1612567820116564,-1.123126930577658,
  1.0992876544297898,-1.0830479804216504,   1.0713113178859344,-1.0624516150969114,
  1.0555359459048546,-1.049992436340172,    1.0454519353364795,-1.0416660012642422,
  1.0384615628127265,-1.035714498198056,    1.0333335508961627,-1.0312501692201073,
  1.0294118808316126,-1.0277778512023086,   1.0263158327950879,-1.0250000237496548,
  1.0238095356162695,-1.0227272776768976,   1.0217391317101858,-1.0208333328204302,
  1.019999998732505,-1.019230767751615,     1.0185185170386828,-1.01785714143099,
  1.0172413778515972,-1.0166666648687321,   1.016131566418359,-1.015697093269014,
  1.0152, -1.0148, 1.0141, -1.014, 1.013, -1.012, 1.011,-1.010, 1.008, -1.006, 1.004, -1.002, 1. };
z_type s=AuZexCo[m];
for(n=m-1;n>0;n--) {s+= AuZexCo[n]; s*=z; }  return s;  }

z_type auzex(z_type z){ int n; 
 for(n=0;n<80;n++){ 
if(abs(z-1.)<.4) return AuZex01(z-1.)+(0.+n); 
if(abs(z)<.1) return AuZexAsy(z)+(0.+n);; 
z=LambertW(z); }
printf("z= %10.5f %10.5f n=%2d what to do?\n",Re(z), Im(z), n); getchar(); 
return 0.;}

//Copyleft 2012 by Dmitrii Kouznetsov // // // // // //