      PROGRAM POLPREC
      IMPLICIT NONE

C     POLPREC is based on POLPRE, which searches triangle 
C     presentations for (4,3)-polyhedra. 

C     POLPREC searches only those presentations, which have 
C     at least one cubic triangle (xj,xj,xj).

C	This is done by going through all possible ways to label the 
C 	rows of a given matrix (containing 
C	permutations of S_6) to obtain incidence tableaus. Each 
C	incidence tableau in turn gives a triangle presentation. 
C	POLPREC gives as a result one representative of each 
C	equivalence class of triangle presentations. 
     
      INTEGER I, J, K, L, P
      INTEGER IN, OUT, ERR
      PARAMETER(IN=5, OUT=6, ERR=0)
      INTEGER N, NMAX
      PARAMETER(NMAX=15)
      INTEGER M(NMAX,3), X(NMAX), NN(NMAX,3)

      INTEGER T(30,3)
      INTEGER KERTOJ, ROWTS, RTS, ART, PLA
      INTEGER UPO(4*NMAX,2), NUPO
      INTEGER POSS, WTS, XX
      REAL CUBES
      
      LOGICAL NEWVN, OK, SRCHVC, THISRR, SNDBCK, NXTPRM, NEWEC

C     Vectors, which tell about the status of the algorithm
      INTEGER C(30), U(30,3)
      
      EXTERNAL NXTPRM, ECHECK

      INTEGER MAXEC

C     Kokeillaan sitten vaihtoehtoista tapaa, luodaan iso matriisi
      INTEGER TRI(30*10000, 3), CUBEL(10000)


      DATA X/1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14 ,15/

C      DATA X/1,2,3,4,10,13,15,8,9,11,7,6,12,14,5/

C      DATA X/1, 2, 3, 4, 6, 15, 10, 12, 9, 14, 5, 13, 8, 11, 7/

C      DATA X/1, 4, 3, 6,7,5, 2, 8, 9,10,11, 12, 13, 14 ,15/

C      DATA X/1,2,10,3,9,15,14,4,13,6,7,8,11,5,12/

C     Case n=7
C      DATA M/2, 3, 2, 1, 4, 5, 3, 0, 0, 0, 0, 0, 0, 0, 0,
C     *     7, 6, 5, 6, 2, 4, 1, 0, 0, 0, 0, 0, 0, 0, 0,
C     *     6, 4, 3, 5, 1, 7, 7, 0, 0, 0, 0, 0, 0, 0, 0/
C      N = 7
      
C     Case n=15
      DATA M/1, 1, 1, 2, 2, 2, 3, 3, 3, 4, 4, 4, 5, 5, 5,
     *     10, 11, 12, 7, 8, 9, 6, 8, 9, 6, 7, 9, 6, 7, 8,
     *     15, 14, 13, 15, 14, 13, 15, 12, 11, 14, 12, 10, 13, 11, 10/
      N = 15

C     The number of the equivalence classes found so far
      MAXEC=0

      NEWVN = .TRUE.
      KERTOJ = 0
      XX = 0

C      OPEN(UNIT=8, FILE='results.txt')
      OPEN(UNIT=9, FILE='x_list.txt')
      OPEN(UNIT=7, FILE='T_list.txt')
      OPEN(UNIT=8, FILE='cubes.txt')

C testi:
CC      P=0

 10   IF (NEWVN) THEN
         KERTOJ = KERTOJ + 1

C testi:
CC         IF (P .GE. 100000000) THEN
CC            GOTO 9999
CC         ELSE
CC            P=P+1
CC         END IF

         OK = .TRUE.
 20      IF (OK) THEN
            OK = .FALSE.
            IF (.NOT. NXTPRM(X,N)) GOTO 9999
            IF (X(4) .NE. XX) THEN
               XX = X(4)
C               WRITE(UNIT=8,FMT=*) 'X(1:4) =', X(1), X(2), X(3), X(4)
               WRITE(OUT,*) 'X=', X(1),X(2),X(3),X(4)
            END IF
            GOTO 20
         END IF

CC         WRITE(*,*) X

C	When vectors beginning with 1 or 2 are gone through, 
C	jump to those beginning with 11.

CC	IF (X(1) .EQ. 3) THEN
CC	    X(1)=11
CC	    X(2)=1
CC	    X(3)=2
CC	    X(4)=3
CC	    X(5)=4
CC	    X(6)=5
CC	    X(7)=6
CC	    X(8)=7
CC	    X(9)=8
CC	    X(10)=9
CC	    X(11)=10
CC	    X(12)=12
CC	    X(13)=13
CC	    X(14)=14
CC	    X(15)=15
CC	END IF

C	When those beginning with 11 are tested, end the program.
CC	IF (X(1) .EQ. 12) GOTO 9999



C        Order the rows of matrix M the way x describes
         DO 36 I = 1, N
            C(I) = 0
            DO 35 J = 1, 3
               U(I,J) = 0 
               T(I,J) = 0
               NN(X(I),J) = M(I,J)
 35         CONTINUE
 36      CONTINUE

         DO 40 I = N+1, 30
            C(I) = 0
            DO 39 J = 1, 3
               U(I,J) = 0 
               T(I,J) = 0
 39         CONTINUE
 40      CONTINUE

C        Search for the triplets:
        
         SRCHVC = .TRUE.
         ROWTS = 1
C        "already tried place on last obtained row"
         ART = 0
         NUPO = 0
         CUBES = 0
 
 50      IF (SRCHVC) THEN
            IF (ROWTS.GE.1.AND.ROWTS.LE.(N+(FLOOR((CUBES+2)/3))*2))THEN


CC               WRITE(OUT, '(3I3)') ((T(I,J), J=1, 3), I=1,20)
CC               WRITE(*,*) ROWTS
CC               WRITE(*,*) ART
CC               WRITE(*,*) CUBES
CC               WRITE(OUT, '(3I3)') ((U(I,J), J=1, 3), I=1,15)

               RTS = ROWTS

C              choose first unused place from the matrix
               PLA = 1
 60            I = PLA / 3
               IF (MOD(PLA,3) .GT. 0) I = I + 1
               J = MOD(PLA-1,3) + 1
               IF (U(I,J) .EQ. 1 .AND. PLA .LT. 3*N) THEN
                  PLA = PLA + 1
                  GOTO 60
               END IF
               IF (PLA .EQ. (3*N+1)) THEN
                  NEWVN = .FALSE.
                  SRCHVC = .FALSE.
                  GOTO 10
               END IF

C Unused element (i,j) found.
C Continue to find a triangle with this element.
               
               T(RTS,1) = NN(I,J)
               U(I,J) = 1
               NUPO = NUPO + 1
               UPO(NUPO,1) = I
               UPO(NUPO,2) = J
               T(RTS,3) = I

               J = 1
 70            IF (C(RTS) .LE. 4) THEN
                  THISRR = .FALSE.
                  POSS = NN(T(RTS,1),J-1)
                  IF (J .GT. ART) THEN
C Search for cubic triangle
                     IF( J .EQ. 1) THEN
                        IF(T(RTS,1) .EQ. T(RTS,3)) THEN
                           CUBES=CUBES+1
CC                           WRITE(*,*) 'Cube'
                           T(RTS,2)=T(RTS,1)
                           U(T(RTS,1),1) = 1
                           NUPO = NUPO + 1
                           UPO(NUPO,1) = UPO(NUPO-1,1)
                           UPO(NUPO,2) = UPO(NUPO-1,2)
                           NUPO = NUPO + 1
                           UPO(NUPO,1) = UPO(NUPO-1,1)
                           UPO(NUPO,2) = UPO(NUPO-1,2)
                           C(RTS) = J
                        ELSE
                           C(RTS)=0 
                        END IF
                        
                     ELSE                     
C Search for non-cubic triangle
                        IF (NN(POSS,1) .EQ. T(RTS,3) 
     *                       .AND. U(T(RTS,1),J-1) .EQ. 0
     *                       .AND. U(POSS,1) .EQ. 0) THEN
                           T(RTS,2) = POSS
                           U(T(RTS,2),1) = 1
                           U(T(RTS,1),J-1) = 1
                           NUPO = NUPO + 1
                           UPO(NUPO,1) = T(RTS,1)
                           UPO(NUPO,2) = J-1
                           NUPO = NUPO + 1
                           UPO(NUPO,1) = T(RTS,2)
                           UPO(NUPO,2) = 1
                           C(RTS) = J
                        ELSE IF (NN(POSS,2) .EQ. T(RTS,3) 
     *                          .AND. U(T(RTS,1),J-1) .EQ. 0
     *                          .AND. U(POSS,2) .EQ. 0) THEN
                           T(RTS,2) = POSS
                           U(T(RTS,2),2) = 1
                           U(T(RTS,1),J-1) = 1                     
                           NUPO = NUPO + 1
                           UPO(NUPO,1) = T(RTS,1)
                           UPO(NUPO,2) = J-1
                           NUPO = NUPO + 1
                           UPO(NUPO,1) = T(RTS,2)
                           UPO(NUPO,2) = 2
                           C(RTS) = J
                        ELSE IF (NN(POSS,3) .EQ. T(RTS,3)
     *                          .AND. U(T(RTS,1),J-1) .EQ. 0
     *                          .AND. U(POSS,3) .EQ. 0) THEN
                           T(RTS,2) = POSS
                           U(T(RTS,2),3) = 1
                           U(T(RTS,1),J-1) = 1
                           NUPO = NUPO + 1
                           UPO(NUPO,1) = T(RTS,1)
                           UPO(NUPO,2) = J-1
                           NUPO = NUPO + 1
                           UPO(NUPO,1) = T(RTS,3)
                           UPO(NUPO,2) = 2
                           C(RTS) = J
                        ELSE
                           C(RTS) = 0
                        END IF
                     END IF
                  END IF
                  
                  IF (C(RTS) .NE. 0) THEN
                     THISRR = .TRUE.
                     ROWTS = ROWTS + 1
                     GOTO 200 
                  END IF

                  IF (J .EQ. 4 .AND. C(RTS) .EQ. 0) THEN
                     
C                    check if it is possible to return to an earlier row
                     SNDBCK = .FALSE.
                     WTS = RTS - 1

 100                 IF (WTS .GE. 1) THEN
                        IF (C(WTS) .LE. 2) THEN
C     If deleting a cube
                           IF(T(WTS,2) .EQ. T(WTS,3) .AND. 
     +                          T(WTS,1) .EQ. T(WTS,2)) THEN
                              CUBES=CUBES-1
CC                              WRITE(*,*) ' Deleting Cube'
                              IF(FLOOR(CUBES) .LT. 0) THEN
                                 WRITE(*,*) 'Negative cube'
                                 WRITE(OUT, '(3I3)') ((T(I,J), J=1, 3), 
     +                                I=1,30)
                                 WRITE(*,*) WTS
                                 WRITE(*,*) X
                              END IF
                           END IF
   
                           SNDBCK = .TRUE.
                           GOTO 110
                        ELSE
                           WTS = WTS - 1
                        END IF
                        GOTO 100
                     END IF
                     
C                    returning to an earlier row
 110                 IF (SNDBCK) THEN
                           
                        ROWTS = WTS
                        ART = C(WTS)
                        J = 1
                        DO 120 K = 1, WTS-1
                           U(K,1) = 0
                           U(K,2) = 0
                           U(K,3) = 0
 120                    CONTINUE
                        DO 130 K = WTS, (N+(FLOOR((CUBES+2)/3))*2)
                           C(K) = 0
                           U(K,1) = 0
                           U(K,2) = 0
                           U(K,3) = 0
 130                    CONTINUE
                        NUPO = 3 * (WTS-1)
                        IF (NUPO .GT. 1) THEN
                           DO 140 L = 1, NUPO
                              U(UPO(L,1), UPO(L,2)) = 1
 140                       CONTINUE
                        END IF
                        RTS = ROWTS
                        
C                    not able to return to an earlier row
                     ELSE
                        SRCHVC = .FALSE.
                     END IF
                     GOTO 200
                  ELSE
                     J = J + 1
                  END IF
                  GOTO 70
               END IF
            END IF

C           check whether a solution is found
 200        IF (ROWTS .EQ. (N+(FLOOR((CUBES+2)/3))*2+1)) THEN 
CC     +           .AND. MOD(FLOOR(CUBES),3) .EQ. 0) THEN

               SRCHVC = .FALSE.

               IF (C(ROWTS-1) .NE. 0) THEN 

C	Check if whether we have a new equivalence class or not
		  NEWEC=.TRUE.

                  CUBEL(MAXEC+1)=CUBES

                  IF (MAXEC .GT. 0) THEN
                     WRITE(OUT,*) 'checking for equivalence'
                     CALL ECHECK(N,T,TRI,CUBEL,MAXEC,NEWEC)
                  END IF

                IF (NEWEC .EQV. .TRUE.) THEN
CCCCCCC
C Found cubic case
                   IF (CUBES .GT. 0) THEN


                      WRITE(*,*) 'new equivalence class!'
                      MAXEC=MAXEC+1
C                      WRITE(UNIT=8,FMT=400) MAXEC, CUBES 
                      WRITE(UNIT=8,FMT=401) FLOOR(CUBES)
                      
                      WRITE(*,*) MAXEC
                      
                      WRITE(*,*) CUBES
                      
C     print the result
                      WRITE(OUT,'(A, 15I3)') 'x:', (X(I), I=1,N)
                      WRITE(OUT,'(A)') 'T:'
                      WRITE(OUT, '(3I3)') ((T(I,J), J=1,3), 
     +                     I=1,N+(FLOOR(CUBES/3))*2)
                      WRITE(UNIT=9, FMT=500) (X(I), I=1,N)
                      WRITE(UNIT=7,FMT=300) ((T(I,J),J=1,3),
     +                     I=1,N+(FLOOR(CUBES/3))*2)
                      
                      DO 270 J=1, (N+(FLOOR(CUBES/3))*2)
                         DO 260 I=1,3
                            TRI(30*(MAXEC-1)+J,I)=T(J,I)
 260                     CONTINUE
 270                  CONTINUE
                      DO 290 J=(N+(FLOOR(CUBES/3))*2+1),30
                         DO 280 I=1,3
                            TRI(30*(MAXEC-1)+J,I)=0
 280                     CONTINUE
 290                  CONTINUE
                      CUBEL(MAXEC)=CUBES
                   END IF
CCCCCCCCCCCCC
                   GOTO 10
                ELSE
                   WRITE(*,*) 'old equivalence class'
C                   WRITE(UNIT=8,FMT=*) MAXEC
                   GOTO 10
                END IF
             END IF
          END IF
          
          GOTO 50
       END IF
       
       GOTO 10
      END IF
      
 300  FORMAT (3I3)
 400  FORMAT (I4,F8.5)
 401  FORMAT (I4)
 500  FORMAT (15I3)

 9999 END











