#Copyright since 2002 by Bin Han, University of Alberta,
#All Rights Reserved 

#The use of all the following software is permitted for
#non-commercial, educational, and research use only. The software
#and/or related materials are provided "as-is" without warranty of any
#kind including any warranties of performance or merchantability or
#fitness for a particular use or purpose or for any purpose
#whatsoever, for the licensed product, however used. In no event shall
#University of Alberta and/or Bin Han be liable for any damages and/or
#costs, including but not limited to incidental or consequential
#damages of any kind, including economic damage or injury to property
#and lost profits, regardless of whether University of Alberta shall
#be advised, have reason to know, or in fact shall know of the
#possibility. User bears all risk relating to quality and performance
#of the software and/or related materials. 

#Any use other than non-commercial, educational, or research, or any
#redistribution in original or modified form requires prior written
#authorization from the copyright holder.
 
#Report errors, mistakes, and bugs to Bin Han at bhan@ualberta.ca.
#Send comments/suggestions to Bin Han at bhan@ualberta.ca
#====================================================================

#MAPLE 7.0 routines for algorithms on constructing 1D pairs
#of dual d-wavelet frames with any given dilation factor d(>=2)
#The algorithms are based on results described in the following paper:
#    Author: Ingrid Daubechies and Bin Han
#    Title:  Pairs of dual wavelet frames from any two
#            refinable functions
#    Submitted to: Constructive Approximation.
#===================================================================
# Explanation of main routines
#
#DwfByCorollary3_5(dilation,apoly,bpoly,primalvm,dualvm,choice)
#
#This routine implements Corollary 3.5 in the above paper
#
#when primalvm<=0 or dualvm<=0, then we take the highest possible orders;
#that is (1+x+...+x^(dilation-1))^primalvm | b(x)
#and (1+x+...+x^(dilation-1))^dualvm | a(x)
#with highest possible primalvm and dualvm.
#when the given primalvm and dualvm are great than the best possible ones,
#we truncate it to be equal to the best possible order
#when choice=0, output all filters in exact format (may be complex-valued)
#when choice=1, output all filters in exact format with real coefficients.
#othewise, output in numerical format with 12 digits and real part of coeffs.
#This routine uses g(x)=1.
#Output all filters format:
#[1,1]=apoly, [1,2..dilation+1]=primal wavelet filters
#[2,1]=bpoly, [2,2..dilation+1]=dual wavelet filters
#[3,1]:=Theta, [3,2..dilation+1]=cpoly vector
#[4,1]=dilation, [4,2]=g(x), [4,3]=h(x)

#DwfByTheorem4_3(dilation,apoly,bpoly,primalvm,dualvm,N,choice)
#
#This routine implements Theorem 4.3 in the above paper
#
#when primalvm<=0 or dualvm<=0, then we take the highest possible orders;
#that is (1+x+...+x^(dilation-1))^primalvm | b(x)
#and (1+x+...+x^(dilation-1))^dualvm | a(x)
#with highest possible primalvm and dualvm.
#when the given primalvm and dualvm are great than the best possible ones,
#we truncate it to be equal to the best possible order
#when choice=0, output all filters in exact number format
#othewise, output in numerical format with 12 digits.
#when choice=0, output all filters in exact format (may be complex-valued)
#when choice=1, output all filters in exact format with real coefficients.
#othewise, output in numerical format with 12 digits and real part of coeffs.
#This routine use g(x)=1.
#Output all filters format:
#[1,1]=apoly, [1,2..dilation+1]=primal wavelet filters
#[2,1]=bpoly, [2,2..dilation+1]=dual wavelet filters
#[3,1]:=Theta, [3,2..dilation+1]=cpoly vector
#[4,1]=dilation, [4,2]=g(x), [4,3]=h(x)

# There are also many other useful routines for outputing
# and check the results. See the documentation for each routine
# in the program.

#Developed by Bin Han on April 2, 2002 (First version).
#Tests have been performed by Bin Han on April 2, 2002.

#Report any bugs and mistakes to me at  bhan@ualberta.ca
#Visit my web page for more detail: http://www.ualberta.ca/~bhan
###################################################################
interface(rtablesize=infinity,quiet=true):
with(LinearAlgebra):

######The following are some output and auxiliary routines======

#Convert a poly in x into a set or an Array
#When choice=0, output=a set of all coefficients of poly
#when choice =1, output=[-2=low, -1=high, 0..(high-low)=coefficients of poly]:
#     poly is supported on [low, high]
#when choice=2, then output=[1..(high-low+3)=[low, high, coefficients ofpoly]
#when choice=3, then output=[1..(high-low+3)=[low, high, realpart of coefficients ofpoly]
#when choice=4, output=[1..(high-low+3)=[low,high, numer coefficients]
#when choice>4 output=[1..(high-low+3)=[low,high,real part of coefficients]
# In this routine, we also try to simplify the coefficients
PolyToSeq:=proc(poly,choice)
  local low,high,k,result:
   high:=degree(poly,x):
   low:=-degree(subs({x=1/x},poly),x):
   if(high=-infinity) then high:=0: low:=0: fi:
   if(choice=0) then
     result:={}:
     for k from low to high do
       result:=result union {simplify(coeff(poly,x,k))}:
     od:
   elif (choice=1) then
     result:=Array((-2)..(high-low)):
     result[-2]:=low:  result[-1]:=high:
     for k from low to high do
       result[k-low]:=simplify(coeff(poly,x,k)):
     od:
   elif(choice<0) then
     result:=Vector[row](high-low+4):
     result[1]:=low: result[2]:=high: low:=1:
     for k from result[1] to result[2] do
         high:=coeff(poly,x,k):
         result[k+4-result[1]]:=high:
         low:=lcm(low, denom(high)):
     end:
     for k from 4 to (result[2]-result[1]+4) do
         result[k]:=result[k]*low:
     end:
     result[3]:=low:
   else
     result:=Vector[row](high-low+3):
     result[1]:=low: result[2]:=high:
     for k from low to high do
      if(choice=2 ) then
         result[k-low+3]:=simplify(coeff(poly,x,k)):
      elif (choice=3) then
         result[k-low+3]:=Re(coeff(poly,x,k)):
      elif (choice=4) then
         result[k-low+3]:=evalf(coeff(poly,x,k),12):
      else
         result[k-low+3]:=Re(evalf(coeff(poly,x,k),12)):
      fi:
     od:
   fi:
   result:=result:
end proc:

#Check whether a sequence is symmetric or antisymmetric
#if so, output its symmetry center: seq(c-k)=seq(k) for all k, that is,
#poly=x^c*subs({x=1/x},poly).
#output[1]=-1 (antisymmetric), 0 (nonsymmetric), 1 (symmetric)
#output[2]:=c (the symmetry center)
IsSym:=proc(poly)
  local k,low,high,seq,total,result:
  seq:=PolyToSeq(poly,1): low:=seq[-2]: high:=seq[-1]:
  result:=Vector(2):
  result[1]:=0: result[2]:=-999999: total:=0:
  if ((seq[0]+seq[high-low])=0) then  result[1]:=-1: fi:
  if ((seq[0]-seq[high-low])=0) then  result[1]:=1:  fi:
  if (result[1]<>0) then
     for k from 0 to trunc((high-low+1)/2) do
        total:=total+abs(seq[k]-result[1]*seq[high-low-k]):
     od:
  fi:
  if (total=0) then
     result[2]:=high+low:
  else
     result[1]:=0:
  fi:
  result:=result:
end proc:

#check whether a sequence is real-valued or not.
#if it is real-valued, output=1, otherwise, 0
IsRealSeq:=proc(poly)
  local k,seq,result:
  seq:=PolyToSeq(poly,1):
  result:=1:
  for k from 0 to (seq[-1]-seq[-2]) do
      if (Im(seq[k])<>0) then
         result:=0:
         break:
      fi:
  od:
  result:=result:
end proc:

#When Bpoly | Apoly, output Apoly/Bpoly in a polynomial form,
#otherwise, output 0 and a warning message
ApolyOverBpoly:=proc(Apoly,Bpoly)
  local ndegA,ndegB,APOLY,BPOLY,result:
  ndegA:=degree(subs({x=1/x},Apoly),x):
  ndegB:=degree(subs({x=1/x},Bpoly),x):
  APOLY:=expand(simplify(x^ndegA*Apoly)):
  BPOLY:=expand(simplify(x^ndegB*Bpoly)):
  if( rem(APOLY,BPOLY,x,'Qpoly')=0 ) then
      result:=Qpoly*x^(ndegB-ndegA):
  else
      result:=0:
      print("Bpoly can not divide Apoly in rotuine ApolyOverBpoly!"):
  fi:
  result:=expand(result):
end proc:

#find the highest order of sum rules satisfied by a sequence
#or find the highest order of vanishing moments satisfied by a sequence
#choice=0, find the order of vanishing moments
#choice<>0, find the order of sum rules
SrOrVm:=proc(dilation,poly,choice)
  local k,absdilation,denompoly,ndegpoly,POLY,result:
  absdilation:=abs(dilation):
  ndegpoly:=degree(subs({x=1/x},poly),x):
  POLY:=expand(x^ndegpoly*poly):
  result:=0: ndegpoly:=degree(POLY,x):
  if(choice<>0) then
     denompoly:=1:
     for k from 1 to (absdilation-1) do
         denompoly:=denompoly+x^k:
     od:
     for k from 1 to ndegpoly do
       if(rem(POLY,denompoly^k,x)=0) then:
         result:=result+1:
       else
         break;
       fi:
    od:
  else
     for k from 0 to ndegpoly do
         if(subs({x=1}, POLY)=0) then:
            result:=result+1:
            POLY:=diff(POLY,x):
         else
            break;
         fi:
       od:
  fi:
  result:=result:
end proc:

#Output all the necessary filters in the construction of DWF
#When apoly and bpoly are real-valued sequence, output
#all the data in real numbers
#when choice=0, exact number form
#when choice<>0, numerical number form with 12 digits:
Output:=proc(allfilters,choice)
   local tmp,k,format,dilation,result:
   if ((IsRealSeq(allfilters[1,1])<>0)and(IsRealSeq(allfilters[2,1])<>0)) then
      if(choice=0) then format:=3: else format:=5: fi:
   else
      if(choice=0) then format:=2: else format:=4: fi:
   fi:
   dilation:=allfilters[4,1]:
   printf("The dilation factor = %d\n", dilation):
   tmp:=PolyToSeq(allfilters[3,1],format):
   printf("The Theta(x) sequence is\n"): print(tmp):
   tmp:=PolyToSeq(allfilters[4,2],format):
   print("The g(x) sequence is"): print(tmp):
   tmp:=PolyToSeq(allfilters[4,3],format):
   print("The h(x) sequence is"): print(tmp):
   for k from 2 to (abs(dilation)+1) do
       tmp:=PolyToSeq(allfilters[3,k],format):
       print("cpoly No.", k-1,  "is"):
       print(tmp):
   od:
   tmp:=PolyToSeq(allfilters[1,1],format):
   print("The primal mask with sum rules", SrOrVm(dilation,allfilters[1,1],1), "is"):
   print(tmp):
   for k from 2 to (abs(dilation)+1) do
       tmp:=PolyToSeq(allfilters[1,k],format):
       print("Primal wavelet filter No.", k-1, "with vanishing moments",SrOrVm(dilation,allfilters[1,k],0) \, "is"):
       print(tmp):
   od:

   tmp:=PolyToSeq(allfilters[2,1],format):
   print("The dual mask with sum rules", SrOrVm(dilation,allfilters[2,1],1), "is"):
   print(tmp):
   for k from 2 to (abs(dilation)+1) do
       tmp:=PolyToSeq(allfilters[2,k],format):
       print("Dual wavelet filter No.", k-1, "with vanishing moments", SrOrVm(dilation,allfilters[2,k],0), "is"):
   print(tmp):
   od:
end proc:

#use the right order of vanishing moments
PossibleVM:=proc(dilation,apoly,bpoly,primalvm,dualvm)
  local bpvm,bdvm,result:
  result:=Vector(2):
  bpvm:=SrOrVm(dilation,bpoly,1):
  bdvm:=SrOrVm(dilation,apoly,1):
  if( (primalvm<=0) or (primalvm>=bpvm) ) then
    result[1]:=bpvm:
    print("Using the best possible order of vanishing moments for primal!",bpvm):
  else
    result[1]:=primalvm:
  fi:
  if( (dualvm<=0) or (dualvm>=bdvm)) then
    result[2]:=bdvm:
    print("Using the best possible order of vanishing moments for dual!",bdvm):
  else
    result[2]:=primalvm:
  fi:
  result:=result:
end proc:

#convert a poly into the real part of the poly
PolyToRealpoly:=proc(poly)
  local low,high,k,result:
   high:=degree(poly,x):
   low:=-degree(subs({x=1/x},poly),x):
   result:=0:
   for k from low to high do
       result:=result+Re(evalc(coeff(poly,x,k)))*x^k:
   od:
   result:=result:
end proc:


#Condition: Apoly(1)=Bpoly(1)=1.
#Find a polynomial Theta with shortest support such that
#Theta(1)=1 and
# (1-x)^vmorder | DiffPoly:=[Theta(x)Apoly(x)-Theta(x^dilation)Bpoly(x)].
#When Apoly and Bpoly are symmetric with symmetry centers cA and cB such that
#     cT:=(cA-cB)/(dilation-1) is an integer,
#then we can find a Theta such that the symmetry center of Theta is cT
#and the symmetry center of DiffPoly is cT+cA=cB+dilation*cT.
FindTheta:=proc(dilation,Apoly,Bpoly,vmorder)
  local k,symA,symB,cA,cB,cT,truncord,halfcT,EQ,VAR,POLY,result:
  result:=0:
  if ( (subs({x=1},Apoly)<>1) and (subs({x=1},Bpoly)<>1)) then
     print("Conditions are not satisfied in routine FindTheta, set Theta=0!"):
     quit:
  fi:
  symA:=IsSym(Apoly):  symB:=IsSym(Bpoly):
  cA:=symA[2]: cB:=symB[2]: cT:=trunc( (cA-cB)/(dilation-1)):
  halfcT:=trunc((cT+1)/2):  truncord:=trunc((vmorder-1)/2): VAR:={}:  
  if( (symA[1]=1) and (symB[1]=1) and (((dilation-1)*cT-cA+cB)=0) ) then
     #find a symmetric Theta with symmetry center cT
     print("Finding a symmetric Theta with symmetry center", cT):
     for k from halfcT to halfcT+truncord do
         result:=result+evaln(cc||(k-halfcT))*(x^k+x^(cT-k)):
         VAR:=VAR union { evaln(cc||(k-halfcT))}:
     od:
   else
     #find a nonsymmetric shortest Theta
     print("Finding a NON-symmetric Theta with symmetry center", cT):
     for k from (halfcT-truncord) to (halfcT-truncord+vmorder-1) do
         result:=result+evaln(cc||(k+truncord-halfcT))*(x^k):
         VAR:=VAR union { evaln(cc||(k+truncord-halfcT))}:
     od:
   fi:
   POLY:=expand(result*Apoly-subs({x=(x^dilation)},result)*Bpoly):
   EQ:={subs({x=1},result)-1}:
   for k from 1 to (vmorder-1) do
       POLY:=diff(POLY,x):
       EQ:=EQ union { subs({x=1},POLY) }:
   od:
   solt:=solve(EQ,VAR); print("have a solution?", solt);
   assign(solt):
   result:=simplify(2*result/2):
   result:=result:
end proc:



#Find the right side of the construction of Dual Wavelet Frames.
#output[1]:=Theta(x)-Theta(x^d)*\olverline{apoly(x)}*bpoly
#output[j+1]:=-Theta(x^d)*\olverline{apoly(exp(-i(\xi+2\pi j/dilation))}*bpoly
#x:=exp(-i\xi)
Findrpoly:=proc(dilation,apoly,bpoly,Theta)
  local APOLY,dTheta,k,result:
  result:=Vector[column](abs(dilation),0):
  dTheta:=subs({x=x^dilation},Theta):
  for k from 0 to (abs(dilation)-1) do
result[k+1]:=expand(-dTheta*subs({x=(exp(I*2*Pi*k/dilation)/x)},apoly)*bpoly):
  od:
  result[1]:=Theta+result[1]:
  result:=result:
end proc:

#Find the polynomial f as given in the paper.
#fpoly[j]:=rpoly/(1-exp(i(\xi+2\pi j/dilation))^primavm\overline(g(x)h(x)}
Findfpoly:=proc(dilation,apoly,bpoly,Theta,gpoly,hpoly,primalvm)
  local numerpoly,denompoly,cgpoly,chpoly,k,result:
  result:=Vector[column](abs(dilation),0):
  chpoly:=subs({x=1/x},hpoly):
  numerpoly:=Findrpoly(dilation,apoly,bpoly,Theta):
  for k from 0 to (abs(dilation)-1) do
denompoly:=expand((1-exp(I*2*Pi*k/dilation)/x)^primalvm*subs({x=(exp(I*2*Pi*k/dilation)/x)},gpoly)*chpoly):
      result[k+1]:=ApolyOverBpoly(numerpoly[k+1],denompoly):
  od:
  result:=result:
end proc:

#given cpoly, a vector of polynomials, generate the matrix cmatrix
#which is defined to be cmatrix[j,k]:=\overline{c^j(\xi+2*\pi*k/dilation)}
GenCmatrix:=proc(dilation,cpoly)
   local j,k,absdilation,result:
   absdilation:=abs(dilation):
   result:=Matrix(absdilation,absdilation,0):
   for j from 1 to absdilation do
      for k from 0 to (absdilation-1) do
          result[k+1,j]:=subs({x=(exp(I*2*Pi*k/dilation)/x)},cpoly[j]):
      od:
   od:
   result:=result:
end proc:

#compute the bpolyvec for wavelet masks for dual part
Findbpolyvec:=proc(dilation,cpoly,fpoly)
  local cmatrix,result:
  cmatrix:=GenCmatrix(dilation,cpoly):
  cmatrix:=Adjoint(cmatrix):
  result:=cmatrix.fpoly:
end proc:

#compute the apolyvect for wavelet masks for primal part 
Findapolyvec:=proc(dilation,cpoly,gpoly,primalvm)
  local k,result:
  result:=Vector[column](abs(dilation)):
  for k from 1 to abs(dilation) do
      result[k]:=(1-x)^primalvm*gpoly*cpoly[k]:
  end:
  result:=result:
end proc:

#Find a desirable polynomial C as givne in Proposition 2.4
#such that (1-e^{-i\xi})^(2N) | c^j(\xi), for j=2,...,dilation
#and cpoly has the require symmetry center.
Findcpoly:=proc(dilation,N,s0)
  local absdilation,m,j,k,result:
  absdilation:=abs(dilation):
  m:=trunc((absdilation+2)/2):
  result:=Vector[column](absdilation,0):
  for j from 0 to m-1 do
      result[j+1]:=1:
      for k from 0 to m-1 do
          if (k<>j) then
              result[j+1]:=result[j+1]*(x+1/x-2*cos(2*Pi*k/dilation)):
          elif((j<>0) and (N>0)) then
              result[j+1]:=result[j+1]*(2-x-1/x)^(N-1):
          fi:
      od:
  od:
  result[1]:=1:
  if ((abs(s0) mod 2)=0) then
     for j from 1 to m do
         result[j]:=x^(trunc(s0/2))*result[j]:
     od:
     for j from 1 to (absdilation-m) do
         result[j+m]:=simplify(result[j+1]*(x-1/x)):
     od:
  elif ( (absdilation mod 2)=0 ) then
     for j from 1 to (m-1) do
         result[j+m-1]:=x^(trunc((s0-1)/2))*(1-x)*result[j+1]:
         result[j]:=x^(trunc((s0-1)/2))*(1+x)*result[j]:
     od:
  else
     for j from 1 to (absdilation-m) do
         result[j+m]:=x^(trunc((s0-1)/2))*(1-x)*result[j+1]:
     od:
     for j from 1 to m do
         result[j]:=x^(trunc((s0-1)/2))*(1+x)*result[j]:
     od:
   fi:
  result:=result:
end proc:


#Find a desirable polynomial C as givne in Lemma ??
#such that (1-e^{-i\xi})^(2N) | c^j(\xi), for j=2,...,dilation
#This program uses a short cpoly comparing with Findcpoly().
Findcpoly2:=proc(dilation,N)
  local absdilation,m,j,k,result:
  absdilation:=abs(dilation):
  m:=trunc((absdilation+2)/2):
  result:=Vector[column](absdilation,0):
  result[1]:=1:
  if ( (absdilation mod 2)=0) then
     if (N>0) then
        result[m]:=(2-x-1/x)^(N-1)
     else
        result[m]:=(2-x-1/x):
     fi:
  fi:
  for j from 1 to absdilation-m do
      result[j+1]:=1:
      for k from j+1 to m do
          result[j+1]:=result[j+1]*(x+1/x-2*cos(2*Pi*k/dilation)):
          if( (N>0)) then
              result[j+1]:=result[j+1]*(2-x-1/x)^(N-1):
          fi:
      od:
  od:
  for j from 1 to (absdilation-m) do
      result[j+m]:=simplify(result[j+1]*(x-1/x)):
  od:
  result:=result:
end proc:

#check the original conditions for pairs of dual wavelet frames
CheckDWF:=proc(dilation,apoly,bpoly,Theta,apolyvec,bpolyvec)
   local k,absdilation,amatrix,rpoly,result:
   amatrix:=GenCmatrix(dilation,apolyvec):
   rpoly:=Findrpoly(dilation,apoly,bpoly,Theta):
   rpoly:=Amatrix.bpolyvec-rpoly:
   result:={}:
   for k from 1 to abs(dilation) do
       result:=result union PolyToSeq(rpoly[k],0):
   od:
   result:=result:
end proc:


#This routine implements Corollary 3.5 in
#  I. Daubechies and B. Han, Pairs of dual wavelet frames from any two
#  refinable functions
#when primalvm<=0 or dualvm<=0, then we take the highest possible orders;
#that is (1+x+...+x^(dilation-1))^primalvm | b(x)
#and (1+x+...+x^(dilation-1))^dualvm | a(x)
#with highest possible primalvm and dualvm.
#when the given primalvm and dualvm are great than the best possible ones,
#we truncate it to be equal to the best possible order
#when choice=0, output all filters in exact number format
#when choice=1, output all filters in exact number format with real coeffs.
#when otherwise, output in numerical format with 12 digits and real part.
#This routine uses g(x)=1.
#Output all filters format as follows:
#[1,1]=apoly, [1,2..dilation+1]=primal wavelet filters
#[2,1]=bpoly, [2,2..dilation+1]=dual wavelet filters
#[3,1]:=Theta, [3,2..dilation+1]=cpoly vector
#[4,1]=dilation, [4,2]=g(x), [4,3]=h(x)

DwfByCorollary3_5:=proc(dilation,apoly,bpoly,primalvm,dualvm,choice)
  local j,k,d,vm,Theta,cpoly,cmatrix,hpoly,fpoly,tmppoly,result:
  d:=abs(dilation):
  result:=Matrix(4,d+1,0):
  vm:=PossibleVM(dilation,apoly,bpoly,primalvm,dualvm):
  cpoly:=Vector[column](d):
  for k from 1 to d do  cpoly[k]:=x^(k-1):  od:
  cmatrix:=GenCmatrix(dilation,cpoly):
  hpoly:=subs({x=1/x}, Determinant(GenCmatrix(dilation,cpoly))):
  result[1,1]:=apoly:  result[2,1]:=bpoly:
  tmppoly:=Findapolyvec(dilation,cpoly,1,vm[1]):
  for k from 1 to d do result[1,k+1]:=tmppoly[k]: od:
  Theta:=FindTheta(dilation,1,subs({x=1/x},apoly)*bpoly,vm[1]+vm[2]):
  fpoly:=Findfpoly(dilation,apoly,bpoly,Theta,1,hpoly,vm[1]):
  tmppoly:=Findbpolyvec(dilation,cpoly,fpoly):
  for k from 1 to d do result[2,k+1]:=tmppoly[k]: result[3,k+1]:=cpoly[k]: od:
  result[4,1]:=dilation:
  result[3,1]:=Theta: result[4,2]:=1: result[4,3]:=hpoly:
  for j from 1 to 3 do
  for k from 1 to (d+1) do
      if(choice=1) then
         result[j,k]:=PolyToRealpoly(result[j,k]):
      elif(choice<>0) then
         result[j,k]:=PolyToRealpoly(evalf(result[j,k], 32)):
      fi:
  od:od:
  result:=result:
end proc:

#This routine implements Theorem 4.3 in
#  I. Daubechies and B. Han, Pairs of dual wavelet frames from any two
#  refinable functions, submitted to Constructive Approximation
#when primalvm<=0 or dualvm<=0, then we take the highest possible orders;
#Noet the order of sum rules means 
#     (1+x+...+x^(dilation-1))^primalvm | b(x)
# and (1+x+...+x^(dilation-1))^dualvm   | a(x)
#with highest possible primalvm and dualvm.
#when the given primalvm and dualvm are greater than the best possible ones,
#we truncate it to be equal to the best possible order
#when choice=0, output all filters in exact number format
#othewise, output in numerical format with 12 digits.
#This routine use g(x)=1.
#Output all filters format:
#[1,1]=apoly, [1,2..dilation+1]=primal wavelet filters
#[2,1]=bpoly, [2,2..dilation+1]=dual wavelet filters
#[3,1]:=Theta, [3,2..dilation+1]=cpoly vector
#[4,1]=dilation, [4,2]=g(x), [4,3]=h(x)
DwfByTheorem4_3:=proc(dilation,apoly,bpoly,primalvm,dualvm,N,choice)
  local j,k,d,vm,SymA,SymB,Theta,Apoly,Bpoly,cpoly,cmatrix,hpoly,fpoly,tmppoly,result:
  d:=abs(dilation):
  result:=Matrix(4,d+1,0):
  vm:=PossibleVM(dilation,apoly,bpoly,primalvm,dualvm):
  SymA:=IsSym(apoly):
  SymB:=IsSym(bpoly):
  #if the symmetry conditions are satisfied, use symmetric cpoly
  if((SymA[1]=1)and(SymB[1]=1)and((abs(SymA[2]-SymB[2]) mod (d-1))=0)) then
       cpoly:=Findcpoly(dilation,N,SymA[2]-vm[1]):
  else
       cpoly:=Findcpoly(dilation,N,0):
  fi:
  result[4,1]:=dilation:
  for k from 1 to d do result[3,k+1]:=cpoly[k]: od:
  hpoly:=subs({x=1/x}, Determinant(GenCmatrix(dilation,cpoly))):
  if(hpoly=0) then
    print("There are something wrong with the program, check it!"):
    quit:
  fi:
  result[4,2]:=1: result[4,3]:=hpoly:
  result[1,1]:=apoly:  result[2,1]:=bpoly:
  tmppoly:=Findapolyvec(dilation,cpoly,1,vm[1]):
  for k from 1 to d do result[1,k+1]:=tmppoly[k]: od:
  if( (d mod 2)=0 ) then
      Apoly:=hpoly*subs({x=x^(1/d)}, expand( x^(trunc(dilation/2))*hpoly)):
   else
      Apoly:=hpoly*subs({x=x^(1/d)}, expand(hpoly)):
   fi:
  Apoly:=Apoly/subs({x=1},Apoly):
  Bpoly:=subs({x=x^dilation},Apoly)*subs({x=1/x},apoly)*bpoly:
  Theta:=FindTheta(dilation,Apoly,Bpoly,vm[1]+vm[2]+2*N):
  Theta:=Theta*Apoly:
  fpoly:=Findfpoly(dilation,apoly,bpoly,Theta,1,hpoly,vm[1]):
  tmppoly:=Findbpolyvec(dilation,cpoly,fpoly):
  for k from 1 to d do result[2,k+1]:=tmppoly[k]: od:
  result[3,1]:=Theta:
  for k from 1 to d do result[3,k+1]:=cpoly[k]: od:
  for j from 1 to 3 do
  for k from 1 to (d+1) do
      if(choice=1) then
         result[j,k]:=PolyToRealpoly(result[j,k]):
      elif(choice<>0) then
         result[j,k]:=PolyToRealpoly(evalf(result[j,k], 32)):
      fi:
  od:od:
  result:=result:
end proc:


#check the original conditions for pairs of dual wavelet frames
CheckDWF1:=proc(dilation,apoly,bpoly,Theta,apolyvec,bpolyvec)
   local k,absdilation,amatrix,rpoly,result:
   amatrix:=GenCmatrix(dilation,apolyvec):
   rpoly:=Findrpoly(dilation,apoly,bpoly,Theta):
   rpoly:=amatrix.bpolyvec-rpoly:
   result:={}:
   for k from 1 to abs(dilation) do
       result:=result union PolyToSeq(rpoly[k],0):
   od:
   result:=result:
end proc:

#the input is in one set
CheckDWF2:=proc(allfilters)
   local k,d,apoly,bpoly,Theta,apolyvec,bpolyvec,result:
   d:=abs(allfilters[4,1]):
   apolyvec:=Vector[column](d):
   bpolyvec:=Vector[column](d):
   apoly:=allfilters[1,1]:
   bpoly:=allfilters[2,1]:
   Theta:=allfilters[3,1]:
   for k from 1 to d do
       apolyvec[k]:=allfilters[1,k+1]:
       bpolyvec[k]:=allfilters[2,k+1]:
   od:
   result:=CheckDWF1(dilation,apoly,bpoly,Theta,apolyvec,bpolyvec):
   result:=result minus {0}:
   if (result={}) then
      print("Filters satisfy the Dual Wavelet Frame Equations!"):
   else
      print("Filters DO NOT satisfy Dual Wavelet Frame Equations!"):
      print("the set of residual is", result):
   fi:
end proc:


OutputInC:=proc(poly,filename)
   local low,high,k:
   high:=degree(poly,x):
   low:=-degree(subs({x=1/x},poly),x):
   writeto(filename):
   printf("%d %d\n", low, high):
   for k from low to high do
       printf("%16.12f\n", evalf(Re(coeff(poly,x,k)), 32)):
   od:
   writeto(terminal):
end proc: