eprb_signal_processing.adb

--********************************************************************
-- This is free and unencumbered software released into the public domain.
--
-- Anyone is free to copy, modify, publish, use, compile, sell, or
-- distribute this software, either in source code form or as a compiled
-- binary, for any purpose, commercial or non-commercial, and by any
-- means.
--
-- In jurisdictions that recognize copyright laws, the author or authors
-- of this software dedicate any and all copyright interest in the
-- software to the public domain. We make this dedication for the benefit
-- of the public at large and to the detriment of our heirs and
-- successors. We intend this dedication to be an overt act of
-- relinquishment in perpetuity of all present and future rights to this
-- software under copyright law.
--
-- THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
-- EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
-- MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
-- IN NO EVENT SHALL THE AUTHORS BE LIABLE FOR ANY CLAIM, DAMAGES OR
-- OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
-- ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
-- OTHER DEALINGS IN THE SOFTWARE.
--
-- For more information, please refer to <https://unlicense.org>
----------------------------------------------------------------------

-- A ‘local realistic’ simulation of the ‘two-channel Bell test’ of
-- https://en.wikipedia.org/w/index.php?title=CHSH_inequality&oldid=1170465048#Experiments,
-- demonstrating that ‘Bell’s theorem’, which is described in the
-- following reference, is wrong:
--
-- [1] J. S. Bell, ‘Bertlmann’s socks and the nature of reality’,
--     preprint, CERN-TH-2926 (1980).
--     http://cds.cern.ch/record/142461/ (Open access, CC BY 4.0)
--
-- Bell commits the infamous fallacy of assuming that causal
-- independence implies logical independence, and asserts this
-- assumption as an axiom in lieu of Bayes’ rule. In other words, he
-- confuses correlation and causation. Therefore his ‘theorem’ is
-- without foundation. One must never, ever ‘separate a and b’ the way
-- Bell did. But a simulation that contradicts Bell’s conclusion would
-- be an even stronger demonstration that he was wrong.
--
----------------------------------------------------------------------
--
-- The simulation is written in Ada. A free software Ada compiler is
-- widely available: GCC. Many readers can compile this program,
-- optimized and with runtime checks, by saving it in a file called
-- ‘eprb_wikipedia.adb’ and then running the command
--
--   gnatmake -O2 -gnata eprb_signal_processing
--
-- which will create an executable program called
-- ‘eprb_signal_processing’.  Alternatively, translate the program
-- into the language of your choice.
--
----------------------------------------------------------------------

pragma ada_2022;
pragma wide_character_encoding (utf8);

with ada.assertions;
with ada.wide_wide_text_io;
with ada.containers;
with ada.containers.doubly_linked_lists;
with ada.numerics;
with ada.numerics.generic_elementary_functions;
with ada.numerics.generic_complex_types;

----------------------------------------------------------------------

--
-- In this simulation, the problem is spoken of as one of ‘signal
-- processing’ rather than as an experiment with particles. The
-- author’s educational background is in electrical engineering, and
-- he recognizes that the issues at hand really have nothing to do
-- with particle physics, and are concerned with random signal
-- analysis. The simulation thus is formulated as a problem in that
-- domain.
--
-- Nevertheless, the simulation is suggestive of what actually IS
-- going on with physical entities. Electrons, photons, etc., probably
-- are never in the least bit ‘entangled’, but instead are correlated
-- by phase relationships of sinusoids.
--

procedure eprb_signal_processing is

  -- A ‘scalar’ is a double precision floating point number.
  type scalar is digits 15;

  subtype scalar_in_0_1 is scalar range 0.0 .. 1.0;
  subtype correlation_coefficient is scalar range -1.0 .. 1.0;

  use ada.assertions;
  use ada.wide_wide_text_io;
  use ada.numerics;
  use ada.containers;

  package scalar_elementary_functions is
    new ada.numerics.generic_elementary_functions (scalar);
  use scalar_elementary_functions;

  package scalar_complexes is
    new ada.numerics.generic_complex_types (scalar);
  use scalar_complexes;

  package scalar_io is new float_io (scalar);
  use scalar_io;

  π     : constant scalar := pi;
  π_2   : constant scalar := π / 2.0;
  π_3   : constant scalar := π / 3.0;
  π_4   : constant scalar := π / 4.0;
  π_6   : constant scalar := π / 6.0;
  π_8   : constant scalar := π / 8.0;
  π_180 : constant scalar := π / 180.0;
  two_π : constant scalar := 2.0 * π;

  subtype tuple_range is integer range 1 .. 2;

----------------------------------------------------------------------

  -- For the sake of reproducibility, let us write our own random
  -- number generator. It will be a simple linear congruential
  -- generator. The author has used one like it, in quicksorts and
  -- quickselects to select the pivot. It is good enough for our
  -- purpose.

  type uint64 is mod 2 ** 64;

  -- The multiplier lcg_a comes from Steele, Guy; Vigna, Sebastiano
  -- (28 September 2021). ‘Computationally easy, spectrally good
  -- multipliers for congruential pseudorandom number generators’.
  -- arXiv:2001.05304v3 [cs.DS]

  lcg_a : constant uint64 := 16#F1357AEA2E62A9C5#;

  -- The value of lcg_c is not critical, but should be odd.

  lcg_c : constant uint64 := 1;

  seed  : uint64 := 0;

  --
  -- random_scalar: returns a non-negative scalar less than 1.
  --
  function random_scalar
  return scalar_in_0_1
  with post => random_scalar'result < 1.0 is
    randval : scalar;
  begin
    -- Take the high 48 bits of the seed and divide it by 2**48.
    randval := scalar (seed / (2**16)) / scalar (2**48);

    -- Update the seed.
    seed := (lcg_a * seed) + lcg_c;

    return randval;
  end random_scalar;

----------------------------------------------------------------------

  --
  -- A SIGNAL is a complex number on the unit circle.
  --

  subtype SIGNAL is complex;

  --
  -- The SIGNAL_SOURCE transmits a random signal.
  --

  function SIGNAL_SOURCE
  return SIGNAL is
  begin
    return compose_from_polar (1.0, random_scalar * two_π);
  end SIGNAL_SOURCE;

  --
  -- A SIGNAL_EXCHANGE receives a SIGNAL and, if it can decide which
  -- of two channels to re-transmit the SIGNAL along, does so. The
  -- channel is indicated indicated by a function output of -1 or
  -- 1. If the SIGNAL_EXCHANGE cannot decide on a channel, the
  -- function output is 0. In that case, the SIGNAL can be considered
  -- lost, as unintelligible due to ‘noise’, ‘uncontrolled variables’,
  -- etc.
  --
  -- A SIGNAL_EXCHANGE has a complex number parameter ζ (zeta).
  --

  function SIGNAL_EXCHANGE (ζ : complex;
                            a : SIGNAL)
  return integer
  with post => (SIGNAL_EXCHANGE'result = -1 or
                SIGNAL_EXCHANGE'result = 0 or
                SIGNAL_EXCHANGE'result = 1) is

    α1   : constant scalar := re (ζ) * re (a);
    α2   : constant scalar := im (ζ) * im (a);
    α    : constant scalar := α1 + α2;

    β1   : constant scalar := -re (ζ) * im (a);
    β2   : constant scalar := im (ζ) * re (a);
    β    : constant scalar := β1 + β2;

    p_α  : constant boolean := (random_scalar < α ** 2);
    p_β  : constant boolean := (random_scalar < β ** 2);

    xchg : integer;
  begin
    if p_α then
      if p_β then
        xchg := 0;
      else
        xchg := 1;
      end if;
    else
      if p_β then
        xchg := -1;
      else
        xchg := 0;
      end if;
    end if;
    return xchg;
  end SIGNAL_EXCHANGE;

  --
  -- The SIGNAL is sent through two different SIGNAL_EXCHANGE and thus
  -- is looked for at four different receivers. We keep looking until
  -- a signal is received coincidentally at two of the
  -- receivers. Which two receivers got the signal coincidentally is
  -- recorded. For the sake of completeness, we also assume the
  -- receivers are very reliable, and so record what the signal was.
  --
  -- (A record of what the signals were can be used for more
  -- sophisticated mathematical analyses that are not yet undertaken
  -- in this program.)
  --

  type EVENT_RECORD is
    record
      r1 : integer;             -- Receiver 1: either -1 or 1.
      r2 : integer;             -- Receiver 2: either -1 or 1.
      a  : SIGNAL;
    end record;

  package EVENT_RECORD_LISTS is
    new ada.containers.doubly_linked_lists
      (element_type => EVENT_RECORD);

  function SIMULATE_EVENT (ζ1, ζ2 : complex)
  return EVENT_RECORD
  with post => ((SIMULATE_EVENT'result.r1 = -1 or
                 SIMULATE_EVENT'result.r1 = 1) and
                (SIMULATE_EVENT'result.r2 = -1 or
                 SIMULATE_EVENT'result.r2 = 1)) is
    ev : EVENT_RECORD;

    procedure fill_ev is
    begin
      ev.a := SIGNAL_SOURCE;
      ev.r1 := SIGNAL_EXCHANGE (ζ1, ev.a);
      ev.r2 := SIGNAL_EXCHANGE (ζ2, ev.a);
    end fill_ev;

  begin
    fill_ev;
    while ev.r1 = 0 or ev.r2 = 0 loop
      fill_ev;
    end loop;
    return ev;
  end SIMULATE_EVENT;

  function SIMULATE_RUN (ζ1, ζ2     : complex;
                         run_length : count_type)
  return EVENT_RECORD_LISTS.list is
    use EVENT_RECORD_LISTS;
    events : list := empty_list;
  begin
    for i in 1 .. run_length loop
      append (events, SIMULATE_EVENT (ζ1, ζ2));
    end loop;
    return events;
  end SIMULATE_RUN;

-- ----------------------------------------------------------------------

  --
  -- A basic function in analyzing run data is to count how many times
  -- a particular receiver pair occurs.
  --

  function count (events : EVENT_RECORD_LISTS.list;
                  r1, r2 : integer)
  return count_type
  with pre => (r1 = -1 or r1 = 1) and (r2 = -1 or r2 = 1) is
    use EVENT_RECORD_LISTS;
    n    : count_type := 0;
    curs : cursor := first (events);
  begin
    while has_element (curs) loop
      declare
        ev : EVENT_RECORD := element (curs);
      begin
        if ev.r1 = r1 and ev.r2 = r2 then
          n := n + 1;
        end if;
      end;
      curs := next (curs);
    end loop;
    return n;
  end count;

-- ----------------------------------------------------------------------

  --
  -- CHSH compute a primitive kind of correlation coefficient that, in
  -- our case, is the difference between the frequency of EVENT_RECORD
  -- with r1=r2 and the frequency of EVENT_RECORD with r1≠r2.
  --

  function chsh_correlation (events : EVENT_RECORD_LISTS.list)
  return correlation_coefficient is
    n1 : constant count_type := count (events, -1, -1)
                                  + count (events, 1, 1);
    n2 : constant count_type := count (events, -1, 1)
                                  + count (events, 1, -1);
  begin
    return scalar (n1 - n2) / scalar (n1 + n2);
  end chsh_correlation;

-- ----------------------------------------------------------------------

  --
  -- The following array of complex number tuples are SIGNAL_EXCHANGE
  -- ζ settings corresponding to the ‘Bell-test angles’.
  --

  type complex_tuple is array (tuple_range) of complex;

  subtype bell_test_settings_range is integer range 1 .. 4;
  type bell_test_settings_array is
    array (bell_test_settings_range) of complex_tuple;

  bell_test_settings : constant bell_test_settings_array :=

    ((compose_from_polar (1.0, 0.0),
      compose_from_polar (1.0, π_8)),

     (compose_from_polar (1.0, 0.0),
      compose_from_polar (1.0, 3.0 * π_8)),

     (compose_from_polar (1.0, π_4),
      compose_from_polar (1.0, π_8)),

     (compose_from_polar (1.0, π_4),
      compose_from_polar (1.0, 3.0 * π_8)));

----------------------------------------------------------------------

begin
  declare
    κ : array (bell_test_settings_range) of correlation_coefficient;
  begin
    for i in bell_test_settings_range loop
      κ(i) := chsh_correlation (SIMULATE_RUN
                                  (bell_test_settings(i)(1),
                                   bell_test_settings(i)(2),
                                   1e6));
    end loop;
    put ("   κ₁ : "); put (κ(1), 2, 5, 0); new_line;
    put ("   κ₂ : "); put (κ(2), 2, 5, 0); new_line;
    put ("   κ₃ : "); put (κ(3), 2, 5, 0); new_line;
    put ("   κ₄ : "); put (κ(4), 2, 5, 0); new_line;
    put ("    S : ");
    put (κ(1) - κ(2) + κ(3) + κ(1), 2, 5, 0);
    new_line;
  end;
end eprb_signal_processing;

----------------------------------------------------------------------
--
-- Output:
--
--     κ₁ :  0.62836
--     κ₂ : -0.62756
--     κ₃ :  0.62693
--     κ₄ :  0.62903
--      S :  2.51122
--
--
--********************************************************************
-- Some instructions for the Emacs text editor.
-- local variables:
-- mode: indented-text
-- tab-width: 2
-- end: