R<I>:=ComplexField(30);
LP<Z> := LaurentSeriesRing(R);
T:=(-Z+2-Z^-1)/4;

TR<U>:=PowerSeriesRing(R);
Pol<X>:=PolynomialRing(R);

readi A,"What accuracy should the wavelets have?";

PolarFactorization:=function(A)
	pp:=Truncate( (1-U+O(U^A))^-A );
	pp:=Evaluate(pp,X);
	p:=LP!1;
	for rr in Roots(pp) do
		r:=rr[1];r;
		for ff in Roots(X^2+(4*r-2)*X+1) do
			f:=ff[1];
			if Abs(f) ge 1 then p*:=Z-f; end if;
		end for;
	end for;
	p1:=Evaluate(p,1);
	return LP![ Real(c): c in Eltseq(p/p1) ],pp;
end function;

p,p2:=PolarFactorization(A);p;a:=2^(1-A)*(1+Z)^A*p;a2:=(1-X)^A*p2;
"Scaling sequence",Coefficients(a),"product filter", Coefficients(a2);
	

Decimation:= function(c)
	return LP![ Coefficient(c,2*k): k in [0..Degree(c) div 2+2] ];
end function;

// Power iteration for the values of the scaling function at integer locations 
shape := LP!(Z^A);
for i := 1 to 10 do
	for k := 1 to 10 do
		shape := Decimation(a*shape); shape := shape/Evaluate(shape,1);
	end for;
	shape:=LP![ R!(1+Coefficient(shape,k))-1: k in [0..Degree(shape)] ];
end for;

wave := LP![ (-1)^k*Coefficient(a,2*A-k): k in [1..2*A] ];
scal := shape;

dx := 1; supp := 2*A-1; pow:=1;
for i := 1 to 4 do
	scal := scal*LP!Evaluate(a,Z^pow);
	dx /:=2; supp *:=2; pow*:=2;
end for;

wave := scal*LP!Evaluate(wave,Z^pow);	
scal := scal*LP!Evaluate(a,Z^pow);
dx /:=2; supp *:=2; pow*:=2;

a2;

fp := Open(Sprintf("daub%o-scal.dat",A),"w");
RO:=RealField(12); 
for k := 1 to supp do
	fprintf fp, "%o\t%o\t%o\t%o\t%o\n",RO!(k*dx),
	RO!Coefficient(scal,k),
	RO!Coefficient(wave,k),
	RO!Abs(Evaluate(scal,Exp(I*2*Pi(RO)*dx^2*k))*dx),
	RO!Abs(Evaluate(wave,Exp(I*2*Pi(RO)*dx^2*k))*dx);
end for;
Flush(fp);