Archivos metodistas

Codigo en C++ para el método de Lin Bairstow

#include <stdlib.h>

#include <stdio.h>

#include <math.h>

#define MAX_TERMS 21

int precision_error_flag;

double a[MAX_TERMS],b[MAX_TERMS],c[MAX_TERMS],d[MAX_TERMS];

void find_poly_roots(int n)

{

                               double r,s,dn,dr,ds,drn,dsn,eps;

                int i,iter;

                r = s = 0;

                dr = 1.0;

                ds = 0;

                eps = 1e-14;

                iter = 1;

                while ((fabs(dr)+fabs(ds)) > eps) {

                               if ((iter % 200) == 0) {

                                               r=(double)rand()/16000.;

                               }

                               if ((iter % 500) == 0) {

                                               eps*=10.0;

                                               precision_error_flag=1;

                                               printf("Loss of precision\n");

                               }

                               b[1] = a[1] - r;

                               c[1] = b[1] - r;

                               for (i=2;i<=n;i++){

                                               b[i] = a[i] - r * b[i-1] - s * b[i-2];

                                               c[i] = b[i] - r * c[i-1] - s * c[i-2];

                               }

                               dn=c[n-1] * c[n-3] - c[n-2] * c[n-2];

                               drn=b[n] * c[n-3] - b[n-1] * c[n-2];

                               dsn=b[n-1] * c[n-1] - b[n] * c[n-2];

                               if (fabs(dn) < 1e-16) {

                                               dn = 1;

                                               drn = 1;

                                               dsn = 1;

                               }

                               dr = drn / dn;

                               ds = dsn / dn;

                               r += dr;

                               s += ds;

                               iter++;

                }

                for (i=0;i<n-1;i++)

                               a[i] = b[i];

                a[n] = s;

                a[n-1] = r;

}

void main()

{

                int i,n,order=0;

                double tmp;

                while ((order < 2) || (order > MAX_TERMS-1)) {

                               printf("Polynomial order (2-20): ");

                               scanf("%d",&order);

                }

                printf("Enter coefficients, high order to low order.\n");

                for (i=0;i<=order;i++){

                               printf("a(%d)=",i);

                               scanf("%lf",&a[i]);

                               if (i==0) {

                                               tmp=a[i];

                                               a[i]=1.0;

                               }

                               else a[i]/=tmp;

                               d[i]=a[i];

                }

                b[0]=c[0]=1.0;

                n=order;

                precision_error_flag=0;

                while (n > 2) {

                               find_poly_roots(n);

                                                               n-=2;

                }

                printf("The quadratic factors are:\n");

               

                for (i=order;i>=2;i-=2) /* print quadratics */

                                printf("t^2  %+.15lg t  %+.15lg\n",a[i-1],a[i]);

                if ((n % 2) == 1)

                               printf("The linear term is: \nt  %+.15lg\n",a[1]);

}