/**************************************************************************
*                  Symbolic Parser with Polynomials                       * 
*                 Example : '(A+B)^2' ==> A^2+2AB+B^2                     *
*             C++ release 2.21 Copyright J-P Moreau, Paris                *
*                 (translated from pascal 2.21 release)                   *
* ----------------------------------------------------------------------- *
*      Release 2.0: program improvement to take care of two digits        *
*                   exponants.                                            *
*      Release 2.1: (A+B)(C+D) ==> (A+B)*(C+D)                            *
*                    25(C-3D)  ==>  25*(C-3D)                             *
*      Release 2.2: 1. Some procedures are transformed into               *
*                   functions to return an index value, in case of        *
*                   table overflows ( SIMPLI, PARE, PUIS, COPIE,          *
*                   IMPR ), that allows to stop the computing             *
*                   without crashing the program.                         *
*                   2. Possibility to read the input string from          *
*                   a text file.                                          *
*     Release 2.21: Minor corrections (Display CHs > 70 characters        * 
*                   and bug in IMPR, line 392).                           *
* ----------------------------------------------------------------------- *
* SAMPLE RUN:                                                             *
*                                                                         *
*                   SYMBOLIC PARSER for POLYNOMIALS                       *
*                                                                         *
*                   Example: (A+B)^2 ==> A^2+2AB+B^2                      *
*                                                                         *
* Inputs (0=SCREEN  1=algebra.dat file): 0                                *
* Input string to evaluate:                                               *
* (A+B+C)^3                                                               *
*                                                                         *
* Detail analysis (0=NO  1=YES): 1                                        *
* Outputs (0=SCREEN  1=algebra.lst file): 0                               *
*                                                                         *
* SIMPLIFICATION:                                                         *
* (A+B+C)^3                                                               *
* POWER:                                                                  *
* (A^3+3A^2B+3A^2C+3AB^2+6ABC+3AC^2+B^3+3B^2C+3BC^2+C^3)                  *
* ADDITION:                                                               *
* A^3+3A^2B+3A^2C+3AB^2+6ABC+3AC^2+B^3+3B^2C+3BC^2+C^3                    *
*                                                                         *
* FINAL SIMPLIFICATION AND RESULT:                                        *
*                                                                         *
* A^3+3A^2B+3A^2C+3AB^2+6ABC+3AC^2+B^3+3B^2C+3BC^2+C^3                    *
*                                                                         *
* Evaluate another string (0=NO  1=YES): 1                                *
*                                                                         *
* (Assuming the input file algabra.dat contains the line:                 *
* (A+B+C+D)^4+(A+B+C-D)^4 )                                               *
*                                                                         *
* Inputs (0=SCREEN  1=algebra.dat file): 1                                *
* Input string to evaluate:                                               *
* (A+B+C+D)^4+(A+B+C-D)^4                                                 *
*                                                                         *
* Detail analysis (0=NO  1=YES): 1                                        *
* Outputs (0=SCREEN  1=algebra.lst file): 1                               *
*                                                                         *
* The file algebra.lst contains:                                          *
*                                                                         *
* STRING TO EVALUATE:                                                     *
*  (A+B+C+D)^4+(A+B+C-D)^4                                                *
*                                                                         *
* SIMPLIFICATION:                                                         *
*  (A+B+C+D)^4+(A+B+C-D)^4                                                *                                               
* POWER:                                                                  *
*  (A^4+4A^3B+4A^3C+4A^3D+6A^2B^2+12A^2BC+12A^2BD+6A^2C^2+12A^2CD+6A^2D^2 *
*  +4AB^3+12AB^2C+12AB^2D+12ABC^2+24ABCD+12ABD^2+4AC^3+12AC^2D+12ACD^2+4A *
*  D^3+B^4+4B^3C+4B^3D+6B^2C^2+12B^2CD+6B^2D^2+4BC^3+12BC^2D+12BCD^2+4BD^ *
*  3+C^4+4C^3D+6C^2D^2+4CD^3+D^4)+(A+B+C-D)^4                             *                             
* ADDITION:                                                               *
*  A^4+4A^3B+4A^3C+4A^3D+6A^2B^2+12A^2BC+12A^2BD+6A^2C^2+12A^2CD+6A^2D^2+ *
*  4AB^3+12AB^2C+12AB^2D+12ABC^2+24ABCD+12ABD^2+4AC^3+12AC^2D+12ACD^2+4AD *
*  ^3+B^4+4B^3C+4B^3D+6B^2C^2+12B^2CD+6B^2D^2+4BC^3+12BC^2D+12BCD^2+4BD^3 *
*  +C^4+4C^3D+6C^2D^2+4CD^3+D^4+(A+B+C-D)^4                               *                              
* POWER:                                                                  *
*  A^4+4A^3B+4A^3C+4A^3D+6A^2B^2+12A^2BC+12A^2BD+6A^2C^2+12A^2CD+6A^2D^2+ *
*  4AB^3+12AB^2C+12AB^2D+12ABC^2+24ABCD+12ABD^2+4AC^3+12AC^2D+12ACD^2+4AD *
*  ^3+B^4+4B^3C+4B^3D+6B^2C^2+12B^2CD+6B^2D^2+4BC^3+12BC^2D+12BCD^2+4BD^3 *
*  +C^4+4C^3D+6C^2D^2+4CD^3+D^4+(A^4+4A^3B+4A^3C-4A^3D+6A^2B^2+12A^2BC-12 *
*  A^2BD+6A^2C^2-12A^2CD+6A^2D^2+4AB^3+12AB^2C-12AB^2D+12ABC^2-24ABCD+12A *
*  BD^2+4AC^3-12AC^2D+12ACD^2-4AD^3+B^4+4B^3C-4B^3D+6B^2C^2-12B^2CD+6B^2D *
*  ^2+4BC^3-12BC^2D+12BCD^2-4BD^3+C^4-4C^3D+6C^2D^2-4CD^3+D^4)            *           
* ADDITION:                                                               *
*  A^4+4A^3B+4A^3C+4A^3D+6A^2B^2+12A^2BC+12A^2BD+6A^2C^2+12A^2CD+6A^2D^2+ *
*  4AB^3+12AB^2C+12AB^2D+12ABC^2+24ABCD+12ABD^2+4AC^3+12AC^2D+12ACD^2+4AD *
*  ^3+B^4+4B^3C+4B^3D+6B^2C^2+12B^2CD+6B^2D^2+4BC^3+12BC^2D+12BCD^2+4BD^3 *
*  +C^4+4C^3D+6C^2D^2+4CD^3+D^4+A^4+4A^3B+4A^3C-4A^3D+6A^2B^2+12A^2BC-12A *
*  ^2BD+6A^2C^2-12A^2CD+6A^2D^2+4AB^3+12AB^2C-12AB^2D+12ABC^2-24ABCD+12AB *
*  D^2+4AC^3-12AC^2D+12ACD^2-4AD^3+B^4+4B^3C-4B^3D+6B^2C^2-12B^2CD+6B^2D^ *
*  2+4BC^3-12BC^2D+12BCD^2-4BD^3+C^4-4C^3D+6C^2D^2-4CD^3+D^4              *             
*                                                                         *
* FINAL SIMPLIFICATION AND RESULT:                                        *
*                                                                         *
*  2A^4+8A^3B+8A^3C+12A^2B^2+24A^2BC+12A^2C^2+12A^2D^2+8AB^3+24AB^2C+24AB *
*  C^2+24ABD^2+8AC^3+24ACD^2+2B^4+8B^3C+12B^2C^2+12B^2D^2+8BC^3+24BCD^2+2 *
*  C^4+12C^2D^2+2D^4                                                      *
*                                                                         *
* ----------------------------------------------------------------------- *
* Reference: "Calcul symbolique et informatique. Du calcul numérique au   *
*             calcul littéral (programmes en BASIC de A. DESHAYES),       *
*             Masson Paris, 1985" [BIBLI 02].                             *
*                                                                         * 
*             C++ Version By Jean-Pierre Moreau, Paris                    *
*                       (www.jpmoreau.fr)                                 * 
**************************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>

#define  MAXSIZE  270                 //max number of terms in development
#define  FALSE  0
#define  TRUE   1

//Labels inherited from Basic (in main):
//debut,e130,e220,e250,e290,e300,e340,e350,e370,e460,e470,e490,e510,e540,e560,
//      e1000,e1465,e1470,e1580,e8000,res,fin

typedef char Str125[126];
typedef char Str20[21];                             //max strlen of a term
typedef char Str1[2];

//global variables
  Str125 CHs;                                         //string to evaluate
  char   RRs[770];                                          //final answer
  Str1   ALFA[27];                               //table of variables A..Z
                                                      //(index 0 not used)   
  Str20 CPs[MAXSIZE];
  Str20 Rs[MAXSIZE], Ps[MAXSIZE],SIs[MAXSIZE];
  float T[MAXSIZE],P[MAXSIZE],R[MAXSIZE];
  Str20 Es[27];
  int   IE[27];
  float CP[MAXSIZE] ,SI[MAXSIZE];
  int   IPI[MAXSIZE],IPO[MAXSIZE];
      
  Str20 Ds,ECs,IEs,Ts,Ys,PRs,NBs,ENs,VLs,Ms,M1s;
  int   I,IE1,IEN,IEX,INB,ISI,IVL,LE1,NB,NC,NL,RE;
  int   IAP,IDI,IDM,IFM,IH,IHC,IHO,IHP,INO,IPF,IPP,JPO,LE;
  int   I4,IDC,IMD,IMG,IO1,IR,ISL,KS,NS,RF,RT;
  unsigned int IC, IC1;

  char  Cs,SUs;
  float EC,EX,M;
  char  E1s[3];
  FILE  *F1;      //pointer to input/output text file
  bool  Ok;	      //if Not ok, something went wrong


  bool Verif(int N) {
  //return TRUE if ( index N <= MAXSIZE
    if (N > MAXSIZE) return FALSE;
    else return TRUE;
  } //Verif()

  void Message() { 
  // Message given if a table index > MAXSIZE
  // To be used after a call to Verif().
    if (RT==0)
      printf("  Too many terms!\n");
    else {
      fprintf(F1,"  Too many terms!\n");
      fclose(F1);
    }
  }

  int SIMPLI(int IS1) {
  // simplify polynomial (put together similar monoms)
    int IS2,NB;
    if (!Verif(IS1)) return IS1;
    if (IS1==1) {NS=1; return 0;} 
    IS2=IS1-1;
    for (NB=1; NB<=IS2; NB++)
      if (SI[NB] != 0.0)
        for (KS=NB+1; KS<=IS1; KS++)
	  if (strcmp(SIs[KS],SIs[NB])==0) {
            SI[NB]=SI[NB]+SI[KS];
            SI[KS]=0.0;
          }

    NS=0;
    for (I=1; I<=IS1; I++)
      if (SI[I] != 0.0) {
        NS++;
        if (!Verif(NS)) return NS;
        strcpy(SIs[NS],SIs[I]);
        SI[NS]=SI[I];
      }

    if (NS==0) {
      NS=1;
      strcpy(SIs[1],"");
    }
    return 0;  //all ok
  } //Simpli()

  void MONOME(Str20 Ts)  {
  //Simplify monom ex.: BAB ==> AB^2
  //Labels: e80, e200
    int I,IVL; unsigned int I2;
    Str20 IEs,SSs;
    char  SCs,Tc1,Tc2;
    for (I=0; I<27; I++) IE[I]=0;
    I2=0;
    SCs=Ts[I2];
e80:if (I2 >= strlen(Ts)) goto e200;
    for (I=1; I<27; I++)
      if (SCs==ALFA[I][0])  NL=I;
    IE[NL]++; I2++;
    SCs=Ts[I2];

    if ((SCs=='^') && (I2<strlen(Ts))) {
      I2++;
      //modif. for release 2.0
      Tc1=Ts[I2]; Tc2=Ts[I2+1];
      if ((Tc2>='0') && (Tc2<='9') && (I2<strlen(Ts)-1)) {
        //Pascal: SSs=Ts[I2]+Ts[I2+1]
	SSs[0]=Tc1; SSs[1]=Tc2; SSs[2]='\0';
      } 
      else {
        //Pascal: SSs=Ts[I2];
        SSs[0]=Tc1; SSs[1]='\0';
      }
      IVL=atoi(SSs);
      IE[NL]=IE[NL]+IVL-1;
      I2 += strlen(SSs);   
      SCs=Ts[I2];
    }
    goto e80;

e200:strcpy(Ts,"");
    for (I2=1; I2<27; I2++) {
      if (IE[I2] != 0) {
        strcat(Ts,ALFA[I2]);
        if (IE[I2] > 1) {
          sprintf(IEs,"%d",IE[I2]);
          //Pascal: Ts=Ts+'^'+IEs
	  strcat(Ts,"^"); strcat(Ts,IEs);
        }
        IE[I2]=0;
      }
    }
  } //Monome()

  int COPIE()  { //copy a monom into another one
				 //used by PARE() and PUIS()
  //Labels: e9620,e9840
    int I,IHR,IZZ;
    float CO;
    if (IDC <= 0) {
      for (I=NS; I>0; I--) {
        IR++;
        strcpy(CPs[IR],SIs[I]);
        CP[IR]=SI[I];
      }
      JPO=IR+1;
      strcpy(CPs[JPO],"(");
      CP[JPO]=(float) IPF;
      CP[IPF]=(float) JPO;
      IDM=JPO-1;
      IAP=IO1+1;
e9620:if (strcmp(CPs[IAP],"$")==0) {
        IAP++;
        goto e9620;
      }
      if (IO1==IDI)  IDI += IDC;
      if (IDC==0) return 0;
      for (I=JPO+1; I<=IO1; I++) {
        strcpy(CPs[I],"$");
        CP[I]=0.0;
      }
      return 0;
    }

    if (IO1!=IDI) {
      IZZ=IO1+1;
      IHR=IHC-1;
      for (I=1; I<=IHR; I++) {
        CO=CP[IPI[I]];
        if (CO > JPO)  CP[IPI[I]]=CO+IDC;
      }
      I=IDI+IDC;
      if (!Verif(I)) return I;
      for (I=IDI; I>=IZZ; I--) {
        strcpy(CPs[I+IDC],CPs[I]);
        if (strcmp(CPs[I],")")==0) {
          CP[I+IDC]=CP[I]+IDC;
          goto e9840;
        }
        if ((strcmp(CPs[I],"(")==0) && (CP[I] > IPF)) {
          CP[I+IDC]=CP[I]+IDC;
          goto e9840;
        }
        CP[I+IDC]=CP[I];
e9840:;} // I loop
	} //if IO1<>IDI

    IDI += IDC;
    if (!Verif(NS)) return NS;
    for (I=NS; I>0; I--) {
      IR++;
      strcpy(CPs[IR],SIs[I]);
      CP[IR]=SI[I];
    }

    JPO=IR+1;
    strcpy(CPs[JPO],"(");
    CP[JPO]=(float) IPF;
    CP[IPF]=(float) JPO;
    IDM=JPO-1;
    IAP=JPO+1;
    while (strcmp(CPs[IAP],"$")==0)  IAP++;
    return 0;
  } // Copie()

  //utility functions used by IMPR()
  //print current result line to screen
  void print_line(int i1, int i2, char *line) {
  int i;
    printf("  ");
    for (i=i1; i<i2; i++) printf("%c",line[i]); printf("%c\n",line[i2]);
  }

  //print current result line to output file
  void fprint_line(int i1, int i2, char *line) {
  int i;
    fprintf(F1,"  ");  
    for (i=i1; i<i2; i++) fprintf(F1,"%c",line[i]); fprintf(F1,"%c\n",line[i2]);
  }

  bool IMPR() {
  //build up and print string RRs of intermediate or final result
  //returns FALSE if result string is too big !
  //Label: e500;
    int  I,IMO,INO,ISP,LR,LR1,LR2,icp;
    char Zs[11];
    ISP=0;
    //Empty work string RRs
	strcpy(RRs,"");

    /*for debugging only
    for (I=IDI; I>2; I--) {
      printf("  i=%d  CP$=%s  CP=%f\n", I,CPs[I],CP[I]);
    }
    scanf("%s", rep); */

    //build up result string RRs from CPs and CP tables
    //CP float values are rounded here to integers (this can be modified by user)
    for (I=IDI; I>2; I--) {
      if (strcmp(CPs[I],"$")==0) goto e500;
      IMO=0; INO=0;
      if ((CPs[I][0]>='A') && (CPs[I][0]<='Z'))  IMO=1;
      if (strlen(CPs[I])<1)  INO=1;
      
      if (ISP==1) {
        if (((IMO==1) || (INO==1)) && (CP[I]>0.0) && (strlen(RRs)>0))
          strcat(RRs,"+");
        if (strcmp(CPs[I],"(")==0) strcat(RRs,"+");
      }

      if (((IMO==1) || (INO==1)) && (floor(fabs(CP[I]))!=1) && (CP[I]!=0.0)) {
        //Pascal: Str(Round(CP[I]),Zs);
        icp=(int) floor(CP[I]);
        sprintf(Zs,"%d", icp);
        strcat(RRs,Zs);
      }

      if ((floor(CP[I])==-1) && (strlen(CPs[I])<1))
        strcat(RRs,"-1");

      if ((floor(CP[I])==-1) && (strlen(CPs[I])>0))
        strcat(RRs,"-");

      if ((floor(CP[I])==1) && (strlen(CPs[I])<1))
        strcat(RRs,"1");

      if (CP[I]!=0.0) strcat(RRs,CPs[I]);

      //only for intermediate results
      //Pascal: if ((CP[I]==0.0) && (CPs[I][1] in ['*','+','-']) then
      if ((CP[I]==0.0) && ((CPs[I][0]=='*') || (CPs[I][0]=='+') || (CPs[I][0]=='-')))
        strcat(RRs,CPs[I]);

      ISP=0;
      if ((IMO==1) || (INO==1) || (strcmp(CPs[I],")")==0))  ISP=1;
e500:;} //I loop

    LR=strlen(RRs);
    //fill eventual uncomplete line with bkanks
    //to avoid unwanted ending characters!
    if (LR<70) {
      LR2=69-LR;
      for (I=LR; I<=LR+LR2; I++) RRs[I]=' ';
    }
    else if (LR<770) {
      LR1=LR / 70;
      LR2=LR-LR1*70;
      for (I=LR; I<=LR+69-LR2; I++)  RRs[I]=' ';
    }

    //now print final result string RRs ! 
    if (RT==0) {  // case print to screen
      if (LR==0) printf("  0");
      else print_line(0,69,RRs);
      if (LR > 69) print_line(70,70+69,RRs); 
      if (LR > 139) print_line(140,140+69,RRs); 
      if (LR > 209) print_line(210,210+69,RRs);
      if (LR > 279) print_line(280,280+69,RRs);
      if (LR > 349) print_line(350,350+69,RRs);
      if (LR > 419) print_line(420,420+69,RRs);
      if (LR > 489) print_line(490,490+69,RRs);
      if (LR > 559) print_line(560,560+69,RRs);
      if (LR > 629) print_line(630,630+69,RRs);
      if (LR > 699) print_line(700,700+69,RRs);
      if (LR > 769) {
        printf("  RESULT STRING TOO BIG (>%d) !\n", LR);
        return FALSE;
      }
    }
    else  {  // case print to output file ALGEBRA.LST
      if (LR==0) fprintf(F1,"  0");
      else fprint_line(0,69,RRs);
      if (LR > 69) fprint_line(70,70+69,RRs); 
      if (LR > 139) fprint_line(140,140+69,RRs); 
      if (LR > 209) fprint_line(210,210+69,RRs);
      if (LR > 279) fprint_line(280,280+69,RRs);
      if (LR > 349) fprint_line(350,350+69,RRs);
      if (LR > 419) fprint_line(420,420+69,RRs);
      if (LR > 489) fprint_line(490,490+69,RRs);
      if (LR > 559) fprint_line(560,560+69,RRs);
      if (LR > 629) fprint_line(630,630+69,RRs);
      if (LR > 699) fprint_line(700,700+69,RRs);
      if (LR > 769) {
        fprintf(F1,"  RESULT STRING TOO BIG (>%d) !\n", LR);
        return FALSE;
      }
    }
    return TRUE; //print ok.
  }

  int PARE() {
  //MULTIPLY PARENTHESES
  //Labels: e5010,e5054,e5075
    int ID1,ID2,IF1,IF2,IO2,IS1,IUM,L1,L2,N1,N2;
    int I,J,MA,MB,NM;
    Ok=TRUE;
    for (I=1; I<=MAXSIZE; I++) strcpy(SIs[I],"");
    IO2=JPO;
    IF2=IPF;
    IF1=IAP+1;
e5010:if ((strcmp(CPs[IF1],"$")==0) && (IF1 < IDI)) {
      IF1++;
      goto e5010;
    }
    IO1=(int) floor(CP[IF1]);
    ID1=IO1-1;
    L1=IF1+1;
    ID2=IO2-1;
    L2=IF2+1;
    if (ID1==L1) {
      N1=1;
      goto e5054;
    }
    IUM=0;

    for (I=ID1; I>=L1; I--) {
      IUM++;
      if (!Verif(IUM)) return IUM;
      strcpy(SIs[IUM],CPs[I]);
      SI[IUM]=CP[I];
    }

    IS1=IUM;
    I=SIMPLI(IS1);
    if (I>0) return I;
    N1=NS;
    IUM=L1-1;
    for (I=N1; I>0; I--) {
      IUM++;
      if (!Verif(IUM)) return IUM;
      strcpy(CPs[IUM],SIs[I]);
      CP[IUM]=SI[I];
    }

e5054:if (ID2==L2) {
      N2=1;
      goto e5075;
	}
    IUM=0;
    for (I=ID2; I>=L2; I--) {
      IUM++;
      if (!Verif(IUM)) return IUM;
      strcpy(SIs[IUM],CPs[I]);
      SI[IUM]=CP[I];
    }

    IS1=IUM;
    I=SIMPLI(IS1);
    if (I>0) return I;
    N2=NS;
    IUM=L2-1;
    for (I=N2; I>0; I--) {
      IUM++;
      if (!Verif(IUM)) return IUM;
      strcpy(CPs[IUM],SIs[I]);
      CP[IUM]=SI[I];
    }

e5075:MA=N1+L1-1;
    MB=N2+L2-1;
    NM=0;
    for (I=MA; I>=L1; I--)
      for (J=MB; J>=L2; J--) {
        //Pascal: Ts=CPs[I]+CPs[J];
        strcpy(Ts,CPs[I]); strcat(Ts,CPs[J]);

        MONOME(Ts); //simplify monom Ts if necessary

        NM++;
        if (!Verif(NM)) return NM;
        strcpy(SIs[NM],Ts);
        SI[NM]=CP[I]*CP[J];
      }

    IS1=NM;
    I=SIMPLI(IS1);
    if (I>0) return I;
    IR=IF2;
    IDC=NS-IO1+L2;

    I=COPIE();
    if (I>0) return I;
    if (RE==1) {
      if (RT==0)
        printf("  PARENTHESES MULTIPLICATION:\n");
      else
        fprintf(F1,"  PARENTHESES MULTIPLICATION:\n");

      if (!IMPR()) { //string already too big ?
        Ok=FALSE;
        return -1;
      }
    }
    return 0;
  } // PARE()

  int PUIS() {
  // Integer power of a polynomial (max. 2 digits power) 
  // Labels: e9080,e9480,e9485
    int I,IS1,ITP,IU,IUP,IX1,IZP,JP,KP;
    Ok=TRUE;
    //current exponant is put in EX
    if (EC==0) {
      EX=CP[IPP-1];
      IO1=JPO;
    }
    if (EX==0.0) {
      strcpy(SIs[1],"");
      SI[1]=1.0;
      NS=1;
      goto e9480;
    }
    if (EX==1.0) {
      if (EC==0.0) goto e9485;
      IUP=0;
      for (I=IDM; I>=IFM; I--) {
        IUP++;
        if (!Verif(IUP)) return IUP;
        strcpy(SIs[IUP],CPs[I]);
        SI[IUP]=CP[I];
      }
      NS=IDM-IFM+1;

      goto e9480;
    }

    ITP=IDM-IFM+1;
    if (ITP==1) {
      NS=1;
      strcpy(Ps[1],CPs[IFM]);
      P[1]=CP[IFM];
      goto e9080;
    }
    I=0;
    for (JP=IDM; JP>=IFM; JP--)  {
      I++;
      if (!Verif(I)) return I;
      strcpy(SIs[I],CPs[JP]);
      SI[I]=CP[JP];
    }

    IS1=ITP;
    I=SIMPLI(IS1);
    if (I>0) return I;
    
    IX1=IFM+NS-1;
    IUP=IX1+1;
    if (!Verif(IUP)) return IUP;
    for (I=1; I<=NS; I++) {
      IUP--;
      strcpy(CPs[IUP],SIs[I]);
      CP[IUP]=SI[I];
      strcpy(Ps[I],SIs[I]);
      P[I]=SI[I];
    }

e9080:IEX=(int) floor(EX)-1;
    for (KP=1; KP<=IEX; KP++) {
      IZP=0;
      for (I=IX1; I>=IFM; I--) 
        for (IU=1; IU<=NS; IU++) {
          IZP++;
          if (!Verif(IZP)) return IZP;
          //Pascal: Ts=CPs[I]+Ps[IU];
	  strcpy(Ts,CPs[I]); strcat(Ts,Ps[IU]);

          MONOME(Ts);   //simplify monom Ts if necessary

          strcpy(SIs[IZP],Ts);
          SI[IZP]=CP[I]*P[IU];
        }

      IS1=IZP;
      I=SIMPLI(IS1);
      if (I>0) return I;

      if (!Verif(NS)) return NS;
      for (I=1; I<=NS; I++) {
        strcpy(Ps[I],SIs[I]);
        P[I]=SI[I];
      }    
    } // end of KP loop

e9480:IDC=NS-IO1+IPF+1;

    IR=IPF;
    COPIE();

    if (EC==1.0) IFM=IPF+1;
e9485:if (EC==0.0) {
      strcpy(CPs[IPP],"$");
      CP[IPP]=0.0;
      strcpy(CPs[IPP-1],"$");
      CP[IPP-1]=0.0;
      IPP=IPP-2;
    }

    if (RE==1) {
      if (RT==0) printf("  POWER:\n");
      else fprintf(F1,"  POWER:\n");
      if (!IMPR()) { //string already too big ?
        Ok=FALSE;
        return -1;
      }
    }
    EC=0.0;
    return 0; //power ok.
  } //PUIS()

  void ADDI() {
  // suppress matching parentheses with a preceding + sign or no +/- sign
    CP[JPO]=0.0;
    CP[IPF]=0.0;
    strcpy(CPs[JPO],"$");
    strcpy(CPs[IPF],"$");
    if (strcmp(CPs[IAP],"+")==0) {
      strcpy(CPs[IAP],"$");
      CP[IAP]=0.0;
    }
    if (RE==1) {
      if (RT==0)
        printf("  ADDITION:\n");
      else
        fprintf(F1,"  ADDITION:\n");
      IMPR();
    }
  } // Addi()

  void SOUS() {
  // suppress matching parentheses with a preceding - sign
      for (I=IDM; I>=IFM; I--)  CP[I]=-CP[I];
      CP[JPO]=0.0;
      CP[IPF]=0.0;
      strcpy(CPs[JPO],"$");
      strcpy(CPs[IPF],"$");
      strcpy(CPs[IAP],"$");
      CP[IAP]=0.0;
      if (RE==1) {
        if (RT==0)
          printf("  SUBSTRACTION:\n");
        else
          fprintf(F1,"  SUBSTRACTION:\n");
        IMPR();
      }
  } // Sous()

  void OPERA(Str125 CHs, unsigned int *IC, char *Cs, int *NC, Str20 Ds) {
    char s[2];
  // Look for a derivative £ operator (anyway not implemented here)
  // Label e670,e680,e700;
	if (*Cs=='£') {
e670: strcpy(Ds,"");
e680: if (*Cs != '(') {
        //Pascal: Ds=Ds+Cs;
        s[0]=*Cs; s[1]='\0';
	strcat(Ds,s);
        *IC=*IC+1;
        *Cs=CHs[*IC];
        goto e680;
      }
      *NC=*NC+1;
      strcpy(CPs[*NC],Ds);
    }
e700:if (*Cs=='(') {
      *NC=*NC+1;
      strcpy(CPs[*NC],"(");
      *IC=*IC+1;
      *Cs=CHs[*IC];
      goto e700;
    }
    if (*Cs=='£') goto e670;
  } //Opera()

  void DERI() {
    // derivation not implemented
  }

  //initialize table of uppercase latters (index 0 not used)
  void Fill_ALFA() {
    strcpy(ALFA[1],"A"); strcpy(ALFA[2],"B"); strcpy(ALFA[3],"C"); strcpy(ALFA[4],"D");
    strcpy(ALFA[5],"E"); strcpy(ALFA[6],"F"); strcpy(ALFA[7],"G"); strcpy(ALFA[8],"H");
    strcpy(ALFA[9],"I"); strcpy(ALFA[10],"J"); strcpy(ALFA[11],"K"); strcpy(ALFA[12],"L");
    strcpy(ALFA[13],"M"); strcpy(ALFA[14],"N"); strcpy(ALFA[15],"O"); strcpy(ALFA[16],"P");
    strcpy(ALFA[17],"Q"); strcpy(ALFA[18],"R"); strcpy(ALFA[19],"S"); strcpy(ALFA[20],"T");
    strcpy(ALFA[21],"U"); strcpy(ALFA[22],"V"); strcpy(ALFA[23],"W"); strcpy(ALFA[24],"X");
    strcpy(ALFA[25],"Y"); strcpy(ALFA[26],"Z");
  }

  // main section (as adapted from original BASIC)
  void main()  {
      unsigned char i, ic, l; char s[2]; int i4;
      RF=0;
      printf("\n                     SYMBOLIC PARSER for POLYNOMIALS\n");
      printf("\n                     Example:  (A+B)^2 ==> A^2+2AB+B^2\n\n");
      Fill_ALFA(); //ALFA[1]='A', ALFA[2]='B', etc. 
debut:if (RF==1) printf("\n\n");
      printf("  Inputs (0=SCREEN  1=algebra.dat file): "); scanf("%d", &RT);
      if (RT==0) {  //read input from screen
        printf("  Input string to evaluate:\n");
        printf("  "); scanf("%s", CHs);
        printf("\n");
      }
      else {  //read from input text file
	F1=fopen("algebra.dat","r");
        fscanf(F1,"%s", CHs);
        fclose(F1);
        printf("  String to evaluate:\n");
	//build up and print echo of input string to evaluate
        IC=strlen(CHs);
        if (IC<70) {
          printf("  ");for (i=0; i<IC; i++) printf("%c",CHs[i]); printf("%c\n",CHs[IC]);
        }
        else if (IC<140) {
          printf("  ");for (i=0; i<69; i++)  printf("%c",CHs[i]); printf("%c\n",CHs[69]);
          printf("  ");for (i=70; i<IC; i++) printf("%c",CHs[i]); printf("%c\n",CHs[IC]);
        }
        else if (IC<210) {
          printf("  ");for (i=0; i<69; i++)   printf("%c",CHs[i]); printf("%c\n",CHs[69]);
          printf("  ");for (i=70; i<139; i++) printf("%c",CHs[i]); printf("%c\n",CHs[139]);
          printf("  ");for (i=140; i<IC; i++) printf("%c",CHs[i]); printf("%c\n",CHs[IC]);
        }
        printf("\n");
      }
      //RE=0 --> no intermediate results
      printf("  Detail Analysis (0=NO  1=YES): "); scanf("%d", &RE);
      //RT=1 --> results are sent to a text output file "algebra.lst"
      printf("  Outputs (0=SCREEN  1=algebra.lst file): "); scanf("%d", &RT);
      if (RT==0) printf("\n");
      if (RT==1) {
	F1=fopen("algebra.lst","w");
        fprintf(F1,"\n  STRING TO EVALUATE:\n");
	//build up and print echo of input string to evaluate
        IC=strlen(CHs);
        if (IC<70) {
          fprintf(F1,"  ");for (i=0; i<IC; i++) fprintf(F1,"%c",CHs[i]); fprintf(F1,"%c\n",CHs[IC]);
        }
        else if (IC<140) {
          fprintf(F1,"  ");for (i=0; i<69; i++)  fprintf(F1,"%c",CHs[i]); fprintf(F1,"%c\n",CHs[69]);
          fprintf(F1,"  ");for (i=70; i<IC; i++) fprintf(F1,"%c",CHs[i]); fprintf(F1,"%c\n",CHs[IC]);
        }
        else if (IC<210) {
          fprintf(F1,"  ");for (i=0; i<69; i++)   fprintf(F1,"%c",CHs[i]); fprintf(F1,"%c\n",CHs[69]);
          fprintf(F1,"  ");for (i=71; i<139; i++) fprintf(F1,"%c",CHs[i]); fprintf(F1,"%c\n",CHs[139]);
          fprintf(F1,"  ");for (i=140; i<IC; i++) fprintf(F1,"%c",CHs[i]); fprintf(F1,"%c\n",CHs[IC]);
        }
        fprintf(F1,"\n");
      }

      // preliminary loop to add a * sign between
      // two parentheses )( or a digit and (.
	  // from release 2.1
	  l=strlen(CHs);
	  for (ic=0; ic<l-1; ic++) {
        Cs=CHs[ic];
        if ((Cs==')') && (CHs[ic+1]=='(')) {
          //Pascal: CHs=Copy(CHs,1,ic)+'*'+Copy(CHs,ic+1,strlen(CHs));
	  for (i=l+1; i>ic; i--) CHs[i]=CHs[i-1];
	  CHs[ic+1]='*';
        }
        if ((Cs>='0') && (Cs<='9') && (CHs[ic+1]=='(')) {
          //Pascal: Chs=Copy(Chs,1,ic)+'*'+Copy(CHs,ic+1,strlen(CHs));
          for (i=l+1; i>ic; i--) CHs[i]=CHs[i-1];
		  CHs[ic+1]='*';          
        }
      }

      IC=0;
      Cs=CHs[IC];
      NC=0;

      OPERA(CHs,&IC,&Cs,&NC,Ds);  //look for derivative £ signs

      // string analysis main loop
	  // -------------------------
e130: if (IC >= strlen(CHs)) goto e1000;
      if (Cs=='-')  ISI=-1; else ISI=1;
      if ((Cs=='+') || (Cs=='-')) {
        IC++;
        Cs=CHs[IC];
      }
      if (Cs=='(') goto e220;
      if (Cs!='£') goto e250;
      // opening parenthesis
      // -------------------
e220: // Pascal: PRs=CHs[IC-1]);
      // (C language does not allow to directly add
      //  a single character to a string)
      s[0]=CHs[IC-1]; s[1]='\0'; strcpy(PRs,s); 
      NC++;
      strcpy(CPs[NC],PRs);

      OPERA(CHs,&IC,&Cs,&NC,Ds);

      goto e130;
e250: if (Cs!=')') goto e460;
      // closing parenthesis
	  // -------------------
e290: if (Cs!=')') goto e300;
      NC++;
      strcpy(CPs[NC],")");
      IC++;
      Cs=CHs[IC];
      goto e290;
e300: SUs=CHs[IC+1];
      if ((SUs=='(') || (SUs=='£')) {
        IC++;
        Cs=SUs;                       
        goto e220;
      }
      if (Cs!='^') goto e370;
      // power
      NC++;
      strcpy(CPs[NC],"^");
      EX=1.0;
e340: if (Cs!='^') goto e350;
        NC--;
        if ((CHs[IC+2]>='0') && (CHs[IC+2]<='9'))  {
          //E1s=CHs[IC+1]+CHs[IC+2]    2 digits exponant
          E1s[0]=CHs[IC+1]; E1s[1]=CHs[IC+2]; E1s[2]='\0';
        } 
        else {
          //E1s=CHs[IC+1];             1 digit exponant
          E1s[0]=CHs[IC+1]; E1s[1]='\0';
        }
        //VAL(E1s,IE1,err);
	IE1=atoi(E1s);
        EX=EX*IE1;
        IC += strlen(E1s);
        Cs=CHs[IC];
      goto e340;

e350: NC++;
      CP[NC]=EX;
      SUs=CHs[IC+1];
      if ((SUs=='(') || (SUs=='£')) {  // £ for derivative
        IC++;
        Cs=SUs;
        goto e220;
      }

e370: if (Cs!='*') goto e130;
  // Begin mult. monom
  // -----------------
      NC++;
      strcpy(CPs[NC],"*");
      NC++;
      strcpy(CPs[NC],"");
      CP[NC]=1.0;
      goto e470;
  // End mult. monom
  // ---------------
e460: NC++;
      CP[NC]=(float) ISI;
      strcpy(CPs[NC],"");
  // Convert/Compact monom
e470: if ((CHs[IC]<'0') || (CHs[IC]>'9') || (IC>=strlen(CHs))) goto e510;
      strcpy(NBs,"");
  // Modif. for release 2.0  (2 digits exponant allowed)
      if (strlen(E1s)==2) IC1=IC-1; else IC1=IC;
  // Read an integer of unknown length
      while ((CHs[IC1]>='0') && (CHs[IC1]<='9') && (IC1<strlen(CHs))) {
        //Pascal: NBs=NBs+CHs[IC1];
        s[0]=CHs[IC1]; s[1]='\0'; strcat(NBs,s);
        IC1++;
      }
      //Pascal: VAL(NBs,NB,err);
      NB=atoi(NBs); //caution: no error handling here
					//the integers must be valid !!!
      IC += strlen(NBs);
      Cs=CHs[IC];
e490: if (Cs=='^') {
        IC++;
        //ENs=CHs[IC];
        s[0]=CHs[IC]; s[1]='\0'; strcat(ENs,s);
        //Pascal: VAL(ENs,IEN,err);
	IEN=atoi(ENs); //no error handling

        // NB (integer) power IEN
        INB=NB;
        NB=1;
        for (i=1; i<=IEN; i++)  NB *= INB;
        IC++;
        Cs=CHs[IC];
        goto e490;
      }
      CP[NC] *= NB;
e510: //current Cs is not an uppercase letter A..Z
      if ((Cs < 'A') || (Cs > 'Z') || (IC>=strlen(CHs))) goto e560;
      //current Cs is an uppercase letter A..Z
      //store rank of letter in NL
      for (i=1; i<27; i++)  if (ALFA[i][0]==Cs) NL=i;
      IE[NL]++; IC++;
      Cs=CHs[IC];
      IEX=1;

e540: if (Cs=='^') {
        IC++;
        if ((CHs[IC+1]>='0') && (CHs[IC+1]<='9')) {
          //Pascal: ECs=CHs[IC]+CHs[IC+1]
          ECs[0]=CHs[IC]; ECs[1]=CHs[IC+1]; ECs[2]='\0';
        }
        else
          ECs[0]=CHs[IC]; ECs[1]='\0';
        //Pascal: Val(ECs,IVL,err);
	IVL=atoi(ECs);
        IEX *= IVL;
        IC += strlen(ECs);
        Cs=CHs[IC];
        goto e540;
      }
      IE[NL] += (IEX-1);
e560: if (Cs=='*') {
        IC++;
        Cs=CHs[IC];
        if (Cs=='-') {
          CP[NC]=-CP[NC];
          IC++;
          Cs=CHs[IC];
        }
      }
      //Cs is still alphanumeric
      if ((Cs>='A') && (Cs<='Z')) goto e470;
      if ((Cs>='0') && (Cs<='9')) goto e470;
      for (i=1; i<27; i++) {
	if (IE[i]!=0) {
          //Pascal: CPs[NC]=CPs[NC]+ALFA[I];
          strcat(CPs[NC],ALFA[i]);
          if (IE[i]>1) {
            //Pascal: Str(IE[i],IEs);
            sprintf(IEs,"%d",IE[i]);
            //Pascal: CPs[NC]=CPs[NC]+'^'+IEs
	    strcat(CPs[NC],"^"); strcat(CPs[NC],IEs);
          }
          IE[i]=0;
        } // if IE(I)<>0
      } // i loop
      // end of closing parenthesis
	  // --------------------------.
      goto e130; //continue string analysis
      // end of string analysis main loop
	  // --------------------------------

      // simplify monoms if necessary
e1000:IDI=NC+2;
      if (!Verif(IDI)) {
        Message();   // too many terms
        goto fin;
      }
      IH=0;
      INO=0;
      LE=IDI / 2;
      for (NC=1; NC<=LE; NC++) {
        strcpy(ECs,CPs[NC]);
        strcpy(CPs[NC],CPs[IDI-NC+1]);
        strcpy(CPs[IDI-NC+1],ECs);
        EC=CP[NC];
        CP[NC]=CP[IDI-NC+1];
        CP[IDI-NC+1]=EC;
      } // end of NC loop
      for (i=3; i<=IDI; i++) {
	if (strcmp(CPs[i],")")==0) {
          IH++;
          IPI[IH]=i;
          INO++;
          IPO[INO]=i;
        }
        if (strcmp(CPs[i],"(")==0) {
          IHO=IPO[INO];
          CP[IHO]=i;
          CP[i]=(float) IHO;
          INO--;
        }
      } // end of i loop
	  //print intermediate result only if required (RE=1) 
      if (RE==1)  {
        if (RT==0)
          printf("  SIMPLIFICATION:\n");
        else
          fprintf(F1,"  SIMPLIFICATION:\n");
        IMPR();
      }

      // operations on monoms
      IHP=IH;  IHC=IH;   EC=0.0;
e1465:if (IHP==0) goto res;
e1470:IPF=IPI[IHC];
      JPO=(int) floor(CP[IPF]);
      IDM=JPO-1;
      IFM=IPF+1;
      IPP=IPF-1;
      IAP=JPO+1;
      if ((strcmp(CPs[IPP],"^")==0) && (strcmp(CPs[IPP-1],"(")==0)) {
        IHC--;
        goto e1470;
      }
      //call PUIS function if a "^" sign is detected
       if (strcmp(CPs[IPP],"^")==0)
       if (PUIS()!=0) {
	  if (!Ok) {
            Message();   //Power failed
            goto fin;
          }
          else RE=0;
        }

      if (strcmp(CPs[IAP],"^")==0) {
        EX=0.0;
        for (i4=IFM; i4<=IDM; i4++)  EX += CP[i4];
        IFM=IAP+2;
        IDM=(int) floor(CP[IAP+1])-1;
        IO1=(int) floor(CP[IAP+1]);
        IHP--;
        EC=1.0;
        if (PUIS()!=0) {
          if (!Ok) {
            Message();    //Power failed
            goto fin;
          }
          else RE=0;
        }
      } // end of if CPs[IAP]="^"
      Cs=CPs[IAP][1];
      if (Cs=='£') DERI();  //no action here
      strcpy(Ms,"");
      M=1.0; IMG=0; IMD=0;
      if ((strcmp(CPs[IAP],"*")==0) && (strcmp(CPs[IAP+1],")")!=0) && (strcmp(CPs[IAP+1],"$")!=0))
        IMG=1;
      Cs=CPs[IPP-1][0];
      if ((strcmp(CPs[IPP],"*")==0) && (strcmp(CPs[IPP-1],"(")!=0) && (Cs!='£'))
        IMD=1;
      if ((IMG!=0) || (IMD!=0)) goto e8000;
e1580:if (strcmp(CPs[IAP],"*")==0) {
	if (PARE()!=0) {
	  if (!Ok) {
            Message();     //PARE() failed
            goto fin;
          }
          else RE=0;
        }
        IHP--;
      }
      
      if (strcmp(CPs[IPP],"*")==0) {
        IHC--;
        goto e1470;
      }
      if (strcmp(CPs[IAP],"-")==0) SOUS(); else ADDI();
      IHP--; IHC--;
      goto e1465;

      // Multiplication of monoms
e8000:if (IMG==1) {
        strcat(Ms,CPs[IAP+1]);
        M=M*CP[IAP+1];
        strcpy(CPs[IAP],"$");
        CP[IAP]=0.0;
        strcpy(CPs[IAP+1],"$");
        CP[IAP+1]=0.0;
        IAP=IAP+2;
      } //if IMG=1

      if (IMD==1) {
        strcat(Ms,CPs[IPP-1]);
        M=M*CP[IPP-1];
        strcpy(CPs[IPP],"$");
        CP[IPP]=0.0;
        strcpy(CPs[IPP-1],"$");
        CP[IPP-1]=0.0;
        IPP=IPP-2;
      } //if IMD=1

      for (i=IDM; i>=IFM; i--) {
        if (CP[i]!=0.0) {
          CP[i]=CP[i]*M;
          //Ts=CPs[i]+Ms;
	  strcpy(Ts,CPs[i]); strcat(Ts,Ms);
          MONOME(Ts);
          strcpy(CPs[i],Ts);
        }
      } //i loop
      if (RE==1) {
        if (RT==0)
          printf("  MULTIPLICATION OF MONOMS:\n");
        else
          fprintf(F1,"  MULTIPLICATION OF MONOMS:\n");
        if (!IMPR()) goto fin;
      }
      goto e1580;

      // final simplifications and  result
res:  if (RT==0) printf("\n  FINAL SIMPLIFICATION AND RESULT:\n\n");
      else fprintf(F1,"\n  FINAL SIMPLIFICATION AND RESULT:\n\n");
      ISL=IDI-1;
	  //final simplification (eliminate zero valued monoms)
      for (i=3; i<=ISL; i++)
        if (CP[i]!=0.0)
          for (KS=i+1; KS<=IDI; KS++)
	    if (strcmp(CPs[KS],CPs[i])==0) {
              CP[i]=CP[i]+CP[KS];
              CP[KS]=0.0;
            }
      //final print (error message is string is too big!) 
      IMPR();
      if (RT==1) fclose(F1);
fin:  //want another algebra computing ?
      printf("\n  Evaluate another string (0=NO  1=YES): ");
      scanf("%d", &RF);
      printf("\n");
      if (RF==1) goto debut; //go to beginning...
  } // end of main program

//  Last update 02/02/2003 (11/07/1994 for pascal french edition)
//  J-P MOREAU, PARIS  (available fortran 77 and BASIC releases)

// End of file algebra.cpp