Marsyas  0.6.0-alpha
/usr/src/RPM/BUILD/marsyas-0.6.0/src/marsyas/basis.cpp
Go to the documentation of this file.
00001 /**********************************************************
00002 *  Reduced and modified version of basis.c by J-P Moreau  *
00003 * ------------------------------------------------------- *
00004 * Reference for original basis.c:                         *
00005 *                                                         *
00006 * "Numerical Algorithms with C, by Gisela Engeln-Muellges *
00007 *  and Frank Uhlig, Springer-Verlag, 1996".               *
00008 **********************************************************/
00009 #include <marsyas/basis.h>
00010 
00011 REAL sqr(REAL x)                    /* square a floating point number */
00012 /***********************************************************************
00013 * Compute the square of a floating point number.                       *
00014 *                                                                      *
00015 * global names used:                                                   *
00016 * ==================                                                   *
00017 * REAL                                                                 *
00018 ***********************************************************************/
00019 {
00020   return x * x;
00021 }
00022 
00023 REAL norm_max      /* Find the maximum norm of a REAL vector .........*/
00024 (
00025   REAL vektor[],               /* vector .................*/
00026   int  n                       /* length of vector .......*/
00027 )                             /* Maximum norm ...........*/
00028 
00029 /***********************************************************************
00030 * Return the maximum norm of a [0..n-1] vector  v.                     *
00031 *                                                                      *
00032 * global names used:                                                   *
00033 * ==================                                                   *
00034 * REAL, FABS, ZERO                                                     *
00035 ***********************************************************************/
00036 {
00037   REAL norm,                                             /* local max */
00038        betrag;                            /* magnitude of a component */
00039 
00040   for (n--, norm = ZERO; n >= 0; n--, vektor++)
00041     if ((betrag = FABS(*vektor)) > norm)
00042       norm = betrag;
00043 
00044   return norm;
00045 }
00046 
00047 void copy_vector        /* copy a REAL vector ........................*/
00048 (
00049   REAL ziel[],            /* copied vector ............*/
00050   REAL quelle[],          /* original vector ..........*/
00051   int  n                  /* length of vector .........*/
00052 )
00053 /***********************************************************************
00054 * copy the n elements of the vector quelle into the vector ziel.       *
00055 *                                                                      *
00056 * global names used:                                                   *
00057 * ==================                                                   *
00058 * REAL                                                                 *
00059 ***********************************************************************/
00060 {
00061   for (n--; n >= 0; n--)
00062     *ziel++ = *quelle++;
00063 }
00064 
00065 /***********************************************************************
00066 * read a vector x of size n from input file fp (index begins at zero)  *
00067 *                                                                      *
00068 * global names used:                                                   *
00069 * ==================                                                   *
00070 * REAL                                                                 *
00071 ***********************************************************************/
00072 int ReadVec (FILE *fp, int n, REAL x[])
00073 {
00074   int i;
00075   REAL tmp;
00076 
00077   for (i = 0; i < n; ++i)
00078   {
00079     if (fscanf (fp,FORMAT_IN, &tmp) <= 0) return (-1);
00080     x[i] = (REAL) tmp;
00081   }
00082 
00083   return (0);
00084 }
00085 
00086 /***********************************************************************
00087 * read a vector x of size n from input file fp (index begins at one)   *
00088 *                                                                      *
00089 * global names used:                                                   *
00090 * ==================                                                   *
00091 * REAL                                                                 *
00092 ***********************************************************************/
00093 int ReadVec1 (FILE *fp, int n, REAL x[])
00094 {
00095   int i;
00096   REAL tmp;
00097 
00098   for (i = 1; i < n+1; ++i)
00099   {
00100     if (fscanf (fp,FORMAT_IN, &tmp) <= 0) return (-1);
00101     x[i] = (REAL) tmp;
00102   }
00103 
00104   return (0);
00105 }
00106 
00107 void SetVec (int n, REAL x[], REAL val)
00108 {
00109   int i;
00110 
00111   for (i = 0; i < n; ++i)
00112     x[i] = val;
00113 }
00114 
00115 int WriteVec (FILE *fp, int n, REAL x[])
00116 /*====================================================================*
00117  *                                                                    *
00118  *  Put out vector of length n to output file. Index begins at zero.  *
00119  *                                                                    *
00120  *====================================================================*
00121  *                                                                    *
00122  *  Input parameters:                                                 *
00123  *  ================                                                  *
00124  *                                                                    *
00125  *      n        int n;                                               *
00126  *               lenght of vector                                     *
00127  *      x        REAL x[];                                            *
00128  *               vector                                               *
00129  *                                                                    *
00130  *   Return value :                                                   *
00131  *   =============                                                    *
00132  *      =  0     All ok                                               *
00133  *      = -1     Error putting out to stdout                          *
00134  *                                                                    *
00135  *====================================================================*/
00136 {
00137   int i;
00138 
00139   for (i = 0; i < n; ++i)
00140     if (fprintf (fp,FORMAT_126LF, x[i]) <= 0) return (-1);
00141   if (fprintf (fp,"\n") <= 0) return (-1);
00142 
00143   return 0;
00144 }
00145 
00146 int WriteVec1 (FILE *fp, int n, REAL x[])
00147 /*====================================================================*
00148  *                                                                    *
00149  *  Put out vector of length n to output file. Index begins at one.   *
00150  *                                                                    *
00151  *====================================================================*
00152  *                                                                    *
00153  *  Input parameters:                                                 *
00154  *  ================                                                  *
00155  *                                                                    *
00156  *      n        int n;                                               *
00157  *               lenght of vector                                     *
00158  *      x        REAL x[];                                            *
00159  *               vector                                               *
00160  *                                                                    *
00161  *   Return value :                                                   *
00162  *   =============                                                    *
00163  *      =  0     All ok                                               *
00164  *      = -1     Error putting out to stdout                          *
00165  *                                                                    *
00166  *====================================================================*/
00167 {
00168   int i;
00169 
00170   for (i = 1; i < n+1; ++i)
00171     if (fprintf (fp,FORMAT_126LF, x[i]) <= 0) return (-1);
00172   if (fprintf (fp,"\n") <= 0) return (-1);
00173 
00174   return 0;
00175 }
00176 
00177 REAL pi() {
00178   return 4.0*atan(1.0);
00179 }
00180 
00181 void CopyMat (int m, int n, REAL * source[], REAL * dest[])
00182 /*====================================================================*
00183  *                                                                    *
00184  *  Copy the m x n matrix source to the  m x n matrix dest.           *
00185  *                                                                    *
00186  *====================================================================*
00187  *                                                                    *
00188  *  Input parameters:                                                 *
00189  *  ================                                                  *
00190  *      m        int m; ( m > 0 )                                     *
00191  *               numnber of rows of matrix                            *
00192  *      n        int n; ( n > 0 )                                     *
00193  *               number of columns of matrix                          *
00194  *      source   REAL * source[];                                     *
00195  *               matrix                                               *
00196  *      dest     REAL * dest[];                                       *
00197  *               matrix to be copied to                               *
00198  *                                                                    *
00199  *   Output parameter:                                                *
00200  *   ================                                                 *
00201  *      dest     same as above                                        *
00202  *                                                                    *
00203  *   ATTENTION : WE do not allocate storage for dest here.            *
00204  *                                                                    *
00205  *====================================================================*/
00206 {
00207   int i, j;
00208 
00209   for (i = 0; i < m; ++i)
00210     for (j = 0; j < n; j++)
00211       dest[i][j] = source[i][j];
00212 }
00213 
00214 REAL maxroot(void)    /* Root of the largest representable number ....*/
00215 /***********************************************************************
00216 * Compute the square root of the largest machine number 2 ^ (MAX_EXP/2)*
00217 * if not already done                                                  *
00218 *                                                                      *
00219 * global names used:                                                   *
00220 * ==================                                                   *
00221 * REAL, boolean, FALSE, TRUE, SQRT, POSMAX                             *
00222 ***********************************************************************/
00223 {
00224   static REAL       save_maxroot;
00225   static boolean    schon_berechnet = FALSE;
00226   REAL              faktor;
00227   unsigned long int n;
00228 
00229   if (! schon_berechnet)
00230   {
00231     save_maxroot = ONE;
00232     faktor       = TWO;
00233     for (n = MAX_EXP / 2; n > 1; n /= 2, faktor *= faktor)
00234       if (n % 2 != 0)
00235         save_maxroot *= faktor;
00236     save_maxroot    *= faktor;
00237     schon_berechnet  = TRUE;
00238   }
00239 
00240   return save_maxroot;
00241 }
00242 
00243 
00244 void quadsolv           /* Complex quadratic equation ................*/
00245 (
00246   REAL    ar,        /* second degree coefficient .......*/
00247   REAL    ai,
00248   REAL    br,        /* linear coefficient ..............*/
00249   REAL    bi,
00250   REAL    cr,        /* polynomial constant .............*/
00251   REAL    ci,
00252   REAL *  tr,        /* solution ........................*/
00253   REAL *  ti
00254 )
00255 /*====================================================================*
00256  *                                                                    *
00257  *  Compute the least magnitude solution of the quadratic equation    *
00258  *  a * t**2 + b * t + c = 0. Here a, b, c and t are complex.         *
00259  *                                         2                          *
00260  *  Formeula used: t = 2c / (-b +/- sqrt (b  - 4ac)).                 *
00261  *  This formula is valid for a=0 .                                   *
00262  *                                                                    *
00263  *====================================================================*
00264  *                                                                    *
00265  *  Input parameters:                                                 *
00266  *  ================                                                  *
00267  *      ar, ai   coefficient of t**2             REAL   ar, ai;       *
00268  *      br, bi   coefficient of t                REAL   br, bi;       *
00269  *      cr, ci   constant term                   REAL   cr, ci;       *
00270  *                                                                    *
00271  *  Output parameter:                                                 *
00272  *  ================                                                  *
00273  *      tr, ti   complex solution of minimal magnitude                *
00274  *                                               REAL   *tr, *ti;     *
00275  *                                                                    *
00276  *  Macro used :     SQRT                                             *
00277  *  ============                                                      *
00278  *                                                                    *
00279  *====================================================================*/
00280 {
00281   REAL pr, pi, qr, qi, h;
00282 
00283   pr = br * br - bi * bi;
00284   pi = TWO * br * bi;                       /*  p = b * b             */
00285 
00286   qr = ar * cr - ai * ci;
00287   qi = ar * ci + ai * cr;                   /*  q = a * c             */
00288 
00289   pr = pr - (REAL)4.0 * qr;
00290   pi = pi - (REAL)4.0 * qi;                 /* p = b * b - 4 * a * c  */
00291 
00292   h  = SQRT (pr * pr + pi * pi);            /* q = sqrt (p)           */
00293 
00294   qr = h + pr;
00295   if (qr > ZERO)
00296     qr = SQRT (qr * HALF);
00297   else
00298     qr = ZERO;
00299 
00300   qi = h - pr;
00301   if (qi > ZERO)
00302     qi = SQRT (qi * HALF);
00303   else
00304     qi = ZERO;
00305 
00306   if (pi < ZERO) qi = -qi;
00307 
00308   h = qr * br + qi * bi;     /* p = -b +/- q, choose sign for large  */
00309   /* magnitude  p                         */
00310   if (h > ZERO)
00311   {
00312     qr = -qr;
00313     qi = -qi;
00314   }
00315 
00316   pr = qr - br;
00317   pi = qi - bi;
00318   h = pr * pr + pi * pi;                      /* t = (2 * c) / p      */
00319 
00320   if (h == ZERO)
00321   {
00322     *tr = ZERO;
00323     *ti = ZERO;
00324   }
00325   else
00326   {
00327     *tr = TWO * (cr * pr + ci * pi) / h;
00328     *ti = TWO * (ci * pr - cr * pi) / h;
00329   }
00330 } //quadsolv
00331 
00332 
00333 REAL comabs             /* Complex absolute value ....................*/
00334 (
00335   REAL  ar,          /* Real part .......................*/
00336   REAL  ai           /* Imaginary part ..................*/
00337 )
00338 /*====================================================================*
00339  *                                                                    *
00340  *  Complex absolute value of   a                                     *
00341  *                                                                    *
00342  *====================================================================*
00343  *                                                                    *
00344  *   Input parameters:                                                *
00345  *   ================                                                 *
00346  *      ar,ai    REAL   ar, ai;                                       *
00347  *               Real, imaginary parts of  a                          *
00348  *                                                                    *
00349  *   Return value :                                                   *
00350  *   =============                                                    *
00351  *      Absolute value of a                                           *
00352  *                                                                    *
00353  *   Macros used :    SQRT, ABS, SWAP                                 *
00354  *   =============                                                    *
00355  *                                                                    *
00356  *====================================================================*/
00357 {
00358   if (ar == ZERO && ai == ZERO) return (ZERO);
00359 
00360   ar = ABS (ar);
00361   ai = ABS (ai);
00362 
00363   if (ai > ar)                                  /* Switch  ai and ar */
00364     SWAP (REAL, ai, ar)
00365 
00366     return ((ai == ZERO) ? (ar) : (ar * SQRT (ONE + ai / ar * ai / ar)));
00367 }
00368 
00369 
00370 int comdiv              /* Complex division ..........................*/
00371 (
00372   REAL   ar,            /* Real part of numerator ..........*/
00373   REAL   ai,            /* Imaginary part of numerator .....*/
00374   REAL   br,            /* Real part of denominator ........*/
00375   REAL   bi,            /* Imaginary part of denominator ...*/
00376   REAL * cr,            /* Real part of quotient ...........*/
00377   REAL * ci             /* Imaginary part of quotient ......*/
00378 )
00379 /*====================================================================*
00380  *                                                                    *
00381  *  Complex division  c = a / b                                       *
00382  *                                                                    *
00383  *====================================================================*
00384  *                                                                    *
00385  *   Input parameters:                                                *
00386  *   ================                                                 *
00387  *      ar,ai    REAL   ar, ai;                                       *
00388  *               Real, imaginary parts of numerator                   *
00389  *      br,bi    REAL   br, bi;                                       *
00390  *               Real, imaginary parts of denominator                 *
00391  *                                                                    *
00392  *   Output parameters:                                               *
00393  *   ==================                                               *
00394  *      cr,ci    REAL   *cr, *ci;                                     *
00395  *               Real , imaginary parts of the quotient               *
00396  *                                                                    *
00397  *   Return value :                                                   *
00398  *   =============                                                    *
00399  *      = 0      ok                                                   *
00400  *      = 1      division by 0                                        *
00401  *                                                                    *
00402  *   Macro used :     ABS                                             *
00403  *   ============                                                     *
00404  *                                                                    *
00405  *====================================================================*/
00406 {
00407   REAL tmp;
00408 
00409   if (br == ZERO && bi == ZERO) return (1);
00410 
00411   if (ABS (br) > ABS (bi))
00412   {
00413     tmp  = bi / br;
00414     br   = tmp * bi + br;
00415     *cr  = (ar + tmp * ai) / br;
00416     *ci  = (ai - tmp * ar) / br;
00417   }
00418   else
00419   {
00420     tmp  = br / bi;
00421     bi   = tmp * br + bi;
00422     *cr  = (tmp * ar + ai) / bi;
00423     *ci  = (tmp * ai - ar) / bi;
00424   }
00425 
00426   return (0);
00427 }
00428 
00429 int ReadMat (FILE *fp, int m, int n, REAL * a[])
00430 /*====================================================================*
00431  *                                                                    *
00432  *  Read an m x n matrix from input file. Index starts at zero.       *
00433  *                                                                    *
00434  *====================================================================*
00435  *                                                                    *
00436  *  Input parameters :                                                *
00437  *  ==================                                                *
00438  *      m        int m; ( m > 0 )                                     *
00439  *               number of rows of matrix                             *
00440  *      n        int n; ( n > 0 )                                     *
00441  *               column number of  matrix                             *
00442  *      a        REAL * a[];                                          *
00443  *               matrix                                               *
00444  *                                                                    *
00445  *   Output parameter:                                                *
00446  *   ================                                                 *
00447  *      a        matrix                                               *
00448  *                                                                    *
00449  *   ATTENTION : WE do not allocate storage for a here.               *
00450  *                                                                    *
00451  *====================================================================*/
00452 {
00453   int i, j;
00454   double x;
00455 
00456   for (i = 0; i < m; ++i)
00457     for (j = 0; j < n; j++)
00458     {
00459       if (fscanf (fp,FORMAT_IN, &x) <= 0) return (-1);
00460       a[i][j] = (REAL) x;
00461     }
00462 
00463   return (0);
00464 }
00465 
00466 //same as ReadMat, but beginning at index 1
00467 int ReadMat1 (FILE *fp, int m, int n, REAL * a[])
00468 {
00469   int i, j;
00470   double x;
00471 
00472   for (i = 1; i < m+1; ++i)
00473     for (j = 1; j < n+1; j++)
00474     {
00475       if (fscanf (fp,FORMAT_IN, &x) <= 0) return (-1);
00476       a[i][j] = (REAL) x;
00477     }
00478 
00479   return (0);
00480 }
00481 
00482 int WriteMat (FILE *fp, int m, int n, REAL * a[])
00483 /*====================================================================*
00484  *                                                                    *
00485  *  Put out an m x n matrix in output file. Index starts at zero.     *
00486  *                                                                    *
00487  *====================================================================*
00488  *                                                                    *
00489  *  Input parameters:                                                 *
00490  *  ================                                                  *
00491  *      m        int m; ( m > 0 )                                     *
00492  *               row number of matrix                                 *
00493  *      n        int n; ( n > 0 )                                     *
00494  *               column number of matrix                              *
00495  *      a        REAL * a[];                                          *
00496  *               matrix                                               *
00497  *                                                                    *
00498  *   Return value :                                                   *
00499  *   =============                                                    *
00500  *      =  0      all put out                                         *
00501  *      = -1      Error writing onto stdout                           *
00502  *                                                                    *
00503  *====================================================================*/
00504 {
00505   int i, j;
00506 
00507   if (fprintf (fp,"\n") <= 0) return (-1);
00508 
00509   for (i = 0; i < m; ++i)
00510   {
00511     for (j = 0; j < n; j++)
00512       if (fprintf (fp,FORMAT_126LF, a[i][j]) <= 0) return (-1);
00513 
00514     if (fprintf (fp,"\n") <= 0) return (-1);
00515   }
00516   if (fprintf (fp,"\n") <= 0) return (-1);
00517 
00518   return (0);
00519 }
00520 
00521 //same as WriteMat, but beginning at index 1
00522 int WriteMat1 (FILE *fp, int m, int n, REAL * a[])
00523 {
00524   int i, j;
00525 
00526   if (fprintf (fp,"\n") <= 0) return (-1);
00527 
00528   for (i = 1; i < m+1; ++i)
00529   {
00530     for (j = 1; j < n+1; j++)
00531       if (fprintf (fp,FORMAT_126LF, a[i][j]) <= 0) return (-1);
00532 
00533     if (fprintf (fp,"\n") <= 0) return (-1);
00534   }
00535   if (fprintf (fp,"\n") <= 0) return (-1);
00536 
00537   return (0);
00538 }
00539 
00540 
00541 void SetMat (int m, int n, REAL * a[], REAL val)
00542 /*====================================================================*
00543  *                                                                    *
00544  *  Initialize an m x n matrix with a constant value val .            *
00545  *                                                                    *
00546  *====================================================================*
00547  *                                                                    *
00548  *  Input parameters:                                                 *
00549  *  ================                                                  *
00550  *      m        int m; ( m > 0 )                                     *
00551  *               row number of matrix                                 *
00552  *      n        int n; ( n > 0 )                                     *
00553  *               column number of matrix                              *
00554  *      a        REAL * a[];                                          *
00555  *               matrix                                               *
00556  *      val      constant value                                       *
00557  *                                                                    *
00558  *   Output parameter:                                                *
00559  *   ================                                                 *
00560  *      a        matrix with constant value val in every position     *
00561  *                                                                    *
00562  *====================================================================*/
00563 {
00564   int i, j;
00565 
00566   for (i = 0; i < m; ++i)
00567     for (j = 0; j < n; j++)
00568       a[i][j] = val;
00569 }
00570 
00571 static char Separator[] =
00572   "--------------------------------------------------------------------";
00573 
00574 int WriteHead (FILE *fp, char * string)
00575 /*====================================================================*
00576  *                                                                    *
00577  *  Put out header with text in string in file fp.                    *
00578  *                                                                    *
00579  *====================================================================*
00580  *                                                                    *
00581  *  Input parameters:                                                 *
00582  *  ================                                                  *
00583  *      string   char *string;                                        *
00584  *               text of headertext (last byte is a 0)                *
00585  *                                                                    *
00586  *   Return value :                                                   *
00587  *   =============                                                    *
00588  *      =  0      All ok                                              *
00589  *      = -1      Error writing onto stdout                           *
00590  *      = -2      Invalid text for header                             *
00591  *                                                                    *
00592  *====================================================================*/
00593 {
00594   if (string == NULL) return (-2);
00595 
00596   if (fprintf (fp,"\n%s\n%s\n%s\n\n", Separator, string, Separator) <= 0)
00597     return (-1);
00598 
00599   return 0;
00600 }
00601 
00602 int WriteHead1( char * string)
00603 /*====================================================================*
00604  *                                                                    *
00605  *  Put out header with text in string in stdout.                     *
00606  *                                                                    *
00607  *====================================================================*
00608  *                                                                    *
00609  *  Input parameters:                                                 *
00610  *  ================                                                  *
00611  *      string   char *string;                                        *
00612  *               text of headertext (last byte is a 0)                *
00613  *                                                                    *
00614  *   Return value :                                                   *
00615  *   =============                                                    *
00616  *      =  0      All ok                                              *
00617  *      = -1      Error writing onto stdout                           *
00618  *      = -2      Invalid text for header                             *
00619  *                                                                    *
00620  *====================================================================*/
00621 {
00622   if (string == NULL) return (-2);
00623 
00624   if (printf ("\n%s\n%s\n%s\n\n", Separator, string, Separator) <= 0)
00625     return (-1);
00626 
00627   return 0;
00628 }
00629 
00630 int WriteEnd (FILE *fp)
00631 /*====================================================================*
00632  *                                                                    *
00633  *  Put out end of writing onto file fp.                              *
00634  *                                                                    *
00635  *====================================================================*
00636  *                                                                    *
00637  *   Return value :                                                   *
00638  *   =============                                                    *
00639  *      =  0      All ok                                              *
00640  *      = -1      error writing onto stdout                           *
00641  *                                                                    *
00642  *====================================================================*/
00643 {
00644   if (fprintf (fp,"\n%s\n\n", Separator) <= 0) return (-1);
00645   return 0;
00646 }
00647 
00648 int WriteEnd1(void)
00649 /*====================================================================*
00650  *                                                                    *
00651  *  Put out end of writing onto  stdout.                              *
00652  *                                                                    *
00653  *====================================================================*
00654  *                                                                    *
00655  *   Return value :                                                   *
00656  *   =============                                                    *
00657  *      =  0      All ok                                              *
00658  *      = -1      error writing onto stdout                           *
00659  *                                                                    *
00660  *====================================================================*/
00661 {
00662   if (printf ("\n%s\n\n", Separator) <= 0) return (-1);
00663   return 0;
00664 }
00665 
00666 void LogError (char * string, int rc, char * file, int line)
00667 /*====================================================================*
00668  *                                                                    *
00669  *  Put out error message onto  stdout.                               *
00670  *                                                                    *
00671  *====================================================================*
00672  *                                                                    *
00673  *  Input parameters:                                                 *
00674  *  ================                                                  *
00675  *      string   char *string;                                        *
00676  *               text of error massage (final byte is 0)              *
00677  *      rc       int rc;                                              *
00678  *               error code                                           *
00679  *      file     char *file;                                          *
00680  *               name of C file in which error was encountered        *
00681  *      line     int line;                                            *
00682  *               line number of C file with error                     *
00683  *                                                                    *
00684  *====================================================================*/
00685 {
00686   if (string == NULL)
00687   {
00688     printf ("Unknown ERROR in file %s at line %d\n", file, line);
00689     return;
00690   }
00691 
00692   if (rc == 0)
00693     printf ("ERROR: %s, File %s, Line %d\n", string, file, line);
00694   else
00695     printf ("ERROR: %s, rc = %d, File %s, Line %d\n",
00696             string, rc, file, line);
00697   return;
00698 }
00699 
00700 
00701 REAL epsroot(void)  /* Compute square root of the machine constant ...*/
00702 /***********************************************************************
00703 * Compute square root of the machine constant, if not already done.    *
00704 *                                                                      *
00705 * global names used:                                                   *
00706 * ==================                                                   *
00707 * REAL, boolean, FALSE, TRUE, SQRT, MACH_EPS                           *
00708 ***********************************************************************/
00709 
00710 {
00711   static REAL    save_mach_eps_root;
00712   static boolean schon_berechnet     = FALSE;
00713 
00714   if (! schon_berechnet)
00715     schon_berechnet    = TRUE,
00716     save_mach_eps_root = SQRT(MACH_EPS);
00717 
00718   return save_mach_eps_root;
00719 }
00720 
00721 
00722 REAL epsquad(void)      /* Find the machine constant squared .........*/
00723 /***********************************************************************
00724 * Compute the square of the machine constant, if not already done.     *
00725 *                                                                      *
00726 * global names used:                                                   *
00727 * ==================                                                   *
00728 * REAL, boolean, FALSE, TRUE, MACH_EPS                                 *
00729 ***********************************************************************/
00730 
00731 {
00732   static REAL    save_mach_eps_quad;
00733   static boolean schon_berechnet     = FALSE;
00734 
00735   if (! schon_berechnet)
00736     schon_berechnet    = TRUE,
00737     save_mach_eps_quad = MACH_EPS * MACH_EPS;
00738 
00739   return save_mach_eps_quad;
00740 }
00741 
00742 // End of file basis_r.cpp