!
!                                **********************************************
!                                *             MODULE pcomb                   *
!                                *  R. J. Purser, NOAA/NCEP/EMC               *
!                                *  jim.purser@noaa.gov                       *
!                                *                                            *
!                                **********************************************
!
! A variety of handy utilities related to the general theme of
! "Combinatorics".
!
! Last minor updates:
! Nov 2012 (simplify arg lists for some next_... routines and allow
!           option for the number sequences to begin with 0. Also
!           include a new "next_signs" variant of next_bins, since a
!           frequent application is to supply sign switches +1 and -1)
!
! DIRECT DEPENDENCIES
! Modules: pietc, pkind
!
!
!============================================================================
module pcomb
!============================================================================
use pkind, only: dp,dpi
use pietc, only: F,T
implicit none
private
public:: kilo,mega,giga,plus1,minus1,flag, &        ! <- the module variables
     factorial,binomial,multinomial,pmultinomial, &            !
     bernoulli,stirling,eulert,partnum, &                      ! <- functions
     abernoulli,dpartnum,dfactorial,rfactorial,lfactorial, &
     compare_bin,next_bin,next_num,int_to_num,num_to_int,next_comb, &
     loc_comb,next_perm,loc_perm,next_boxnum,next_dist,loc_dist,next_mdist, &
     next_signs,next_perm0,next_perm1,next_comb0,next_comb1, &
     compare_part,partable,partnums,loc_table_part,dualpart, & ! <-partitions
     next_sectorvector,count_sectorvector,sectormixers, & ! <- Hyperspheres
     qprime, & ! <- query whether number is prime
     taypow,taymul,tayval,taycom, & ! <- taylor series manipulations
     mulrat,invrat,divrat,addrat,subrat,euclid ! <- rational manipulations
integer,parameter                :: kilo=1000,mega=1000000,giga=1000000000
integer, allocatable,dimension(:):: kv,mv   ! <- used by next_sectorvector
real(dp),allocatable,dimension(:):: xs,dxs  ! <- used by next_sectorvector
integer,dimension(2)             :: plus1,minus1
logical                          :: flag=.true.
data plus1 / 1,1/
data minus1/-1,1/

interface bernoulli;     module procedure bernoulli;              end interface
interface abernoulli;    module procedure abernoulli;             end interface
interface factorial;     module procedure factorial4,factorial8;  end interface
interface dfactorial;    module procedure dfactorial;             end interface
interface rfactorial;    module procedure rfactorial;             end interface
interface lfactorial;    module procedure lfactorial;             end interface
interface binomial;      module procedure binomial4,binomial8;    end interface
interface multinomial;   module procedure multinomial;            end interface
interface pmultinomial;  module procedure pmultinomial;           end interface
interface stirling;      module procedure stirling4,stirling8;    end interface
interface eulert;        module procedure eulert4,eulert8;        end interface
interface partnum;       module procedure partnum;                end interface
interface dpartnum;      module procedure dpartnum;               end interface
interface compare_bin;   module procedure compare_bin;            end interface
interface next_bin;      module procedure next_bin,next_bins;     end interface
interface next_signs;    module procedure next_signs;             end interface
interface next_num;      module procedure next_num,next_nums;     end interface
interface int_to_num;    module procedure int_to_num;             end interface
interface num_to_int;    module procedure num_to_int;             end interface
interface next_comb;     module procedure next_comb;              end interface
interface next_comb0;    module procedure next_comb0;             end interface
interface next_comb1;    module procedure next_comb1;             end interface
interface loc_comb;      module procedure loc_comb;               end interface
interface next_perm;     module procedure next_perm,next_mperm;   end interface
interface next_perm0;    module procedure next_perm0;             end interface
interface next_perm1;    module procedure next_perm1;             end interface
interface loc_perm;      module procedure loc_perm;               end interface
interface next_boxnum;   module procedure next_boxnum;            end interface
interface next_dist;     module procedure next_dist;              end interface
interface loc_dist;      module procedure loc_dist;               end interface
interface next_mdist;    module procedure next_mdist;             end interface
interface loc_mdist;     module procedure loc_mdist;              end interface
interface compart_part;  module procedure compare_part;           end interface
interface partable;      module procedure partable;               end interface
interface partnums;      module procedure ipartnums4,ipartnums8;  end interface
interface loc_table_part;module procedure loc_table_part;         end interface
interface dualpart;      module procedure dualpart;               end interface
interface next_sectorvector;  module procedure next_sectorvector; end interface
interface count_sectorvector; module procedure count_sectorvector;end interface
interface sectormixers;  module procedure sectormixers;           end interface
interface qprime;        module procedure qprime,qprime_d;        end interface
interface taypow;        module procedure taypow;                 end interface
interface taymul;        module procedure taymul;                 end interface
interface tayval;        module procedure tayval;                 end interface
interface taycom;        module procedure taycom;                 end interface
interface mulrat;        module procedure mulrat,mulrat8;         end interface
interface invrat;        module procedure invrat,invrat8;         end interface
interface divrat;        module procedure divrat,divrat8;         end interface
interface addrat;        module procedure addrat,addrat8;         end interface
interface subrat;        module procedure subrat,subrat8;         end interface
interface euclid;        module procedure euclid,euclid8;         end interface

contains
!=============================================================================
function bernoulli(n)result(bn)!                                   [bernoulli]
!=============================================================================
! Deliver the positive even Bernoulli number of index 2n as a real(dp)
! using Seidel's method.
!=============================================================================
integer,intent(in ):: n
real(dp)           :: bn
!-----------------------------------------------------------------------------
real(dp),dimension(n):: b
!=============================================================================
call abernoulli(n,b)
bn=b(n)
end function bernoulli

!=============================================================================
subroutine abernoulli(n,b)!                                       [abernoulli]
!=============================================================================
! Use Seidel's method to compute the first n even Bernoulli
! numbers in an array, b.
!=============================================================================
integer,              intent(in ):: n
real(dp),dimension(n),intent(out):: b
!-----------------------------------------------------------------------------
real(dp),dimension(-n:n):: s
real(dp)                :: p4
integer                 :: i,j
!=============================================================================
s=0; s(0)=1
p4=4
b(1)=(2*s(0))/(p4*(p4-1))! = +1/6
do j=1,n-1
   p4=p4*4
   do i=1-j,j
      s(i)=s(i-1)+s(i)
   enddo
   do i=j-1,-j,-1
      s(i)=s(i+1)+s(i)
   enddo
   b(j+1)=(2*(j+1)*s(-j))/(p4*(p4-1))
   if(mod(j+1,2)==0)b(j+1)=-b(j+1)
enddo
end subroutine abernoulli

!============================================================================
function factorial4(n)result(f)!                                   [factorial]
!============================================================================
integer           :: f
integer,intent(IN):: n
!----------------------------------------------------------------------------
integer           :: i
!============================================================================
if(n<0)stop 'In factorial4; negative value of n'
if(n>12)stop &
  'In factorial4; n>12 too large; call dfactorial or rfactorial instead'
f=1; do i=1,n; f=f*i; enddo
end function factorial4
!============================================================================
function factorial8(n)result(f)!                                  [factorial]
!============================================================================
! Note: function name overloading with "factorial" is not possible owing to
! ambiguity, since all versions need the same single-integer argument.
integer(dpi)            :: f
integer(dpi),intent(IN ):: n
!----------------------------------------------------------------------------
integer                 :: i
!============================================================================
if(n<0)stop 'In factorial8; negative value of n'
if(n>20)stop &
  'In factorial8; n>20 too large for integer function; call rfactorial instead'
f=1; do i=1,n; f=f*i; enddo
end function factorial8


!============================================================================
function dfactorial(n)result(f)!                                 [dfactorial]
!============================================================================
! Note: function name overloading with "factorial" is not possible owing to
! ambiguity, since all versions need the same single-integer argument.
integer(dpi)       :: f
integer,intent(IN ):: n
!----------------------------------------------------------------------------
integer            :: i
!============================================================================
if(n<0)stop 'In dfactorial; negative value of n'
if(n>20)stop &
  'In dfactorial; n>20 too large for integer function; call rfactorial instead'
f=1; do i=1,n; f=f*i; enddo
end function dfactorial

!============================================================================
function rfactorial(n)result(f)!                                 [rfactorial]
!============================================================================
! Note: function name overloading with "factorial" is not possible owing to
! ambiguity, since all versions need the same single-integer argument.
real(dp)           :: f
integer,intent(IN ):: n
!----------------------------------------------------------------------------
integer            :: i
!============================================================================
if(n<0)stop 'In rfactorial; negative value of n'
f=1; do i=1,n; f=f*i; enddo
end function rfactorial

!=============================================================================
function Lfactorial(n)result(L)!                                  [lfactorial]
!=============================================================================
! Return the log of n! for integer n >= 0
!=============================================================================
use pietc, only: one,half,pi2
integer,intent(in ):: n
real(dp)           :: L
!-----------------------------------------------------------------------------
integer,parameter    :: nn=10
integer,dimension(nn):: num,den ! Coefficients in rising-exponent Stirling
integer              :: i,k
real(dp)             :: z,zp,pz
data num/   1,       1,      59,      29,     533,  &
         1577,  280361,   69311,36226519, 7178335/
data den/  12,      12,     360,      60,     280,  &
          168,    5040,     180,   11880,     264/
!=============================================================================
if(n<=50)then
   if(n<0)stop 'In Lfactorial; n<0 is out of bounds'
! Just sum the logs of the positive integers up to n
   L=0; do i=1,n; L=L+log(one*i); enddo
else
! For large n, apply Stirling, using 10 terms of convergent asymptotic series
   z=n+1
   L=(z-half)*log(z)-z+half*log(pi2)
   zp=z
   pz=one
   do k=1,nn
      zp=zp+one
      pz=pz*zp ! Current "rising exponent" of z
      L=L+num(k)/(den(k)*pz)
   enddo
endif
end function lfactorial

!============================================================================
function binomial4(n,m)result(f)!                                  [binomial]
!============================================================================
integer           :: f
integer,intent(IN):: n,m
!----------------------------------------------------------------------------
integer           :: i,h
!============================================================================
if(n<0 .or. m<-1 .or. m>n+1)stop 'In binomial4; invalid n or m'
if(m==n+1 .or. m==-1)then
   f=0
   return
endif
h=huge(1)
f=1
do i=1,min(m,n-m)
   if(f>h/(n+1-i))stop 'In binomial4; numbers are becoming too large'
   f=(f*(n+1-i))/i
enddo
end function binomial4

!============================================================================
function binomial8(n,m)result(f)!                                  [binomial]
!============================================================================
integer(dpi)           :: f
integer(dpi),intent(IN):: n,m
!----------------------------------------------------------------------------
integer(dpi)           :: i,h
!============================================================================
if(n<0 .or. m<-1 .or. m>n+1)stop 'In binomial8; invalid n or m'
if(m==n+1 .or. m==-1)then
   f=0
   return
endif
h=huge(1_dpi)
f=1
do i=1,min(m,n-m)
   if(f>h/(n+1-i))stop 'In binomial8; numbers are becoming too large'
   f=(f*(n+1-i))/i
enddo
end function binomial8

!=============================================================================
integer function multinomial(nm,m)!                              [multinomial]
!=============================================================================
! Return the generalized multinomial coefficient,
! c={m(1)+..m(nm); m(1), ...., m(nm) }
! in which all the arguments, m(i), need not be non-negative. When more
! one of them are negative, the result c is 0. Also, when one of them is
! negative but sum(m) is non-negative, the result c is again 0. However,
! when just one of m(i) and sum(m) are both negative numbers, the value of
! c returned is f times the c that one would obtain by replacing the
! negative m(i) by (-1-sum(m)), and f is +1 or -1 according to whether
! the original sum(m) and the original negative m(i) had same or different
! parity. This is the natural and consistent way to generalize the
! multinomial coefficients to all the integer arguments.
!=============================================================================
integer,              intent(IN ):: nm
integer,dimension(nm),intent(IN ):: m
!-----------------------------------------------------------------------------
integer                          :: i,f,s,cneg
integer,dimension(nm)            :: mp
!=============================================================================
multinomial=0
cneg=0
s=sum(m)
mp=m
do i=1,nm
   if(m(i)<0)then
      if(cneg>0)return
      cneg=i
   endif
enddo
f=1
if(cneg>0)then
   if(s>=0)return
   f=(-1)**(s-m(cneg))
   mp(cneg)=-s-1
endif
multinomial=f*pmultinomial(nm,mp)
end function multinomial

!=============================================================================
integer function pmultinomial(nm,m)!                            [pmultinomial]
!=============================================================================
! A version of the multinomial coefficient evaluation in which all the
! arguments, m, are non-negative.
!=============================================================================
integer,              intent(IN ):: nm
integer,dimension(nm),intent(IN ):: m
!-----------------------------------------------------------------------------
integer                          :: i,j,s
!=============================================================================
s=sum(m)
pmultinomial=1
do j=1,nm-1
   do i=1,m(j); pmultinomial=(pmultinomial*s)/i; s=s-1; enddo
enddo
end function pmultinomial

!=============================================================================
function stirling4(n,m)result(s)!                                   [stirling]
!=============================================================================
integer           :: s
integer,intent(in):: n,m
!-----------------------------------------------------------------------------
integer:: k,b
!=============================================================================
s=0
do k=0,m
   b=binomial(m,k)
   s=s+(-1)**(m-k)*b*k**n
enddo
s=s/factorial(m)
end function stirling4
!=============================================================================
function stirling8(n,m)result(s)!                                   [stirling]
!=============================================================================
integer(dpi)           :: s
integer(dpi),intent(in):: n,m
!-----------------------------------------------------------------------------
integer(dpi):: k,b
!=============================================================================
s=0
do k=0,m
   b=binomial(m,k)
   s=s+(-1)**(m-k)*b*k**n
enddo
s=s/factorial(m)
end function stirling8

!=============================================================================
function eulert4(n,m)result(e)!                                       [eulert]
!=============================================================================
! Evaluate the Euler triangle number E(n,m) indexed as in Riordan
!=============================================================================
integer           :: e
integer,intent(in):: n,m
!-----------------------------------------------------------------------------
integer:: k,b
!=============================================================================
e=0
do k=0,m-1
   b=binomial(n+1,k)
   e=e+b*(-1)**k * (m-k)**n
enddo
end function eulert4
!=============================================================================
function eulert8(n,m)result(e)!                                       [eulert]
!=============================================================================
integer(dpi)           :: e
integer(dpi),intent(in):: n,m
!-----------------------------------------------------------------------------
integer(dpi):: k,b
!=============================================================================
e=0
do k=0,m-1
   b=binomial(n+1,k)
   e=e+b*(-1)**k * (m-k)**n
enddo
end function eulert8

!=============================================================================
function partnum(n)result(p)!                                       [partnum]
!=============================================================================
integer,       intent(IN):: n
integer(4),dimension(0:n):: pn
integer(4)               :: c,p
!=============================================================================
call partnums(n,pn,c)
p=pn(n)
end function partnum

!=============================================================================
function dpartnum(n)result(p)!                                      [dpartnum]
!=============================================================================
integer,         intent(IN):: n
integer(dpi),dimension(0:n):: pn
integer(dpi)               :: c,p
!=============================================================================
call partnums(n,pn,c)
p=pn(n)
end function dpartnum

!============================================================================
subroutine compare_bin(n,bin1,bin2,leg)!                        [compare_bin]
!============================================================================
! Compare two binary numbers, bin1 and bin2, of n bits. 
! If bin1 <  bin2,  return leg=-1
! If bin1 == bin2,  return leg= 0
! If bin1 >  bin2,  return leg=+1
!
! The most significant bit is assumed to be bin(n)
! (Actually, this routine will also work with bases other than 2.)
!============================================================================
integer,             intent(IN ):: n
integer,dimension(n),intent(IN ):: bin1,bin2
integer,             intent(OUT):: leg
!---------------------------------------------------------------------------
integer :: i
!============================================================================
leg=0
do i=n,1,-1
   if(bin1(i)<bin2(i))then
      leg=-1
      return
   elseif(bin1(i)>bin2(i))then
      leg=1
      return
   endif
enddo
end subroutine compare_bin
     
!============================================================================
subroutine next_bin(n,bin) !                                       [next_bin]
!============================================================================
! Replace bin with the NEXT binary number of n bits
! If there is no "next", return a "TRUE" flag, otherwise it remains .FALSE.
! !!!!! PHASE OUT IN FAVOR OF EITHER NEXT_BINS OR NEXT_SIGNS
!============================================================================
integer,             intent(IN   ):: n
integer,dimension(n),intent(INOUT):: bin
!----------------------------------------------------------------------------
integer                           :: i
!============================================================================
print '(''IN NEXT_BIN; REPLACE BY SIMPLER NEXT_BINS WITH SIMPLER ARG LIST'')'
if(n<=0)then; flag=T; return; endif
if(flag)then; bin=0; flag=F
else
   do i=1,n
      if(bin(i)==1)cycle
      bin(i)=1
      bin(1:i-1)=0
      return
   enddo
   flag=T
endif
end subroutine next_bin
!============================================================================
subroutine next_bins(bin) !                                       [next_bin]
!============================================================================
! Replace bin with the NEXT binary number of n bits
! If there is no "next", return a "TRUE" flag, otherwise it remains .FALSE.
!============================================================================
integer,dimension(:),intent(INOUT):: bin
!----------------------------------------------------------------------------
integer                           :: i,n
!============================================================================
n=size(bin)
if(n<=0)then; flag=T; return; endif
if(flag)then; bin=0; flag=F
else
   do i=1,n
      if(bin(i)==1)cycle
      bin(i)=1
      bin(1:i-1)=0
      return
   enddo
   flag=T
endif
end subroutine next_bins

!===========================================================================+
subroutine next_signs(s)!                                        [next_signs]
!===========================================================================+
! Adapted from next_bin, this routine delivers the next sequence of "signs"
! in the form of a string of +1 and -1 (in place of the binary digits of 
! next_bin). The array size is obtaied from the SIZE intrinsic function,
! which simpifies calling this routine.
!============================================================================
integer,dimension(:),intent(inout):: s
!----------------------------------------------------------------------------
integer:: i,n
!============================================================================
n=size(s)
if(n==0)then; flag=T; return; endif
if(flag)then; s=1; flag=F
else;do i=1,n; if(s(i)==-1)cycle; s(i)=-1; s(1:i-1)=1; return; enddo; flag=T
endif
end subroutine next_signs

!============================================================================
subroutine next_num(n,r,digits)!                                   [next_num]
!============================================================================
! Replace digits with the NEXT radix-r number of n digits
!============================================================================
integer,             intent(IN   ):: n,r
integer,dimension(n),intent(INOUT):: digits
!----------------------------------------------------------------------------
integer                           :: i
!============================================================================
print'('' IN NEXT_NUM; REPLACE WITH NEXT_NUMS FOR SIMPLER ARG LIST'')'
if(n<=0)then
   flag=T
   return
endif
if(flag)then
   digits=0
   flag=F
else
   do i=1,n
      if(digits(i)==r-1)cycle
      digits(i)=digits(i)+1
      digits(1:i-1)=0
      return
   enddo
   flag=T
endif
end subroutine next_num
!============================================================================
subroutine next_nums(r,digits)!                                    [next_num]
!============================================================================
! Replace digits with the NEXT radix-r number of n digits
!============================================================================
integer,             intent(IN   ):: r
integer,dimension(:),intent(INOUT):: digits
!----------------------------------------------------------------------------
integer                           :: i,n
!============================================================================
n=size(digits)
if(n<=0)then
   flag=T
   return
endif
if(flag)then
   digits=0
   flag=F
else
   do i=1,n
      if(digits(i)==r-1)cycle
      digits(i)=digits(i)+1
      digits(1:i-1)=0
      return
   enddo
   flag=T
endif
end subroutine next_nums

!============================================================================
subroutine int_to_num(int,n,r,digits)!                           [int_to_num]
!============================================================================
! Convert a positive integer, int,to its equivalent radix-r n-digit number in
! the array, digits. If this fails, set flag to .TRUE.
!============================================================================
integer,             intent(IN ):: int,n,r
integer,dimension(n),intent(OUT):: digits
!----------------------------------------------------------------------------
integer                         :: i,c
!============================================================================
flag=(r<=0 .or. int<0 .or. n<=0); if(flag)return
c=int
do i=1,n
   digits(i)=mod(c,r); c=c/r
enddo
flag=c>0 ! Set flag when there is still a residual (conversion incomplete).
end subroutine int_to_num

!============================================================================
subroutine num_to_int(int,n,r,digits)!                           [num_to_int]
!============================================================================
! Convert a positive radix-r n-digit number in the array, digits, to its 
! equivalent integer, int. If this fails, set the flag to .TRUE.
!============================================================================
integer,             intent(OUT):: int
integer,             intent(IN ):: n,r
integer,dimension(n),intent(IN ):: digits
!----------------------------------------------------------------------------
integer                         :: i,h
!============================================================================
flag=(r<=0 .or. n<=0); if(flag)return
h=huge(1)/r
int=0
flag=T
do i=n,1,-1
   if(int+1>h)return
   int=int*r+digits(i)
enddo
flag=F
end subroutine num_to_int

!============================================================================
subroutine next_comb(nc,n,comb)!!!                                [next_comb]
!============================================================================
! If possible, replace comb with the next of the ordered combinations of n 
! integers from the set, {1,...,nc} and return a "false" flag. If there is 
! no next one, return a "true" flag.
! !!!!!!!!!!! PHASE OUT IN FAVOR OF NEXT_COMB0 AND NEXT_COMB1
!============================================================================
integer,             intent(IN   ):: n,nc
integer,dimension(n),intent(INOUT):: comb
!----------------------------------------------------------------------------
integer                           :: i,j
!============================================================================
if(flag)then
   do i=1,n; comb(i)=i; enddo
   flag=F
else
   do i=n,1,-1
      if(comb(i)+n-I==nc)cycle
      comb(i)=comb(i)+1
      do j=i+1,n
         comb(j)=comb(j-1)+1
      enddo
      return
   enddo
   flag=T
endif
end subroutine next_comb

!============================================================================
subroutine next_comb0(nc,comb)!                                  [next_comb0]
!============================================================================
! If possible, replace comb with the next of the ordered combinations of n 
! integers from the set, {0,...,nc-1} and return a "false" flag. If there is 
! no next one, return a "true" flag.
!============================================================================
integer,             intent(IN   ):: nc
integer,dimension(:),intent(INOUT):: comb
!----------------------------------------------------------------------------
integer                           :: i,j,n,ncm
!============================================================================
n=size(comb); ncm=nc-1
if(flag)then
   do i=1,n; comb(i)=i-1; enddo
   flag=F
else
   do i=n,1,-1
      if(comb(i)+n-I==ncm)cycle
      comb(i)=comb(i)+1
      do j=i+1,n
         comb(j)=comb(j-1)+1
      enddo
      return
   enddo
   flag=T
endif
end subroutine next_comb0

!============================================================================
subroutine next_comb1(nc,comb)!                                  [next_comb1]
!============================================================================
! If possible, replace comb with the next of the ordered combinations of n 
! integers from the set, {1,...,nc} and return a "false" flag. If there is 
! no next one, return a "true" flag.
!============================================================================
integer,             intent(IN   ):: nc
integer,dimension(:),intent(INOUT):: comb
!----------------------------------------------------------------------------
integer                           :: i,j,n
!============================================================================
n=size(comb)
if(flag)then
   do i=1,n; comb(i)=i; enddo
   flag=F
else
   do i=n,1,-1
      if(comb(i)+n-I==nc)cycle
      comb(i)=comb(i)+1
      do j=i+1,n
         comb(j)=comb(j-1)+1
      enddo
      return
   enddo
   flag=T
endif
end subroutine next_comb1

!============================================================================
subroutine loc_comb(n,m,comb,index)!                               [loc_comb]
!============================================================================
! Locate the index assigned to this combination of m integers from {1,..,n}
! in the standard ordering, as produced by iterated next_comb.
!============================================================================
integer,             intent(IN ):: n,m
integer,dimension(m),intent(IN ):: comb
integer,             intent(OUT):: index
!---------------------------------------------------------------------------
integer                         :: i
!============================================================================
index=binomial(n,m)
do i=1,m
   index=index-binomial(n-comb(i),m+1-i)
enddo
end subroutine loc_comb

!============================================================================
subroutine next_perm(n,b)!                                        [next_perm]
!============================================================================
! If possible, replace the permutation sequence, b, of the first n integers
! by the next one in the standard ordering of the list, and return a "false"
! flag. If there is no next one, (b is already the last one) return a "true"
! flag.
! !!!!!!!!!!! PHASE OUT IN FAVOR OF NEXT_PERM0 AND NEXT_PERM1
!============================================================================
integer,             intent(IN   ):: n
integer,dimension(n),intent(INOUT):: b
!----------------------------------------------------------------------------
integer                           :: i,j,tt
!============================================================================
if (flag)then
   do i=1,n; b(i)=i; enddo
   flag=F
else
   do i=n,2,-1
      if(b(i-1)>b(i))cycle
      tt=b(i-1) ! This element will be replaced by the next higher from tail:
      do j=n,i,-1
         if(b(j)<tt)cycle
         b(i-1)=b(j)
         b(j)=tt
         b(i:n)=b(n:i:-1) ! <- Flip the order of new tail
         return
      enddo
   enddo
   flag=T  ! To signify that all permutations have already been exhausted
endif
end subroutine next_perm

!============================================================================
subroutine next_mperm(nn,b,dom)!                                  [next_perm]
!============================================================================
integer,              intent(IN   ):: nn
integer,dimension(nn),intent(INOUT):: b
integer,dimension(nn),intent(IN   ):: dom
!----------------------------------------------------------------------------
integer                            :: k,l,m,n,iseg,k1,k2,tt
!============================================================================
if(flag)then
   if(maxval(dom)/=1 .or. minval(dom)<0)stop 'In next_mperm; invalid dom'
   if(dom(nn)/=1)stop 'In next_mperm; dom must terminate in a "1")'
   call next_perm(nn,b)
else
   k=nn
   do iseg=1,nn
      flag=F
      if(k<1)exit
      m=k
      do k=m-1,1,-1
         if(dom(k)==1)exit
      enddo
      L=k+1
      n=m+1-L
      call next_perm(n,b(l:m))
      if(.not.flag)return
      do k1=L,(L+m)/2; k2=L+m-k1
         tt   =b(k1)
         b(k1)=b(k2)
         b(k2)=tt
      enddo
   enddo
   flag=T
endif
end subroutine next_mperm
 
!============================================================================
subroutine next_perm0(b)!                                        [next_perm0]
!============================================================================
! For n=size(n)
! If possible, replace the permutation sequence, b, of the n integers,
! 0,1,....n-1, by the next permutation in the standard ordering of the list, 
! and return a "false" (F) flag. If there is no next one, (b is already the 
! last one) return a "true" (T) flag.
! This routine and the almost identical next_perm1 routine, are derived from
! next_perm, but are a little simpler to call.
!============================================================================
integer,dimension(:),intent(INOUT):: b
!----------------------------------------------------------------------------
integer                           :: i,j,n,tt
!============================================================================
n=size(b)
if(flag)then
   do i=1,n; b(i)=i-1; enddo
   flag=F
else
   do i=n,2,-1
      if(b(i-1)>b(i))cycle
      tt=b(i-1) ! This element will be replaced by the next higher from tail:
      do j=n,i,-1
         if(b(j)<tt)cycle
         b(i-1)=b(j)
         b(j)=tt
         b(i:n)=b(n:i:-1) ! <- Flip the order of new tail
         return
      enddo
   enddo
   flag=T  ! To signify that all permutations have already been exhausted
endif
end subroutine next_perm0
     
!============================================================================
subroutine next_perm1(b)!                                        [next_perm1]
!============================================================================
! Trivial variant of NEXT_PERM0, but permuting 1,2,...n.
!============================================================================
integer,dimension(:),intent(INOUT):: b
!----------------------------------------------------------------------------
integer                           :: i,j,n,tt
!============================================================================
n=size(b)
if(flag)then
   do i=1,n; b(i)=i; enddo
   flag=F
else
   do i=n,2,-1
      if(b(i-1)>b(i))cycle
      tt=b(i-1) ! This element will be replaced by the next higher from tail:
      do j=n,i,-1
         if(b(j)<tt)cycle
         b(i-1)=b(j)
         b(j)=tt
         b(i:n)=b(n:i:-1) ! <- Flip the order of new tail
         return
      enddo
   enddo
   flag=T  ! To signify that all permutations have already been exhausted
endif
end subroutine next_perm1
     
!============================================================================
subroutine loc_perm(n,perm,index)!                                 [loc_perm]
!============================================================================
! Locate the index assigned to this permutation of integers from {1,..,n}
! in the standard ordering, as produced by iterated next-perm.
!============================================================================
integer,             intent(IN ):: n
integer,dimension(n),intent(IN ):: perm
integer,             intent(OUT):: index
!----------------------------------------------------------------------------
integer                         :: i,j,k
!============================================================================
index=1
do i=1,n
   k=perm(i)
   do j=1,i
      if(perm(j)<=perm(i))k=k-1
   enddo
   index=index+k*factorial(n-i)
enddo
end subroutine loc_perm

!============================================================================
subroutine next_boxnum(m,n,num)!                                 [next_boxnum]
!============================================================================
! Initialize (flag==.true. on input) or increment (flag==.false. on input)
! the n-vector of "box-numbers" of norm, m. If successfully incremented (or
! initialized) return flag==.false. but if the series of box numbers is
! already ended, return flag==.true. as an escape flag. The "box numbers" 
! are the integer vectors with L_infty norm of m, i.e., they form a hollow
! n-dimensional "box" in their totality.
!============================================================================
integer,             intent(IN   ):: m,n
integer,dimension(n),intent(INOUT):: num
!----------------------------------------------------------------------------
integer                           :: i
!============================================================================
! num is of vanishing length:
if(n<=0)then
   flag=.true.
   return
endif

! n>0 and flag is true; initialize num:
if(flag)then
   num=-m
   flag=.false.
   return
endif

! flag is false, n==1:
if(n==1)then
   flag=(num(1)==m)
   num(1)=m
   return
endif

! flag is false, n>1:
do i=1,n
   if(num(i)==m)cycle
   if(i==1 .and. minval(num(2:n))>-m .and. maxval(num(2:n))<m)then
      num(i)=m
   else
      num(i)=num(i)+1
      num(:i-1)=-m
   endif
   return
enddo

! All num(i) have already attained the max allowed value, m; terminate series:
flag=.true.
end subroutine next_boxnum

!============================================================================
subroutine next_dist(nc,n,dist)!                                  [next_dist]
!============================================================================
! Find the next distribution in the list of distributions of nc objects into
! n boxes.
!============================================================================
integer,             intent(IN   ):: nc,n
integer,dimension(n),intent(INOUT):: dist
!----------------------------------------------------------------------------
integer                           :: s,i
!============================================================================
if(flag)then
  if(n>0)then; dist=0; dist(1)=nc; endif
  flag=.false.
else
   flag=(nc==0); if(flag)return
   s=dist(n)+1
   dist(n)=0
   do i=n-1,1,-1
      if(dist(i)>0)then
         dist(i)=dist(i)-1
         dist(i+1)=s
         return
      endif
   enddo
   flag=.true.
endif
end subroutine next_dist

!============================================================================
subroutine loc_dist(nc,n,dist,index)!                              [loc_dist]
!============================================================================
! Locate the index assigned to this distribution of integers summing to nc
! over n boxes, where the first index of the list is 1.
!============================================================================
integer,             intent(IN ):: nc,n
integer,dimension(n),intent(IN ):: dist
integer,             intent(OUT):: index
!----------------------------------------------------------------------------
integer                         :: i,ic,r
!===========================================================================
index=1
r=nc
do i=1,n-1
   ic=n-i
   r=r-dist(i)
   index=index+binomial(ic+r-1,ic)
enddo
end subroutine loc_dist

!============================================================================
subroutine next_mdist(nnc,n,dist)!                               [next_mdist]
!============================================================================
! Find the next member of the list of multi-sum distributions of the
! nonnegative integers with sum not exceeding nnc distributed over
! n boxes. If no further members exist, set the escape flag to .true.
!============================================================================
integer,             intent(IN   ):: nnc,n
integer,dimension(n),intent(INOUT):: dist
!----------------------------------------------------------------------------
integer                           :: nc
!============================================================================
if(flag)then
   dist=0
   flag=.false.
else
   nc=sum(dist)
   if(dist(n)==nc)then
      flag=(nc==nnc); if(flag)return
      nc=nc+1
      dist=0; dist(1)=nc
      return
   endif
   call next_dist(nc,n,dist)
   if(flag)stop 'anomalous results in next_mdist'
endif
end subroutine next_mdist

!============================================================================
subroutine loc_mdist(n,dist,index)!                               [loc_mdist]
!============================================================================
integer,             intent(IN ):: n
integer,dimension(n),intent(IN ):: dist
integer,             intent(OUT):: index
!----------------------------------------------------------------------------
integer                         :: nc
!============================================================================
nc=sum(dist)
call loc_dist(nc,n,dist,index)
index=index+binomial(n+nc-1,n)
end subroutine loc_mdist

!=============================================================================
subroutine compare_part(n,part1,part2,leg)!                     [compare_part]
!=============================================================================
integer,             intent(IN ):: n
integer,dimension(n),intent(IN ):: part1,part2
integer,             intent(OUT):: leg
!-----------------------------------------------------------------------------
integer            :: i,m1,m2,n1,n2
!=============================================================================
n1=sum(part1); n2=sum(part2)
if    (n1<n2)then;                leg=-1; return
elseif(n1>n2)then;                leg= 1; return;  endif
do m1=n,1,-1; if(part1(m1)>0)exit;  enddo
do m2=n,1,-1; if(part2(m2)>0)exit;  enddo
if    (m1<m2)then;                leg=-1; return
elseif(m1>m2)then;                leg= 1; return;  endif
do i=1,m1
   if    (part1(i)>part2(i))then; leg=-1; return
   elseif(part1(i)<part2(i))then; leg= 1; return; endif
enddo
leg=0
end subroutine compare_part

!=============================================================================
subroutine partable(nn,nk, pn,ls,ns,pk,ok,mk)
!=============================================================================
! List and index the partition table for all integers up to nn.
! nk     : (found from subroutine partnums) the size of the cumulative table.
! pn(0:n): The partition number for n items (eg, pn(0)=1, pn(5)=7)
! ls(m,0:n): pointer to index k of 1st entry for partitions of n into m sets.
! ns(m,0:n): number of such entries
! pk(:,k): the table of partitions for eack k
! ok(:,k): the occupancy table for each k
! mk(k)  : The m of this k.
!=============================================================================
integer,                     intent(IN ):: nn,nk
integer,dimension(0:nn),     intent(OUT):: pn
integer,dimension(0:nn,0:nn),intent(OUT):: ls,ns
integer,dimension(nn,nk),    intent(OUT):: pk,ok
integer,dimension(nk),       intent(OUT):: mk
!-----------------------------------------------------------------------------
integer                                 :: j,k,l,m,n,p,r,k1
!============================================================================= 
pk=0; ok=0; ls=0; ns=0
ls(:,0)=1; ns(0,0)=1; mk(1)=0; pn(0)=1
k=2
do n=1,nn
   ls(0,n)=k; ns(0,n)=0
   ls(1,n)=k; ns(1,n)=1
   pk(1,k)=n; ok(n,k)=1
   mk(k)=1
   do m=2,n
      ls(m,n)=k+1
      do l=1,n
         if(m*l>n)exit
         r=n-l
         do k1=ls(m-1,r),ls(m-1,r)+ns(m-1,r)-1
            if(pk(1,k1)<l)cycle
            k=k+1
            pk(1,k)=l
            pk(2:nn,k)=pk(1:nn-1,k1)
            do j=1,m
               p=pk(j,k)
               ok(p,k)=ok(p,k)+1
            enddo
            mk(k)=m
         enddo
      enddo
      ns(m,n)=k+1-ls(m,n)
   enddo
   k=k+1
   ls(n+1:nn,n)=k
   pn(n)=k-ls(1,n)
enddo
k=k-1
if(k/=nk)stop 'In partable; k /= nk'
! Reverse the order of nonzero digits in each pk(:,k) to put them in 
! decreasing order (i.e., most significant first).
do k=1,nk
   do l=nn,1,-1
      if(pk(l,k)/=0)then
         pk(1:L,k)=pk(L:1:-1,k)
         exit
      endif
   enddo
enddo

end subroutine partable

!=============================================================================
subroutine ipartnums8(nn,pn,clast)!                                 [partnums]
!=============================================================================
! For a given nn, return the partition numbers in array pn for integers up
! to and including nn. Also return the size k required for a cumulative list
! of all these partitions, clast= sum_{i=1}^{nn} pn(i). 
!
! This is believed to be an extremely efficient routine, scaling as nn**(3/2).
! But note that the recurrence involved in this algorithm is numerically
! unstable if done using real numbers (with round-off error) so, if pn is
! required for numbers nn larger than can be accommodated by integer(dpi),
! one should use the "lpartnums" routine in conjunction with the "plint"
! long-integer utilities (not combined with the present module)., or use a
! less efficient (nn**2) but numerically stable routine like logparts to
! obtain a floating-point approximation to the logarithm of pn.
!=============================================================================
integer,                     intent(IN ):: nn
integer(dpi),dimension(0:nn),intent(OUT):: pn
integer(dpi),                intent(OUT):: clast
!-----------------------------------------------------------------------------
integer                     :: j,k,l,n,s
integer(dpi),dimension(0:nn):: c
!=============================================================================
if(nn>360)then
  print '('' WARNING IN PARTNUMS: nn is too large; results not guaranteed!!'')'
  read(*,*)
endif
c(0)=1
clast=1
pn(0)=1
do n=1,nn
   s=1
   j=n
   c(n)=1
   do k=1,n
      j=j-k
      if(j<0)exit
      l=j-k
      if(l>=0)c(j)=c(j)+c(L)
      c(n)=c(n)+s*c(j)
      s=-s
   enddo
   pn(n)=c(n)-clast
   clast=c(n)
enddo
end subroutine ipartnums8
!=============================================================================
subroutine ipartnums4(nn,pn,clast)!                                 [partnums]
!=============================================================================
integer,                   intent(IN ):: nn
integer(4),dimension(0:nn),intent(OUT):: pn
integer(4),                intent(OUT):: clast
!-----------------------------------------------------------------------------
integer                   :: j,k,l,n,s
integer(4),dimension(0:nn):: c
!=============================================================================
c(0)=1
clast=1
pn(0)=1
do n=1,nn
   s=1
   j=n
   c(n)=1
   do k=1,n
      j=j-k
      if(j<0)exit
      l=j-k
      if(l>=0)c(j)=c(j)+c(L)
      c(n)=c(n)+s*c(j)
      s=-s
   enddo
   pn(n)=c(n)-clast
   clast=c(n)
enddo
end subroutine ipartnums4

!=============================================================================
subroutine loc_table_part(n,nk,pk,item,index)
!=============================================================================
integer,                intent(IN ):: nk,n
integer,dimension(n,nk),intent(IN ):: pk
integer,dimension(n),   intent(IN ):: item
integer,                intent(OUT):: index
!-----------------------------------------------------------------------------
integer                            :: indexa,indexb,L,leg
!=============================================================================
index=0; flag=(nk<=0); if(flag)return
indexa=1; indexb=nk
call compare_part(n,item,pk(:,indexa),leg); flag=(leg<0)
if(flag)then
   index=indexa-1
   return
elseif(leg==0)then
   index=indexa
   return
endif

call compare_part(n,item,pk(:,indexb),leg); flag=(leg>0)
if(flag .or. leg==0)then
   index=indexb
   return
endif
do L=1,n
   index=(indexa+indexb)/2
   flag=(index==indexa); if(flag)return
   call compare_part(n,item,pk(:,index),leg); if(leg==0)return
   if(leg<0)then; indexb=index
            else; indexa=index;  endif
enddo
stop 'In loc_table_part; full number of iterations used - check the logic'
end subroutine loc_table_part

!=============================================================================
subroutine dualpart(n,p,q)
!=============================================================================
! Given a partition, p, return the dual partition,q
!=============================================================================
integer,             intent(IN ):: n
integer,dimension(n),intent(IN ):: p
integer,dimension(n),intent(OUT):: q
!-----------------------------------------------------------------------------
integer:: i,j,k
!=============================================================================
k=0
do i=n,1,-1
   if(p(i)>k)then
      do j=k+1,p(i)
         q(j)=i
      enddo
      k=p(i)
   endif
enddo
do j=k+1,n
   q(j)=0
enddo
end subroutine dualpart

!=============================================================================
subroutine next_sectorvector(m,n,  iv, xv,ws)
!=============================================================================
! Routine for delivering the next unit vector on a unit hypersphere embedded
! in a euclidean space of m >= 2 amongst a set that occupies a sector that
! serves as a fundamental region of the group of relection and rotation 
! symmetries of the concentric embedded hypercube aligned with the coordinates.
! When FLAG is given as .TRUE> (T), initialize the status arrays and deliver
! the first vector, resetting the flag to 'FALSE' (F). At subsequent calls,
! deliver the next vector of the fundamental set, or, if there are no more to
! be delivered, set the flag to T. The vectors a found by gridding the 
! positive orthogonal sector (quadrant, octant, etc ) of the hypercube of
! half-width n with a cartesian grid staggered wrt to the zero-coordinates
! in each of the m-1 directions. The indices of such points run from 1 to n
! towards edges, and the fundamental sector is chosen to be the one where
! the grid indices satisfy i1>=i2>=i3 etc. The m-dimension vector on the
! fundamental sector of the face with x(m)=n is projected back onto the
! unit hypersphere and the jacobian of this mapping is used to obtain the
! appropriate quadrature weight for using these vectors in the integration of
! a function over the sphere's image of the fundamental sector. The weights
! assigned to grid points lying on the diagonal boundaries of the sector
! are appropriately reduced to make them correct for whole-sphere 
! integrations obtained by replication of the fundamental domain's
! sampling vectors (2**m)*m! times necessary to complete the full coverage.
!
! The grid separation on the hypercube is not uniform, but is chosen to have
! uniform angular increments along the cartesian center lines. The real arrays 
! xs and dxdts are tables that define the positions and 1-D jacobian 
! of the grid respectively, where the jacobian is with respect to the 
! central angle coordinate. These table are computed on first call to
! NEXT_SECTORVECTOR and must not be overwritten between calls. The status
! integer vectors IV (the grid indices at the completiong of the previous
! step), KV (running counts of the numbers of consecutive iv that are
! equal), and MV (the multiplicity of overlapping replicates of the grid
! point under the application of the hypercube's symmetry group) are 
! both initialized here and updated here at each call but, again, must not
! be overwritten between calls to this routine. XV is the delivered next
! vector, and WS is its associated quadrature weight that takes care of
! the projection jacobian and the overlapping replicates mentioned above.
! Thus, if f(x) is a function defined over the m-dimensional space x, then the
! integral of f over the "surface" of the unit hypersphere is approximated
! by summing ws*f(G(xv)), where G denotes the set of all group elements
! applied to xv to construct the symmetry replicates needed to cover the
! whole sphere, and not just the fundamental secto covered by xv by itself.
!=============================================================================
use pietc, only: F,one,half,pih
integer,               intent(in   ):: m,n
integer,dimension(m-1),intent(inout):: iv
real(dp),dimension(m), intent(  out):: xv
real(dp),              intent(  out):: ws
!-----------------------------------------------------------------------------
real(dp)          :: t2,dis
integer           :: i,i1,mm,L1,L1m,L2
!=============================================================================
mm=m-1
if(flag)then
   if(allocated(xs))deallocate(xs,dxs,kv,mv)
   allocate( xs(n),dxs(n),kv(m-1),mv(m-1))
   do i1=1,n; t2=(i1-half)*pih/n; xs(i1)=tan(t2*half)
      dxs(i1)=pih/(n*(one+cos(t2)))
   enddo
   iv=1; i1=1; do i=1,mm; kv(i)=i; i1=i1*i; mv(i)=i1; enddo
   flag=F
else
! Assume that xs and dxs to be given with correct entries and assume that
! the status vectors iv, kv, mv are also valid and mutually consistent.
! If possible, increment index vector iv by one step, updating kv and mv
! as necessary; if not possible, set flag and return.
   flag=(iv(mm)==n)
   if(flag)then; deallocate(xs,dxs,kv,mv); return; endif
   do L1=mm,1,-1
      iv(L1)=iv(L1)+1
      L1m=L1-1
      if(L1>1 .and. iv(L1)> iv(L1m))cycle
      if(L1>1 .and. iv(L1)==iv(L1m))then
         kv(L1)=kv(L1m)+1; mv(L1)=mv(L1m)*kv(L1)
      endif
      do L2=L1+1,mm; iv(L2)=1; kv(L2)=L2-L1; mv(L2)=mv(L2-1)*kv(L2); enddo
      exit
   enddo
endif
xv(m)=one; ws=one
do i=1,mm; i1=iv(i); xv(i)=xs(i1); ws=ws*dxs(i1); enddo
dis=sqrt(dot_product(xv,xv)); xv=xv/dis; ws=ws/(mv(mm)*dis**m)

end subroutine next_sectorvector

!=============================================================================
subroutine count_sectorvector(m,n,count)
!=============================================================================
integer,intent(in ):: m,n
integer,intent(out):: count
count=binomial(m+n-2,m-1)
end subroutine count_sectorvector

!=============================================================================
subroutine sectormixers(m,nperm,nbin,perms,bins)
!=============================================================================
! Construct the arrays PERMS and BINS that can be used to shuffle the 
! components of a set of vectors covering the symmetry group's fundamental
! sector of the surface of a hypercube in m dimensions, so that the resulting
! signed component-permuted vectors collectively cover the entire surface,
! or do the same thing with the vectors normalized to |v|=1 in order to
! cover the surface of the unit hypersphere embedded in this space.
! See also subroutine next_sectorvector which can fill the fundamental
! sector with the requisite vectors.
!=============================================================================
integer,                   intent(in ):: m,nperm,nbin
integer,dimension(m,nperm),intent(OUT):: perms
integer,dimension(m,nbin) ,intent(OUT):: bins
!-----------------------------------------------------------------------------
integer,dimension(m):: bin,perm
integer             :: k
!=============================================================================
if(m>10 .or.m<2)stop 'in sectormixers; m is outside the reasonable range'
if(nperm/=factorial(m))stop 'In sectormixers; nperm is not consistent with m'
if(nbin /=2**m)stop 'In sectormixers; nbin is not consistent with m'
flag=T
do k=1,giga; call next_signs(bin);   if(flag)exit; bins(:,k)=bin; enddo
do k=1,giga; call next_perm1(perm); if(flag)exit; perms(:,k)=perm;   enddo
flag=F
end subroutine sectormixers

!=============================================================================
subroutine qprime(m,prime) !                                          [qprime]
!=============================================================================
! Determine whether or not the positive integer m is prime.
! PRIME==.TRUE. when m is prime.
!=============================================================================
integer,intent(IN ):: m
logical,intent(OUT):: prime
!-----------------------------------------------------------------------------
integer            :: n,i
real               :: r
!=============================================================================
prime=(m<4 .or. m==5)
if(prime .or. mod(m,2)==0 .or. mod(m,3)==0)return
r=m
n=sqrt(r+.5)
do i=5,n,6
   if(mod(m,i)  ==0)return
   if(mod(m,i+2)==0)return
enddo
prime=.true.
end subroutine qprime

!=============================================================================
subroutine qprime_d(mm,prime) !                                       [qprime]
!=============================================================================
! Determine whether or not the positive integer m is prime.
! PRIME==.TRUE. when m is prime.
!=============================================================================
integer(dpi),intent(IN ):: mm
logical,     intent(OUT):: prime
!-----------------------------------------------------------------------------
integer                :: n,i
integer(dpi)           :: ii,iim,iip
real(dp)               :: r
!=============================================================================
prime=(mm<4 .or. mm==5)
if(prime .or. mod(mm,2_dpi)==0 .or. mod(mm,3_dpi)==0)return
r=mm
n=sqrt(r+.5)
do i=5,n,6
   iim=i
   iip=i+2
   if(mod(mm,iim)  ==0)return
   if(mod(mm,iip)==0)return
enddo
prime=.true.
end subroutine qprime_d

!=============================================================================
subroutine taypow(n,pow,q,qp)
!=============================================================================
! (+/-Integer power of function)
! Use partitions and generalized multinomomial coefficients to derive the
! Taylor expansion qp that represents a function Q raised to the integer power,
! pow, where the Taylor expansion of Q itself has coefficients q. pow need not
! be positive.
!=============================================================================
integer,                intent(IN ):: n,pow
real(dp),dimension(0:n),intent(IN ):: q
real(dp),dimension(0:n),intent(OUT):: qp
!-----------------------------------------------------------------------------
integer,dimension(0:n,0:n)        :: ls,ns
integer,dimension(0:n)            :: pn
integer,dimension(0:n)            :: oka
integer                           :: i,jn,k,nk
integer,allocatable,dimension(:,:):: pk,ok
integer,allocatable,dimension(:)  :: mk
real(dp)                          :: s,tt
!=============================================================================
call partnums(n,pn,nk)
allocate(  pk(n,nk), ok(n,nk), mk(nk) )
call partable(n,nk, pn,ls,ns, pk,ok,mk)
qp(0)=q(0)**pow
do jn=1,n
   s=0
   do k=ls(1,jn),ls(1,jn)+pn(jn)-1
      oka(1:jn)=ok(1:jn,k); oka(0)=pow-sum(ok(1:jn,k))
      tt=1
      do i=0,jn
         if(oka(i)/=0)tt=tt*q(i)**oka(i)
      enddo
      tt=tt*multinomial(jn+1,oka(0:jn))
      s=s+tt
   enddo
   qp(jn)=s
enddo
deallocate(  pk, ok, mk )
end subroutine taypow

!=============================================================================
subroutine taymul(n,qa,qb,qc)
!=============================================================================
! (Multiplication of functions)
! Multiply Taylor expansions qa and qb to get qc
!=============================================================================
integer,                intent(IN ):: n
real(dp),dimension(0:n),intent(IN ):: qa,qb
real(dp),dimension(0:n),intent(OUT):: qc
!-----------------------------------------------------------------------------
integer                            :: i,j,k
real(dp)                           :: s
!=============================================================================
do k=0,n
   s=0
   do i=0,k
      j=k-i
      s=s+qa(i)*qb(j)
   enddo
   qc(k)=s
enddo
end subroutine taymul
   
!=============================================================================
subroutine tayval(n,q,z,qz)
!=============================================================================
! (Evaluation of function)
! Evaluate the degree-n Taylor series, q, for Q(z)at a given z.
!=============================================================================
integer,                intent(IN ):: n
real(dp),dimension(0:n),intent(IN ):: q
real(dp),               intent(IN ):: z
real(dp),               intent(OUT):: qz
!-----------------------------------------------------------------------------
integer                            :: i
real(dp)                           :: zp
!=============================================================================
qz=q(0); zp=1; do i=1,n; zp=zp*z; qz=qz+q(i)*zp; enddo
end subroutine tayval

!=============================================================================
subroutine taycom(n,f,g,h)
!=============================================================================
! (Composition of functions)
! Derive the Taylor expansion h that represents the function H(z)=F(G(z))
! when given the expansions f and g of F(z) and G(z).
!=============================================================================
integer,                intent(IN ):: n
real(dp),dimension(0:n),intent(IN ):: f,g
real(dp),dimension(0:n),intent(OUT):: h
!-----------------------------------------------------------------------------
integer                            :: pow
real(dp),dimension(0:n)            :: gp,ggp
!=============================================================================
h=0
h(0)=f(0)
gp=g
do pow=1,n
   h=h+f(pow)*gp
   call taymul(n,g,gp,ggp); gp=ggp
enddo
end subroutine taycom

!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
! Some routines for elementary operations with rational coefficients:

!=============================================================================
subroutine mulrat(rata,ratb,ratc)!                                   [mulrat]
!=============================================================================
integer,dimension(2),intent(IN ):: rata,ratb
integer,dimension(2),intent(OUT):: ratc
!-----------------------------------------------------------------------------
ratc(1)=rata(1)*ratb(1)
ratc(2)=rata(2)*ratb(2)
call euclid(ratc)
end subroutine mulrat
!=============================================================================
subroutine mulrat8(rata,ratb,ratc)!                                   [mulrat]
!=============================================================================
integer(dpi),dimension(2),intent(IN ):: rata,ratb
integer(dpi),dimension(2),intent(OUT):: ratc
!-----------------------------------------------------------------------------
ratc(1)=rata(1)*ratb(1)
ratc(2)=rata(2)*ratb(2)
call euclid(ratc)
end subroutine mulrat8

!=============================================================================
subroutine invrat(rata,ratb)!                                        [invmat]
!=============================================================================
integer,dimension(2),intent(IN ):: rata
integer,dimension(2),intent(OUT):: ratb
!-----------------------------------------------------------------------------
ratb(1)=rata(2)
ratb(2)=rata(1)
call euclid(ratb)
end subroutine invrat
!=============================================================================
subroutine invrat8(rata,ratb)!                                        [invmat]
!=============================================================================
integer(dpi),dimension(2),intent(IN ):: rata
integer(dpi),dimension(2),intent(OUT):: ratb
!-----------------------------------------------------------------------------
ratb(1)=rata(2)
ratb(2)=rata(1)
call euclid(ratb)
end subroutine invrat8

!=============================================================================
subroutine divrat(rata,ratb,ratc)!                                   [divrat]
!=============================================================================
integer,dimension(2),intent(IN ):: rata,ratb
integer,dimension(2),intent(OUT):: ratc
!-----------------------------------------------------------------------------
ratc(1)=rata(1)*ratb(2)
ratc(2)=rata(2)*ratb(1)
call euclid(ratc)
end subroutine divrat
!=============================================================================
subroutine divrat8(rata,ratb,ratc)!                                   [divrat]
!=============================================================================
integer(dpi),dimension(2),intent(IN ):: rata,ratb
integer(dpi),dimension(2),intent(OUT):: ratc
!-----------------------------------------------------------------------------
ratc(1)=rata(1)*ratb(2)
ratc(2)=rata(2)*ratb(1)
call euclid(ratc)
end subroutine divrat8

!=============================================================================
subroutine addrat(rata,ratb,ratc)!                                   [addrat]
!=============================================================================
integer,dimension(2),intent(IN ):: rata,ratb
integer,dimension(2),intent(OUT):: ratc
!-----------------------------------------------------------------------------
ratc(1)=rata(1)*ratb(2) + ratb(1)*rata(2)
ratc(2)=rata(2)*ratb(2)
call euclid(ratc)
end subroutine addrat
!=============================================================================
subroutine addrat8(rata,ratb,ratc)!                                   [addrat]
!=============================================================================
integer(dpi),dimension(2),intent(IN ):: rata,ratb
integer(dpi),dimension(2),intent(OUT):: ratc
!-----------------------------------------------------------------------------
ratc(1)=rata(1)*ratb(2) + ratb(1)*rata(2)
ratc(2)=rata(2)*ratb(2)
call euclid(ratc)
end subroutine addrat8

!=============================================================================
subroutine subrat(rata,ratb,ratc)!                                   [subrat]
!=============================================================================
integer,dimension(2),intent(IN ):: rata,ratb
integer,dimension(2),intent(OUT):: ratc
!-----------------------------------------------------------------------------
ratc(1)=rata(1)*ratb(2) - ratb(1)*rata(2)
ratc(2)=rata(2)*ratb(2)
call euclid(ratc)
end subroutine subrat
!=============================================================================
subroutine subrat8(rata,ratb,ratc)!                                   [subrat]
!=============================================================================
integer(dpi),dimension(2),intent(IN ):: rata,ratb
integer(dpi),dimension(2),intent(OUT):: ratc
!-----------------------------------------------------------------------------
ratc(1)=rata(1)*ratb(2) - ratb(1)*rata(2)
ratc(2)=rata(2)*ratb(2)
call euclid(ratc)
end subroutine subrat8

!=============================================================================
subroutine euclid(rat)!                                              [euclid]
!=============================================================================
integer,dimension(2),intent(INOUT):: rat
integer                           :: kdum,a,b,tt
!=============================================================================
if(rat(2)<0)rat=-rat
if(rat(1)==0)then
   rat(2)=1
   return
endif
b=rat(2)
a=mod(abs(rat(1)),b)

do kdum=1,huge(1)
   if(a==0)exit
   tt=mod(b,a)
   b=a
   a=tt
enddo
rat=rat/b
end subroutine euclid
!=============================================================================
subroutine euclid8(rat)!                                              [euclid]
!=============================================================================
integer(dpi),dimension(2),intent(INOUT):: rat
integer(dpi)                           :: a,b,tt
integer                                :: kdum
!=============================================================================
if(rat(2)<0)rat=-rat
if(rat(1)==0)then
   rat(2)=1
   return
endif
b=rat(2)
a=mod(abs(rat(1)),b)

do kdum=1,huge(1)
   if(a==0)exit
   tt=mod(b,a)
   b=a
   a=tt
enddo
rat=rat/b
end subroutine euclid8

end module pcomb