eprb_simulation.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 simulation of an EPR-B experiment, demonstrating that ‘Bell’s
-- theorem’ assumes its conclusion and therefore is fallacious.
-- (https://en.wikipedia.org/w/index.php?title=Begging_the_question&oldid=1163995745)

-- Bell’s conclusion is not proven. Indeed, because the proof is by
-- counterexample, the contrary IS proven: hidden variable theories
-- ARE possible. Therefore ‘Bell’s theorem’ is shown to be worthless,
-- even if it were not begging the question: its premises are
-- demonstrably false.

-- Let us list the references here, near the top:

-- [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)

-- [2] A. F. Kracklauer, ‘EPR-B correlations: non-locality or
--     geometry?’, J. Nonlinear Math. Phys. 11 (Supp.) 104–109
--     (2004). https://doi.org/10.2991/jnmp.2004.11.s1.13 (Open
--     access, CC BY-NC)

-- The simulation is written in Ada—which, if carefully used, seems to
-- me a good language for conveying scientific algorithms. Also, 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_simulation.adb’ and then running
-- the command

--   gnatmake -O2 -gnata eprb_simulation

-- which will create an executable program called ‘eprb_simulation’.

-- Ada, simply put, seemed the best choice, among the programming
-- languages I felt comfortable using.

-- However, if you know another language, such as Fortran, Python, or
-- Common Lisp, it would be very useful for the reader’s understanding
-- to translate the program from Ada into the other language.

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

pragma ada_2012;
pragma wide_character_encoding (utf8);

with ada.exceptions;
with ada.text_io;
with ada.containers.ordered_sets;
with ada.numerics;
with ada.numerics.generic_elementary_functions;

procedure eprb_simulation is

  subtype range_1_to_2 is integer range 1 .. 2;

  type scalar is digits 15;     -- double precision
  type scalar_pair is array (range_1_to_2) of scalar;

  use ada.text_io;
  use ada.numerics;

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

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

  π   : constant scalar := pi;
  π_2 : constant scalar := π / 2.0;
  π_4 : constant scalar := π / 4.0;
  π_8 : constant scalar := π / 8.0;

  programmer_mistake : exception;

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

  -- 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;

  function random_scalar
  return scalar
  with post => 0.0 <= random_scalar'result
                and 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;

  function random_1_or_2
  return range_1_to_2 is
  begin
    return (if random_scalar < 0.5 then 1 else 2);
  end random_1_or_2;

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

  -- The simulation is of a two-arm Clauser-Aspect type
  -- experiment. See [2] for description of a closely related
  -- simulation. We will be referring to equation numbers in that
  -- paper.

  -- The objects of our simulation are as follows.

  -- There are light pulses that we shall call photons. These are
  -- randomly generated in pairs, with either a vertically polarized
  -- pulse to the left and a horizontally polarized pulse to the
  -- right, or vice versa. Call the polarization state φ. In our
  -- simulation, these photons are assumed also to carry ‘hidden
  -- variables’ called ξ. A photon generated with vertical
  -- polarization has ξ=ξa, and a photon with horizontal polarization
  -- has ξ=ξb. Thus the following is our representation:

  type ξ_value is (ξa, ξb, ξ_quiescent);
  type photon is
    record
      φ : scalar;       -- The current polarization angle, in radians.
      ξ : ξ_value;      -- The hidden variable. Either ξa or ξb.
    end record;

  -- The following are the only possible initial photons:

  photon_vertical   : constant photon := (φ => π_2, ξ => ξa);
  photon_horizontal : constant photon := (φ => 0.0, ξ => ξb);

  -- Later we will introduce special photodetectors that can measure
  -- the value of a photon’s ξ. The quiescent state of such a
  -- photodetector will be ξ_quiescent.

  -- Please note something carefully here: by including a ‘hidden
  -- variable’ in the simulation of a photon, we are contradicting an
  -- assumption in [1], which Bell attempts to justify on pages 15 and
  -- 16. To quote [1]: ‘[I]t may be that it is not permissible to
  -- regard the experimental settings a and b in the analyzers as
  -- independent variables, as we did.’ Merely imagining hidden
  -- variables, along with photodetectors able to measure the hidden
  -- variables, directly contradicts ‘regarding the experimental
  -- settings a and b ... as independent variables’. That we can
  -- imagine such hidden variables therefore proves beyond doubt that
  -- Bell’s argument is founded on a false assumption! Thus his
  -- argument as a whole fails.

  -- So we have done enough already. But let us proceed with the
  -- simulation, nonetheless. Doing so will help illustrate why simply
  -- assuming the existence of a photon with hidden variables
  -- contradicts Bell’s assumption—or, in other words, how Bell’s
  -- argument is circular. Bell does the following:

  --   1. He assumes that a particle must be an ‘independent
  --   variable’.

  --   2. This is equivalent to assuming that a particle has no
  --   ‘hidden variables’.

  --   3. He performs various mathematical operations from this
  --   assumption, concluding that a particle cannot have any hidden
  --   variables.

  -- It is a classic example of someone accidentally assuming the
  -- conclusion, which is an informal fallacy. There is no deductive
  -- error, and yet the conclusion remains unproven.

  -- (Side note: the contrast is striking, between our simple disproof
  -- by counterexample, and the sophistries at the top of page 16 of
  -- [1]. Notable is Bell’s reliance on his physical intuitions,
  -- employing no principles of logic and mathematics. Also he
  -- overcomplicates the matter, seemingly in a struggle to visualize
  -- the situation. All that one needs to visualize is a hidden
  -- variable and a detector capable of measuring it! Bell obviously
  -- did not realize he had to visualize a photon WITH a hidden
  -- variable, rather than the photon WITHOUT a hidden variable that
  -- he was used to imagining.)

  -- Next after photons, there are two polarizing beam splitters (PBS)
  -- (https://en.wikipedia.org/w/index.php?title=Polarizer&oldid=1152014034#Beam-splitting_polarizers).
  -- These redirect photons in proportion to the Law of Malus
  -- (https://en.wikipedia.org/w/index.php?title=Polarizer&oldid=1152014034#Malus's_law_and_other_properties),
  -- altering the photons’ polarization angles φ, but not their hidden
  -- variables ξ.

  type pbs is
    record
      θ : scalar;               -- The angle of beam splitting.
    end record;

  -- These are the different PBS settings, in radians, making up four
  -- pairs of (left,right) settings, chosen for maximum nominal CHSH
  -- contrast (=2√2):

  θ_left  : constant scalar_pair := (π_4, 0.0);
  θ_right : constant scalar_pair := (π_8, 3.0 * π_8);

  -- Without loss of generality, we have specified the PBS angles in
  -- the first quadrant. PBS geometry is actually a sort of ×, where
  -- one of the two slopes represents one of the channels, and the
  -- crossing slope represents the other.

  -- Finally, there are four photodetectors
  -- (https://en.wikipedia.org/w/index.php?title=Photodetector&oldid=1168137030).
  -- These are fancy new-model photodetectors: not only do they detect
  -- incident photons without fail, but also they tell you the value
  -- of the photon’s ξ. Their output is a ξ_value: the quiescent state
  -- is ξ_quiescent, and they go to ξa or ξb to register a photon
  -- detection. They can detect only one photon at a time, but the
  -- experimental arrangement is such that they will encounter at most
  -- one at a time.

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

  -- Here is a procedure that generates a (left,right) pair of
  -- photons.

  procedure generate_photon_pair (photon_left  : out photon;
                                  photon_right : out photon) is
  begin
    if random_1_or_2 = 1 then
      photon_left := photon_vertical;
      photon_right := photon_horizontal;
    else
      photon_left := photon_horizontal;
      photon_right := photon_vertical;
    end if;
  end generate_photon_pair;

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

  -- The following procedure decides towards which of two
  -- photodetectors a PBS will transmit a given photon. The output of
  -- the procedure is the outputs of the two photodetectors. The
  -- photon will have its angle of polarization altered, but the
  -- photodetectors cannot measure angle of polarization. The outputs
  -- are of what the photodetectors CAN measure: hidden variable
  -- values.

  procedure split_beam (photon0   : photon;
                        pbs0      : pbs;
                        detector1 : out ξ_value;
                        detector2 : out ξ_value) is
  begin
    detector1 := ξ_quiescent;
    detector2 := ξ_quiescent;

    -- We implement the Law of Malus this way, which I think is
    -- clearer than if I wrote it in terms of differences of angles—

    -- If the incident photon is polarized horizontally, then the
    -- probability of transmission along channel 1 will be cos²(θ),
    -- where θ is the pbs setting. Thus the probability it will
    -- instead be transmitted along channel 2 is sin²(θ), and the
    -- total probability of transmission is cos²(θ)+sin²(θ)=1.

    -- On the other hand, if the incident photon is polarized
    -- vertically, the roles of cosine and sine will be reversed.
    -- (This role reversal results from application of the identities
    -- cos(π/2-x)≡sin(x) and sin(π/2-x)≡cos(x).) The probability of
    -- transmission along channel 1 is sin²(θ), along channel 2 is
    -- cos²(θ).

    if abs (photon0.φ) < π_4 * scalar'model_epsilon then
      -- The photon is polarized horizontally.
      if random_scalar < cos (pbs0.θ) ** 2 then
        detector1 := photon0.ξ;
      else
        detector2 := photon0.ξ;
      end if;
    elsif abs (photon0.φ - π_2) < π_4 * scalar'model_epsilon then
      -- The photon is polarized vertically.
      if random_scalar < sin (pbs0.θ) ** 2 then
        detector1 := photon0.ξ;
      else
        detector2 := photon0.ξ;
      end if;
    else
      -- There should be only horizontally or vertically polarized
      -- photons.
      raise programmer_mistake;
    end if;
  end split_beam;

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

  -- The function one_event, below, simulates an event of the
  -- simulation: from generation of a pair of photons to their
  -- measurement by four photodetectors. The settings of the two PBS
  -- are chosen randomly and included in the output. The output is
  -- gathered into a record and given a unique identifier (so it can
  -- be stored in set objects).

  current_identifier : long_integer := 1;

  type event_record is
    record
      identifier        : long_integer;
      pbs_left_setting  : range_1_to_2;
      pbs_right_setting : range_1_to_2;
      detector_left_1   : ξ_value;
      detector_left_2   : ξ_value;
      detector_right_1  : ξ_value;
      detector_right_2  : ξ_value;
    end record;

  function one_event
  return event_record is
    ev           : event_record;
    pbs_left     : pbs;
    pbs_right    : pbs;
    photon_left  : photon;
    photon_right : photon;
  begin
    ev.identifier := current_identifier;
    current_identifier := current_identifier + 1;

    ev.pbs_left_setting := random_1_or_2;
    pbs_left.θ := θ_left(ev.pbs_left_setting);

    ev.pbs_right_setting := random_1_or_2;
    pbs_right.θ := θ_right(ev.pbs_right_setting);

    generate_photon_pair (photon_left, photon_right);

    split_beam (photon_left, pbs_left,
                ev.detector_left_1, ev.detector_left_2);
    split_beam (photon_right, pbs_right,
                ev.detector_right_1, ev.detector_right_2);

    return ev;
  end one_event;

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

  -- We wish to simulate some number of events, and to gather the
  -- output records in a set. Here is where the set type is created.

  function "<" (a, b : event_record)
  return boolean is
  begin
    return (a.identifier < b.identifier);
  end "<";

  package event_record_sets is
    new ada.containers.ordered_sets (element_type => event_record);

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

  -- The following function simulates a given number of events and
  -- gathers the output records into a set.

  function simulate_events (number_of_events : natural)
  return event_record_sets.set is
    events : event_record_sets.set;
  begin
    for i in 1 .. number_of_events loop
      event_record_sets.include (events, one_event);
    end loop;
    return events;
  end simulate_events;

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

  -- Now comes analysis.

  procedure analyze (events : event_record_sets.set) is

    use event_record_sets;

    events_L1R1 : set;
    events_L1R2 : set;
    events_L2R1 : set;
    events_L2R2 : set;

    procedure separate_into_four_subsets_by_pbs_settings is
      curs              : cursor;
      ev                : event_record;
    begin
      curs := first (events);
      while has_element (curs) loop
        ev := element (curs);
        case ev.pbs_left_setting is
          when 1 =>
            case ev.pbs_right_setting is
              when 1 => include (events_L1R1, ev);
              when 2 => include (events_L1R2, ev);
            end case;
          when 2 =>
            case ev.pbs_right_setting is
              when 1 => include (events_L2R1, ev);
              when 2 => include (events_L2R2, ev);
            end case;
        end case;
        curs := next (curs);
      end loop;
    end separate_into_four_subsets_by_pbs_settings;

    function estimate_correlation_coefficient (subset : set)
    return scalar is
      curs                   : cursor;
      ev                     : event_record;
      n, nLc, nLs, nRc, nRs  : scalar;
      cosL, sinL, cosR, sinR : scalar;
      cosLR, sinLR           : scalar;
    begin

      -- We use equation (2.4) of [2] to estimate cosines and sines
      -- for calculation of the correlation coefficient by equations
      -- (2.3) and (2.5).
      --
      -- We take account of the hidden variable in this way: if ξ=ξa,
      -- this indicates that the PBS followed the sine-squared law for
      -- channel 1, rather than the cosine-squared law. For the hidden
      -- variable tells us that the photon had a vertical angle of
      -- polarization when it entered the PBS. We must, therefore,
      -- count a ξ=ξa detection towards the sine estimate. By
      -- contrast, and symmetrically, a ξ=ξb detection goes towards
      -- the cosine estimate.
      --
      -- If ξ=ξ_quiescent then a detector did not receive a photon.

      n := 0.0;
      nLc := 0.0;
      nLs := 0.0;
      nRc := 0.0;
      nRs := 0.0;

      curs := first (subset);

      while has_element (curs) loop

        ev := element (curs);
        if ev.detector_left_1 /= ξ_quiescent then
          if ev.detector_left_1 = ξa then
            nLs := nLs + 1.0;
          else
            nLc := nLc + 1.0;
          end if;
        elsif ev.detector_left_2 /= ξ_quiescent then
          if ev.detector_left_2 = ξb then
            nLs := nLs + 1.0;
          else
            nLc := nLc + 1.0;
          end if;
        else
          -- The way the program currently works, this should not
          -- happen.
          raise programmer_mistake;
        end if;

        if ev.detector_right_1 /= ξ_quiescent then
          if ev.detector_right_1 = ξa then
            nRs := nRs + 1.0;
          else
            nRc := nRc + 1.0;
          end if;
        elsif ev.detector_right_2 /= ξ_quiescent then
          if ev.detector_right_2 = ξb then
            nRs := nRs + 1.0;
          else
            nRc := nRc + 1.0;
          end if;
        else
          -- The way the program currently works, this should not
          -- happen.
          raise programmer_mistake;
        end if;

        n := n + 1.0;

        curs := next (curs);

      end loop;

      cosL := sqrt (nLc / n);
      sinL := sqrt (nLs / n);
      cosR := sqrt (nRc / n);
      sinR := sqrt (nRs / n);

      -- Reference [2], Equation (2.3).
      cosLR := (cosR * cosL) + (sinR * sinL);
      sinLR := (sinR * cosL) - (cosR * sinL);

      -- Reference [2], Equation (2.5).
      return (cosLR ** 2) - (sinLR ** 2);

    end estimate_correlation_coefficient;

    κ_L1R1, κ_L1R2      : scalar;
    κ_L2R1, κ_L2R2      : scalar;
    κ_nominal           : constant scalar := sqrt (0.5);
    chsh_contrast       : scalar;
    chsh_nominal        : constant scalar := 2.0 * sqrt (2.0);
    chsh_difference     : scalar;
    relative_difference : scalar;

  begin
    separate_into_four_subsets_by_pbs_settings;
    κ_L1R1 := estimate_correlation_coefficient (events_L1R1);
    κ_L1R2 := estimate_correlation_coefficient (events_L1R2);
    κ_L2R1 := estimate_correlation_coefficient (events_L2R1);
    κ_L2R2 := estimate_correlation_coefficient (events_L2R2);

    chsh_contrast := κ_L1R1 + κ_L1R2 + κ_L2R1 - κ_L2R2;
    chsh_difference := chsh_contrast - chsh_nominal;
    relative_difference := chsh_difference / chsh_nominal;

    put_line ("  =========================");
    put_line ("   no. of events:  " &
              event_record_sets.length (events)'image);
    put_line ("  =========================");
    put_line ("   correlation coefficients");
    put ("          kappa11");
    put (κ_L1R1, 4, 5, 0);
    put_line ("");
    put ("          kappa12");
    put (κ_L1R2, 4, 5, 0);
    put_line ("");
    put ("          kappa21");
    put (κ_L2R1, 4, 5, 0);
    put_line ("");
    put ("          kappa22");
    put (κ_L2R2, 4, 5, 0);
    put_line ("");
    put ("          nominal");
    put (κ_nominal, 4, 5, 0);
    put_line ("");
    put_line ("  =========================");
    put ("    CHSH contrast");
    put (chsh_contrast, 4, 5, 0);
    put_line ("");
    put ("          nominal");
    put (chsh_nominal, 4, 5, 0);
    put_line ("");
    put_line ("       --------------------");
    put ("       difference");
    put (chsh_difference, 4, 5, 0);
    put_line ("");
    put ("    relative diff");
    put (relative_difference, 4, 5, 0);
    put_line ("");
    put_line ("  =========================");
  end analyze;

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

-- Finally, the main program.

  number_of_events : natural := 0;

begin

  -- Simulate many events.

  seed := 12345;
  number_of_events := 1000000;
  analyze (simulate_events (number_of_events));

end eprb_simulation;

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

-- The condition for supposed demonstration of ‘non-locality’ is CHSH
-- contrast greater than 2. The nominal value predicted by quantum
-- mechanics is 2√2≅2.82843. Here is the program’s output for
-- seed=12345, number_of_events=1000000:

--  =========================
--   no. of events:   1000000
--  =========================
--   correlation coefficients
--          kappa11   0.70739
--          kappa12   0.70775
--          kappa21   0.70710
--          kappa22  -0.70741
--          nominal   0.70711
--  =========================
--    CHSH contrast   2.82965
--          nominal   2.82843
--       --------------------
--       difference   0.00123
--    relative diff   0.00043
--  =========================

-- As you can see, the CHSH contrast is essentially that which is
-- predicted by quantum mechanics, and certainly exceeds the value of
-- 2, which supposedly is impossible for a ‘local realistic’ model to
-- achieve. We have disproven that ‘impossibility’, by counterexample.

-- (One might object this is a random simulation, not closed
-- expressions. But closed expressions can be written by analogy to
-- the simulation. One might write them as an exercise.)

-- However, let us go deduce further. Let us imagine that the photons
-- actually do NOT have ξ values, and that the photodetectors merely
-- detect photons without fail. Does anything change in the pattern of
-- their detections? Not one bit. A PBS reacts only to a photon’s φ,
-- its polarization angle. Therefore the correlations are exactly the
-- same. We get the exact same results we did above.

-- It turns out that ξ was merely a fiction introduced so we could
-- include the initial angles of polarization in the analysis stage,
-- and yet make it so we were using only information actually measured
-- in the experiment. We never actually needed that fiction at all,
-- but could simply have included the initial angles of polarization
-- in the event records.

-- Recall the quote from [1]: ‘[I]t may be that it is not permissible
-- to regard the experimental settings a and b in the analyzers as
-- independent variables, as we did.’ In other words, it may be that
-- one must regard Bell’s a and b as functions of the initial angles
-- of polarization. This is exactly what we have shown: that one must
-- regard a and b as FUNCTIONS, and that the initial angle of
-- polarization of one of the photons of the pair may be regarded as
-- the sole independent variable.

-- Furthermore, we showed this by inventing a method whereby the
-- ACTUAL independent variable is converted into a fictitious
-- experimental ‘hidden variable’. Thus it is shown that assuming a
-- and b must be the independent variables is equivalent to assuming
-- there can be no such hidden variable. Bell has assumed away the
-- ability to construct what (by ignoring his assumption) we succeeded
-- in constructing: a ‘locally realistic’ hidden variable theory.

-- Note one thing more: we have shown that there is nothing necessary
-- to obtain the correlations but that the light be
-- polarized. ‘Entanglement’, at least in this case, is simply an
-- artifact of quantum mechanics notation.

-- (Indeed, classical electromagnetic wave theory predicts the same
-- correlations. See [2].)

--********************************************************************
-- Some instructions for the Emacs text editor.
-- local variables:
-- mode: indented-text
-- tab-width: 2
-- end: