# include <stdio.h>
# define MSP    100000
# define  MM      200
# define  MV      1000
#define MDIM      100
# define MPR      2000
char opt[5];
int msp=MSP;
short mv=MV,mm=MM,mdim=MDIM,mpr=MPR,prime,dim,*spv,**spm,
      mspace[MSP],*vec[MV],**mat[MM],cp[500],pinv[MPR];
FILE *fopen(),*ip,*op;

main()
{ short i,j,l,m,n,*p,**dp,maxv,maxm,ord;
  int quot;
  char c,fault,flnm[80],f;
  printf("Do you wish to read matrices from a file(y/n)?  ");
  if ((c=getchar())=='y')
  { f=1; seeknln();
reenter:
    printf("Input filename:  "); scanf("%s",flnm);
    if ((ip=fopen(flnm,"r"))==0)
    { fprintf(stderr,"Cannot open %s.\n",flnm);
      seeknln(); goto reenter;
    }
    fscanf(ip,"%hd%hd%hd",&prime,&dim,&n);
    if (prime>mpr) {fprintf(stderr,"prime too big.\n"); exit(1);}
  }
  else
  {f=0;printf("Input modulus,dim:    "); scanf("%hd%hd",&prime,&dim);}
  if (prime>mpr) {fprintf(stderr,"modulus too big.\n"); exit(1);}
  if (dim>mdim) {fprintf(stderr,"dim too big.\n"); exit(1);}
  seeknln(); setpinv();
/* Set matrix and vector pointers as in mcd.c  */
  quot=msp/dim; if (quot>=mv) quot=mv-1; maxv=quot;
  for (i=1;i<=maxv;i++) vec[i]= mspace-1+(i-1)*dim;
  maxm= (maxv-10)/dim-1; if (maxm>=mm) maxm=mm-1;
  for (i=0;i<=maxm;i++) mat[i]=vec+10+i*dim;
  dp=mat[0]; p= *dp;
  for (i=1;i<=dim*dim;i++) *(++p)=0;
  for (i=1;i<=dim;i++) dp[i][i]=1;
  spm=mat[maxm]; spv=spm[1];
  printf("10 free vectors and %d matrices can be defined.\n",maxm-1);
  printf("Mat 0 is the identity.\n");
  if (f) { for (i=1;i<=n;i++) readmat(mat[i]);  fclose(ip); }
  while(1)
  { printf("Choose option. Print 'l' for list.\n");
    scanf("%s",opt);
    if (strcmp(opt,"l")==0)
    { seeknln();
      printf("rv n  or  rm n        = read vec or mat no. n.\n");
      printf("pv  n  or  pm  n      = write vec or mat no. n.\n");
      printf("tr m n                = calc transpose m[n] of m[m].\n");
      printf("inv m n               = calc inverse m[n] of m[m].\n");
      printf("im l m n              = calc image v[m] of v[l] under m[n].\n");
      printf("c  l m n              = calc comm  v[m] of v[l] under m[n].\n");
      printf("pr n l i(1)..i(l)     = calc prod m[n] of m[i(1)]...m[i(l)].\n");
      printf("sum n l i(1)...i(l)   = calc sum m[n] of m[i(1)]...m[i(l)].\n");
      printf("                        (mat nos may be negative for difference.)\n");
      printf("ch n i j k            = change m[n][i][j] to k.\n");
      printf("eq m n                = test equality of m[m],m[n].\n");
      printf("ord m n               = calc order of m[m] and put inverse in m[n].\n");
      printf("op filename           = output some matrices to filename.\n");
      printf("q                     = quit.\n");
    }
    else if (strcmp(opt,"rv")==0)
    { scanf("%hd",&n); seeknln();
      if (n<=0 || n>=10)
      { fprintf(stderr,"Invalid vector number %d\n",n); continue;}
      for (i=1;i<=dim;i++)
      { scanf("%hd",&m); m%=prime; if (m<0) m+=prime; vec[n][i]=m; }
      seeknln();
    }
    else if (strcmp(opt,"rm")==0)
    { scanf("%hd",&n); seeknln();
      if (n<=0 || n>=maxm)
      { fprintf(stderr,"Invalid matrix number %d\n",n); continue;}
      for (i=1;i<=dim;i++) for (j=1;j<=dim;j++)
      { scanf("%hd",&m); m%=prime; if (m<0) m+=prime; mat[n][i][j]=m; }
      seeknln();
    }
    else if (strcmp(opt,"pv")==0)
    { scanf("%hd",&n); seeknln();
      if (n<=0 || n>=10)
      { fprintf(stderr,"Invalid vector number %d\n",n); continue;}
      for (i=1;i<=dim;i++) printf("%3d",vec[n][i]); printf("\n");
    }
    else if (strcmp(opt,"pm")==0)
    { scanf("%hd",&n); seeknln();
      if (n<0 || n>=maxm)
      { fprintf(stderr,"Invalid matrix number %d\n",n); continue;}
      for (i=1;i<=dim;i++)
      { for (j=1;j<=dim;j++) printf("%3d",mat[n][i][j]); printf("\n"); }
    }
    else if (strcmp(opt,"tr")==0)
    { scanf("%hd%hd",&m,&n); seeknln();
      if (n<=0 || n>=maxm)
      { fprintf(stderr,"Invalid matrix number %d\n",n); continue;}
      if (m<0 || m>=maxm)
      { fprintf(stderr,"Invalid matrix number %d\n",m); continue;}
      trans(mat[m],mat[n]);
    }
    else if (strcmp(opt,"inv")==0)
    { scanf("%hd%hd",&m,&n); seeknln();
      if (n<=0 || n>=maxm)
      { fprintf(stderr,"Invalid matrix number %d\n",n); continue;}
      if (m<0 || m>=maxm)
      { fprintf(stderr,"Invalid matrix number %d\n",m); continue;}
      inv(mat[m],mat[n]);
    }
    else if (strcmp(opt,"im")==0 || strcmp(opt,"c")==0)
    { scanf("%hd%hd%hd",&l,&m,&n); seeknln();
      if (l<=0 || l>=10)
      { fprintf(stderr,"Invalid vector number %d\n",l); continue;}
      if (m<=0 || m>=10)
      { fprintf(stderr,"Invalid vector number %d\n",m); continue;}
      if (n<0 || n>=maxm)
      { fprintf(stderr,"Invalid matrix number %d\n",n); continue;}
      if (strcmp(opt,"im")==0) im(vec[l],vec[m],mat[n]);
      else comm(vec[l],vec[m],mat[n]);
    }
    else if (strcmp(opt,"pr")==0)
    { scanf("%hd%hd",&n,cp);
      if (n<=0 || n>=maxm)
      { fprintf(stderr,"Invalid matrix number %d\n",n); seeknln(); continue;}
      fault=0;
      for (i=1;i<=*cp;i++)
      { scanf("%hd",&m);
        if (m<0 || m>=maxm)
        { fprintf(stderr,"Invalid matrix number %d\n",m); fault=1; break;}
        cp[i]=m;
      }
      seeknln();
      if (fault) continue;
      prod(cp,mat[n]);
    }
    else if (strcmp(opt,"sum")==0)
    { scanf("%hd%hd",&n,cp);
      if (n<=0 || n>=maxm)
      { fprintf(stderr,"Invalid matrix number %d\n",n); seeknln(); continue;}
      fault=0;
      for (i=1;i<=*cp;i++)
      { scanf("%hd",&m);
        if (m<0 || m>=maxm)
        { fprintf(stderr,"Invalid matrix number %d\n",m); fault=1; break;}
        cp[i]=m;
      }
      seeknln();
      for (i=1;i<=dim;i++) for (j=1;j<=dim;j++)
      { m=0; for (l=1;l<= *cp;l++)
        if (cp[l]>0) m+=mat[cp[l]][i][j]; else m-=mat[-cp[l]][i][j];
        m%=prime; if (m<0) m+=prime; mat[n][i][j]=m;
      }
    }
    else if (strcmp(opt,"ch")==0)
    { scanf("%hd%hd%hd%hd",&n,&i,&j,&l); seeknln();
      if (n<=0 || n>=maxm)
      { fprintf(stderr,"Invalid matrix number %d\n",n); continue;}
      if (i<=0 || j<=0 || i>dim || j>dim)
      { fprintf(stderr,"Invalid row or column no %d,%d\n",i,j); continue;}
      mat[n][i][j]=l;
    }
    else if (strcmp(opt,"eq")==0)
    { scanf("%hd%hd",&m,&n); seeknln();
      if (n<0 || n>=maxm)
      { fprintf(stderr,"Invalid matrix number %d\n",n); continue;}
      if (m<0 || m>=maxm)
      { fprintf(stderr,"Invalid matrix number %d\n",m); continue;}
      f=1;
      for (i=1;i<=dim;i++)
      { for (j=1;j<=dim;j++) if (mat[m][i][j]!=mat[n][i][j]) { f=0; break; }
        if (f==0) break;
      }
      if (f) printf("Matrices are equal.\n");
      else printf("Matrices differ in place  %d,%d.\n",i,j);
    }
    else if (strcmp(opt,"ord")==0)
    { scanf("%hd%hd",&m,&n); seeknln();
      if (n<=0 || n>=maxm)
      { fprintf(stderr,"Invalid matrix number %d\n",n); continue;}
      if (m<0 || m>=maxm)
      { fprintf(stderr,"Invalid matrix number %d\n",m); continue;}
      *cp=1; cp[1]=m; prod(cp,mat[n]); prod (cp,spm); ord=1;
      while(1)
      { f=1;
        for (i=1;i<=dim;i++)
        { for (j=1;j<=dim;j++) if (mat[0][i][j]!=spm[i][j]) { f=0; break; }
          if (f==0) break;
        }
        if (f) { printf("Order=%d\n",ord); break;}
        *cp=1; cp[1]=maxm;  prod(cp,mat[n]);
        ord++; *cp=2; cp[1]=m; cp[2]=n; prod(cp,spm);
      }
    }
    else if (strcmp(opt,"op")==0)
    { scanf("%s",flnm); op=fopen(flnm,"w");
      printf("Input matrix nos to be saved preceded by their number\n");
      scanf("%hd",&n); fprintf(op,"%4d%4d%4d\n",prime,dim,n);
      fault=0;
      for (i=1;i<=n;i++)
      { scanf("%hd",&j);
        if (j<0 || j>=maxm)
        { fprintf(stderr,"Invalid matrix number %d\n",j); fault=1; break;}
        printmat(mat[j]);
      }
      seeknln();
      if (fault) { fclose(op); unlink(flnm); continue;}
      fclose(op);
    }
    else if (strcmp(opt,"q")==0) {seeknln(); break;}
    else { seeknln(); printf("Invalid option.\n");}
  }
  exit(0);
}

seeknln() { char c; while ((c=getchar())!='\n'); }

setpinv()
{ int i,j; int sum;
  for (i=0;i<prime;i++) pinv[i]=0;
  for (i=1;i<prime;i++) if (pinv[i]==0) for (j=1;j<prime;j++)
  if ((sum=i*j) % prime ==1) { pinv[i]=j; pinv[j]=i; break; }
}