clog.c

gridgen是一款强大的网格生成程序

 * Also tested by csin(casin(z)) = z.
 */

double complex
casin (z)
     double complex z;
{
  double complex w;
  static double complex ca, ct, zz, z2;
  double x, y;

  x = creal (z);
  y = cimag (z);

  if (y == 0.0)
    {
      if (fabs(x) > 1.0)
	{
	  w = PIO2 + 0.0 * I;
	  mtherr ("casin", DOMAIN);
	}
      else
	{
	  w = asin (x) + 0.0 * I;
	}
      return (w);
    }

/* Power series expansion */
/*
b = cabs(z);
if( b < 0.125 )
{
z2.r = (x - y) * (x + y);
z2.i = 2.0 * x * y;

cn = 1.0;
n = 1.0;
ca.r = x;
ca.i = y;
sum.r = x;
sum.i = y;
do
	{
	ct.r = z2.r * ca.r  -  z2.i * ca.i;
	ct.i = z2.r * ca.i  +  z2.i * ca.r;
	ca.r = ct.r;
	ca.i = ct.i;

	cn *= n;
	n += 1.0;
	cn /= n;
	n += 1.0;
	b = cn/n;

	ct.r *= b;
	ct.i *= b;
	sum.r += ct.r;
	sum.i += ct.i;
	b = fabs(ct.r) + fabs(ct.i);
	}
while( b > MACHEP );
w->r = sum.r;
w->i = sum.i;
return;
}
*/


  ca = x + y * I;
  ct = ca * I;
	/* sqrt( 1 - z*z) */
  /* cmul( &ca, &ca, &zz ) */
  /*x * x  -  y * y */
  zz = (x - y) * (x + y) + (2.0 * x * y) * I;

  zz = 1.0 - creal(zz) - cimag(zz) * I;
  z2 = csqrt (zz);

  zz = ct + z2;
  zz = clog (zz);
  /* multiply by 1/i = -i */
  w = zz * (-1.0 * I);
  return (w);
}
/*							cacos()
 *
 *	Complex circular arc cosine
 *
 *
 *
 * SYNOPSIS:
 *
 * double complex cacos();
 * double complex z, w;
 *
 * w = cacos (z);
 *
 *
 *
 * DESCRIPTION:
 *
 *
 * w = arccos z  =  PI/2 - arcsin z.
 *
 *
 *
 *
 * ACCURACY:
 *
 *                      Relative error:
 * arithmetic   domain     # trials      peak         rms
 *    DEC       -10,+10      5200      1.6e-15      2.8e-16
 *    IEEE      -10,+10     30000      1.8e-14      2.2e-15
 */

double complex
cacos (z)
     double complex z;
{
  double complex w;

  w = casin (z);
  w = (PIO2  -  creal (w)) - cimag (w) * I;
  return (w);
}
/*							catan()
 *
 *	Complex circular arc tangent
 *
 *
 *
 * SYNOPSIS:
 *
 * double complex catan();
 * double complex z, w;
 *
 * w = catan (z);
 *
 *
 *
 * DESCRIPTION:
 *
 * If
 *     z = x + iy,
 *
 * then
 *          1       (    2x     )
 * Re w  =  - arctan(-----------)  +  k PI
 *          2       (     2    2)
 *                  (1 - x  - y )
 *
 *               ( 2         2)
 *          1    (x  +  (y+1) )
 * Im w  =  - log(------------)
 *          4    ( 2         2)
 *               (x  +  (y-1) )
 *
 * Where k is an arbitrary integer.
 *
 * catan(z) = -i catanh(iz).
 *
 * ACCURACY:
 *
 *                      Relative error:
 * arithmetic   domain     # trials      peak         rms
 *    DEC       -10,+10      5900       1.3e-16     7.8e-18
 *    IEEE      -10,+10     30000       2.3e-15     8.5e-17
 * The check catan( ctan(z) )  =  z, with |x| and |y| < PI/2,
 * had peak relative error 1.5e-16, rms relative error
 * 2.9e-17.  See also clog().
 */

double complex
catan (z)
     double complex z;
{
  double complex w;
  double a, t, x, x2, y;

  x = creal (z);
  y = cimag (z);

  if ((x == 0.0) && (y > 1.0))
    goto ovrf;

  x2 = x * x;
  a = 1.0 - x2 - (y * y);
  if (a == 0.0)
    goto ovrf;

  t = 0.5 * atan2 (2.0 * x, a);
  w = redupi (t);

  t = y - 1.0;
  a = x2 + (t * t);
  if (a == 0.0)
    goto ovrf;

  t = y + 1.0;
  a = (x2 + (t * t))/a;
  w = w + (0.25 * log (a)) * I;
  return (w);

ovrf:
  mtherr ("catan", OVERFLOW);
  w = MAXNUM + MAXNUM * I;
  return (w);
}


/*							csinh
 *
 *	Complex hyperbolic sine
 *
 *
 *
 * SYNOPSIS:
 *
 * double complex csinh();
 * double complex z, w;
 *
 * w = csinh (z);
 *
 * DESCRIPTION:
 *
 * csinh z = (cexp(z) - cexp(-z))/2
 *         = sinh x * cos y  +  i cosh x * sin y .
 *
 * ACCURACY:
 *
 *                      Relative error:
 * arithmetic   domain     # trials      peak         rms
 *    IEEE      -10,+10     30000       3.1e-16     8.2e-17
 *
 */

double complex
csinh (z)
     double complex z;
{
  double complex w;
  double x, y;

  x = creal(z);
  y = cimag(z);
  w = sinh (x) * cos (y)  +  (cosh (x) * sin (y)) * I;
  return (w);
}


/*							casinh
 *
 *	Complex inverse hyperbolic sine
 *
 *
 *
 * SYNOPSIS:
 *
 * double complex casinh();
 * double complex z, w;
 *
 * w = casinh (z);
 *
 *
 *
 * DESCRIPTION:
 *
 * casinh z = -i casin iz .
 *
 * ACCURACY:
 *
 *                      Relative error:
 * arithmetic   domain     # trials      peak         rms
 *    IEEE      -10,+10     30000       1.8e-14     2.6e-15
 *
 */

double complex
casinh (z)
     double complex z;
{
  double complex w;

  w = -1.0 * I * casin (z * I);
  return (w);
}


/*							ccosh
 *
 *	Complex hyperbolic cosine
 *
 *
 *
 * SYNOPSIS:
 *
 * double complex ccosh();
 * double complex z, w;
 *
 * w = ccosh (z);
 *
 *
 *
 * DESCRIPTION:
 *
 * ccosh(z) = cosh x  cos y + i sinh x sin y .
 *
 * ACCURACY:
 *
 *                      Relative error:
 * arithmetic   domain     # trials      peak         rms
 *    IEEE      -10,+10     30000       2.9e-16     8.1e-17
 *
 */

double complex
ccosh (z)
     double complex z;
{
  double complex w;
  double x, y;

  x = creal(z);
  y = cimag(z);
  w = cosh (x) * cos (y)  +  (sinh (x) * sin (y)) * I;
  return (w);
}



/*							cacosh
 *
 *	Complex inverse hyperbolic cosine
 *
 *
 *
 * SYNOPSIS:
 *
 * double complex cacosh();
 * double complex z, w;
 *
 * w = cacosh (z);
 *
 *
 *
 * DESCRIPTION:
 *
 * acosh z = i acos z .
 *
 * ACCURACY:
 *
 *                      Relative error:
 * arithmetic   domain     # trials      peak         rms
 *    IEEE      -10,+10     30000       1.6e-14     2.1e-15
 *
 */

double complex
cacosh (z)
     double complex z;
{
  double complex w;

  w = I * cacos (z);
  return (w);
}


/*							ctanh
 *
 *	Complex hyperbolic tangent
 *
 *
 *
 * SYNOPSIS:
 *
 * double complex ctanh();
 * double complex z, w;
 *
 * w = ctanh (z);
 *
 *
 *
 * DESCRIPTION:
 *
 * tanh z = (sinh 2x  +  i sin 2y) / (cosh 2x + cos 2y) .
 *
 * ACCURACY:
 *
 *                      Relative error:
 * arithmetic   domain     # trials      peak         rms
 *    IEEE      -10,+10     30000       1.7e-14     2.4e-16
 *
 */

double complex
ctanh (z)
     double complex z;
{
  double complex w;
  double x, y, d;

  x = creal(z);
  y = cimag(z);
  d = cosh (2.0 * x) + cos (2.0 * y);
  w = sinh (2.0 * x) / d  +  (sin (2.0 * y) / d) * I;
  return (w);
}


/*							catanh
 *
 *	Complex inverse hyperbolic tangent
 *
 *
 *
 * SYNOPSIS:
 *
 * double complex catanh();
 * double complex z, w;
 *
 * w = catanh (z);
 *
 *
 *
 * DESCRIPTION:
 *
 * Inverse tanh, equal to  -i catan (iz);
 *
 * ACCURACY:
 *
 *                      Relative error:
 * arithmetic   domain     # trials      peak         rms
 *    IEEE      -10,+10     30000       2.3e-16     6.2e-17
 *
 */

double complex
catanh (z)
     double complex z;
{
  double complex w;

  w = -1.0 * I * catan (z * I);
  return (w);
}

/*							cpow
 *
 *	Complex power function
 *
 *
 *
 * SYNOPSIS:
 *
 * double complex cpow();
 * double complex a, z, w;
 *
 * w = cpow (a, z);
 *
 *
 *
 * DESCRIPTION:
 *
 * Raises complex A to the complex Zth power.
 * Definition is per AMS55 # 4.2.8,
 * analytically equivalent to cpow(a,z) = cexp(z clog(a)).
 *
 * ACCURACY:
 *
 *                      Relative error:
 * arithmetic   domain     # trials      peak         rms
 *    IEEE      -10,+10     30000       9.4e-15     1.5e-15
 *
 */

double complex
cpow (a, z)
     double complex a, z;
{
  double complex w;
  double x, y, r, theta, absa, arga;

  x = creal (z);
  y = cimag (z);
  absa = cabs (a);
  if (absa == 0.0)
    {
      return (0.0 + 0.0 * I);
    }
  arga = carg (a);
  r = pow (absa, x);
  theta = x * arga;
  if (y != 0.0)
    {
      r = r * exp (-y * arga);
      theta = theta + y * log (absa);
    }
  w = r * cos (theta) + (r * sin (theta)) * I;
  return (w);
}