extern int seriesdeclarations();
extern int smult(int sn1,unsigned int * *sin1,unsigned int * *sin2);
extern int sadd(int sn1,unsigned int * *sin1,unsigned int * *sin2);
extern int ssub(int sn1,unsigned int * *sin1,unsigned int * *sin2);
extern int stonth(int sn3,unsigned int * *sin4,unsigned int *sp);
extern int smultfac(int sn2,unsigned int * *sin3,unsigned int *sfac);
extern int sstorit(int sn4,unsigned int * *sinstore,unsigned int * *sstore); //<
extern int sder(int sn5,unsigned int * *sin5);
extern int slogder(int sn6,unsigned int * *sin6);
//extern scontoc(int sn6, unsigned int * *sin6);  //moved to oversere.hc
extern int seval(int sn6, unsigned int * *sin6, unsigned int * s77);
extern int sinvert(int sn1,unsigned int * * sin1);            
extern int ssubstitute(int sn1, unsigned int * * sin1, unsigned int * * sin2);

int sn1,sn2,sn3,sn4,sn5,sn5snmax;

realscalar sp_ram;unsigned int *sp=sp_ram.s;
realscalar stemp_ram;unsigned int *stemp=stemp_ram.s;
realscalar stemp2_ram;unsigned int *stemp2=stemp2_ram.s;
realscalar s77_ram;unsigned int *s77=s77_ram.s;
realscalar sfac_ram;unsigned int *sfac=sfac_ram.s;

realvector sin1_ram;unsigned int * *sin1=sin1_ram.v;
realvector sin2_ram;unsigned int * *sin2=sin2_ram.v;
realvector sin3_ram;unsigned int * *sin3=sin3_ram.v;
realvector sin4_ram;unsigned int * *sin4=sin4_ram.v;
realvector sin5_ram;unsigned int * *sin5=sin5_ram.v;
realvector sin6_ram;unsigned int * *sin6=sin6_ram.v;
realvector sout_ram;unsigned int * *sout=sout_ram.v;
realvector sout1_ram;unsigned int * *sout1=sout1_ram.v;
realvector sX5_ram;unsigned int * *sX5=sX5_ram.v;
realvector sX6_ram;unsigned int * *sX6=sX6_ram.v;
realvector sstore_ram;unsigned int * *sstore=sstore_ram.v;
realvector sinstore_ram;unsigned int * *sinstore=sinstore_ram.v;
realvector sc_ram;unsigned int * *sc=sc_ram.v;
realvector sd_ram;unsigned int * *sd=sd_ram.v;

realmatrix sm1_ram;unsigned int * * *sm1=sm1_ram.m;
realmatrix sm2_ram;unsigned int * * *sm2=sm2_ram.m;

//unsigned int * * *soutc;  
//unsigned int * *spval;       //a complex number:: used for ?
/*
int siseriesdeclarations=seriesdeclarations();
seriesdeclarations()
  {int si22,si33;
   spval = new unsigned int * [2];      //complex variable
   for (si22=0;si22<2;si22++) spval[si22] = new unsigned int [zprecpad];
//   soutc = new unsigned int * * [zmaxsizepad];    //complex series
//   for (si22=0;si22<zmaxsizepad;si22++)
//     {soutc[si22] = new unsigned int * [2];          //complex series
//      for (si33=0;si33<2;si33++) soutc[si22][si33]= new unsigned int [zprecpad];
//     }
   return 0;
   }
*/

int sadd(int sn1,unsigned int * *sin1,unsigned int * *sin2) 
  {int si;
   for (si=0;si<sn1+1;si++)
     {zadd(sin1[si],sin2[si]);              //output is sout
      zstore(sout[si]);
     }
   return 0;
  }

int ssub(int sn1,unsigned int * *sin1,unsigned int * *sin2) 
  {int si;
   for (si=0;si<sn1+1;si++)
     {zsub(sin1[si],sin2[si]);
      zstore(sout[si]);                    //output is sout
     }
   return 0;
  }

int smult(int sn1,unsigned int * *sin1,unsigned int * *sin2)  
  {int si,sj;                          
   for (si=0;si<sn1+1;si++)
     {zzero(sX6[si]);
      for (sj=0;sj<si+1;sj++) 
        {zmult(sin1[sj],sin2[si-sj]);
         zadd(z,sX6[si]);
         zstore(sX6[si]);
        }                                     //output is sout
     }
   for (si=0;si<sn1+1;si++)
     {zrecall(sX6[si]); 
      zstore(sout[si]);
     }
   return 0;
  }

int smultfac(int sn2,unsigned int * *sin3,unsigned int *sfac)
  {int si;
   for (si=0;si<sn2+1;si++)
     {zmult(sin3[si],sfac);
      zstore(sout[si]);                      //output is sout
     }
   return 0;
  }

int sstorit(int sn4,unsigned int * *sinstore,unsigned int * *sstore)
  {int si;
   for (si=0;si<sn4+1;si++)
     {zrecall(sinstore[si]);
      zstore(sstore[si]);
     }
   return 0;
  }

int stonth(int sn3,unsigned int * *sin4,unsigned int *sp)
  {int si55,si5,sk;
   for (si55=0;si55<sn3+1;si55++)
     {zrecall(sin4[si55]);
      zstore(sX5[si55]);
     }
   if (sX5[0][prec+1]==0 || sX5[0][prec]==0)  
     {zmessage="first term in a series raised to a power must be >0";
      zbreak();
     }
   for (si55=1;si55<sn3+1;si55++) 
     {zdiv(sX5[si55],sX5[0]);
      zstore(sX5[si55]);
     }
   zlog(sX5[0]);
   zmult(z,sp);
   zexp(z);
   zstore(stemp);
   zzero(sX5[0]);                    //series is 0+sX5[1]*x+sX5[2]*x^2 etc
   for (si55=1;si55<sn3+1;si55++)
     {zzero(sout1[si55]);}
   zconinttomp(1);
   zstore(sout1[0]);                 //starts with 1
   for (si55=sn3;si55>=1;si55--)
     {sk=si55-1;
      zconinttomp(sk);
      zsub(sp,z);
      zdivint(z,si55);
      zstore(stemp2);
      smult(sn3,sout1,sX5);
      smultfac(sn3,sout,stemp2);             //output is sout
      sstorit(sn3,sout,sout1);
      zconinttomp(1);
      zstore(sout1[0]);
     }
   for (si55=0;si55<sn3+1;si55++)
     {zrecall(sout1[si55]);
      zmult(z,stemp);
      zstore(sout[si55]);}
   return 0;
  }

int sder(int sn5, unsigned int * *sin5)
  {int si;
   for (si=0;si<sn5;si++)
     {zmultint(sin5[si+1],si+1);             //output is sout
      zstore(sout[si]);
     }
   return 0;
  }

int slogder(int sn6, unsigned int * *sin6)
  {int si6;
   si6=sn6-1;
   sder(sn6,sin6);       //derivative of sin6
   sstorit(si6,sout,sc);
   zconinttomp(1);                       //output is sout
   zchs(z);
   zstore(sp);
   stonth(si6,sin6,sp);     //sd=sin6^-1
   sstorit(si6,sout,sd);
   smult(si6,sc,sd);    //this should put answer in sout
   return 0;
  }

int seval(int sn6, unsigned int * *sin6, unsigned int * s77)
  {int si6;
   zrecall(s77);
   zstore(stemp2);
   zzero(stemp);
   for (si6=sn6;si6>=0;si6--)
     {zmult(stemp,stemp2);
      zadd(z,sin6[si6]);
      zstore(stemp);                   //output is z
      }
   zrecall(stemp);
   return 0;
  }

int sinvert(int sn1,unsigned int * * sin1)            //nmax,input series=sn1
  {int si,sj,sk;   //input series assumed to start with x
                   //as in log(1-x). See note at end.
   for (si=1;si<sn1+1;si++) 
     {for (sj=1;sj<sn1+1;sj++) zzero(sm1[si][sj]);}

   for (si=1;si<sn1+1;si++)
     {zrecall(sin1[si]);
      zstore(sm1[1][si]);              //use sm1 for matrix x
     }

   for (si=2;si<sn1+1;si++)                
     {cerr<<".";
      for (sk=si;sk<sn1+1;sk++)
       {zzero(sm1[si][sk]);            //x[i][k]=0
        for (sj=1;sj<sk-si+2;sj++)        //x(i,k)=x(i,k)+x(i-1,k-j)*x(1,j) 
          {zmult(sm1[si-1][sk-sj],sm1[1][sj]);
           zadd(z,sm1[si][sk]);
           zstore(sm1[si][sk]);
          }
       }
     }
   cerr<<'\n';

   for (si=1;si<sn1+1;si++) 
     {for (sj=1;sj<sn1+1;sj++) zzero(sm2[si][sj]);}

   for (si=sn1;si>=1;si--)    //for i=nmax to 1 step -1
     {cerr<<".";
      zrecip(sm1[si][si]);      //c(i,i)=1/x(i,i)
      zstore(sm2[si][si]);      //use sm2 for matrix c
      for (sk=si+1;sk<sn1+1;sk++)               //for k=i+1 to nmax
        {zzero(sm2[si][sk]);                   //c(i,k)=0
         for (sj=si+1;sj<sn1+1;sj++)           //for j=i+1 to nmax
           {zmult(sm1[si][sj],sm2[sj][sk]);   //c(i,k)=c(i,k)+x(i,j)*c(j,k)
            zadd(z,sm2[si][sk]);
            zstore(sm2[si][sk]);
           }
         zdiv(sm2[si][sk],sm1[si][si]);   //c(i,k)=-c(i,k)/x(i,i)
         zchs(z);                          
         zstore(sm2[si][sk]);
        }
     }
     cerr<<'\n';

   zzero(sout[0]);  //added October 20,99--the inverse series must
                    //also start with zero.
   for (si=1;si<sn1+1;si++) 
      {zrecall(sm2[1][si]);       //inverted series   //output is sout
       zstore(sout[si]);
      }
   cerr<<"inversion complete"<<'\n';
   return 0;
  }                                                  
                                                   
int ssubstitute(int sn1, unsigned int * * sin1, unsigned int * * sin2)
  {int si,sj,sk;                 //sin1 into sin2
                                 //this assumes both start with x^1
                                 //really, only sin1 needs to start
for (si=1;si<sn1+1;si++)         //with x
  {for (sj=1;sj<sn1+1;sj++) zzero(sm1[si][sj]);}

   for (si=1;si<sn1+1;si++) zzero(sout1[si]);

   for (si=1;si<sn1+1;si++)
     {zrecall(sin1[si]);
      zstore(sm1[1][si]);              //use sm1 for matrix x
      zrecall(sin2[si]);
      zstore(sout1[si]);
     }
   for (si=2;si<sn1+1;si++)                
     {cerr<<".";
      for (sk=si;sk<sn1+1;sk++)
       {zzero(sm1[si][sk]);            //x[i][k]=0
        for (sj=1;sj<sk-si+2;sj++)        //x(i,k)=x(i,k)+x(i-1,k-j)*x(1,j) 
          {zmult(sm1[si-1][sk-sj],sm1[1][sj]);
           zadd(z,sm1[si][sk]);
           zstore(sm1[si][sk]);
          }
       }
     }
   cerr<<'\n';
   sout[0]=sin2[0];    //added Oct.30, 1999 so that sin2 does not
                       //need to begin with x^1.
   for (si=1;si<sn1+1;si++)
     {zzero(sout[si]);        //why?
      for (sj=1;sj<si+1;sj++)                    //output is sout
        {zmult(sout1[sj],sm1[sj][si]);
         zadd(z,sout[si]);
         zstore(sout[si]);
        }
     }
   return 0;
  }