home *** CD-ROM | disk | FTP | other *** search
/ The C Users' Group Library 1994 August / wc-cdrom-cusersgrouplibrary-1994-08.iso / listings / v_09_01 / 9n01125a < prev    next >
Text File  |  1990-11-28  |  3KB  |  143 lines
  1. //#include <stream.hpp>    //for Zortech C++
  2. #include <iostream.h> //for TURBO C++
  3. #include <math.h>
  4. #include "complex.hpp"
  5. #define errorcode -1.e60
  6. #define ABS(x) ((x)>0.?(x):-(x))
  7. #define max(a,b) ((a)>(b)?(a):(b))
  8. #define pi 3.141592653589793238462643383279
  9.  
  10.  
  11. void complex::print( char *ahead,char *behind)
  12.    {char *between="";
  13.    if(y>=0.)between="+";
  14.    printf("%s %e%s%e i %s",ahead,x,between,y,behind);}
  15.  
  16. double complex::abs()
  17.    { return sqrt(x*x+y*y);}
  18.  
  19. complex complex::cexp()
  20.    {double scale;
  21.    scale= exp(((double)x));
  22.    return complex(scale*cos(y),scale*sin(y));
  23.    }
  24. complex complex::clog()
  25.    {double mant,arg,mag;
  26.    mant = log((*this).abs());
  27.    arg= atan2(y,x);
  28.    return complex(mant,arg);
  29.    }
  30.  
  31. complex complex::operator^(double expon)
  32.    {
  33.    complex z;
  34.    z= (*this).clog( ) * expon;
  35.    return z.cexp();
  36.    }
  37.  
  38. complex complex::operator^(complex expon)
  39.    {
  40.    complex z;
  41.    z= (*this).clog( ) * expon;
  42.    return z.cexp();
  43.    }
  44.  
  45.  
  46. double Rzeta2( int k)
  47. // Riemann Zeta function of 2*(k+1) returned
  48. {double z[6]={1.64493406684822643647,1.08232323371113819152,
  49.    1.01734306198444913971, 1.00407735619794433938,
  50.    1.00099457512781808534,1.00024608655330804830};
  51. if(k<6)return z[k]; 
  52. // return 1+ 2^(-2k)+...+ 6^(-2k) in simplified form;
  53. if(k>=30)return 1.;
  54. double fk=-((k+1)<<1); double p2=pow(2.,fk);
  55. return 1.+((p2+pow(3.,fk))*(1.+p2)+pow(5.,fk));
  56. }
  57.  
  58. void prep( complex& u, complex& v, complex& f)
  59. {
  60. f=0.;
  61. do   {
  62.    f+=u^(-v);u++;
  63.    if(v.real()<1.)
  64.       { while(u.real()>2.){u-=1.;f-=u^(-v);}}
  65.    }while(u.real()<=v.real());
  66. }
  67.  
  68. #define tolabs 1.e-10
  69. #define tolrel 1.e-5
  70.  
  71.  
  72. int size=80,szused,jused;
  73.  
  74. complex complex::hurwitz(complex u)
  75. {
  76. complex b,c,f,t,p,q,*h,hold,tpu,ev;int i,j,k; double diff,tk,z,scale;
  77. if( *this==1. || ((u-*this+(*this).abs())==0.&&
  78. (*this).imaginary()==0.))
  79.   return complex(errorcode,0.);
  80. prep(u,*this,f);
  81. //u.print("modified a=","");f.print(" modifed sum=","\n");
  82. ev= -(*this);
  83. j=4;scale=.25;
  84. jused=j=max(j, (int)(((*this).abs()+.999)*scale));
  85. if(x<0.)jused=j=0;
  86. b=u+ ((double)j); c=u;
  87. for(i=0;i<j;i++){f+= c^ev;c+=1.;} 
  88. //f.print(" modifed sum=","\n");
  89. if( size%2)size++;//size must be even integer
  90. h=new complex[size+1] ;
  91. t=complex(-2.,0.);
  92. h[0]=1.;k=0; 
  93. tpu=2.*pi*b;
  94. tpu=1./(tpu*tpu);
  95. while(k<size){
  96.    tk=k<<1;
  97.    t*=(*this+tk)*(1.-*this-tk)*tpu;
  98.    z=Rzeta2(k);
  99.    k++;
  100.    h[k]=h[k-1]+t*z;
  101.    q=0.;j=k;
  102.    do   {
  103.       p=h[j-1];
  104.       if(h[j]==p)
  105.          {
  106.          h[j-1]=-errorcode;
  107.          if(p==-errorcode)h[j-1]=q;
  108.          }
  109.       else h[j-1]= q+1./(h[j]-p);
  110.       q=p;j--;
  111.       }while(j);
  112.    if(!(k%2))
  113.       {
  114.       diff=(hold-h[0]).abs();
  115.       if(diff<tolabs || diff<tolrel * h[0].abs())break;
  116.       hold=h[0];
  117.       }
  118.    };
  119. szused=k;
  120. ev= -(*this);
  121. tpu= 1.+ev;
  122. q=-h[0]*(b^tpu)/tpu;
  123. delete h;
  124. return f+.5*(b^ev)+q;
  125. }
  126.  
  127. main()
  128. {complex u,v,w;double a,b,c,d;
  129.  
  130. while(1)
  131.    {
  132.    cout << " enter v,u" ;
  133.        //scanf("%le%le%le%le",&a,&b,&c,&d);
  134.       cin >> a >> b >> c >> d ;
  135.    if(a==0.&&b==0.&&c==0.&&d==0.)break;
  136.       u=complex(c,d);v=complex(a,b);
  137.       w=v.hurwitz(u);
  138.        w.print(" Hurwitz=","\n");
  139.    cout  << " max k used=" << szused << ", J used: " << jused << "\n" ;
  140.    };
  141. }
  142.  
  143.