#pragma mode( separator(.,;) integer(h32) )

// (C) 2026 Roger D Metcalf
// Program concept and rudimentary code by Roger D. Metcalf
// Code corrected, expanded, and vastly improved by Claude Code Fable 5
// For academic use only – do not rely on results for anything important
// Results have not been rigorously verified and/or validated
//
// FREE FOR EDUCATIONAL AND PERSONAL USE.  Keep this header with the
//   file and credit the author if you pass it on.  Not for sale or
//   commercial redistribution.  Written with AI assistance (Anthropic
//   Claude); concept, direction and all hardware testing by the author.
// DISCLAIMER: a teaching and demonstration aid. Not validated for
//   forensic, clinical or legal use; do not rely on its results for
//   legal purposes.


// Sampling Distribution Calculator — Roger D Metcalf
//   Sampling distribution of the sample mean, x-bar ~ Normal(μ, σ/√n).
// Launch:  SampDist()   (zero-argument, menu-first — BayesTree pattern)
//   The Prime numbers the menu items itself:
//   1  Find probability      left / right / between / outside
//   2  Find critical value   lower / upper / two-tailed
//   3  Compare sample sizes  three curves, shared axes
//   4  Watch sample means    40 simulated means fall onto the axis
//   5  Help
// Exports EP_SD_* so results are readable in Home after the menu is cancelled.
//
// Written for the 2026-09-09 firmware (HP Prime G2 and G1).  v2.3b has
//   NOT yet been run on hardware or on the emulator, only in the PC
//   harness (ppl_run_v3e, 19 Oct 2026); v2.3 passed the emulator's round
//   (17 Oct 2026).  See the checklist at the end, SampDist_v2.3b_CHANGES.md,
//   SampDist_v2.3a_CHANGES.md and SampDist_v2.3_CHANGES.md.
//
// v2.3b (19 Oct 2026): S1, S2 (fleet review 19 Oct 2026, findings_B).
//   Line numbers below are v2.3a's.
//   * S1: in critical-value mode α is read back from the dragged line and
//     kept off 0 and 1 for NORMALD_ICDF (sd_fix, 594).  The floor was
//     0.000001, which the panel's six decimals print as an exact
//     "= 0.000001": a red line dragged to the right edge read
//     "P(x-bar > 144.000) = 0.000001" while its own z line said
//     "< 0.001%" (α at z = 16 is about 1E-58).  The floor is now 1E-7
//     (and the ceiling 1 - 1E-7), below what six decimals show, so the
//     panel reads "< 0.000001".
//   * S2: sd_fix's right tail (592) was the one 1 - NORMALD_CDF the v2.3a
//     pass did not turn round; it is now NORMALD_CDF(μ, σx-bar, 2μ - a),
//     as in sd_prob and Mode 1.  Nothing on the panel changes (the clamp
//     hid the difference); the code is now the same everywhere.
//   Nothing else changed.
//
// v2.3a (17 Oct 2026, evening): no "=" before a relation, and the right
//   tail kept.  Roger approved the pass, 17 Oct 2026, 21:18.  Line
//   numbers below are v2.3's.
//   * sd_fmtp writes "< 0.000001" below its floor, and the result panel
//     put that after "= ": "P(x-bar > 1.000) = < 0.000001" and, on the z line,
//     "= < 0.001%".  New sd_Eq (after sd_fmtp) returns " = 0.158655", or
//     " < 0.000001" when the relation is already there; the panel's four
//     probability lines and its percent use it (641, 645, 650, 655, 658).
//   * The right tail is the lower tail NORMALD_CDF(μ, σ(x-bar), 2μ - a), the
//     same number by symmetry, in sd_prob (509, 514) and Mode 1 (1229,
//     1234), for "Right tail" and the right half of "Outside".  1 -
//     NORMALD_CDF is exactly 0 once the tail is below the CDF's last digit,
//     so a far right tail read "= 0.000000" where the same tail on the left
//     read "< 0.000001".
//   Nothing else changed.
//
// v2.3 (16 Oct 2026): the touch and graphics version, built on v2.2.
//   Line numbers below are v2.2's.
//   * The parent population N(μ, σ) is drawn as a faint grey curve
//     behind x-bar's curve on the probability and critical-value
//     graphs, labelled "population (σ = 15)" near its right shoulder.
//     Both curves share ONE y-scale, fitted to x-bar's curve, so the
//     population's peak is 1/√n of x-bar's (18 % at n = 30).  The
//     window is fitted once, on entry, to show both (μ ± 16σx-bar,
//     but never wider than μ ± 3.5σ or narrower than μ ± 4σx-bar),
//     and then FROZEN, so changing n moves x-bar's curve and leaves
//     the population where it is (Lessons Learned section 6).
//   * sd_drawDist (227-387) is split into a compose step
//     (sd_compose) and a live loop (sd_live).  The result graph now
//     carries soft keys  n−  n+  n×4  μ,σ  Tail  Menu :  n− / n+ step
//     n through 1 2 4 5 10 16 20 25 30 40 50 64 100 160 200 400,
//     n×4 multiplies n by 4 (up to 1600) so x-bar's curve snaps to
//     half its width, μ,σ reopens the form, Tail cycles the tail
//     type; the probability, the σx-bar line and the shading are
//     recomputed after each.  In critical-value mode a new n moves
//     the critical values (α is kept); in probability mode the
//     bounds stay and P changes.
//   * Live bound: a finger in the plot box moves the NEARER bound and
//     the probability follows it (two NORMALD_CDF calls per frame;
//     the curves are drawn from cached pixel heights).  On lift the
//     bound snaps to a step of σx-bar/4 rounded to 1, 2 or 5.  Left /
//     Right nudge the bound by the same step, Up / Down choose which
//     bound in a two-bound region.  In critical-value mode dragging
//     the critical value changes α (both values move together in the
//     two-tailed case).  + / − = n+ / n−, Enter = a menu with the
//     other soft-key actions, Esc = back to the menu.
//   * The footer text moved into the result panel (the soft-key
//     strip occupies y >= 218); the plot box is now x 40-314,
//     y 20-150, the panel y 166-216.
//   * y-axis density labels in grey on the result graph.
//   * Compare sample sizes (sd_drawCmp, 392-464): the three curves
//     are filled (widest first, light red / blue / green) and then
//     outlined; on opening, n grows from 1 to its value over 8
//     frames (any key or tap skips); a finger on the plot is a
//     vertical cursor showing P(x-bar < x) for each n in its colour;
//     a tap on a curve picks it (drawn thicker, n and σx-bar
//     printed).  The static picture is composed once into G2 and
//     copied into G1 per cursor frame (GaltonBoard's idiom).
//   * New menu item "Watch sample means": 40 samples of n draws are
//     simulated (Marsaglia polar normals, PPL, sd_Norm), each mean
//     drops as a dot onto the axis under x-bar's curve, then a
//     histogram of the 40 means (bars one σx-bar wide, area 1) is
//     drawn under the curve with their mean and SD beside μ and
//     σx-bar.  Any key or tap skips the drop.
//   * Help is two pages and documents the soft keys and gestures;
//     pages are held by sd_WaitAny (a key OR a tap).
//   * New tagged helpers: sd_Drain, sd_WaitAny, sd_Norm, sd_pdfList,
//     sd_pyList, sd_poly, sd_fill, sd_quick, sd_grid, sd_curves,
//     sd_axes, sd_prob, sd_crit, sd_fix, sd_compose, sd_live,
//     sd_watch.  sd_drawOne and sd_drawDist are gone.
//   * The numbers are unchanged: sd_prob and sd_crit hold the same
//     NORMALD_CDF / NORMALD_ICDF formulas as modes 1 and 2.
//
// v2.2 (16 Oct 2026): fixes from the fleet review of 16 Oct 2026.
//   Line numbers below are v2.1's.
//   * Drain-wait-drain round all three graph holds (sd_drawDist 329,
//     sd_drawCmp 400, sd_showHelp 574).  Each was a bare
//     `UNTIL kk > -1` straight after INPUT or CHOOSE, so the OK press
//     that closed the form could satisfy it and the graph vanished
//     before it was seen (Lessons Learned section 7).
//   * The comparison title printed σ as semax·√n₁ (385), which is
//     only σ when n₁ is the SMALLEST n: n = 100, 30, 5 showed
//     "σ = 67.08" for σ = 15.  sd_drawCmp now receives σ itself.
//   * Menu strings no longer carry their own numbers (588-591): the
//     Prime numbers CHOOSE items, so each line showed two.  The help
//     page headings (556, 560, 564) are the item names, unnumbered.
//   * Compare mode (526) now rejects a non-integer n like modes 1 and
//     2 do (418, 481); the v2.0 note above claimed it already did.
//   * sd_lblChars (151-166) is gone.  It existed because SIZE() on a
//     string was unverified; NormExplorer's tick labels use
//     3*SIZE(label) and have run on the G2 and G1, and DIM on strings
//     is on the confirmed list.  Tick labels now centre on
//     3*SIZE(label) directly.
//   * One footer wording on all three screens, in the house form and
//     colour: "any key = menu", dodger blue RGB(30,144,255), at y = 226.
//   * The 7-name LOCAL line (198) is split; its bare x, y, i are now
//     ddX, ddY, ddI (a bare x is the CAS's symbolic variable).
//   * NORMALD is evaluated once per pixel column into a list (ddPdf);
//     the shading loop (240) and the curve loop (252) both read it
//     instead of each calling NORMALD at the same x.
//   * The result panel gains the line "σ = 15.00, σx-bar = 15.00/√30 =
//     2.7386" — σ against σx-bar is THE sampling-distribution
//     misconception, so the two now sit side by side with the
//     arithmetic between them.  sd_drawDist receives σ for it.
//   * Compare-mode defaults 5 / 30 / 100 (519) become 10 / 40 / 160,
//     so the "4× n halves σx-bar" line under the plot describes the
//     picture it sits under.
//   * sd_fmtp's `pv != 0` is written `pv <> 0` (corpus-confirmed form,
//     Lessons Learned section 3) so the lint gate's E8 check passes;
//     `!=` has also run on hardware (EpiStats), so this is cosmetic.
//
// v2.1 (10 Oct 2026): every helper routine now carries the program
//   tag "sd_" (fmtp -> sd_fmtp, drawDist -> sd_drawDist).
//   The Prime keeps the routines of ALL programs in one shared
//   pool, so two programs that define the same helper name
//   (NiceNum, fmtp and lngamma were each shared by up to eight
//   programs) can run each other's code.  That is what sent
//   StudRange into a Burr3 screen and forced a reset.
//   Rule: every non-EXPORT routine gets the program's tag.
// v2.0 CHANGES (29 Jul 2026) — three latent syntax errors, four silent
// wrong-answer paths, and the fleet-wide formatting pass.
//
//  SYNTAX (would not compile, or would fail the moment the branch ran)
//   * `col` was a PARAMETER of DrawCurve — `col` is a built-in name, and the
//     parameter parser is stricter than the LOCAL parser (the `fp` lesson).
//     Renamed clr. This is the one that would have bitten first.
//   * SEVEN `LOCAL` declarations sat inside IF blocks (tmp in Mode_FindProb,
//     msg ×3 and tailp in Mode_CriticalVal, pxa and pxb in DrawGraph). LOCAL
//     inside IF/CASE is rejected; all hoisted to the top of their function.
//   * DrawGraph declared `LOCAL stp, gmin, gmax, xval, tickx, nfrac` halfway
//     down the body, after executable statements. Hoisted.
//
//  SILENT WRONG ANSWERS (compiled fine, produced a plausible number)
//   * All five INPUT() returns were discarded. Cancelling the bounds form in
//     Mode 1 left a := μ−SE and b := μ+SE — the pre-seeded defaults — so the
//     program drew a graph and reported P = 0.682689 as though the user had
//     typed those bounds. Every INPUT now tests its return.
//   * The three CHOOSE index variables were never initialised to 0, and none
//     was reset before its call. CHOOSE does not touch its index on CANCEL,
//     so a cancelled menu silently re-used the previous answer. modeM is now
//     reset INSIDE the menu loop, not once before it (the FDAES lesson).
//   * `INPUT(prob, "…", "α:", "…")` used the scalar form; the first argument
//     must be a list of variable NAMES. Now {prob} with list labels/help.
//   * n was never checked for integrality. n = 2.5 gave σ/√2.5, a standard
//     error for a sample size that cannot exist. Now rejected.
//
//  DISPLAY
//   * 33 of the 41 STRING() calls were mode 1 (Standard), which drops
//     trailing zeros regardless of the digit count: α = 0.05 printed "5%"
//     not "5.00%", and a p of exactly 0.5 printed "0.5" where six decimals
//     were asked for. All converted to fixed STRING(x,2,d); the 7 genuine
//     integer sites keep STRING(n,1,0).
//   * Probabilities now route through sd_fmtp() — the same helper as FEP /
//     StatDist / Norm3 — so a far-tail p prints "< 0.000001" instead of
//     "0.000000", which reads as exactly zero.
//   * x-bar was written as X + U+0304 COMBINING MACRON in 14 display strings.
//     The macron does not stack on the device; x-bar is CHAR(57344).
//   * The readout block (μ, σ of the mean, n, P) was drawn at x=5, y=45..90
//     — INSIDE the plot, whose left edge is x=30 and top is y=40. The SE line
//     is 11 characters ≈ 66 px, so it ran to x≈71 and sat on top of the curve
//     and the shading. Results moved to an orange panel below the axis.
//   * x tick labels were centred with a flat `tickx-9` regardless of label
//     width; now centred on the actual character count (sd_lblChars in
//     v2.0-2.1; 3*SIZE(label) from v2.2).
//   * No double-buffering: both graphs were built in view on G0. Now drawn
//     to G1 and BLIT once.
//   * Polyline gap bug: when a point fell outside the y-window the code
//     skipped it WITHOUT updating lastpy, but left firstpoint at 0 — so the
//     next in-range point drew a segment from a stale y at px-1. Points are
//     now clamped, not skipped.
//   * 19 EXPORTs, including six colour globals (BLACK/WHITE/RED/BLUE/
//     LTBLUE/GRAY) squatting on names any other program might use. Only the
//     entry point and the six EP_SD_* results are exported now; the eleven
//     helpers are LOCAL, so they are defined before first use.
//   * `WAIT;` replaced with the device-verified REPEAT/GETKEY hold.
//   * Colours moved to Roger's house palette.
//
//  NEW
//   * Tail-first flow: the tail type is chosen BEFORE the bounds, so the
//     form asks for one boundary or two as appropriate. v1 inferred the tail
//     from ABS(b-a) < 0.0001, which meant "one-tailed" was expressed by
//     typing the same number twice.
//   * Fourth region: Outside, P(x-bar < a or x-bar > b) — the two-tailed rejection
//     region, which Mode 1 could not express at all.
//   * Menu loops instead of exiting after one calculation.
//
// VERIFIED (verify_SampDist.py — 308 checks, 0 failures; v2.3 adds
//   verify_SampDist_v2.3.py, which runs this file's own routines in
//   hp-prime-kit's interpreter against SciPy; v2.3a: verify_SampDist_v2.3a.py;
//   v2.3b: verify_SampDist_v2.3b.py):
//   SE = σ/√n; all four region probabilities against scipy AND against an
//   independent quadrature of the pdf; critical values round-tripped through
//   the CDF; the complement identity P(between) + P(outside) = 1; the
//   two-tailed region symmetric about μ; the μ±4σ window covering 99.9937%;
//   the tick ladder landing inside a negative range; and (v2.0-2.1 only)
//   sd_lblChars matching the true printed width of STRING(x,2,d).
//
// NOT DEVICE-VERIFIED — see the checklist at the end of this file.

EXPORT EP_SD_SE;
EXPORT EP_SD_P;
EXPORT EP_SD_ZLO;
EXPORT EP_SD_ZHI;
EXPORT EP_SD_LO;
EXPORT EP_SD_HI;

// ---- fixed-decimal formatter with a floor ----
// Identical to the helper in FEP / StatDist / BayesTree / Norm3. Mode 2 keeps
// trailing zeros; the floor stops a small tail probability printing as a row
// of zeros. Params pv/dg, not p/d — short parameter names collide with
// builtins. Non-EXPORT so it stays out of the catalog.
sd_fmtp(pv, dg)
BEGIN
  LOCAL thr;
  thr := 10^(-dg);
  IF pv <> 0 AND ABS(pv) < thr THEN
    IF pv > 0 THEN RETURN "< " + STRING(thr,2,dg); END;
    RETURN "> -" + STRING(thr,2,dg);
  END;
  RETURN STRING(pv,2,dg);
END;

// " = 0.158655", or " < 0.000001" when sd_fmtp has already written the
// relation (v2.3a): "P(x-bar > a) = < 0.000001" said two things at once;
// "P(x-bar > a) < 0.000001" says the one that is true.  The first character
// decides, read only when there is one (DIM first: an empty string has
// no first character).
sd_Eq(eqS)
BEGIN
  IF DIM(eqS) > 0 THEN
    IF LEFT(eqS, 1) == "<" OR LEFT(eqS, 1) == ">" THEN RETURN " " + eqS; END;
  END;
  RETURN " = " + eqS;
END;

// ---- "nice" number for axis ticks (Heckbert) ----
sd_NiceNum(x, roundf)
BEGIN
  LOCAL expv, f, nf;
  IF x <= 0 THEN RETURN 1; END;
  expv := FLOOR(LOG(x));
  f := x / 10^expv;
  IF roundf THEN
    CASE
      IF f < 1.5 THEN nf := 1; END;
      IF f < 3 THEN nf := 2; END;
      IF f < 7 THEN nf := 5; END;
      DEFAULT nf := 10;
    END;
  ELSE
    CASE
      IF f <= 1 THEN nf := 1; END;
      IF f <= 2 THEN nf := 2; END;
      IF f <= 5 THEN nf := 5; END;
      DEFAULT nf := 10;
    END;
  END;
  RETURN nf * 10^expv;
END;

// (v2.2: sd_lblChars, which worked out a tick label's printed width from
// its value because SIZE() on a string was unverified, is gone.  SIZE on a
// string has since run on the G2 and G1 in NormExplorer's tick labels, so
// the labels below centre on 3*SIZE(label) directly.)

// ---- the house touch kit (ChiTable beta 1k's csDrain / csWaitAny) ----
// Wait for the finger to lift -- but never for more than a second.
// MOUSE(1) is -1 while nothing touches the screen.  Without the time
// limit a stuck or phantom touch would hold the program here for ever,
// and the only way out would be a reset.
sd_Drain()
BEGIN
  LOCAL drT;
  drT := TICKS;                                  // milliseconds
  REPEAT UNTIL MOUSE(1) == -1 OR TICKS - drT > 1000;
END;

// Wait for any key OR a tap (drain, wait, drain).  Returns the key code,
// or -1 for a tap.  A GETKEY drain alone does not drain touch, so the
// finger that tapped to open a page could otherwise close it at once.
sd_WaitAny()
BEGIN
  LOCAL wk, wm, wm1;
  REPEAT wk := GETKEY; UNTIL wk == -1;
  sd_Drain();
  REPEAT
    wk := GETKEY;
    wm := MOUSE;
    wm1 := wm(1);
  UNTIL wk > -1 OR SIZE(wm1) > 0;
  IF SIZE(wm1) > 0 THEN wk := -1; END;
  sd_Drain();
  REPEAT UNTIL GETKEY == -1;
  RETURN wk;
END;

// ---- one standard normal random number, Marsaglia's polar method ----
// Pick a random point in the square [-1,1] x [-1,1], keep it only if it
// falls inside the unit circle, then transform.  No sine or cosine, so
// the angle mode (degrees or radians) cannot affect it (Lessons Learned
// section 11; CentralLimit's cltNorm is the same routine).
sd_Norm()
BEGIN
  LOCAL nv1, nv2, nsq;
  REPEAT
    nv1 := 2*RANDOM - 1;
    nv2 := 2*RANDOM - 1;
    nsq := nv1*nv1 + nv2*nv2;
  UNTIL nsq > 0 AND nsq < 1;
  RETURN nv1 * (-2*LN(nsq)/nsq)^0.5;
END;

// ---- the density ONCE per pixel column ----
// Returns a list of nc+1 values: element j+1 is the Normal(mu, sgm) density
// at the x of pixel column j (0 .. nc) of a window xlo .. xhi.  Every
// curve, shading pass and hit test reads this list (or the pixel heights
// made from it) instead of calling NORMALD again (Lessons Learned
// section 6: cache the density for EVERY curve).
sd_pdfList(mu, sgm, xlo, xhi, nc)
BEGIN
  LOCAL plV, plJ;
  plV := MAKELIST(0, K, 0, nc, 1);
  FOR plJ FROM 0 TO nc DO
    plV(plJ + 1) := NORMALD(mu, sgm, xlo + plJ/nc*(xhi - xlo));
  END;
  RETURN plV;
END;

// ---- densities -> screen rows ----
// The y-flip: row = ybot - density/ymx*yh, clamped at the top of the box
// (clamped, never skipped, so a polyline cannot jump from a stale y).
sd_pyList(vals, ymx, ybot, ytop, yh)
BEGIN
  LOCAL pyL, pyJ, pyN, pyY;
  pyN := SIZE(vals);
  pyL := MAKELIST(0, K, 1, pyN, 1);
  FOR pyJ FROM 1 TO pyN DO
    pyY := IP(ybot - vals(pyJ)/ymx*yh);
    IF pyY < ytop THEN pyY := ytop; END;
    pyL(pyJ) := pyY;
  END;
  RETURN pyL;
END;

// ---- a curve from cached rows ----
// Element j of yl is the row at column lft + j - 1.  A segment lying flat
// on the axis (both ends at row ybot) is skipped: the axis line drawn
// later covers it, and on a narrow curve that saves most of the calls.
// yo shifts the line up or down a pixel, to draw a thick curve.
sd_poly(poG, yl, lft, ybot, clr, yo)
BEGIN
  LOCAL poJ, poN;
  poN := SIZE(yl);
  FOR poJ FROM 2 TO poN DO
    IF yl(poJ - 1) < ybot OR yl(poJ) < ybot THEN
      LINE_P(poG, lft + poJ - 2, yl(poJ - 1) + yo, lft + poJ - 1, yl(poJ) + yo, clr);
    END;
  END;
END;

// ---- fill under a curve, one vertical line per pixel column ----
// (pixel columns, not math samples, so there are no white gap columns)
sd_fill(fiG, yl, lft, ybot, clr)
BEGIN
  LOCAL fiJ, fiN;
  fiN := SIZE(yl);
  FOR fiJ FROM 1 TO fiN DO
    IF yl(fiJ) < ybot THEN
      LINE_P(fiG, lft + fiJ - 1, ybot, lft + fiJ - 1, yl(fiJ), clr);
    END;
  END;
END;

// ---- a quick curve for animation frames: every 2nd column, uncached ----
// qg holds the window and the box: {xlo, xhi, ymx, LFT, RGT, TOPy, BOT}.
// Used only by the 8-frame opening of Compare, where each frame has a
// different n and a cache would be used once.
sd_quick(quG, mu, sgm, qg, clr)
BEGIN
  LOCAL quX, quY, quPX, quPY, quW, quH;
  quW := qg(5) - qg(4); quH := qg(7) - qg(6);
  quPX := -1;
  FOR quX FROM qg(4) TO qg(5) STEP 2 DO
    quY := IP(qg(7) - NORMALD(mu, sgm, qg(1) + (quX - qg(4))/quW*(qg(2) - qg(1)))/qg(3)*quH);
    IF quY < qg(6) THEN quY := qg(6); END;
    IF quPX >= 0 THEN LINE_P(quG, quPX, quPY, quX, quY, clr); END;
    quPX := quX; quPY := quY;
  END;
END;

// ============================================================
// THE RESULT GRAPH (probability, critical value, Watch)
// Plot box x 40..314, y 20..150 (W = 274, H = 130).  These constants are
// repeated in sd_grid, sd_curves, sd_axes, sd_compose, sd_live and
// sd_watch, so the drawing and the touch code cannot disagree; change
// them together.
// ============================================================

// white page and the vertical grid lines at the x ticks
sd_grid(grG, xlo, xhi)
BEGIN
  LOCAL LFT, RGT, TOPy, BOT, W;
  LOCAL stp, xv, tickx;
  LFT := 40; RGT := 314; TOPy := 20; BOT := 150; W := RGT - LFT;
  RECT_P(grG, 0, 0, 319, 239, RGB(255,255,255));
  stp := sd_NiceNum((xhi - xlo)/6, 1);
  xv := CEILING(xlo/stp)*stp;
  WHILE xv <= xhi + stp/1000 DO
    tickx := LFT + IP((xv - xlo)/(xhi - xlo)*W);
    IF tickx > LFT AND tickx <= RGT THEN
      LINE_P(grG, tickx, TOPy, tickx, BOT-1, RGB(225,228,235));
    END;
    xv := xv + stp;
  END;
END;

// the population (faint grey, with its label), x-bar's curve (blue) over
// it, and the dotted green mean marker
sd_curves(cuG, xlo, xhi, mu, sg, yX, yP)
BEGIN
  LOCAL LFT, RGT, TOPy, BOT, W;
  LOCAL cPop, cuL, cuX, cuY, cuI, tickx;
  LOCAL cuJ;
  LFT := 40; RGT := 314; TOPy := 20; BOT := 150; W := RGT - LFT;
  cPop := RGB(150,150,150);
  // the parent population N(μ, σ): same y-scale as x-bar's curve, so its
  // peak is 1/√n of x-bar's.  Drawn first, so x-bar's curve sits on top.
  // Every 2nd column: it is a faint guide line, and it spans the whole
  // window, so this halves the biggest cost of a drag frame.
  cuJ := 3;
  WHILE cuJ <= SIZE(yP) DO
    IF yP(cuJ - 2) < BOT OR yP(cuJ) < BOT THEN
      LINE_P(cuG, LFT + cuJ - 3, yP(cuJ - 2), LFT + cuJ - 1, yP(cuJ), cPop);
    END;
    cuJ := cuJ + 2;
  END;
  // its label, just above the curve at the right shoulder (μ + 1.2σ);
  // when the shoulder is off the window the label goes to the right
  // end, where the population is still visible as a low line
  cuL := "population (σ = " + IFTE(FP(sg) == 0, STRING(sg,2,0), STRING(sg,2,2)) + ")";
  cuX := mu + 1.2*sg;
  IF cuX > xlo AND cuX < xhi THEN
    cuI := LFT + IP((cuX - xlo)/(xhi - xlo)*W) + 3;
  ELSE
    cuI := RGT - 6*SIZE(cuL);
  END;
  IF cuI + 6*SIZE(cuL) > RGT THEN cuI := RGT - 6*SIZE(cuL); END;
  IF cuI < LFT + 2 THEN cuI := LFT + 2; END;
  cuY := yP(cuI - LFT + 1) - 12;
  IF cuY < TOPy + 1 THEN cuY := TOPy + 1; END;
  IF cuY > BOT - 12 THEN cuY := BOT - 12; END;
  TEXTOUT_P(cuL, cuG, cuI, cuY, 1, cPop);
  // x-bar's curve over everything else
  sd_poly(cuG, yX, LFT, BOT, RGB(0,60,200), 0);
  // dotted mean marker, so it reads differently from the solid red bounds
  tickx := LFT + IP((mu - xlo)/(xhi - xlo)*W);
  IF tickx > LFT AND tickx < RGT THEN
    FOR cuI FROM TOPy TO BOT STEP 4 DO
      PIXON_P(cuG, tickx, cuI, RGB(0,140,60));
    END;
  END;
END;

// axes; x ticks with labels centred BELOW the axis; y ticks with the
// density labels in grey to the LEFT of the axis (font 1 is about 6 px
// per character, so half a label is 3*SIZE and a whole one 6*SIZE)
sd_axes(axG, xlo, xhi, ymx)
BEGIN
  LOCAL LFT, RGT, TOPy, BOT, W, H;
  LOCAL stp, nfrac, xv, tickx, axL, axY;
  LFT := 40; RGT := 314; TOPy := 20; BOT := 150;
  W := RGT - LFT; H := BOT - TOPy;
  LINE_P(axG, LFT, BOT, RGT, BOT, RGB(0,0,0));
  LINE_P(axG, LFT, BOT, LFT, TOPy, RGB(0,0,0));
  stp := sd_NiceNum((xhi - xlo)/6, 1);
  nfrac := MAX(0, -FLOOR(LOG(stp)));
  xv := CEILING(xlo/stp)*stp;
  WHILE xv <= xhi + stp/1000 DO
    tickx := LFT + IP((xv - xlo)/(xhi - xlo)*W);
    IF tickx >= LFT AND tickx <= RGT THEN
      LINE_P(axG, tickx, BOT, tickx, BOT+3, RGB(0,0,0));
      axL := STRING(xv,2,nfrac);
      // centred on the tick, but never past the right edge of the screen
      TEXTOUT_P(axL, axG, MIN(tickx - 3*SIZE(axL), 318 - 6*SIZE(axL)), BOT+4, 1, RGB(90,90,90));
    END;
    xv := xv + stp;
  END;
  // density ticks: about four, on the nice 1/2/5 ladder; none in the top
  // 5 % of the box, where a label would touch the header
  stp := sd_NiceNum(ymx/4, 1);
  nfrac := MAX(0, -FLOOR(LOG(stp)));
  xv := stp;
  WHILE xv <= ymx*0.95 DO
    axY := BOT - IP(xv/ymx*H);
    LINE_P(axG, LFT - 3, axY, LFT, axY, RGB(0,0,0));
    axL := STRING(xv,2,nfrac);
    TEXTOUT_P(axL, axG, MAX(0, LFT - 5 - 6*SIZE(axL)), axY - 5, 1, RGB(90,90,90));
    xv := xv + stp;
  END;
END;

// ---- the probability for a region (the formulas of mode 1, unchanged) ----
// tailK: 1 = left of lo, 2 = right of lo, 3 = between lo and hi,
//        4 = outside lo..hi
sd_prob(mu, se, lo, hi, tailK)
BEGIN
  LOCAL pv;
  CASE
    IF tailK == 1 THEN pv := NORMALD_CDF(mu, se, lo); END;
    // v2.3a: a right tail is the LOWER tail at 2μ - a (the same number by
    // symmetry): 1 - NORMALD_CDF would be exactly 0 for a tail below the
    // CDF's last digit, where the left tail still reads "< 0.000001"
    IF tailK == 2 THEN pv := NORMALD_CDF(mu, se, 2*mu - lo); END;
    IF tailK == 3 THEN
      pv := NORMALD_CDF(mu, se, hi) - NORMALD_CDF(mu, se, lo);
    END;
    DEFAULT
      pv := NORMALD_CDF(mu, se, lo) + NORMALD_CDF(mu, se, 2*mu - hi);
  END;
  RETURN pv;
END;

// ---- critical values for a tail area (the formulas of mode 2) ----
// tailK: 1 lower, 2 upper, 3 two-tailed (central 1 - α).  Returns {lo, hi}.
sd_crit(mu, se, alp, tailK)
BEGIN
  LOCAL zlo, zhi, lo, hi;
  CASE
    IF tailK == 1 THEN
      zlo := NORMALD_ICDF(0,1,alp);
      lo := mu + zlo*se; hi := lo;
    END;
    IF tailK == 2 THEN
      zlo := NORMALD_ICDF(0,1,1-alp);
      lo := mu + zlo*se; hi := lo;
    END;
    DEFAULT
      zlo := NORMALD_ICDF(0,1,alp/2);
      zhi := NORMALD_ICDF(0,1,1-alp/2);
      lo := mu + zlo*se; hi := mu + zhi*se;
  END;
  RETURN {lo, hi};
END;

// ---- tidy the bounds after one has moved ----
// act says which bound moved (1 = lo, 2 = hi).  Critical-value mode,
// two-tailed: the other bound mirrors it about μ.  Two-bound regions keep
// lo <= hi (swapping which bound is active if the finger crossed over).
// Critical-value mode: α is read back from the moved bound, and kept
// inside 0.0000001 .. 0.9999999 so NORMALD_ICDF never sees 0 or 1.
// v2.3b (S1): the floor was 0.000001, which the panel's six decimals
// print as an exact "= 0.000001" -- a red line dragged to z = 16 read
// "P(x-bar > 144.000) = 0.000001" though α there is about 1E-58.  At
// 1E-7 the floor is below what six decimals show, so the panel says
// "< 0.000001", which is true.  NORMALD_ICDF is fine at 1E-7.
// Returns {lo, hi, act, alp}.
sd_fix(mu, se, lo, hi, tailK, crit, act, alp)
BEGIN
  LOCAL fxL, fxH, fxA, fxP, fxT;
  fxL := lo; fxH := hi; fxA := act; fxP := alp;
  IF crit == 1 AND tailK == 3 THEN
    IF fxA == 1 THEN fxH := 2*mu - fxL; ELSE fxL := 2*mu - fxH; END;
  END;
  IF tailK >= 3 AND fxL > fxH THEN
    fxT := fxL; fxL := fxH; fxH := fxT; fxA := 3 - fxA;
  END;
  IF crit == 1 THEN
    IF tailK == 1 THEN fxP := NORMALD_CDF(mu, se, fxL); END;
    // v2.3b (S2): the right tail as the LOWER tail at 2μ - a, as sd_prob
    // and Mode 1 have done since v2.3a (1 - NORMALD_CDF keeps only the
    // digits left after the 9s of a CDF near 1)
    IF tailK == 2 THEN fxP := NORMALD_CDF(mu, se, 2*mu - fxL); END;
    IF tailK == 3 THEN fxP := 2*NORMALD_CDF(mu, se, fxL); END;
    fxP := MAX(0.0000001, MIN(0.9999999, fxP));
  END;
  RETURN {fxL, fxH, fxA, fxP};
END;

// ---- compose one frame of the result graph in G1, blit, soft keys ----
// st  = {mu, sg, se, nn, lo, hi, tailK, pv, crit, act}
// geo = {xlo, xhi, ymx};  yX, yP = cached rows of x-bar's curve and of the
// population;  note = a one-line message for the panel ("" = the key help)
sd_compose(st, geo, yX, yP, note)
BEGIN
  LOCAL LFT, RGT, TOPy, BOT, W, H;
  LOCAL mu, sg, se, nn, lo, hi;
  LOCAL tailK, pv, act, xlo, xhi, ymx;
  LOCAL pxx, ddX, shade, tickx, cMk, cpN;
  LOCAL zlo, zhi, s1, s2, s3;

  LFT := 40; RGT := 314; TOPy := 20; BOT := 150;
  W := RGT - LFT; H := BOT - TOPy;
  mu := st(1); sg := st(2); se := st(3); nn := st(4);
  lo := st(5); hi := st(6); tailK := st(7); pv := st(8);
  act := st(10);
  xlo := geo(1); xhi := geo(2); ymx := geo(3);
  cMk := RGB(255,0,0);

  DIMGROB_P(G1,320,240);
  sd_grid(G1, xlo, xhi);

  // shading: iterate PIXEL COLUMNS so there are no gap columns; the
  // height of each column comes from the cache (no NORMALD here)
  FOR pxx FROM LFT+1 TO RGT DO
    ddX := xlo + (pxx - LFT)/W*(xhi - xlo);
    shade := 0;
    IF tailK == 1 AND ddX <= lo THEN shade := 1; END;
    IF tailK == 2 AND ddX >= lo THEN shade := 1; END;
    IF tailK == 3 AND ddX >= lo AND ddX <= hi THEN shade := 1; END;
    IF tailK == 4 AND (ddX <= lo OR ddX >= hi) THEN shade := 1; END;
    IF shade THEN
      IF yX(pxx - LFT + 1) < BOT THEN
        LINE_P(G1, pxx, BOT, pxx, yX(pxx - LFT + 1), RGB(170,205,255));
      END;
    END;
  END;

  // the population ghost, x-bar's curve and the mean marker
  sd_curves(G1, xlo, xhi, mu, sg, yX, yP);
  sd_axes(G1, xlo, xhi, ymx);

  // boundary markers; the ACTIVE bound (the one the arrows move) is
  // two pixels wide with a small red handle at the top
  tickx := LFT + IP((lo - xlo)/(xhi - xlo)*W);
  IF tickx >= LFT AND tickx <= RGT THEN
    LINE_P(G1, tickx, TOPy, tickx, BOT, cMk);
    IF act == 1 THEN
      LINE_P(G1, tickx + 1, TOPy, tickx + 1, BOT, cMk);
      RECT_P(G1, tickx - 2, TOPy, tickx + 3, TOPy + 4, cMk);
    END;
  END;
  IF tailK == 3 OR tailK == 4 THEN
    tickx := LFT + IP((hi - xlo)/(xhi - xlo)*W);
    IF tickx >= LFT AND tickx <= RGT THEN
      LINE_P(G1, tickx, TOPy, tickx, BOT, cMk);
      IF act == 2 THEN
        LINE_P(G1, tickx + 1, TOPy, tickx + 1, BOT, cMk);
        RECT_P(G1, tickx - 2, TOPy, tickx + 3, TOPy + 4, cMk);
      END;
    END;
  END;

  // header
  TEXTOUT_P(CHAR(57344) + " ~ Normal(μ = " + STRING(mu,2,2) + ", σ" +
            CHAR(57344) + " = " + STRING(se,2,4) + "),  n = " +
            STRING(nn,2,0), G1, 6, 4, 1, RGB(0,60,200), 308);

  // result panel (orange frame, light orange fill — house style)
  zlo := (lo - mu)/se;
  zhi := (hi - mu)/se;
  CASE
  // v2.3a: sd_Eq -- "P(x-bar > 1.000) < 0.000001", never "= < 0.000001"
    IF tailK == 1 THEN
      s1 := "P(" + CHAR(57344) + " < " + STRING(lo,2,3) + ")" + sd_Eq(sd_fmtp(pv,6));
      s2 := "z = " + STRING(zlo,2,4);
    END;
    IF tailK == 2 THEN
      s1 := "P(" + CHAR(57344) + " > " + STRING(lo,2,3) + ")" + sd_Eq(sd_fmtp(pv,6));
      s2 := "z = " + STRING(zlo,2,4);
    END;
    IF tailK == 3 THEN
      s1 := "P(" + STRING(lo,2,3) + " < " + CHAR(57344) + " < " +
            STRING(hi,2,3) + ")" + sd_Eq(sd_fmtp(pv,6));
      s2 := "z = " + STRING(zlo,2,4) + " to " + STRING(zhi,2,4);
    END;
    DEFAULT
      s1 := "P(outside " + STRING(lo,2,3) + " .. " + STRING(hi,2,3) +
            ")" + sd_Eq(sd_fmtp(pv,6));
      s2 := "z = " + STRING(zlo,2,4) + " / " + STRING(zhi,2,4);
  END;
  s2 := s2 + "  " + sd_Eq(sd_fmtp(pv*100,3)) + "%";
  // σ and σx-bar side by side, with the arithmetic between them: that σx-bar
  // is σ/√n and NOT σ is the one thing this program exists to teach (v2.2)
  s3 := "σ = " + STRING(sg,2,2) + ",  σ" + CHAR(57344) + " = " +
        STRING(sg,2,2) + "/√" + STRING(nn,2,0) + " = " + STRING(se,2,4);
  // the fourth line: a message if there is one, else the key help (the
  // footer lives here now -- the soft-key strip takes y 218-239)
  cpN := note;
  IF DIM(cpN) == 0 THEN
    cpN := "drag or " + CHAR(8646) + " = bound, Ent = more, Esc = menu";
  END;

  RECT_P(G1, 4, 166, 315, 216, RGB(255,69,0), RGB(255,240,220));
  TEXTOUT_P(s1, G1, 8, 169, 1, RGB(0,0,0), 305);
  TEXTOUT_P(s2, G1, 8, 180, 1, RGB(200,0,200), 305);
  TEXTOUT_P(s3, G1, 8, 191, 1, RGB(0,140,60), 305);
  TEXTOUT_P(cpN, G1, 8, 202, 1, RGB(30,144,255), 305);

  BLIT_P(G0, 0, 0, G1);
  // the soft keys are painted on the SCREEN, after the blit, every frame
  DRAWMENU({"n" + CHAR(8722), "n+", "n" + CHAR(215) + "4",
            CHAR(956) + "," + CHAR(963), "Tail", "Menu"});
END;

// ---- the live result graph ----
// Probability mode (crit = 0): the bounds are the student's; a new n
// changes P.  Critical-value mode (crit = 1): α is the student's; a new
// n moves the critical values.  All state lives in this routine's
// LOCALs (a PPL routine cannot change its caller's variables), and the
// parameters are copied first because INPUT must never write into a
// parameter (Lessons Learned section 2).
sd_live(mu0, sg0, nn0, lo0, hi0, tk0, crit, alp0)
BEGIN
  LOCAL LFT, RGT, TOPy, BOT, W, H;
  LOCAL mu, sg, nn, lo, hi, tailK;
  LOCAL alp, se, pv, xlo, xhi, ymx;
  LOCAL ww, vX, vP, yX, yP, note;
  LOCAL act, stp, fresh, refit, go, rd;
  LOCAL kk, mm, mm1, mx, dn, op;
  LOCAL ix, tmp, nl, iok, ch, cb;
  LOCAL mu2, sg2, nn2, al2, lo2, hi2;
  LOCAL zlo, zhi, mtx;

  LFT := 40; RGT := 314; TOPy := 20; BOT := 150;
  W := RGT - LFT; H := BOT - TOPy;
  mu := mu0; sg := sg0; nn := nn0; lo := lo0; hi := hi0;
  tailK := tk0; alp := alp0;
  // the n ladder for n− / n+
  nl := {1, 2, 4, 5, 10, 16, 20, 25, 30, 40, 50, 64, 100, 160, 200, 400};
  act := 1; dn := 0; go := 1; rd := 1; note := "";
  fresh := 1; refit := 1; ww := 0;

  // drain the OK that closed the last form, and any touch (section 7)
  REPEAT kk := GETKEY; UNTIL kk == -1;
  sd_Drain();

  WHILE go DO
    // ---- first frame, or a new μ / σ: fit the window and FREEZE it.
    // μ ± 16σx-bar shows the population's shoulders at moderate n and
    // keeps x-bar's curve a quarter of the plot wide; never wider than
    // μ ± 3.5σ (the whole population), never narrower than μ ± 4σx-bar.
    IF fresh THEN
      se := sg / SQRT(nn);
      ww := MAX(4*se, MIN(3.5*sg, 16*se));
      xlo := mu - ww; xhi := mu + ww;
      IF crit THEN
        cb := sd_crit(mu, se, alp, tailK);
        lo := cb(1); hi := cb(2);
      END;
      // stretched if a boundary falls outside it (as v2.2 did)
      IF lo < xlo THEN xlo := lo - 0.5*se; END;
      IF hi > xhi THEN xhi := hi + 0.5*se; END;
      // the population's density depends on μ and σ only: cached here,
      // reused for every n
      vP := sd_pdfList(mu, sg, xlo, xhi, W);
      fresh := 0; refit := 1;
    END;
    // ---- a new n: x-bar's curve, the y-scale (refitted, because n sets
    // the peak height), the arrow step, and in critical-value mode the
    // critical values for the same α
    IF refit THEN
      se := sg / SQRT(nn);
      IF crit THEN
        cb := sd_crit(mu, se, alp, tailK);
        lo := cb(1); hi := cb(2);
      END;
      vX := sd_pdfList(mu, se, xlo, xhi, W);
      ymx := NORMALD(mu, se, mu)*1.15;
      yX := sd_pyList(vX, ymx, BOT, TOPy, H);
      yP := sd_pyList(vP, ymx, BOT, TOPy, H);
      stp := sd_NiceNum(se/4, 1);
      refit := 0; rd := 1;
    END;
    // ---- the numbers and the frame (only when something changed: a
    // finger resting still must not redraw on every poll)
    IF rd THEN
      IF crit THEN
        // α is the student's number: P is α, or 1 - α for the central region
        IF tailK == 3 THEN pv := 1 - alp; ELSE pv := alp; END;
      ELSE
        pv := sd_prob(mu, se, lo, hi, tailK);
      END;
      zlo := (lo - mu)/se;
      zhi := (hi - mu)/se;
      EP_SD_SE := se; EP_SD_P := pv;
      EP_SD_LO := lo; EP_SD_ZLO := zlo;
      IF tailK <= 2 THEN
        EP_SD_HI := lo; EP_SD_ZHI := zlo;
      ELSE
        EP_SD_HI := hi; EP_SD_ZHI := zhi;
      END;
      sd_compose({mu, sg, se, nn, lo, hi, tailK, pv, crit, act},
                 {xlo, xhi, ymx}, yX, yP, note);
      rd := 0;
    END;

    // ---- wait for a key or a touch.  While a finger is down (dn = 1)
    // the loop comes straight back, so a lift is seen at once.
    REPEAT
      kk := GETKEY;
      mm := MOUSE;
      mm1 := mm(1);
    UNTIL kk > -1 OR SIZE(mm1) > 0 OR dn == 1;

    op := 0;
    IF SIZE(mm1) > 0 THEN
      kk := -1;                              // a touch is not a keypress
      mx := mm1(1);
      IF dn == 0 AND mm1(2) >= 218 THEN
        // the soft-key strip: six 53-px slots, left to right
        sd_Drain();
        op := MIN(6, IP(mx / 53.4) + 1);
      ELSE
        IF dn == 1 OR (mx >= LFT AND mx <= RGT AND mm1(2) >= TOPy AND mm1(2) <= BOT) THEN
          // the finger is in the plot: the x it points at, through the
          // SAME box constants the drawing uses
          mx := MAX(LFT, MIN(RGT, mx));
          tmp := xlo + (mx - LFT)/W*(xhi - xlo);
          IF dn == 0 THEN
            // the press: pick the NEARER bound and keep it for the drag
            act := 1;
            IF tailK >= 3 AND ABS(tmp - hi) < ABS(tmp - lo) THEN act := 2; END;
            dn := 1;
          END;
          IF act == 1 THEN lo := tmp; ELSE hi := tmp; END;
          cb := sd_fix(mu, se, lo, hi, tailK, crit, act, alp);
          lo := cb(1); hi := cb(2); act := cb(3); alp := cb(4);
          note := ""; rd := 1;
        END;
      END;
    ELSE
      IF dn == 1 THEN
        // the finger has lifted: snap the bound that moved onto the step
        dn := 0;
        IF act == 1 THEN lo := ROUND(lo/stp, 0)*stp; ELSE hi := ROUND(hi/stp, 0)*stp; END;
        cb := sd_fix(mu, se, lo, hi, tailK, crit, act, alp);
        lo := cb(1); hi := cb(2); act := cb(3); alp := cb(4);
        rd := 1;
      END;
    END;

    // ---- keys: shortcuts for everything a finger does
    IF kk == 4 THEN go := 0; END;                      // Esc = menu
    IF kk == 50 THEN op := 2; END;                     // +   = n+
    IF kk == 45 THEN op := 1; END;                     // -   = n-
    IF kk == 7 OR kk == 8 THEN                         // Left / Right
      // move the active bound to the next multiple of the step, so an odd
      // value is pulled onto the grid and stays round from then on
      IF act == 1 THEN tmp := lo; ELSE tmp := hi; END;
      IF kk == 7 THEN tmp := (CEILING(tmp/stp - 0.000001) - 1)*stp; END;
      IF kk == 8 THEN tmp := (FLOOR(tmp/stp + 0.000001) + 1)*stp; END;
      IF act == 1 THEN lo := tmp; ELSE hi := tmp; END;
      cb := sd_fix(mu, se, lo, hi, tailK, crit, act, alp);
      lo := cb(1); hi := cb(2); act := cb(3); alp := cb(4);
      note := ""; rd := 1;
    END;
    IF (kk == 2 OR kk == 12) AND tailK >= 3 THEN       // Up / Down
      act := 3 - act; rd := 1;                         // the other bound
    END;
    IF kk == 30 THEN                                   // Enter = more
      ch := 0;
      IF crit THEN mtx := "Edit μ, σ, n, α"; ELSE mtx := "Edit μ, σ, n"; END;
      CHOOSE(ch, "More", {"n " + CHAR(215) + " 4", mtx, "Next tail type", "Back to menu"});
      REPEAT kk := GETKEY; UNTIL kk == -1;
      sd_Drain();
      IF ch > 0 THEN op := ch + 2; END;                // 3, 4, 5 or 6
      rd := 1;                                         // CHOOSE drew over us
    END;

    // ---- the actions behind the soft keys
    // 1 n−, 2 n+, 3 n×4: a new n.  The window stays frozen (the population
    // stays put) unless x-bar's curve would no longer fit in it.
    IF op >= 1 AND op <= 3 THEN
      note := "";
      nn2 := nn;
      IF op == 1 THEN                                  // the next rung down
        FOR ix FROM 1 TO 16 DO
          IF nl(ix) < nn THEN nn2 := nl(ix); END;
        END;
        IF nn2 == nn THEN note := "n is already at its smallest, 1"; END;
      END;
      IF op == 2 THEN                                  // the next rung up
        tmp := 0;
        FOR ix FROM 1 TO 16 DO
          IF tmp == 0 AND nl(ix) > nn THEN nn2 := nl(ix); tmp := 1; END;
        END;
        IF nn2 == nn THEN note := "n+ stops at 400; n" + CHAR(215) + "4 goes on to 1600"; END;
      END;
      IF op == 3 THEN
        IF nn*4 <= 1600 THEN
          nn2 := nn*4;
        ELSE
          note := "n" + CHAR(215) + "4 stops at n = 1600";
        END;
      END;
      IF nn2 <> nn THEN
        nn := nn2;
        IF 2.5*sg/SQRT(nn) > ww THEN
          fresh := 1;
          note := "window refitted: the curve outgrew it";
        ELSE
          refit := 1;
        END;
      END;
      rd := 1;
    END;
    // 4 μ,σ: the form again (with α in critical-value mode); then, in
    // probability mode, the bounds form, so the bounds can follow a new μ
    IF op == 4 THEN
      note := "";
      mu2 := mu; sg2 := sg; nn2 := nn; al2 := alp;
      IF crit THEN
        iok := INPUT({mu2, sg2, nn2, al2}, "Critical Values",
                     {"μ:", "σ:", "n:", "α:"},
                     {"Population mean", "Population SD (> 0)",
                      "Sample size, whole number >= 1", "Tail area, e.g. 0.05"});
      ELSE
        iok := INPUT({mu2, sg2, nn2}, "Population Parameters",
                     {"μ:", "σ:", "n:"},
                     {"Population mean", "Population SD (> 0)",
                      "Sample size, whole number >= 1"});
      END;
      REPEAT kk := GETKEY; UNTIL kk == -1;
      sd_Drain();
      IF iok THEN
        CASE
          IF sg2 <= 0 THEN MSGBOX("σ must be > 0; nothing changed"); END;
          IF nn2 < 1 THEN MSGBOX("n must be >= 1; nothing changed"); END;
          IF nn2 - IP(nn2) > 0 THEN MSGBOX("n must be a whole number; nothing changed"); END;
          IF al2 <= 0 OR al2 >= 1 THEN MSGBOX("α must be between 0 and 1; nothing changed"); END;
          DEFAULT
            mu := mu2; sg := sg2; nn := nn2; alp := al2;
            IF crit == 0 THEN
              lo2 := lo; hi2 := hi;
              IF tailK <= 2 THEN
                iok := INPUT({lo2}, "Enter the boundary", {"a:"},
                             {"Boundary value on the " + CHAR(57344) + " scale"});
              ELSE
                iok := INPUT({lo2, hi2}, "Enter the bounds", {"a:", "b:"},
                             {"Lower bound", "Upper bound"});
              END;
              REPEAT kk := GETKEY; UNTIL kk == -1;
              sd_Drain();
              IF iok THEN
                IF tailK >= 3 AND hi2 < lo2 THEN tmp := lo2; lo2 := hi2; hi2 := tmp; END;
                lo := lo2;
                IF tailK >= 3 THEN hi := hi2; END;
              END;
            END;
            fresh := 1;
        END;
      END;
      rd := 1;
    END;
    // 5 Tail: the next tail type.  Probability mode cycles left, right,
    // between, outside; critical-value mode lower, upper, two-tailed (α
    // kept).  Entering a two-bound region from a one-bound one makes the
    // pair symmetric about μ at the old distance (at least σx-bar).
    IF op == 5 THEN
      note := "";
      IF crit THEN
        IF tailK >= 3 THEN tailK := 1; ELSE tailK := tailK + 1; END;
        refit := 1;
      ELSE
        IF tailK >= 4 THEN tailK := 1; ELSE tailK := tailK + 1; END;
        IF tailK == 3 AND hi <= lo THEN
          tmp := ABS(lo - mu);
          IF tmp < se/2 THEN tmp := se; END;
          lo := mu - tmp; hi := mu + tmp;
        END;
      END;
      act := 1; rd := 1;
    END;
    IF op == 6 THEN go := 0; END;                      // Menu
  END;

  REPEAT kk := GETKEY; UNTIL kk == -1;
  sd_Drain();
END;

// ---- comparison view: three sample sizes on shared axes ----
// sg is the population σ, for the title.  v2.1 reconstructed it as
// semax*√n1, which is only σ when n1 happens to be the smallest n.
// v2.3: filled curves (widest first), an 8-frame opening in which n
// grows from 1, a finger cursor with P(x-bar < x) for each n, and a tap
// that picks a curve.  The finished picture is composed once into G2;
// each cursor frame copies G2 into G1 and adds the cursor (GaltonBoard's
// scene + frame idiom), so a drag costs two blits, not 1,700 lines.
sd_drawCmp(mu, sg, se1, se2, se3, n1, n2, n3)
BEGIN
  LOCAL LFT, RGT, TOPy, BOT, W, H;
  LOCAL xlo, xhi, ymx, semax, semin, kk;
  LOCAL stp, nfrac, xv, tickx, dcLbl, cAx;
  LOCAL cGrid, ses, nns, cls, lts, odr;
  LOCAL ys, yy, fr, ix, jx, tmp;
  LOCAL mm, mm1, dn, cx, hl, bg;
  LOCAL x0, y0, mv, skip, go, rd;
  LOCAL qg, pc, dd, best;

  LFT := 30; RGT := 314; TOPy := 22; BOT := 172;
  W := RGT - LFT; H := BOT - TOPy;
  cAx := RGB(0,0,0); cGrid := RGB(225,228,235);

  // widest curve (largest SE) sets the x-range; tallest (smallest SE) sets
  // the y-scale, so no curve is clipped and the equal-area lesson survives
  semax := MAX(se1, MAX(se2, se3));
  semin := MIN(se1, MIN(se2, se3));
  xlo := mu - 4*semax;
  xhi := mu + 4*semax;
  ymx := NORMALD(mu, semin, mu)*1.15;
  stp := sd_NiceNum((xhi - xlo)/6, 1);
  nfrac := MAX(0, -FLOOR(LOG(stp)));

  // the three curves as lists, in legend order: n1 red, n2 blue, n3 green
  ses := {se1, se2, se3};
  nns := {n1, n2, n3};
  cls := {RGB(255,0,0), RGB(0,60,200), RGB(0,140,60)};
  lts := {RGB(255,235,235), RGB(235,242,255), RGB(232,246,235)};
  // fill order: widest (largest SE) first, so each narrower, taller curve
  // is painted over the wider one and all three areas stay visible
  odr := {1, 2, 3};
  FOR ix FROM 1 TO 2 DO
    FOR jx FROM 1 TO 3 - ix DO
      IF ses(odr(jx)) < ses(odr(jx + 1)) THEN
        tmp := odr(jx); odr(jx) := odr(jx + 1); odr(jx + 1) := tmp;
      END;
    END;
  END;
  qg := {xlo, xhi, ymx, LFT, RGT, TOPy, BOT};

  // ---- the opening: n grows from 1 to its value over 8 frames, so all
  // three curves start as the population and squeeze.  Any key or tap
  // skips to the finished picture.
  REPEAT kk := GETKEY; UNTIL kk == -1;
  sd_Drain();
  skip := 0;
  FOR fr FROM 1 TO 8 DO
    IF skip == 0 THEN
      DIMGROB_P(G1,320,240);
      RECT_P(G1,0,0,319,239,RGB(255,255,255));
      xv := CEILING(xlo/stp)*stp;
      WHILE xv <= xhi + stp/1000 DO
        tickx := LFT + IP((xv - xlo)/(xhi - xlo)*W);
        IF tickx > LFT AND tickx <= RGT THEN
          LINE_P(G1, tickx, TOPy, tickx, BOT-1, cGrid);
        END;
        xv := xv + stp;
      END;
      FOR ix FROM 1 TO 3 DO
        // n^(fr/8) runs from about 1 to n; σx-bar = σ/√(that)
        sd_quick(G1, mu, sg/SQRT(nns(ix)^(fr/8)), qg, cls(ix));
      END;
      LINE_P(G1, LFT, BOT, RGT, BOT, cAx);
      LINE_P(G1, LFT, BOT, LFT, TOPy, cAx);
      TEXTOUT_P("Effect of sample size on σ" + CHAR(57344) + "   (μ = " +
                STRING(mu,2,2) + ", σ = " + STRING(sg,2,2) + ")",
                G1, 6, 4, 1, RGB(0,60,200), 308);
      TEXTOUT_P("n growing:  " + STRING(nns(1)^(fr/8),2,0) + " / " +
                STRING(nns(2)^(fr/8),2,0) + " / " + STRING(nns(3)^(fr/8),2,0),
                G1, 8, 186, 1, RGB(90,90,90), 300);
      TEXTOUT_P("any key or tap = skip", G1, 6, 226, 1, RGB(30,144,255));
      BLIT_P(G0, 0, 0, G1);
      kk := GETKEY;
      mm := MOUSE;
      mm1 := mm(1);
      IF kk > -1 OR SIZE(mm1) > 0 THEN skip := 1; END;
    END;
  END;
  REPEAT kk := GETKEY; UNTIL kk == -1;
  sd_Drain();

  // ---- the finished curves, cached as rows (one list per curve)
  ys := {sd_pyList(sd_pdfList(mu, se1, xlo, xhi, W), ymx, BOT, TOPy, H),
         sd_pyList(sd_pdfList(mu, se2, xlo, xhi, W), ymx, BOT, TOPy, H),
         sd_pyList(sd_pdfList(mu, se3, xlo, xhi, W), ymx, BOT, TOPy, H)};

  hl := 0; cx := -1; dn := 0; bg := 1; go := 1; rd := 1;
  x0 := 0; y0 := 0; mv := 0;
  WHILE go DO
    // ---- the scene into G2 (again only when the picked curve changes)
    IF bg THEN
      DIMGROB_P(G2,320,240);
      RECT_P(G2,0,0,319,239,RGB(255,255,255));
      xv := CEILING(xlo/stp)*stp;
      WHILE xv <= xhi + stp/1000 DO
        tickx := LFT + IP((xv - xlo)/(xhi - xlo)*W);
        IF tickx > LFT AND tickx <= RGT THEN
          LINE_P(G2, tickx, TOPy, tickx, BOT-1, cGrid);
        END;
        xv := xv + stp;
      END;
      // fills, widest first; then every outline over the fills
      FOR ix FROM 1 TO 3 DO
        jx := odr(ix);
        yy := ys(jx);
        sd_fill(G2, yy, LFT, BOT, lts(jx));
      END;
      FOR ix FROM 1 TO 3 DO
        yy := ys(ix);
        sd_poly(G2, yy, LFT, BOT, cls(ix), 0);
      END;
      // the picked curve, three pixels thick, with its n and σx-bar
      IF hl > 0 THEN
        yy := ys(hl);
        sd_poly(G2, yy, LFT, BOT, cls(hl), -1);
        sd_poly(G2, yy, LFT, BOT, cls(hl), 1);
        TEXTOUT_P("n = " + STRING(nns(hl),2,0) + ":  σ" + CHAR(57344) + " = " +
                  STRING(ses(hl),2,4), G2, LFT + 4, TOPy + 2, 1, cls(hl));
      END;
      LINE_P(G2, LFT, BOT, RGT, BOT, cAx);
      LINE_P(G2, LFT, BOT, LFT, TOPy, cAx);
      xv := CEILING(xlo/stp)*stp;
      WHILE xv <= xhi + stp/1000 DO
        tickx := LFT + IP((xv - xlo)/(xhi - xlo)*W);
        IF tickx >= LFT AND tickx <= RGT THEN
          LINE_P(G2, tickx, BOT, tickx, BOT+3, cAx);
          dcLbl := STRING(xv,2,nfrac);
          TEXTOUT_P(dcLbl, G2, tickx - 3*SIZE(dcLbl), BOT+5, 1, RGB(90,90,90));
        END;
        xv := xv + stp;
      END;

      TEXTOUT_P("Effect of sample size on σ" + CHAR(57344) + "   (μ = " +
                STRING(mu,2,2) + ", σ = " + STRING(sg,2,2) + ")",
                G2, 6, 4, 1, RGB(0,60,200), 308);

      RECT_P(G2, 4, 182, 315, 220, RGB(255,69,0), RGB(255,240,220));
      TEXTOUT_P("n = " + STRING(n1,1,0) + ": σ" + CHAR(57344) + " = " +
                STRING(se1,2,4), G2, 8, 186, 1, RGB(255,0,0));
      TEXTOUT_P("n = " + STRING(n2,1,0) + ": σ" + CHAR(57344) + " = " +
                STRING(se2,2,4), G2, 8, 197, 1, RGB(0,60,200));
      TEXTOUT_P("n = " + STRING(n3,1,0) + ": σ" + CHAR(57344) + " = " +
                STRING(se3,2,4), G2, 170, 186, 1, RGB(0,140,60));
      TEXTOUT_P("4× n halves σ" + CHAR(57344), G2, 170, 197, 1, RGB(90,90,90));
      TEXTOUT_P("touch the plot: P(" + CHAR(57344) + " < x) for each n",
                G2, 8, 208, 1, RGB(90,90,90), 300);

      // house footer: dodger blue, "key = action"
      TEXTOUT_P("drag = cursor, tap a curve = pick, any key = menu",
                G2, 6, 226, 1, RGB(30,144,255), 308);
      bg := 0; rd := 1;
    END;

    // ---- one frame: the scene, the cursor and its three readings
    IF rd THEN
      BLIT_P(G1, 0, 0, G2);
      IF cx >= LFT AND cx <= RGT THEN
        LINE_P(G1, cx, TOPy, cx, BOT, RGB(0,0,0));
        tmp := xlo + (cx - LFT)/W*(xhi - xlo);
        RECT_P(G1, 6, 207, 313, 218, RGB(255,240,220));
        TEXTOUT_P("P(" + CHAR(57344) + " < " + STRING(tmp,2,nfrac+1) + "):",
                  G1, 8, 208, 1, RGB(0,0,0));
        FOR ix FROM 1 TO 3 DO
          pc := NORMALD_CDF(mu, ses(ix), tmp);
          TEXTOUT_P(sd_fmtp(pc,4), G1, 128 + 60*(ix - 1), 208, 1, cls(ix));
          yy := ys(ix);
          RECT_P(G1, cx - 1, yy(cx - LFT + 1) - 1, cx + 1, yy(cx - LFT + 1) + 1, cls(ix));
        END;
      END;
      BLIT_P(G0, 0, 0, G1);
      rd := 0;
    END;

    // ---- wait for a key or a touch (straight back while a finger is down)
    REPEAT
      kk := GETKEY;
      mm := MOUSE;
      mm1 := mm(1);
    UNTIL kk > -1 OR SIZE(mm1) > 0 OR dn == 1;
    IF SIZE(mm1) > 0 THEN
      IF dn == 0 THEN
        // the press: remember where, to tell a tap from a drag on lift
        x0 := mm1(1); y0 := mm1(2); mv := 0; dn := 1;
      END;
      tmp := MAX(ABS(mm1(1) - x0), ABS(mm1(2) - y0));
      IF tmp > mv THEN mv := tmp; END;
      IF mm1(1) >= LFT AND mm1(1) <= RGT AND mm1(1) <> cx THEN
        cx := mm1(1); rd := 1;
      END;
    ELSE
      IF dn == 1 THEN
        dn := 0;
        // moved no more than 8 px: a tap.  Pick the curve whose row at the
        // tapped column is within 6 px of the finger (and above the axis);
        // a tap away from every curve clears the pick.
        IF mv <= 8 AND x0 >= LFT AND x0 <= RGT THEN
          best := 0; dd := 7;
          FOR ix FROM 1 TO 3 DO
            yy := ys(ix);
            IF yy(x0 - LFT + 1) < BOT - 2 AND ABS(y0 - yy(x0 - LFT + 1)) < dd THEN
              dd := ABS(y0 - yy(x0 - LFT + 1)); best := ix;
            END;
          END;
          IF best <> hl THEN hl := best; bg := 1; END;
        END;
      ELSE
        IF kk > -1 THEN go := 0; END;                  // any key = menu
      END;
    END;
  END;
  REPEAT kk := GETKEY; UNTIL kk == -1;
  sd_Drain();
END;

// ---- Mode 1: find a probability ----
sd_modeProb()
BEGIN
  LOCAL mu, sg, nn, se, iok;
  LOCAL lo, hi, tailK, pv;
  LOCAL zlo, zhi, tmp;

  mu := 0; sg := 1; nn := 30;
  iok := INPUT({mu, sg, nn}, "Population Parameters",
               {"μ:", "σ:", "n:"},
               {"Population mean", "Population SD (> 0)",
                "Sample size, whole number >= 1"});
  IF iok == 0 THEN RETURN; END;
  IF sg <= 0 THEN MSGBOX("σ must be > 0"); RETURN; END;
  IF nn < 1 THEN MSGBOX("n must be >= 1"); RETURN; END;
  IF nn - IP(nn) > 0 THEN MSGBOX("n must be a whole number"); RETURN; END;

  se := sg / SQRT(nn);

  tailK := 0;
  CHOOSE(tailK, "Which probability?",
    {"Left tail  P(" + CHAR(57344) + " < a)",
     "Right tail  P(" + CHAR(57344) + " > a)",
     "Between  P(a < " + CHAR(57344) + " < b)",
     "Outside  P(" + CHAR(57344) + " < a or > b)"});
  IF tailK == 0 THEN RETURN; END;

  lo := mu - se; hi := mu + se;
  IF tailK == 1 OR tailK == 2 THEN
    iok := INPUT({lo}, "Enter the boundary", {"a:"},
                 {"Boundary value on the " + CHAR(57344) + " scale"});
    IF iok == 0 THEN RETURN; END;
    hi := lo;
  ELSE
    iok := INPUT({lo, hi}, "Enter the bounds", {"a:", "b:"},
                 {"Lower bound", "Upper bound"});
    IF iok == 0 THEN RETURN; END;
    IF hi < lo THEN tmp := lo; lo := hi; hi := tmp; END;
  END;

  // NORMALD_CDF is evaluated on the x-bar scale directly — standardising first
  // and calling the standard normal would give the same number by one more
  // subtraction. z is reported, not used.
  zlo := (lo - mu)/se;
  zhi := (hi - mu)/se;
  CASE
    IF tailK == 1 THEN pv := NORMALD_CDF(mu, se, lo); END;
    // v2.3a: the right tail as the lower tail at 2μ - a, as in sd_prob
    IF tailK == 2 THEN pv := NORMALD_CDF(mu, se, 2*mu - lo); END;
    IF tailK == 3 THEN
      pv := NORMALD_CDF(mu, se, hi) - NORMALD_CDF(mu, se, lo);
    END;
    DEFAULT
      pv := NORMALD_CDF(mu, se, lo) + NORMALD_CDF(mu, se, 2*mu - hi);
  END;

  EP_SD_SE := se; EP_SD_P := pv;
  EP_SD_ZLO := zlo; EP_SD_ZHI := zhi;
  EP_SD_LO := lo; EP_SD_HI := hi;

  // the live graph recomputes the same numbers with sd_prob (the same
  // formulas) every time n, a bound or the tail changes
  sd_live(mu, sg, nn, lo, hi, tailK, 0, 0.05);
END;

// ---- Mode 2: find a critical value ----
sd_modeCrit()
BEGIN
  LOCAL mu, sg, nn, se, iok;
  LOCAL alp, tailK, pv;
  LOCAL zlo, zhi, lo, hi;

  mu := 0; sg := 1; nn := 30; alp := 0.05;
  // ONE form rather than two sequential prompts: fewer cancel paths.
  iok := INPUT({mu, sg, nn, alp}, "Critical Values",
               {"μ:", "σ:", "n:", "α:"},
               {"Population mean", "Population SD (> 0)",
                "Sample size, whole number >= 1", "Tail area, e.g. 0.05"});
  IF iok == 0 THEN RETURN; END;
  IF sg <= 0 THEN MSGBOX("σ must be > 0"); RETURN; END;
  IF nn < 1 THEN MSGBOX("n must be >= 1"); RETURN; END;
  IF nn - IP(nn) > 0 THEN MSGBOX("n must be a whole number"); RETURN; END;
  IF alp <= 0 OR alp >= 1 THEN MSGBOX("α must be between 0 and 1"); RETURN; END;

  se := sg / SQRT(nn);

  tailK := 0;
  CHOOSE(tailK, "Critical value type",
    {"Lower tail", "Upper tail", "Two-tailed (central)"});
  IF tailK == 0 THEN RETURN; END;

  CASE
    IF tailK == 1 THEN
      zlo := NORMALD_ICDF(0,1,alp);
      lo := mu + zlo*se; hi := lo; zhi := zlo; pv := alp;
    END;
    IF tailK == 2 THEN
      zlo := NORMALD_ICDF(0,1,1-alp);
      lo := mu + zlo*se; hi := lo; zhi := zlo; pv := alp;
    END;
    DEFAULT
      zlo := NORMALD_ICDF(0,1,alp/2);
      zhi := NORMALD_ICDF(0,1,1-alp/2);
      lo := mu + zlo*se; hi := mu + zhi*se; pv := 1 - alp;
  END;

  EP_SD_SE := se; EP_SD_P := pv;
  EP_SD_ZLO := zlo; EP_SD_ZHI := zhi;
  EP_SD_LO := lo; EP_SD_HI := hi;

  // the live graph keeps α and recomputes the critical values (sd_crit,
  // the same formulas) whenever n changes
  sd_live(mu, sg, nn, lo, hi, tailK, 1, alp);
END;

// ---- Mode 3: compare three sample sizes ----
sd_modeCompare()
BEGIN
  LOCAL mu, sg, n1, n2, n3, iok;
  LOCAL se1, se2, se3;

  // Defaults 10, 40, 160: each n is 4x the one before, so the picture
  // shows exactly the "4x n halves σx-bar" rule printed under it.  The old
  // 5 / 30 / 100 showed a 6x and a 3.3x, which the slogan did not describe.
  mu := 100; sg := 15; n1 := 10; n2 := 40; n3 := 160;
  iok := INPUT({mu, sg, n1, n2, n3}, "Compare Sample Sizes",
               {"μ:", "σ:", "n₁:", "n₂:", "n₃:"},
               {"Population mean", "Population SD (> 0)",
                "Small n", "Medium n", "Large n"});
  IF iok == 0 THEN RETURN; END;
  IF sg <= 0 THEN MSGBOX("σ must be > 0"); RETURN; END;
  IF n1 < 1 OR n2 < 1 OR n3 < 1 THEN
    MSGBOX("every n must be >= 1"); RETURN;
  END;
  // same whole-number check as modes 1 and 2: n = 2.5 has no meaning
  IF n1 - IP(n1) > 0 OR n2 - IP(n2) > 0 OR n3 - IP(n3) > 0 THEN
    MSGBOX("every n must be a whole number"); RETURN;
  END;

  se1 := sg / SQRT(n1);
  se2 := sg / SQRT(n2);
  se3 := sg / SQRT(n3);

  EP_SD_SE := se2;

  sd_drawCmp(mu, sg, se1, se2, se3, n1, n2, n3);
END;

// ---- Mode 4: watch sample means (v2.3) ----
// Forty samples of n draws from N(μ, σ) are simulated in PPL (sd_Norm);
// each sample's mean drops as a dot onto the axis under x-bar's curve,
// stacking where means pile up.  Then the forty become a histogram (bars
// one σx-bar wide, scaled to area 1 so they sit on the same scale as the
// curve), with their mean and SD printed beside μ and σx-bar.  Above
// n = 100 each mean is drawn directly from N(μ, σx-bar) -- the same
// distribution, without 40 x n draws -- and the screen says so.
sd_watch()
BEGIN
  LOCAL LFT, RGT, TOPy, BOT, W, H;
  LOCAL mu, sg, nn, se, iok, ww;
  LOCAL xlo, xhi, ymx, yX, yP, ms;
  LOCAL stk, ix, jx, sm, mv, px;
  LOCAL bn, yl, fr, nf, skip, kk;
  LOCAL mm, mm1, t0, cnt, e1, e2;
  LOCAL hy, bx1, bx2, mbar, sd2, drct;
  LOCAL s1, s2, s3, s4, yf;

  LFT := 40; RGT := 314; TOPy := 20; BOT := 150;
  W := RGT - LFT; H := BOT - TOPy;

  mu := 100; sg := 15; nn := 30;
  iok := INPUT({mu, sg, nn}, "Watch Sample Means",
               {"μ:", "σ:", "n:"},
               {"Population mean", "Population SD (> 0)",
                "Draws per sample, whole number >= 1"});
  IF iok == 0 THEN RETURN; END;
  IF sg <= 0 THEN MSGBOX("σ must be > 0"); RETURN; END;
  IF nn < 1 THEN MSGBOX("n must be >= 1"); RETURN; END;
  IF nn - IP(nn) > 0 THEN MSGBOX("n must be a whole number"); RETURN; END;

  se := sg / SQRT(nn);
  drct := 0;
  IF nn > 100 THEN drct := 1; END;
  // the same window and scale as the result graph
  ww := MAX(4*se, MIN(3.5*sg, 16*se));
  xlo := mu - ww; xhi := mu + ww;
  ymx := NORMALD(mu, se, mu)*1.15;
  yX := sd_pyList(sd_pdfList(mu, se, xlo, xhi, W), ymx, BOT, TOPy, H);
  yP := sd_pyList(sd_pdfList(mu, sg, xlo, xhi, W), ymx, BOT, TOPy, H);

  // ---- the scene, once, into G2: each animation frame copies it to G1
  // and adds the falling dot, and each landed dot is added to G2
  DIMGROB_P(G2,320,240);
  sd_grid(G2, xlo, xhi);
  sd_curves(G2, xlo, xhi, mu, sg, yX, yP);
  sd_axes(G2, xlo, xhi, ymx);
  TEXTOUT_P("40 means of n = " + STRING(nn,2,0) + " draws from N(" +
            STRING(mu,2,2) + ", " + STRING(sg,2,2) + ")", G2, 6, 4, 1, RGB(0,60,200), 308);
  RECT_P(G2, 4, 166, 315, 216, RGB(255,69,0), RGB(255,240,220));
  TEXTOUT_P("each dot is ONE sample's mean " + CHAR(57344), G2, 8, 191, 1, RGB(0,0,0), 305);
  IF drct THEN
    TEXTOUT_P("n > 100: each mean drawn from N(μ, σ" + CHAR(57344) + ")",
              G2, 8, 202, 1, RGB(90,90,90), 305);
  END;
  TEXTOUT_P("any key or tap = skip", G2, 6, 226, 1, RGB(30,144,255));

  stk := MAKELIST(0, K, 1, 70, 1);       // dots per 4-px column, for stacking
  ms := MAKELIST(0, K, 1, 40, 1);        // the forty means
  skip := 0;
  REPEAT kk := GETKEY; UNTIL kk == -1;
  sd_Drain();

  FOR ix FROM 1 TO 40 DO
    // one sample of n draws from N(μ, σ), and its mean
    IF drct THEN
      mv := mu + se*sd_Norm();
    ELSE
      sm := 0;
      FOR jx FROM 1 TO nn DO
        sm := sm + mu + sg*sd_Norm();
      END;
      mv := sm/nn;
    END;
    ms(ix) := mv;
    // where it lands: its column, stacked on the dots already there
    px := LFT + IP((mv - xlo)/(xhi - xlo)*W);
    px := MAX(LFT + 1, MIN(RGT - 2, px));
    bn := IP((px - LFT)/4) + 1;
    yl := BOT - 3 - 4*stk(bn);
    IF yl < TOPy + 2 THEN yl := TOPy + 2; END;
    stk(bn) := stk(bn) + 1;
    // the fall: five frames for the first six dots, two after that
    IF skip == 0 THEN
      nf := 2;
      IF ix <= 6 THEN nf := 5; END;
      FOR fr FROM 1 TO nf DO
        IF skip == 0 THEN
          yf := TOPy + IP((yl - TOPy)*fr/nf);
          BLIT_P(G1, 0, 0, G2);
          RECT_P(G1, px - 1, yf - 1, px + 2, yf + 2, RGB(255,0,0));
          TEXTOUT_P("sample " + STRING(ix,2,0) + " of 40:  " + CHAR(57344) +
                    " = " + STRING(mv,2,3), G1, 8, 169, 1, RGB(200,0,200), 305);
          BLIT_P(G0, 0, 0, G1);
          t0 := TICKS;
          REPEAT UNTIL TICKS - t0 > 25;      // about 40 frames a second at most
          kk := GETKEY;
          mm := MOUSE;
          mm1 := mm(1);
          IF kk > -1 OR SIZE(mm1) > 0 THEN skip := 1; END;
        END;
      END;
    END;
    // landed: the dot joins the scene (magenta, CentralLimit's colour
    // for a pile of means)
    RECT_P(G2, px - 1, yl - 1, px + 2, yl + 2, RGB(200,0,200));
  END;
  REPEAT kk := GETKEY; UNTIL kk == -1;
  sd_Drain();

  // ---- the forty means as a histogram under the curve.  Bars one
  // σx-bar wide, edges at μ + kσx-bar; height = count / (40 σx-bar), a
  // density, so the bars and the curve share the y-scale (area 1).
  DIMGROB_P(G1,320,240);
  sd_grid(G1, xlo, xhi);
  FOR ix FROM -8 TO 7 DO
    e1 := mu + ix*se; e2 := e1 + se;
    cnt := 0;
    FOR jx FROM 1 TO 40 DO
      IF ms(jx) >= e1 AND ms(jx) < e2 THEN cnt := cnt + 1; END;
    END;
    IF cnt > 0 THEN
      hy := BOT - IP(cnt/(40*se)/ymx*H);
      IF hy < TOPy THEN hy := TOPy; END;
      bx1 := MAX(LFT + 1, LFT + IP((e1 - xlo)/(xhi - xlo)*W));
      bx2 := MIN(RGT, LFT + IP((e2 - xlo)/(xhi - xlo)*W) - 1);
      IF bx2 > bx1 THEN
        RECT_P(G1, bx1, hy, bx2, BOT - 1, RGB(200,0,200), RGB(250,235,250));
      END;
    END;
  END;
  sd_curves(G1, xlo, xhi, mu, sg, yX, yP);
  sd_axes(G1, xlo, xhi, ymx);
  // a rug: each mean as a short magenta tick on the axis
  FOR ix FROM 1 TO 40 DO
    px := LFT + IP((ms(ix) - xlo)/(xhi - xlo)*W);
    IF px > LFT AND px <= RGT THEN LINE_P(G1, px, BOT - 3, px, BOT - 1, RGB(200,0,200)); END;
  END;
  // the forty means against the theory: their mean beside μ, their SD
  // (n - 1 divisor) beside σx-bar = σ/√n
  sm := 0;
  FOR ix FROM 1 TO 40 DO sm := sm + ms(ix); END;
  mbar := sm/40;
  sm := 0;
  FOR ix FROM 1 TO 40 DO sm := sm + (ms(ix) - mbar)^2; END;
  sd2 := (sm/39)^0.5;
  s1 := "mean of the 40 " + CHAR(57344) + " = " + STRING(mbar,2,3) +
        "   (μ = " + STRING(mu,2,2) + ")";
  s2 := "SD of the 40 " + CHAR(57344) + " = " + STRING(sd2,2,3) +
        "   (σ" + CHAR(57344) + " = " + STRING(se,2,4) + ")";
  s3 := "one draw: σ = " + STRING(sg,2,2) + ";  a mean of " +
        STRING(nn,2,0) + ": σ" + CHAR(57344) + " = " + STRING(se,2,4);
  s4 := "bars one σ" + CHAR(57344) + " wide, total area 1";
  IF drct THEN s4 := "n > 100: each mean drawn from N(μ, σ" + CHAR(57344) + ")"; END;
  TEXTOUT_P("40 means of n = " + STRING(nn,2,0) + " draws from N(" +
            STRING(mu,2,2) + ", " + STRING(sg,2,2) + ")", G1, 6, 4, 1, RGB(0,60,200), 308);
  RECT_P(G1, 4, 166, 315, 216, RGB(255,69,0), RGB(255,240,220));
  TEXTOUT_P(s1, G1, 8, 169, 1, RGB(200,0,200), 305);
  TEXTOUT_P(s2, G1, 8, 180, 1, RGB(200,0,200), 305);
  TEXTOUT_P(s3, G1, 8, 191, 1, RGB(0,140,60), 305);
  TEXTOUT_P(s4, G1, 8, 202, 1, RGB(90,90,90), 305);
  TEXTOUT_P("any key or tap = menu", G1, 6, 226, 1, RGB(30,144,255));
  BLIT_P(G0, 0, 0, G1);
  sd_WaitAny();
END;

// ---- Help: two drawn pages, each held by a key OR a tap ----
sd_showHelp()
BEGIN
  LOCAL hK, hBk, hBl, hMg;
  hBk := RGB(0,0,0); hBl := RGB(0,60,200); hMg := RGB(200,0,200);

  // ---- page 1: the idea and the menu items
  DIMGROB_P(G1,320,240);
  RECT_P(G1,0,0,319,239,RGB(255,255,255));
  TEXTOUT_P("Sampling Distribution — Help (1 of 2)", G1, 6, 4, 1, hBl);
  LINE_P(G1, 6, 17, 313, 17, RGB(90,90,90));
  TEXTOUT_P("The sample mean " + CHAR(57344) + " of n draws from a", G1, 6, 24, 1, hBk);
  TEXTOUT_P("Normal(μ, σ) population is itself Normal,", G1, 6, 35, 1, hBk);
  TEXTOUT_P("centred on μ with standard error", G1, 6, 46, 1, hBk);
  TEXTOUT_P("σ" + CHAR(57344) + " = σ / SQRT(n).", G1, 6, 57, 1, hMg);

  // the menu items by name (the Prime numbers them on the menu)
  TEXTOUT_P("Find probability", G1, 6, 74, 1, hBl);
  TEXTOUT_P("   Choose left, right, between or outside,", G1, 6, 85, 1, hBk);
  TEXTOUT_P("   then enter the bounds.", G1, 6, 96, 1, hBk);
  TEXTOUT_P("Find critical value", G1, 6, 110, 1, hBl);
  TEXTOUT_P("   Enter α; get the " + CHAR(57344) + " cutting off that tail.", G1, 6, 121, 1, hBk);
  TEXTOUT_P("Compare sample sizes", G1, 6, 135, 1, hBl);
  TEXTOUT_P("   Three curves, shared axes: 4× n halves σ" + CHAR(57344) + ".", G1, 6, 146, 1, hBk);
  TEXTOUT_P("Watch sample means", G1, 6, 160, 1, hBl);
  TEXTOUT_P("   40 simulated means fall onto the axis,", G1, 6, 171, 1, hBk);
  TEXTOUT_P("   then become a histogram under the curve.", G1, 6, 182, 1, hBk);

  RECT_P(G1, 4, 196, 315, 214, RGB(255,69,0), RGB(255,240,220));
  TEXTOUT_P("Results persist in Home as EP_SD_*", G1, 8, 200, 1, hBk);
  TEXTOUT_P("any key or tap = next", G1, 6, 226, 1, RGB(30,144,255));
  BLIT_P(G0, 0, 0, G1);
  hK := sd_WaitAny();

  // ---- page 2: the live graph and the comparison view
  DIMGROB_P(G1,320,240);
  RECT_P(G1,0,0,319,239,RGB(255,255,255));
  TEXTOUT_P("Sampling Distribution — Help (2 of 2)", G1, 6, 4, 1, hBl);
  LINE_P(G1, 6, 17, 313, 17, RGB(90,90,90));
  TEXTOUT_P("The result graph", G1, 6, 22, 1, hBl);
  TEXTOUT_P("   grey = the population N(μ, σ); blue = " + CHAR(57344) + ",", G1, 6, 33, 1, hBk);
  TEXTOUT_P("   on ONE scale: at n = 30 the population's", G1, 6, 44, 1, hBk);
  TEXTOUT_P("   peak is 18% of " + CHAR(57344) + "'s.", G1, 6, 55, 1, hBk);
  TEXTOUT_P("Soft keys", G1, 6, 69, 1, hBl);
  TEXTOUT_P("   n" + CHAR(8722) + " n+  step n through 1, 2, 4 ... 400", G1, 6, 80, 1, hBk);
  TEXTOUT_P("   n" + CHAR(215) + "4  quadruples n: σ" + CHAR(57344) + " halves", G1, 6, 91, 1, hBk);
  TEXTOUT_P("   μ,σ  edit;  Tail  next region;  Menu", G1, 6, 102, 1, hBk);
  TEXTOUT_P("Touch", G1, 6, 116, 1, hBl);
  TEXTOUT_P("   drag in the plot: the nearer bound moves", G1, 6, 127, 1, hBk);
  TEXTOUT_P("   " + CHAR(8646) + " nudge it, " + CHAR(8645) + " the other bound, +/- = n,", G1, 6, 138, 1, hBk);
  TEXTOUT_P("   Ent = more, Esc = menu", G1, 6, 149, 1, hBk);
  TEXTOUT_P("Compare sample sizes", G1, 6, 163, 1, hBl);
  TEXTOUT_P("   drag = a cursor with P(" + CHAR(57344) + " < x) for each n;", G1, 6, 174, 1, hBk);
  TEXTOUT_P("   tap a curve to pick it.", G1, 6, 185, 1, hBk);
  TEXTOUT_P("Critical values: a new n moves them (α kept).", G1, 6, 199, 1, hMg);
  TEXTOUT_P("any key or tap = menu", G1, 6, 226, 1, RGB(30,144,255));
  BLIT_P(G0, 0, 0, G1);
  hK := sd_WaitAny();
END;

// ---- entry point ----
EXPORT SampDist()
BEGIN
  LOCAL modeM;

  REPEAT
    // Reset INSIDE the loop. CHOOSE does not touch its index on CANCEL, so
    // initialising once before the loop would let the second pass silently
    // re-run the first pass's choice.
    modeM := 0;
    // No numbers in the strings: the Prime numbers CHOOSE items itself
    // (1-9, then A, B, C), so a "1." here showed as a double number.
    CHOOSE(modeM, "Sampling Distribution",
      {"Find probability",
       "Find critical value",
       "Compare sample sizes",
       "Watch sample means",
       "Help"});
    CASE
      IF modeM == 1 THEN sd_modeProb(); END;
      IF modeM == 2 THEN sd_modeCrit(); END;
      IF modeM == 3 THEN sd_modeCompare(); END;
      IF modeM == 4 THEN sd_watch(); END;
      IF modeM == 5 THEN sd_showHelp(); END;
    END;
  UNTIL modeM == 0;

  RETURN "EP_SD_SE, EP_SD_P, EP_SD_ZLO, EP_SD_ZHI, EP_SD_LO, EP_SD_HI";
END;

// ============================================================
// DEVICE CHECKLIST — untested on hardware, check these first
// (the full ordered list is in SampDist_v2.3_CHANGES.md)
// ============================================================
//  1. CHOOSE items built by concatenation: the Mode 1 tail menu builds its
//     four item strings with "…" + CHAR(57344) + "…" inside the list literal.
//     Every other CHOOSE in the fleet uses plain literals. If the menu comes
//     up blank or errors, replace CHAR(57344) with a plain "xbar" there.
//  2. ₃ in the "n₃:" INPUT label. Only ₀ ₁ ₂ are on the verified subscript
//     list; ₃ is assumed to come from the same HP palette but has not been
//     seen on screen. If it renders as a box, drop all three to n1/n2/n3.
//  3. Tick centring by 3*SIZE(label) — cosmetic only; check labels sit under
//     their ticks at both a wide (σ = 100) and a narrow (σ = 0.5) scale.
//  4. Panel width at the extremes: the Mode 1 "Between" result line is the
//     longest string on screen; every panel line now passes a 305-px width,
//     so a long line is clipped rather than run off the edge.
//  5. The soft keys (v2.3): DRAWMENU with labels built by CHAR(8722),
//     CHAR(215), CHAR(956), CHAR(963) concatenation; a tap on each slot.
//  6. The live bound (v2.3): drag inside the plot box, lift, check the snap.
//  7. G2 as a scene buffer (v2.3, Compare and Watch): BLIT_P(G1, 0, 0, G2)
//     is GaltonBoard's idiom, not yet run on hardware by any fleet program.
//  8. √ inside a TEXTOUT_P string (the panel's "/√30").  It is the Prime's own
//     keyboard glyph, but no device-run program has printed it yet.  If it
//     shows as a box, replace "/√" with "/sqrt(" and close the bracket.
//  9. CHAR(8646) / CHAR(8645) (the ⇆ ⇅ arrows) in the panel and help text:
//     Lessons Learned section 13 records both as rendering.
