.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "api/gallery/numerical_analysis/floating_point/plot_01_kahan_summation.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_numerical_analysis_floating_point_plot_01_kahan_summation.py>`
        to download the full example code.

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

.. _sphx_glr_api_gallery_numerical_analysis_floating_point_plot_01_kahan_summation.py:


Kahan's compensated summation
================================

Adding :math:`n` floating-point numbers one at a time can accumulate
error proportional to :math:`n`. William Kahan's 1965 algorithm carries
the rounding error of each addition forward in a correction term, so
the error stays at a few units in the last place however many terms
are added. This script sums :math:`0.1` repeatedly and compares the
naive loop, Kahan's algorithm, and NumPy's pairwise summation against
the exactly rounded result from :func:`math.fsum`.

.. GENERATED FROM PYTHON SOURCE LINES 15-30

.. code-block:: Python

    import math

    import matplotlib.pyplot as plt
    import numpy as np

    from mathematicskit.numerical_analysis import kahan_sum


    def naive_sum(values):
        s = 0.0
        for v in values:
            s += v
        return s









.. GENERATED FROM PYTHON SOURCE LINES 31-33

Relative error against the number of terms
----------------------------------------------

.. GENERATED FROM PYTHON SOURCE LINES 33-57

.. code-block:: Python


    sizes = np.unique(np.logspace(1, 6, 16).astype(int))
    eps = np.finfo(float).eps
    errors = {"naive loop": [], "numpy.sum (pairwise)": [], "Kahan": []}
    for n in sizes:
        values = [0.1] * int(n)
        exact = math.fsum(values)
        for name, total in (("naive loop", naive_sum(values)), ("numpy.sum (pairwise)", float(np.sum(values))), ("Kahan", kahan_sum(values))):
            errors[name].append(max(abs(total - exact) / exact, eps / 10))

    fig, ax = plt.subplots(figsize=(7, 4.5))
    for name, marker in zip(errors, "os^"):
        ax.loglog(sizes, errors[name], marker + "-", label=name)
    ax.loglog(sizes, sizes * eps, "--", color="gray", label=r"$n\,u$")
    ax.set_xlabel("number of terms $n$")
    ax.set_ylabel("relative error")
    ax.set_title("Summing 0.1 repeatedly")
    ax.legend(fontsize=8)
    fig.tight_layout()

    for name in errors:
        print(f"{name:22s} n = {sizes[-1]}: relative error {errors[name][-1]:.1e}")

    plt.show()



.. image-sg:: /api/gallery/numerical_analysis/floating_point/images/sphx_glr_plot_01_kahan_summation_001.png
   :alt: Summing 0.1 repeatedly
   :srcset: /api/gallery/numerical_analysis/floating_point/images/sphx_glr_plot_01_kahan_summation_001.png
   :class: sphx-glr-single-img


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

 .. code-block:: none

    naive loop             n = 1000000: relative error 1.3e-11
    numpy.sum (pairwise)   n = 1000000: relative error 2.9e-16
    Kahan                  n = 1000000: relative error 2.2e-17





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

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


.. _sphx_glr_download_api_gallery_numerical_analysis_floating_point_plot_01_kahan_summation.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_kahan_summation.ipynb <plot_01_kahan_summation.ipynb>`

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

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

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

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


.. only:: html

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

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