# include <stdio.h>
extern short  npt,nf,cp[],orb[], pno[],fp[],tsv1[],tsv2[],tsv3[],orep[],
       *pptr[],mb;
short cpno;

addperm()
/* This computes the perm cp as perm no cpno, and adds it to the list pno.
   Externals: nf,cpno,fp,pno,pptr.
*/
{ short k,m,n;
  if (nf== -1)
  { fprintf(stderr,"Out of space. Increase PSP.\n"); return(-1);}
  cpno=nf; nf=fp[nf]; (*pno)++; pno[*pno]=cpno; m=cpno+1;
  for (n=1;n<=npt;n++) {k=image(n); pptr[cpno][n]=k;pptr[m][k]=n;}
  return(0);
}

intbase(pt,pos,stad,nbad,b,lorb,svptr) short pt,pos,*stad,*nbad,
        *b,*lorb,**svptr;
/* This changes base no pos to pt (unless it is already an earlier base point.
   b is the base, lorb, svptr as usual. nbad is tha address of nb (no of
   base points). It is assumed that the perm nos are linked by the extern
   array fp, starting at perm no *stad.
   Externals:npt,mb,tsv1,tsv2,tsv3,pno,pptr,svptr,orb,orep,nf,cp.
*/
{ char recalc,sk1,sk2;
  short i,j,k,l,m,n,v,z,b1,b2,lo,lo1,lo2,nlo1,nlo2,np1,np2,
        fpno,lpno,nr,stpos,npt1;
  npt1=npt+1; recalc=0; stpos=pos-1;
  if (pos== *nbad+1) stpos=0;
  else for (i=1;i<= *nbad;i++)
  { if (i==pos) stpos=0;
    if (b[i]==pt)
    { if (stpos!=0) {fprintf(stderr,"Base element there.\n"); return(1);}
      stpos=i;
    }
  }

/* If stpos=0, then bno is not at present a base point, so we introduce it
   as base[nb+1], and put stpos=nb+1. Otherwise, it is already present as
   base[stpos]. The idea is to swap base points l,l+1  for l=stpos-1,...pos
   successively.
*/
  if (stpos==0)
  { (*nbad)++;
    if (*nbad>mb)
    { fprintf(stderr,"nb too big. Increase SVSP (or MB).\n"); return(-1);}
    stpos= *nbad; b[*nbad]=pt; lorb[*nbad]=1;
    for (n=1;n<=npt;n++) svptr[*nbad][n]=0; svptr[*nbad][pt]= -1;
  }
  for (l=stpos-1;l>=pos;l--)
  { b1=b[l];b2=b[l+1]; lo1=lorb[l]; lo2=lorb[l+1]; b[l]=b2; b[l+1]=b1;
/* Now we list the permutations that may need to be changed */
    *pno=0; sk1=1; sk2=1;
    for (i= *stad;i!= -1;i=fp[i])
    { k=pptr[i][npt1];
      if (sk1 && k<=l+1)  {np1= *pno; sk1=0;}
      if (sk2 && k<=l)    {np2= *pno; sk2=0;}
      if (k<l) break; (*pno)++; pno[*pno]=i;
    }
    if (sk1) np1= *pno; if (sk2) np2= *pno;
    if (lo1==1)
    { for (n=1;n<=npt;n++) {svptr[l][n]=svptr[l+1][n];svptr[l+1][n]=0;}
      svptr[l+1][b1]= -1; lorb[l]=lo2; lorb[l+1]=1;
      for (i=np1+1;i<=np2;i++) pptr[pno[i]][npt1]=l;
      continue;
    }
    lo = orbitsv(b2,tsv1,0);
    if (lo==1)
    { for (n=1;n<=npt;n++) {svptr[l+1][n]=svptr[l][n];svptr[l][n]=0;}
      svptr[l][b2]= -1; lorb[l]=1; lorb[l+1]=lo1;
      for (i=np2+1;i<= *pno;i++) pptr[pno[i]][npt1]=l+1;
      continue;
    }
    nlo1=lo;lorb[l]=lo; nlo2= lo1*lo2/lo; lorb[l+1]=nlo2;
    if (nlo2==1)
    { for (n=1;n<=npt;n++) {svptr[l][n]=tsv1[n];svptr[l+1][n]=0;}
      svptr[l+1][b1]= -1;
      for (i=np1+1;i<=np2;i++) pptr[pno[i]][npt1]=l;
      continue;
    }
/* The easy cases have now been handled! */
    if (recalc) orbitsv(b1,svptr[l],0); else recalc=1;
    lpno=pno[*pno]; *pno=np2; orbitsv(b2,tsv2,0); *pno=np1;
    for (n=1;n<=npt;n++) tsv3[n]=1; nr=0;
    for (n=1;n<=npt;n++) if (svptr[l][n]>0 && tsv3[n])
    { nr++; orb[1]=n; orep[nr]=n; tsv3[n]=0; m=1;
      for (j=1;j<=m;j++)
      { z=orb[j]; for (k=1;k<= *pno;k++)
        { v=pptr[pno[k]][z];
          if (tsv3[v])
          { tsv3[v]=0; m++; orb[m]=v; }
        }
      }
    }
    for (n=1;n<=npt;n++) tsv3[n]=0;
    tsv3[b1]= -1; lo=1; orb[1]=b1; fpno= -1;
    for (i=1;i<=nr;i++)
    { *cp=0; addsv(orep[i],svptr[l]); j=image(b2); if (tsv2[j]!=0)
      { addsv(j,tsv2); j=image(b1); if (tsv3[j]==0)
        { if (fpno== -1) fpno=nf; else fp[cpno]=nf;
          if ((addperm())== -1) return(-1);
          pptr[cpno][npt1]=l+1; if ((lo=orbitsv(b1,tsv3,lo))==nlo2) break;
        }
      }
    }
    for (n=1;n<=npt;n++) {svptr[l+1][n]=tsv3[n]; tsv3[n]=0; }
    tsv3[b2] = -1; lo=1; orb[1]=b2;
    for (n=1;n<=npt;n++) if (tsv1[n]>0 && tsv3[n]==0)
    { *cp=0;
      addsv(n,tsv1); fp[cpno]=nf;
      if (addperm()== -1) return(-1);
      pptr[cpno][npt1]=l;
      if ((lo=orbitsv(b2,tsv3,lo))==nlo1) break;
    }
    for (n=1;n<=npt;n++) svptr[l][n]=tsv3[n];
    fp[cpno]=fp[lpno]; fp[lpno]=nf;
    if (np1==0) { nf= *stad; *stad=fpno; }
    else { j=pno[np1]; nf=fp[j]; fp[j]=fpno; }
  }
  while (lorb[*nbad]==1 && b[*nbad]!=pt) (*nbad)--;
  if (recalc && pos>1)
  { *pno=0; j=pos-1;
    for (i= *stad;i!= -1;i=fp[i])
    { k=pptr[i][npt1]; if (k<j) {orbitsv(b[j],svptr[j],0); j=k; }
      (*pno)++; pno[*pno]=i;
    }
    orbitsv(b[j],svptr[j],0);
  }
  return(0);
}