自学教程：C++ GSL_SET_COMPLEX函数代码示例

gsl_complexgsl_complex_sub (gsl_complex a, gsl_complex b){                               /* z=a-b */  double ar = GSL_REAL (a), ai = GSL_IMAG (a);  double br = GSL_REAL (b), bi = GSL_IMAG (b);  gsl_complex z;  GSL_SET_COMPLEX (&z, ar - br, ai - bi);  return z;}

gsl_complexgsl_complex_add (gsl_complex a, gsl_complex b){                               /* z=a+b */  double ar = GSL_REAL (a), ai = GSL_IMAG (a);  double br = GSL_REAL (b), bi = GSL_IMAG (b);  gsl_complex z;  GSL_SET_COMPLEX (&z, ar + br, ai + bi);  return z;}

 /** Subtraction operator (vector) */ vector<complex> vector<complex>::operator-(const vector<complex>& v) const {     vector<complex> v1(_vector);     gsl_complex z1;     GSL_SET_COMPLEX(&z1, -1., 0.);     if (gsl_blas_zaxpy(z1, v.as_gsl_type_ptr(), v1.as_gsl_type_ptr())) {         std::cout << "/n Error in vector<complex> -" << std::endl;         exit(EXIT_FAILURE);     }     return v1; }

/* ------------------------------------------------------ */gsl_complex gsl_complex_min_real (gsl_complex a, gsl_complex b){  gsl_complex z;  double min;	  /* just consider real parts */  min = GSL_REAL(a) < GSL_REAL(b) ? GSL_REAL(a) : GSL_REAL(b);  GSL_SET_COMPLEX (&z, min, 0);	  return z;}

gsl_complexgsl_complex_cos (gsl_complex a){				/* z = cos(a) */  double R = GSL_REAL (a), I = GSL_IMAG (a);  gsl_complex z;  if (I == 0.0)     {      /* avoid returing negative zero (-0.0) for the imaginary part  */      GSL_SET_COMPLEX (&z, cos (R), 0.0);      }   else     {      GSL_SET_COMPLEX (&z, cos (R) * cosh (I), sin (R) * sinh (-I));    }  return z;}

gsl_complexgsl_complex_pow_real (gsl_complex a, double b){				/* z=a^b */  gsl_complex z;  if (GSL_REAL (a) == 0 && GSL_IMAG (a) == 0)    {      GSL_SET_COMPLEX (&z, 0, 0);    }  else    {      double logr = gsl_complex_logabs (a);      double theta = gsl_complex_arg (a);      double rho = exp (logr * b);      double beta = theta * b;      GSL_SET_COMPLEX (&z, rho * cos (beta), rho * sin (beta));    }  return z;}

/* ------------------------------------------------------ */gsl_complex gsl_complex_max_real (gsl_complex a, gsl_complex b){  gsl_complex z;  double max;	  /* just consider real parts */  max = GSL_REAL(a) > GSL_REAL(b) ? GSL_REAL(a) : GSL_REAL(b);  GSL_SET_COMPLEX (&z, max, 0);	  return z;}

complex_spinor complex_spinor::operator+(const complex_spinor & complex_spinor_param)const{	gsl_complex aux_UP;	gsl_complex aux_DOWN;	GSL_SET_COMPLEX(&aux_UP, 0, 0);	GSL_SET_COMPLEX(&aux_DOWN, 0, 0);		int i;		int nSitios = this->_numSites;	complex_spinor complex_spinor_res(nSitios);		for(i=0; i<nSitios ; i++){	aux_UP = gsl_complex_add(complex_spinor_param.complex_spinor_get(i, UP) , this->complex_spinor_get(i, UP));   	aux_DOWN = gsl_complex_add(complex_spinor_param.complex_spinor_get(i, DOWN) , this->complex_spinor_get(i,DOWN));	complex_spinor_res.complex_spinor_set(i, UP, aux_UP);	complex_spinor_res.complex_spinor_set(i, DOWN, aux_DOWN);	}	return complex_spinor_res;}

gsl_complex gsl_complex_pow(gsl_complex a, gsl_complex b){				/* z=a^b */    gsl_complex z;    if (GSL_REAL(a) == 0 && GSL_IMAG(a) == 0.0) {	GSL_SET_COMPLEX(&z, 0.0, 0.0);    } else {	double logr = gsl_complex_logabs(a);	double theta = gsl_complex_arg(a);	double br = GSL_REAL(b), bi = GSL_IMAG(b);	double rho = exp(logr * br - bi * theta);	double beta = theta * br + bi * logr;	GSL_SET_COMPLEX(&z, rho * cos(beta), rho * sin(beta));    }    return z;}

gsl_complexgsl_complex_sqrt (gsl_complex a){				/* z=sqrt(a) */  gsl_complex z;  if (GSL_REAL (a) == 0.0 && GSL_IMAG (a) == 0.0)    {      GSL_SET_COMPLEX (&z, 0, 0);    }  else    {      double x = fabs (GSL_REAL (a));      double y = fabs (GSL_IMAG (a));      double w;      if (x >= y)	{	  double t = y / x;	  w = sqrt (x) * sqrt (0.5 * (1.0 + sqrt (1.0 + t * t)));	}      else	{	  double t = x / y;	  w = sqrt (y) * sqrt (0.5 * (t + sqrt (1.0 + t * t)));	}      if (GSL_REAL (a) >= 0.0)	{	  double ai = GSL_IMAG (a);	  GSL_SET_COMPLEX (&z, w, ai / (2.0 * w));	}      else	{	  double ai = GSL_IMAG (a);	  double vi = (ai >= 0) ? w : -w;	  GSL_SET_COMPLEX (&z, ai / (2.0 * vi), vi);	}    }  return z;}

static intdo_ccmul(lua_State *L, mMatComplex *a, double b_re, double b_im){    mMatComplex *r = qlua_newMatComplex(L, a->l_size, a->r_size);    gsl_complex z;    gsl_matrix_complex_memcpy(r->m, a->m);    GSL_SET_COMPLEX(&z, b_re, b_im);    gsl_matrix_complex_scale(r->m, z);    return 1;}

 void complex::assign(const double& real=0., const double& imag=0.,                 bool polar=false) {   if (polar)     _complex = gsl_complex_polar(real,imag);   else     {       GSL_SET_COMPLEX(&_complex, real, imag);       // *_complex = gsl_complex_rect(real,imag);     } }

static intcm_mul_cm(lua_State *L){    mMatComplex *a = qlua_checkMatComplex(L, 1);    mMatComplex *b = qlua_checkMatComplex(L, 2);    int al = a->l_size;    int ar = a->r_size;    int bl = b->l_size;    int br = b->r_size;    mMatComplex *r = qlua_newMatComplex(L, al, br);    gsl_complex z1, z0;    if (ar != bl)        return luaL_error(L, "matrix sizes mismatch in m * m");    GSL_SET_COMPLEX(&z1, 1.0, 0.0);    GSL_SET_COMPLEX(&z0, 0.0, 0.0);    gsl_blas_zgemm(CblasNoTrans, CblasNoTrans, z1, a->m, b->m, z0, r->m);    return 1;}

void CModel::AllFacetHarmonics(int bw, double r, gsl_complex *coeff){  for(int i=0; i<bw*bw; i++) GSL_SET_COMPLEX(&coeff[i], 0.0, 0.0);    gsl_complex *temp = new gsl_complex[bw*bw];    for(int i=0; i<f_no; i++){    FacetHarmonics(i, bw, r, temp);    for(int j=0; j<bw*bw; j++) coeff[j] = gsl_complex_add(coeff[j],temp[j]);  }	  delete [] temp;}

gsl_complex gsl_complex_arccos(gsl_complex a){				/* z = arccos(a) */    double R = GSL_REAL(a), I = GSL_IMAG(a);    gsl_complex z;    if (I == 0) {	z = gsl_complex_arccos_real(R);    } else {	double x = fabs(R), y = fabs(I);	double r = hypot(x + 1, y), s = hypot(x - 1, y);	double A = 0.5 * (r + s);	double B = x / A;	double y2 = y * y;	double real, imag;	const double A_crossover = 1.5, B_crossover = 0.6417;	if (B <= B_crossover) {	    real = acos(B);	} else {	    if (x <= 1) {		double D =		    0.5 * (A + x) * (y2 / (r + x + 1) + (s + (1 - x)));		real = atan(sqrt(D) / x);	    } else {		double Apx = A + x;		double D = 0.5 * (Apx / (r + x + 1) + Apx / (s + (x - 1)));		real = atan((y * sqrt(D)) / x);	    }	}	if (A <= A_crossover) {	    double Am1;	    if (x < 1) {		Am1 = 0.5 * (y2 / (r + (x + 1)) + y2 / (s + (1 - x)));	    } else {		Am1 = 0.5 * (y2 / (r + (x + 1)) + (s + (x - 1)));	    }	    imag = log1p(Am1 + sqrt(Am1 * (A + 1)));	} else {	    imag = log(A + sqrt(A * A - 1));	}	GSL_SET_COMPLEX(&z, (R >= 0) ? real : M_PI - real,			(I >= 0) ? -imag : imag);    }    return z;}

static VALUE Vector_complex_set_all(VALUE obj, VALUE va) {  gsl_vector_complex * v;  gsl_complex z;  double d0, d1;  d0 = NUM2DBL(rb_ary_entry(va,0));  d1 = NUM2DBL(rb_ary_entry(va,1));  GSL_SET_COMPLEX(&z, d0, d1);  Data_Get_Struct(obj, gsl_vector_complex, v);  gsl_vector_complex_set_all(v, z);  return obj;}

gsl_complex gsl_complex_arccot(gsl_complex a){				/* z = arccot(a) */    gsl_complex z;    if (GSL_REAL(a) == 0.0 && GSL_IMAG(a) == 0.0) {	GSL_SET_COMPLEX(&z, M_PI_2, 0);    } else {	z = gsl_complex_inverse(a);	z = gsl_complex_arctan(z);    }    return z;}

//go through all sets of Bethe rapidities and calculate g1 for all identical sets. return value is sum of all these terms.gsl_complex exact_integrationDE(double orteins, double ortzwei, double** kEp_plusminus, int anzahl_sets, double** coeffoverlap){  //variables to make loop easily readable, could delete from final version  double c, Energie, Energieprime, normierung, normierungprime;  double k[N];  double kprime[N];//corresponds to primed Bethe rapidity set in papers  gsl_complex overlap, overlapprime, integral;  GSL_SET_COMPLEX(&integral, 0.0, 0.0);//sum of all values, technically no needed but good for quick check of initial shape of g1    int zaehler_integrale = 0;//counter of total sets of integrals    double Nminusone_fak = 1.0;//(N-1) factorial, see definition of integrals on restricted spatial domain  for(int j=1; j<=(N-1); j++)    Nminusone_fak = Nminusone_fak * ((double) j);  for(int zaehler_eigenstate = 0 ; zaehler_eigenstate < anzahl_sets; zaehler_eigenstate++){//loop over all Bethe rapidity sets        c = kEp_plusminus[zaehler_eigenstate][0];//interaction strength    for(int bb=0; bb<N; bb++)      k[bb]=kEp_plusminus[zaehler_eigenstate][bb+1];//Bethe rapidities    Energie = kEp_plusminus[zaehler_eigenstate][N+1];//energy    normierung=kEp_plusminus[zaehler_eigenstate][N+6];//norm        GSL_SET_COMPLEX(&overlap, coeffoverlap[zaehler_eigenstate][0], coeffoverlap[zaehler_eigenstate][1]);//overlap of set with initial state      //calculate g1_DE for this Bethe set    gsl_complex integralwert = integraleigenfunction(k, k, c, normierung, normierung, orteins, ortzwei);        integralwert = gsl_complex_mul_real(integralwert, (Nminusone_fak*Laenge));//integral evaluated in Rp only -> missing factor of (N-1)!. Laenge (=system length): definition of g1 has prefactor of N!/(n*(N-1)!) = Laenge    integral = gsl_complex_add(integral, gsl_complex_mul(gsl_complex_mul(overlap, gsl_complex_conjugate(overlap)), integralwert));       }//end for loop zaehler_eigenstate     return integral;};

gsl_complexgsl_complex_arccos_real (double a){				/* z = arccos(a) */  gsl_complex z;  if (fabs (a) <= 1.0)    {      GSL_SET_COMPLEX (&z, acos (a), 0);    }  else    {      if (a < 0.0)	{	  GSL_SET_COMPLEX (&z, M_PI, -acosh (-a));	}      else	{	  GSL_SET_COMPLEX (&z, 0, acosh (a));	}    }  return z;}

gsl_complexgsl_complex_arccsc_real (double a){				/* z = arccsc(a) */  gsl_complex z;  if (a <= -1.0 || a >= 1.0)    {      GSL_SET_COMPLEX (&z, asin (1 / a), 0.0);    }  else    {      if (a >= 0.0)	{	  GSL_SET_COMPLEX (&z, M_PI_2, -acosh (1 / a));	}      else	{	  GSL_SET_COMPLEX (&z, -M_PI_2, -acosh (-1 / a));	}    }  return z;}

gsl_complex gsl_complex_tan(gsl_complex a){				/* z = tan(a) */    double R = GSL_REAL(a), I = GSL_IMAG(a);    gsl_complex z;    if (fabs(I) < 1) {	double D = pow(cos(R), 2.0) + pow(sinh(I), 2.0);	GSL_SET_COMPLEX(&z, 0.5 * sin(2 * R) / D, 0.5 * sinh(2 * I) / D);    } else {	double u = exp(-I);	double C = 2 * u / (1 - pow(u, 2.0));	double D = 1 + pow(cos(R), 2.0) * pow(C, 2.0);	double S = pow(C, 2.0);	double T = 1.0 / tanh(I);	GSL_SET_COMPLEX(&z, 0.5 * sin(2 * R) * S / D, T / D);    }    return z;}

gsl_complexgsl_complex_tanh (gsl_complex a){				/* z = tanh(a) */  double R = GSL_REAL (a), I = GSL_IMAG (a);  gsl_complex z;  if (fabs(R) < 1.0)     {      double D = pow (cos (I), 2.0) + pow (sinh (R), 2.0);            GSL_SET_COMPLEX (&z, sinh (R) * cosh (R) / D, 0.5 * sin (2 * I) / D);    }  else    {      double D = pow (cos (I), 2.0) + pow (sinh (R), 2.0);      double F = 1 + pow (cos (I) / sinh (R), 2.0);      GSL_SET_COMPLEX (&z, 1.0 / (tanh (R) * F), 0.5 * sin (2 * I) / D);    }  return z;}

gsl_complexgsl_complex_arccosh_real (double a){				/* z = arccosh(a) */  gsl_complex z;  if (a >= 1)    {      GSL_SET_COMPLEX (&z, acosh (a), 0);    }  else    {      if (a >= -1.0)	{	  GSL_SET_COMPLEX (&z, 0, acos (a));	}      else	{	  GSL_SET_COMPLEX (&z, acosh (-a), M_PI);	}    }  return z;}

symrec *putsym (char *sym_name, symrec_type sym_type){	symrec *ptr;	ptr = (symrec *)malloc(sizeof(symrec));	/* names are always lowercase */	ptr->name = strdup(sym_name);	str_tolower(ptr->name);	ptr->type = sym_type;	GSL_SET_COMPLEX(&ptr->value.c, 0, 0); /* set value to 0 even if fctn.  */	ptr->next = (struct symrec *)sym_table;	sym_table = ptr;	return ptr;}

void sym_init_table ()  /* puts arithmetic functions in table. */{	int i;	symrec *ptr;	for (i = 0; arith_fncts[i].fname != 0; i++){		ptr = putsym (arith_fncts[i].fname, S_FNCT);		ptr->value.fnctptr = arith_fncts[i].fnct;	}	/* now the constants */	for (i = 0; arith_cnts[i].fname != 0; i++){		ptr = putsym(arith_cnts[i].fname, S_CMPLX);		GSL_SET_COMPLEX(&ptr->value.c, arith_cnts[i].c, 0);	}}

gsl_complex gsl_complex_div(gsl_complex a, gsl_complex b){				/* z=a/b */    double ar = GSL_REAL(a), ai = GSL_IMAG(a);    double br = GSL_REAL(b), bi = GSL_IMAG(b);    double s = 1.0 / gsl_complex_abs(b);    double sbr = s * br;    double sbi = s * bi;    double zr = (ar * sbr + ai * sbi) * s;    double zi = (ai * sbr - ar * sbi) * s;    gsl_complex z;    GSL_SET_COMPLEX(&z, zr, zi);    return z;}

static intmc_conj(lua_State *L){    mMatComplex *a = qlua_checkMatComplex(L, 1);    mMatComplex *b = qlua_newMatComplex(L, a->l_size, a->r_size);    int l, r;    for (l = 0; l < b->l_size; l++) {        for (r = 0; r < b->r_size; r++) {            gsl_complex q;            gsl_complex z = gsl_matrix_complex_get(a->m, l, r);            GSL_SET_COMPLEX(&q, GSL_REAL(z), -GSL_IMAG(z));            gsl_matrix_complex_set(b->m, l, r, q);        }    }    return 1;}