extern int csolvpoly(int cynpoly,unsigned int * * *cyapoly,unsigned int *cy);  //error
extern int cyfindmax();
extern int cdivideoutroot();
extern int cymeshpatch();
extern int cyrealroot();
extern int cysaveistep();
extern int cyrecoveristep();
extern int cytestroot();
extern int cywhateverf();
extern int cyiterate();
extern int cyfunc();
extern int cyro();
extern int cyiterq();
extern int cyoriglq();
extern int cpolydeclarations();   
extern int cyspecialcases();
extern int cmultiplyout();
extern int cmultpolys(int cyn1in,unsigned int * * *cypoly1in,unsigned int cyn2in,unsigned int * * *cypoly2in); //error
extern int crotatepoly(int cynpoly,unsigned int * * *cyapoly,unsigned int *cy);  //cy=angle
extern int crotateroots(int cynpoly, unsigned int * * *croot,unsigned int *cy);


int cyfindflag;
int cykk; 
int cyj9,cynnn;
int cypp,cyk8; 
int cyi7,cyj7,cyk7,cym7;
int cyisavezit;
int cyj2,cyn1in,cyn2in; 

int cyniter;
int cynten1=100;
int cynitermax;
int cyn2,cyi,cyj,cyk,cynpasses,cynpolymax; 
int cynisnow,cynpassq,cymessage2,cymessage;
int cyimax,cynpassq2,cymessagesave5;
int cyerrorflag;
int cyponk;
int cynpoly;

realscalar cyxoldmax_ram;unsigned int *cyxoldmax=cyxoldmax_ram.s;   //ok realscalar
realscalar cyxmax_ram;unsigned int *cyxmax=cyxmax_ram.s;
realscalar cyxleft_ram;unsigned int *cyxleft=cyxleft_ram.s;
realscalar cyfmax_ram;unsigned int *cyfmax=cyfmax_ram.s;          //ok
realscalar cyx_ram;unsigned int *cyx=cyx_ram.s;
realscalar cyw_ram;unsigned int *cyw=cyw_ram.s;
realscalar cy_ram;unsigned int *cy=cy_ram.s;
realscalar cywlower_ram;unsigned int *cywlower=cywlower_ram.s;     //ok
realscalar cydel_ram;unsigned int *cydel=cydel_ram.s;
realscalar cyf_ram;unsigned int *cyf=cyf_ram.s;                    //ok
realscalar cyfprev_ram;unsigned int *cyfprev=cyfprev_ram.s;         //ok
realscalar cyfsave5_ram;unsigned int *cyfsave5=cyfsave5_ram.s;
realscalar cyfprevsave5_ram;unsigned int *cyfprevsave5=cyfprevsave5_ram.s;  //ok
realscalar cyxprev_ram;unsigned int *cyxprev=cyxprev_ram.s;       //ok
realscalar cyx2save_ram;unsigned int *cyx2save=cyx2save_ram.s;    //ok
realscalar cy1D16_ram;unsigned int *cy1D16=cy1D16_ram.s;
realscalar cy1Dm5_ram;unsigned int *cy1Dm5=cy1Dm5_ram.s;
realscalar cytemp_ram;unsigned int *cytemp=cytemp_ram.s;
realscalar cytemp1_ram;unsigned int *cytemp1=cytemp1_ram.s;
realscalar cyxsave_ram;unsigned int *cyxsave=cyxsave_ram.s;
realscalar cyfsave_ram;unsigned int *cyfsave=cyfsave_ram.s;   //21 realscalar

unsigned int * *cyq;        //6 real vectors specially declared so as not
unsigned int * *cyr;        //to need zmaxsizepad in cpoly.h
unsigned int * *cyqpr;
unsigned int * *cyrpr;
unsigned int * *cyxrot;      
unsigned int * *cyyrot;
int cyispecialcases=cyspecialcases();
int cyspecialcases()
  {cyxrot = new unsigned int * [cmaxsizepad];
   cyyrot = new unsigned int * [cmaxsizepad];
   cyq = new unsigned int * [cmaxsizepad];                 
   cyr = new unsigned int * [cmaxsizepad];
   cyrpr = new unsigned int * [cmaxsizepad];
   cyqpr = new unsigned int * [cmaxsizepad];
   for (cyi=0;cyi<cmaxsizepad;cyi++)  
     {cyxrot[cyi]=new unsigned int [zprecpad];
      cyyrot[cyi]=new unsigned int [zprecpad];
      cyq[cyi]=new unsigned int [zprecpad]; 
      cyr[cyi]=new unsigned int [zprecpad];
      cyrpr[cyi]=new unsigned int [zprecpad]; 
      cyqpr[cyi]=new unsigned int [zprecpad];
     }
   return 0;
  }

complexscalar cylam_ram;unsigned int * *cylam=cylam_ram.cs;
complexscalar cylamroot_ram;unsigned int * *cylamroot=cylamroot_ram.cs;
complexscalar cyuf_ram;unsigned int * *cyuf=cyuf_ram.cs;
complexscalar crot_ram;unsigned int * *crot=crot_ram.cs;      //4 complexscalar

complexvector cyuq_ram;unsigned int * * *cyuq=cyuq_ram.cv;
complexvector cya_ram;unsigned int * * *cya=cya_ram.cv;
complexvector cyb_ram;unsigned int * * *cyb=cyb_ram.cv;
complexvector cypoly1in_ram;unsigned int * * *cypoly1in=cypoly1in_ram.cv;
complexvector cypoly2in_ram;unsigned int * * *cypoly2in=cypoly2in_ram.cv;
complexvector cyc_ram;unsigned int * * *cyc=cyc_ram.cv;
complexvector cyapoly_ram;unsigned int * * *cyapoly=cyapoly_ram.cv;
complexvector croot_ram;unsigned int * * *croot=croot_ram.cv;       //8 complexvector
complexvector ccpolyout_ram;unsigned int * * *ccpolyout=ccpolyout_ram.cv;

unsigned int *cymess;  //special cmaxsizepad
int ccpolydeclarations=cpolydeclarations();
int cpolydeclarations()   
  {cymess = new unsigned int [cmaxsizepad];
   return 0;
  }

int csolvpoly(int cynpoly,unsigned int * * *cyapoly,unsigned int *cy)
    {
   int zmultiplicity=1;
   cynpasses=0;
   cyn2=cynpoly;
   cconinttomp(1,0);
   cstore(cylamroot);           //cylamroot starts at 1+i0
   for (cyi=0;cyi<cyn2+1;cyi++)        
     {crecall(cyapoly[cyi]);          
      cstore(cya[cyi]);}
   cynisnow=3;

//rotate the poly:

//zconinttomp(2);     //set cy=0.2 radians
//zdivint(z,10);
//zstore(cy);

crotatepoly(cyn2,cya,cy);  //cy=angle,output is cyapoly

//end of rotate the poly

cywayback:

   for (cyi=0;cyi<cyn2;cyi++)
     {cdiv(cya[cyi],cya[cyn2]);   //dividing out coeff of x^n
      cstore(cya[cyi]);}
   cconinttomp(1,0);
   cstore(cya[cyn2]);
   cstore(cylam);        //for use when no scaling is done as below
                         //(the code between /* and */ omitted
/*
   cln(cya[0]);          //log of last coeff
   cdivint(cz,cyn2);     //divided by order
   cexp(cz);             //exp of that cylam=nth root of last coeff=cyn
   cstore(cylam);
   cmult(cylamroot,cylam);
   cstore(cylamroot);
   for (cyi=1;cyi<cyn2+1;cyi++)   //this will make first coeff and
     {crectopol(cylam);           //last coeff 1.
      cxtonthpol(cz,cyi);    
      cpoltorec(cz);
      cmult(cya[cyi],cz);
      cdiv(cz,cya[0]);       //divides out last root
      cstore(cya[cyi]);
     }
   cconinttomp(1,0);
   cstore(cya[0]);
*/

   cyi=cyn2/2;
   cyi=2*cyi;
   if (cyi!=cyn2)
     {for (cyi=0;cyi<cyn2+1;cyi++)
       {cconinttomp(0,1);
        cmult(cz,cya[cyi]);
        cstore(cya[cyi]);}
     }
   for (cyi=0;cyi<cyn2+1;cyi++)
     {zzero(cyq[cyi]);
      zzero(cyr[cyi]);
      zzero(cyqpr[cyi]);
      zzero(cyrpr[cyi]);        
     }
   for (cyi=0;cyi<prec+1;cyi++)
     {cyxoldmax[cyi]=0;}
   cyxoldmax[prec+1]=2;
   cyxoldmax[prec]=1;              
   cyxoldmax[0]=0x40000000+9999;
   cynpassq=0;
   cyerrorflag=0;
   cynisnow=3;                           
   cyerrorflag=cyerrorflag+1;
   cyfindmax();                    
   zsub(cyxmax,cyxoldmax);
   if (z[prec+1]==2 && z[prec]>0)
     {zrecall(cyxoldmax);
      zstore(cyxmax);
      goto cybhyu;
     }
   zrecall(cyxmax);
   zstore(cyxoldmax);
cybhyu:
    if (cynpassq==0)
      {zrecall(cyxmax);
       zchs(z);
       zstore(cyxleft);
       cynpassq=1;
      }
    cymessage2=3;
    cyrealroot();
    if (cymessage==1)
      {zrecall(cyxmax);
       zchs(z);
       zmult(z,cyxmax);
       zstore(cywlower);
       cymeshpatch();         
       cerr<<'\n';
       if (cymessage==9)      //at this point could take sqrt of poly.
         {cerr<<"cannot find all roots"<<'\n';
          cerr<<"number found "<<cynpasses<<'\n';
          cerr<<"continuing search for multiple roots"<<'\n';
          cyi=cyn2/2;
          cyi=2*cyi;
          if (cyi!=cyn2)
             {zmessage="not an even number of multiple roots";
              zbreak();}

//          for (cyi=0;cyi<cyn2+1;cyi++) cprnt(cya[cyi]);
//          zmessage="remaining series--sqrt will be taken";
//          zbreak();

          cconinttomp(1,0);
          cdivint(cz,2);
          cstore(csp);
          cstonth(cyn2,cya,csp);
          cyn2=cyn2/2;

          for (cyi=0;cyi<cyn2+1;cyi++)
             {crecall(csout[cyi]);
              cstore(cya[cyi]);
             }

//          for (cyi=0;cyi<cyn2+1;cyi++) cprnt(cya[cyi]);
//          zmessage="sqrt of remaining series";
//          zbreak();
          
//          cmultpolys(cyn2,cya,cyn2,cya);  //square for check

//          for (cyi=0;cyi<2*cyn2+1;cyi++) cprnt(cyb[cyi]);
//          zmessage="sqrt squared";         //checked perfectly
//          zbreak();

//now: what goes here?  Need to go to cywayback with a message that
//                      a new regime has been established.
//                      Introduce a multiplicity variable

          zmultiplicity=2*zmultiplicity;
          cerr<<"going to cywayback after doubling zmultiplicity"<<'\n';
          goto cywayback; }
//          goto cyahead12;}

       cerr<<"going to cywayback"<<'\n';
       goto cywayback;
      }                      
    if (cymessage==3)           //in case a root is found
      {
//       cynpasses=cynpasses+1;
//       cymess[cynpasses]=3;
       zrecall(cyx);
       zstore(crot[0]);
       zrecall(cyw);
       zstore(crot[1]);
       cmult(crot,cylamroot);
//       cstore(croot[cynpasses]);
       for (cyi=1;cyi<zmultiplicity+1;cyi++)
         {cynpasses=cynpasses+1;
          cymess[cynpasses]=3;
          cstore(croot[cynpasses]);}
       cerr<<"r";
       cdivideoutroot();      //divide out root only once
       if (cyn2==0)
         {cerr<<'\n';
          goto cyahead12;}
       zrecall(cyx);
       zstore(cyxleft);
       goto cywayback;
      }
cyahead2:
      zmessage="error in main before yahead12";
      zbreak();

cyahead12:
     crotateroots(cynpoly,croot,cy);    //can't be cyn2
  return 0;
  }

int cymeshpatch()
  {cynten1=2*cynten1;
   if (cynten1>=10000)
     {cymessage=9;}
   return 0;
  }

int cyfindmax()                  
  {zconinttomp(1);
   zdivint(z,10);                         
   zstore(cydel);
   cyfindflag=0;
   zconinttomp(2);
   zsqrt(z);
   zsqrt(z);
   zsqrt(z);
   zstore(cytemp1);
cyback:
   if (cyfindflag==0)
     {zconinttomp(1);
      zstore(cyx);
     }
   if (cyfindflag==1)
     {zdiv(cyx,cytemp1);
      zstore(cyx);
     }
   if (cyfindflag==2)
     {zmult(cyx,cytemp1);
      zstore(cyx);
     }
   cabs(cya[cyn2]);
   zstore(cyf);
   for (cyi=cyn2-1;cyi>=0;cyi--)
     {zmult(cyf,cyx);
      zstore(cyf);
      cabs(cya[cyi]);
      zsub(cyf,z);
      zstore(cyf);
     }
   if (cyf[prec+1]==0)
     {if (cyfindflag==0)
        {cyfindflag=2;             
         goto cyback;
        }
      if (cyfindflag==1)           
        {zmultint(cyx,2);              
         zstore(cyx);
         goto cyfinishit;}
      if (cyfindflag==2)          
        {goto cyback;}
      }
   if (cyf[prec]==0)
      {goto cyfinishit;}
   if (cyf[prec+1]==2 && cyf[prec]>0)
     {if (cyfindflag==2)           
        {goto cyfinishit;}
      if (cyfindflag==0)
        {cyfindflag=1;             
         goto cyback;
        }
      if (cyfindflag==1)          
        {goto cyback;}            
     }
cyfinishit: zrecall(cyx);
   zstore(cyxmax);
   return 0;
  }

int cdivideoutroot()
  {for (cyj=0;cyj<cyn2+1;cyj++)
     {czero(cyc[cyj]);}
   crecall(cya[cyn2]);
   cstore(cyc[cyn2-1]);
   for (cykk=cyn2-2;cykk>=0;cykk--)
     {cmult(crot,cyc[cykk+1]);
      cadd(cz,cya[cykk+1]);
      cstore(cyc[cykk]);
     }
   for (cykk=0;cykk<cyn2;cykk++)
     {crecall(cyc[cykk]);
      cstore(cya[cykk]);
     }
   cyn2=cyn2-1;
   return 0;
  }

int cyrealroot()
  {int cyhelp;
   cyimax=cynten1;
   zconinttomp(cyimax);       
   zdiv(cyxmax,z);
   zmultint(z,100);
   zdivint(z,99);
   zstore(cydel);
   cynpassq2=0;
   for (cyi=-cyimax;cyi<cyimax+1;cyi++)
     {cyhelp=abs(cyi);
      zmultint(cydel,cyhelp);
      if (cyi<0) zchs(z);
      zstore(cyx);
      zconinttomp(1);
      zdivint(z,10);
      zadd(z,cyx);
      zsub(z,cyxleft);
      if (z[prec+1]==0)
        {goto cyni5;}
      cywhateverf();
      if (cynpassq2==0)
        {cynpassq2=1;
         goto cyni5;
        }
      if (cyi==-cyimax)
        {goto cyni5;}
      if (cyf[prec]==0)                                  
        {cysaveistep();
         cyponk=3;
         goto cyahead500;          
        }                                                          
      zmult(cyfprev,cyf);                                            
      if (z[prec+1]==0)
        {cysaveistep();
         cyiterate();
         cyponk=4;
         goto cyahead500;            
        }
      goto cyni5;
      zmessage="it got below goto cyni5 somehow";
      zbreak();
cyahead500: zchs(z);                     
      cytestroot();                        
      if (cymessage==5)
           {cymessage=3;
            goto cyahead200;}
      cyrecoveristep();
cyni5: zrecall(cyf);
      zstore(cyfprev);
      zrecall(cyx);
      zstore(cyxprev);
     }
   cymessage=1;
cyahead200:
   return 0;
  }

int cywhateverf()
  {if (cymessage2==3)
     {cyro();} 
   return 0;
  }

int cysaveistep()
  {zrecall(cyf);
   zstore(cyfsave5);
   zrecall(cyfprev);
   zstore(cyfprevsave5);
   cymessagesave5=cymessage;
   return 0;
  }

int cyiterate() 
  {cynitermax=24*prec;    
   cyniter=0;
cyback1:
   cyniter=cyniter+1;
   zrecall(cyx);
   zstore(cyxsave);
   zrecall(cyf);
   zstore(cyfsave);
   zadd(cyx,cyxprev);
   zdivint(z,2);
   zstore(cyx); 
   cywhateverf();               
   if (cyf[prec]==0)
     {goto cyxit;}
   zmult(cyfprev,cyf);
   if (z[prec+1]==0)
     {if (cyniter>cynitermax)
       {goto cyxit;}
      goto cyback1;
     }
   if (cyniter>cynitermax)
     {zrecall(cyxprev);
      zstore(cyx);
      goto cyxit;
     }
   zrecall(cyf);
   zstore(cyfprev);
   zrecall(cyx);
   zstore(cyxprev);
   zrecall(cyfsave);
   zstore(cyf);
   zrecall(cyxsave);
   zstore(cyx);
   goto cyback1;
   zmessage="error in yiterate";
   zbreak();
cyxit: return 0;
  }

int cytestroot()
   {if (cymessage2==3)
     {
      if (cyq[1][prec]!=0)
        {zdiv(cyq[0],cyq[1]);
         zchs(z);
         zstore(cyw);
         cymessage=5;
         goto cyahead21;}
     }
    cymessage=6;
cyahead21:
   return 0;
  }

int cyrecoveristep() 
  {zrecall(cyfsave5);
   zstore(cyf);
   zrecall(cyfprevsave5);
   zstore(cyfprev);
   cymessage=cymessagesave5;
   return 0;
  }

int cyfunc()
   {crecall(cya[cyn2]);
    cstore(cyuf);
    for (cyj2=cyn2-1;cyj2>=0;cyj2--)
      {cmultcbyr(cyuf,cyx);
       cadd(cz,cya[cyj2]);
       cstore(cyuf);
      }
    return 0;
   }

int cyro()
   {cyisavezit=cyi; 
    cyoriglq();
    cyiterq();
    zrecall(cyr[0]);
    zstore(cyf);
    cyi=cyisavezit;
    return 0;
   }
      
int cyoriglq()         
   {cyfunc();  
    crecall(cyuf);
    cstore(cyuq[cyn2]);               
    zrecall(cyuq[cyn2][0]);
    zstore(cyq[cyn2]);
    zrecall(cyuq[cyn2][1]);
    zstore(cyr[cyn2]);
    for (cyi7=0;cyi7<cyn2+1;cyi7++)  
      {crecall(cya[cyi7]);
       cstore(cyb[cyi7]);
      }
    for (cyk7=1;cyk7<cyn2+1;cyk7++)    
      {for (cyj7=1;cyj7<cyn2-cyk7+2;cyj7++)    
          {cmultint(cyb[cyj7],cyj7);               
           cdivint(cz,cyk7);
           cstore(cyb[cyj7-1]);
           cconinttomp(0,1);
           cmult(cz,cyb[cyj7-1]);
           cstore(cyb[cyj7-1]);
          }                  
       czero(cyuq[cyn2-cyk7]);
       for (cyi7=0;cyi7<cyn2-cyk7+1;cyi7++)     
          {zxtonth(cyx,cyi7);                 
           zstore(crfac);
           cmultcbyr(cyb[cyi7],crfac);                
           cadd(cz,cyuq[cyn2-cyk7]);
           cstore(cyuq[cyn2-cyk7]);         
          }
       zrecall(cyuq[cyn2-cyk7][0]); 
       zstore(cyq[cyn2-cyk7]);
       zrecall(cyuq[cyn2-cyk7][1]);
       zstore(cyr[cyn2-cyk7]);
      }                                        
    for (cykk=0;cykk<cyn2+1;cykk++)       
      {zrecall(cyq[cyn2-cykk]);               
       zstore(cyqpr[cykk]);                 
      }
    for (cykk=0;cykk<cyn2;cykk++)
      {zrecall(cyr[cyn2-cykk]);
       zstore(cyrpr[cykk]);
      }
    for (cykk=0;cykk<cyn2+1;cykk++)           
      {zrecall(cyqpr[cykk]);
       zstore(cyq[cykk]);
       zrecall(cyrpr[cykk]);
       zstore(cyr[cykk]);
      } 
    return 0;
   }

int cyiterq() 
   {for (cypp=cyn2;cypp>=2;cypp--)        
      {for (cyk8=0;cyk8<cypp;cyk8++)
        {if (cyk8==0)
           {zmult(cyq[cyk8],cyr[cypp-1]);
            zstore(cyqpr[cyk8]);
            goto cyapril;
           }
         zmult(cyq[cyk8],cyr[cypp-1]);
         zstore(cytemp);
         zmult(cyq[cypp],cyr[cyk8-1]);
         zchs(z);
         zadd(z,cytemp);
         zstore(cyqpr[cyk8]);
cyapril:  zrecall(cyr[cyk8]);
         zstore(cyrpr[cyk8]);
        }                          
       for (cyk8=0;cyk8<cypp;cyk8++)
        {zrecall(cyrpr[cyk8]);  
         zstore(cyq[cyk8]);
        }
       for (cyk8=0;cyk8<cypp-1;cyk8++)
        {zmult(cyqpr[cypp-1],cyrpr[cyk8]); 
         zstore(cytemp);
         zmult(cyqpr[cyk8],cyrpr[cypp-1]);
         zsub(z,cytemp);
         zstore(cyr[cyk8]);
        }
      }
    return 0;
   }

int cmultpolys(int cyn1in,unsigned int * * *cypoly1in,unsigned int cyn2in,unsigned int * * *cypoly2in)
   {int yj9,yj8;
   for (yj9=0;yj9<cyn1in+cyn2in+1;yj9++)
     {czero(cyb[yj9]);}
   for (yj8=0;yj8<cyn1in+cyn2in+1;yj8++) 
     {for (yj9=0;yj9<cyn1in+1;yj9++)
        {if (yj8-yj9>cyn2in) goto cybkyl; 
         if (yj8-yj9<0) goto cybkyl;
         cmult(cypoly1in[yj9],cypoly2in[yj8-yj9]);
         cadd(cz,cyb[yj8]);
         cstore(cyb[yj8]);
cybkyl: continue;                        //output is cyb
        }
     }
   return 0;
  }


int cmultiplyout()      
  {for (cyj9=0;cyj9<cynpasses+1;cyj9++)  
     {czero(cyb[cyj9]);}
   cynnn=0;   
   cconinttomp(1,0);
   cstore(cyb[0]);
   for (cyi=1;cyi<cynpasses+1;cyi++)  
     {if (cymess[cyi]==3)
        {cynnn=cynnn+1;
         cconinttomp(1,0);
         cstore(cyb[cynnn]);
         for (cyj=cynnn-1;cyj>=1;cyj--)
           {cmult(croot[cyi],cyb[cyj]);
            cchs(cz);                                            
            cadd(cz,cyb[cyj-1]);
            cstore(cyb[cyj]);
           }
         cmult(croot[cyi],cyb[0]);
         cchs(cz);
         cstore(cyb[0]);
        }
     }
cyfred: zchs(z);         
   return 0;
  }                     //why not unsigned int

int crotatepoly(int cynpoly,unsigned int * * *cyapoly,unsigned int *cy)  //cy=angle
  {cexpitheta(cy);
   cstore(cyr);
   cconinttomp(1,0);
   cstore(cyrpr);
   for (cykk=0;cykk<cynpoly+1;cykk++)
     {cmult(cyapoly[cykk],cyrpr);
      cstore(cya[cykk]);       //store rotated poly in cya
      cmult(cyrpr,cyr);
      cstore(cyrpr);
     }
   return 0;
  }

int crotateroots(int cynpoly, unsigned int * * *croot,unsigned int *cy)
  {cexpitheta(cy);
   cstore(cyr);
   for (cykk=1;cykk<cynpoly+1;cykk++)
     {cmult(croot[cykk],cyr);
      cstore(croot[cykk]);
     }
   return 0;
  }

int crotaterealpoly(int cynpoly, unsigned int * *cyxrot,unsigned int *cy)
  {cexpitheta(cy);
   cstore(cyr);
   cconinttomp(1,0);
   cstore(cyrpr);
   for (cykk=0;cykk<cynpoly+1;cykk++)
     {cmultcbyr(cyrpr,cyxrot[cykk]);      //result is complex
      cstore(ccpolyout[cykk]);           //output
      cmult(cyr,cyrpr);
      cstore(cyrpr);
     }
   return 0;
  }