/********************************************************************************/ /* Copyright (C) 2014 Christelle Guichard */ /* 2014 Vincent Delecroix */ /* */ /* */ /* Distributed under the terms of the GNU General Public License (GPL) */ /* http://www.gnu.org/licenses/ */ /********************************************************************************/ /*****************************************************************/ /* Roots localization with respect to the unit disk */ /* ------------------------------------------------ */ /* */ /* The main function of this module is polrootslocalization(p). */ /* It returns a vector [kout,kon,kin] where kout,kon and kin */ /* are the number of roots of the polynomial outside the unit */ /* disk, on the unit disk and inside. */ /* */ /* see also the following related gp functions: */ /* polsturm(p,{ab}) : localization on the real axis */ /* poliscyclo(f) : check if f is cyclotomic */ /* poliscycloprod(f): check if the f is a product of cyclotomic */ /*****************************************************************/ polschur(p)= { my(pstar); pstar = conj(polrecip(p)); return(conj(polcoeff(p,0)) * p - conj(polcoeff(pstar,0)) * pstar); } addhelp(polschur,"polschur(p): The Schur transform of the polynomial p.") polchebyshevtransform(p)= { my(d,sol=0,t0,t1,x=variable(p)); if(polrecip(p) != p, error("p is not reciprocal")); t0 = 1; t1 = x; d = poldegree(p); if(d%2 == 0, sol = polcoeff(p,d/2), sol=0); d = floor((d+1)/2); sol += 2*polcoeff(p,d-1) * t1; for(i=2,d, [t0,t1] = [t1,2*x*t1 - t0]; sol += 2* polcoeff(p,d-i) * t1); sol; } addhelp(polchebyshevtransform,"polchebyshevtransform(p,{v}): given a reciprocal polynomial p"\ " returns the polynomial whose roots are the (a+1/a)/2 where a are the roots of a. The"\ " result has half degree as the input.") polrootslocalization(p)= { my(q,qq,l,n,kin=0,kkon,kon=0,x=variable(p)); if(p == 0, error("zero polynomial in localization.")); if(poldegree(p)==0, return([0,0,0])); /* remove roots 0, 1 and -1 */ while(polcoeff(p,0)==0, p/=x; kin+=1); while(subst(p,x,-1)==0, p/=(x+1); kon+=1); while(subst(p,x,1)==0, p/=(x-1); kon+=1); /* degree 0 is fine */ if(poldegree(p) == 0, return([0,kon,kin])); /* now we consider the case Tp = 0 (<=> p is reciprocal) */ /* we apply Chebychev transform */ /* (the polynomial can not be anti-reciprocal since -1 */ /* is not a root anymore) */ if(polrecip(p) == p, q = polchebyshevtransform(p); n = poldegree(q); if((subst(q,x,-1) == 0) || (subst(q,x,1) == 0), error("hum hum...")); kkon = 0; /* polsturm count roots without multiplicity */ /* we hence need to rerun it several times */ while(poldegree(q) >= 1, qq = gcd(q,q'); kkon += polsturm(q/qq,[-1,1]); q = qq; ); return([n-kkon, kon+2*kkon, kin+n-kkon]); ); /* now consider Schur transform... if |a_0| = |a_n| where p = a_n x^n + ... + a_0 */ /* we can not conclude. In that case, we transform the polynomial into (x-2) p */ /* and it might be enough!? */ n = poldegree(p); if(abs(polcoeff(p,0)) == abs(polcoeff(p,n)), l = polrootslocalization(p*(x-2)); return([l[1]-1, kon+l[2], kin+l[3]]); ); /* Tp(0)>0, same number of zeros inside */ /* Tp(0)<0, same number of zeros outside */ q = polschur(p); l = polrootslocalization(q); if(polcoeff(q,0) > 0, return([n-l[2]-l[3], kon+l[2], kin+l[3]]), return([l[3], kon+l[2], kin+n-l[2]-l[3]])); } addhelp(polrootslocalization,"polrootslocalization(p): return a vector [k_out, k_on, k_in] where k_out, k_on, k_in are respectively\n"\ "the number of roots out of the disk, on the circle and inside the disk. If we were unable to compute these\n"\ "numbers then -1 is returned") /**************/ /* some tests */ /**************/ testlocalization()= { (localization_via_approx(p)= my(a,k,ipos,ineg,izero); a=apply(abs,polroots(p)); for(k=1,length(a),ipos+=a[k]>1.00000001; ineg+=a[k]<0.99999999); return([ipos,length(a)-ipos-ineg,ineg]); ); /* test Chebyshev desymmetrization */ random_pols = [\ 1,\ x,\ x^3-x,\ x^2 - 1,\ x^3 + x - 2,\ x^4 + x^2 - 5,\ x^5 + 3*x^3 - x^2 + 1\ ]; for(i=1,length(random_pols), p=random_pols[i]; n=poldegree(p); q=x^n*subst(p,x,(x+1/x)/2); r=polchebyshevtransform(q); if(p!=r, printf("ERROR Chebyshev transform p=%s q=%s r=%s\n",p,q,r)); ); /* now test localization */ pisot_polynomials = [\ x^2 - 3*x - 2,\ x^2 - 1*x - 1,\ x^2 - 2*x - 1,\ x^2 - 3*x - 1,\ x^2 - 2*x - 2,\ x^2 - 4*x + 1,\ x^3 - x - 1,\ x^3 - x^2 - x - 1,\ x^3 - 3*x^2 + 1,\ x^5 - x^4 - x^3 - x^2 - 1,\ x^5 - 10*x^4 - 15*x^3 - 3*x^2 + 3*x + 1,\ x^6 - 15*x^5 - 20*x^4 + 6*x^3 + 18*x^2 + 8*x + 1,\ x^6 - 30*x^5 - 60*x^4 - 32*x^3 + 3*x^2 + 6*x + 1,\ x^8 - x^7 - x^6 + x^2 - 1\ ]; if(polrootslocalization(x^2*(x-1)^3*(x+1)^4) != [0,7,2],print("ERROR x^2(x-1)^3(x+1)^4")); for(i=1,length(pisot_polynomials), p = pisot_polynomials[i]; if(polrootslocalization(p) != [1,0,poldegree(p)-1], printf("ERROR with Pisot polynomial p=%s\n",p);); q = x^12 * p - polrecip(p); /* a Salem built from Pisot */ if(polrootslocalization(q) != [1,poldegree(q)-2,1], printf("ERROR with Salem from Pisot polynomial q=%s\n",q);); q *= x^2*(x-1)^2*(x+1)^2; if(polrootslocalization(q) != [1,poldegree(q)-4,3], printf("ERROR with Salem from Pisot polynomial q=%s\n",q);); q = x^12 * p + polrecip(p); /* a Salem built from Pisot */ if(polrootslocalization(q) != [1,poldegree(q)-2,1], printf("ERROR with Salem from Pisot polynomial q=%s\n",q);); ); (exampleSalem1(i,n,v='x)= return(v^(n+1) - 2*v^n - v^(n+1-i) + v^(n-i) + v^(i+1) - v^i - 2*v + 1); ); for(i=3,6, for(n=i,10, p = exampleSalem1(i,n); if(polrootslocalization(p) != [1,poldegree(p)-2,1], printf("ERROR with Salem polynomial p=%s\n",p)); ); ); test_polynomials = [x^3-2, x^6 - 4*x^4 + 4*x - 2, x^6 + x - 4, 2*x^6 - x - 5, 4*x^5 - 3*x + 5]; for(i=1, length(test_polynomials), p = test_polynomials[i]; if(localization_via_approx(p) != polrootslocalization(p), printf("ERROR with p=%s\n",p)); p *= (x-1)^2*(x+1)^2; if(localization_via_approx(p) != polrootslocalization(p), printf("ERROR with p=%s\n",p)); ); for(k=2, 30, p = polcyclo(k); if(localization_via_approx(p) != [0,eulerphi(k),0], printf("ERROR cyclo(%d)\n",k)); if(localization_via_approx(p^2) != [0,2*eulerphi(k),0], printf("ERROR cyclo(%d)\n",k)); p = x^2*(x^3 - 2) * (x^3 - 1/4) * polcyclo(k); if(localization_via_approx(p) != [3,eulerphi(k),5], printf("ERROR x^2(x^3-2)(x^3-1/4)cyclo(%d)\n",k)); ); more_examples=[\ [x^6 + 4*x^5 - x^4 + 15*x^3 - x^2 + 4*x + 1, [3,0,3]],\ [x^11 - 3*x^10 - 4*x^9 + 12*x^7 + x^6 - x^5 - 12*x^4 + 4*x^2 + 3*x - 1, [3,5,3]]\ ]; for(i=1,length(more_examples), [p,l]=more_examples[i]; if(polrootslocalization(p) != l, printf("ERROR with p=%s\n",p);); if(polrootslocalization(p^2) != 2*l, printf("ERROR with p=%s\n",p);); if(polrootslocalization(x*p^2) != 2*l+[0,0,1], printf("ERROR with p=%s\n",p);); ); }