speed up logarithms modulo composite integers

[yyyyx4]

The current logic for IntegerModRing_abstract.log() is as follows:

1. If logarithm_exists=True was passed: Call Pari's znlog with just the elements.
2. Else, if the modulus is prime: Compute the orders of a,b and (if the logarithm exists) call Pari's znlog with the elements and the order of b.
3. Else, use generic discrete_log (i.e., Pohlig–Hellman).

The main reason things are done like this appears to be that znlog is not guaranteed to terminate if the logarithm doesn't exist, and checking if it exists is not straightforward when the ambient group is non-cyclic. Since we're already at it, here's what Pari does in znlog:

1. If the order of the base element is given, factor it and use Pohlig–Hellman.
2. Else, factor the modulus and solve each prime-power part separately using either index calculus (if p is a suitable prime) or Pohlig–Hellman with either BSGS or Pollard ρ for the base cases (depending on size).

We really shouldn't resort to generic discrete_log() right away if the modulus isn't prime: While indeed (ℤ/n)* is usually non-cyclic, we can typically reduce to the cyclic case using CRT, which we can then pass on to Pari after checking that the logarithm exists (by comparing orders). This patch takes care of this.

Summary of the new algorithm:

1. Factor the modulus.
2. Solve the problem modulo each prime power separately:
   1. Modulo powers of two, use Sage's discrete_log.
   2. Modulo odd prime powers, check that the logarithm exists and call Pari's znlog.
3. Use CRT to combine the local solutions.

The reason to treat powers of two differently is that this is the only case where the unit group is not always cyclic, hence we might run into Pari's infinite loop when no solution exists (side note: I'm not sure if this can actually happen, but of course we should rely only on documented behaviour). Since 2-groups are a very easy case for generic discrete_log(), we simply deal with that in Sage instead of risking the infinite loop. In all other cases, unit groups modulo prime powers are cyclic, hence we can safely call Pari after checking the orders.

I know factoring the modulus looks scary, but let me point out that we're already doing that in all cases:

• If logarithm_exists=True is given, the order is factored in Pari first thing.
• Else, we are computing the orders of elements in Sage, which requires first factoring the group order.

The new doctests introduced in the patch are cases where the old algorithm would take a very long time (if it managed to finish at all without getting killed for lack of memory), whereas the new algorithm takes much less than a second. I also tested the old and new code with thousands of random DLPs modulo (1) random integers, and (2) primes, of sizes up to 2^64, and there didn't appear to be any noticeable performance degradation even in the previous best case of prime moduli. For composite integers, the new code is unsurprisingly a lot faster on average.

I found out about the current, suboptimal implementation here: https://web.archive.org/web/20210813155145/https://cryptohack.gitbook.io/cryptobook/abstract-algebra/groups/untitled

CC: @JohnCremona @kwankyu @defeo @edgarcosta @kedlaya

Component: number theory

Keywords: discrete logarithms, finite rings, composite characteristic

Author: Lorenz Panny

Branch/Commit: d866118

Reviewer: Travis Scrimshaw, Edgar Costa

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

[yyyyx4]

Commit: 7559691

[yyyyx4]

Branch: public/faster_logarithms_in_ZZ_mod_n

[yyyyx4]

Author: Lorenz Panny

[tscrim]

comment:2

I am cc-ing some people who might have an opinion on this. (Sorry if you don't care.)

Some quick comments:

I don't think we need the assert b**n == self check to be run. I would comment it out.

I would just copy the error message where appropriate rather than put it in a lambda function. It obfuscates the code a bit for something effectively trivial (having an easy way to change an error message). Plus you can then have more nuanced error messages explaining where the computation failed.

Split these lines:

-n, m = 0, 1
+n = 0
+m = 0

-na, nb = a_red.multiplicative_order(), b_red.multiplicative_order()
+na = a_red.multiplicative_order()
nb = b_red.multiplicative_order()

It makes a difference in Cython code. (Related question: do we know anything about the types of some of these variables? This could make the code even faster.)

Move the import statement from sage.arith.all import crt, lcm outside of the for loop. It only needs to be done once in the function.

It is only the crt that could fail here, right? So

             try:
                 n = crt(n, v, m, nb)
-                m = lcm(m, nb)
             except ValueError:  # contradictory partial solutions
                 raise exc(self.modulus())
             m = lcm(m, nb)

[edgarcosta]

comment:3

Overall it looks like a solid improvement, and I agree with Travis's comments.

Another minor comment, I would consider storing a list with nbs and vs and at the end calling CRT_list which is literally doing the same as you

                n = crt(n, v, m, nb)
                m = lcm(m, nb)

but now the code would be more readable and modular.

[sagetrac-git]

Changed commit from 7559691 to 161082a

[sagetrac-git]

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

db23291	typos; add tests for failure modes
161082a	incorporate feedback from #32375

[yyyyx4]

comment:5

Thanks for the comments! All done (I think), except for this:

Replying to @edgarcosta:

I would consider storing a list with nbs and vs and at the end calling CRT_list which is literally doing the same

— that's what I had at first, but the current version aborts earlier in some cases. (We don't want to continue running an expensive DLP computation when it's already clear that there will be no solution.)

[edgarcosta]

comment:6

Makes sense! I will wait for the tests to pass to take a more careful look into it.

[tscrim]

comment:7

Green patchbot. Now to some technical comments:

-for kw in kwargs.keys():
+for kw in kwargs:

However, while I somewhat understand why you are using kwargs input, you could just set logarithm_exists=None and check if None. I doubt anyone is passing None here and meant to.

Trivial formatting point: We try to follow Python's convention and have error messages starting with a lower case letter.

[edgarcosta]

comment:8

In terms of math, the code looks good! Once Travis is happy, I will give a positive review.

I agree that logarithm_exists=None sounds ideal to me, and less opaque why we still have it, no one needs to wonder what goes in **kwargs, and we can just mark it as deprecated in the documentation.

Optional formatting:

• delete the empty lines at the beginning of if, for, try... blocks
• for p,e in ... -> for p, e in ...
• extra minor, I'm a strong fan of https://docs.python.org/3/tutorial/inputoutput.html#fancier-output-formatting

[tscrim]

comment:9

Replying to @edgarcosta:

I agree that logarithm_exists=None sounds ideal to me, and less opaque why we still have it, no one needs to wonder what goes in **kwargs, and we can just mark it as deprecated in the documentation.

I would just leave it out altogether from the documentation since it is deprecated and the warning/code makes it clear that it is. No need to tell a more casual user about something that is disappearing.

[sagetrac-git]

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

99e216f	more feedback from #32375
e6f06b7	slight documentation tweaks while we're at it

[sagetrac-git]

Changed commit from 161082a to e6f06b7

[yyyyx4]

comment:11

Thanks both of you! Done.

[tscrim]

comment:12

Thank you. I am setting this to a positive review from comment:8.

Edgar, if there is something else you feel needs to be addressed, feel free to revert the positive review.

[tscrim]

Reviewer: Travis Scrimshaw, Edgar Costa

[edgarcosta]

comment:13

Thank you, Travis!

[vbraun]

comment:14

sage -t --long --warn-long 42.3 --random-seed=0 src/sage/functions/log.py
**********************************************************************
File "src/sage/functions/log.py", line 385, in sage.functions.log.log
Failed example:
    log(a,3)
Expected:
    Traceback (most recent call last):
    ...
    ValueError: No discrete log of 8 found to base 3 modulo 13
Got:
    <BLANKLINE>
    Traceback (most recent call last):
      File "/home/release/Sage/local/lib64/python3.9/site-packages/sage/functions/log.py", line 451, in log
        return args[0].log(args[1])
      File "sage/rings/finite_rings/integer_mod.pyx", line 777, in sage.rings.finite_rings.integer_mod.IntegerMod_abstract.log (build/cythonized/sage/rings/finite_rings/integer_mod.c:9308)
        raise ValueError(f"no logarithm of {self} found to base {b} modulo {self.modulus()}" \
    ValueError: no logarithm of 8 found to base 3 modulo 13
    <BLANKLINE>
    During handling of the above exception, another exception occurred:
    <BLANKLINE>
    Traceback (most recent call last):
      File "sage/symbolic/function.pyx", line 573, in sage.symbolic.function.Function.__call__ (build/cythonized/sage/symbolic/function.cpp:6731)
        args = [SR.coerce(a) for a in args]
      File "sage/structure/parent.pyx", line 1177, in sage.structure.parent.Parent.coerce (build/cythonized/sage/structure/parent.c:11028)
        cpdef coerce(self, x):
      File "sage/structure/parent.pyx", line 1209, in sage.structure.parent.Parent.coerce (build/cythonized/sage/structure/parent.c:10974)
        return (<map.Map>mor)._call_(x)
      File "sage/structure/coerce_maps.pyx", line 161, in sage.structure.coerce_maps.DefaultConvertMap_unique._call_ (build/cythonized/sage/structure/coerce_maps.c:4638)
        raise
      File "sage/structure/coerce_maps.pyx", line 156, in sage.structure.coerce_maps.DefaultConvertMap_unique._call_ (build/cythonized/sage/structure/coerce_maps.c:4530)
        return C._element_constructor(x)
      File "sage/symbolic/ring.pyx", line 402, in sage.symbolic.ring.SymbolicRing._element_constructor_ (build/cythonized/sage/symbolic/ring.cpp:6851)
        raise TypeError('positive characteristic not allowed in symbolic computations')
    TypeError: positive characteristic not allowed in symbolic computations
    <BLANKLINE>
    During handling of the above exception, another exception occurred:
    <BLANKLINE>
    Traceback (most recent call last):
      File "/home/release/Sage/local/lib64/python3.9/site-packages/sage/doctest/forker.py", line 704, in _run
        self.compile_and_execute(example, compiler, test.globs)
      File "/home/release/Sage/local/lib64/python3.9/site-packages/sage/doctest/forker.py", line 1098, in compile_and_execute
        exec(compiled, globs)
      File "<doctest sage.functions.log.log[32]>", line 1, in <module>
        log(a,Integer(3))
      File "/home/release/Sage/local/lib64/python3.9/site-packages/sage/functions/log.py", line 455, in log
        return logb(args[0], args[1])
      File "sage/symbolic/function.pyx", line 1179, in sage.symbolic.function.BuiltinFunction.__call__ (build/cythonized/sage/symbolic/function.cpp:12857)
        res = super(BuiltinFunction, self).__call__(
      File "sage/symbolic/function.pyx", line 586, in sage.symbolic.function.Function.__call__ (build/cythonized/sage/symbolic/function.cpp:6893)
        raise TypeError("cannot coerce arguments: %s" % (err))
    TypeError: cannot coerce arguments: positive characteristic not allowed in symbolic computations
**********************************************************************
1 item had failures:
   1 of  50 in sage.functions.log.log
    [330 tests, 1 failure, 1.97 s]
----------------------------------------------------------------------
sage -t --long --warn-long 42.3 --random-seed=0 src/sage/functions/log.py  # 1 doctest failed
----------------------------------------------------------------------

[sagetrac-git]

Changed commit from e6f06b7 to d866118

[sagetrac-git]

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

d866118

the error message has changed

[yyyyx4]

comment:16

This error is a result of this absolutely hilarious piece of code in src/sage/functions/log.py:

    except ValueError as ex:
        if repr(ex)[12:27] == "No discrete log":
            raise

The "fix" is easy, but the underlying issue isn't: Has there been any effort to systematize exception types in Sage so that this could be replaced by something along the lines of except NoSolutionFound:?

[tscrim]

comment:17

That is quite impressive of a failure. The fix LGTM.

I don't think there has been, although it could be good to have a NoSolutionError, perhaps as a subclass of ArithmeticError. This likely would need a discussion on sage-devel.

[vbraun]

Changed branch from public/faster_logarithms_in_ZZ_mod_n to d866118