{
  "nbformat": 4,
  "nbformat_minor": 5,
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "name": "python",
      "version": "3.12.14"
    },
    "execution": {
      "method": "Python subprocess, exact code cell source",
      "deterministic_replay_passed": true
    },
    "runtime": {
      "packages": [],
      "local_assets": []
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "id": "question",
      "metadata": {},
      "source": [
        "# Why parallel sums disagree: tracing FP32 reduction order\n",
        "Two correct addition orders return zero and one for inputs whose exact sum is two. Enumerating every input order reveals the numerical contract.\n",
        "\nCPU FP32 rounding illustration with exact chosen inputs; not a device race test or GPU benchmark.\n"
      ]
    },
    {
      "cell_type": "code",
      "id": "experiment",
      "metadata": {},
      "execution_count": 1,
      "source": [
        "\"\"\"Why parallel sums disagree: tracing FP32 reduction order\n",
        "Deterministic, standard-library experiment. RESULTS contains the full data.\n",
        "\"\"\"\n",
        "import math\n",
        "import struct\n",
        "\n",
        "def f32(value):\n",
        "    return struct.unpack('f', struct.pack('f', value))[0]\n",
        "\n",
        "def serial(values):\n",
        "    total = f32(0)\n",
        "    for value in values:\n",
        "        total = f32(total + value)\n",
        "    return total\n",
        "\n",
        "def tree(values):\n",
        "    work = [f32(value) for value in values]\n",
        "    while len(work) > 1:\n",
        "        if len(work) % 2:\n",
        "            work.append(0.0)\n",
        "        work = [f32(work[i] + work[i + 1])\n",
        "                for i in range(0, len(work), 2)]\n",
        "    return work[0] if work else 0.0\n",
        "\n",
        "for values in ([1e8, 1, -1e8, 1], [1e8, -1e8, 1, 1]):\n",
        "    reference = math.fsum(values)\n",
        "    print(f'serial={serial(values):.1f} tree={tree(values):.1f} reference={reference:.1f}')\n",
        "assert serial([1e8, 1, -1e8, 1]) == 1\n",
        "assert tree([1e8, 1, -1e8, 1]) == 0\n",
        "assert tree([]) == 0\n",
        "print('Rounding-order demonstration, not a GPU reduction benchmark.')\n",
        "\n",
        "\n",
        "from itertools import permutations\n",
        "rows=[]\n",
        "for values in sorted(set(permutations([1e8,1.,-1e8,1.]))):\n",
        "    reference=math.fsum(values); a=serial(values); b=tree(values)\n",
        "    rows.append(dict(inputs=list(values),serial=a,pairwise=b,reference=reference,serial_abs_error=abs(a-reference),pairwise_abs_error=abs(b-reference),sum_abs=sum(abs(v) for v in values)))\n",
        "assert len(rows)==12\n",
        "RESULTS=dict(permutations=rows,serial_results=sorted(set(r['serial'] for r in rows)),pairwise_results=sorted(set(r['pairwise'] for r in rows)),empty_tree=tree([]))\n",
        "print('Distinct input permutations:',len(rows))\n",
        "print('Serial outcomes:',RESULTS['serial_results'],'pairwise outcomes:',RESULTS['pairwise_results'])\n"
      ],
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "serial=1.0 tree=0.0 reference=2.0\n",
            "serial=2.0 tree=2.0 reference=2.0\n",
            "Rounding-order demonstration, not a GPU reduction benchmark.\n",
            "Distinct input permutations: 12\n",
            "Serial outcomes: [0.0, 1.0, 2.0] pairwise outcomes: [0.0, 2.0]\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "from pathlib import Path\n",
        "import base64\n",
        "print('__NOTEBOOK_IMAGE__svg+xml:' + base64.b64encode(Path('figure.svg').read_bytes()).decode())\n"
      ]
    }
  ]
}