bairstow method code

#include<iostream>
#include<cstdlib>
#include<math.h>
using namespace std;
int main()
{
    int n,i,N1,e;
          float p,q,r,pr,qr,ps,qs,D,jam1,jam2,s,k;
    cin>>N1>>n>>p>>q>>e;
    float a[n],b[n-1],dpb[n-1],dqb[n-1];
    for(i=0;i<=n;i++)
        {
            cin>>a[i];
    }
    while(N1--){
            b[n-2]=a[n];
            b[n-3]=a[n-1]-p*b[n-2];
            for(i=2;i<=n-2;i++)
            {
                b[n-(i+2)]=a[n-i]-p*b[n-(i+1)]-q*b[n-i];
            }

            r=a[1]-p*b[0]-q*b[1];
            s=a[0]-q*b[0];
            k=fabs(r);
          if(fabs(r)<=e && fabs(s)<=e){goto jam;}
            dpb[n-2]=0;
            dpb[n-3]=-b[n-2];
             for(i=2;i<=n-2;i++)
            {
                dpb[n-(i+2)]=-b[n-i-1]-p*dpb[n-i-1]-q*dpb[n-i];
            }
            dqb[n-2]=0;
            dqb[n-3]=0;
            dqb[n-4]=-b[n-2];
             for(i=3;i<=n-2;i++)
            {
                dqb[n-(i+2)]=-p*dqb[n-i-1]-q*dpb[n-i]-b[n-i];
            }
            pr=-b[0]-p*dpb[0]-q*dpb[1];
            qr=-b[1]-p*dqb[0]-q*dqb[1];
            ps=-q*dpb[0];
            qs=-b[0]-q*dqb[0];
         D=pr*qs-(ps*qr);
            jam1=(-r*qs+s*qr)/D;
            jam2=(-s*pr+r*ps)/D;
            p=p+jam1;
            q=q+jam2;
}
    jam:
        cout<<r<<"\t"<<s<<endl<<p<<"\t"<<q;

}