1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
|
package
global(polyval,polyadd,roots)
{
/*
A polynomial p(x) = p[1]*x^n + p[2]*x^(n-1) + ... + p[n-1]*x + p[n]
is represented by vector #(p[1],...,p[n]). That is, the first element
of the vector is the higher order coefficient.
*/
function y = polyval(p,x)
/* polyval(p,x) evaluates polynomial p at argument x.
How polynonials are represented, see roots. */
{
n = length(p);
y = p[1];
for (i=2; i<=n; i++)
y = y*x + p[i];
};
function y = polyadd(p,q)
/* polyadd(p,q) forms the sum of two polynomials (which are represented
in vector form).
How polynonials are represented, see roots. */
{
Lp = length(p);
Lq = length(q);
L = max(Lp,Lq);
y = izeros(L) + 0*p[1] + 0*q[1];
y[L-Lp+1:L] = p;
y[L-Lq+1:L] = y[L-Lq+1:L] + q;
};
function y = roots(p)
/* roots(p) returns the vector of zeros of polynomial p.
It uses the companion matrix method.
A polynomial p(x) = p[1]*x^n + p[2]*x^(n-1) + ... + p[n]
is represented by vector #(p[1],...,p[n]). That is, the first
element of the vector is the higher order coefficient
See also: fsolve*/
{
n = length(p);
n1 = n - 1;
p1 = p[1];
A = zeros(n1,n1) + 0*p1;
if (n > 2) {
A[n:n:n1^2] = 1;
A[1,:] = -p[2:n]/p1;
} else if (n == 2) {
A[1,1] = -p[2]/p1;
} else {
y = #();
return;
};
y = eig(A);
};
}
|
