Ccluster.jl
Ccluster.jl is a Julia wrapper for Ccluster (https://github.com/rimbach/Ccluster.git) that implements a clustering algorithm for univariate polynomials whose coefficients are complex numbers.
Ccluster.jl also provides a clustering function for triangular systems of polynomial equations.
The most recent release (0.3.0) is only com compatible with julia >= 1.6.0. For a 1.0 <= julia version < 1.6.0, see release 0.2.0 which is not intended to be maintained. The Branch compat-julia-v0.6 is compatible with julia 0.6, but is not intended to be maintained.
Brief description
Univariate solvers
Complex root clustering
The main function provided by Ccluster.jl is ccluster. It takes as input a polynomial P in C[z], a square complex box B and a precision eps.
It outputs a set of natural clusters of roots together with the sum of multiplicities of the roots in each cluster. A cluster is a complex disc D containing at least one root, and it is natural when 3D contains the same roots than D. Each root of P in B is in exactly one cluster of the output, and clusters may contain roots of P in 2B.
Input polynomial P can be given as an oracle polynomial, that is a function returning arbitrary precision approximations of coefficients of P. In case where P has rational coefficients, it can be given exactly.
To cluster all the roots of P, a bound on their modulus (e.g. Fujiwara bound) to find an initial box B containing all the roots. This is done in ccluster when it is called with no input box. Notice that input oracle for P must have no zero leading coefficients.
The implemented algorithm is described here: https://dl.acm.org/citation.cfm?id=2930939
Please cite: https://link.springer.com/chapter/10.1007/978-3-319-96418-8_28 if you use it in your research.
Real root isolator
In the case where P has integer or rational coefficients, one may want to find isolating intervals for its real roots.
Ccluster.jl provides a real root isolator called risolate.
The implemented algorithm is described here: https://arxiv.org/abs/2102.10821
Solver for triangular systems
The notion of natural clusters is straightforwardly extended to the multivariate case. Our function tcluster (t for triangular) takes as input a triangular polynomial system P which polynomials have rational coefficients, a vector of square complex boxes B and a precision eps.
It outputs a set of natural clusters of solutions of P together with the sum of multiplicities of the solutions in each cluster. Each solution of P in B is in exactly one cluster of the output, and clusters may contain solutions of P in 2B.
Let z1<z2<...<zn. Input triangular system P: f1=f2=...=fn=0 must satisfy:
- S has a finite number of solutions,
- greatest variable in fi is zi and deg(fi,zi)>0.
Notice that we do not require simple solutions, and we do not require that the leading coefficient of fi has no common factor with the fj's for j<i.
In the special case where the leading coefficient of fi has no common factor with the fj's for j<i (for instance P is a regular chain), tcluster returns all solutions of P when no initial vector of square complex boxes B is given in input.
The implemented algorithm is described here: https://arxiv.org/abs/1806.10164
Installation
Enter the packages manager with
]
then install the last registered version with:
add Ccluster.jl
Alternatively, install the last version with:
add https://github.com/rimbach/Ccluster.jl
Ccluster depends on Nemo that will be automatically installed.
For graphical outputs, install the package CclusterPlot with
add https://github.com/rimbach/CclusterPlot.jl
in the packages manager.
CclusterPlot depends on PyCall and PyPlot, and requires that matplotlib is installed on your system. It is heavy both to install and to load.
Usage: univariate solver
Simple example: clustering the complex roots of a Mignotte-like polynomial
See the file examples/mignotte.jl
using Nemo
R, x = PolynomialRing(QQ, "x")
d=64
a=14
P = x^d - 2*((2^a)*x-1)^2 #mignotte polynomial
using Ccluster
bInit = [fmpq(0,1),fmpq(0,1),fmpq(4,1)] #box centered in 0 + sqrt(-1)*0 with width 4
precision = 53 #get clusters of size 2^-53
Res = ccluster(P, bInit, precision, verbosity="silent");
#verbosity can take value "silent" (default value),
# "brief" (brief report),
# "results" (brief report and clusters are printed as complex balls
# with precision bits mantissa )
For computing all the roots of P, do:
Res = ccluster(P, precision, verbosity="silent");
Res is an array of couples (sum of multiplicity, disc):
63-element Array{Any,1}:
Any[1, Nemo.fmpq[975//1024, 1025//1024, 15//2048]]
⋮
Any[1, Nemo.fmpq[-2995//4096, 4805//4096, 15//8192]]
Any[2, Nemo.fmpq[0, 0, 15//16384]] # the cluster with sum of multiplicity 2
Any[1, Nemo.fmpq[6935//8192, -8955//8192, 15//16384]]
Any[1, Nemo.fmpq[6935//8192, 8955//8192, 15//16384]]
each element of Res being an array which
- second element is a complex disc (defined by the real and imaginary parts of its center and its radius)
- first element is the sum of multiplicities of the roots in the disk.
You can print the clusters with
printClusters(stdout, Res, precision)
If you care about geometry, so do we. If you have installed CclusterPlot.jl, you can plot the clusters with:
using CclusterPlot
plotCcluster(Res)
or
plotCcluster(Res, bInit, focus=false)
The last argument is a flag telling the function wether to focus on clusters (when true) or not (when false). You can also add markers=false as an optional argument to avoid plotting approximations of the roots with markers.
Real roots: isolating the roots of the same Mignotte polynomials
using Nemo
R, x = PolynomialRing(QQ, "x")
d=64
a=14
P = x^d - 2*((2^a)*x-1)^2 #mignotte polynomial
using Ccluster
bInit = [fmpq(0,1),fmpq(1,1)] #interval centered in 0 with width 1: [-1/2,1/2]
eps = 2
Res = risolate(P, bInit, eps); # each solution is isolated in an interval of radius less than 2^-eps
For computing all the real roots of P, do:
Res = risolate(P, eps);
You can print the clusters with
printClusters(stdout, Res, precision)
Other example: clustering the complex roots of a polynomial whose coefficients are roots of polynomials
See the file examples/coeffsBernoulli.jl
Find the 64 roots of the Bernoulli polynomial of degree 64
using Nemo
RR, x = PolynomialRing(Nemo.QQ, "x")
n = 64 #degree
P = zero(RR)
bernoulli_cache(n)
for k = 0:n
global P
coefficient = (Nemo.binomial(n,k))*(bernoulli(n-k))
P = P + coefficient*x^k
end #P is now the Bernoulli polynomial of degree 64
using Ccluster
bInit = [fmpq(0,1),fmpq(0,1),fmpq(100,1)] #box centered in 0 + sqrt(-1)*0 with width 100
precision = 53 #get clusters of size 2^-53
Coeffs = ccluster(P, bInit, precision)
Define an approximation function for the polynomial whose coefficients are the found roots
function getApproximation( dest::Ptr{acb_poly}, preci::Int )
function getApp(prec::Int)::Nemo.acb_poly
eps=fmpq(1,fmpz(2)^prec)
R = Nemo.RealField(prec)
C = Nemo.ComplexField(prec)
CC, y = PolynomialRing(C, "y")
res = zero(CC)
for i=1:n
btemp = [ Coeffs[i][2][1], Coeffs[i][2][2], 2*Coeffs[i][2][3] ]
temp = ccluster(P, btemp, prec)
approx::Nemo.acb = C( Nemo.ball(R(temp[1][2][1]),R(eps)), Nemo.ball(R(temp[1][2][2]),R(eps)))
res = res + approx*y^(i-1)
end
return res
end
precTemp::Int = 2*preci
poly = getApp(precTemp)
while Ccluster.checkAccuracy( poly, preci ) == 0
precTemp = 2*precTemp
poly = getApp(precTemp)
end
Ccluster.ptr_set_acb_poly(dest, poly)
end
Cluster the roots
bInit = [fmpq(0,1),fmpq(0,1),fmpq(100,1)] #box centered in 0 + sqrt(-1)*0 with width 100
precision = 53 #get clusters of size 2^-53
Roots = ccluster(getApproximation, bInit, precision, verbosity="brief")
Output (total time in s on a Intel(R) Core(TM) i7-7600U CPU @ 2.80GHz):
-------------------Ccluster: ----------------------------------------
-------------------Input: ----------------------------------------
|box: cRe: 0 cIm: 0 wid: 100 |
|eps: 1/100 |
|strat: newton tstarOpt predPrec anticip |
-------------------Output: ----------------------------------------
|number of clusters: 63 |
|number of solutions: 63 |
-------------------Stats: ----------------------------------------
|total time: 1.802102 |
---------------------------------------------------------------------
63-element Array{Any,1}:
Any[1, Nemo.fmpq[-3125//32768, 5125//8192, 375//65536]]
⋮
Any[1, Nemo.fmpq[211625//262144, -105125//262144, 375//524288]]
Defining an approximation function
ccluster takes as input a function prototyped as:
function getApproximation( dest::Ptr{Nemo.acb_poly}, p::Int )
Here is an example for a polynomial with complex coefficients (see also the file examples/spiral.jl)
degr=64
function getApproximation( dest::Ptr{acb_poly}, precision::Int )
function getAppSpiral( degree::Int, prec::Int )::Nemo.acb_poly
CC = ComplexField(prec)
R2, y = PolynomialRing(CC, "y")
res = R2(1)
for k=1:degree
modu = fmpq(k,degree)
argu = fmpq(4*k,degree)
root = modu*Nemo.exppii(CC(argu))
res = res * (y-root)
end
return res
end
precTemp::Int = 2*precision
poly = getAppSpiral( degr, precTemp)
while Ccluster.checkAccuracy( poly, precision ) == 0
precTemp = 2*precTemp
poly = getAppSpiral(degr, precTemp)
end
Ccluster.ptr_set_acb_poly(dest, poly)
end
Defining an initial box
The initial box is an array of three Nemo.fmpq defining respectively the real part of the center, the imaginary part of the center and the width of the box.
The following code:
bInit = [fmpq(0,1),fmpq(0,1),fmpq(150,1)]
defines a box centered in 0+i0 with width 150.
The precision
The precision is an integer p. Ccluster computes clusters of size eps=2^-p.
The verbosity flag
The last, optional, argument of ccluster is a verbosity flag. When no verbosity is given, ccluster is silent. Values can be "brief" and "results".
Usage: solver for triangular systems
Simple example: clustering the roots of a Mignotte-like polynomial
See the file examples/triangularSys.jl
using Nemo
using Ccluster
Rx, x = PolynomialRing(QQ, "x") #Ring of polynomials in x with rational coefficients
Syx, y = PolynomialRing(Rx, "y") #Ring of polynomials in y with coefficients that are in Rx
d1=30
c = 10
delta=128
d2=10
twotodelta = fmpq(2,1)^(delta)
f = Rx( x^d1 - (twotodelta*x-1)^(c) )
g = Syx( y^d2 - x^d2 )
precision = 53
bInitx = [fmpq(0,1),fmpq(0,1),fmpq(10,1)^40]
nbSols, clusters, ellapsedTime = tcluster( [f,g], [bInitx], precision, verbosity = "silent" );
Here g is monic and its leading coefficient has no common factor with f. In general, you can check this property using Nemo.gcd:
Nemo.gcd(f,Nemo.lead(g))
Providing this is 1, compute clusters containing all the solutions with:
nbSols, clusters, ellapsedTime = tcluster( [f,g], precision, verbosity = "silent" );
nbSols is the total number of solutions (counted with multiplicity), clusters is an array of clusters of solutions and ellapsedTime is the time spent to solve the system.
Each element in clusters is an array which first element is the sum of multiplicities of the solution in the cluster and second element is a precision-bit approximation of the solutions (of type Nemo::acb).
You can print the clusters with:
printClusters(stdout, nbSols, clusters)