.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "api/gallery/number_theory/primality/plot_01_miller_rabin.py"
.. LINE NUMBERS ARE GIVEN BELOW.

.. only:: html

    .. note::
        :class: sphx-glr-download-link-note

        :ref:`Go to the end <sphx_glr_download_api_gallery_number_theory_primality_plot_01_miller_rabin.py>`
        to download the full example code.

.. rst-class:: sphx-glr-example-title

.. _sphx_glr_api_gallery_number_theory_primality_plot_01_miller_rabin.py:


The Miller-Rabin probabilistic primality test
===================================================

Rabin's 1980 randomization of Miller's test declares :math:`n`
"probably prime" after :math:`k` random witnesses; for a composite
:math:`n` each witness is fooled with probability at most :math:`1/4`,
so the error shrinks like :math:`4^{-k}`. This script shows the test
exposing a Carmichael number that fools Fermat's test, measures the
error rate against the number of rounds for two hard composites, and
certifies a 39-digit prime far beyond trial division's reach.

.. GENERATED FROM PYTHON SOURCE LINES 15-20

.. code-block:: Python

    import matplotlib.pyplot as plt
    import numpy as np

    from mathematicskit.number_theory import fast_mod_pow, is_prime_miller_rabin, prime_factorization








.. GENERATED FROM PYTHON SOURCE LINES 21-23

A Carmichael number fools Fermat but not Miller-Rabin
-----------------------------------------------------

.. GENERATED FROM PYTHON SOURCE LINES 23-29

.. code-block:: Python


    n = 3215031751
    print(f"{n} = {' * '.join(map(str, prime_factorization(n)))}")
    print("Fermat test a^(n-1) = 1 for a = 2, 3, 5, 7:", all(fast_mod_pow(a, n - 1, n) == 1 for a in (2, 3, 5, 7)))
    print("Miller-Rabin (40 random witnesses) says prime:", is_prime_miller_rabin(n))





.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    3215031751 = 151 * 751 * 28351
    Fermat test a^(n-1) = 1 for a = 2, 3, 5, 7: True
    Miller-Rabin (40 random witnesses) says prime: False




.. GENERATED FROM PYTHON SOURCE LINES 30-34

The error probability shrinks like 4^-k
-----------------------------------------------------
Run the test with ``k`` rounds on many random seeds and record how often
a composite slips through.

.. GENERATED FROM PYTHON SOURCE LINES 34-50

.. code-block:: Python


    trials = 4000
    ks = np.arange(1, 7)
    fig, ax = plt.subplots()
    for composite in (3215031751, 9080191):
        rates = [np.mean([is_prime_miller_rabin(composite, k=int(k), seed=s) for s in range(trials)]) for k in ks]
        print(f"n = {composite}: false 'prime' rate by rounds = {np.round(rates, 4).tolist()}")
        seen = np.array(rates) > 0  # rates below 1/trials are unresolved
        ax.semilogy(ks[seen], np.array(rates)[seen], "o-", label=f"n = {composite}")
    ax.semilogy(ks, 4.0 ** (-ks), "k--", label="Rabin's bound 4^-k")
    ax.set_xlabel("rounds k (random witnesses)")
    ax.set_ylabel("probability a composite passes")
    ax.set_title("Miller-Rabin error rate vs. number of witnesses")
    ax.legend()
    plt.show()




.. image-sg:: /api/gallery/number_theory/primality/images/sphx_glr_plot_01_miller_rabin_001.png
   :alt: Miller-Rabin error rate vs. number of witnesses
   :srcset: /api/gallery/number_theory/primality/images/sphx_glr_plot_01_miller_rabin_001.png
   :class: sphx-glr-single-img


.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    n = 3215031751: false 'prime' rate by rounds = [0.2458, 0.0595, 0.013, 0.0038, 0.0, 0.0]
    n = 9080191: false 'prime' rate by rounds = [0.2502, 0.0635, 0.0145, 0.0032, 0.001, 0.0002]




.. GENERATED FROM PYTHON SOURCE LINES 51-55

Certifying a large prime
-----------------------------------------------------
Trial division would need about :math:`10^{19}` divisions for
:math:`2^{127} - 1`; Miller-Rabin needs 40 modular exponentiations.

.. GENERATED FROM PYTHON SOURCE LINES 55-58

.. code-block:: Python


    print("2^127 - 1 is prime:", is_prime_miller_rabin(2**127 - 1))
    print("2^127 + 1 is prime:", is_prime_miller_rabin(2**127 + 1))




.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    2^127 - 1 is prime: True
    2^127 + 1 is prime: False





.. rst-class:: sphx-glr-timing

   **Total running time of the script:** (0 minutes 0.516 seconds)


.. _sphx_glr_download_api_gallery_number_theory_primality_plot_01_miller_rabin.py:

.. only:: html

  .. container:: sphx-glr-footer sphx-glr-footer-example

    .. container:: sphx-glr-download sphx-glr-download-jupyter

      :download:`Download Jupyter notebook: plot_01_miller_rabin.ipynb <plot_01_miller_rabin.ipynb>`

    .. container:: sphx-glr-download sphx-glr-download-python

      :download:`Download Python source code: plot_01_miller_rabin.py <plot_01_miller_rabin.py>`

    .. container:: sphx-glr-download sphx-glr-download-zip

      :download:`Download zipped: plot_01_miller_rabin.zip <plot_01_miller_rabin.zip>`


.. only:: html

 .. rst-class:: sphx-glr-signature

    `Gallery generated by Sphinx-Gallery <https://sphinx-gallery.github.io>`_