Doya.cin

// doya.cin is the C++ complex(double) implementation of the Tania function and the Doya function.

/*

Header

The typical header required for the Doya function is below: */

#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 "doya.cin"

/* After such a header, the functions can be defined as follows:

Tania and Doya

They are short; so, let them stay together */

z_type ArcTania(z_type z) {return z + log(z) - 1. ;}
z_type ArcTaniap(z_type z) {return 1. + 1./z ;}

z_type TaniaTay(z_type z) { int n; z_type s;
s=1.+z*(.5+z*(1./16.+z*(-1./192.+z*(-1./3072.+z*(1.3/6144.+z*(-4.7/147456.
//+z*(7.3/4128768.) //reserve term for testing
)))))); DO(n,2) s+=(z-ArcTania(s))/ArcTaniap(s); return s ; }

z_type TaniaNega(z_type z){int n;z_type s=exp(z-exp(z)+1.); 
DO(n,3) s+=(z-ArcTania(s))/ArcTaniap(s); return s ; }

z_type TaniaBig0(z_type z){int n;z_type  L=log(z), s=z-L+1.;  s-=(1.-L)/z;  
DO(n,3) s+=(z-ArcTania(s))/ArcTaniap(s);  return s ; }

z_type TaniaS(z_type z){int n; z_type s,t=z+z_type(2.,-M_PI);t*=2./9.; t=I*sqrt(t);
s=-1.+t*(3.+t*(-3.+t*(.75+t*(.3+t*(.9/16.+t*(-.3/7.+t*(-12.51/224. //+t*(-.9/28.)
)))))));
DO(n,3) s+=(z-ArcTania(s))/ArcTaniap(s); return s ; }

z_type Tania(z_type z){ z_type t;
if( fabs(Im(z))< M_PI && Re(z)<-2.51) return TaniaNega(z);
if( abs(z)>7. || Re(z)>3.8 ) return TaniaBig0(z);
if( Im(z) > .7 ) return TaniaS(z);
if( Im(z) < -.7) return conj(TaniaS(conj(z)));
return TaniaTay(z);
}

z_type Doya(z_type t, z_type z){ return Tania(t+ArcTania(z)) ;}

/*

Algorithm and the use

The Doya function is implemented through the Tania function. Both return complex values; the input may be also long(double), but not of the "int" type. If need to use with other type, the conversion is necessary. It is sufficient to add ".0" to the integer parameter to convert it to "double". The same about the output; even if the real value is expected, the real part of the result should be used.

The TORI commands are screened; the C++ is not supposed to complain if loaded as is */