Opened 4 years ago

Closed 3 years ago

# Numerical glitch causes error in computation of fundamental group of curves

Reported by: Owned by: mmarco major sage-8.8 algebraic geometry braids, numerical, QQbar, 32bit tscrim, slelievre, tmonteil Miguel Marco Travis Scrimshaw, Samuel Lelièvre, Frédéric Chapoton, Thierry Monteil N/A b91818d b91818d769ea9bfd81b37d46a50d8dcc6d04a3eb

In Sage 8.3 we get the following (wrong) answer:

```sage: var('t')
sage: wp = QQ[t](t^2 + t + 1).roots(QQbar)[0][0]
sage: Kw.<wp> = NumberField(wp.minpoly(), embedding=wp)
sage: R.<x,y> = Kw[]
sage: z = 1
sage: f = y * (y + (-wp - 1)) * x * (x - 1) * (x - y) * (x + (-wp - 1)*y - 1) * (x + (-wp - 1)*y + (wp))
sage: from sage.schemes.curves import zariski_vankampen as zvk
sage: zvk.fundamental_group(f)
Finitely presented group < x0, x1, x2, x3, x5 | x2*x1*x2^-1*x1^-1, x1*x0*x1^-1*x0^-1, x5*x2*x5^-1*x2^-1, x0^-1*x3*x0*x3^-1, x0^-1*x2*x0*x2^-1, x1*x3*x1^-1*x3^-1, x5^-1*x0*x5*x0^-1, x3^-1*x5*x3*x5^-1, x1*x5*x1^-1*x5^-1 >
```

It must be wrong because the abelianized must have rank 7, so it must have at least 7 generators.

Note that Sage 8.1 used to give right answers.

Anyways, I think I did find the problem. When computing the braids corresponding to each segment, the endpoints are represented by floating point numbers, which may cause some errors. Using algebraic numbers for this seems to solve the issue.

### comment:1 Changed 4 years ago by mmarco

• Branch set to u/mmarco/numerical_glitch_causes_error_in_computation_of_fundamental_group_of_curves

### comment:2 Changed 4 years ago by mmarco

• Branch u/mmarco/numerical_glitch_causes_error_in_computation_of_fundamental_group_of_curves deleted
• Milestone changed from sage-8.4 to sage-8.5

### comment:3 Changed 4 years ago by mmarco

• Status changed from new to needs_review

### comment:4 Changed 4 years ago by slelievre

• Branch set to u/mmarco/numerical_glitch_causes_error_in_computation_of_fundamental_group_of_curves
• Commit set to 7b55efda8aa69d6cfa2d3eaf5493c84e031bdc2b
• Description modified (diff)

Trying to see if Miguel Marco's branch is still there.

New commits:

 ​d0ee394 `Use algebraic numbers to pin the beginning and end of braids` ​cb4c624 `Added docytest` ​7b55efd `Merge branch 'fixzvk' into t/26503/numerical_glitch_causes_error_in_computation_of_fundamental_group_of_curves`

### comment:5 follow-up: ↓ 7 Changed 4 years ago by tscrim

So this change will likely slow things down, but it will make things more numerically stable (as evidenced by this example).

One thing that should be changed:

```-    TESTS::
+    TESTS:

-        # Check that # 26503 is fixed
+    Check that :trac:`26503` is fixed::
```

### comment:6 Changed 4 years ago by slelievre

Pushing the suggestion by Travis Scrimshaw a bit further, I would change

```    TESTS::

# Check that # 26503 is fixed
sage: var('t')
sage: wp = QQ[t](t^2+t+1).roots(QQbar)[0][0]
sage: Kw.<wp>=NumberField(wp.minpoly(), embedding=wp)
sage: R.<x,y> = Kw[]
sage: z = 1
sage: f = y * (y + (-wp - 1)) * x * (x - 1) * (x - y) * (x + (-wp - 1)*y - 1) * (x + (-wp - 1)*y + (wp))
sage: from sage.schemes.curves import zariski_vankampen as zvk
sage: g = f.subs({x:x+y})
sage: g = g.subs({x:x+y})
sage: disc = zvk.discrim(g)
sage: segs = zvk.segments(disc)
sage: segs[16]
(0.577350269189626*I, 0.500000000000000 + 0.288675134594813*I)
sage: zvk.braid_in_segment(g, *segs[16])  # optional - sirocco
s0*s1*s3*s5*s1^-1*s0^-1*s3^-1
```

to the following, if you agree it means the same (please double check):

```    TESTS:

Check that :trac:`26503` is fixed::

sage: wp = QQ['t']([1, 1, 1]).roots(QQbar)[0][0]
sage: Kw.<wp> = NumberField(wp.minpoly(), embedding=wp)
sage: R.<x, y> = Kw[]
sage: z = -wp - 1
sage: f = y*(y + z)*x*(x - 1)*(x - y)*(x + z*y - 1)*(x + z*y + wp)
sage: from sage.schemes.curves import zariski_vankampen as zvk
sage: g = f.subs({x: x + 2*y})
sage: disc = zvk.discrim(g)
sage: segs = zvk.segments(disc)
sage: segs[16]
(0.577350269189626*I, 0.500000000000000 + 0.288675134594813*I)
sage: zvk.braid_in_segment(g, *segs[16])  # optional - sirocco
s0*s1*s3*s5*s1^-1*s0^-1*s3^-1
```

This combines some pep8 fixes, getting rid of a symbolic variable, using `z` to clarify `f`, and combining two `.subs` into a single one.

### comment:7 in reply to: ↑ 5 Changed 4 years ago by mmarco

So this change will likely slow things down, but it will make things more numerically stable (as evidenced by this example).

Yes, but on the other hand, the other minor change I made that computes radical of the polynomial before computing the discriminant would make things faster, so overall the change could actually be faster.

Samuel: I think your change is better. I will push a commit with it tomorrow.

### comment:8 Changed 4 years ago by git

• Commit changed from 7b55efda8aa69d6cfa2d3eaf5493c84e031bdc2b to 5ba2b90026614dad76caf4bb3406942529e2a9dd

Branch pushed to git repo; I updated commit sha1. New commits:

 ​5ba2b90 `Improve doctests with reviewers suggestions`

### comment:9 Changed 4 years ago by slelievre

A blank line is needed between `::` and the code block.

```        sage: g = f.subs({x:x+y})
sage: g = g.subs({x:x+y})
```

into

```        sage: g = f.subs({x: x + 2*y})
```
Last edited 4 years ago by slelievre (previous) (diff)

### comment:10 Changed 4 years ago by git

• Commit changed from 5ba2b90026614dad76caf4bb3406942529e2a9dd to 6e966770c36f215128d9362ee6e3b48908bf1ac4

Branch pushed to git repo; I updated commit sha1. New commits:

 ​6e96677 `Minor fix in the doctest`

### comment:11 Changed 4 years ago by chapoton

• Branch changed from u/mmarco/numerical_glitch_causes_error_in_computation_of_fundamental_group_of_curves to public/ticket/26503
• Commit changed from 6e966770c36f215128d9362ee6e3b48908bf1ac4 to 67f1188713d9f2ba52f649301073beb262008e08
• Reviewers set to Travis Scrimshaw, Samuel Lelièvre, Frédéric Chapoton
• Status changed from needs_review to positive_review

I have made some trivial changes, and I allow myself to turn this to positive.

New commits:

 ​fd8864e `Merge branch 'u/mmarco/numerical_glitch_causes_error_in_computation_of_fundamental_group_of_curves' in 8.5.b5` ​67f1188 `trac 26503 some minor details`

### comment:12 Changed 4 years ago by vbraun

• Status changed from positive_review to needs_work

On 32-bit:

```Doctesting 3781 files using 2 threads.
**********************************************************************
File "src/sage/schemes/curves/zariski_vankampen.py", line 333, in sage.schemes.curves.zariski_vankampen.braid_in_segment
Failed example:
segs[16]
Expected:
(0.577350269189626*I, 0.500000000000000 + 0.288675134594813*I)
Got:
(3.82181265801140e-17 + 0.577350269189626*I,
-0.500000000000000 + 0.288675134594813*I)
**********************************************************************
1 of  15 in sage.schemes.curves.zariski_vankampen.braid_in_segment
[35 tests, 1 failure, 0.91 s]
----------------------------------------------------------------------
sage -t --long src/sage/schemes/curves/zariski_vankampen.py  # 1 doctest failed
----------------------------------------------------------------------
```

### comment:13 Changed 3 years ago by git

• Commit changed from 67f1188713d9f2ba52f649301073beb262008e08 to 1269f8c4d1c93a255c45d492ff7e717111f13518

Branch pushed to git repo; I updated commit sha1. New commits:

 ​1269f8c `#26503 Add abs tol to one doctest`

### comment:14 Changed 3 years ago by slelievre

• Milestone changed from sage-8.5 to sage-8.7
• Status changed from needs_work to needs_review

Added `# abs tol 1e-16`.

### comment:15 Changed 3 years ago by chapoton

• Status changed from needs_review to positive_review

ok

### comment:16 Changed 3 years ago by vbraun

• Status changed from positive_review to needs_work

On 32-bit:

```File "src/sage/schemes/curves/zariski_vankampen.py", line 333, in sage.schemes.curves.zariski_vankampen.braid_in_segment
Failed example:
segs[16]  # abs tol 1e-16
Expected:
(0.577350269189626*I, 0.500000000000000 + 0.288675134594813*I)
Got:
(3.82181265801140e-17 + 0.577350269189626*I,
-0.500000000000000 + 0.288675134594813*I)
**********************************************************************
1 of  15 in sage.schemes.curves.zariski_vankampen.braid_in_segment
[35 tests, 1 failure, 0.79 s]
```

### comment:17 Changed 3 years ago by chapoton

There is a sign problem, not only a matter of tolerance..

### comment:18 follow-up: ↓ 19 Changed 3 years ago by slelievre

Seems the order of the elements in `disc` can vary.

Maybe exactifying the elements of disc and sorting them would fix that:

```        sage: wp = QQ['t']([1, 1, 1]).roots(QQbar)[0][0]
sage: Kw.<wp> = NumberField(wp.minpoly(), embedding=wp)
sage: R.<x, y> = Kw[]
sage: z = -wp - 1
sage: f = y*(y + z)*x*(x - 1)*(x - y)*(x + z*y - 1)*(x + z*y + wp)
sage: from sage.schemes.curves import zariski_vankampen as zvk
sage: g = f.subs({x: x + 2*y})
sage: disc = zvk.discrim(g)
sage: for a in disc: a.exactify()
sage: disc = sorted(disc)
sage: segs = zvk.segments(disc)
sage: segs[16]  # abs tol 1e-16
(0.577350269189626*I, 0.500000000000000 + 0.288675134594813*I)
```

### comment:19 in reply to: ↑ 18 Changed 3 years ago by mmarco

Seems the order of the elements in `disc` can vary.

Maybe exactifying the elements of disc and sorting them would fix that:

```        sage: wp = QQ['t']([1, 1, 1]).roots(QQbar)[0][0]
sage: Kw.<wp> = NumberField(wp.minpoly(), embedding=wp)
sage: R.<x, y> = Kw[]
sage: z = -wp - 1
sage: f = y*(y + z)*x*(x - 1)*(x - y)*(x + z*y - 1)*(x + z*y + wp)
sage: from sage.schemes.curves import zariski_vankampen as zvk
sage: g = f.subs({x: x + 2*y})
sage: disc = zvk.discrim(g)
sage: for a in disc: a.exactify()
sage: disc = sorted(disc)
sage: segs = zvk.segments(disc)
sage: segs[16]  # abs tol 1e-16
(0.577350269189626*I, 0.500000000000000 + 0.288675134594813*I)
```

Yes, that is what I think too. The failing doctest relies on the order of the roots of the discriminant being fixed, which apparently might not be depending on the architecture (as the failing doctest shows). Samuel's proposal makes more sense.

Last edited 3 years ago by slelievre (previous) (diff)

### comment:20 Changed 3 years ago by git

• Commit changed from 1269f8c4d1c93a255c45d492ff7e717111f13518 to 98c921d1de31190676c339133c57b91fedd4cb74

Branch pushed to git repo; I updated commit sha1. New commits:

 ​98c921d `Sort to make braid_in_segment doctest pass on 32-bit`

### comment:21 Changed 3 years ago by slelievre

• Status changed from needs_work to needs_review

Pushed a commit to add this exactifying and sorting. Can someone test on 32-bit?

### comment:22 Changed 3 years ago by slelievre

• Keywords braids numerical QQbar 32bit added

### comment:23 Changed 3 years ago by tscrim

Patchbot seems to have trouble with it, plus one of them is getting an error with `sirocco`:

```File "src/sage/schemes/curves/zariski_vankampen.py", line 337, in sage.schemes.curves.zariski_vankampen.braid_in_segment
Failed example:
zvk.braid_in_segment(g, *segs[16])  # optional - sirocco
Expected:
s0*s1*s3*s5*s1^-1*s0^-1*s3^-1
Got:
s5^-1*s3*s0*s1*s0*s3^-1*s4^-1*s3^-1*s2^-1
```

### comment:24 Changed 3 years ago by tmonteil

Let me run a 32-bit patchbot on it.

### comment:25 Changed 3 years ago by tmonteil

• Status changed from needs_review to needs_work

My 32-bit patchbot is unhappy :(

### comment:26 Changed 3 years ago by slelievre

Order seems to depend on the platform, so we can't rely on `segs[16]` and we need a better doctest.

### comment:27 Changed 3 years ago by mmarco

Can you check if the discrepancy on the order appears on `disc` or `segs` ?

Anyways, maybe the right thing to do about this doctest is to be explicit and give the specific segment to test?

### comment:28 Changed 3 years ago by tscrim

What precisely is that part of the doctest suppose to be testing?

### comment:29 Changed 3 years ago by mmarco

It is testing the computation of the braid along the segment. The segments are created from the Voronoi diagram of the points in disc. We compute them all, and pick a specific one (where the bug shows up) and test it. The problem seems to be that the orderof the computed segments seems to deppend on the platform, so we are picking a different one on 32 bits, and get different results.

### comment:30 Changed 3 years ago by tscrim

So if the bug is showing up in that specific segment, perhaps the thing to do is isolate that segment, maybe by checking if that segment is in `segs` and then passing it to `braid_in_segments`?

### comment:31 Changed 3 years ago by tscrim

Ah, wait, that won't work outright because of the numerical noise.

### comment:32 Changed 3 years ago by tscrim

Perhaps instead of that specific segment, we check

```s0*s1*s3*s5*s1^-1*s0^-1*s3^-1 in [zvk.braid_in_segment(g,s) for s in segs]
```

I guess we could then get the index of that braid in the corresponding list and then check the segment itself (if you also want to show that output).

### comment:33 Changed 3 years ago by git

• Commit changed from 98c921d1de31190676c339133c57b91fedd4cb74 to 5012a7ab8e44b38aa55a9266f89d10a7f6773f53

Branch pushed to git repo; I updated commit sha1. New commits:

 ​a348797 `Use algebraic numbers to pin the beginning and end of braids` ​b3cd564 `Added docytest` ​446a5fb `Improve doctests with reviewers suggestions` ​666941e `Minor fix in the doctest` ​581cf35 `trac 26503 some minor details` ​2f081a6 `#26503 Add abs tol to one doctest` ​5c74d7a `Sort to make braid_in_segment doctest pass on 32-bit` ​1a1af79 `Explicit path in doctest for braid_in_segment` ​5012a7a `Merge branch 'public/ticket/26503' of git://trac.sagemath.org/sage into t/26503/public/ticket/26503`

### comment:34 Changed 3 years ago by mmarco

I changed the doctest by giving explicitely the endpoints of the path. It still detects the bug.

I didn't test it in a 32 bits system though.

### comment:35 Changed 3 years ago by mmarco

• Status changed from needs_work to needs_review

### comment:36 Changed 3 years ago by git

• Commit changed from 5012a7ab8e44b38aa55a9266f89d10a7f6773f53 to bfa164b0e61abd54a3fb87fc366ce0b1ff470169

Branch pushed to git repo; I updated commit sha1. New commits:

 ​bfa164b `Use left normal form of braid in doctest`

### comment:37 Changed 3 years ago by tscrim

The changes are good with me, but you need to mark 2 doctests with `# optional - sirocco` (see pathcbot). Thierry, can you run your 32-bit patchbot on this?

Last edited 3 years ago by tscrim (previous) (diff)

### comment:38 Changed 3 years ago by git

• Commit changed from bfa164b0e61abd54a3fb87fc366ce0b1ff470169 to b91818d769ea9bfd81b37d46a50d8dcc6d04a3eb

Branch pushed to git repo; I updated commit sha1. New commits:

 ​0a95947 `Merge branch 'public/ticket/26503' in 8.7.b7` ​b91818d `adding the optional - sirocco tags`

### comment:39 Changed 3 years ago by chapoton

I have added the 2 missing sirocco tags

### comment:40 Changed 3 years ago by tscrim

Sorry, this fell off my radar by the time I got to a place where I could have made those modifications. Thank you Frédéric for the fix and the reminder.

Can someone with a 32-bit machine retest this for completeness?

### comment:42 follow-up: ↓ 43 Changed 3 years ago by tmonteil

OK. I have to reboostrap one patchbot, so do not expect some result before this evening.

### comment:43 in reply to: ↑ 42 Changed 3 years ago by tscrim

OK. I have to reboostrap one patchbot, so do not expect some result before this evening.

No problem. Thank you for doing that.

### comment:44 Changed 3 years ago by tmonteil

• Reviewers changed from Travis Scrimshaw, Samuel Lelièvre, Frédéric Chapoton to Travis Scrimshaw, Samuel Lelièvre, Frédéric Chapoton, Thierry Monteil

This looks good now on 32bit, see https://patchbot.sagemath.org/log/26503/debian/9.7/i686/4.9.0-7-686/tmonteil-debian-stretch-32/2019-03-21%2013:03:07.

There is an error for the file `src/sage/tests/books/computational-mathematics-with-sagemath/graphique_doctest.py` but it also exists for ticket 0, see https://patchbot.sagemath.org/ticket/0/ so it is unrelated to this ticket.

### comment:45 Changed 3 years ago by tmonteil

The error is actually tracked at #27250.

### comment:46 Changed 3 years ago by tscrim

• Status changed from needs_review to positive_review

Thank you. I am now setting this to positive review.

### comment:47 Changed 3 years ago by slelievre

There's been a lot of trial and error on this, do we want to keep all these commits? We could squash to a single commit instead -- I'm fine either way.

### comment:48 Changed 3 years ago by embray

• Milestone changed from sage-8.7 to sage-8.8

Ticket retargeted after milestone closed (if you don't believe this ticket is appropriate for the Sage 8.8 release please retarget manually)

### comment:49 Changed 3 years ago by vbraun

• Branch changed from public/ticket/26503 to b91818d769ea9bfd81b37d46a50d8dcc6d04a3eb
• Resolution set to fixed
• Status changed from positive_review to closed
Note: See TracTickets for help on using tickets.