flatsurf / sage-flatsurf

Flat surfaces in Sage
https://flatsurf.github.io/sage-flatsurf/
GNU General Public License v2.0
10 stars 10 forks source link

Spurious doctest error in rotation_matrix_error #151

Closed saraedum closed 2 years ago

saraedum commented 2 years ago

In some setups we get a

File "flatsurf/geometry/matrix_2x2.py", line 208, in flatsurf.geometry.matrix_2x2.rotation_matrix_angle
Failed example:
    for _ in range(100):
        r = QQ.random_element(x=0,y=500)
        r -= r.floor()
        m = rot_matrix(r.numerator(), r.denominator())
        assert rotation_matrix_angle(m) == r
Exception raised:
    Traceback (most recent call last):
      File "/usr/share/miniconda/envs/test/lib/python3.7/site-packages/sage/doctest/forker.py", line 694, in _run
        self.compile_and_execute(example, compiler, test.globs)
      File "/usr/share/miniconda/envs/test/lib/python3.7/site-packages/sage/doctest/forker.py", line 1088, in compile_and_execute
        exec(compiled, globs)
      File "<doctest flatsurf.geometry.matrix_2x2.rotation_matrix_angle[5]>", line 5, in <module>
        assert rotation_matrix_angle(m) == r
      File "/home/runner/work/sage-flatsurf/sage-flatsurf/flatsurf/geometry/matrix_2x2.py", line 224, in rotation_matrix_angle
        raise RuntimeError
    RuntimeError
saraedum commented 2 years ago

Currently, the implementation reads

    e0,e1 = r.change_ring(CDF).eigenvalues()
    m0 = (e0.log() / 2 / CDF.pi()).imag()
    m1 = (e1.log() / 2 / CDF.pi()).imag()
    r0 = RR(m0).nearby_rational(max_denominator=10000)
    r1 = RR(m1).nearby_rational(max_denominator=10000)
    if r0 != -r1:
        raise RuntimeError
    r0 = r0.abs()
    if r[0][1] > 0:
        return QQ.one() - r0
    else:
        return r0

    if check:
        e = r.change_ring(AA).eigenvalues()[0]
        if e.minpoly() != ZZ['x'].cyclotomic_polynomial()(r.denominator()):
            raise RuntimeError
        z = QQbar.zeta(r.denominator())
        if z**r.numerator() != e:
            raise RuntimeError

    return r

Curiously, the check part is never executed.

saraedum commented 2 years ago

I think I prefer to use tr(R) = 1 + 2cos α. I.e., compute an arccos, pick the sign and approximate as a rational. Any reason the current implementation is better? @videlec

saraedum commented 2 years ago

Hm…more importantly, this method is not used anywhere.