Initially you said you wanted to solve
(3*x)*R**2 - (2*x)*R + y =0, which has polynomial coefficients {3*x, -2*x, y}
But later you use
1*R**2 + 1*R + 2*e1 = 0
Which is it?
In any case, to get the roots of a low-degree polynomial, use the POLYROOT function in SAS IML. Create a SAS data set that has the polynomial coefficients ordered from highest degree to lowest (the constant term). Read the coefficients into IML and loop over the coefficients.
A polynomial of degree d always has d complex roots. If you are interested only in the real roots, then find the roots for which the imaginary part is 0 (or very tiny).
The following code assumes that you want to solve the equation in your original post. It then creates a SAS data set that has ONLY the polynomial coefficients. The IML program extracts the REAL roots.
data have;
input y x;
datalines;
-1 2
-2 3
-5 6
-8 9
;
/* root of
(3*x)*R**2 - 2*x*R + y = 0
Let a = 3*x, b= -2*x, and c=y and use the quadratic equation
*/
data Coef;
keep c2 c1 c0;
set have;
c2 = 3*x;
c1 = -2*x;
c0 = y;
run;
proc iml;
delta = 1E-14; /* small number to use as a custoff value */
use Coef;
read all var _num_ into coeff; /* coefficients of R^2, R, constant */
close;
numPolys = nrow(coeff);
deg = ncol(coeff)-1;
roots = j(numPolys, deg, .); /* poly of degree d has d complex roots */
do i = 1 to numPolys;
r = polyroot( coeff[i,] ); /* returns all complex roots */
print i r;
/* store only the real roots for which r[,2]=0 */
realIdx = loc( abs(r[,2]) < delta );
numRealRts = ncol(realIdx);
if numRealRts > 0 then
roots[i, 1:numRealRts] = rowvec( r[realIdx, 1] ); /* get the real parts */
end;
PolyNumber = T(1:numPolys);
RootNumber = 1:deg;
print PolyNumber roots[c=RootNumber];
You can use the same IML program to find roots of higher-degree polynomials. For example, if you want to find cubic roots, use the following DATA step to create the coefficients and use the same IML code:
data Coef;
input c3 c2 c1 c0;
datalines;
1 -1 2 1
1 -1 -4 4
1 1 -3 -3
;
... View more
