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’, 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)

-- 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_2022;
pragma wide_character_encoding (utf8);

with ada.exceptions;
with ada.wide_wide_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.wide_wide_text_io;
  use ada.numerics;

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

  package range_1_to_2_io is new integer_io (range_1_to_2);
  use range_1_to_2_io;

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

  programmer_mistake : exception;

  function cos2 (x : scalar)
  return scalar is
  begin
    return cos (x) ** 2;
  end cos2;

  function sin2 (x : scalar)
  return scalar is
  begin
    return sin (x) ** 2;
  end sin2;

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

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

  --
  -- random_1_or_2: returns the integer 1 or the integer 2.
  --
  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 optical experiments, involving
  -- sources of plane-polarized photon pairs, polarizing beam
  -- splitters, and photodetectors. There will be two versions of the
  -- simulation.

  -- Simulation A will be a ‘local realistic’ experiment.

  -- Simulation B will be a ‘quantum entanglement’ experiment.

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

  -- In more detail, the objects of Simulation A, the ‘local
  -- realistic’ experiment, are as follows.

  -- There are light pulses that we shall call photons. These are
  -- randomly generated in pairs, which are oppositely
  -- plane-polarized, either horizontally or vertically. Call the
  -- polarization state φ. In our simulation, a photon also has a
  -- state called ρ, which is a non-negative number less than
  -- one. Thus the following is our representation:

  type φ_value is (horizontal, vertical);
  subtype ρ_value is scalar range 0.0 .. 1.0;

  type photon_A is
    record
      φ : φ_value;
      ρ : ρ_value;
    end record;

  -- 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 φ. However, we will not
  -- care about the new value of φ.

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

  -- Finally, there are four photodetectors
  -- (https://en.wikipedia.org/w/index.php?title=Photodetector&oldid=1168137030).
  -- These are ideal photodetectors, that detect incident photons
  -- without fail. 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
  -- at random. The correlated polarization angles are set by a system
  -- of randomly positioned polarizing filters, so angles can be
  -- recorded. The two polarizing filters can be separated from each
  -- other by any distance whatsoever, and there is no physical
  -- influence of one filter setting upon the other. The value of ρ is
  -- randomly generated separately for each photon, but not recorded.

  procedure generate_photon_A_pair (photon_left  : out photon_A;
                                    photon_right : out photon_A) is
  begin
    if random_1_or_2 = 1 then
      photon_left.φ := horizontal;
      photon_right.φ := vertical;
    else
      photon_left.φ := vertical;
      photon_right.φ := horizontal;
    end if;
    photon_left.ρ := random_scalar;
    photon_right.ρ := random_scalar;
  end generate_photon_A_pair;

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

  -- The following procedure decides which of two photodetectors will
  -- detect the photon that enters a given PBS. 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, so that
  -- information is simply lost.

  procedure split_beam_A (photon    : photon_A;
                          pbs       : pbs_A;
                          detector1 : out boolean;
                          detector2 : out boolean) is
    φ : scalar;
  begin

    -- The Law of Malus tells us that the the square of the cosine of
    -- the angle between the plane of polarization and the plane of
    -- the PBS setting gives us the probability of transmission to
    -- detector1.

    -- The photon’s hidden variable ρ is used as the random value.

    φ := (if photon.φ = horizontal then 0.0 else π_2);
    if photon.ρ < cos2 (φ - pbs.θ) then
      detector1 := true;
      detector2 := false;
    else
      detector1 := false;
      detector2 := true;
    end if;

  end split_beam_A;

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

  -- The function one_event_A, 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 included in the output. The polarization angles of the
  -- photons are included, as well. The output is gathered into a
  -- record and given a unique identifier (so it can be stored in set
  -- objects).

  type identifier_tag is mod 2 ** 128; -- (Set to 2**64 if necessary.)
  current_identifier : identifier_tag := 0;

  type event_record is
    record
      identifier        : identifier_tag;
      θ_left, θ_right   : scalar := 0.0;
      φ_left, φ_right   : φ_value := horizontal;
      detector_left_1   : boolean := false;
      detector_left_2   : boolean := false;
      detector_right_1  : boolean := false;
      detector_right_2  : boolean := false;
    end record;

  function one_event_A (θ_L, θ_R : scalar)
  return event_record is
    ev           : event_record;
    pbs_left     : pbs_A;
    pbs_right    : pbs_A;
    photon_left  : photon_A;
    photon_right : photon_A;
  begin
    ev.identifier := current_identifier;
    current_identifier := current_identifier + 1;

    ev.θ_left := θ_L;
    pbs_left.θ := θ_L;

    ev.θ_right := θ_R;
    pbs_right.θ := θ_R;

    generate_photon_A_pair (photon_left, photon_right);
    ev.φ_left := photon_left.φ;
    ev.φ_right := photon_right.φ;

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

    return ev;
  end one_event_A;

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

  -- 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_A (number_of_events : natural;
                              θ_L, θ_R         : scalar)
  return event_record_sets.set is
    use event_record_sets;
    events : set;
  begin
    for i in 1 .. number_of_events loop
      include (events, one_event_A (θ_L, θ_R));
    end loop;
    return events;
  end simulate_events_A;

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

  -- For Simulation B, we must totally ignore the operations of the
  -- apparatus. Quantum mechanics follows Bohr’s doctrine that eschews
  -- cause-and-effect, as outside the realm of interest. Instead
  -- quantum mechanics simply tells us what the probabilities of
  -- different measurements are for a given experimental
  -- arrangement. This is all we are supposed to care about. Thus we
  -- do the simulation by calculating those probabilities and then
  -- rolling dice.

  -- The apparatus for Simulation B will be the same as for Simulation
  -- A, and it will be operated in the same way. However, we make no
  -- assumptions that photons have ρ values. Strictly speaking, we
  -- also do not assume they have φ values. Rather, it is the
  -- polarizing filters that have φ values. We know those values, and
  -- record them in the event records.

  -- (In some experiments people have actually run, there may be
  -- analogs of our φ values that are measured at the far ends of the
  -- experimental arms, rather than at the light source. Perhaps
  -- someone will object to Simulation B on grounds of this
  -- difference. Is Bell’s argument supposed to apply only to
  -- particular experimental arrangements? Or perhaps it is that
  -- photon pairs that have passed through polarizing filters have
  -- their ‘magic entanglement’ robbed from them? My guess is the
  -- objection would be a Schrödinger’s cat of all possible
  -- challenges. However, an actually curious but doubtful reader can
  -- modify this program to handle the different experimental
  -- arrangement. Doing so would be much more convincing than if I did
  -- the modification for them.)

  -- The only other measurements in the experiment are the PBS
  -- settings and the photodetector outputs. These, too, are recorded
  -- in the event records.

  -- Our task now is to calculate probabilities of different event
  -- records. That is what the formulas of quantum mechanics provide
  -- means to do in such a situation. It is what Niels Bohr considered
  -- all that physicists ought to expect of their work.

  -- But I do not know how to do the formulas of quantum mechanics. In
  -- my education, we got an introduction to quantum mechanics, but
  -- not one deep enough to learn the special language that is used in
  -- cases such as this. Instead I will use general methods of logical
  -- inference. These will solve the problem entirely, without any
  -- specialty knowledge.

  -- (That may come as a surprise to some quantum physicists, who are
  -- not so deep inside the circle that they cannot see what is going
  -- on outside it. They may have expected it impossible to solve a
  -- quantum mechanics problem without its special language and a
  -- benediction.)

  -- First of all, we must disabuse ourselves of Bell’s immensely
  -- ignorant assertions in [1], that it is possible to declare two
  -- variables ‘statistically independent’ on grounds that one cannot
  -- have a causal effect on the other. Nor can one do so on grounds
  -- of vague intuitions. Bell is doing nothing less than displaying
  -- to us his ignorance of the methods of scientific inference. To
  -- declare ‘statistical independence’, one must do nothing less than
  -- present a mathematical proof. In our case, there is no possible
  -- such proof, because the two variables are indubitably correlated:
  -- if one is a vertical polarization, then the other is certainly a
  -- horizontal polarization, and vice versa.

  -- Thus I have already shown that Bell is wrong, and thus the entire
  -- opus based on his claims is discredited. CHSH is wrong, Aspect is
  -- wrong, and so on, and so on. But let us continue nonetheless.

  -- We have disabused ourselves of Bell’s assertions, and instead
  -- will use the correct expressions. A scientist or engineer in
  -- practically any field except quantum physics would do the same,
  -- because these are general methods of logical inference.

  -- We will use notation as follows (and analogously throughout):

  --    P(hL1)    probability of horizontal left polarizing filter
  --              and a detection in left channel 1
  --    P(vR2)    probability of a vertical right polarizing filter
  --              and a detection in right channel 2
  --    P(vL2,vR1)  joint probability of vL2 and vR1
  --    P(vL2|vR1)  conditional probability of vL2, given vR1

  -- Calculations will be done separately for each pair of PBS
  -- settings, and thus the PBS settings need not be considered in the
  -- probabilities.

  -- To do LOGICAL INFERENCE correctly, one must use MATHEMATICAL
  -- THEOREMS, not physical or intuitive justifications such as Bell
  -- attempts. Thus it is not legal to write something like

  --    P(hL1,hR1) = P(hL1) × P(hR1)

  -- as Bell does. It is simply not allowed. There does not exist, in
  -- mathematics, any theorem one can cite to justify such a step. One
  -- MUST instead write one or the other of

  --    P(hL1,hR1) = P(hL1) × P(hR1|hL1)
  --    P(hL1,hR1) = P(hL1|hR1) × P(hR1)

  -- This is a form of Bayes’ theorem
  -- (https://en.wikipedia.org/w/index.php?title=Bayes%27_theorem&oldid=1171865590). (See
  -- my footnote on ‘Bayes’ theorem’.)

  -- Notice that the value of P(hL1,hR1) according to Bell’s incorrect
  -- calculation is a positive number. On the other hand, according to
  -- the correct expressions, the value is zero: it is impossible to
  -- have horizontally polarized filters simultaneously on both sides,
  -- and therefore the conditional probability factors are zero.

  -- Already we see a tangible demonstration that we cannot expect to
  -- get correct answers, if we reason as Bell does. Obviously there
  -- is no causal effect of one polarizing filter on the other. Their
  -- angles are correlated solely because they are set that way by a
  -- third device, whose nature is unimportant (though it probably
  -- would involve a system of shutters). Even so, we cannot write the
  -- expression Bell claims we can. We MUST use Bayes’ theorem.

  -- (Here is another example of the need to use Bayes’ theorem rather
  -- than arguments about causal influences—

  -- My father and I both uploaded some data about our DNAs to the
  -- same ‘family tree’ service. The service confirmed that our DNAs
  -- were about 50 percent similar, in whatever it was they
  -- measured. I am not a molecular biologist and do not understand
  -- these things. All I know is the two DNAs proved to be extremely
  -- alike. Now think about that coincidence. There was this powerful
  -- correlation between two samples of deoxyribonucleic acid. The
  -- only connection between them is they had a common origin in one
  -- organism, my father. There was no possible causal influence of
  -- one sample on the other.

  -- By Bell’s reasoning, the only explanation for the DNA service’s
  -- result is a spooky kind of action at a distance!

  -- Of course, in the field of genomics, such a conclusion would not
  -- be publishable. If a renowned genomicist gave a lecture and came
  -- to such a conclusion, it might be wise to call for
  -- paramedics. However, strange beliefs may be due to indoctrination
  -- rather than organic or conventional psychiatric causes, so the
  -- situation in quantum theory is different.)

  -- The experiment, then, is a very simple one. The correlations may
  -- be characterized thus: any conditional probability that has a ‘v’
  -- on both sides of the ‘|’ must equal zero. The same is true if
  -- instead of ‘v’ on both sides there is an ‘h’. Because excluding
  -- the possibility of either ‘h’ or ‘v’ from one arm removes half of
  -- the possible cases to choose from, the other conditional
  -- probabilities must each equal twice their equivalents without the
  -- condition. (That is, the denominator of the fraction is halved.)

  -- So let us write out the full set of probabilities, in the P
  -- notation:

  --    P(hL1,hR1) = 0
  --    P(hL1,hR2) = 0
  --    P(hL2,hR1) = 0
  --    P(hL2,hR2) = 0
  --    P(vL1,vR1) = 0
  --    P(vL1,vR2) = 0
  --    P(vL2,vR1) = 0
  --    P(vL2,vR2) = 0
  --    P(hL1,vR1) = P(hL1) × P(vR1|hL1) = P(hL1) × (2 × P(vR1))
  --    P(hL1,vR2) = P(hL1) × P(vR2|hL1) = P(hL1) × (2 × P(vR2))
  --    P(hL2,vR1) = P(hL2) × P(vR1|hL2) = P(hL2) × (2 × P(vR1))
  --    P(hL2,vR2) = P(hL2) × P(vR2|hL2) = P(hL2) × (2 × P(vR2))
  --    P(vL1,hR1) = P(vL1) × P(hR1|vL1) = P(vL1) × (2 × P(hR1))
  --    P(vL1,hR2) = P(vL1) × P(hR2|vL1) = P(vL1) × (2 × P(hR2))
  --    P(vL2,hR1) = P(vL2) × P(hR1|vL2) = P(vL2) × (2 × P(hR1))
  --    P(vL2,hR2) = P(vL2) × P(hR2|vL2) = P(vL2) × (2 × P(hR2))

  -- It remains to replace the P notations with more useful
  -- expressions. A PBS is still presumed to follow the cosine-squared
  -- rule, as in Simulation A. Also, the probability of either
  -- polarization filter setting is one half, resulting in a factor of
  -- ¼ in the joint probability. (This can be proved with Bayes’
  -- theorem.) Therefore:

  --    P(hL1,vR1) = 2 × P(hL1) × P(vR1) = ½ × cos²(θ_L) × sin²(θ_R)
  --    P(hL1,vR2) = 2 × P(hL1) × P(vR2) = ½ × cos²(θ_L) × cos²(θ_R)
  --    P(hL2,vR1) = 2 × P(hL2) × P(vR1) = ½ × sin²(θ_L) × sin²(θ_R)
  --    P(hL2,vR2) = 2 × P(hL2) × P(vR1) = ½ × sin²(θ_L) × cos²(θ_R)
  --    P(vL1,hR1) = 2 × P(vL1) × P(hR1) = ½ × sin²(θ_L) × cos²(θ_R)
  --    P(vL1,hR2) = 2 × P(vL1) × P(hR2) = ½ × sin²(θ_L) × sin²(θ_R)
  --    P(vL2,hR1) = 2 × P(vL2) × P(hR1) = ½ × cos²(θ_L) × cos²(θ_R)
  --    P(vL2,hR2) = 2 × P(vL2) × P(vR1) = ½ × cos²(θ_L) × sin²(θ_R)

  -- We should expect certain symmetries:

  --    P(hL1,vR2) = P(vL2,hR1) = ½ × cos²(θ_L) × cos²(θ_R)
  --    P(vL1,hR2) = P(hL2,vR1) = ½ × sin²(θ_L) × sin²(θ_R)
  --    P(hL2,vR2) = P(vL1,hR1) = ½ × sin²(θ_L) × cos²(θ_R)
  --    P(hL1,vR1) = P(vL2,hR2) = ½ × cos²(θ_L) × sin²(θ_R)

  -- Plugging in some numbers for θ_L and θ_R confirms that the sum of
  -- these joint probabilities is one.

  -- A good test that have indeed characterized the supposed ‘quantum
  -- entanglement’ will be our ability to ‘clamp’ probabilities to
  -- zero, by adjusting the PBS settings. It can be deduced, from the
  -- (incorrect) assumption in [1] of ‘statistical independence of a
  -- and b’, that such clamping ought to be impossible without ‘action
  -- at a distance’.

  -- Set θ_L=0 and θ_R=π_2. Then

  --    P(hL1,vR2) = P(vL2,hR1) = ½ × cos²(θ_L) × cos²(θ_R) = 0
  --    P(vL1,hR2) = P(hL2,vR1) = ½ × sin²(θ_L) × sin²(θ_R) = 0
  --    P(hL2,vR2) = P(vL1,hR1) = ½ × sin²(θ_L) × cos²(θ_R) = 0
  --    P(hL1,vR1) = P(vL2,hR2) = ½ × cos²(θ_L) × sin²(θ_R) = ½

  -- Thus there is definitely the supposed ‘quantum entanglement’.

  -- The same methods could be used to solve for the probabilities of
  -- the outcomes in the experiment described in the Wikipedia
  -- article. No knowledge of ‘Pauli matrices’ or ‘tensor products of
  -- vector spaces’ is needed. One need not know what a ‘bra-ket’ is.

  -- Here, then, is Simulation B. It is merely random generation of
  -- possible event records. However, that is all one should expect of
  -- a ‘simulation’ of quantum theory.

  function one_event_B (θ_L, θ_R : scalar)
  return event_record is
    ev        : event_record;
    P_hL1_vR2 : scalar;
    P_vL2_hR1 : scalar;
    P_vL1_hR2 : scalar;
    P_hL2_vR1 : scalar;
    P_hL2_vR2 : scalar;
    P_vL1_hR1 : scalar;
    P_hL1_vR1 : scalar;
    P_vL2_hR2 : scalar;
    r1, r2    : scalar;
    r3, r4    : scalar;
    r5, r6    : scalar;
    r7, r8    : scalar;
    r         : scalar;
  begin
    ev.identifier := current_identifier;
    current_identifier := current_identifier + 1;

    ev.θ_left := θ_L;
    ev.θ_right := θ_R;

    P_hL1_vR1 := 0.5 * cos2 (θ_L) * sin2 (θ_R);
    P_hL1_vR2 := 0.5 * cos2 (θ_L) * cos2 (θ_R);
    P_hL2_vR1 := 0.5 * sin2 (θ_L) * sin2 (θ_R);
    P_hL2_vR2 := 0.5 * sin2 (θ_L) * cos2 (θ_R);
    P_vL1_hR1 := 0.5 * sin2 (θ_L) * cos2 (θ_R);
    P_vL1_hR2 := 0.5 * sin2 (θ_L) * sin2 (θ_R);
    P_vL2_hR1 := 0.5 * cos2 (θ_L) * cos2 (θ_R);
    P_vL2_hR2 := 0.5 * cos2 (θ_L) * sin2 (θ_R);

    r1 := P_hL1_vR1;
    r2 := r1 + P_hL1_vR2;
    r3 := r2 + P_hL2_vR1;
    r4 := r3 + P_hL2_vR2;
    r5 := r4 + P_vL1_hR1;
    r6 := r5 + P_vL1_hR2;
    r7 := r6 + P_vL2_hR1;
    r8 := r7 + P_vL2_hR2;

    -- Sanity check.
    if abs (r8 - 1.0) > 50.0 * scalar'model_epsilon then
      raise programmer_mistake;
    end if;

    ev.detector_left_1 := false;
    ev.detector_left_2 := false;
    ev.detector_right_1 := false;
    ev.detector_right_2 := false;

    r := random_scalar;
    if r < r1 then
      -- P_hL1_vR1
      ev.φ_left := horizontal;
      ev.detector_left_1 := true;
      ev.detector_right_1 := true;
    elsif r < r2 then
    -- P_hL1_vR2
      ev.φ_left := horizontal;
      ev.detector_left_1 := true;
      ev.detector_right_2 := true;
    elsif r < r3 then
    -- P_hL2_vR1
      ev.φ_left := horizontal;
      ev.detector_left_2 := true;
      ev.detector_right_1 := true;
    elsif r < r4 then
    -- P_hL2_vR2
      ev.φ_left := horizontal;
      ev.detector_left_2 := true;
      ev.detector_right_2 := true;
    elsif r < r5 then
    -- P_vL1_hR1
      ev.φ_left := vertical;
      ev.detector_left_1 := true;
      ev.detector_right_1 := true;
    elsif r < r6 then
    -- P_vL1_hR2
      ev.φ_left := vertical;
      ev.detector_left_1 := true;
      ev.detector_right_2 := true;
    elsif r < r7 then
    -- P_vL2_hR1
      ev.φ_left := vertical;
      ev.detector_left_2 := true;
      ev.detector_right_1 := true;
    else
    -- P_vL2_hR2
      ev.φ_left := vertical;
      ev.detector_left_2 := true;
      ev.detector_right_2 := true;
    end if;

    if ev.φ_left = vertical then
      ev.φ_right := horizontal;
    else
      ev.φ_right := vertical;
    end if;

    return ev;
  end one_event_B;

  function simulate_events_B (number_of_events : natural;
                              θ_L, θ_R         : scalar)
  return event_record_sets.set is
    use event_record_sets;
    events : set;
  begin
    for i in 1 .. number_of_events loop
      include (events, one_event_B (θ_L, θ_R));
    end loop;
    return events;
  end simulate_events_B;

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

  -- The following function calculates the frequencies of detections
  -- in a set of event records, so these frequencies can be compared
  -- to the probabilities.

  type detections_array is
    array (φ_value,     -- Setting of the left-side polarizing filter.
           range_1_to_2, -- Which detector was activated on the left.
           range_1_to_2) -- Which detector was activated on the right.
           of scalar;

  function detections_frequencies (events : event_record_sets.set)
  return detections_array is
    use event_record_sets;

    a    : detections_array;

    procedure initialize is
    begin
      for i in φ_value loop
        for j in range_1_to_2 loop
          for k in range_1_to_2 loop
            a(i, j, k) := 0.0;
          end loop;
        end loop;
      end loop;
    end initialize;

    procedure count is
      curs : cursor;
      ev   : event_record;

      procedure check_record is
      begin
        if ev.φ_left = horizontal and ev.φ_right = horizontal then
          raise programmer_mistake;
        end if;
        if ev.φ_left = vertical and ev.φ_right = vertical then
          raise programmer_mistake;
        end if;
        if not (ev.detector_left_1 xor ev.detector_left_2) then
          raise programmer_mistake;
        end if;
        if not (ev.detector_right_1 xor ev.detector_right_2) then
          raise programmer_mistake;
        end if;
      end check_record;

    begin
      curs := first (events);
      while has_element (curs) loop
        ev := element (curs);
        check_record;
        if ev.detector_left_1 then
          if ev.detector_right_1 then
            a(ev.φ_left, 1, 1) := @ + 1.0;
          else
            a(ev.φ_left, 1, 2) := @ + 1.0;
          end if;
        else
          if ev.detector_right_1 then
            a(ev.φ_left, 2, 1) := @ + 1.0;
          else
            a(ev.φ_left, 2, 2) := @ + 1.0;
          end if;
        end if;
        curs := next (curs);
      end loop;
    end count;

    procedure normalize is
      n : scalar;
    begin
      n := scalar (length (events));
      for i in φ_value loop
        for j in range_1_to_2 loop
          for k in range_1_to_2 loop
            a(i, j, k) := a(i, j, k) / n;
          end loop;
        end loop;
      end loop;
    end normalize;

  begin
    initialize;
    count;
    normalize;
    return a;
  end detections_frequencies;

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

  -- The following function fills in a detections_array with the
  -- probabilities themselves, rather than simulated frequencies.

  function detections_probabilities (θ_L, θ_R : scalar)
  return detections_array is
    a : detections_array;
  begin
    a(horizontal, 1, 1) := 0.5 * cos2 (θ_L) * sin2 (θ_R);
    a(horizontal, 1, 2) := 0.5 * cos2 (θ_L) * cos2 (θ_R);
    a(horizontal, 2, 1) := 0.5 * sin2 (θ_L) * sin2 (θ_R);
    a(horizontal, 2, 2) := 0.5 * sin2 (θ_L) * cos2 (θ_R);
    a(vertical, 1, 1) := 0.5 * sin2 (θ_L) * cos2 (θ_R);
    a(vertical, 1, 2) := 0.5 * sin2 (θ_L) * sin2 (θ_R);
    a(vertical, 2, 1) := 0.5 * cos2 (θ_L) * cos2 (θ_R);
    a(vertical, 2, 2) := 0.5 * cos2 (θ_L) * sin2 (θ_R);
    return a;
  end detections_probabilities;

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

begin

  -- The main program.

  declare
    number_of_events   : natural;
    θ_L, θ_R           : scalar;
    events_A, events_B : event_record_sets.set;
    probabilities      : detections_array;
    frequencies_A      : detections_array;
    frequencies_B      : detections_array;
    φ                  : φ_value;
    i, j               : range_1_to_2;
  begin
    seed := 1234;
    number_of_events := 1000000;

    for i in 1 .. 4 loop
      if i = 1 then
        θ_L := π / 4.0;
        θ_R := 2.0 * π / 3.0;
      else
        θ_L := π_2 * random_scalar;
        θ_R := π_2 * random_scalar;
      end if;

      events_A := simulate_events_A (number_of_events, θ_L, θ_R);
      events_B := simulate_events_B (number_of_events, θ_L, θ_R);

      probabilities := detections_probabilities (θ_L, θ_R);
      frequencies_A := detections_frequencies (events_A);
      frequencies_B := detections_frequencies (events_B);

      new_line;

      set_col (2);
      put ("PBS left  θ");
      set_col (15);
      put (θ_L / π_180, 3, 5, 0);
      new_line;

      set_col (2);
      put ("PBS right θ");
      set_col (15);
      put (θ_R / π_180, 3, 5, 0);
      new_line;

      set_col (5);
      put ("φ");

      set_col (10);
      put ("detec L");

      set_col (20);
      put ("detec R");

      set_col (30);
      put ("probability");

      set_col (44);
      put ("freq classical");

      set_col (60);
      put ("freq quantum");

      for φ in φ_value loop
        for i in range_1_to_2 loop
          for j in range_1_to_2 loop

            set_col (4);
            if φ in horizontal then
              put ("H/V");
            else
              put ("V/H");
            end if;

            set_col (13);
            put (i, 1);

            set_col (23);
            put (j, 1);

            set_col (30);
            put (probabilities(φ, i, j), 3, 5, 0);

            set_col (45);
            put (frequencies_A(φ, i, j), 3, 5, 0);

            set_col (60);
            put (frequencies_B(φ, i, j), 3, 5, 0);

          end loop;
        end loop;
      end loop;

      new_line;

    end loop;

    new_line;
  end;

end eprb_simulation;

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

-- The output of the program should look something like the following,
-- demonstrating that both the classical and quantum simulations
-- approximate the quantum predictions.


--  PBS left  θ   45.00000
--  PBS right θ  120.00000
--     φ    detec L   detec R   probability   freq classical  freq quantum
--    H/V      1         1        0.18750        0.18709        0.18705
--    H/V      1         2        0.06250        0.06252        0.06234
--    H/V      2         1        0.18750        0.18646        0.18788
--    H/V      2         2        0.06250        0.06223        0.06256
--    V/H      1         1        0.06250        0.06263        0.06275
--    V/H      1         2        0.18750        0.18838        0.18733
--    V/H      2         1        0.06250        0.06258        0.06256
--    V/H      2         2        0.18750        0.18809        0.18752

--  PBS left  θ   35.40353
--  PBS right θ   58.09312
--     φ    detec L   detec R   probability   freq classical  freq quantum
--    H/V      1         1        0.23939        0.23912        0.23933
--    H/V      1         2        0.09280        0.09297        0.09271
--    H/V      2         1        0.12093        0.12105        0.12094
--    H/V      2         2        0.04688        0.04711        0.04710
--    V/H      1         1        0.04688        0.04661        0.04698
--    V/H      1         2        0.12093        0.12076        0.12101
--    V/H      2         1        0.09280        0.09265        0.09259
--    V/H      2         2        0.23939        0.23973        0.23933

--  PBS left  θ   25.34591
--  PBS right θ   89.05482
--     φ    detec L   detec R   probability   freq classical  freq quantum
--    H/V      1         1        0.40826        0.40835        0.40873
--    H/V      1         2        0.00011        0.00011        0.00010
--    H/V      2         1        0.09160        0.09127        0.09122
--    H/V      2         2        0.00002        0.00003        0.00002
--    V/H      1         1        0.00002        0.00002        0.00002
--    V/H      1         2        0.09160        0.09147        0.09224
--    V/H      2         1        0.00011        0.00011        0.00010
--    V/H      2         2        0.40826        0.40864        0.40758

--  PBS left  θ   75.33916
--  PBS right θ   76.26485
--     φ    detec L   detec R   probability   freq classical  freq quantum
--    H/V      1         1        0.03022        0.03033        0.03039
--    H/V      1         2        0.00181        0.00183        0.00181
--    H/V      2         1        0.44159        0.44165        0.44100
--    H/V      2         2        0.02638        0.02656        0.02618
--    V/H      1         1        0.02638        0.02653        0.02636
--    V/H      1         2        0.44159        0.44117        0.44220
--    V/H      2         1        0.00181        0.00178        0.00180
--    V/H      2         2        0.03022        0.03015        0.03026

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

-- Afterword:

-- Of course, Simulation A also is soluble in closed form by
-- probability theory. And you will find that its closed form solution
-- is exactly the same as that of Simulation B. As, of course, it must
-- be.

-- The problem with ‘Bell’s theorem’ is that Bell obviously did not
-- know any probability theory. He only pretended to know
-- probability theory, and then he solved our Simulation A, and
-- every problem in its class, incorrectly. He did not know how to
-- do scientific inference, and he found an audience that also did
-- not know how to do it.

-- Perhaps it is the case to this day that quantum physicists do not
-- realize ‘quantum’ problems can be solved without quantum
-- theory. They believe in the magic of their incantations. In any
-- case, methods of scientific inference are not taught to quantum
-- physics students. What can be solved by scientific inference
-- probably involves cause and effect by contact action, not
-- ‘entanglement’ and ‘non-locality’. Restricting the teaching to the
-- special incantations of quantum theory, and to the set of magic
-- doctrines that accompany the incantations, neatly does its
-- work. The students come out of their lessons believing in magic and
-- telling the public that it is ‘science’.

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

-- A footnote on ‘Bayes’ theorem’:

-- We should ignore all the stuff on the ‘Bayes’ theorem’ Wikipedia
-- page about ‘Interpretations’. Such disputes are a topic separate
-- from what we are dealing with here. What matters to us is that
-- probability theory done properly is FORMAL, not merely a batch of
-- intuitions.

-- By far most specialists in mathematics regard probability theory as
-- a subset of ‘measure theory’, which is a kind of generalization of
-- the notion of areas. Problems concerning various sets of interest
-- to mathematicians tend to be discussed.

-- Alternatively, probability theory can be built up as a very
-- carefully devised generalization of boolean logic, specifically
-- designed for the task of logical inference. This form of
-- probability theory is probably more useful to empirical scientists
-- and engineers than is the other.

-- Either form of probability works for us, because in either form of
-- probability theory the theorems that we use are the same.

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