Loading...
Searching...
No Matches
qreval.c
Go to the documentation of this file.
/* qreval.c CCMATH mathematics library source code.
*
* Copyright (C) 2000 Daniel A. Atkinson All rights reserved.
* This code may be redistributed under the terms of the GNU library
* public license (LGPL). ( See the lgpl.license file for details.)
* ------------------------------------------------------------------------
*/
#include "ccmath.h"
int qreval(double *ev, double *dp, int n)
{
double cc, sc = 0.0, d, x, y, h, tzr = 1.e-15;
int j, k, m, mqr = 8 * n;
for (j = 0, m = n - 1;; ++j) {
while (1) {
if (m < 1)
return 0;
k = m - 1;
if (fabs(dp[k]) <= fabs(ev[m]) * tzr)
--m;
else {
x = (ev[k] - ev[m]) / 2.;
h = sqrt(x * x + dp[k] * dp[k]);
if (m > 1 && fabs(dp[m - 2]) > fabs(ev[k]) * tzr)
break;
x += ev[m];
ev[m--] = x - h;
ev[m--] = x + h;
}
}
if (j > mqr)
return -1;
if (x > 0.)
d = ev[m] + x - h;
else
d = ev[m] + x + h;
cc = 1.;
y = 0.;
ev[0] -= d;
for (k = 0; k < m; ++k) {
x = ev[k] * cc - y;
y = dp[k] * cc;
h = sqrt(x * x + dp[k] * dp[k]);
if (k > 0)
dp[k - 1] = sc * h;
ev[k] = cc * h;
cc = x / h;
sc = dp[k] / h;
ev[k + 1] -= d;
y *= sc;
ev[k] = cc * (ev[k] + y) + ev[k + 1] * sc * sc + d;
}
ev[k] = ev[k] * cc - y;
dp[k - 1] = ev[k] * sc;
ev[k] = ev[k] * cc + d;
}
return 0;
}
ccmath.h
qreval
int qreval(double *ev, double *dp, int n)
Definition qreval.c:9
x
#define x