e.cpp

一个非常有用的开源代码

}           

  
void quadsolv           /* Complex quadratic equation ................*/
             (
               REAL    ar,        /* second degree coefficient .......*/
               REAL    ai,
               REAL    br,        /* linear coefficient ..............*/
               REAL    bi,
               REAL    cr,        /* polynomial constant .............*/
               REAL    ci,
               REAL *  tr,        /* solution ........................*/
               REAL *  ti
             )
/*====================================================================*
 *                                                                    *
 *  Compute the least magnitude solution of the quadratic equation    *
 *  a * t**2 + b * t + c = 0. Here a, b, c and t are complex.         *
 *                                         2                          *
 *  Formeula used: t = 2c / (-b +/- sqrt (b  - 4ac)).                 *
 *  This formula is valid for a=0 .                                   *
 *                                                                    *
 *====================================================================*
 *                                                                    *
 *  Input parameters:                                                 *
 *  ================                                                  *
 *      ar, ai   coefficient of t**2             REAL   ar, ai;       *
 *      br, bi   coefficient of t                REAL   br, bi;       *
 *      cr, ci   constant term                   REAL   cr, ci;       *
 *                                                                    *
 *  Output parameter:                                                 *
 *  ================                                                  *
 *      tr, ti   complex solution of minimal magnitude                *
 *                                               REAL   *tr, *ti;     *
 *                                                                    *
 *  Macro used :     SQRT                                             *
 *  ============                                                      *
 *                                                                    *
 *====================================================================*/
{
  REAL pr, pi, qr, qi, h;

  pr = br * br - bi * bi;
  pi = TWO * br * bi;                       /*  p = b * b             */

  qr = ar * cr - ai * ci;
  qi = ar * ci + ai * cr;                   /*  q = a * c             */

  pr = pr - (REAL)4.0 * qr;
  pi = pi - (REAL)4.0 * qi;                 /* p = b * b - 4 * a * c  */

  h  = SQRT (pr * pr + pi * pi);            /* q = sqrt (p)           */

  qr = h + pr;
  if (qr > ZERO)
    qr = SQRT (qr * HALF);
  else
    qr = ZERO;

  qi = h - pr;
  if (qi > ZERO)
    qi = SQRT (qi * HALF);
  else
    qi = ZERO;

  if (pi < ZERO) qi = -qi;

  h = qr * br + qi * bi;     /* p = -b +/- q, choose sign for large  */
                             /* magnitude  p                         */
  if (h > ZERO)
  {
    qr = -qr;
    qi = -qi;
  }

  pr = qr - br;
  pi = qi - bi;
  h = pr * pr + pi * pi;                      /* t = (2 * c) / p      */

  if (h == ZERO)
  {
    *tr = ZERO;
    *ti = ZERO;
  }
  else
  {
    *tr = TWO * (cr * pr + ci * pi) / h;
    *ti = TWO * (ci * pr - cr * pi) / h;
  }
} //quadsolv 
   
   
REAL comabs             /* Complex absolute value ....................*/
              (
               REAL  ar,          /* Real part .......................*/
               REAL  ai           /* Imaginary part ..................*/
              )
/*====================================================================*
 *                                                                    *
 *  Complex absolute value of   a                                     *
 *                                                                    *
 *====================================================================*
 *                                                                    *
 *   Input parameters:                                                *
 *   ================                                                 *
 *      ar,ai    REAL   ar, ai;                                       *
 *               Real, imaginary parts of  a                          *
 *                                                                    *
 *   Return value :                                                   *
 *   =============                                                    *
 *      Absolute value of a                                           *
 *                                                                    *
 *   Macros used :    SQRT, ABS, SWAP                                 *
 *   =============                                                    *
 *                                                                    *
 *====================================================================*/
{
  if (ar == ZERO && ai == ZERO) return (ZERO);

  ar = ABS (ar);
  ai = ABS (ai);

  if (ai > ar)                                  /* Switch  ai and ar */
    SWAP (REAL, ai, ar)

  return ((ai == ZERO) ? (ar) : (ar * SQRT (ONE + ai / ar * ai / ar)));
}   
 
 
int comdiv              /* Complex division ..........................*/
           (
            REAL   ar,            /* Real part of numerator ..........*/
            REAL   ai,            /* Imaginary part of numerator .....*/
            REAL   br,            /* Real part of denominator ........*/
            REAL   bi,            /* Imaginary part of denominator ...*/
            REAL * cr,            /* Real part of quotient ...........*/
            REAL * ci             /* Imaginary part of quotient ......*/
           )
/*====================================================================*
 *                                                                    *
 *  Complex division  c = a / b                                       *
 *                                                                    *
 *====================================================================*
 *                                                                    *
 *   Input parameters:                                                *
 *   ================                                                 *
 *      ar,ai    REAL   ar, ai;                                       *
 *               Real, imaginary parts of numerator                   *
 *      br,bi    REAL   br, bi;                                       *
 *               Real, imaginary parts of denominator                 *
 *                                                                    *
 *   Output parameters:                                               *
 *   ==================                                               *
 *      cr,ci    REAL   *cr, *ci;                                     *
 *               Real , imaginary parts of the quotient               *
 *                                                                    *
 *   Return value :                                                   *
 *   =============                                                    *
 *      = 0      ok                                                   *
 *      = 1      division by 0                                        *
 *                                                                    *
 *   Macro used :     ABS                                             *
 *   ============                                                     *
 *                                                                    *
 *====================================================================*/
{
  REAL tmp;

  if (br == ZERO && bi == ZERO) return (1);

  if (ABS (br) > ABS (bi))
  {
    tmp  = bi / br;
    br   = tmp * bi + br;
    *cr  = (ar + tmp * ai) / br;
    *ci  = (ai - tmp * ar) / br;
  }
  else
  {
    tmp  = br / bi;
    bi   = tmp * br + bi;
    *cr  = (tmp * ar + ai) / bi;
    *ci  = (tmp * ai - ar) / bi;
 }

 return (0);
}       
   
int ReadMat (FILE *fp, int m, int n, REAL * a[])
/*====================================================================*
 *                                                                    *
 *  Read an m x n matrix from input file. Index starts at zero.       *
 *                                                                    *
 *====================================================================*
 *                                                                    *
 *  Input parameters :                                                *
 *  ==================                                                *
 *      m        int m; ( m > 0 )                                     *
 *               number of rows of matrix                             *
 *      n        int n; ( n > 0 )                                     *
 *               column number of  matrix                             *
 *      a        REAL * a[];                                          *
 *               matrix                                               *
 *                                                                    *
 *   Output parameter:                                                *
 *   ================                                                 *
 *      a        matrix                                               *
 *                                                                    *
 *   ATTENTION : WE do not allocate storage for a here.               *
 *                                                                    *
 *====================================================================*/
{
  int i, j;
  double x;

  for (i = 0; i < m; i++)
    for (j = 0; j < n; j++)
    {
      if (fscanf (fp,FORMAT_IN, &x) <= 0) return (-1);
      a[i][j] = (REAL) x;
    }

  return (0);
}

//same as ReadMat, but beginning at index 1
int ReadMat1 (FILE *fp, int m, int n, REAL * a[])
{
  int i, j;
  double x;

  for (i = 1; i < m+1; i++)
    for (j = 1; j < n+1; j++)
    {
      if (fscanf (fp,FORMAT_IN, &x) <= 0) return (-1);
      a[i][j] = (REAL) x;
    }

  return (0);
}

int WriteMat (FILE *fp, int m, int n, REAL * a[])
/*====================================================================*
 *                                                                    *
 *  Put out an m x n matrix in output file. Index starts at zero.     *
 *                                                                    *
 *====================================================================*
 *                                                                    *
 *  Input parameters:                                                 *
 *  ================                                                  *
 *      m        int m; ( m > 0 )                                     *
 *               row number of matrix                                 *
 *      n        int n; ( n > 0 )                                     *
 *               column number of matrix                              *
 *      a        REAL * a[];                                          *
 *               matrix                                               *
 *                                                                    *
 *   Return value :                                                   *
 *   =============                                                    *
 *      =  0      all put out                                         *
 *      = -1      Error writing onto stdout                           *
 *                                                                    *
 *====================================================================*/
{
  int i, j;

  if (fprintf (fp,"\n") <= 0) return (-1);

  for (i = 0; i < m; i++)
  {
    for (j = 0; j < n; j++)
      if (fprintf (fp,FORMAT_126LF, a[i][j]) <= 0) return (-1);

    if (fprintf (fp,"\n") <= 0) return (-1);
  }
  if (fprintf (fp,"\n") <= 0) return (-1);

  return (0);
}   

//same as WriteMat, but beginning at index 1
int WriteMat1 (FILE *fp, int m, int n, REAL * a[])
{
  int i, j;

  if (fprintf (fp,"\n") <= 0) return (-1);

  for (i = 1; i < m+1; i++)
  {
    for (j = 1; j < n+1; j++)
      if (fprintf (fp,FORMAT_126LF, a[i][j]) <= 0) return (-1);

    if (fprintf (fp,"\n") <= 0) return (-1);
  }
  if (fprintf (fp,"\n") <= 0) return (-1);

  return (0);
}