# A. Galligo and M. van Hoeij # # Version 0.02, Jan 2007. # Weaknesses in current version: # # The settings of "Digits" inside the program probably needs some fine # tuning, in particular for higher degree inputs. # # Termination is not guaranteed. # # Easy to add features (like separate error bound for each coefficient) # have not yet been implemented. # # Inputs in C[x,y] not yet implemented, only R[x,y] implemented so far. # # TrivialFactoring at the moment only works for 2 factors # Terminology: An Npoly is a polynomial in R[x,y] or in C[x,y] that # has precisely N singularities, all of which are ordinary double points # Assumption: F is close to an Npoly that has the same degree as F. # # The output of this program should be an Npoly that, among all Npolys # of the same degree as F, has a distance from F that is essentially minimal. # If F is too far away from an Npoly then this program may fail. # # Here distance is measured as follows: # distance(F,G) = sqrt( sum of squares of abs values of coeffs of F-G ). # # Note that it is not difficult to modify this program to allow a more # general distance function, of the type: distance(F,G) # = sqrt( sum of w_{i,j} * abs(coeff(F-G, x^i y^j))^2 ) # # where the weights w_{i,j} may be any positive real number but may also # be infinite (indicating that the coefficient of that x^i y^j is exact). # # If the error in the input F is of size O(epsilon) then make sure # that the number Digits is well above log10( 1/epsilon ). # # This first program works over R, the second over C. # NearestNpoly := proc(F, N, x, y) local degY, degX, Fr, Fx, Fy, Ry, v, V, absF, var, cf, xy, eqV, newf, i, j, N2, flag; if nargs=2 then return procname(F, N, op(indets(F, name))) elif hastype(F, nonreal) then return NearestNpolyC(F, N, x, y) fi; degY := degree(F,y); for j from 0 to degY do degX[j] := max(-1, degree(coeff(F,y,j),x)) od; Digits := Digits + 5; Fr := primpart(collect(F*10^(Digits+10), {x,y}, round), {x,y}); Fx,Fy := diff(Fr,x), diff(Fr,y); Ry := evalf(resultant(Fx,Fy,y)); v := [fsolve(Ry, x, complex)]; v := RemoveConjugates(v); V := [seq([i, Y_value(Fx,Fy,x,y,i)], i=v)]; V := [seq( [abs(eval(F, {x=i[1],y=i[2]})),i], i=V)]; # V := sort(V, proc(i,j) evalb(i[1]= N then flag := v-N; if flag=1 then WARNING("found one too many points") fi; N2 := i; break fi od; # If F is small on the first N solutions of Fx=Fy=0 and not small # on the N+1'th solution, then we can be reasonably confident that # we have correctly located (at least within some reasonable error # bound) the N singularities. # # However, if the N+1'th is not much larger than the N'th, then # we're not confident that this program will work, and we generate # a warning message: (note: N'th solution is N2'th element of V) v := V[N2+1, 1]/V[N2, 1]; if v < 10 then WARNING(cat("low level of confidence: ", convert(evalf(v,2),string))) fi; if nargs=5 and args[5]=`return V` then return V[1..N2] fi; cf := add(add(c[i,j]*x^i*y^j, i=0..degX[j]), j=0..degY): var := [seq(seq(c[i,j], i=0..degX[j]), j=0..degY)]: xy := [seq(seq(x^i*y^j, i=0..degX[j]), j=0..degY)]: eqV := {seq(eval(cf, {x=V[i,2,1], y=V[i,2,2]}), i=1..N2)}: eqV := [seq([seq(coeff(i,j),j=var)], i=eqV)]: # Make the equations real: eqV := [seq(map(Re,i),i=eqV), seq(map(Im,i),i=select(hastype,eqV,nonreal))]; V := LinearAlgebra:-NullSpace( Matrix(eqV) ); if nops(V) <> add(degX[j]+1, j=0..degY) - (N+flag) then error "unexpected dimension" fi; newf := ProjectOnV(CoeffVector(F,x,y,degX,degY), V, "GS"); newf := LinearAlgebra:-DotProduct(newf, Vector(xy)); i := Distance(F, newf, x, y); lprint(i); if i > 1/10^iquo(Digits-2,2) then if assigned(_EnvSmallStep) then # This option doesn't seem to help. newf := _EnvSmallStep * F + (1-_EnvSmallStep) * newf fi; # Restore value of Digits: Digits := Digits - 5; procname(newf, N, x, y) else newf fi end: RemoveConjugates := proc(v) local S,i; S := NULL; for i in v do if not member(conjugate(i),[S]) then S := S, i fi od; S := [S]; if add(`if`(type(i,nonreal),2,1),i=S) <> nops(v) then error "failed to remove conjugates" fi; S end: CoeffVector := proc(F,x,y,degX,degY) local i,j; if nargs=4 then return procname(F,x,y,table([seq(i=degX-i,i=0..degX)]), degX) fi; Vector([seq(seq(coeff(coeff(F,x,i),y,j), i=0..degX[j]), j=0..degY)]) end: ProjectOnV := proc(v, V) local i; if nargs=3 then # Use a 3'rd argument if V is already GramSchmidt normalized eval(add(LinearAlgebra:-DotProduct(v,i)*i, i=V)) else procname(v, LinearAlgebra:-GramSchmidt(V,normalized), "GS") fi end: Distance := proc(F, G, x, y) local DEG, h; if nargs=2 then procname(F, G, op(indets([F,G], name))) elif nargs=4 then h := F-G; DEG := degree(h, {x,y}); sqrt(add(add(coeff(coeff(h,y,j),x,i)^2,i=0..DEG-j),j=0..DEG)) else error "wrong number of arguments" fi end: Y_value := proc(Fx, Fy, x, y, X_value) local v; v := eval([Fx,Fy], x=X_value); while has(v[2],y) do v := [v[2],rem(v[1],v[2],y)] od; if degree(v[1],y)<>1 then error fi; solve(v[1],y) end: NearestNpolyC := proc(F, x, y, N) error "not yet implemented" end: # Assumption: just two factors, that intersect transversally, and no other # singularities, giving N = n1*n2 double points, where n1, n2 are the # degrees of the factors (so F has degree n1+n2). # # The input polynomial F should be accurate up to the setting of Digits (if it # is not, then first use the NearestNpoly program to increase its accuracy). # # The method used is a trivial one: just interpolate the singularities. # It also assumes that degree(F,x)=degree(F,y)=degree(F,{x,y}). TrivialFactoring := proc(F, N, x, y) local DEG,n1,n2,n,V,cf,var,eqV,xy; if nargs=2 then return procname(F, N, op(indets(F, name))) elif hastype(F, nonreal) then return TrivialFactoringC(F, N, x, y) fi; DEG := degree(F, {x,y}); n1, n2 := solve(x^2 - DEG*x + N, x); if not type([n1,n2],list(posint)) then error "input inconsistent with assumptions" fi; n := min(n1,n2); V := NearestNpoly(F, N, x, y, `return V`); Digits := Digits + 3; cf := add(add(c[i,j]*x^i*y^j, i=0..n-j), j=0..n): var := [seq(seq(c[i,j], i=0..n-j), j=0..n)]: xy := [seq(seq(x^i*y^j, i=0..n-j), j=0..n)]: eqV := {seq(eval(cf, {x=i[2,1], y=i[2,2]}), i=V)}: eqV := [seq([seq(coeff(i,j),j=var)], i=eqV)]: # Make the equations real: eqV := [seq(map(Re,i),i=eqV), seq(map(Im,i),i=select(hastype,eqV,nonreal))]; # V := LinearAlgebra:-NullSpace( Matrix(eqV) ); # if nops(V) <> `if`(n1=n2,2,1) then # error "unexpected dimension" # fi; V := LinearAlgebra:-SingularValues( Matrix(eqV),'output'=['Vt'] ); V := {seq(Vector[column](V[i, 1..-1]), i = -`if`(n1=n2,2,1) .. -1)}; V := seq(LinearAlgebra:-DotProduct(i, Vector(xy)), i=V); if n1=n2 then FindCombination(F, x, y, V) else Digits := Digits + 10; n := expand(quo(F,V,x)); Digits := Digits - 10; evalf(V) * evalf(n) fi end: FindCombination := proc(F, x, y, f1, f2) local sc,a,b,g,cg,S; sc := [op(expand(F))]; # list of terms of F sc := [seq([eval(i,{x=1,y=1}),degree(i,x),degree(i,y)], i = sc)]; sc := sort(sc, (i,j) -> abs(i[1]) > abs(i[2]) ); # Now sc = terms of F sorted w.r.t. size g := expand( (f1+b*f2)*(f2+c*f1) ); cg := coeff(coeff(g,x,sc[1,2]),y,sc[1,3]); g := expand(sc[1,1]*g - cg*F ); S := {seq(coeff(coeff(g,x,i[2]),y,i[3]), i=sc[2..3])}; S := solve(S, {b,c}); # S contains 2 solutions, pick one: g := abs(b)+abs(c); if eval(g,S[1])