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 // // // // // //