!PolyGen.f95
!2021.04.06.0610cdt JMS- Traveler2/Athlon64/WinXPPro/APF9.0/Ogl1.2.1
!
!              Up to 8th order and seven variables this program
!                generates N-Variable polynomial coefficients
!                spanning a [-1.,1.]**#ofVariables hyperspace
!            and computes scaled exact integer matrix coefficients
!                             
!    The user specifies the number & order of each variable in `PolyGen.ini`

!Table of Contents:
!Module PolyGenDef
!Program Main
! Subroutine ReadPolyGenIni(iP)
! Subroutine Allocator(iP)
! Subroutine PowersMatrixGen(iP)
! Subroutine ComputeCoe(iP)
! Subroutine MatrixExportPrep( nTot,A,ADiv,AMin,AMax,AMaxScaled &
!                             ,iAout,iADigits,iADeciPl,iP)
! Subroutine ExportBinaryData(iP)
! Subroutine Deallocator(iP)
! Subroutine Invertr16(N,A,ValMin,iRank,DetN,Vused,iP)
! Subroutine Printr16(N,A,Noise,iRank,DetN,iRu,jCu,iP)
! Subroutine DisclaimerEtc(iP)
! Subroutine CoefficientCount(iP)
!---------------------------------------------------------------------------7 9

Module PolyGenDef
!use PolyGenDef                                                 !JMS 2021.03.31
! -- Polynomial generator data structures
! -- Specify the multi-variable polynomial:
  integer*4,public            ::nVT       !Variables- total number [1,2,...,10]
  integer*4,public            ::nOrd(10)  !Order of polynomial of each variable
!                                          1's are forced to be 2's.
!                                          for symmetric & on-center relevance.
  integer*4,public            ::nOrdMax   !Largest of the nOrd()'s
  integer*4,public            ::nCT       !Coefficients total number
  integer*4,public            ::MemoryEst !Real*16 memory estimate:
! -- Specify the polynomial's domain & output characteristics:
  character,public            ::cFileOut*79!Filename of the binary matrices
  character*1    ::AugHex(0:20)=(/"0","1","2","3","4","5","6","7","8","9","a" &
                                     ,"b","c","d","e","f","g","h","i","j","k"/)

  integer*4,public            ::iFileBytes !Binary file size in bytes

  integer*4,public            ::iOutFiles  &!Output type:
                                           =1!:"Comma separated values" (.csv)
!                                          =2 : Binary                  (.bin)
!                                          =3 : Both .csv & .bin
  integer*4,public            ::iPolyClass &!Coefficients type:
!                                          =1 :Exact
                                           =2!:Unscaled
  integer*4,public            ::iDomain    &
                                           =1!:[ -1.  , 1. ]
!                                          =2 :[-N/2. ,N/2.]
  integer*4,public            ::iReInSm    &!Real or integer or small arrays
                                           =1!:real
!                                          =2 :integer
!                                          =3 :most compact representation
  integer*4,public            ::iOutCoes   &!Binary data type to use
!                                          =0 :By default- none
!                                          =1 :i* 1-   8 bits ~  2 digits
!                                          =2 :i* 2-  16  "   ~  4   "
!                                          =3 :r* 4-  32  "   ~  6   "
!                                          =4 :i* 4-  32  "   ~  9   "
!                                          =5 :r* 8-  64  "   ~ 15   "
!                                          =6 :i* 8-  64  "   ~ 18   "
                                           =7!:r*16- 128  "   ~ 31   "
  integer*4,public            ::iDecPlaces!Number of decimal places in .csv
!                                          =[0,1,2,3,4,...,31,32]
! -- Refine the runtime estimate:
  real*16                     ::TEC       &!Time Estimate Coefficient
                                           = 208408._16
!                                           units: nCT's squared per minute 
  real*16                     ::TimeEstMin!= nCT*nCT/iToc ~Runtime in minutes
!                                             ...
  integer*4,public            ::iPolyAlloc =0!:unAllocated
!                                          =1 :  Allocated
! --Data samplepoint values:
  real*16  ,public,allocatable::X( :,:)   !Coefficient values in [-1,1]
!                         (0:nOrdMax,nVT)

!  --Sequence of matrix values for both rows and columns of both matrices:
  integer*4,public,allocatable::xRC(:,:)  !(_,Repeats & Cycles) counts
!                                (2,nVT)   = product(nOrd(i)+1),i=1,nVT)
!                                         This bypasses the need for recursion!
  integer*4,public,allocatable::nVnC(:,:) !Indexes & powers lookup table
!                                 (nVT,nCT)
! --Dataset powers matrix:
  real*16  ,public,allocatable:: Pow(:,:) !The matrix to be inverted
!                                 (nCT,nCT)
  real*16  ,public            :: PowDiv   !Divisor of Coe()- real
  real*16  ,public            :: PowMin   !MinAbs Pow() coefficient >1.e-20
  real*16  ,public            :: PowMax   !MaxAbs coefficient-
  real*16  ,public            :: PowMaxScaled !              - scaled
  integer*4,public            ::iPowOut   !Binary reporting Compression:
!                                          =0: By default- none
!                                          =1: i* 1-   8 bits ~  2 digits
!                                          =2: i* 2-  16  "   ~  4   "
!                                          =3: r* 4-  32  "   ~  6   "
!                                          =4: i* 4-  32  "   ~  9   "
!                                          =5: r* 8-  64  "   ~ 15   "
!                                          =6: i* 8-  64  "   ~ 18   "
!                                          =7: r*16- 128  "   ~ 31   "
  integer*4,public            ::iPowDigits!Number of digits for exact values
  integer*4,public            ::iPowDeciPl!Number of decimal places

! --Polynomial coefficients matrix:
  real*16  ,public,allocatable:: Coe(:,:) !The matrix to be inverted
!                                 (nCT,nCT)
  real*16  ,public            :: CoeMult  !Additional Coe() multiplier
  real*16  ,public            :: CoeDiv   !Divisor of Coe()- real
  real*16  ,public            :: CoeMin   !MinAbs Coe() coefficient >1.e-20
  real*16  ,public            :: CoeMax   !MaxAbs coefficient-
  real*16  ,public            :: CoeMaxScaled !              - scaled
  integer*4,public            ::iCoeOut   !Binary reporting Compression:
!                                          =0: By default- none
!                                          =1: i* 1-   8 bits ~  2 digits
!                                          =2: i* 2-  16  "   ~  4   "
!                                          =3: r* 4-  32  "   ~  6   "
!                                          =4: i* 4-  32  "   ~  9   "
!                                          =5: r* 8-  64  "   ~ 15   "
!                                          =6: i* 8-  64  "   ~ 18   "
!                                          =7: r*16- 128  "   ~ 31   "
  integer*4,public            ::iCoeDigits!Number of digits for exact values
  integer*4,public            ::iCoeDeciPl!Number of decimal places

  real*16  ,public,allocatable:: Vused(:) !Matrix inverter reporting.
!                                    (nCT)
  real*16  ,public            :: Frac     !The division fraction accumulator

  real*16  ,public            ::CoePowMax ! AbsMax value in [Coe*Pow-I]

  character,public            ::Chrono(0:9)*22!TimeStamps:
!                                     (0):                for differential use
!                                     (1):                Start
!                                     (2):Data entry    - done
!                                     (3):Powers matrix - done
!                                     (4):Powers matrix - inverted

!                                     (5):Coeffs matrix - done
!                                     (6):Binary file   - exported
!                                     (7):CSV file      - exported
!                                     (8):                Done
!                                     (9):Elapsed runtime

! --The external runtime environment variables:
  character::cIniFIle*79  ='PolyGen.ini'c
  character::cOutFile*79  ='PCoes.txt'c
  integer*4::uD  =14 !: Open(uD...) for reading & writing documentation
  integer*4::uP  =13 !: write(uP,...) Print to- file
  integer*4::uS  = 6 !: write(uS,...) Print to- screen
  integer*4::uT  = 0 ! [0,uP, or uS] temporarily requested outputs
!                     if(uT>5) write(uT,...)
  logical  ::LExists
  integer*4::iAlloc     !Allocation flag: used with allocatable arrays     flag
!Contains
End Module PolyGenDef
!---------------------------------------------------------------------------7 9
Program Main                                                   !2021.03.31.1710
!2021.03.31.1710cdt JMS- Opens an output file & calls PolyGen()
 use PolyGenDef,only: &
   cIniFile,cOutFile,iAlloc,LExists,uP,uS,Chrono,TimeEstMin,nCT,iOutFiles,TEC &
  ,Coe,CoeDiv,CoeMin,CoeMax,CoeMaxScaled,iCoeOut,iCoeDigits,iCoeDeciPl        &
  ,Pow,PowDiv,PowMin,PowMax,PowMaxScaled,iPowOut,iPowDigits,iPowDeciPl
!--End Globals
 implicit none
 !--Arguments
 !--Internals
  character::cL*80
  !--EndDefs---------------------------
   uP = 13  ! = 6 ...for screen only (e.g. for no writes to a hard drive.)
                                                  call DateAndTime22(Chrono(1))
   !--- Open output
   if(uP.gt.12) then
     open( uP,file=cOutFile,action='write',err=10)
       cL='Opening output file: '//cOutFile    ;call TimeStamp(cL,chrono(1),uP)
       goto 20
   endif !uP>12
10     uP=6
       cL='The computer monitor will be the output. 'c ;call WSF(cL,uP)
       cL=''c                                                 ;call WSF(cL,uP)
20 continue
   cL='This is:    "PolyGen.exe",         By: Jeff Setterholm'c
                                                               ;call WSF(cL,uP)
   cL='         Version 1. 2021.04.06    Lakeville, MN 55044 USA'c
                                                                call WSF(cL,uP)
   cL=''c                                                      ;call WSF(cL,uP)
   cL=' Output: a large binary file of very accurate matrix coefficients.'c
                                                                call WSF(cL,uP)
   cL='       : Up to 8th order the nearest integers are scaled exact values'c
                                                                call WSF(cL,uP)
   cL='           & can be written to a .csv file by `Polyize.exe`.'c
                                                                call WSF(cL,uP)
 cL='       : 2000+ coefficient polynomials up to 20th order can be computed.'c
                                                                call WSF(cL,uP)
   cL=''c                                                      ;call WSF(cL,uP)
   cL='"IntroToPolyGen.txt" - has some background information.'c
                                                                call WSF(cL,uP)
   !--- Checking for PolyGen.ini
   Inquire(file=cIniFIle,exist=LExists)
     if(LExists)then
       cL='       "PolyGen.ini" - has your program control inputs.'c
                                                                call WSF(cL,uP)
      else; cL='       "PolyGen.ini" - wasn`t found.'c         ;call WSF(cL,uP)
            cL='          The file is self-documented.'c       ;call WSF(cL,uP)
       pause"               Press return to exit PolyGen..."          ;goto 100
     endif !LExists
   !-- end Inquiry.

!  call CoefficientCount(uP)

   call ReadPolyGenIni(uP) !/////////////////////
   call DateAndTime22(Chrono(2)) !Data entry - done
   write(cL,"('Data entry - done')");call TimeStamp(cL,Chrono(2),uP)
   call ElapsedTime(Chrono(1),Chrono(2),Chrono(9),TimeEstMin,0)
   call TimeStamp( & !                                  |
   "                                       task`s time: "c       ,Chrono(9),uP)

   call Allocator(uP) !//////////////////////////

   call PowersMatrixGen(uP) !////////////////////
   call MatrixExportPrep( nCT,Pow,PowDiv,PowMin,PowMax,PowmaxScaled &
                            ,iPowout,iPowDigits,iPowDeciPl,uP) !////////

   call DateAndTime22(Chrono(3)) !Powers matrix - done
   write(cL,"('Powers matrix - done')")        ;call TimeStamp(cL,Chrono(3),uP)
   call ElapsedTime(Chrono(2),Chrono(3),Chrono(9),TimeEstMin,0)
   call TimeStamp( & !                                  |
   "                                       task`s time: "c       ,Chrono(9),uP)

   call ComputeCoe(uP) !/////////////////////////
   call DateAndTime22(Chrono(4)) !Powers matrix - inverted
   write(cL,"('Powers matrix - inverted')")    ;call TimeStamp(cL,Chrono(4),uP)
   call ElapsedTime(Chrono(3),Chrono(4),Chrono(9),TimeEstMin,0)
   call TimeStamp( & !                                  |
   "                                       task`s time: "c       ,Chrono(9),uP)


   call MatrixExportPrep( nCT,Coe,CoeDiv,CoeMin,CoeMax,CoeMaxScaled &
                            ,iCoeout,iCoeDigits,iCoeDeciPl,uP) !///////

   call DateAndTime22(Chrono(5)) !Coeffs matrix - done
   write(cL,"('Coeffs matrix - done')")        ;call TimeStamp(cL,Chrono(5),uP)
   call ElapsedTime(Chrono(4),Chrono(5),Chrono(9),TimeEstMin,0)
   call TimeStamp( & !                                  |
   "                                       task`s time: "c       ,Chrono(9),uP)

   Chrono(6)=Chrono(5)
   call ExportBinaryData(uP) !///////////////////
   call DateAndTime22(Chrono(6)) !Binary file   - exported
   write(cL,"('Binary file - exported')")    ;call TimeStamp(cL,Chrono(6),uP)
   call ElapsedTime(Chrono(5),Chrono(6),Chrono(9),TimeEstMin,0)
   call TimeStamp( & !                                  |
   "                                       task`s time: "c     ,Chrono(9),uP)

   Chrono(7)=Chrono(6) !No .csv export.

   call Deallocator(uP) !////////////////////////
   call DateAndTime22(Chrono(8)) !Done
     write(cL,"('Done.')")                     ;call TimeStamp(cL,Chrono(7),uP)
     call ElapsedTime(Chrono(7),Chrono(8),Chrono(9),TimeEstMin,0)
     call TimeStamp( & !                                  |
   "                                         task`s time: "c     ,Chrono(9),uP)
   write(cL,"('')")                                ;call WSF(cL,uP)
   call ElapsedTime(Chrono(1),Chrono(8),Chrono(9),TimeEstMin,uP)

   TEC = nCT*nCT/TimeEstMin
   write(cL,"('')")                                            ;call WSF(cL,uP)
   write(cL,"(f15.3,',:TEC- measured matrix coefficients per minute')") TEC
                                                                call WSF(cL,uP)
   write(cL,"('')")                                            ;call WSF(cL,uP)
   call DisclaimerEtc(uP) !//////////////////////
   cL=''c                                                      ;call WSF(cL,uP)
100 cL='Ending "PolyGen.exe"'c  
   call DateAndTime22(Chrono(9))                  ;call TimeStamp(cL,Chrono,uP)
   if(uP.gt.12) close(uP)
   pause 'Press enter to continue - thus closing this output screen & exiting.'
End Program Main
!---------------------------------------------------------------------------7 9

Subroutine ReadPolyGenIni(iP)
!2021.03.31.1710cdt JMS- Print out the run setup
 use PolyGenDef
!--End Globals
 implicit none
 !--Arguments
  integer*4::iP           !Write enable>5: write(iP,...)
!--Internals
  integer*4::N,nV,i,nCh
  character::cL*80,cString*80
  real*16  ::TECin
  !--EndDefs---------------------------
!  cL=''c                                                      ;call WSF(cL,iP)
!  cL='Subroutine ReadPolyGenIni():'c                          ;call WSF(cL,iP)

!  --Default values:
   nVT         = 4             !#1
   nOrd        = 0
   nOrd(1:2)   = (/2,2/)       !#2
   iPolyClass  = 1             !#3 exact
   TEC         = 208408._16    !#4 #of matrix coefficients computed per minute
   iOutFiles   = 2             !#5 .bin
   iDecPlaces  = 0             !#8 no trailing zeros after the decimal point
!  --Read the PolyGen.ini:
   cL='                       ... opening it now...'c          ;call WSF(cL,uP)
   Open(unit=uD,file=cIniFIle,action='read')
10   read(uD,"(a80)") cString
     read(cString(1:1),"(i1)") i
     select case(i)
      case(1); read(cString(3:),*) nVT                          !nVT
                                if(nVT< 1) nVT  =  1
                                if(nVT>10) nVT  = 10
      case(2);                             nOrd =  0            !nOrd(1:nVT)
               read(cString(3:),*) nOrd(1:nVT)

       do nV = 1,nVT ;if(mod(nOrd(nV),2)== 1) nOrd(nV) = nOrd(nV)- 1
                      if(    nOrd(nV)   == 0) nOrd(nV) =           2
                      if(    nOrd(nV)    >20) nOrd(nV) =          20
       enddo!nV
      case(3); read(cString(3:),*) iPolyClass                   !iPolyClass
                                if(iPolyClass.ne.2) iPolyClass = 1
      case(4); read(cString(3:),*,err=10) TECin                 !TEC
                                if(TECin<  1000._16) TEC = 200000._16
                                if(TECin>= 1000._16) Tec = TecIn
      case(0);  goto 20!                        ///Ends the .ini file read///
     end select!(i)
     goto 10
20 close(uD)

!  --Find the max order, count coefficients, & precompute the exact divisor:
   nOrdMax = maxval(nOrd); if(nOrdMax>8) iPolyClass = 2 !Unscaled
                 PowDiv = 1._16;            CoeDiv = 1._16;  nCt =  1
   do nV = 1,nVT;                                      nCT = nCT*(nOrd(nV)+1)
     if(iPolyClass==1) then
       select case(nOrd(nV))
        case(2); PowDiv = PowDiv*    1._16; CoeDiv = CoeDiv*  2._16
        case(4); PowDiv = PowDiv*   16._16; CoeDiv = CoeDiv*  6._16
        case(6); PowDiv = PowDiv*  729._16; CoeDiv = CoeDiv* 80._16
        case(8); PowDiv = PowDiv*65536._16; CoeDiv = CoeDiv*630._16
       end select!nOrd(nV)
     endif!(iPolyClass==1)
   enddo!nV

!  --Name the binary output file:
   cL=char(0)
   if(iPolyClass==1) then; write(cL(1:7),"('PCoesE',i1,\)")nVT; nCh = 7
                     else; write(cL(1:7),"('PCoesU',i1,\)")nVT; nCh = 7
   endif!(iPolyClass==1)
     do nV=nVT,1,-1;write(cL(nCh+1:nCh+1),"(a1,\)")AugHex(nOrd(nV))
                                                            nCh = nCh+1;enddo!nV
                    write(cL(nCh+1:nCh+4),"('.bin')")
               cFileOut = cL(    1:   79)

!  --Report the run configuration:
   cL='!--- Run configuration:'                                ;call WSF(cL,iP)
   cL='Values  #'c                                             ;call WSF(cL,iP)
   write(cL,"( i6,',:1 nVT             -independent variables')") nVT
                                                                call WSF(cL,iP)
   write(cL,"('      ,:2`order` =power of each variable- evens only')")
                                                                call WSF(cL,iP)
   do N = 1,nVT; write(cL,"(i6,','\)") nOrd(n); call WSF(cL,iP)        ;enddo!N
!    cL=''                                                     ;call WSF(cL,iP)
   write(cL,"( i6,',:3 iPolyClass       - 1:exact, 2:unscaled'   )") iPolyClass
                                                                call WSF(cL,iP)
   write(cL,"('      ,:4 TEC             - Time Estimate Coefficient')")
                                                                call WSF(cL,iP)
   write(cL,"('=',f13.1)")TEC                                  ;call WSF(cL,iP)
   cL='!--- Assigned internally:'c                             ;call WSF(cL,iP)
   write(cL,"( i6,',:5 iOutFiles Output- 1:.csv, 2:.bin, 3:both')")iOutFiles 
                                                                call WSF(cL,iP)
   write(cL,"( 6x,',:6 cFileOut        - `PolyGen`+hex+`.bin`'\)")
                                                                call WSF(cL,iP)
                                                                call WSF(cL,iP)
   write(cL,"('=',28x,79(a1,:))") (cFileOut(i:i),i=1,len_trim(cFileOut))
                                                                call WSF(cL,iP)
   write(cL, &
   "(22x,'order indices ^ are augmented hexadecimals:{ 0,15}->{0:f}')")
                                                                call WSF(cL,iP)
   write(cL,"(39x,'-&- 16->g 17->h 18->i 19->j   20 ->   k')")
                                                                call WSF(cL,iP)
   cL='!--- Overview:'c                                        ;call WSF(cL,iP)

   write(cL,"(i13,', :Highest `order`' )") nOrdMax             ;call WSF(cL,iP)
   write(cL,"(i13,', :# of variables')") nVT                   ;call WSF(cL,iP)
   do nV=1,nVT; write(cL,"(i13,', :Order of X(',i2,')')") nOrd(nV),nV
                                                      call WSF(cL,iP) ;enddo!nV
   write(cL,"(i13,', :# of  Poly  Coeff`s')") nCT              ;call WSF(cL,iP)
   write(cL,"(i13,', :# of matrix Coeff`s')") nCT*nCT          ;call WSF(cL,iP)
   write(cL,"(f17.3,',  estimated runtime (minutes)')") nCT*nCT/TEC
                                                                call WSF(cL,iP)
!  pause"in ReadPolyGenIni"
                                                                         return
End Subroutine ReadPolyGenIni

!---------------------------------------------------------------------------7 9

Subroutine Allocator(iP)
!2021.03.24.0815cdt JMS- Deallocates PolyGen arrays.
!                      - Traveler2/Athlon64/WinXPPro/APF9.0/Ogl1.2.1
 use PolyGenDef
!--End Globals
 implicit none
 !--Arguments
  integer*4::iP           !Write enable>5: write(iP,...)
 !--Internals
  character::cL*80
  !--EndDefs-----------------
!  --Deallocate any allocated arrays:
   if(iPolyAlloc.ne.0) call Deallocator(iP)
!  --Allocate the arrays:
   allocate(X(0:nOrdMax,nVT),stat=iAlloc);         if(iAlloc /= 0) then
       pause 'X(0:nOrdMax,nVT) alloc. error. `Enter` to exit.'    ;stop; endif
   X        = 0._16                                  ;iPolyAlloc = iPolyAlloc+1
   allocate(xRC(2,nVT),stat=iAlloc);               if(iAlloc /= 0) then
       pause 'XRC(2,nVT) alloc. error. `Enter` to exit.'          ;stop; endif
   XRC      = 0                                      ;iPolyAlloc = iPolyAlloc+1
   allocate(nVnC(nVT,nCT),stat=iAlloc);            if(iAlloc /= 0) then
       pause 'nVnC(nVT,nCT) alloc. error. `Enter` to exit.'       ;stop; endif
   nVnC      = 0                                     ;iPolyAlloc = iPolyAlloc+1
   allocate(Pow(nCT,nCT),stat=iAlloc);             if(iAlloc /= 0) then
       pause 'Pow(nCT,nCT) alloc. error. `Enter` to exit.'        ;stop; endif
   Pow       = 1._16                                 ;iPolyAlloc = iPolyAlloc+1
   allocate(Coe(nCT,nCT),stat=iAlloc);             if(iAlloc /= 0) then
       pause 'Coe(nCT,nCT) alloc. error. `Enter` to exit.'        ;stop; endif
   Coe      = 1._16                                  ;iPolyAlloc = iPolyAlloc+1
   allocate(Vused(nCT),stat=iAlloc);               if(iAlloc /= 0) then
       pause 'Vused(nCT) alloc. error. `Enter` to exit.'   ;stop; endif
   Vused    = 0._16                                  ;iPolyAlloc = iPolyAlloc+1
   write(cL,"('All arrays needed have been successfully allocated')")
                                                        call WSF(cL,iP)
   return
End Subroutine Allocator
!---------------------------------------------------------------------------7 9

Subroutine PowersMatrixGen(iP)
!2021.03.24.0815cdt JMS- Populates PolyGenDef
 use PolyGenDef
!--End Globals
 implicit none
 !--Arguments
  integer*4::iP           !Write enable>5: write(iP,...)
!--Internals
  integer*4::nV,nC,nCTlocal,N,nCyc,nRep,nRow,nCol
  real*16  ::n16,PowUse
  character::cL*80
  !--EndDefs---------------------------
     write(cL,"('')")                              ;call WSF(cL,iP)
     write(cL,"('Subroutine PowersMatrixGen():')")
                                                        call WSF(cL,iP)
!  --Compute the sample points of the variables:
   do nV  = 1,nVT;                   n16 = quad(nOrd(nV))
     do N = 0,nOrd(nV)
       if(iDomain==1) &
       X(N,nV) = ((2._16*quad(N) -n16)/n16  )  ![ -1. ,+1. ] ///////////
       if(iDomain==2) &
       X(N,nV) = ((2._16*quad(N) -n16)/2._16)  ![ -N/2,+N/2] ///////////
     enddo!N
   enddo!nV
                 nCTlocal = 1
   do nV = 1,nVT
     xRC(1,nV) = nCTlocal                          !# of repeats
     nCTlocal  = nCTlocal*(nOrd(nV)+1)
   enddo!nV
   do nV = 1,nVT
     xRC(2,nV) = nCTlocal/(xRC(1,nV)*(nOrd(nV)+1)) !# of cycles of the repeat
   enddo!nV

!  --Generate the indices & powers lookup table:
   do nV               = 1,nVT       ; nC = 0
     do nCyc           = 1,xRC(2,nV)              !Number of repeats
       do N            = 0,nOrd(nV)
         do nRep       = 1,xRC(1,nV) ; nC = nC +1 !Number of repetitions
           nVnC(nV,nC) = N                        !Indices = Powers
         enddo!nRep
       enddo!N
     enddo!nCyc
   enddo!nV

!  --Compute the Powers Matrix:
                        call Beamer( 0  ,nCT,6)
   do nRow   = 1,nCT   ;call Beamer(nRow,nCT,6)
     do nCol = 1,nCT
       do nV=1,nVT
         PowUse=X(nVnC(nV,nRow),nV)**nVnC(nV,nCol)
         if(nVnC(nV,nCol)==0) PowUse = 1._16
         Pow(nrow,nCol)=Pow(nrow,nCol)*PowUse
       enddo!nV
     enddo!nCol
   enddo!nRow
                                                                         return
End Subroutine PowersMatrixGen
!---------------------------------------------------------------------------7 9
Subroutine ComputeCoe(iP)
!2021.03.24.0815cdt JMS- Populates PolyGenDef
 use PolyGenDef
!--End Globals
 implicit none
 !--Arguments
  integer*4::iP           !Write enable>5: write(iP,...)
!--Internals
  integer*4::nC,nRow,nCol
! Inverter arguments:
  real*16  :: ValMin   = 0.d0
  integer*4:: iRank    = 0
  real*16  :: DetN    !Abs(determinant product) for each iteration
  real*16  :: CoePow
  character::cL*80

  !--EndDefs---------------------------
!  --Computing the coefficient matrix by inversion of Pow():
   Coe=Pow

   write(cL,"('')")                              ;call WSF(cL,iP)
   write(cL,"('Coe[,] = inverse(Pow[,]), inversion progress bar...')")
                                                        call WSF(cL,iP)

   call Invertr16(nCT,Coe,ValMin,iRank,DetN,Vused,0)  !iP!/////////////////////
   write(cL,"(e54.30,', :Determinant' )") DetN     ;call WSF(cL,iP)
   write(cL,"(f50.30,','              )") DetN     ;call WSF(cL,iP)

   if(iRank.ne.nCT) &
   write(cL,"('Pseudoinverse result. Missed ',i6,' dimensions')") nCT-iRank
   if(iRank.ne.nCT)                                     call WSF(cL,iP)

   call DateAndTime22(Chrono(9)) !Inversion - done
   write(cL,"('')")                            ;call TimeStamp(cL,Chrono(9),uP)

!  --Confirm inversion accuracy:
   write(cL,"('Evaluating the inversion accuracy...')")
                                                        call WSF(cL,iP)
   CoePowMax=0._16                                     ;call Beamer(0   ,nCT,6)
   do nRow=1,nCT                                       ;call Beamer(nRow,nCT,6)
     do nCol=1,nCT; CoePow=0._16
       do nC=1,nCT; CoePow=CoePow+Coe(nRow,nC)*Pow(nC,nCol)         ;enddo!nC
     if(nRow==nCol) CoePow=CoePow-1._16
     if(abs(CoePow)>CoePowMax) CoePowMax=abs(CoePow)          ;enddo;enddo!nRow
   write(cL,"(e12.3,', :AbsMax r*16 inversion error in [Coe*Pow-I]')") &
                                CoePowMax              ;call WSF(cL,iP)
                                                                         return
End Subroutine ComputeCoe
!---------------------------------------------------------------------------7 9

Subroutine MatrixExportPrep( nTot,A,ADiv,AMin,AMax,AMaxScaled &
                            ,iAout,iADigits,iADeciPl,iP)
!2021.04.01.0835cdt JMS- Preparing the Coe() matrix for export
 use PolyGenDef     , only:iPolyClass,iDecPlaces
!--End Globals
 implicit none
 !--Arguments
  integer*4::nTot
  real*16  ::A(nTot,nTot) !Pow(,) or Coe(,)
  real*16  ::ADiv
  real*16  ::AMin
  real*16  ::AMax         !MaxAbs coefficient-
  real*16  ::AMaxScaled   !              - scaled
  integer*4::iAOut        !Binary reporting Compression: Not yet implemented.
!                         =0: By default- none
!                         =1: i* 1-   8 bits ~  2 digits
!                         =2: i* 2-  16  "   ~  4   "
!                         =3: r* 4-  32  "   ~  6   "
!                         =4: i* 4-  32  "   ~  9   "
!                         =5: r* 8-  64  "   ~ 15   "
!                         =6: i* 8-  64  "   ~ 18   "
!                         =7: r*16- 128  "   ~ 31   "
  integer*4::iADigits     !Number of digits for exact values
  integer*4::iADeciPl     !Number of decimal places
  integer*4::iP           !Write enable>5: write(iP,...)
!--Internals
  integer*4::nRow,nCol
  real*16  ::AUse
  character::cL*80
  !--EndDefs---------------------------
   write(cL,"('!--- MatrixExportPrep():')");       ;call WSF(cL,iP)

   write(cL,"(e54.30,', Division factor')") ADiv   ;call WSF(cL,iP)
   write(cL,"(f50.30,','                )") ADiv   ;call WSF(cL,iP)
   write(cL,"('')")                                ;call WSF(cL,iP)

   AMax = maxval(abs(A))
   write(cL,"(e54.30,', matrix MaxAbs')") AMax     ;call WSF(cL,iP)
   write(cL,"(f50.30,','              )") AMax     ;call WSF(cL,iP)
                     AMin = 1._16
   do nRow   = 1,nTot
     do nCol = 1,nTot
       AUse = abs(A(nRow,nCol))
       if(AUse<1.e-20_16) cycle
       if(AUse<AMin) AMin=AUse
     enddo!nCol
   enddo!nRow
   write(cL,"('')")                                ;call WSF(cL,iP)
   write(cL,"(e54.30,', matrix MinAbs')") AMin     ;call WSF(cL,iP)
   write(cL,"(f50.30,','              )") AMin     ;call WSF(cL,iP)
   write(cL,"('')")                                ;call WSF(cL,iP)

!  --Exporting the coefficient matrix:
   iADeciPl=0
   if(iPolyClass==2) then
     write(cL,"('Unscaled matrix coefficients are needed.')")
                                                        call WSF(cL,iP)
     if(iDecPlaces>0) then; iADeciPl = iDecPlaces
      else
       write(cL,"('  0 decimal place results were also requested.')")
                                                        call WSF(cL,iP)
       write(cL,"('  Your .ini file input #8 can override this guess:')")
                                                        call WSF(cL,iP)

       iADeciPl = int(log10(1._16/AMin))+3
       if(iADeciPl>24) iADeciPl = 24
       write(cL,"( i19,', :decimal places used')") iADeciPl
                                                        call WSF(cL,iP)
     endif!(iDecPlaces>0)
   endif!(iPolyClass==2)

   AMaxScaled = ADiv*AMax
   write(cL,"(f36.3,', :Matrix - biggest integer- decimal')") AMaxScaled
                                                        call WSF(cL,iP)
   if(     AMaxScaled>100000000000000._16) then; iAOut = 16
    elseif(AMaxScaled>        1000000._16) then; iAOut =  8
    else                                       ; iAOut =  4
   endif!AMaxScaled>1.e15_16

   write(cL,"(f36.3,,', :storage space needed- i*1,i*2,i*4,i*8')") &
               AMaxScaled                              ;call WSF(cL,iP)
   write(cL,"('21098765432109876543210987654321, :digits index')")
                                                        call WSF(cL,iP)
   write(cL,"('  3         2         1         ')");call WSF(cL,iP)
   iADigits = int(log10(AMaxScaled))+3
   write(cL,"( i19,', : digits recommended w/`.`')") iADigits  ;call WSF(cL,iP)

!  if(iPolyClass==1) then  !Multiply A*ADiv during binary write
!   !Exactly rescale the Coe Matrix:
!    write(cL,"(8x,f50.30,',Matrix division factor'/)") ADiv   ;call WSF(cL,iP)
!    do nRow   = 1,nTot
!      do nCol = 1,nTot
!        AUse = A(nRow,nCol)*ADiv
!       !write(cL,"(2i4,f50.30)") nRow,nCol,AUse              ;call WSF(cL,iP)
!        if(AUse<1.e-5_16) then; A(nRow,nCol) = 0._16
!                          else; A(nRow,nCol) = floor(AUse+.5_16)
!        endif!(AUse<1.e-5_16)
!       !write(cL,"(8x,f50.30)")             AUse              ;call WSF(cL,iP)
!       !cL=""c                                                ;call WSF(cL,iP)
!      enddo!nCol
!    enddo!nRow
!  endif!iPolyClass==1)
                                                                         return
End Subroutine MatrixExportPrep
!---------------------------------------------------------------------------7 9

Subroutine ExportBinaryData(iP)
!2021.04.01.0725cdt JMS- Exports the Coe() matrix & supporting variables
 use PolyGenDef
!--End Globals
 implicit none
 !--Arguments
  integer*4::iP           !Write enable>5: write(iP,...)
!--Internals
  integer*4::nB
  character::cL*80

  !--EndDefs---------------------------
   write(cL,"('')")                              ;call WSF(cL,iP)
   write(cL,"('Subroutine ExportBinaryData():')")

 !  --Export the Binary matrices
   write(cL,*) cFileOut                                ;call WSF(cL,iP)
                                                             write(cL, &
   "('Binary file data write/read sequence- integer*4`s & real*16`s:')")
                                                        call WSF(cL,iP)
   nB = sizeof(nVT);     iFileBytes =            nB         ;write(  cL,  &
   "(i12,', i* 4 nVT       -Number of independent variables  fixed size' )") nB
                                                        call WSF(cL,iP)
   nB = sizeof(nOrdMax); iFileBytes = iFileBytes+nB         ;write(  cL,  &
   "(i12,', i* 4 nOrdMax   -Highest order variable             `    `'   )") nB
                                                        call WSF(cL,iP)
   nB = sizeof(nCT);     iFileBytes = iFileBytes+nB         ;write(  cL,  &
   "(i12,', i* 4 nCT       -Number of polynomial coefficients  `    `'   )") nB
                                                        call WSF(cL,iP)
   nB = sizeof(nOrd);     iFileBytes = iFileBytes+nB        ;write(  cL,  &
   "(i12,', i* 4 nOrd(20)  -Order of each variable             `    `'   )") nB
                                                        call WSF(cL,iP)
   nB = sizeof(iPolyClass);iFileBytes = iFileBytes+nB        ;write(  cL,  &
   "(i12,', i* 4 iPolyClass-=1:exact, =2:unscaled              `    `'   )") nB
                                                        call WSF(cL,iP)
   nB = sizeof(iCoeDigits);iFileBytes = iFileBytes+nB       ;write(  cL,  &
   "(i12,', i* 4 iCoeDigits- +***. <-count                     `    `'   )") nB
                                                        call WSF(cL,iP)
   nB = sizeof(iCoeDeciPl);iFileBytes = iFileBytes+nB       ;write(  cL,  &
   "(i12,', i* 4 iCoeDeciPl-number of digits after .           `    `'   )") nB
                                                        call WSF(cL,iP)
   nB = sizeof(CoeDiv);   iFileBytes = iFileBytes+nB        ;write(  cL,  &
   "(i12,', r*16 CoeDiv           -Coe[,] dividing factor      `    `'   )") nB
                                                        call WSF(cL,iP)
   nB = sizeof(CoeMaxScaled); iFileBytes = iFileBytes+nB    ;write(  cL,  &
   "(i12,', r*16 CoeMaxScaled     -AbsMax Coe[,]               `    `'   )") nB
                                                        call WSF(cL,iP)
   nB = sizeof(iPowDigits);iFileBytes = iFileBytes+nB       ;write(  cL,  &
   "(i12,', i* 4 iPowDigits- +***. <-count                     `    `'   )") nB
                                                        call WSF(cL,iP)
   nB = sizeof(iPowDeciPl);iFileBytes = iFileBytes+nB       ;write(  cL,  &
   "(i12,', i* 4 iPowDeciPl-number of digits after .           `    `'   )") nB
                                                        call WSF(cL,iP)
   nB = sizeof(PowDiv);   iFileBytes = iFileBytes+nB        ;write(  cL,  &
   "(i12,', r*16 PowDiv           -Pow[,] dividing factor      `    `'   )") nB
                                                        call WSF(cL,iP)
   nB = sizeof(PowMaxScaled); iFileBytes = iFileBytes+nB    ;write(  cL,  &
   "(i12,', r*16 PowMaxScaled     -AbsMax Pow[,]               `    `'   )") nB
                                                        call WSF(cL,iP)
   nB = sizeof(X);        iFileBytes = iFileBytes+nB        ;write(  cL,  &
   "(i12,', r*16 X[0:nOrdMax,nVT] -Variable values to use    allocated'  )") nB
                                                        call WSF(cL,iP)
   nB = sizeof(nVnC);     iFileBytes = iFileBytes+nB        ;write(  cL,  &
   "(i12,', i* 4 nVnC[nVT   ,nCT] -Values combinations sequence  `'      )") nB
                                                        call WSF(cL,iP)
   nB = sizeof(Coe);      iFileBytes = iFileBytes+nB        ;write(  cL,  &
   "(i12,', r*16 Coe[ nCT   ,nCT] -The coefficients solver       `'      )") nB
                                                        call WSF(cL,iP)
   nB = sizeof(Pow);      iFileBytes = iFileBytes+nB        ;write(  cL,  &
   "(i12,', r*16 Pow[ nCT   ,nCT]      Coe[,] = inverse(Pow[,])  `'      )") nB
                                                        call WSF(cL,iP)
   write(  cL,"(i12,', Total Bytes')") iFileBytes;  call WSF(cL,iP)
  write(cL,"(f15.3,12x,':write-time-minutes ballpark')")iFileBytes/1.e10_16
                                                        call WSF(cL,iP)
   open(unit=14,file=cFileOut,form='binary',action='write')
     write(14)  nVT, nOrdMax, nCT,nOrd, iPolyClass            &
               , iCoeDigits, iCoeDeciPl,CoeDiv, CoeMaxScaled  &
               , iPowDigits, iPowDeciPl,PowDiv, PowMaxScaled  &
               , X, nVnC, (Coe*CoeDiv), (Pow*PowDiv)      !<- Verify this!
   close(14)
                                                                         return
End Subroutine ExportBinaryData
!---------------------------------------------------------------------------7 9

Subroutine Invertr16(N,A,ValMin,iRank,DetN,Vused,iP)           !2021.03.22.0535
!2021.03.22.0535cdt JMS- Vused is an array returning the division sequence.
!2010.01.06.1205cst JMS- Traveler2/Athlon64/WinXPPro-32/APF9.0-32
!2006.05.25.0715cdt JMS- An instantiation of the "Setterholm Matrix Inverter"
!                          Expects,uses, & returns real*16's (~QUAD precision).
!---                                                                 No globals
 implicit none                                                       !arguments
  integer*4:: N
  real*16  :: A(N,N)
  real*16  :: ValMin
  integer*4:: iRank
  real*16  :: DetN    !Abs(determinant product) for each iteration
! integer*4:: iUsed(N)
  real*16  :: Vused(N)                                              !2021.03.14
  integer*4:: iP
  !---                                                               !internals
  integer*4,allocatable::iRu(:)
  integer*4,allocatable::jCu(:)
  real*16  ,allocatable::Temp(:)
  integer*4::i,iu,iu2,j,ju,ju2,L
  real*16  ::Amult,ValMax
  integer*8::iBinSum
  integer*4::iAlloc
   !-----                                                            !end defs
   if(iP.gt.5) write(iP,"('Overwriting Matrix Inversion:')")
   DetN=0._16 ;Vused(1:N)=0._16 ;iRank=0
   !---
   allocate(iRu(N),stat=iAlloc)
   if(iAlloc.ne.0) stop 'Invert: iRu(N) allocation error. Halt.'
   iRu=0
   allocate(jCu(N),stat=iAlloc)
   if(iAlloc.ne.0) stop 'Invert: jCu(N) allocation error. Halt.'
   jCu=0
   allocate(Temp(N),stat=iAlloc)
   if(iAlloc.ne.0) stop 'Invert: Temp(N) allocation error. Halt.'
   Temp=0._16
   !---   
     Temp=0._16; iRu   =0 ;jCu   =0 ;ValMin=0._16
   do i=1,N    ;iRu(i)=i ;jCu(i)=i                ;end do
           if(iP.gt.5) call Printr16(N,A,ValMin,iRank,DetN,iRu,jCu,iP) !*******
!  Inversion Loop: ----------------------------------------
                                                call Beamer(0,N,6)
   do L=1,N ;iU=0     ;jU=0     ;ValMax =0._16 ;call Beamer(L,N,6)
     do i=1,N  ;if(iRu(i).lt.0) cycle
       do j=1,N ;if(jCu(j).lt.0) cycle
         if(ValMax.ge.abs(A(i,j))) cycle
           ValMax=abs(A(i,j)) ;iu=i ;ju=j
       enddo!j
     enddo!i
     if(L==1) ValMin=ValMax*1.e-20_16!Establishes the "noise floor".
     if(L==1) DetN=ValMax
     if(iU.eq.0) Exit
     if(ValMax.lt.ValMin) Exit
     iRank=iRank+1
     if(L>1) DetN=DetN*ValMax
     ValMax=A(iu,ju) ;iu2=jCu(ju) ;ju2=iRu(iu)
    !if(iP.gt.5) write(iP,"(3i3,f15.6)")  L,iu,ju,ValMax
     if(iP.gt.5) write(iP,"(1x,2i4,f15.9,'   :iu ju pivot, & value')") &
                                              iu,ju,        ValMax
     Vused(L)=ValMax                                                !2021.03.14

     Temp(1:N)=A(1:N,ju)                            ;i=iRu(iu)
               A(1:N,ju)=A(1:N,ju2)                 ;  iRu(iu)=iRu(iu2)
                         A(1:N,ju2)=Temp(1:N)       ;          iRu(iu2)=-abs(i)
     Temp(1:N)=A(iu,1:N)                            ;j=jCu(ju)
               A(iu,1:N)=A(iu2,1:N)                 ;  jCu(ju)=jCu(ju2)
                         A(iu2,1:N)=Temp(1:N)/ValMax ;         jCu(ju2)=-abs(j)
     do i=1,N ;if(i.eq.iu2) cycle ;Amult=A(i,ju2)
       do j=1,N ;A(i,j)=A(i,j)-A(iu2,j)*Amult ;enddo!j
       A(i,ju2)=-Amult/ValMax
     enddo; A(iu2,ju2)=1._16/ValMax
           if(iP.gt.5) call Printr16(N,A,ValMin,iRank,DetN,iRu,jCu,iP) !*******
   enddo!L  Inversion loop- done. -------------------------

   iBinSum=0
   do i=1,n
     if((iRu(i).lt.0).and.(i.lt.63)) iBinsum=iBinSum+2**(i-1)
!!!! if(iRu(i).lt.0) iUsed(i)=1                                     !2021.03.14
!    Zero linearly-dependent rows and columns, if any:
     if(iRu(i).gt.0) A(i,1:N)=0._16
     if(jCu(i).gt.0) A(1:N,i)=0._16
   enddo!i
   if((iP >5).and.(iu == 0)) then
     write(iP,"('Linear dependency dealt with. Partial inverse results:')")
     write(iP,"(b63.2,' -Kp unused flag')") iBinSum !For the first 62 Kp's
                       call Printr16(N,A,ValMin,iRank,DetN,iRu,jCu,iP) !*******
   endif!iP>5 & iu==0
   if(iP.gt.5) write(iP,"('Overwriting Inverter- done.')")
   deallocate(iRu)
   deallocate(jCu)
   deallocate(Temp)
   return
End Subroutine Invertr16
!---------------------------------------------------------------------------7-9

Subroutine Printr16(N,A,Noise,iRank,DetN,iRu,jCu,iP)
!2018.09.04.1340cdt JMS- Prints Invert's progress.
!                          Expects,uses, & returns real*16's (~QUAD precision).
!                      - Traveler2/Athlon64/WinXPPro-32/APF9.0-32
 !---                                                                No globals
 implicit none                                                       !arguments
  integer*4:: N
  real*16  :: A(N,N)
  real*16  :: Noise
  integer*4:: iRank
  real*16  :: DetN
  integer*4:: iRu(N)
  integer*4:: jCu(N)
  integer*4:: iP
  !---                                                               !internals
  integer*4::i,j
   !-----                                                            !end defs
                                                              if(iP < 6) return
   write(iP,"('Inverter: @ Rank =',i5,8x,'abs(Det)=',e39.30/2x\)") iRank,DetN
   do j=1,N;   write(iP,"(i18,\)") jCu(j) ; enddo; write(iP,"('')")
   do i=1,N;   write(iP,"(i5,' ' ,\)") iRu(i)
     do j=1,N; write(iP,"(sp,f18.9\)") A(i,j); enddo; write(iP,"('')")
   enddo!i
   write(iP,"('')");                                                     return
End Subroutine Printr16
!---------------------------------------------------------------------------7-9

Subroutine Deallocator(iP)
!2021.03.24.0815cdt JMS- Deallocates PolyGen arrays.
!                      - Traveler2/Athlon64/WinXPPro/APF9.0/Ogl1.2.1
 use PolyGenDef
!--End Globals
 implicit none
 !--Arguments
  integer*4::iP           !Write enable>5: write(iP,...)
 !--Internals
  !--EndDefs-----------------
   if(iPolyAlloc>0) then
     deallocate(X)
     deallocate(xRC)
     deallocate(nVnC)
     deallocate(Pow)
     deallocate(Coe)
     deallocate(Vused)
     iPolyAlloc = 0
     if(iP>5) write(iP,"(1x,'Arrays deallocated')")
   endif!iPolyAlloc>0
   return
 End Subroutine Deallocator
!---------------------------------------------------------------------------7 9

Subroutine DisclaimerEtc(iP)
!2021.04.06.0610cdt JMS- Disclaimer, etc.
!Elapsed is updated first,
!       so using the same argument twice results in returning the DaTime value.
!--Arguments
 implicit none
  integer*4::iP           !Write enable>5: write(iP,...)
 !--Internals
  character::cL*80
  !--EndDefs-----------------------------------------------
 cL=""                                                         ;call WSF(cL,iP)
 cL="   Link: http://ftp.setterholm.com/ExactInversePolynomials/PolyGen.zip"c
                                                                call WSF(cL,iP)
 cL="Analyst: Jeff Setterholm  Lakeville MN 55044 USA    Tuesday 2020.04.06"c
                                                                call WSF(cL,iP)
 cL=""                                                         ;call WSF(cL,iP)
 cL="  ********, PolyGen.exe is post-copyright i.e.:FREE ,********"c
                                                                call WSF(cL,iP)
 cL=""                                                         ;call WSF(cL,iP)
 cL="  My legal disclaimer:"c                                  ;call WSF(cL,iP)
 cL="  ***********************************************************"c
                                                                call WSF(cL,iP)
 cL="  ********,  Individual cognition is always flawed  ,********"c
                                                                call WSF(cL,iP)
 cL="  ********,          including yours & mine.        ,********"c
                                                                call WSF(cL,iP)
 cL="  ********,                 - So: -                 ,********"c
                                                                call WSF(cL,iP)
 cL="  ********,    Use these results at your own risk.  ,********"c
                                                                call WSF(cL,iP)
 cL="  ***********************************************************"c
                                                                call WSF(cL,iP)
 cL="           Mitigate malfunctions. Nurture synergies."c
                                                                call WSF(cL,iP)
                                                                         return
End Subroutine DisclaimerEtc
!---------------------------------------------------------------------------7 9
Subroutine CoefficientCount(iP)
!2021.03.27.1035cdt JMS- # of coeff`s of even-order multi-variable Polynomials
!                                        up to 20th order.
 implicit none
 !--Arguments
  integer*4 ::iP           !Write enable>5: write(iP,...)
 !--Internals
  integer*4 ::nMax=7,j,nOrder(10)=1,nCh,iUsed,nLine
  integer*4 ::i1,i2,i3,i4,i5,i6,i7
  integer*4 ::j1,j2,j3,j4,j5,j6,j7
  character::cL*100
  !--EndDefs-----------------------------------------------
                                                                if(iP<6) return
   cL=''c ;write(iP,*) CL
   cL= &
   '    Number of polynomial coefficients <= 9999'// &
   '  for 7 variables up to order 20'c
   write(iP,*) CL
   cL= &
   '            Order of v1 = 0    2    4    6    8|'// &
   '  10   12   14   16   18   20'c
   write(iP,*) CL
   cL= &
   '                      ----------------- <-exact|unscaled->'// &
   ' ------------------'c
   write(iP,*) CL
                              nLine=1
   do i1=0,20,2                ; j1 = i1+1;
     do i2=i1,20,2             ; j2 = i2+1;
       do i3=i2,20,2           ; j3 = i3+1;
         do i4=i3,20,2         ; j4 = i4+1;
           do i5=i4,20,2       ; j5 = i5+1;
             do i6=i5,20,2     ; j6 = i6+1;
                                j  = j1*j2*j3*j4*j5*j6
                             if(J>9999) exit
                 cL=char(0);  iUsed = 0;           nCh = 22
                 write(cL(1:22),"(i3,1x,6i3.0\)") nLine,i1,i2,i3,i4,i5,i6
                                 if(nLine==  1) cL(17:22) = 'V2 = 0'
                                 if(nLine== 11) cL(14:19) = 'V3 = 0'
                                 if(nLine== 66) cL(11:16) = 'V4 = 0'
                                 if(nLine==132) cL( 8:13) = 'V5 = 0'
                                 if(nLine==161) cL( 5:10) = 'V6 = 0'
                                 if(nLine==171) cL( 4: 7) = 'V7=0'
               do i7=2,i6,2;              write(cL(nCh+1:nCh+5),"('     '\)")
                                                   nCh = nCh+5        ;enddo!i7
               do i7=i6,20,2;  j7 = i7+1;
                                  j  = j1*j2*j3*j4*j5*j6*j7
                               if(J>9999) exit
                               iUsed = iUsed+1
                                          write(cL(nCh+1:nCh+5),"(i5\)") j
                                                   nCh = nCh+5
               enddo!i7
                            if(iUsed==0) exit
                              nLine=nLine+1
                                                      cL(48:48)='|'
                                          write(iP,*) CL
             enddo!i6
           enddo!i5
         enddo!i4
       enddo!i3
     enddo!i2
   enddo!i1
   cL= &
   '                      ----------------- <-exact|unscaled->'// &
   ' ------------------'c
   write(iP,*) CL
   cL=''c ;write(iP,*) CL
                                                                         return
End Subroutine CoefficientCount
!---------------------------------------------------------------------------7 9

