SuZex.cin

From TORI
Jump to navigation Jump to search

// SuZex.cin is the C++ implementation of function SuZex
// which is superfunction for the transfer function T(z)=zex(z)= z exp(z)
// See SuZex approximation for the details.

// Copyleft 2012 by Dmitrii Kouznetsov;

 //Example of the header:
 /*
 #include <math.h>
 #include <stdio.h>
 #include <stdlib.h>
 #define DB double
 #define DO(x,y) for(x=0;x<y;x++)
 using namespace std;
 #include<complex>
 typedef complex<double> z_type;
 #define Re(x) x.real()
 #define Im(x) x.imag()
 #define I z_type(0.,1.)
 #include "conto.cin"
 //#include "fsexp.cin"
 //#include "fslog.cin"
 */

z_type zex(z_type z) { return z*exp(z);}
z_type suzexo(z_type z){ int n,m; z_type L=log(-z); z_type c[21]; z_type s;
double a[21][21]=
{{-1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.},
{0., 0.5, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.},
{-0.16666666666666666, 0.25, -0.25, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.},
{-0.14583333333333334, 0.375, -0.3125, 0.125, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.},
{-0.16365740740740742, 0.4791666666666667, -0.53125, 0.2708333333333333, -0.0625,
0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.},
{-0.1894675925925926, 0.6487268518518519, -0.8645833333333334, 0.578125, -0.20052083333333334, 0.03125,
0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.},
{-0.22663111772486771, 0.8927662037037037, -1.4053819444444444, 1.1536458333333333, -0.5338541666666666, 0.1359375, -0.015625,
0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.},
{-0.2772791969797178, 1.239592013888889, -2.265031828703704, 2.216435185185185, -1.2763671875, 0.44166666666666665, -0.087109375,
0.0078125, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.},
{-0.3462407775931875, 1.7289127948633156, -3.6116999421296296, 4.1282600308641975, -2.854618778935185,
1.2419270833333333, -0.3379991319444444, 0.05368303571428571, -0.00390625, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.},
{-0.4416612459670082, 2.4225398966010014, -5.695903759507275, 7.481679928626543, -6.071601924189815, 3.1901204427083334,
-1.1004448784722223, 0.24412667410714287, -0.032149832589285714, 0.001953125, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.},
{-0.5750503318438397, 3.4195761781355416, -8.904301621256142, 13.234012896825396, -12.387900872878086,
7.666662145543982, -3.2086561414930554, 0.9080953931051587, -0.16865408761160713, 0.018837580605158732, -0.0009765625,
0., 0., 0., 0., 0., 0., 0., 0., 0., 0.},
{-0.7590483780036865, 4.872564914208889, -13.85598530050081, 22.941559420715624, -24.390718169573965,
17.460022032937886, -8.632101704161844, 2.9751346648685515, -0.7086426265656002, 0.11248517717633928, -0.010848950582837302,
0.00048828125, 0., 0., 0., 0., 0., 0., 0., 0., 0.},
{-1.004225435438346, 6.9905727251265635, -21.545687392877074, 39.182750311359435, -46.60769821586042, 37.9988728199577,
-21.77607288501881, 8.886511650287286, -2.5856723119342138, 0.5286710062985698, -0.07291558159722222, 0.006161750033820346,
-0.000244140625, 0., 0., 0., 0., 0., 0., 0., 0.},
{-1.324367258711613, 10.02275169291253, -33.492205053099866, 66.27369784024671, -86.97581836388929, 79.5894440905974,
-52.05348199746358, 24.66389493266111, -8.513126871825525, 2.1317655062128837, -0.38009394489268145, 0.04616735505163239,
-0.0034596849141526878, 0.0001220703125, 0., 0., 0., 0., 0., 0., 0.},
{-1.765838606129714, 14.281946657437558, -51.825733451743794, 111.28532737735638, -159.46688040237638,
161.5608677547437, -118.88109243776208, 64.38542946379413, -25.83747150199123, 7.687203653415184, -1.6822828267953593,
0.2649616424575226, -0.028660799570528573, 0.0019239424177168317, -0.00006103515625, 0., 0., 0., 0., 0., 0.},
{-2.4453295113397027, 20.384762874691635, -79.47016669126273, 185.20699731803768, -288.39342903373137,
319.9807544809364, -261.3916309123107, 159.56531377235643, -73.28007587330262,
25.374828078366953, -6.606544153459068, 1.2794918394983246, -0.1799314263212159,
0.017497047884163362, -0.0010612010161875885, 0.000030517578125, 0., 0., 0., 0., 0.},
{-3.472868025763881, 29.755017528063433, -121.2741348443979, 304.5239431690528,
-514.610709152941, 621.4198636944384, -557.3368214307372, 378.5159493573904,
-196.20535170900774, 77.82525925989691, -23.603134539423095, 5.44450530408303,
-0.9429602728261576, 0.1194755555243684, -0.010528913584758572, 0.0005812326643625472,
-0.0000152587890625, 0., 0., 0., 0.},
{-4.613171346702781, 44.39688698302471, -187.09589191646856, 495.8720203103204,
-904.4187338107077, 1185.548137407219, -1159.0132176158231, 866.0241149874477,
-500.2758720467312, 224.2176840217891, -77.95303764062392, 20.961038432504814,
-4.328004726805225, 0.6762887253792872, -0.07780318693217438, 0.006257000696877798,
-0.0003164092474738532, 7.62939453125e-6, 0., 0., 0.},
{-5.273050484527783, 63.71698561183738, -293.3339373818455, 809.223685904566,
-1567.9214126035747, 2220.7277895628836, -2357.82881491874, 1923.1719087140677,
-1224.4150653842441, 612.3847140576257, -240.77243443992216, 74.2602772858538,
-17.884781187781225, 3.334455327400953, -0.4736586312099147, 0.049810412507743526,
-0.003677767515730688, 0.00017132547534001787, -3.814697265625e-6, 0., 0.},
{-7.083324624421979, 81.95247240893262, -449.3226503471503, 1333.5026446614604, -2705.8669603251315,
4089.414578728234, -4695.066740933935, 4161.496488889609, -2895.974174290078, 1598.6304816010704,
-702.1516955747055, 245.06996838649602, -67.73177677852487, 14.736875454771372, -2.4994954306270327,
0.3248890060201511, -0.03141381618433806, 0.00214088576116657, -0.00009232912728448859, 1.9073486328125e-6, 0.},
{-18.214476068392667, 111.80948244868608, -634.4236872182383, 2164.4934901545644, -4686.6900918162155,
7456.441210014379, -9163.224335014025, 8787.986445778999, -6649.857698257051, 4017.0643233450655, -1949.7063293884232,
760.8547074429804, -238.0908620446758, 59.469804480097125, -11.776087325864497,
1.828774790094764, -0.21876253685476346, 0.019549158283135085, -0.0012355455420681165,
0.000049547951834558144, -9.5367431640625e-7}};

// #include "SuZexCoef20.cin"
c[0]=-1.; c[1]=.5*L;
for(n=2;n<21;n++){ s=a[n][n];
                   for(m=n-1; m>=0; m--){ s*=L; s+=a[n][m]; }
                   c[n]=s; }
s=c[20]/z; for(n=19;n>=0;n--){ s+=c[n]; s/=z;}
return s;
}
 z_type suzex2008t12(z_type z){ int n,m=80; z_type s;
 double c[112]={
  0.09263113502690032,0.008211206810042051,
  0.000712647069880373,0.000061002792135026616,
  5.168897004880931e-6,4.3445677647721587e-7,
  3.627510329345628e-8,3.011749762480825e-9,
  2.4882932760219414e-10,2.0469528670147353e-11,
  1.6774049075221718e-12,1.3697934125298772e-13,
  1.1150519031183394e-14,9.050499680673825e-16,
  7.32631450169268e-17,5.915874252987347e-18,
  4.765920796716144e-19,3.8311975285996394e-20,
  3.0735487052404278e-21,2.4610209884198943e-22,
  1.967011009647219e-23,1.5694830383038919e-24,
  1.250267839550047e-25,9.944454380674088e-27,
  7.898093109767335e-28,6.264065061257306e-29,
  4.961466114678482e-30,3.9247189079708547e-31,
  3.1008090414779705e-32,2.446983078569256e-33,
  1.9288411564301896e-34,1.5187638797678427e-35,
  1.1946203045539106e-36,9.387107747358053e-38,
  7.369048831404122e-39,5.779407456250524e-40,
  4.5285704598919127e-41,3.5453392367706237e-42,
  2.773227658284088e-43,2.167483778760481e-44,
  1.692699001606891e-45,1.320892724326662e-46,
  1.0299808302576541e-47,8.025533074423089e-49,
  6.249002602542793e-50,4.8623709172587175e-51,
  3.780890726937988e-52,2.9380322310816133e-53,
  2.2816170670960654e-54,1.7707603355215024e-55,
  1.373454866536895e-56,1.0646651593426265e-57,
  8.248252131650093e-59,6.386559879065092e-60,
  4.9423548308245524e-61,3.8226822480200836e-62,
  2.955119712176247e-63,2.2832823166625066e-64,
  1.7633018554609683e-65,1.3610714134613786e-66,
  1.0500908268140331e-67,8.097833892562851e-69,
  6.241817958992842e-70,4.809031344505806e-71,
  3.7034995072112494e-72,2.850879848329411e-73,
  2.193618572531656e-74,1.6871838457734253e-75,
  1.297138043467975e-76,9.968636715255186e-78,
  7.657981460120605e-79,5.8806442027849044e-80,
  4.514093360695293e-81,3.4638103401082574e-82,
  2.6569196029061346e-83,2.0372587207723823e-84,
  1.56156498861741e-85,1.196527295125126e-86,
  9.165081177043452e-88,7.017842204645432e-89,
  5.371886263590005e-90,4.1106285121843686e-91,
  3.1444890288849017e-92,2.404664710358596e-93,
  1.8383308299060688e-94,1.4049457397975093e-95,
  1.0734063622275724e-96,8.19859532315087e-98,
  6.260186874436073e-99,4.778697590779901e-100,
  3.6467668355383486e-101,2.7821746766120975e-102,
  2.1219758899253086e-103,1.6179973905647313e-104,
  1.2333834974885623e-105,9.399461648725002e-107,
  7.161333609088019e-108,5.454719486363905e-109,
  4.153746652896943e-110,3.1622638017912685e-111,
  2.4068449211041933e-112,1.8314346120767418e-113,
  1.3932510247170333e-114,1.059652059254174e-115,
  8.057391237803329e-117,6.125254358912119e-118,
  4.655362973359923e-119,3.5373979441576056e-120,
  2.6873016732269505e-121,2.0410432859962447e-122,
  1.5498598577315127e-123,1.1766252828987956e-124};
 // #include "SuZexTay2008co.cin"
 s=c[m]; for(n=m-1;n>=0;n--){ s*=z; s+=c[n];}
 return s;}

 z_type SuZexTay0(z_type z){ int n,m=96; z_type s;
 // Coefficients of the Taylor expansion at zero of function SuZex. Furst ten values seem to be precize.
 double SuZexTay0co[97]={
   1.,0.7136859972397819,
   0.44760159770758723,0.260172774210734,
   0.14351435561731263,0.07612265909856022,
   0.039149270539685414,0.01963286277085933,
   0.009639736455495637,0.004648327410935626,
   0.002206510699406525,0.0010330188519909768,
   0.000477708623371847,0.00021848101908428043,
   0.00009892678352844481,0.00004438614595334552,
   0.000019748886668485166,8.719325500381778e-6,
   3.822215192158781e-6,1.6643935234951081e-6,
   7.202727821897625e-7,3.098905541329569e-7,
   1.3259961066070563e-7,5.644630910362613e-8,
   2.3911772971848777e-8,1.0082842042825059e-8,
   4.233027550988762e-9,1.769739581178266e-9,
   7.369587989119768e-10,3.0572517487417084e-10,
   1.2637053312940668e-10,5.2053996884250256e-11,
   2.1370670158299365e-11,8.745723776272641e-12,
   3.568130503659941e-12,1.4514540270404245e-12,
   5.887492850084059e-13,2.3815849132007153e-13,
   9.608412431405376e-14,3.866570577462619e-14,
   1.5521210050724034e-14,6.215628529414511e-15,
   2.4833444698339505e-15,9.899474755844535e-16,
   3.937674510545802e-16,1.5629592062539282e-16,
   6.19101603512115e-17,2.4474152450402498e-17,
   9.65625256611268e-18,3.802659539422078e-18,
   1.4947402770926363e-18,5.864947187812793e-19,
   2.297218140115094e-19,8.982535748615876e-20,
   3.506487297990436e-20,1.3665937660186885e-20,
   5.317616148368429e-21,2.065954622336573e-21,
   8.014308213914461e-22,3.104326526535793e-22,
   1.200712142614215e-22,4.637608493477486e-23,
   1.7887347815479634e-23,6.88980682205225e-24,
   2.650265866508781e-24,1.0181336565847141e-24,
   3.906288694553889e-25,1.4968521982261763e-25,
   5.728738794347915e-26,2.189852009234572e-26,
   8.360956138353809e-27,3.188541158152261e-27,
   1.2145969233506357e-27,4.621524566640595e-28,
   1.7565474638078776e-28,6.669056170567843e-29,
   2.5293375859074745e-29,9.582856551540693e-30,
   3.626908015616651e-30,1.371320077917939e-30,
   5.179752729604282e-31,1.954582817311684e-31,
   7.368528020156649e-32,2.7752083037388775e-32,
   1.0442512518764559e-32,3.925682914243774e-33,
   1.4744573317046644e-33,5.533020374821927e-34,
   2.0744771710938354e-34,7.771257125665694e-35,
   2.9082833536311864e-35,1.0883331081199238e-35,
   4.0535616084969196e-36,1.5357809956671497e-36,
   5.370444540799253e-37,2.5415674051776864e-37,
   1.8802501580541168e-38};
// #include "SuZexTay0co.cin"
 s=SuZexTay0co[m]; for(n=m-1;n>=0;n--){ s*=z; s+=SuZexTay0co[n];}
 return s;}

 z_type suzex(z_type z){int m,n; z_type s;
 if( abs(z) < 1.6 ) return SuZexTay0(z) ; // I made the Taylor expansion for this case
 if( Re(z)>0 && fabs(Im(z))<1.5){n=int(Re(z)+.5); s=SuZexTay0(z-(0.+n)); DO(m,n) s=zex(s); return s;}
 z+=-1.1259817765745026; // WARNING! ARGUEMENT CHANGES ITS VALUE!
 if( abs(z+12.) < 8.1 ) return suzex2008t12(z+12.) ; // I made the Taylor expansion for this case
 if( Re(z)<-12. || fabs(Im(z)) > 8. ) return suzexo(z); // definitely, |z|>8
 n= int(Re(z)+12.); s=suzex2008t12(z+12.-(0.+n)); DO(m,n) s=zex(s); return s;
 }

// Copyleft 2012 by Dmitrii Kouznetsov //