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

static voidcalc_primes(char *bitmap, uint64_t max){   uint64_t a,            r;   /*** Main algorithm start ***/   for (a = 3; a <= sqrtl(max); a += 2) {      if (bitmap[a/2/8] & 1 << ((a/2) % 8))         continue;      for (r = a * a; r <= max; r += a * 2)         bitmap[r/2/8] |= 1 << ((r/2) % 8);   }   /*** Main algorithm end ***/   /* mark 1 as not prime */   bitmap[0] |= 0x1;}

long doublehw2p7a (long double h, long double p){  // trough  // h is height of water from top of trough  // volume is 12.4 feet.  // Got the formula for the water volume by height from here:  // http://answers.yahoo.com/question/index?qid=20120329114323AAgspZg  long double r = 1.0L;		// radius  long double l = 10.0L;	// length  long double maxVol = M_PI * r * r * 0.5L * l;  h = 1-h;  long double result =  l*( r*r * ( M_PI/2 -asinl((1-h)/r) ) - (r-h) * sqrtl( h*2*r - h*h) ) ;  //long double result = maxVol - l *(  r*r* acosl((r-h)/r) - (r-h)*powl(2*r*h - h*h, 0.5L) ) ;  //long double result = maxVol - l *( r*r * asinl(h/r - h* sqrtl(r*r - h*h))  ) ;  printf("a h %Lf r %Lf l %Lf max %Lf result %Lf p %Lf/n", h, r, l, maxVol, result, p );  return result-p;}

compl compl_create(const real rl,const real img){    /* DESC: Constructor del tipo */    compl c;     c.mod = sqrtl(rl*rl + img*img);    if (rl > 0){      c.ang = (atanl(img/rl));    }else if ((rl < 0) && (img >= 0)) {            c.ang = atanl((img/rl)+ PI);    }else if ((rl < 0) && (img < 0)) {      c.ang = (atanl (img/rl)) - PI;    }else if ((rl == 0) &&(img > 0)){      c.ang = PI / 2;    }else if ((rl == 0) && (img < 0)){      c.ang = -PI/2;    }else if ((rl == 0) && (img == 0)){      c.ang = 0;    }                                      return c;}

intmain(int argc, char *argv[]){   uint64_t  max = argc > 1 ? atol(argv[1]) : 1000000000ull;   char     *bitmap = calloc(max/30 + sqrtl(max) + 16*8, 1);   uint64_t  adj_max;   int       a_i;   /* adjust max to a multiple of 2,3,5 */   for (adj_max = (max - 1) / 30 * 30 + 1, a_i = 0;        adj_max + diff[a_i%8] <= max;        adj_max += diff[a_i++%8])      ;   calc_primes(bitmap, adj_max);   printf("%ld/n", count_primes(bitmap, adj_max / 30 * 8 + num2bit(adj_max)) + 3);   free(bitmap);   return EXIT_SUCCESS;}

static voidcalc_primes(char *bitmap, uint64_t max){   uint64_t a,            r;   /*** Main algorithm start ***/   for (a = 2; a <= sqrtl(max); a++) {      if (bitmap[a / 8] & 1 << (a % 8))         continue;      for (r = a * a; r <= max; r += a)         bitmap[r / 8] |= 1 << (r % 8);   }   /*** Main algorithm end ***/   /* mark bit 0 and 1 as not prime */   bitmap[0] |= 0x3;}

int main(){	double sp, arie, a, b, c, mid;	scanf("%d %d %d %d %d %d %d", &A.x, &A.y, &B.x, &B.y, &C.x, &C.y, &L);	a = dist(B, C);	b = dist(A, C);	c = dist(A, B);		Mid.x = (A.x + B.x) / 2;	Mid.y = (A.y + B.y) / 2;	mid = dist(C, Mid);	sp = (a + b + c) / 2;	arie = sqrtl(sp * (sp - a) * (sp - b) * (sp - c));	if (c == 0)		Dif = b;	else	if (mid >= min(a, b) && mid <= max(a, b))		Dif = min(a, b);	else		Dif = 2 * arie / c;	Dif = (Dif < L) ? 0 : Dif - L;	LMax = max(a, b);	LMax = (LMax < L) ? 0 : LMax - L;	printf("%.2lf/n%.2lf/n", Dif, LMax);	return 0;}

  void plot_loops (polygon_t * p, poly_dim_t speed, double flowrate, int dir)  {    if (!p)      return;    poly_contour_t *c;    for (c = p->contours; c; c = c->next)      if (c->vertices && (!dir || dir == c->dir))	{	  poly_vertex_t *v = c->vertices;	  poly_dim_t d = sqrtl ((px - v->x) * (px - v->x) + (py - v->y) * (py - v->y));	  if (pe && d > layer * 5)	    {			// hop and pull back extruder while moving	      move (px, py, z + hop, back);	      move (v->x, v->y, z + hop, back);	    }	  move (v->x, v->y, z, 0);	  for (v = c->vertices->next; v; v = v->next)	    extrude (v->x, v->y, z, speed, flowrate);	  v = c->vertices;	  extrude (v->x, v->y, z, speed, flowrate);	}  }

voidSqrt (Token ** Pila){  Token *Operando = EntornoEjecucion_BuscaSimbolo (*Pila);  if (Buzon.GetHuboError ())    return;  if (NoEsReal (Operando))    {      BorrarTokenSiEsVariable (Operando);      return;    }  long double ValorRetorno = sqrtl (Operando->GetDatoReal ());  BorrarTokenSiEsVariable (Operando);  if (Buzon.GetHuboError ())    return;  Token *TokenRetorno = ConsigueToken (ValorRetorno);  delete Desapila (Pila);  Apila (Pila, TokenRetorno);  return;}

int main(void){	pthread_attr_t attr;	pthread_attr_init(&attr);	pthread_attr_setstacksize(&attr, 65536);		uint64_t cpus = 6;		pthread_t threads[cpus];	doRangeArg arg[cpus];		uint64_t max = sqrtl(MAX);	uint64_t *fsms = malloc(max*sizeof(uint64_t));	for(uint64_t i = 0; i < max; ++i){		fsms[i] = 1;	}	for(uint64_t p = 2; p < max; ++p){		if(fsms[p] != 1){//p is not prime			continue;		}//p is prime, mark off its multiples		for(uint64_t m = p; m < max; m += p){			mark(fsms, p, m);		}	}		for(uint64_t i = 0; i < cpus; ++i){		arg[i] = (doRangeArg){.start = i + 1, .step = cpus, .fsms = fsms, .max = max};		if(pthread_create(&threads[i], &attr, doRange, &arg[i])){			return 1;		}	}	for(uint64_t i = 0; i < cpus; ++i){		if(pthread_join(threads[i], NULL)){			return 1;		}	}		free(fsms);	printf("%"PRIu64"/n", s);}

int main(){	long long num,a,b,count,c=1,i,temp;	scanf("%lld",&num);	while(num--)	{		count=0;		scanf("%lld%lld",&a,&b);		for(i=a;i<=b;i++)		{			if(check(i))			{			temp=sqrtl(i);			if(temp*temp==i)				if(check(temp))					count++;			}		}		printf("Case #%lld: %lld/n",c++,count);	}	return 0 ;}

void prepareColCor(double * x, size_t nr, int cosine, double * res, size_t * nNAentries, int * NAmean){  *nNAentries = 0;  size_t count = 0;  LDOUBLE mean = 0, sum = 0;  for (size_t k = 0; k < nr; k++)    if (!ISNAN(x[k]))    {      mean += x[k];      sum += ((LDOUBLE) x[k])*( (LDOUBLE) x[k]);      count ++;    }  if (count > 0)  {    *NAmean = 0;    *nNAentries = nr-count;    if (cosine) mean = 0; else mean = mean/count;    sum = sqrtl(sum - count * mean*mean);    if (sum > 0)    {      // Rprintf("sum: %Le/n", sum);       for (size_t k=0; k<nr; k++)         if (!ISNAN(x[k]))            res[k] = (x[k] - mean)/sum;         else            res[k] = 0;    } else {       // Rprintf("prepareColCor: have zero variance./n");       *NAmean = 1;       for (size_t k=0; k<nr; k++) res[k] = 0;    }  } else {    *NAmean = 1;    *nNAentries = nr;    for (size_t k=0; k<nr; k++)       res[k] = 0;  }}

void fun(long N){	   double *a, *b, *c;   long i;   double t_start, t_end;     a = dvector(N);   b = dvector(N);   c = dvector(N);   if ( !a || !b || !c ) {       printf("Not enough memory!/n");       exit(2);   }   for (i=0; i < N; i++)      a[i] = sqrtl(i);   for (i=0; i < N; i++)      b[i] = 1.0;   t_start = omp_get_wtime();     // #pragma omp parallel for	for(i=0;i<N;i++){		c[i] = a[i]*b[i];	}      t_end = omp_get_wtime();      printf("%10ld %20.15lf/n", N, t_end - t_start);     free_dvector(c);   free_dvector(b);   free_dvector(a);	}

long double asinl(long double x){	union ldshape u = {x};	long double z, r, s;	uint16_t e = u.i.se & 0x7fff;	int sign = u.i.se >> 15;	if (e >= 0x3fff) {   /* |x| >= 1 or nan */		/* asin(+-1)=+-pi/2 with inexact */		if (x == 1 || x == -1)			return x*pio2_hi + 0x1p-120f;		return 0/(x-x);	}	if (e < 0x3fff - 1) {  /* |x| < 0.5 */		if (e < 0x3fff - (LDBL_MANT_DIG+1)/2) {			/* return x with inexact if x!=0 */			FORCE_EVAL(x + 0x1p120f);			return x;		}		return x + x*__invtrigl_R(x*x);	}	/* 1 > |x| >= 0.5 */	z = (1.0 - fabsl(x))*0.5;	s = sqrtl(z);	r = __invtrigl_R(z);	if (CLOSETO1(u)) {		x = pio2_hi - (2*(s+s*r)-pio2_lo);	} else {		long double f, c;		u.f = s;		CLEARBOTTOM(u);		f = u.f;		c = (z - f*f)/(s + f);		x = 0.5*pio2_hi-(2*s*r - (pio2_lo-2*c) - (0.5*pio2_hi-2*f));	}	return sign ? -x : x;}

int main(int argc, char **argv) {    bitmap *sieve = bitmap_new(0, SIEVEFILE);    if (!sieve){        perror(SIEVEFILE);        exit(errno);    }    U64 start = argc > 1 ? atoll(argv[1]) : 0;    U64 size  = argc > 2 ? atoll(argv[2]) : 1000*1000*1000;    bitmap *b = bitmap_new(size >> 1, NULL);    bitmap_fill(b, 1);    if (start == 0) bitmap_set(b, 0, 0);    U64 pmax = (U64)sqrtl(start+size);#ifdef VERBOSE    printf("pmax = %llu/n", pmax);#endif    U64 p, i;    for (p = 3; p < pmax;){#ifdef VERBOSE        printf("sieving %llu/r", p);        fflush(stdout);#endif        for ( i = p + p - (start % p) ; i <= size ; i += p ) {            if ((i & 1) == 0) continue; /* skip even */            bitmap_set( b, i>>1, 0 );        }        for(p += 2; !bitmap_get(sieve, p>>1); p += 2); /* seek next prime */    }    char filename[256];    snprintf(filename, 256, "%llu~%llu.bm", start, start+size);#ifdef VERBOSE    printf("saving %s/n", filename);#endif    bitmap_save(b, filename);    return 0;}

static voidcalc_primes(char *bitmap, uint64_t max){   char     *bmp,            *bitmap_end = bitmap + 32*1024;   uint64_t  a,             cur;   int       i,             a_i,             offsets[8];   /*** Main algorithm start ***/   for (cur = 0; cur <= max; cur += 32*1024*30, bitmap_end += 32*1024) {      for (a = 7, a_i = 1; a <= sqrtl(cur + 32*1024*30); a += diff[a_i++%8]) {         if (bitmap[a / 30] & 1 << (a_i%8))            continue;         for (i = 0; i < 8; i++)            offsets[num2bit(a * bval[i])] = a * bval[i] / 30;         for (bmp = bitmap + MAX(a * (a / 30), (cur / 30) / a * a);              bmp < bitmap_end;              bmp += a) {            for (i = 0; i < 8; i++)               *(bmp + offsets[i]) |= 1 << i;         }      }      /* First byte gets clobbered for the first block */      if (cur == 0)         bitmap[0] = 0x1;   }   /*** Main algorithm end ***/}

int main(void){	int cases;	scanf("%d",&cases);	double a,b,c;	double t3,tn,s;	long long xt3,xtn,xs;	double temp1,temp2;	long long n,first;	long long i,d;	while(cases--)	{		scanf("%lld %lld %lld",&xt3,&xtn,&xs);		t3=xt3;		tn=xtn;		s=xs;		b=(-1*((7*t3)+(5*tn)+(2*s)));		a=(t3+tn);		c=(12*s);		temp1=((b*b)-(4*a*c));		temp1=sqrtl(temp1);		temp2=((-1*b)+temp1)/(2*a);		n=floor(temp2);		printf("%lld/n",n);		d=((xtn-xt3)/(n-6));		first=xt3-(2*d);		printf("%lld",first);		for(i=1;i<n;i++)		{			first+=d;			printf(" %lld",first);		}		printf("/n");	}	return 0;}

int main(){int i,n;double d,l,area;double a,b,pi;pi=2*acos(0);scanf("%d",&n);for(i=0;i<n;i++){   scanf("%lf %lf",&d,&l);        a=l/2;    b=sqrtl(((l/2)*(l/2))-((d/2)*(d/2)));    area=pi*a*b;    printf("%.3lf/n",area);    }return 0;}

int main(){	int test_cases;	scanf("%d",&test_cases);	while(test_cases--)	{		long long int x,y,z,done=0,c,b,ans,root,d,a1,a2,ans1,ans2;		scanf("%lld%lld%lld",&x,&y,&z);		d = x*x+y*y+2*y*z;		root = sqrtl(d);		a1 = x + root;		a2 = x - root;		c=y+2*z;		ans1 = (a1*z)/c;		ans2 = (a2*z)/c;		if(root*root == d && (a1*z)%c==0)			printf("%lld/n",ans1);		else if(root*root == d && (a2*z)%c==0)			printf("%lld/n",ans2);		else			printf("Not this time./n");	}	return 0;}

floatx80 floatx80_sqrt( floatx80 a STATUS_PARAM){    return sqrtl(a);}

ldcomplexcasinl(ldcomplex z) {	long double x, y, t, R, S, A, Am1, B, y2, xm1, xp1, Apx;	int ix, iy, hx, hy;	ldcomplex ans;	x = LD_RE(z);	y = LD_IM(z);	hx = HI_XWORD(x);	hy = HI_XWORD(y);	ix = hx & 0x7fffffff;	iy = hy & 0x7fffffff;	x = fabsl(x);	y = fabsl(y);	/* special cases */	/* x is inf or NaN */	if (ix >= 0x7fff0000) {	/* x is inf or NaN */		if (isinfl(x)) {	/* x is INF */			LD_IM(ans) = x;			if (iy >= 0x7fff0000) {				if (isinfl(y))					/* casin(inf + i inf) = pi/4 + i inf */					LD_RE(ans) = pi_4 + pi_4_l;				else	/* casin(inf + i NaN) = NaN + i inf */					LD_RE(ans) = y + y;			} else	/* casin(inf + iy) = pi/2 + i inf */				LD_RE(ans) = pi_2 + pi_2_l;		} else {		/* x is NaN */			if (iy >= 0x7fff0000) {				/* INDENT OFF */				/*				 * casin(NaN + i inf) = NaN  + i inf				 * casin(NaN + i NaN) = NaN  + i NaN				 */				/* INDENT ON */				LD_IM(ans) = y + y;				LD_RE(ans) = x + x;			} else {				/* INDENT OFF */				/* casin(NaN + i y  ) = NaN  + i NaN */				/* INDENT ON */				LD_IM(ans) = LD_RE(ans) = x + y;			}		}		if (hx < 0)			LD_RE(ans) = -LD_RE(ans);		if (hy < 0)			LD_IM(ans) = -LD_IM(ans);		return (ans);	}	/* casin(+0 + i 0) = 0 + i 0. */	if (x == zero && y == zero)		return (z);	if (iy >= 0x7fff0000) {	/* y is inf or NaN */		if (isinfl(y)) {	/* casin( x + i inf ) = 0 + i inf */			LD_IM(ans) = y;			LD_RE(ans) = zero;		} else {		/* casin( x + i NaN ) = NaN + i NaN */			LD_IM(ans) = x + y;			if (x == zero)				LD_RE(ans) = x;			else				LD_RE(ans) = y;		}		if (hx < 0)			LD_RE(ans) = -LD_RE(ans);		if (hy < 0)			LD_IM(ans) = -LD_IM(ans);		return (ans);	}	if (y == zero) {	/* region 1: y=0 */		if (ix < 0x3fff0000) {	/* |x| < 1 */			LD_RE(ans) = asinl(x);			LD_IM(ans) = zero;		} else {			LD_RE(ans) = pi_2 + pi_2_l;			if (ix >= ip1)	/* |x| >= i386 ? 2**65 : 2**114 */				LD_IM(ans) = ln2 + logl(x);			else if (ix >= 0x3fff8000)	/* x > Acrossover */				LD_IM(ans) = logl(x + sqrtl((x - one) * (x +					one)));			else {				xm1 = x - one;				LD_IM(ans) = log1pl(xm1 + sqrtl(xm1 * (x +					one)));			}		}	} else if (y <= E * fabsl(x - one)) {	/* region 2: y < tiny*|x-1| */		if (ix < 0x3fff0000) {	/* x < 1 */			LD_RE(ans) = asinl(x);			LD_IM(ans) = y / sqrtl((one + x) * (one - x));		} else {			LD_RE(ans) = pi_2 + pi_2_l;			if (ix >= ip1)	/* i386 ? 2**65 : 2**114 */				LD_IM(ans) = ln2 + logl(x);//.........这里部分代码省略.........

long doublehw2p1 (long double x, long double p){  return sqrtl (x) - cosl (x);}

spotsPump * acopleOptico(long double wFuenteLaserTan, long double wFuenteLaserSag, long double divergencia, long double lambda0,/                         long double La, long double tL, long double R11, long double R12,/                         long double nL, long double thetaL, long double Lb, long double tE,/                         long double R21, long double R22, long double nE, long double thetaE,/                          long double nC, long double delta1, matriz * epsilon1)                        {                            long double spotWaistTan=2*lambda0/M_PI/divergencia;                            long double spotWaistSag=2*lambda0/M_PI/divergencia;                            long double LtoInTan=sqrtl(powl(wFuenteLaserTan/spotWaistTan,2)-1)*(M_PI*powl(spotWaistTan,2))/lambda0;                            long double LtoInSag=sqrtl(powl(wFuenteLaserSag/spotWaistSag,2)-1)*(M_PI*powl(spotWaistSag,2))/lambda0;                            long double _Complex qInTan=La+LtoInTan+I*M_PI*1*powl(spotWaistTan,2)/lambda0;                            long double _Complex qInSag=La+LtoInSag+I*M_PI*1*powl(spotWaistSag,2)/lambda0;                            assert(epsilon1->filas==1);                            long double Lc;                            // Prepara estructura de salida                            spotsPump *Salida=nuevoSpotsPump(epsilon1->columnas);                            // Variables de trabajo                            long double tanA11, tanB11, tanC11, tanD11, tanA12, tanB12, tanC12, tanD12,/                             tanA21, tanB21, tanC21, tanD21, tanA22, tanB22, tanC22, tanD22;                            // Entrada y salida de lentes gruesas - tangencial                            // Lente 1                            tanA11=sqrtl(1-powl(sinl(thetaL)/nL,2));                            tanB11=0;                            tanC11=(sqrtl(1-powl(sinl(thetaL)/nL,2))-nL*cosl(thetaL))/(R11*cosl(thetaL)*sqrtl(1-powl(sinl(thetaL)/nL,2)));                            tanD11=1/tanA11;                            tanA12=1/tanA11;                            tanB12=0;                            tanC12=-tanC11;                            tanD12=tanA11;                            //Lente 2 (Espejo cóncavo) -tangencial                            tanA21=sqrtl(1-powl(sinl(thetaL)/nL,2));                            tanB21=0;                            tanC21=(sqrtl(1-powl(sinl(thetaL)/nL,2))-nL*cosl(thetaL))/(R12*cosl(thetaL)*sqrtl(1-powl(sinl(thetaL)/nL,2)));                            tanD21=1/tanA11;                            tanA22=1/tanA11;                            tanB22=0;                            tanC22=-tanC11;                            tanD22=tanA11;                            long double sagA11, sagB11, sagC11, sagD11, sagA12, sagB12, sagC12, sagD12,/                             sagA21, sagB21, sagC21, sagD21, sagA22, sagB22, sagC22, sagD22;                            // Entrada y salida de lentes gruesas - sagital                            // Lente 1                            sagA11=sqrtl(1-powl(sinl(thetaE)/nE,2));                            sagB11=0;                            sagC11=(sqrtl(1-powl(sinl(thetaE)/nE,2))-nE*cosl(thetaE))/(R21*cosl(thetaE)*sqrtl(1-powl(sinl(thetaE)/nE,2)));                            sagD11=1/sagA11;                            sagA12=1/sagA11;                            sagB12=0;                            sagC12=-sagC11;                            sagD12=sagA11;                            //Lente 2 (Espejo cóncavo) - sagital                            sagA21=sqrtl(1-powl(sinl(thetaE)/nE,2));                            sagB21=0;                            sagC21=(sqrtl(1-powl(sinl(thetaE)/nE,2))-nE*cosl(thetaE))/(R22*cosl(thetaE)*sqrtl(1-powl(sinl(thetaE)/nE,2)));                            sagD21=1/sagA11;                            sagA22=1/sagA11;                            sagB22=0;                            sagC22=-sagC11;                            sagD22=sagA11;                            // Declarando matrices                            matriz *Tan1, *Tan2, *Tan3, *Tan4, *Tan5, *Tan6, *Tan7,/                             *Tan8, *Tan9, *Tan10;                            matriz *Sag1, *Sag2, *Sag3, *Sag4, *Sag5, *Sag6, *Sag7,/                             *Sag8, *Sag9, *Sag10;                            // Plano tangencial                            Tan1=llenaMatriz(1, La, 0, 1);                            Tan2=llenaMatriz(tanA11,tanB11,tanC11,tanD11);                            Tan3=llenaMatriz(1,tL/nL,0,1);                            Tan4=llenaMatriz(tanA12,tanB12,tanC12,tanD12);                            Tan5=llenaMatriz(1,Lb,0,1);                            Tan6=llenaMatriz(tanA21,tanB21,tanC21,tanD21);                            Tan7=llenaMatriz(1,tE/nE,0,1);                            Tan8=llenaMatriz(tanA22,tanB22,tanC22,tanD22);                            Tan9=nuevaMatriz(2,2); // Variable                            Tan10=llenaMatriz(nC,0,0,1/nC);                            // Plano tangencial                            Sag1=llenaMatriz(1, La, 0, 1);                            Sag2=llenaMatriz(sagA11,sagB11,sagC11,sagD11);                            Sag3=llenaMatriz(1,tL/nL,0,1);                            Sag4=llenaMatriz(sagA12,sagB12,sagC12,sagD12);                            Sag5=llenaMatriz(1,Lb,0,1);                            Sag6=llenaMatriz(sagA21,sagB21,sagC21,sagD21);                            Sag7=llenaMatriz(1,tE/nE,0,1);                            Sag8=llenaMatriz(sagA22,sagB22,sagC22,sagD22);                            Sag9=nuevaMatriz(2,2); // Variable                            Sag10=llenaMatriz(1,0,0,1);//.........这里部分代码省略.........

long doublehypotl(long double x, long double y){	long double a=x,b=y,t1,t2,y1,y2,w;	int32_t j,k,ha,hb;	GET_HIGH_WORD(ha,x);	ha &= 0x7fff;	GET_HIGH_WORD(hb,y);	hb &= 0x7fff;	if(hb > ha) {a=y;b=x;j=ha; ha=hb;hb=j;} else {a=x;b=y;}	a = fabsl(a);	b = fabsl(b);	if((ha-hb)>DESW(MANT_DIG+7)) {return a+b;} /* x/y > 2**(MANT_DIG+7) */	k=0;	if(ha > ESW(MAX_EXP/2-12)) {	/* a>2**(MAX_EXP/2-12) */	   if(ha >= ESW(MAX_EXP)) {	/* Inf or NaN */	       man_t manh, manl;	       /* Use original arg order iff result is NaN; quieten sNaNs. */	       w = fabsl(x+0.0)-fabsl(y+0.0);	       GET_LDBL_MAN(manh,manl,a);	       if (manh == LDBL_NBIT && manl == 0) w = a;	       GET_LDBL_MAN(manh,manl,b);	       if (hb >= ESW(MAX_EXP) && manh == LDBL_NBIT && manl == 0) w = b;	       return w;	   }	   /* scale a and b by 2**-(MAX_EXP/2+88) */	   ha -= DESW(MAX_EXP/2+88); hb -= DESW(MAX_EXP/2+88);	   k += MAX_EXP/2+88;	   SET_HIGH_WORD(a,ha);	   SET_HIGH_WORD(b,hb);	}	if(hb < ESW(-(MAX_EXP/2-12))) {	/* b < 2**-(MAX_EXP/2-12) */	    if(hb <= 0) {		/* subnormal b or 0 */	        man_t manh, manl;		GET_LDBL_MAN(manh,manl,b);		if((manh|manl)==0) return a;		t1=0;		SET_HIGH_WORD(t1,ESW(MAX_EXP-2));	/* t1=2^(MAX_EXP-2) */		b *= t1;		a *= t1;		k -= MAX_EXP-2;	    } else {		/* scale a and b by 2^(MAX_EXP/2+88) */		ha += DESW(MAX_EXP/2+88);		hb += DESW(MAX_EXP/2+88);		k -= MAX_EXP/2+88;		SET_HIGH_WORD(a,ha);		SET_HIGH_WORD(b,hb);	    }	}    /* medium size a and b */	w = a-b;	if (w>b) {	    t1 = a;	    union IEEEl2bits uv;	    uv.e = t1; uv.bits.manl = 0; t1 = uv.e;	    t2 = a-t1;	    w  = sqrtl(t1*t1-(b*(-b)-t2*(a+t1)));	} else {	    a  = a+a;	    y1 = b;	    union IEEEl2bits uv;	    uv.e = y1; uv.bits.manl = 0; y1 = uv.e;	    y2 = b - y1;	    t1 = a;	    uv.e = t1; uv.bits.manl = 0; t1 = uv.e;	    t2 = a - t1;	    w  = sqrtl(t1*y1-(w*(-w)-(t1*y2+t2*b)));	}	if(k!=0) {	    u_int32_t high;	    t1 = 1.0;	    GET_HIGH_WORD(high,t1);	    SET_HIGH_WORD(t1,high+DESW(k));	    return t1*w;	} else return w;}

voidfun_C(int x, int y, int z){    b=sqrtl(a);    printf("%Lf/n", b);}

// sqrt(r)// uses bntmp1 - bntmp6 - global temp bignumbers//  SIDE-EFFECTS://      n ends up as |n|bn_t sqrt_bn(bn_t r, bn_t n){    int comp, almost_match = 0;    LDBL f;    bn_t orig_r, orig_n;    int  orig_bnlength,         orig_padding,         orig_rlength,         orig_shiftfactor;    // use Newton's recursive method for zeroing in on sqrt(n): r=.5(r+n/r)    if (is_bn_neg(n))    {   // sqrt of a neg, return 0        clear_bn(r);        return r;    }    f = bntofloat(n);    if (f == 0) // division by zero will occur    {        clear_bn(r); // sqrt(0) = 0        return r;    }    f = sqrtl(f); // approximate square root    // no need to check overflow    // With Newton's Method, there is no need to calculate all the digits    // every time.  The precision approximately doubles each iteration.    // Save original values.    orig_bnlength      = bnlength;    orig_padding       = padding;    orig_rlength       = rlength;    orig_shiftfactor   = shiftfactor;    orig_r             = r;    orig_n             = n;    // calculate new starting values    bnlength = intlength + (int)(LDBL_DIG/LOG10_256) + 1; // round up    if (bnlength > orig_bnlength)        bnlength = orig_bnlength;    calc_lengths();    // adjust pointers    r = orig_r + orig_bnlength - bnlength;    floattobn(r, f); // start with approximate sqrt    copy_bn(bntmp4, r);    for (int i = 0; i < 25; i++) // safety net, this shouldn't ever be needed    {        // adjust lengths        bnlength <<= 1; // double precision        if (bnlength > orig_bnlength)            bnlength = orig_bnlength;        calc_lengths();        r = orig_r + orig_bnlength - bnlength;        n = orig_n + orig_bnlength - bnlength;        copy_bn(bntmp6, r);        copy_bn(bntmp5, n);        unsafe_div_bn(bntmp4, bntmp5, bntmp6);        add_a_bn(r, bntmp4);        half_a_bn(r);        if (bnlength == orig_bnlength && (comp = abs(cmp_bn(r, bntmp4))) < 8)  // if match or almost match        {            if (comp < 4  // perfect or near perfect match                    || almost_match == 1) // close enough for 2nd time                break;            else // this is the first time they almost matched                almost_match++;        }    }    // restore original values    bnlength      = orig_bnlength;    padding       = orig_padding;    rlength       = orig_rlength;    shiftfactor   = orig_shiftfactor;    r             = orig_r;    return r;}

voidCheby2D_Test(Cheby2D *cheby, int nx_test, int ny_test,        void        (*func)(long double *x, long double *y,             int nx, int ny, long double *z, void *info),        void *info,        long double *residualRMS, long double *residualMAV){    int ix_test, iy_test;    long double  fprime;    long double sdiff, ssdiff, maxdiff, diff;    long double *f = (long double*)malloc(nx_test*ny_test*sizeof(long double));    long double *x = (long double *)malloc(nx_test*sizeof(long double));    long double *y = (long double *)malloc(ny_test*sizeof(long double));    sdiff=ssdiff=maxdiff=0.0L;    // call func to get actual values on grid    for (iy_test=0; iy_test<ny_test; iy_test++)        y[iy_test] = -1.0L + 2.0L*(iy_test+0.5L)/ny_test;     for (ix_test=0; ix_test<nx_test; ix_test++)        x[ix_test] = -1.0L + 2.0L*(ix_test+0.5L)/nx_test;     func(x, y, nx_test, ny_test, f, info);    // compare    for (iy_test=0; iy_test<ny_test; iy_test++)    {        for (ix_test=0; ix_test<nx_test; ix_test++)        {            fprime = Cheby2D_Evaluate(cheby, x[ix_test], y[iy_test]);            diff = f[iy_test*nx_test+ix_test]-fprime;#if 0            {                // redwards hack to simulate 1-D polyco.. first get prediction from                // centre freq                long double phase1d = Cheby2D_Evaluate(cheby, x[ix_test], 0.0);                // second, get instantaneous freq at band center                long double dx = 1.0e-5L;                long double freq1 =                     (Cheby2D_Evaluate(cheby, x[ix_test]+dx*0.5L, y[iy_test])-                     Cheby2D_Evaluate(cheby, x[ix_test]-dx*0.5L, y[iy_test])) / dx;                long double freq=                    (Cheby2D_Evaluate(cheby, x[ix_test]+dx*0.5L, 0.0L)-                     Cheby2D_Evaluate(cheby, x[ix_test]-dx*0.5L, 0.0L)) / dx;                long double skyfreq = 653.0L+y[iy_test]*32.0L; // MHzx                long double dmdelay = 1.0L/2.41e-4L *(1.0L/ (skyfreq*skyfreq) -1.0L/(653.0L*653.0L))                    *48.901787L; // s                long double dphase = -freq*dmdelay;                long double diff1d = f[iy_test*nx_test+ix_test]-(phase1d+dphase);                freq *= 2.0L/0.03L / 86400.0L; // normalized->day^-1->Hz                freq1 *= 2.0L/0.03L / 86400.0L;                //             printf("YYY %Lf %Lf %Lf/n", freq, dmdelay, dphase);                //       printf("%Lf %Lf %Lg %Lg %Lg %Lg %Lg XXX/n", x[ix_test], y[iy_test],                 //   	     f[iy_test*nx_test+ix_test], fprime, diff, diff1d,                 // 	     (freq1-freq)/freq);            }#endif            if (fabs(diff) > maxdiff)                maxdiff = fabs(diff);            sdiff += diff;            ssdiff += diff*diff;        }        //      printf("/n");    }    sdiff /=nx_test*ny_test;    *residualRMS = sqrtl(ssdiff/(nx_test*ny_test));    *residualMAV = maxdiff;    free(f);    free(x);    free(y);}

/*Function for the gravitational potential to be used. Potential for Plummer model.*/mydbl potential(mydbl r){  return -G*M/(sqrtl(A*A + r*r));}

/*$cp^{0}$ multiplied by this factor allows to obtain the energy for a local inertial observer in this spacetime.*/mydbl energy_factor(mydbl r){  mydbl g = sqrtl(1.0 + 2.0*potential(r)/(C*C));  return g;}

long double complexcsqrtl(long double complex z){	long double complex result;	long double a, b;	long double t;	int scale;	a = creall(z);	b = cimagl(z);	/* Handle special cases. */	if (z == 0)		return (cpackl(0, b));	if (isinf(b))		return (cpackl(INFINITY, b));	if (isnan(a)) {		t = (b - b) / (b - b);	/* raise invalid if b is not a NaN */		return (cpackl(a, t));	/* return NaN + NaN i */	}	if (isinf(a)) {		/*		 * csqrt(inf + NaN i)  = inf +  NaN i		 * csqrt(inf + y i)    = inf +  0 i		 * csqrt(-inf + NaN i) = NaN +- inf i		 * csqrt(-inf + y i)   = 0   +  inf i		 */		if (signbit(a))			return (cpackl(fabsl(b - b), copysignl(a, b)));		else			return (cpackl(a, copysignl(b - b, b)));	}	/*	 * The remaining special case (b is NaN) is handled just fine by	 * the normal code path below.	 */	/* Scale to avoid overflow. */	if (fabsl(a) >= THRESH || fabsl(b) >= THRESH) {		a *= 0.25;		b *= 0.25;		scale = 1;	} else {		scale = 0;	}	/* Algorithm 312, CACM vol 10, Oct 1967. */	if (a >= 0) {		t = sqrtl((a + hypotl(a, b)) * 0.5);		result = cpackl(t, b / (2 * t));	} else {		t = sqrtl((-a + hypotl(a, b)) * 0.5);		result = cpackl(fabsl(b) / (2 * t), copysignl(t, b));	}	/* Rescale. */	if (scale)		return (result * 2);	else		return (result);}

/*To set the initial value of pr, it must hold $g_{/mu/nu}p^{/mu}p^{/nu} = 0$.This factor multiplies p0 to guarantee that p1 fulfill the null geodesic condition.*/mydbl condition_factor(mydbl r, double a){  return (mydbl)(1.0/a)*sqrtl((1.0+2.0*potential(r)/(C*C))/(1.0 - 2.0*potential(r)/(C*C)));}