Characteristic polynomial of central Hyperplane arrangement returns wrong result?

[jplab]

A central hyperplane arrangement must have an even number of regions by central symmetry... yet the one below gets 31 regions(!).

R.<y> = QQ[]
v1 = AA.polynomial_root(AA.common_polynomial(y^2 - 3), RIF(RR(1.7320508075688772), RR(1.7320508075688774)))
v2 = QQbar.polynomial_root(AA.common_polynomial(y^4 - y^2 + 1), CIF(RIF(RR(0.8660254037844386), RR(0.86602540378443871)), RIF(-RR(0.50000000000000011), -RR(0.49999999999999994))))
my_vectors = (vector(AA, [-v1, -1, 1]), vector(AA, [0, 2, 1]), vector(AA,[v1, -1, 1]), vector(AA, [1, 0, 0]), vector(AA, [1/2, AA(-1/2*v2^3 + v2),0]), vector(AA, [-1/2, AA(-1/2*v2^3 + v2), 0]))


sage: H = HyperplaneArrangements(AA,names='xyz')
sage: x,y,z = H.gens()
sage: A = H(backend="normaliz")
sage: for v in my_vectors:
....:     a,b,c = v
....:     A = A.add_hyperplane(a*x + b*y + c*z)
sage: A
Arrangement of 6 hyperplanes of dimension 3 and rank 3
sage: A.n_regions()
31
sage: A.is_central()
True

Here is another failure in characteristic polynomial:

sage: tau = AA((1+sqrt(5))/2)                                                                                                                                                                                      
sage: ncn = [[2*tau+1,2*tau,tau],[2*tau+2,2*tau+1,tau+1],[1,1,1],[tau+1,tau+1,tau],[2*tau,2*tau,tau],[tau+1,tau+1,1],[1,1,0],[0,1,0],[1,0,0],[tau+1,tau,tau]]                                                      
sage: H = HyperplaneArrangements(AA,names='xyz')                                                                                                                                                                   
sage: x,y,z = H.gens()                                                                                                                                                                                             
sage: A = H()                                                                                                                                                                                                      
sage: for v in ncn: 
....:     a,b,c = v 
....:     A = A.add_hyperplane(a*x + b*y + c*z) 
....:                                                                                                                                                                                                              
sage: A.n_regions()
Traceback (most recent call last):
...
ValueError: arrangement cannot simultaneously have h and -h as hyperplane

#30078 fixes this and we add another doctest here.

CC: @kliem @LaisRast @sagetrac-nailuj

Component: geometry

Keywords: hyperplane arrangements, regions

Author: Jonathan Kliem

Branch/Commit: d74929d

Reviewer: Travis Scrimshaw

Issue created by migration from https://trac.sagemath.org/ticket/30749

[fchapoton]

Description changed:

--- 
+++ 
@@ -1,23 +1,18 @@
 A central hyperplane arrangement must have an even number of regions by central symmetry... yet the one below gets 31 regions(!).
 
 ```
-sage: R.<y> = QQ[]
-....: v1 = AA.polynomial_root(AA.common_polynomial(y^2 - 3), RIF(RR(1.7320508075
-....: 688772), RR(1.7320508075688774)))
-....: v2 = QQbar.polynomial_root(AA.common_polynomial(y^4 - y^2 + 1), CIF(RIF(RR
-....: (0.8660254037844386), RR(0.86602540378443871)), RIF(-RR(0.5000000000000001
-....: 1), -RR(0.49999999999999994))))
-....: my_vectors = (vector(AA, [-v1, -1, 1]), vector(AA, [0, 2, 1]), vector(AA, 
-....: [v1, -1, 1]), vector(AA, [1, 0, 0]), vector(AA, [1/2, AA(-1/2*v2^3 + v2), 
-....: 0]), vector(AA, [-1/2, AA(-1/2*v2^3 + v2), 0]))
-....: 
+R.<y> = QQ[]
+v1 = AA.polynomial_root(AA.common_polynomial(y^2 - 3), RIF(RR(1.7320508075688772), RR(1.7320508075688774)))
+v2 = QQbar.polynomial_root(AA.common_polynomial(y^4 - y^2 + 1), CIF(RIF(RR(0.8660254037844386), RR(0.86602540378443871)), RIF(-RR(0.50000000000000011), -RR(0.49999999999999994))))
+my_vectors = (vector(AA, [-v1, -1, 1]), vector(AA, [0, 2, 1]), vector(AA,[v1, -1, 1]), vector(AA, [1, 0, 0]), vector(AA, [1/2, AA(-1/2*v2^3 + v2),0]), vector(AA, [-1/2, AA(-1/2*v2^3 + v2), 0]))
+
+
 sage: H = HyperplaneArrangements(AA,names='xyz')
 sage: x,y,z = H.gens()
 sage: A = H(backend="normaliz")
 sage: for v in my_vectors:
 ....:     a,b,c = v
 ....:     A = A.add_hyperplane(a*x + b*y + c*z)
-....:     
 sage: A
 Arrangement of 6 hyperplanes of dimension 3 and rank 3
 sage: A.n_regions()

[fchapoton]

comment:2

The char poly gives the same answer

sage: p = A.characteristic_polynomial() : p
x^3 - 6*x^2 + 13*x - 11
sage: p(-1)                                                                     
-31

[jplab]

comment:3

Replying to @fchapoton:

The char poly gives the same answer

sage: p = A.characteristic_polynomial() : p
x^3 - 6*x^2 + 13*x - 11
sage: p(-1)                                                                     
-31

Yes, this is how the number of regions is computed up to sign. The hyperplane arrangement seems to be sane: computing .regions() computes the expected (number of) regions.

I asked the number of regions to double-check something and then got this scary answer...

[sagetrac-nailuj]

comment:5

I encountered this problem in some example a long while ago, before I was contributing to Sage myself. A simple workaround is to have n_regions output len(self.regions()), which was good enough for me. Of course, this is much slower if one doesn’t want to compute the regions anyway. I don’t understand the characteristic polynomial well enough to figure out the error in there.

[jplab]

Description changed:

--- 
+++ 
@@ -20,3 +20,21 @@
 sage: A.is_central()
 True
 ```
+
+Here is another failure in characteristic polynomial:
+
+```
+sage: tau = AA((1+sqrt(5))/2)                                                                                                                                                                                      
+sage: ncn = [[2*tau+1,2*tau,tau],[2*tau+2,2*tau+1,tau+1],[1,1,1],[tau+1,tau+1,tau],[2*tau,2*tau,tau],[tau+1,tau+1,1],[1,1,0],[0,1,0],[1,0,0],[tau+1,tau,tau]]                                                      
+sage: H = HyperplaneArrangements(AA,names='xyz')                                                                                                                                                                   
+sage: x,y,z = H.gens()                                                                                                                                                                                             
+sage: A = H()                                                                                                                                                                                                      
+sage: for v in ncn: 
+....:     a,b,c = v 
+....:     A = A.add_hyperplane(a*x + b*y + c*z) 
+....:                                                                                                                                                                                                              
+sage: A.n_regions()
+Traceback (most recent call last):
+...
+ValueError: arrangement cannot simultaneously have h and -h as hyperplane
+```

[mkoeppe]

comment:7

Setting new milestone based on a cursory review of ticket status, priority, and last modification date.

[sagetrac-nailuj]

comment:10

Replying to @sagetrac-nailuj:

I don’t understand the characteristic polynomial well enough to figure out the error in there.

The characteristic polynomial is obtained through recursive deletions and contractions of hyperplane arrangements. During this process, it may occur that some hyperplane h is included in an arrangement multiple times with different scalings of the defining linear expression which are not properly detected as duplicates. This is due to a defect in the method hyperplane.primitive() addressed in #30078.

[mkoeppe]

Dependencies: #30078

[mkoeppe]

comment:12

Is this defect solved by #30078?

[tscrim]

comment:13

It seems like it is. Using the first example, I get

sage: A.n_regions()
24
sage: len(A.regions())
24
sage: A.characteristic_polynomial()(-1)
-24

So it is now consistent. We probably just want to add a doctest for this ticket.