C
C
C     Filenames:
      CHARACTER*80 FILE2
C
C     Logical unit numbers:
      INTEGER LU2
      PARAMETER (LU2=2)
C
C     Array dimensions:
      INTEGER MLE,MH,MSPINE
      PARAMETER (MLE=400,MH=MLE**2+1,MSPINE=MLE*(MLE+1)/2)
C
      INTEGER ISP1(MSPINE),ISP2(MSPINE),ISP3(MSPINE)
      INTEGER IS1(MSPINE),IS2(MSPINE),IS3(MSPINE)
C
C.......................................................................
C
      NH=100
      FILE2='net.fs'
      WRITE(*,'(A)')
     *'+Enter maximum f.s. size (100), and f.s. filename (''net.fs''): '
      READ(*,*) NH,FILE2
      NH=NH*NH+1
      IF(NH.GT.MH) THEN
        STOP 'ERROR: TOO LARGE FORWARD STAR.'
      END IF
      OPEN(LU2,FILE=FILE2)
C
      DH=1./SQRT(12.*FLOAT(NH))
      DH=DH/2.
      NFS=1
      NSPINE=0
      DO 69 IH=1,NH
        DO 58 I1=INT(SQRT(FLOAT(IH)/3.)-DH+1.),INT(SQRT(FLOAT(IH))+DH)
          I=IH-I1*I1
          DO 57 I2=INT(SQRT(FLOAT(I)/2.)-DH+1.),
     *                                   MIN0(INT(SQRT(FLOAT(I))+DH),I1)
            I3I3=I-I2*I2
            IF(I3I3.EQ.0) THEN
                I3=0
                DO 10 J=2,I1
                  IF(MOD(I1,J).EQ.0.AND.MOD(I2,J).EQ.0) THEN
                    GO TO 56
                  END IF
   10           CONTINUE
C
C               New spine
                NSPINE=NSPINE+1
                IF(NSPINE.GT.MSPINE) THEN
                  STOP 'ERROR: TOO MANY SPINES.'
                END IF
                ISP1(NSPINE)=I3
                ISP2(NSPINE)=I2
                ISP3(NSPINE)=I1
C
   56           CONTINUE
            END IF
   57     CONTINUE
   58   CONTINUE
C
        IF(IH.EQ.NFS*NFS+1) THEN
          DO 66 ISPINE=NSPINE,3,-1
C           Looking for two closest spines
            II=ISP2(ISPINE)*ISP2(ISPINE)+ISP3(ISPINE)*ISP3(ISPINE)
            JL=2
            JLL=2
            AL=999999.
            JR=1
            JRR=1
            AR=999999.
            DO 64 JSPINE=1,NSPINE
              IF(ISP3(JSPINE).GT.0) THEN
                IF(JSPINE.NE.ISPINE) THEN
                  IJ=ISP2(ISPINE)*ISP3(JSPINE)-ISP3(ISPINE)*ISP2(JSPINE)
                  JJ=ISP2(JSPINE)*ISP2(JSPINE)+ISP3(JSPINE)*ISP3(JSPINE)
                  A=FLOAT(IJ*IJ)/FLOAT(II*JJ)
                  IF(IJ.LT.0) THEN
                    IF(A.LT.AL) THEN
                      JL=JSPINE
                      AL=A
                      JLR=ISP2(JL)*ISP3(JR)-ISP3(JL)*ISP2(JR)
                      JLL=ISP2(JL)*ISP2(JL)+ISP3(JL)*ISP3(JL)
                      R=FLOAT(JLR*JLR*IH)/FLOAT(JLL*JRR)
                      IF(R.LT.0.999999) THEN
                        ISP3(ISPINE)=-ISP3(ISPINE)
                        GO TO 65
                      END IF
                    END IF
                  ELSE
                    IF(A.LT.AR) THEN
                      JR=JSPINE
                      AR=A
                      JLR=ISP2(JL)*ISP3(JR)-ISP3(JL)*ISP2(JR)
                      JRR=ISP2(JR)*ISP2(JR)+ISP3(JR)*ISP3(JR)
                      R=FLOAT(JLR*JLR*IH)/FLOAT(JLL*JRR)
                      IF(R.LT.0.999999) THEN
                        ISP3(ISPINE)=-ISP3(ISPINE)
                        GO TO 65
                      END IF
                    END IF
                  END IF
                END IF
              END IF
   64       CONTINUE
   65       CONTINUE
   66     CONTINUE
          NS=0
          DO 67 ISPINE=1,NSPINE
            IF(ISP3(ISPINE).GT.0) THEN
              NS=NS+1
              IS1(NS)=ISP1(ISPINE)
              IS2(NS)=ISP2(ISPINE)
              IS3(NS)=ISP3(ISPINE)
            END IF
   67     CONTINUE
          WRITE(*,'(A,I7,2I8)') '+',NFS,NS
          WRITE(LU2,'(I3,I6,I7,E13.6)') NFS,NS
          WRITE(LU2,'(8(3I3,1X))') (IS1(I),IS2(I),IS3(I),I=1,NS)
          NFS=NFS+1
          DO 68 ISPINE=1,NSPINE
            ISP3(ISPINE)=IABS(ISP3(ISPINE))
   68     CONTINUE
        END IF
   69 CONTINUE
      WRITE(LU2,'(I3,I6,I7,E13.6)') 0,0
C
      STOP
      END
C
C=======================================================================
C
C