module polymod
	implicit none
	
	complex, parameter	:: cmplxone=(1.0,0.0), cmplxzero=(0.0,0.0)
	
	contains

	function polymult(p1,p2) result(p)
		real		:: p1(:),p2(:),p(size(p1)+size(p2)-1)
		integer	:: i
		
		p=0
		do i=1,size(p1)
			p(i:i-1+size(p2))=p(i:i-1+size(p2))+p1(i)*p2
		enddo
	end function
	
	function cpolymult(p1,p2) result(p)
		complex	:: p1(:),p2(:),p(size(p1)+size(p2)-1)
		integer	:: i
		
		p=0
		do i=1,size(p1)
			p(i:i-1+size(p2))=p(i:i-1+size(p2))+p1(i)*p2
		enddo
	end function

	function polysqr(p1) result(p)
		real		:: p1(:),p(2*size(p1)-1)
		p=polymult(p1,p1)
	end function

	function cpolysqr(p1) result(p)
		complex	:: p1(:),p(2*size(p1)-1)
		p=cpolymult(p1,p1)
	end function

	function polyderiv(p) result (dp)
		real		:: p(:),dp(size(p)-1)
		integer	:: i
		
		do i=1,size(p)-1
			dp(i)=i*p(i+1)
		enddo
	end function
	
	function cpolyderiv(p) result (dp)
		complex	:: p(:),dp(size(p)-1)
		integer	:: i
		
		do i=1,size(p)-1
			dp(i)=i*p(i+1)
		enddo
	end function

	function polyminus(p) result(p2)
		real		:: p(:),p2(size(p))
		integer	:: i
		do i=1,size(p)
			p2(i)=merge(p(i),-p(i),mod(i,2)==1)
		enddo
	end function
	
	function polyevsqrt(p) result(p2)
		real		:: p(:),p2(size(p)/2+1)
		integer	:: i
		p2(1)=p(1)
		do i=2,size(p2)
			p2(i)=p(2*i-1)
		enddo
	end function

	function polyevsqr(p) result(p2)
		real		:: p(:),p2(2*size(p)-1)
		integer	:: i
		p2(1)=p(1)
		do i=2,size(p2)
			p2(i)=merge(p((i+1)/2),0.0,mod(i,2)==1)
		enddo
	end function

	function polyevalcmplx(p,z) result(val)
		real		:: p(:)
		complex	:: z,val, r
		integer	:: i
		val=0; r=1
		do i=1,size(p)
			val=val+r*p(i)
			r=r*z
		enddo
	end function
	
	function cpolyevalcmplx(p,z) result(val)
		complex	:: p(:)
		complex	:: z,val, r
		integer	:: i
		val=0; r=1
		do i=1,size(p)
			val=val+r*p(i)
			r=r*z
		enddo
	end function

	function realpolyfromroots2(r) result(p)
		complex	:: r(:),pc(size(r)+1)
		real		:: p(size(r)+1)
		integer	:: i
		
		pc(1:2)=[-r(1),cmplxone]
		do i=2,size(r)
			pc(1:i+1)=[-r(i)*pc(1:i),cmplxzero]+[cmplxzero,pc(1:i)]
		enddo
		p=real(pc)
	end function
		
	function realpolyfromroots(r) result(p)
		complex	:: r(:),pc(size(r)+1)
		real		:: p(size(r)+1)
		integer	:: i
		
		pc(1:2)=[-r(1),cmplxone]
		do i=2,size(r)
			pc(1:i+1)=cpolymult(pc(1:i),[-r(i),cmplxone])
		enddo
		p=real(pc)
	end function

	function cpolyfromroots(r) result(pc)
		complex	:: r(:),pc(size(r)+1)
		integer	:: i
		
		pc(1:2)=[-r(1),cmplxone]
		do i=2,size(r)
			pc(1:i+1)=[-r(i)*pc(1:i),cmplxzero]+[cmplxzero,pc(1:i)]
		enddo
	end function

end module
		
		
		
		
		