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

// ============================================================
// STATDIST - v2.7e
// Probabilities, densities and moments for 27 distributions,
// each drawn as a plot you can touch
// ============================================================
// Author: Roger D. Metcalf, UT Arlington, 2026 (framework after
//   S. Micciche).  CRCJ 3340 - Statistics and Introduction to
//   Research Methods.  Written with AI assistance (Anthropic
//   Claude); concept, direction and all hardware testing by the
//   author.
//
// 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.  Provided as is, with no
//   warranty: check any result you intend to rely on.  The
//   results have not been rigorously verified or validated; do
//   not rely on them for anything important.
// DISCLAIMER: a teaching and demonstration aid. Not validated for
//   forensic, clinical or legal use; do not rely on its results
//   for legal purposes.
//
// Firmware: this version has NOT yet been run on hardware or on
//   the emulator.  On the HP Prime Virtual Calculator (build
//   2026-09-09, the calculators' own firmware): v2.6 ran in round 1
//   (16-17 Oct 2026; v2.7b fixed what it found), and v2.7b in round 2
//   (17 Oct 2026): F(2.5, 7.5) and chi-square(2.5) gave 0.800427 and
//   0.876911, and the Formula page showed the CAS's exact answers
//   (15/11, 405/77, 2880/847, "exact (CAS)").  The CAS never failed
//   there, so the fallback fixed in v2.7c has run only in the PC
//   harness (ppl_run_v3d).  v2.7e has run only in the PC harness
//   (ppl_run_v3e, 19 Oct 2026), not on the emulator or a calculator.
//   Written for the 2026-09-09 firmware (2.4.2); none of its commands
//   is new in 2.4.x, so older firmware should also work.
//
// WHAT CHANGED in v2.7e (19 Oct 2026): the right tail of every family
//   computed as a tail, and four smaller items -- D9, B1, B2, B5, B7
//   (fleet review 19 Oct 2026, findings_B).  Line numbers are v2.7d's.
//   - THE RIGHT TAIL, P(X > x), FAMILY BY FAMILY (D9, Roger's decision
//     of 19 Oct; it removes B3 and B4).  v2.7d took every right tail as
//     1 - F(x) (st_Prob 2514, st_CutCalc 2536-2538).  The Prime keeps
//     12 significant digits, so far out F(x) is 0.999999999999... and
//     1 - F(x) is exactly 0 once the tail is below about 5E-13: a far
//     right tail read "= 0.000000" where the same tail on the left read
//     "< 0.000001" (B4), and the Formula page, which checks the CAS's
//     exact answer against the program's number to 6 significant
//     digits, threw away right answers for discrete tails below about
//     1E-6 (Binomial(20, 0.1), P(X > 11): "no closed form shown", B3).
//     New st_sf (after st_cdf) is each family's own right tail, a lower
//     tail or a formula with nothing subtracted from 1 (the 27 formulas
//     of findings_B part I): the Normal mirrored, NORMALD_CDF(mu,
//     sigma, 2 mu - x); e^(-lambda x) for the Exponential; the Binomial
//     turned round, BINOMIAL_CDF(n, 1 - p, n - k - 1); and so on.  New
//     st_gammaQ (after st_gammaP) is the upper incomplete gamma Q(a, x)
//     for the chi-square, Gamma and Erlang, returned without the 1 - Q
//     that st_gammaP ends with.  Three of the 27 go elsewhere at sizes
//     the program's own loops cannot reach (neither form limits lambda
//     or the Erlang's k): above a shape of 5000 the Poisson's P(k + 1,
//     lambda) and the Erlang's Q(k, lambda x) come from the built-ins
//     CHISQUARE_CDF(2k + 2, 2 lambda) and POISSON_CDF(lambda x, k - 1),
//     the same numbers (500 passes are too few there); and the F keeps
//     v2.7d's 1 - F(x) when a df passes 10000 (st_betaI's log-gamma
//     terms lose digits on 12-digit reals there: about 3E-6 of the tail
//     at df 100000, measured in a PC model of the Prime's rounding).  The
//     Triangular below its mode is written so that nothing cancels when
//     the mode sits at the maximum.  Used by: st_Prob for P(X > x) (new
//     st_SfAt, beside st_CdfAt); a Between in the upper half, F(a) >
//     1/2, as S(a) - S(b) -- two small numbers -- instead of F(b) -
//     F(a) (Normal(0, 1), P(8 < X < 9) read "= 0.000000"); the "%
//     beyond >>" note (st_CutCalc); the panel, the Numbers page, the
//     summary returned to Home and the Formula page's check (all
//     through st_Prob).  The drag reads the right tail from a new cache
//     stSv, built in st_Cache: for bars, P(X > k) from the right, the
//     exact tail past the last bar plus the pmf bar by bar; for a curve,
//     st_sf where F(x) > 1/2 and 1 - F(x) elsewhere (exact there).  The
//     Negative binomial's LOWER tail (st_cdf 1364) had the same loss
//     (NegBin(30, 0.3), P(Y <= 0) = 2E-16 read "CDF: 0.000000" on the
//     Numbers page): it is now BINOMIAL_CDF(k + r, 1 - p, k), at most k
//     failures in k + r trials, the same number with no subtraction.
//   - B1: st_gammaP stops after 500 passes, too few from a shape of
//     about 10000 (Gamma(20000, 1) at its mean gave 0.5007, not
//     0.5009).  The Gamma form (2266) refuses k above 5000 and the
//     chi-square form (1896) df above 5000, with a box that says why;
//     their help lines say "0 < k <= 5000", "0 < df <= 5000".
//   - B2: the panel's "poor approximation" note (2883) was 58
//     characters, past the screen's edge at font 1; now "approx. =
//     0.7602: poor here, np or n(1-p) < 10" (46).
//   - B5: the Advanced menu (4208-4210) spells "Chi-square" and
//     "Maxwell-Boltzmann", as the plot and the Numbers page do.
//   - B7: a typed x below 0.0001 put a relation inside the label:
//     "P(X <= < 0.0001) = 0.500020".  New st_XStr (after st_fmtx) gives
//     the value's own digits instead ("P(X <= 0.00005)"); it is used for
//     x, a and b on the panel (2844, 2851, 2860), the Numbers page
//     (3230-3233, 3250) and the Formula page's "at x" line (3719).
//     Parameters keep v2.7d's "σ < 0.0001" (st_ParStr, by choice).
//   - The help page's title reads v2.7e; the file map names the new
//     routines.
//
// WHAT CHANGED in v2.7d (17 Oct 2026, evening): no "=" before a
//   relation.  Roger approved the pass, 17 Oct 2026, 21:18.  Line
//   numbers are v2.7c's.  Nothing else changed.
//   - st_fmtp writes "< 0.000001" (or "> -0.000001") below its floor,
//     and st_fmtx, for a value that needs more than three decimals,
//     "< 0.0001".  The displays put that after "= ": a far tail read
//     "P(X <= -6) = < 0.000001" on the panel, on the Numbers page and
//     in the summary returned to Home; a density "f(x) = < 0.0001"; a
//     sigma of 0.001 "variance = < 0.0001"; a sigma of 0.00005
//     "σ = < 0.0001".  New st_Eq (after st_fmtp) returns " = 0.158655",
//     or " < 0.000001" when the relation is already there.  The sites:
//     st_ParStr (1523-1550: each value carries its own " = "), the
//     panel (2777-2808, 2832), the Numbers page (3156-3196), the
//     Formula page's check lines (3573, 3653).
//   - The panel's normal approximation (Binomial, Poisson) took its
//     right tail as 1 - NORMALD_CDF (2821), which is exactly 0 once the
//     tail is below the CDF's last digit, where the left tail read
//     "< 0.0001".  It is now the lower tail at 2 mean - (x + 0.5), the
//     same number by symmetry.
//   - Seen, not changed: the right tail of every family, P(X > x), is
//     1 - CDF (st_Prob), so a far right tail can read "= 0.000000"
//     where the left reads "< 0.000001"; for a family with a bounded
//     range 0 is the true answer, so the floor cannot be applied there.
//     A fix needs each family's own upper tail; that is Roger's call.
//   - The help page's title reads v2.7d.
//
// WHAT CHANGED in v2.7c (17 Oct 2026): two string rules that the
//   second emulator round found, typed in Home.  Line numbers are
//   v2.7b's.  Nothing else changed.
//   - INSTRING("", "a") is "Error: Invalid input": the program
//     stops.  st_Cas (728-749) sets its answer to "" when the CAS
//     fails and then tested INSTRING(cR, "undef") (747), so a failed
//     CAS call would have stopped the Formula page instead of showing
//     the decimals and the "decimal" badge; the test now runs only
//     when DIM(cR) > 0.  st_Repl (780-797), which turns CAS text into
//     what the page shows, looked for the next match in what was left
//     of the text; when an answer ENDED in a name it replaces -- the
//     Pareto(1, 1) CDF (stqx-1)/stqx, a mean 1/stqa -- what was left
//     was "" and the next INSTRING stopped the page, with the CAS
//     answering.  It now looks only while text is left, and returns
//     an empty pattern's text unchanged.  The other INSTRINGs search a
//     fixed digit string (st_IVal 761, st_Digits 882, st_FmRes 3512)
//     or a text already known not to be empty (st_FmGood 3395-3403,
//     st_FmSub 3416); every LEFT (480, 741, 786, 814) and MID has a
//     count and a start of 1 or more.
//   - STRING(x, 2, d) prints in SCIENTIFIC form when 0 < |x| <
//     10^-d (STRING(0.096, 2, 1) = "9.6ᴇ-2").  st_fmtx (419-428)
//     printed a value within 1E-10 of a whole number with 0 decimals
//     as it was, so a mean of 1E-17 left by rounding error would read
//     "1ᴇ-17"; it now prints the rounded value (0, never "-0").
//     The other 15 fixed-decimal calls cannot reach the case: st_fmtp
//     shows a value below 10^-d as "< 0.000001"; the axis labels are
//     multiples of a 1-2-5 step with d = -FLOOR(LOG(step)); the
//     "% beyond" notes need 0.15 % or more; k, a and b of a discrete
//     family are whole; the "not a simple ratio" line (3706) shows 12
//     decimals, so only a parameter below 1E-12 would turn scientific,
//     the one form that can show it.
//   - The help page's title reads v2.7c.
//
// WHAT CHANGED in v2.7b (17 Oct 2026): what the first emulator round
//   found (Emulator_Round1_Results_2026-10-17_v1.md).  Line numbers are
//   v2.7a's.  Nothing else changed.
//   - A FRACTIONAL df NO LONGER STOPS THE PROGRAM.  On the emulator,
//     Advanced Continuous -> Chi-square with df = 2.5 stopped with
//     "Error: Invalid input": GfxProbe T11 and the Home checks showed
//     that CHISQUARE_CDF, CHISQUARE_ICDF, FISHER_CDF and FISHER_ICDF
//     take WHOLE-number df only (CHISQUARE_CDF(2.5, 3) is "Invalid
//     input"; SciPy gives 0.69415), while the densities and all of
//     Student's t take any df.  Now (st_cdf 1122-1129, st_pdf 997-1000):
//       chi-square, df not whole: F(x) = P(df/2, x/2), the program's own
//         regularized incomplete gamma (st_gammaP, already used for the
//         Gamma family);
//       F, either df not whole: F(x) = I_z(df1/2, df2/2) with
//         z = df1 x / (df1 x + df2), the program's own regularized
//         incomplete beta (st_betaI, already used for the Beta family);
//         and its density by the formula, through st_lngamma (FISHER's
//         density with a fractional df was not among the emulator's
//         checks, so it is not relied on);
//       whole df: the built-ins, exactly as before.
//     The chi-square quantile already used CHISQUARE_ICDF only for a
//     whole df; for any other df its bisection now calls the new CDF.
//     The F quantile always bisected.  The other calls were checked:
//     Erlang's CHISQUARE_CDF(2k, ...) gets a whole 2k (k is rounded
//     down at the form), Maxwell's CHISQUARE_ICDF(3, ...) a whole 3,
//     and POISSON_CDF is never asked for k < 0 (v2.5's guard; T12
//     showed POISSON_CDF(2, -1) is "Invalid input").
//   - SHORTER FORM LABELS.  The INPUT label "Degrees of freedom:" was
//     cut to "Degrees of fr..." on the emulator.  Every label longer
//     than about 14 characters is now short (the help line under each
//     field says the rest): "df:" for t and chi-square, "df1:" and
//     "df2:" for F, "N (pop.):", "K (pos. in N):" and "k (pos. in n):"
//     for the hypergeometric, and "Scale:" for the Cauchy (14
//     characters, but wide ones).
//   - The help page's title reads v2.7b.
//
// WHAT CHANGED in v2.7a (17 Oct 2026): the CAS call form, and the
//   one reader that waited for an infinity.  Line numbers are v2.7's.
//   - st_Cas (651-672): CASProbe beta 1b found that CAS(s) with s a
//     variable hands the text back UNEVALUATED, so v2.7's st_Cas
//     returned its own input; every Formula-page answer then failed
//     st_FmGood (it still said "int(") and the page showed "no closed
//     form shown" throughout.  Its IFERR line is now
//       STRING(EXPR("CAS(" + CHAR(34) + "string(" + s + ")" + CHAR(34) + ")"))
//     (BRIEF_CAS_v2), the quotes STRING() adds are dropped, and an
//     answer holding the infinity sign is refused like undef.
//   - "No mean" / "no variance" (st_FmCont, 3447 and 3468; st_FmInf,
//     3332-3340): the page used to wait for the CAS to answer the
//     divergent integral with an infinity ("+infinity", or the sign).
//     st_Cas now refuses any answer holding the infinity sign, so the
//     page asks a FINITE question instead: exp(-(integral)), which is
//     0 exactly when the integral is +infinity (a finite integral
//     gives a nonzero e^(-v); minus infinity gives an infinity, which
//     st_Cas refuses; an integral left undone stays unevaluated).
//     st_FmInf now reads that "0".
//   - st_FmGood (3292-3300) also refuses the integral and summation
//     signs (CHAR(8747), CHAR(8721), CHAR(931)), in case the Prime
//     writes an integral or a sum it cannot do with its sign rather
//     than as "integrate(" or "sum(" (not measured: CASProbe saw no
//     undone integral).  The value check would not catch it, since an
//     undone integral still evaluates to the right number.
//   - The help page's title reads v2.7a; comments that said the call
//     form was still to be decided say what CASProbe found.  Nothing
//     else changed.
//
// Run StatDist() from the Program Catalog.  All 27 distributions
//   are live.  Continuous: Normal, Uniform, Exponential, Weibull,
//   LogNormal, Laplace, Student t, Chi-square, Fisher F, Burr XII,
//   Gumbel, Pareto, Rayleigh, Logistic, Triangular, Cauchy,
//   Frechet, Gompertz, Maxwell-Boltzmann, Erlang, Gamma, Beta.
//   Discrete: Binomial, Poisson, Geometric, Negative binomial,
//   Hypergeometric.
//
// What changed in v2.7 (16 Oct 2026; the CAS version, built on
//   v2.6 with every v2.6 change kept; line numbers are v2.6's):
//   * A FORMULA PAGE for the distribution on the plot: soft-key
//     page 3 (More, More) -> Formula, or Enter -> "Formula page
//     (CAS)".  The CAS does the calculus with YOUR numbers:
//     - continuous: E[X] = integral of x f(x), E[X^2] = integral
//       of x^2 f(x), and Var(X) = E[X^2] - E[X]^2, as exact
//       expressions (Exponential(1/2): 2, 8, 4; Normal(100, 15):
//       100, 10225, 225; Weibull(2, 1): sqrt(pi)/2 ...);
//     - for the Exponential, Weibull and Pareto, the CDF itself,
//       F(x) = integral of f(t) from the lowest value to x;
//     - discrete (Binomial, Poisson, Geometric, Negative binomial,
//       Hypergeometric, at moderate sizes): the probability on the
//       plot as an exact fraction (Binomial(20, 1/2), P(X <= 10) =
//       308333/524288), the red bar's P(X = k), and for the
//       Binomial and Hypergeometric the mean as an exact sum;
//     - a moment that does not exist (Cauchy, t with df <= 1,
//       Pareto with alpha <= 1 ...) is shown as an integral the CAS
//       finds infinite.
//     The parameters go to the CAS as exact ratios (0.085 is
//     85/10^3).  One that is not a simple ratio (more than six
//     decimals, such as 1/3 typed as 0.333333333333) stays a SYMBOL,
//     so the CAS gives the formula, and the decimal is the
//     program's own.
//   * NEVER A WRONG FORMULA: every CAS answer is checked against
//     the program's own value (st_Mom, st_cdf, st_Prob) before it is
//     shown; an answer that fails the check, or a CAS that fails,
//     shows "no closed form shown" and the program's decimal.  In
//     desktop Giac the check rejects, for example, Logistic E[X^2]
//     (Giac says -pi^2/3: the sign is wrong) and leaves out Gumbel,
//     Frechet and Gompertz moments (Giac cannot do the integral).
//   * The fleet's CAS fence: ONE call routine (st_Cas), whole
//     numbers only in CAS text (st_IStr), tagged CAS variables
//     (stqx, stqt, stqk, stqa, stqb, stqc), the CAS asked only when
//     the Formula page OPENS and only if the plot changed since the
//     last time (the page is kept), a badge "exact (CAS)" or
//     "decimal".  The plot, the drag and the caches are untouched:
//     no frame and no drag ever calls the CAS.
//   * Soft keys: "More" now turns THREE pages; page 3 is {Formula,
//     Help, , , Menu, More}.  The Enter options gain "Formula page
//     (CAS)".  Help covers both.
//
// What changed in v2.6 (16 Oct 2026; built on the v2.5 fix
//   version, and every v2.5 fix is kept -- see the v2.5 list):
//   * Every result is now a DRAWN PLOT (st_Plot): the density in
//     blue over a nice-number grid, the requested probability
//     shaded column by column, x as a red line, the mean as a
//     green dotted line, and an orange-framed panel with the
//     parameters, mean, sd, the density at x ("a density, not a
//     probability") and the probability in magenta.  Discrete
//     distributions draw pmf bars: the counted bars filled, the
//     chosen k in red (an outline only when k itself is NOT
//     counted, as in P(X > k) -- the off-by-one made visible).
//   * The window fits itself to the 0.1 % .. 99.9 % range; heavy
//     tails (Cauchy, Pareto, small-df t ...) are cut at 8
//     interquartile ranges and the plot says how much lies beyond.
//   * TOUCH: tap or drag in the plot to move x (or k); the shading
//     and the readout follow the finger from a cache, and x snaps
//     to a nice grid when the finger lifts.  Soft keys: P<=x,
//     P>x, a<X<b (then tap a and b; later a drag moves the nearer
//     line), Zoom-, Zoom+ (about x) and More.
//     Slot 6 is "More" (the house convention of 16 Oct 2026): it
//     turns to a second page -- Edit, Numbers, CDF, Approx / N(0,1)
//     / Refit, Menu -- so everything works by touch alone.
//     Arrows move x one grid step; Up/Down pick the upper/lower
//     tail (or, in a<X<b, which bound the arrows move); +/- zoom;
//     0 refits; Enter opens the options; Esc returns to the menu.
//   * Enter -> options: the v2.5 NUMBERS PAGE (now built from the
//     live plot, so it follows a dragged x), the CDF as a magenta
//     S-curve on a right-hand 0..1 axis, the normal approximation
//     over Binomial and Poisson bars (with "poor here" when np,
//     n(1-p) or lambda is under 10), N(0,1) in pale grey behind
//     Student t, Between, Edit, refit and help.
//   * The ten v2.5 stubs are implemented: Triangular, Cauchy,
//     Frechet, Gompertz, Maxwell-Boltzmann, Erlang, Gamma (any
//     k > 0: our own regularized incomplete gamma -- series, then
//     a Lentz continued fraction), Beta (our own incomplete beta,
//     Lentz with the symmetry swap), Negative binomial (through
//     BINOMIAL_CDF) and Hypergeometric (COMB in a loop).  Checked
//     against SciPy in verify_StatDist.py.
//   * A "typical use" line for every distribution, chosen for
//     criminal-justice and forensic students, in the INPUT title,
//     on the plot and on the numbers page.
//   * GALLERY, the new first menu item: six thumbnails (Normal,
//     Exponential, Chi-square(4), t(3) over N(0,1), Binomial(10,
//     0.3), Poisson(3)); tap one to open it with those parameters.
//   * Help pages cover every gesture and soft key; st_Page now
//     waits for a key OR a tap (st_WaitAny, with the capped
//     finger-lift drain) and its footer is in the house form.
//   * Structure: the menus call st_Run(id); one dispatch section
//     computes every distribution (st_pdf, st_cdf, st_icdf,
//     st_Mom ...); one input section per distribution (st_ask_*);
//     EXPORT StatDist() moved to the end so every routine is
//     defined before it is used.  The v2.5 st_calc_* routines are
//     gone (their forms live on in st_ask_*, their result text in
//     st_Numbers).  New code writes ^0.5 and number literals, not
//     the root sign and pi.
//   * Small corrections found on the way: Student t's mean is
//     "undef" when df <= 1 (v2.5 printed 0); the a/b bounds of
//     Between are put in order when typed backwards.
//
// What changed in v2.5 (16 Oct 2026, the fix version; line
//   numbers are v2.4's -- all of this is kept in v2.6):
//   * Every result went to MSGBOX as 8-10 lines with the
//     probability LAST, and MSGBOX shows about five lines (Lessons
//     Learned s3).  Results are drawn pages (st_Page): orange title
//     band, font 2 at 14 px, the "P(...) = ..." line in magenta.
//   * The ten stubs (787-796) became st_calc_* (tagged names).
//   * The two header blocks became one; #pragma is line 1.
//   * Poisson "between" with lower k = 0 called POISSON_CDF(l, -1)
//     (629).  v2.6 keeps the fix: st_cdf returns 0 for k < 0 and
//     never asks POISSON_CDF for a negative k.
//   * n and k typed as fractions are rounded down with a one-line
//     box ("n rounded down to 20"), now in st_Int.
//   * The tail type is asked first, then only the fields that tail
//     needs (x for the tails, a and b for Between).
//   * Tagged LOCALs (stC, stK, stS, stN, stP, stZ ...).
//   * Weibull stays on both the Common and the Advanced menu.
// v2.4 (10 Oct 2026): every helper routine carries the tag "st_";
//   all programs share ONE pool of routine names on the Prime.
// v2.3 (29 Sep 2026): CHOOSE index set to 0, called, then tested.
// v2.2: every displayed number goes through a formatter (st_fmtp,
//   st_fmtv); "undef ..." strings for moments that do not exist.
// v2.1: CHOOSE argument order fixed; LOCALs hoisted out of IF;
//   Burr closed-form CDF and Lanczos moments; Geometric by hand.
//
// MAP OF THIS FILE (top to bottom)
//   1. shared state (file-level LOCALs)
//   2. small helpers: st_Init, st_Drain, st_WaitAny, formatters
//      (st_fmtp, st_Eq, st_fmtv, st_fmtx, st_XStr), st_Nice, st_Int,
//      st_Page
//   3. special functions: st_lngamma, st_atan, st_tanr, st_gammaP,
//      st_gammaQ, st_betacf, st_betaI, st_hyp, st_lapCDF, st_gompMom
//   3b. the CAS (v2.7): st_IStr, st_Cas, st_IVal, st_Repl, st_Wrap,
//      st_Scale, st_CasNum, st_Rat, st_RTxt, st_LTxt, st_Digits
//   4. THE DISPATCH -- every distribution's maths in one place:
//      st_Name, st_Use, st_Disc, st_Def, st_Sup, st_pdf, st_cdf,
//      st_sf (the right tail, v2.7e), st_qbis, st_dq, st_icdf,
//      st_Mom, st_ParStr, st_DefAB, st_ABHelp
//   5. the forms: st_get_tail_type, st_AskAB, one st_ask_* per
//      distribution (with its typical use), st_Ask
//   6. help: st_show_help
//   7. the plot: window, caches, drawing, touch, options, st_Plot,
//      st_Run; the Formula page (v2.7): st_Fm* and st_Formula
//   8. the gallery
//   9. the menus and EXPORT StatDist()
//
// THE DISTRIBUTION NUMBERS (the "id" used everywhere)
//    1 Normal        2 Uniform       3 Exponential   4 Weibull
//    5 LogNormal     6 Laplace       7 Student t     8 Chi-square
//    9 Fisher F     10 Burr XII     11 Gumbel       12 Pareto
//   13 Rayleigh     14 Logistic     15 Binomial     16 Poisson
//   17 Geometric    18 Triangular   19 Cauchy       20 Frechet
//   21 Gompertz     22 Maxwell      23 Erlang       24 Gamma
//   25 Beta         26 Neg.binomial 27 Hypergeometric
//   15, 16, 17, 26 and 27 are discrete.
// ============================================================


// ============================================================
// 1. SHARED STATE
// Everything the plot, the touch handler and the gallery share
// lives here at file level, so the drawing code and the touch
// code read the SAME numbers (Lessons Learned s6: keep the plot
// box in shared constants so the two can never disagree).
// ============================================================
// ---- the plot box: x from stPL to stPR, y from stPT down to stPB ----
LOCAL stPL, stPR, stPT, stPB, stW, stH;
// ---- the house colours (set once in st_Init) ----
LOCAL stBLU, stMAG, stRED, stGRN, stORN, stGRY;
LOCAL stGRID, stSHD, stFOOT, stLORN, stPALE, stLast;
// ---- the distribution on the plot ----
//  stId    which distribution (1..27, the table above)
//  stPar   its parameters packed in ONE list {p1, p2, p3}: a PPL
//          call takes at most 16 arguments, so every dispatch
//          routine takes the list instead of loose numbers
//  stTl    the tail: 1 = P(X <= x), 2 = P(X > x), 3 = P(a < X < b)
//  stXv    x (or k);  stAv, stBv  the bounds a and b
//  stPh    a<X<b set-up: 1 = the next tap sets a, 2 = it sets b
//  stAct   which bound the arrows move: 1 = a, 2 = b
//  stDrg   1 while a finger is dragging (readouts come from cache)
//  stDisc  1 for a discrete distribution;  stZm 1 once zoomed
LOCAL stId, stPar, stTl, stXv, stAv, stBv;
LOCAL stPh, stAct, stDrg, stDisc, stZm;
// ---- the window and the caches built from it ----
//  stLo..stHi  the x range on screen (discrete: first and last k)
//  stGs        the grid step: arrows move x by it, a lifted
//              finger snaps x onto it;  stTk the tick spacing
//  stPd        density (or pmf) at every pixel column (or bar)
//  stCd        the CDF at the same places
//  stPy        the pixel height of stPd, so a frame does no maths
//  stNy        pixel heights of the overlay curve (N(0,1) behind
//              t, or the normal approximation over the bars)
//  stYmx       the top of the y axis;  stCut {left, right}: how
//              much probability lies outside the window
//  stMu, stVr  mean and variance (a STRING when they do not exist)
LOCAL stLo, stHi, stGs, stTk, stYmx, stCut;
LOCAL stPd, stCd, stPy, stNy, stMu, stVr;
//  stSv  (v2.7e) the RIGHT tail P(X > x) at the same places as stCd,
//        so a drag reads P(X > x) from a list as it reads F(x); it is
//        built with no "1 - F(x)" where F(x) is near 1 (st_sf, below)
LOCAL stSv;
// ---- overlays, switched from the Enter menu or soft-key page 2;
//      stSKP is the soft-key page showing (1 or 2, "More" turns it) ----
LOCAL stOvN, stOvT, stOvC, stClip, stSKP;
// ---- the gallery: thumbnail heights, tap boxes, selection ----
LOCAL stGH, stGX, stGY, stGSel;
// ---- v2.7: the CAS and the Formula page ----
//  stCasOK  1 while the Formula page being built has shown a CAS
//           answer (each one checked); its badge then reads
//           "exact (CAS)", else "decimal"
//  stFmL, stFmC  the page's lines and their colours, KEPT: the CAS
//           is asked again only when stFmK -- the plot's state when
//           the page was built -- no longer matches the plot
//  stFmOK   the kept page's badge
//  stFmSy   for p1, p2, p3: 1 = not a simple ratio, so a symbol
//           (stqa, stqb, stqc) in the CAS text;  stFmR their values
//           as 12-digit ratios, for checking;  stFmXR x the same way;
//           stFmPT the three parameter texts handed to the CAS
LOCAL stCasOK, stFmL, stFmC, stFmK, stFmOK, stFmSy;
LOCAL stFmR, stFmXR, stFmPT;


// ============================================================
// 2. SMALL HELPERS
// ============================================================

// The house palette and the plot geometry, set once at the start.
st_Init()
BEGIN
  stBLU := #003CC8h;   // blue    RGB(0,60,200)    the distribution
  stMAG := #C800C8h;   // magenta RGB(200,0,200)   probability, CDF
  stRED := #FF0000h;   // red     RGB(255,0,0)     "you are here": x
  stGRN := #008C3Ch;   // green   RGB(0,140,60)    the mean
  stORN := #FF4500h;   // orange  RGB(255,69,0)    frames
  stGRY := #5A5A5Ah;   // grey    RGB(90,90,90)    tick labels
  stGRID := #E1E4EBh;  // grid    RGB(225,228,235)
  stSHD := #AACDFFh;   // shade   RGB(170,205,255) the probability
  stFOOT := #1E90FFh;  // dodger blue               key legends
  stLORN := #FFF0DCh;  // light orange              the panel
  stPALE := #A8A8A8h;  // pale grey                 N(0,1) behind t
  // The plot box: 274 pixel columns (40..314) and 121 rows (17..138).
  // Below it: tick labels at 142, the panel from 152 to 217, and the
  // soft-key strip from 218 to 239.
  stPL := 40; stPR := 314; stPT := 17; stPB := 138;
  stW := stPR - stPL; stH := stPB - stPT;
  stOvN := 0; stOvT := 0; stOvC := 0; stClip := 0;
  stDrg := 0; stZm := 0; stPh := 0; stAct := 2;
  stCut := {0, 0};
  stLast := "no plot drawn yet";
  stGX := MAKELIST(0, K, 1, 6, 1); stGY := MAKELIST(0, K, 1, 6, 1);
  stGSel := 1;
  stSKP := 1;
  // v2.7: no Formula page kept yet
  stCasOK := 0; stFmL := {}; stFmC := {}; stFmK := {}; stFmOK := 0;
  stFmSy := {0, 0, 0}; stFmR := {"(0)", "(0)", "(0)"}; stFmXR := "(0)";
  stFmPT := {"(0)", "(0)", "(0)"};
END;

// 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 (Lessons Learned s7).  The ChiTable shape, unchanged.
st_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 keys and touch, wait, drain.
// A GETKEY drain alone does not drain touch, so the finger lift
// that opened a page could close it (Lessons Learned s7).
st_WaitAny()
BEGIN
  LOCAL waK, waM, waM1;
  REPEAT waK := GETKEY; UNTIL waK == -1;
  st_Drain();
  REPEAT
    waK := GETKEY;
    waM := MOUSE;
    waM1 := waM(1);
  UNTIL waK > -1 OR SIZE(waM1) > 0;
  st_Drain();
  REPEAT waK := GETKEY; UNTIL waK == -1;
END;

// ---- probability / density formatter ----
// v2.2: every result used to be concatenated raw ("  PDF: " + pdf_val), which
// hands formatting to whatever number format the user happens to have set in
// Home -- so the same calculation printed differently on two calculators, and
// often as a 12-digit smear.  st_fmtp() pins the width AND guards the floor: a
// value below the displayable resolution prints "< 0.000001" instead of
// "0.000000", which reads as exactly zero.  Params pv/dg, not p/d -- short
// parameter names are where builtin collisions bite.  Non-EXPORT.
st_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;

// ---- " = value", or " < 0.000001" (v2.7d) ----
// st_fmtp writes the relation itself below its floor ("< 0.000001"),
// so "P(X <= -6) = < 0.000001" said two things at once; "P(X <= -6)
// < 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).  Param eqS: short names are where built-ins collide.
st_Eq(eqS)
BEGIN
  IF DIM(eqS) > 0 THEN
    IF LEFT(eqS, 1) == "<" OR LEFT(eqS, 1) == ">" THEN RETURN " " + eqS; END;
  END;
  RETURN " = " + eqS;
END;

// ---- polymorphic value formatter (numbers OR "undef" strings) ----
// A moment holds a STRING when it does not exist: "undef (K*C<=1)" for
// Burr, "undef (df<=2)" for Student, "undef (no mean)" for Cauchy.
// Formatting those must pass the string through untouched, so dispatch on
// TYPE() first.  TYPE(x)==0 means a real number.
st_fmtv(vv, dg)
BEGIN
  IF TYPE(vv) != 0 THEN RETURN vv; END;
  RETURN st_fmtp(vv, dg);
END;

// ---- a short number for the plot (v2.6) ----
// The fewest decimals (0 to 3) that show the value exactly, else four:
// 2 -> "2", 0.3 -> "0.3", 1.96 -> "1.96", 1/3 -> "0.3333".  A typed
// x = 1.96 must never be shown rounded to "2.0".
// v2.7c: the value is ROUNDed before STRING: a value within 1E-10 of a
// whole number, such as a mean of 1E-17 left by rounding error, printed
// as "1ᴇ-17" (fixed decimals turn scientific below 10^-d on the Prime,
// emulator 17 Oct 2026); rounded, it prints 0.
st_fmtx(fxV)
BEGIN
  LOCAL fxD, fxR;
  IF TYPE(fxV) <> 0 THEN RETURN fxV; END;
  FOR fxD FROM 0 TO 3 DO
    fxR := ROUND(fxV, fxD);
    IF ABS(fxR - fxV) < 0.0000000001*MAX(1, ABS(fxV)) THEN
      IF fxR == 0 THEN fxR := 0; END;     // a rounded -1E-17 is 0, never "-0"
      RETURN STRING(fxR, 2, fxD);
    END;
  END;
  RETURN st_fmtp(fxV, 4);
END;

// ---- x, a or b inside a label such as "P(X <= x)" (v2.7e) ----
// st_fmtx and st_fmtp write "< 0.0001" for a value too small for their
// decimals.  That is right for a probability, but x is a VALUE, not a
// probability: a typed x = 0.00005 read "P(X <= < 0.0001) = 0.500020"
// (fleet review 19 Oct 2026, B7).  xsS is the text the page would
// print (from st_fmtx or st_fmtp), xsV the value itself.  A text that
// starts with a relation is replaced by the value's own digits: the
// fewest decimals, 5 to 11, that show it exactly (0.00005 -> 0.00005),
// else four significant digits.  The value is ROUNDed before STRING:
// fixed decimals turn scientific below 10^-d on the Prime (Lessons
// Learned s3), and a rounded value is either 0 or at least 10^-d.
// Any other text comes back unchanged.  Params xsV, xsS: tagged.
st_XStr(xsV, xsS)
BEGIN
  LOCAL xsD, xsR;
  IF DIM(xsS) > 0 THEN
    IF LEFT(xsS, 1) <> "<" AND LEFT(xsS, 1) <> ">" THEN RETURN xsS; END;
  END;
  IF xsV == 0 THEN RETURN "0"; END;
  // the fewest decimals that hold the value exactly
  FOR xsD FROM 5 TO 11 DO
    xsR := ROUND(xsV, xsD);
    IF ABS(xsR - xsV) < 0.0000000001*ABS(xsV) THEN RETURN STRING(xsR, 2, xsD); END;
  END;
  // otherwise four significant digits (at most 11 decimals)
  xsD := MIN(11, 3 - FLOOR(LOG(ABS(xsV))));
  xsR := ROUND(xsV, xsD);
  RETURN STRING(xsR, 2, xsD);
END;

// A "nice" step (1, 2 or 5 times a power of ten) NEAR the value
// handed in -- Heckbert's rule, for axis ticks and the x grid.
// (The rounding branch of Norm3's NiceNum.)
st_Nice(nnV)
BEGIN
  LOCAL nnE, nnF;
  IF nnV <= 0 THEN RETURN 1; END;
  nnE := FLOOR(LOG(nnV));
  nnF := nnV / 10^nnE;
  IF nnF < 1.5 THEN RETURN 10^nnE; END;
  IF nnF < 3 THEN RETURN 2*10^nnE; END;
  IF nnF < 7 THEN RETURN 5*10^nnE; END;
  RETURN 10^(nnE + 1);
END;

// A whole number from a typed value.  v2.5: a fraction is still
// rounded down with IP(), but a one-line box now says so.
st_Int(rdV, rdNm)
BEGIN
  IF FP(rdV) > 0 THEN
    MSGBOX(rdNm + " rounded down to " + STRING(IP(rdV), 2, 0));
  END;
  RETURN IP(rdV);
END;

// ---- drawn result page (v2.5; v2.6 waits for a key OR a tap) ----
// MSGBOX shows about five lines and simply drops the rest (Lessons
// Learned s3).  st_Page draws the lines on a page of its own instead:
//   pgTtl  the title, white on an orange band
//   pgLst  a list of strings, one per line, drawn in font 2 at 14 px
//          from y = 30; the last line is magenta when it is the
//          "P(...) = ..." answer; "" leaves a blank line
// Up to 13 lines fit above the footer; a longer list is paged.  The
// page is composed in G1 and blitted once (no tearing).
st_Page(pgTtl, pgLst)
BEGIN
  LOCAL pgN, pgI, pgY, pgCol, pgS;
  LOCAL pgFirst, pgLast;
  pgN := SIZE(pgLst);
  pgFirst := 1;
  REPEAT
    pgLast := MIN(pgN, pgFirst + 12);              // at most 13 lines a page
    DIMGROB_P(G1, 320, 240, #FFFFFFh);              // white page
    RECT_P(G1, 0, 0, 319, 21, #FF4500h, #FF4500h);  // orange title band
    TEXTOUT_P(pgTtl, G1, 6, 4, 2, #FFFFFFh, 308);   // white title, font 2
    pgY := 30;
    FOR pgI FROM pgFirst TO pgLast DO
      pgS := pgLst(pgI);
      pgCol := #000000h;                            // black text
      IF pgI == pgN AND LEFT(pgS, 2) == "P(" THEN
        pgCol := #C800C8h;                          // magenta: the answer line
      END;
      IF DIM(pgS) > 0 THEN
        TEXTOUT_P(pgS, G1, 6, pgY, 2, pgCol, 308);
      END;
      pgY := pgY + 14;
    END;
    // the footer, in the house form (comma-separated, spaces round =)
    IF pgLast < pgN THEN
      TEXTOUT_P("any key or tap = next page", G1, 6, 224, 1, #1E90FFh, 308);
    ELSE
      TEXTOUT_P("any key or tap = back", G1, 6, 224, 1, #1E90FFh, 308);
    END;
    BLIT_P(G0, 0, 0, G1);                           // one blit
    st_WaitAny();                                   // drain, wait, drain
    pgFirst := pgLast + 1;
  UNTIL pgFirst > pgN;
END;


// ============================================================
// 3. SPECIAL FUNCTIONS
// Pure arithmetic, built only from commands that have run on the
// device (EXP, LN, ATAN, TAN, COMB, NORMALD_CDF ...).
// ============================================================

// ---- log-gamma (Lanczos) for moments, Gamma and Beta; pure numeric ----
// LNGAMMA does not exist outside the CAS (Lessons Learned s3).
// v2.5: the parameter and the locals carry the st tag.
st_lngamma(stX)
BEGIN
  LOCAL stY, tmp, ser, stJ, cof;
  cof := {76.18009172947146,-86.50532032941677,24.01409824083091,-1.231739572450155,0.001208650973866179,-0.000005395239384953};
  stY := stX; tmp := stX + 5.5; tmp := tmp - (stX + 0.5)*LN(tmp);
  ser := 1.000000000190015;
  FOR stJ FROM 1 TO 6 DO stY := stY + 1; ser := ser + cof(stJ)/stY; END;
  RETURN -tmp + LN(2.5066282746310005*ser/stX);
END;

// ---- arctangent in RADIANS, whatever the angle mode ----
// ATAN answers in the calculator's angle mode: radians, degrees or
// grads.  ATAN(1) is a quarter turn in that same unit (0.785398...,
// 45 or 50), so ATAN(z) * (pi/4) / ATAN(1) is always radians.  The
// Marsaglia rule of Lessons Learned s11 for the same reason: never
// let the angle mode leak into a probability.
st_atan(atZ)
BEGIN
  RETURN ATAN(atZ)*0.785398163397448/ATAN(1);
END;

// ---- tangent of an angle given in RADIANS, whatever the mode ----
st_tanr(tnR)
BEGIN
  RETURN TAN(tnR*ATAN(1)/0.785398163397448);
END;

// ---- regularized lower incomplete gamma P(a, x) ----
// P(a, x) = (1/Gamma(a)) * integral from 0 to x of t^(a-1) e^(-t) dt.
// It is the Gamma distribution's CDF: Gamma(k, theta) has
// F(x) = P(k, x/theta).  Two methods, each where it converges fast
// (Numerical Recipes 6.2):
//   x < a + 1   the series  e^(-x) x^a / Gamma(a+1) * sum of
//               x^n / ((a+1)(a+2)...(a+n))
//   otherwise   a continued fraction for Q = 1 - P, evaluated by
//               the modified Lentz method
// CHISQUARE_CDF would give the same answer for whole 2k, but it takes
// WHOLE-number df only (GfxProbe T11 on the emulator, 17 Oct 2026:
// CHISQUARE_CDF(2.5, 3) is "Invalid input"), so the Gamma here never
// relies on it -- and since v2.7b the chi-square with a fractional df
// comes here too: P(X <= x) = P(df/2, x/2).
st_gammaP(gpA, gpX)
BEGIN
  LOCAL gpS, gpD, gpAp, gpI, gpG, gpB;
  LOCAL gpC, gpH, gpAn, gpDel, gpGo;
  IF gpX <= 0 THEN RETURN 0; END;
  gpG := st_lngamma(gpA);
  IF gpX < gpA + 1 THEN
    // the series: each term is the last one times x/(a+n)
    gpAp := gpA; gpS := 1/gpA; gpD := gpS;
    gpI := 0; gpGo := 1;
    WHILE gpGo == 1 AND gpI < 500 DO
      gpI := gpI + 1;
      gpAp := gpAp + 1;
      gpD := gpD*gpX/gpAp;
      gpS := gpS + gpD;
      IF ABS(gpD) < ABS(gpS)*0.000000000001 THEN gpGo := 0; END;
    END;
    RETURN MIN(1, gpS*EXP(-gpX + gpA*LN(gpX) - gpG));
  END;
  // the continued fraction for Q(a, x), modified Lentz
  gpB := gpX + 1 - gpA; gpC := 1E300; gpD := 1/gpB; gpH := gpD;
  gpI := 0; gpGo := 1;
  WHILE gpGo == 1 AND gpI < 500 DO
    gpI := gpI + 1;
    gpAn := -gpI*(gpI - gpA);
    gpB := gpB + 2;
    gpD := gpAn*gpD + gpB;
    IF ABS(gpD) < 1E-300 THEN gpD := 1E-300; END;
    gpC := gpB + gpAn/gpC;
    IF ABS(gpC) < 1E-300 THEN gpC := 1E-300; END;
    gpD := 1/gpD;
    gpDel := gpD*gpC;
    gpH := gpH*gpDel;
    IF ABS(gpDel - 1) < 0.00000000001 THEN gpGo := 0; END;
  END;
  RETURN MAX(0, 1 - EXP(-gpX + gpA*LN(gpX) - gpG)*gpH);
END;

// ---- regularized UPPER incomplete gamma Q(a, x) = 1 - P(a, x) (v2.7e) ----
// The right tail of the Gamma family: Gamma(k, theta) has P(X > x) =
// Q(k, x/theta), and the chi-square (Gamma(df/2, scale 2)) and the
// Erlang use it too.  WHY A SECOND ROUTINE: for x >= a + 1, st_gammaP
// works out Q by the continued fraction and then returns 1 - Q.  On the
// Prime's 12-digit reals, 1 - Q is exactly 1 once Q is below about
// 5E-13, and 1 - (1 - Q) is then 0: a far right tail would read
// "= 0.000000".  Here Q itself is returned, every digit kept.  For
// x < a + 1 the series gives P and 1 - P is safe: Q there is at least
// Q(a, a + 1), which is 0.14 at a = 1 and 0.48 at a = 1000 (for a
// shape below 1 it is smaller, 0.0002 at a = 0.001, and 1 - P still
// keeps 8 of its 12 digits).  The continued fraction below is
// st_gammaP's, line for line, with the gp names as gq.  At most 500
// passes: about 490 are needed at a = 5000, so the forms refuse a
// Gamma k or a chi-square df above 5000 (B1, fleet review 19 Oct 2026).
st_gammaQ(gqA, gqX)
BEGIN
  LOCAL gqB, gqC, gqD, gqH, gqI, gqAn;
  LOCAL gqDel, gqGo, gqG;
  IF gqX <= 0 THEN RETURN 1; END;
  // below a + 1: the series, through st_gammaP (Q is not small here)
  IF gqX < gqA + 1 THEN RETURN 1 - st_gammaP(gqA, gqX); END;
  gqG := st_lngamma(gqA);
  // the continued fraction for Q(a, x), modified Lentz (as st_gammaP)
  gqB := gqX + 1 - gqA; gqC := 1E300; gqD := 1/gqB; gqH := gqD;
  gqI := 0; gqGo := 1;
  WHILE gqGo == 1 AND gqI < 500 DO
    gqI := gqI + 1;
    gqAn := -gqI*(gqI - gqA);
    gqB := gqB + 2;
    gqD := gqAn*gqD + gqB;
    IF ABS(gqD) < 1E-300 THEN gqD := 1E-300; END;
    gqC := gqB + gqAn/gqC;
    IF ABS(gqC) < 1E-300 THEN gqC := 1E-300; END;
    gqD := 1/gqD;
    gqDel := gqD*gqC;
    gqH := gqH*gqDel;
    IF ABS(gqDel - 1) < 0.00000000001 THEN gqGo := 0; END;
  END;
  // Q = e^(-x) x^a / Gamma(a) times the fraction: NO subtraction
  RETURN MIN(1, EXP(-gqX + gqA*LN(gqX) - gqG)*gqH);
END;

// ---- the continued fraction inside the incomplete beta ----
// Numerical Recipes 6.4 (betacf), by the modified Lentz method: each
// pass takes one even and one odd step of the fraction.
st_betacf(bcA, bcB, bcX)
BEGIN
  LOCAL bcQab, bcQap, bcQam, bcC, bcD, bcH;
  LOCAL bcM, bcM2, bcAa, bcDel, bcGo;
  bcQab := bcA + bcB; bcQap := bcA + 1; bcQam := bcA - 1;
  bcC := 1; bcD := 1 - bcQab*bcX/bcQap;
  IF ABS(bcD) < 1E-300 THEN bcD := 1E-300; END;
  bcD := 1/bcD; bcH := bcD;
  bcM := 0; bcGo := 1;
  WHILE bcGo == 1 AND bcM < 300 DO
    bcM := bcM + 1; bcM2 := 2*bcM;
    // the even step
    bcAa := bcM*(bcB - bcM)*bcX/((bcQam + bcM2)*(bcA + bcM2));
    bcD := 1 + bcAa*bcD;
    IF ABS(bcD) < 1E-300 THEN bcD := 1E-300; END;
    bcC := 1 + bcAa/bcC;
    IF ABS(bcC) < 1E-300 THEN bcC := 1E-300; END;
    bcD := 1/bcD; bcH := bcH*bcD*bcC;
    // the odd step
    bcAa := -(bcA + bcM)*(bcQab + bcM)*bcX/((bcA + bcM2)*(bcQap + bcM2));
    bcD := 1 + bcAa*bcD;
    IF ABS(bcD) < 1E-300 THEN bcD := 1E-300; END;
    bcC := 1 + bcAa/bcC;
    IF ABS(bcC) < 1E-300 THEN bcC := 1E-300; END;
    bcD := 1/bcD; bcDel := bcD*bcC; bcH := bcH*bcDel;
    IF ABS(bcDel - 1) < 0.00000000001 THEN bcGo := 0; END;
  END;
  RETURN bcH;
END;

// ---- regularized incomplete beta I_x(a, b): the Beta(a, b) CDF ----
// I_x(a,b) = x^a (1-x)^b / (a B(a,b)) * betacf(a, b, x).  The
// fraction converges fast only for x < (a+1)/(a+b+2); above that,
// use the SYMMETRY  I_x(a, b) = 1 - I_(1-x)(b, a).
st_betaI(biA, biB, biX)
BEGIN
  LOCAL biL;
  IF biX <= 0 THEN RETURN 0; END;
  IF biX >= 1 THEN RETURN 1; END;
  // ln of x^a (1-x)^b / B(a,b), with B(a,b) = G(a)G(b)/G(a+b)
  biL := st_lngamma(biA + biB) - st_lngamma(biA) - st_lngamma(biB) + biA*LN(biX) + biB*LN(1 - biX);
  IF biX < (biA + 1)/(biA + biB + 2) THEN
    RETURN EXP(biL)*st_betacf(biA, biB, biX)/biA;
  END;
  RETURN 1 - EXP(biL)*st_betacf(biB, biA, 1 - biX)/biB;
END;

// ---- hypergeometric probability by COMB ----
// Draw hyM items without replacement from hyN, of which hyK are
// "successes": P(X = hyJ) = C(K, j) C(N-K, M-j) / C(N, M).
// N is capped at 1000 in the form, so C(N, M) stays far inside the
// Prime's range (C(1000, 500) is about 2.7E299).
st_hyp(hyN, hyK, hyM, hyJ)
BEGIN
  IF hyJ < 0 OR hyJ > hyM OR hyJ > hyK OR hyM - hyJ > hyN - hyK THEN RETURN 0; END;
  RETURN COMB(hyK, hyJ)*COMB(hyN - hyK, hyM - hyJ)/COMB(hyN, hyM);
END;

// Laplace CDF helper (piecewise).  v2.5: parameters x, b -> stX, stB.
st_lapCDF(stX, mu, stB)
BEGIN
  IF stX < mu THEN RETURN 0.5*EXP((stX-mu)/stB); END;
  RETURN 1 - 0.5*EXP(-(stX-mu)/stB);
END;

// ---- Gompertz mean and variance, by Simpson's rule ----
// For X >= 0, E[X] = integral of S(x) and E[X^2] = integral of
// 2x S(x), where S(x) = exp(-eta (e^(bx) - 1)) is the survival
// function.  S falls faster than exponentially, so 400 Simpson
// steps out to the point where S = e^-30 are plenty.  (There is no
// elementary closed form: the mean needs the exponential integral.)
st_gompMom(gmE, gmB)
BEGIN
  LOCAL gmX, gmH, gmI, gmS, gmW, gmM1;
  LOCAL gmM2, gmXi;
  gmX := LN(1 + 30/gmE)/gmB;
  gmH := gmX/400;
  gmM1 := 0; gmM2 := 0;
  FOR gmI FROM 0 TO 400 DO
    gmXi := gmI*gmH;
    gmS := EXP(-gmE*(EXP(gmB*gmXi) - 1));
    // Simpson weights 1, 4, 2, 4, ..., 2, 4, 1
    gmW := 2;
    IF gmI - 2*IP(gmI/2) == 1 THEN gmW := 4; END;
    IF gmI == 0 OR gmI == 400 THEN gmW := 1; END;
    gmM1 := gmM1 + gmW*gmS;
    gmM2 := gmM2 + gmW*2*gmXi*gmS;
  END;
  gmM1 := gmM1*gmH/3; gmM2 := gmM2*gmH/3;
  RETURN {gmM1, MAX(0, gmM2 - gmM1^2)};
END;


// ============================================================
// 3b. THE CAS (v2.7)
// The fleet's CAS kit (BRIEF_CAS, copied exactly, xx -> st) and the
// helpers around it.  The fence:
//   * st_Cas is the ONE place this program calls the CAS;
//   * CAS text holds whole numbers only, written by st_IStr: a decimal
//     goes in as a ratio of whole numbers (0.085 -> (85/10^3)), never
//     as STRING of a real (Home's number format can add digit marks
//     or a decimal comma the CAS cannot read);
//   * CAS variables carry the tag: stqx, stqt, stqk, stqa, stqb, stqc
//     (a student may own x, t, k, a, b or c in the CAS view);
//   * the CAS is asked only when the Formula page opens, never in a
//     frame or a drag, and the page it built is kept;
//   * every CAS answer is checked against the program's own number
//     (st_CasNum) before it is shown -- never a wrong formula;
//   * the numbers the plot draws come from the program, never the CAS.
// ============================================================

// A whole number as plain digits: no digit-grouping marks, no decimal
// point, no exponent.  Home's STRING() follows the Home number format
// (PyProbe saw grouping marks and a special minus sign), which the
// CAS cannot read, so CAS text is built only from this routine.
st_IStr(stN)
BEGIN
  LOCAL dS, dV, dD;
  dV := ROUND(ABS(stN), 0);
  IF dV == 0 THEN RETURN "0"; END;
  dS := "";
  WHILE dV > 0 DO
    dD := dV MOD 10;
    dS := CHAR(48 + dD) + dS;
    dV := (dV - dD) / 10;
  END;
  IF stN < 0 THEN dS := "-" + dS; END;
  RETURN dS;
END;

// The ONE place this program calls the CAS.  It returns the CAS answer as text,
// or "" if the CAS failed, gave undef, an error or an infinity.
// The call form is the one CASProbe beta 1b found on the HP Prime Virtual
// Calculator (build 2026-09-09, 17 Oct 2026): CAS(s) with s a variable hands
// the text back UNEVALUATED, so the text is built into a literal CAS("...")
// that EXPR reads at run time, and the CAS is asked for string(...) of the
// answer, so a big integer crosses into PPL whole, not as a 12-digit real.
st_Cas(stS)
BEGIN
  LOCAL cR;
  cR := "";
  IFERR
    cR := STRING(EXPR("CAS(" + CHAR(34) + "string(" + stS + ")" + CHAR(34) + ")"));
  THEN
    cR := "";
  END;
  // The answer is a PPL string, and STRING() of it adds one pair of quotes
  // (CASProbe beta 1b T12, emulator): drop that pair.
  IF DIM(cR) >= 2 THEN
    IF LEFT(cR, 1) == CHAR(34) AND MID(cR, DIM(cR), 1) == CHAR(34) THEN
      cR := MID(cR, 2, DIM(cR) - 2);
    END;
  END;
  // undef, an error, or an infinity (CHAR(8734) is the infinity sign; it
  // also catches the Prime's plus-minus infinity, 1/0): no usable answer
  // an empty answer is already a failure; INSTRING("", ...) is "Invalid input" on
  // the Prime (emulator, 17 Oct 2026), so the checks run only on a non-empty answer
  IF DIM(cR) > 0 THEN
    IF INSTRING(cR, "undef") > 0 OR INSTRING(cR, "rror") > 0 OR INSTRING(cR, CHAR(8734)) > 0 THEN cR := ""; END;
  END;
  RETURN cR;
END;

// A whole number from CAS text (the CAS's round(...) answer), read
// digit by digit, so no Home number format, decimal separator or EXPR
// is involved.  Returns {value, 1}, or {0, 0} if there is no number.
st_IVal(ivS)
BEGIN
  LOCAL ivI, ivC, ivD, ivV, ivSg, ivN;
  LOCAL ivGo;
  ivV := 0; ivSg := 1; ivN := 0; ivI := 1; ivGo := 1;
  WHILE ivGo AND ivI <= DIM(ivS) DO
    ivC := MID(ivS, ivI, 1);
    ivD := INSTRING("0123456789", ivC);     // 1 for "0" ... 10 for "9"
    IF ivD > 0 THEN
      ivV := ivV*10 + ivD - 1;
      ivN := ivN + 1;
    ELSE
      IF ivN == 0 AND (ivC == "-" OR ivC == CHAR(8722)) THEN
        ivSg := -1;
      ELSE
        IF ivN > 0 OR ivC <> " " THEN ivGo := 0; END;
      END;
    END;
    ivI := ivI + 1;
  END;
  IF ivN == 0 THEN RETURN {0, 0}; END;
  RETURN {ivSg*ivV, 1};
END;

// Every rpO in the text rpS replaced by rpN.  LEFT is never asked for
// 0 characters (LEFT(s, 0) is the WHOLE string on the Prime).  v2.7c:
// INSTRING("", ...) is "Invalid input" on the Prime (emulator, 17 Oct
// 2026), so the position is looked for only while text is left: a CAS
// answer ENDING in a replaced name, such as Pareto(1, 1)'s CDF
// (stqx-1)/stqx, emptied rpT and stopped the Formula page.  An empty
// pattern returns the text unchanged (INSTRING(s, "") is 1, which
// would loop for ever).
st_Repl(rpS, rpO, rpN)
BEGIN
  LOCAL rpR, rpT, rpP, rpLn;
  IF DIM(rpO) == 0 THEN RETURN rpS; END;
  rpR := ""; rpT := rpS; rpLn := DIM(rpO);
  rpP := 0; IF DIM(rpT) > 0 THEN rpP := INSTRING(rpT, rpO); END;
  WHILE rpP > 0 DO
    IF rpP > 1 THEN rpR := rpR + LEFT(rpT, rpP - 1); END;
    rpR := rpR + rpN;
    IF rpP + rpLn > DIM(rpT) THEN
      rpT := "";
    ELSE
      rpT := MID(rpT, rpP + rpLn, DIM(rpT) - rpP - rpLn + 1);
    END;
    rpP := 0; IF DIM(rpT) > 0 THEN rpP := INSTRING(rpT, rpO); END;
  END;
  RETURN rpR + rpT;
END;

// Long text cut into screen lines of at most wrW characters, broken
// just before a + or - (or at a space) near the end of a line when
// there is one, so a formula breaks between its terms.
st_Wrap(wrS, wrW)
BEGIN
  LOCAL wrL, wrR, wrB, wrI, wrC;
  wrL := {}; wrR := wrS;
  WHILE DIM(wrR) > wrW DO
    wrB := 0;
    wrI := wrW;
    WHILE wrB == 0 AND wrI > wrW - 16 DO
      wrC := MID(wrR, wrI, 1);
      IF wrC == "+" OR wrC == "-" OR wrC == " " THEN wrB := wrI - 1; END;
      wrI := wrI - 1;
    END;
    IF wrB < 1 THEN wrB := wrW; END;
    wrL := CONCAT(wrL, {LEFT(wrR, wrB)});
    wrR := MID(wrR, wrB + 1, DIM(wrR) - wrB);
  END;
  RETURN CONCAT(wrL, {wrR});
END;

// The power of ten that makes scV a whole number of about ten digits.
st_Scale(scV)
BEGIN
  IF ABS(scV) < 1E-300 THEN RETURN 9; END;
  RETURN 9 - FLOOR(LOG(ABS(scV)));
END;

// The CAS's own decimal for the CAS text cnT, asked as the whole
// number round(value x 10^cnE).  Returns {value, 1}, or {0, 0}.
st_CasNum(cnT, cnE)
BEGIN
  LOCAL cnR, cnV;
  cnR := st_Cas("round(evalf(" + cnT + ")*10^(" + st_IStr(cnE) + "))");
  cnV := st_IVal(cnR);
  IF cnV(2) == 0 THEN RETURN {0, 0}; END;
  RETURN {cnV(1)/10^cnE, 1};
END;

// Is raV a SIMPLE ratio -- a whole number over 10^d, d at most 6?
// Typed values nearly always are (2.5, 0.085, 1.96).  Returns
// {numerator, d}, or {} when it is not (1/3 typed as 0.333333333333).
st_Rat(raV)
BEGIN
  LOCAL raD, raM, raR, raGo;
  raD := 0; raGo := 1; raR := 0;
  WHILE raGo == 1 AND raD <= 6 DO
    raM := raV*10^raD;
    raR := ROUND(raM, 0);
    IF ABS(raM - raR) <= 0.000000001*MAX(1, ABS(raM)) AND ABS(raR) < 1E11 THEN
      raGo := 0;
    ELSE
      raD := raD + 1;
    END;
  END;
  IF raGo == 1 THEN RETURN {}; END;
  RETURN {raR, raD};
END;

// The CAS text of a simple ratio from st_Rat: (85/10^3), or (3).
st_RTxt(rtR)
BEGIN
  IF rtR(2) == 0 THEN RETURN "(" + st_IStr(rtR(1)) + ")"; END;
  RETURN "(" + st_IStr(rtR(1)) + "/10^" + st_IStr(rtR(2)) + ")";
END;

// ANY real as a ratio of whole numbers, to all 12 of its digits:
// 0.333333333333 -> (333333333333/10^(12)).  Used only to CHECK a
// symbolic answer at the user's value, never for display.
st_LTxt(ltV)
BEGIN
  LOCAL ltD;
  IF ABS(ltV) < 1E-300 THEN RETURN "(0)"; END;
  ltD := 11 - FLOOR(LOG(ABS(ltV)));
  RETURN "(" + st_IStr(ROUND(ltV*10^ltD, 0)) + "/10^(" + st_IStr(ltD) + "))";
END;

// How many digits a text holds (for "an exact fraction of 41 digits").
st_Digits(dgS)
BEGIN
  LOCAL dgI, dgN;
  dgN := 0;
  FOR dgI FROM 1 TO DIM(dgS) DO
    IF INSTRING("0123456789", MID(dgS, dgI, 1)) > 0 THEN dgN := dgN + 1; END;
  END;
  RETURN dgN;
END;


// ============================================================
// 4. THE DISPATCH
// Every routine here takes the distribution number dId (1..27) and,
// where it needs them, the parameters packed in one list dPl =
// {p1, p2, p3}.  Each CASE line is one distribution, so the
// formulas for a distribution can be read across st_pdf, st_cdf,
// st_icdf and st_Mom side by side.
// ============================================================

// The short name, used on the plot and on the numbers page.
st_Name(dId)
BEGIN
  LOCAL nmL;
  nmL := {"Normal", "Uniform", "Exponential", "Weibull", "LogNormal",
    "Laplace", "Student t", "Chi-square", "Fisher F", "Burr XII",
    "Gumbel (max)", "Pareto", "Rayleigh", "Logistic", "Binomial",
    "Poisson", "Geometric", "Triangular", "Cauchy", "Frechet",
    "Gompertz", "Maxwell-Boltzmann", "Erlang", "Gamma", "Beta",
    "Negative binomial", "Hypergeometric"};
  RETURN nmL(dId);
END;

// A typical use for criminal-justice and forensic students.
st_Use(dId)
BEGIN
  LOCAL usL;
  usL := {"test scores, errors", "time of a crime in a shift",
    "time to the next incident", "time to re-offence",
    "sentence lengths", "age-estimate errors", "small-sample means",
    "goodness of fit, tables", "ANOVA: ratio of variances",
    "time to re-arrest", "worst day of the year",
    "few offenders, most crime", "shot distance from aim",
    "risk-assessment scores", "convictions in n cases",
    "calls per shift", "searches to the first hit",
    "guess: min, likeliest, max", "ratio of errors (no mean!)",
    "largest losses, extremes", "years lived past age 20",
    "3-D miss distance", "time to the k-th call",
    "case processing time", "clearance rates",
    "stops before r-th arrest", "sampling seized packets"};
  RETURN usL(dId);
END;

// 1 for the discrete distributions (bars, k), 0 for the continuous.
st_Disc(dId)
BEGIN
  IF dId == 15 OR dId == 16 OR dId == 17 OR dId == 26 OR dId == 27 THEN RETURN 1; END;
  RETURN 0;
END;

// The starting values {p1, p2, p3, x, a, b} for each form -- the
// v2.5 form defaults (a and b are filled in by st_DefAB's rule).
st_Def(dId)
BEGIN
  CASE
    IF dId == 1 THEN RETURN {0, 1, 0, 0, -1, 1}; END;
    IF dId == 2 THEN RETURN {0, 1, 0, 0.5, 0, 1}; END;
    IF dId == 3 THEN RETURN {1, 0, 0, 0, 0, 1}; END;
    IF dId == 4 THEN RETURN {2, 1, 0, 1, 0, 1}; END;
    IF dId == 5 THEN RETURN {0, 1, 0, 1, 0.5, 2}; END;
    IF dId == 6 THEN RETURN {0, 1, 0, 0, -1, 1}; END;
    IF dId == 7 THEN RETURN {10, 0, 0, 0, -2, 2}; END;
    IF dId == 8 THEN RETURN {10, 0, 0, 5, 5, 20}; END;
    IF dId == 9 THEN RETURN {10, 15, 0, 1, 0.5, 2}; END;
    IF dId == 10 THEN RETURN {2, 2, 1, 1, 0.5, 2}; END;
    IF dId == 11 THEN RETURN {0, 1, 0, 0, -1, 2}; END;
    IF dId == 12 THEN RETURN {1, 3, 0, 2, 1, 2}; END;
    IF dId == 13 THEN RETURN {1.5, 0, 0, 2, 0, 2}; END;
    IF dId == 14 THEN RETURN {0, 1, 0, 0, -1, 1}; END;
    IF dId == 15 THEN RETURN {20, 0.5, 0, 10, 7, 12}; END;
    IF dId == 16 THEN RETURN {5, 0, 0, 3, 0, 9}; END;
    IF dId == 17 THEN RETURN {0.3, 0, 0, 3, 1, 8}; END;
    IF dId == 18 THEN RETURN {2, 5, 12, 5, 2, 5}; END;
    IF dId == 19 THEN RETURN {0, 1, 0, 1, -1, 1}; END;
    IF dId == 20 THEN RETURN {3, 1, 0, 1, 0, 2}; END;
    IF dId == 21 THEN RETURN {0.01, 0.085, 0, 50, 40, 60}; END;
    IF dId == 22 THEN RETURN {1, 0, 0, 1.5, 0, 2}; END;
    IF dId == 23 THEN RETURN {3, 1, 0, 2, 0, 3}; END;
    IF dId == 24 THEN RETURN {2.5, 2, 0, 4, 0, 5}; END;
    IF dId == 25 THEN RETURN {2, 5, 0, 0.3, 0.1, 0.4}; END;
    IF dId == 26 THEN RETURN {3, 0.4, 0, 4, 0, 4}; END;
    IF dId == 27 THEN RETURN {50, 40, 6, 5, 2, 4}; END;
  END;
  RETURN {0, 1, 0, 0, 0, 1};
END;

// The support {lowest, highest} value X can take; 1E300 stands for
// "no limit".  A discrete support is in whole numbers.
st_Sup(dId, dPl)
BEGIN
  LOCAL pr1, pr2, pr3;
  pr1 := dPl(1); pr2 := dPl(2); pr3 := dPl(3);
  CASE
    IF dId == 2 THEN RETURN {pr1, pr2}; END;
    IF dId == 3 OR dId == 4 OR dId == 5 OR dId == 8 OR dId == 9 OR dId == 10 THEN RETURN {0, 1E300}; END;
    IF dId == 13 OR dId == 21 OR dId == 22 OR dId == 23 OR dId == 24 THEN RETURN {0, 1E300}; END;
    IF dId == 12 THEN RETURN {pr1, 1E300}; END;
    IF dId == 15 THEN RETURN {0, pr1}; END;
    IF dId == 16 OR dId == 26 THEN RETURN {0, 1E300}; END;
    IF dId == 17 THEN RETURN {1, 1E300}; END;
    IF dId == 18 THEN RETURN {pr1, pr3}; END;
    IF dId == 20 THEN RETURN {pr3, 1E300}; END;
    IF dId == 25 THEN RETURN {0, 1}; END;
    IF dId == 27 THEN RETURN {MAX(0, pr3 - (pr1 - pr2)), MIN(pr3, pr2)}; END;
  END;
  RETURN {-1E300, 1E300};
END;

// ---- the density f(x), or for a discrete distribution P(X = x) ----
// Outside the support every density is 0; the guards come first so
// that no LN or power of a negative number is ever taken (a zoomed-out
// window reaches x < 0 for distributions that live on x >= 0).
st_pdf(dId, dPl, dX)
BEGIN
  LOCAL pr1, pr2, pr3, zz;
  pr1 := dPl(1); pr2 := dPl(2); pr3 := dPl(3);
  CASE
    // 1 Normal(mu, sigma): the built-in
    IF dId == 1 THEN RETURN NORMALD(pr1, pr2, dX); END;
    // 2 Uniform(a, b): flat at 1/(b-a) between a and b
    IF dId == 2 THEN
      IF dX < pr1 OR dX > pr2 THEN RETURN 0; END;
      RETURN 1/(pr2 - pr1);
    END;
    // 3 Exponential(rate lambda): lambda e^(-lambda x), x >= 0
    IF dId == 3 THEN
      IF dX < 0 THEN RETURN 0; END;
      RETURN pr1*EXP(-pr1*dX);
    END;
    // 4 Weibull(shape k, scale lambda): (k/l)(x/l)^(k-1) e^(-(x/l)^k)
    IF dId == 4 THEN
      IF dX <= 0 THEN RETURN 0; END;
      RETURN (pr1/pr2)*(dX/pr2)^(pr1 - 1)*EXP(-(dX/pr2)^pr1);
    END;
    // 5 LogNormal(mu, sigma of ln X): phi(z)/(x sigma), z = (ln x - mu)/sigma
    IF dId == 5 THEN
      IF dX <= 0 THEN RETURN 0; END;
      zz := (LN(dX) - pr1)/pr2;
      RETURN NORMALD(0, 1, zz)/(dX*pr2);
    END;
    // 6 Laplace(mu, b): e^(-|x - mu|/b) / (2b)
    IF dId == 6 THEN RETURN EXP(-ABS(dX - pr1)/pr2)/(2*pr2); END;
    // 7 Student t(df): the built-in
    IF dId == 7 THEN RETURN STUDENT(pr1, dX); END;
    // 8 Chi-square(df): the built-in density
    IF dId == 8 THEN
      IF dX <= 0 THEN RETURN 0; END;
      RETURN CHISQUARE(pr1, dX);
    END;
    // 9 Fisher F(df1, df2): the built-in for whole df.  v2.7b: with a
    //    fractional df, the formula, by logs (FISHER's density with a
    //    fractional df was not among the emulator's checks):
    //    f(x) = d1^(d1/2) d2^(d2/2) x^(d1/2 - 1)
    //           / ( B(d1/2, d2/2) (d1 x + d2)^((d1 + d2)/2) )
    IF dId == 9 THEN
      IF dX <= 0 THEN RETURN 0; END;
      IF FP(pr1) == 0 AND FP(pr2) == 0 THEN RETURN FISHER(pr1, pr2, dX); END;
      RETURN EXP((pr1/2)*LN(pr1) + (pr2/2)*LN(pr2) + (pr1/2 - 1)*LN(dX) - ((pr1 + pr2)/2)*LN(pr1*dX + pr2) + st_lngamma((pr1 + pr2)/2) - st_lngamma(pr1/2) - st_lngamma(pr2/2));
    END;
    // 10 Burr XII(C, K, S): K C (x/S)^C / (x (1 + (x/S)^C)^(K+1))
    IF dId == 10 THEN
      IF dX <= 0 THEN RETURN 0; END;
      RETURN (pr2*pr1*(dX/pr3)^pr1)/(dX*(1 + (dX/pr3)^pr1)^(pr2 + 1));
    END;
    // 11 Gumbel max(mu, beta): (1/beta) e^(-(z + e^-z)), z = (x - mu)/beta
    IF dId == 11 THEN
      zz := (dX - pr1)/pr2;
      IF zz < -600 THEN RETURN 0; END;   // e^-z would overflow; f is 0
      RETURN EXP(-(zz + EXP(-zz)))/pr2;
    END;
    // 12 Pareto(x_m, alpha): alpha x_m^alpha / x^(alpha+1), x >= x_m
    IF dId == 12 THEN
      IF dX < pr1 THEN RETURN 0; END;
      RETURN pr2*pr1^pr2/dX^(pr2 + 1);
    END;
    // 13 Rayleigh(sigma): (x/sigma^2) e^(-x^2/(2 sigma^2))
    IF dId == 13 THEN
      IF dX < 0 THEN RETURN 0; END;
      RETURN (dX/pr1^2)*EXP(-dX^2/(2*pr1^2));
    END;
    // 14 Logistic(mu, s): e^-z / (s (1 + e^-z)^2), z = (x - mu)/s
    IF dId == 14 THEN
      zz := (dX - pr1)/pr2;
      IF ABS(zz) > 600 THEN RETURN 0; END;
      RETURN EXP(-zz)/(pr2*(1 + EXP(-zz))^2);
    END;
    // 15 Binomial(n, p): the built-in, as v2.5 used it
    IF dId == 15 THEN
      IF dX < 0 OR dX > pr1 THEN RETURN 0; END;
      RETURN BINOMIAL(pr1, pr2, dX);
    END;
    // 16 Poisson(lambda): the built-in
    IF dId == 16 THEN
      IF dX < 0 THEN RETURN 0; END;
      RETURN POISSON(pr1, dX);
    END;
    // 17 Geometric(p), k = the trial of the first success: p(1-p)^(k-1)
    IF dId == 17 THEN
      IF dX < 1 THEN RETURN 0; END;
      RETURN pr1*(1 - pr1)^(dX - 1);
    END;
    // 18 Triangular(min a, mode c, max b): two straight sides
    IF dId == 18 THEN
      IF dX <= pr1 OR dX >= pr3 THEN RETURN 0; END;
      IF dX < pr2 THEN RETURN 2*(dX - pr1)/((pr3 - pr1)*(pr2 - pr1)); END;
      IF dX == pr2 THEN RETURN 2/(pr3 - pr1); END;
      RETURN 2*(pr3 - dX)/((pr3 - pr1)*(pr3 - pr2));
    END;
    // 19 Cauchy(x0, scale g): 1 / (pi g (1 + z^2)), z = (x - x0)/g
    IF dId == 19 THEN
      zz := (dX - pr1)/pr2;
      RETURN 1/(3.14159265358979*pr2*(1 + zz^2));
    END;
    // 20 Frechet(alpha, s, m): (alpha/s) z^(-1-alpha) e^(-z^-alpha), z = (x - m)/s
    IF dId == 20 THEN
      IF dX <= pr3 THEN RETURN 0; END;
      zz := (dX - pr3)/pr2;
      RETURN (pr1/pr2)*zz^(-1 - pr1)*EXP(-zz^(-pr1));
    END;
    // 21 Gompertz(eta, b): b eta e^(eta + b x - eta e^(b x)), x >= 0
    IF dId == 21 THEN
      IF dX < 0 OR pr2*dX > 700 THEN RETURN 0; END;
      RETURN pr2*pr1*EXP(pr1 + pr2*dX - pr1*EXP(pr2*dX));
    END;
    // 22 Maxwell-Boltzmann(a): sqrt(2/pi) z^2 e^(-z^2/2) / a, z = x/a
    IF dId == 22 THEN
      IF dX <= 0 THEN RETURN 0; END;
      zz := dX/pr1;
      RETURN 0.797884560802865*zz^2*EXP(-zz^2/2)/pr1;
    END;
    // 23 Erlang(k, lambda): 2 lambda * chi-square(2k) density at 2 lambda x
    //    (the sum of k exponential waits IS a scaled chi-square, 2k df)
    IF dId == 23 THEN
      IF dX <= 0 THEN RETURN 0; END;
      RETURN 2*pr2*CHISQUARE(2*pr1, 2*pr2*dX);
    END;
    // 24 Gamma(k, theta): x^(k-1) e^(-x/theta) / (Gamma(k) theta^k), by logs
    IF dId == 24 THEN
      IF dX <= 0 THEN RETURN 0; END;
      RETURN EXP((pr1 - 1)*LN(dX) - dX/pr2 - st_lngamma(pr1) - pr1*LN(pr2));
    END;
    // 25 Beta(alpha, beta): x^(a-1) (1-x)^(b-1) / B(a, b), by logs
    IF dId == 25 THEN
      IF dX <= 0 OR dX >= 1 THEN RETURN 0; END;
      RETURN EXP((pr1 - 1)*LN(dX) + (pr2 - 1)*LN(1 - dX) + st_lngamma(pr1 + pr2) - st_lngamma(pr1) - st_lngamma(pr2));
    END;
    // 26 Negative binomial(r, p), k = failures before the r-th success:
    //    C(k+r-1, k) p^r (1-p)^k = (r/(k+r)) * BINOMIAL(k+r, p, r)
    IF dId == 26 THEN
      IF dX < 0 THEN RETURN 0; END;
      RETURN pr1/(dX + pr1)*BINOMIAL(dX + pr1, pr2, pr1);
    END;
    // 27 Hypergeometric(N, K, n): COMB, in st_hyp
    IF dId == 27 THEN RETURN st_hyp(pr1, pr2, pr3, dX); END;
  END;
  RETURN 0;
END;

// ---- the CDF F(x) = P(X <= x) ----
st_cdf(dId, dPl, dX)
BEGIN
  LOCAL pr1, pr2, pr3, zz, hk;
  pr1 := dPl(1); pr2 := dPl(2); pr3 := dPl(3);
  CASE
    IF dId == 1 THEN RETURN NORMALD_CDF(pr1, pr2, dX); END;
    IF dId == 2 THEN RETURN MIN(1, MAX(0, (dX - pr1)/(pr2 - pr1))); END;
    IF dId == 3 THEN
      IF dX <= 0 THEN RETURN 0; END;
      RETURN 1 - EXP(-pr1*dX);
    END;
    IF dId == 4 THEN
      IF dX <= 0 THEN RETURN 0; END;
      RETURN 1 - EXP(-(dX/pr2)^pr1);
    END;
    IF dId == 5 THEN
      IF dX <= 0 THEN RETURN 0; END;
      RETURN NORMALD_CDF(0, 1, (LN(dX) - pr1)/pr2);
    END;
    IF dId == 6 THEN RETURN st_lapCDF(dX, pr1, pr2); END;
    IF dId == 7 THEN RETURN STUDENT_CDF(pr1, dX); END;
    // 8 Chi-square(df).  v2.7b: CHISQUARE_CDF takes a WHOLE df only
    //    (GfxProbe T11: CHISQUARE_CDF(2.5, 3) is "Invalid input"), so a
    //    fractional df uses the incomplete gamma: a chi-square with df
    //    degrees of freedom is a Gamma(df/2, scale 2), F(x) = P(df/2, x/2)
    IF dId == 8 THEN
      IF dX <= 0 THEN RETURN 0; END;
      IF FP(pr1) == 0 THEN RETURN CHISQUARE_CDF(pr1, dX); END;
      RETURN st_gammaP(pr1/2, dX/2);
    END;
    // 9 Fisher F(df1, df2).  v2.7b: FISHER_CDF takes WHOLE df only (T11),
    //    so a fractional df uses the incomplete beta: if X is F(d1, d2),
    //    then d1 X / (d1 X + d2) is Beta(d1/2, d2/2), so
    //    F(x) = I_z(d1/2, d2/2) with z = d1 x / (d1 x + d2)
    IF dId == 9 THEN
      IF dX <= 0 THEN RETURN 0; END;
      IF FP(pr1) == 0 AND FP(pr2) == 0 THEN RETURN FISHER_CDF(pr1, pr2, dX); END;
      RETURN st_betaI(pr1/2, pr2/2, pr1*dX/(pr1*dX + pr2));
    END;
    // 10 Burr XII: 1 - (1 + (x/S)^C)^(-K)
    IF dId == 10 THEN
      IF dX <= 0 THEN RETURN 0; END;
      RETURN 1 - (1 + (dX/pr3)^pr1)^(-pr2);
    END;
    // 11 Gumbel: e^(-e^(-z))
    IF dId == 11 THEN
      zz := (dX - pr1)/pr2;
      IF zz < -600 THEN RETURN 0; END;
      RETURN EXP(-EXP(-zz));
    END;
    // 12 Pareto: 1 - (x_m/x)^alpha
    IF dId == 12 THEN
      IF dX <= pr1 THEN RETURN 0; END;
      RETURN 1 - (pr1/dX)^pr2;
    END;
    // 13 Rayleigh: 1 - e^(-x^2/(2 sigma^2))
    IF dId == 13 THEN
      IF dX <= 0 THEN RETURN 0; END;
      RETURN 1 - EXP(-dX^2/(2*pr1^2));
    END;
    // 14 Logistic: 1/(1 + e^-z)
    IF dId == 14 THEN
      zz := (dX - pr1)/pr2;
      IF zz < -600 THEN RETURN 0; END;
      IF zz > 600 THEN RETURN 1; END;
      RETURN 1/(1 + EXP(-zz));
    END;
    // 15 Binomial: the built-in BINOMIAL_CDF(n, p, k), v2.5's order
    IF dId == 15 THEN
      IF dX < 0 THEN RETURN 0; END;
      IF dX >= pr1 THEN RETURN 1; END;
      RETURN BINOMIAL_CDF(pr1, pr2, dX);
    END;
    // 16 Poisson: 0 below k = 0 -- never POISSON_CDF(lambda, -1)
    //    (v2.5's fix; what the built-in does there is unverified)
    IF dId == 16 THEN
      IF dX < 0 THEN RETURN 0; END;
      RETURN POISSON_CDF(pr1, dX);
    END;
    // 17 Geometric: 1 - (1-p)^k
    IF dId == 17 THEN
      IF dX < 1 THEN RETURN 0; END;
      RETURN 1 - (1 - pr1)^dX;
    END;
    // 18 Triangular: a parabola up to the mode, then one down
    IF dId == 18 THEN
      IF dX <= pr1 THEN RETURN 0; END;
      IF dX >= pr3 THEN RETURN 1; END;
      IF dX <= pr2 THEN RETURN (dX - pr1)^2/((pr3 - pr1)*(pr2 - pr1)); END;
      RETURN 1 - (pr3 - dX)^2/((pr3 - pr1)*(pr3 - pr2));
    END;
    // 19 Cauchy: 1/2 + arctan(z)/pi (arctan in radians, any angle mode)
    IF dId == 19 THEN RETURN 0.5 + st_atan((dX - pr1)/pr2)/3.14159265358979; END;
    // 20 Frechet: e^(-z^-alpha)
    IF dId == 20 THEN
      IF dX <= pr3 THEN RETURN 0; END;
      RETURN EXP(-((dX - pr3)/pr2)^(-pr1));
    END;
    // 21 Gompertz: 1 - e^(-eta (e^(b x) - 1))
    IF dId == 21 THEN
      IF dX <= 0 THEN RETURN 0; END;
      IF pr2*dX > 700 THEN RETURN 1; END;
      RETURN 1 - EXP(-pr1*(EXP(pr2*dX) - 1));
    END;
    // 22 Maxwell-Boltzmann: erf(x/(a sqrt 2)) - sqrt(2/pi)(x/a) e^(-x^2/(2a^2)),
    //    with erf(z/sqrt 2) = 2 Phi(z) - 1
    IF dId == 22 THEN
      IF dX <= 0 THEN RETURN 0; END;
      zz := dX/pr1;
      RETURN 2*NORMALD_CDF(0, 1, zz) - 1 - 0.797884560802865*zz*EXP(-zz^2/2);
    END;
    // 23 Erlang(k, lambda): CHISQUARE_CDF(2k, 2 lambda x) -- whole df only;
    //    k is a whole number (st_ask_erlang rounds it down), so 2k is too
    IF dId == 23 THEN
      IF dX <= 0 THEN RETURN 0; END;
      RETURN CHISQUARE_CDF(2*pr1, 2*pr2*dX);
    END;
    // 24 Gamma(k, theta): our own P(k, x/theta)
    IF dId == 24 THEN RETURN st_gammaP(pr1, dX/pr2); END;
    // 25 Beta(alpha, beta): our own I_x(alpha, beta)
    IF dId == 25 THEN RETURN st_betaI(pr1, pr2, dX); END;
    // 26 Negative binomial: P(Y <= k) = P(at least r successes in k+r
    //    trials) = P(at most k FAILURES in k+r trials), and the failures
    //    are Binomial(k+r, 1-p): BINOMIAL_CDF(k+r, 1-p, k).  v2.7e: it
    //    was 1 - BINOMIAL_CDF(k+r, p, r-1), the same number, but on the
    //    Prime's 12 digits a tiny lower tail (NegBin(30, 0.3), P(Y <= 0)
    //    = 2E-16) came out exactly 0 (fleet review 19 Oct 2026, B4)
    IF dId == 26 THEN
      IF dX < 0 THEN RETURN 0; END;
      RETURN BINOMIAL_CDF(dX + pr1, 1 - pr2, dX);
    END;
    // 27 Hypergeometric: add up P(X = j) from the smallest possible j
    IF dId == 27 THEN
      zz := 0;
      FOR hk FROM MAX(0, pr3 - (pr1 - pr2)) TO MIN(dX, MIN(pr3, pr2)) DO
        zz := zz + st_hyp(pr1, pr2, pr3, hk);
      END;
      RETURN MIN(1, zz);
    END;
  END;
  RETURN 0;
END;

// ---- the RIGHT tail S(x) = P(X > x), family by family (v2.7e) ----
// Why not simply 1 - F(x)?  The Prime keeps 12 significant digits.  Far
// out on the right F(x) is 0.999999999999..., and 1 - F(x) keeps only
// the few digits left after the 9s: 1 - NORMALD_CDF(0, 1, 7.5) is
// exactly 0 on a 12-digit machine, though the tail is 3.19E-14.  So a
// far right tail read "= 0.000000" where the same tail on the left read
// "< 0.000001", and the Formula page threw away exact answers that were
// right (fleet review 19 Oct 2026, B3, B4; Roger's decision D9).  Here
// each family's right tail is a LOWER tail of something, or a formula
// with no subtraction from 1 (the 27 formulas of findings_B part I,
// each checked against SciPy):
//   * a symmetric family is mirrored about its centre: P(X > x) =
//     P(X < 2 mu - x) (Normal, Laplace, t, LogNormal through ln x);
//   * a CDF of the form 1 - g(x) gives g(x) itself (Exponential,
//     Weibull, Burr, Pareto, Rayleigh, Gompertz, Geometric ...);
//   * a count is turned round: X > k when the n - X failures number
//     fewer than n - k (Binomial); "fewer than r successes in the
//     first k + r trials" (Negative binomial);
//   * the incomplete gamma and beta give the upper part directly
//     (chi-square, Gamma, Erlang, F, Beta).
// Left of the support the tail is 1, right of it 0.
st_sf(dId, dPl, dX)
BEGIN
  LOCAL pr1, pr2, pr3, zz, hk, sfU;
  LOCAL sfLo, sfHi;
  pr1 := dPl(1); pr2 := dPl(2); pr3 := dPl(3);
  CASE
    // 1 Normal: P(X > x) = P(X < 2 mu - x), the mirror image
    IF dId == 1 THEN RETURN NORMALD_CDF(pr1, pr2, 2*pr1 - dX); END;
    // 2 Uniform: the share of [a, b] above x
    IF dId == 2 THEN RETURN MIN(1, MAX(0, (pr2 - dX)/(pr2 - pr1))); END;
    // 3 Exponential: e^(-lambda x)
    IF dId == 3 THEN
      IF dX <= 0 THEN RETURN 1; END;
      RETURN EXP(-pr1*dX);
    END;
    // 4 Weibull: e^(-(x/lambda)^k)
    IF dId == 4 THEN
      IF dX <= 0 THEN RETURN 1; END;
      RETURN EXP(-(dX/pr2)^pr1);
    END;
    // 5 LogNormal: ln X is Normal(mu, sigma): the mirrored normal tail
    IF dId == 5 THEN
      IF dX <= 0 THEN RETURN 1; END;
      RETURN NORMALD_CDF(0, 1, -(LN(dX) - pr1)/pr2);
    END;
    // 6 Laplace: mirrored about mu
    IF dId == 6 THEN RETURN st_lapCDF(2*pr1 - dX, pr1, pr2); END;
    // 7 Student t: symmetric about 0, P(T > t) = P(T < -t)
    IF dId == 7 THEN RETURN STUDENT_CDF(pr1, -dX); END;
    // 8 Chi-square(df) is Gamma(df/2, scale 2): Q(df/2, x/2), our own
    //    (the form keeps df at most 5000, within st_gammaQ's 500 passes)
    IF dId == 8 THEN
      IF dX <= 0 THEN RETURN 1; END;
      RETURN st_gammaQ(pr1/2, dX/2);
    END;
    // 9 Fisher F(d1, d2): z = d2/(d1 x + d2) is Beta(d2/2, d1/2), and z
    //    is SMALL when x is large, so P(X > x) = I_z(d2/2, d1/2), a lower
    //    tail.  With a df above 10000, v2.7d's 1 - F(x) is kept (for a
    //    whole df the built-in FISHER_CDF): st_betaI's log-gamma terms are
    //    each about a ln a, and on 12-digit reals their difference loses
    //    digits -- about 3E-6 of the tail at df 100000, visible in the
    //    sixth decimal (PC model of the Prime's rounding) -- and above df
    //    1E6 its 300 passes run out.  The far tail of so large an F is no
    //    classroom case.
    IF dId == 9 THEN
      IF dX <= 0 THEN RETURN 1; END;
      IF MAX(pr1, pr2) > 10000 THEN RETURN 1 - st_cdf(9, dPl, dX); END;
      RETURN st_betaI(pr2/2, pr1/2, pr2/(pr1*dX + pr2));
    END;
    // 10 Burr XII: (1 + (x/S)^C)^(-K)
    IF dId == 10 THEN
      IF dX <= 0 THEN RETURN 1; END;
      RETURN (1 + (dX/pr3)^pr1)^(-pr2);
    END;
    // 11 Gumbel: 1 - e^(-u) with u = e^(-z); for u below 1E-4 the first
    //    terms of its series, u - u^2/2 + u^3/6 (1 - e^(-u) would keep
    //    only 8 of the 12 digits there, and none at all below 5E-13)
    IF dId == 11 THEN
      zz := (dX - pr1)/pr2;
      IF zz < -600 THEN RETURN 1; END;
      sfU := EXP(-zz);
      IF sfU < 0.0001 THEN RETURN sfU - sfU^2/2 + sfU^3/6; END;
      RETURN 1 - EXP(-sfU);
    END;
    // 12 Pareto: (x_m/x)^alpha
    IF dId == 12 THEN
      IF dX <= pr1 THEN RETURN 1; END;
      RETURN (pr1/dX)^pr2;
    END;
    // 13 Rayleigh: e^(-x^2/(2 sigma^2))
    IF dId == 13 THEN
      IF dX <= 0 THEN RETURN 1; END;
      RETURN EXP(-dX^2/(2*pr1^2));
    END;
    // 14 Logistic: 1/(1 + e^z); beyond 600 scale units 0, as the CDF's
    //    left side is (st_cdf), so the two sides read alike
    IF dId == 14 THEN
      zz := (dX - pr1)/pr2;
      IF zz > 600 THEN RETURN 0; END;
      IF zz < -600 THEN RETURN 1; END;
      RETURN 1/(1 + EXP(zz));
    END;
    // 15 Binomial(n, p): X > k exactly when the n - X failures number
    //    at most n - k - 1, and n - X is Binomial(n, 1 - p)
    IF dId == 15 THEN
      IF dX < 0 THEN RETURN 1; END;
      IF dX >= pr1 THEN RETURN 0; END;
      RETURN BINOMIAL_CDF(pr1, 1 - pr2, pr1 - dX - 1);
    END;
    // 16 Poisson(lambda): X > k exactly when the (k + 1)-th event comes
    //    before time 1 -- the lower incomplete gamma P(k + 1, lambda), a
    //    lower tail: our own st_gammaP up to k + 1 = 5000; beyond, where
    //    its 500 passes are too few, the same number from the built-in
    //    CHISQUARE_CDF (2 x Gamma(k + 1) is a chi-square with the WHOLE
    //    df 2k + 2, the identity the Erlang below uses)
    IF dId == 16 THEN
      IF dX < 0 THEN RETURN 1; END;
      IF dX + 1 <= 5000 THEN RETURN st_gammaP(dX + 1, pr1); END;
      RETURN CHISQUARE_CDF(2*dX + 2, 2*pr1);
    END;
    // 17 Geometric: no success in the first k trials, (1 - p)^k
    IF dId == 17 THEN
      IF dX < 1 THEN RETURN 1; END;
      RETURN (1 - pr1)^dX;
    END;
    // 18 Triangular(min a, mode c, max b): above the mode the mirror
    //    quadratic (b - x)^2/((b - a)(b - c)).  Below it, 1 - (x - a)^2/
    //    ((b - a)(c - a)) would cancel when the mode sits at or near the
    //    maximum (the tail can then be tiny); the same number written as
    //    ((c - a)(b - c) + (c - x)(c + x - 2a)) / ((b - a)(c - a)) adds
    //    two parts that are never negative, so nothing cancels.
    IF dId == 18 THEN
      IF dX <= pr1 THEN RETURN 1; END;
      IF dX >= pr3 THEN RETURN 0; END;
      IF dX <= pr2 THEN RETURN ((pr2 - pr1)*(pr3 - pr2) + (pr2 - dX)*(pr2 + dX - 2*pr1))/((pr3 - pr1)*(pr2 - pr1)); END;
      RETURN (pr3 - dX)^2/((pr3 - pr1)*(pr3 - pr2));
    END;
    // 19 Cauchy: 1/2 - arctan(z)/pi; for z > 1 the same number as
    //    arctan(1/z)/pi (arctan(z) + arctan(1/z) is a right angle), which
    //    keeps its digits far out where 1/2 - arctan(z)/pi would not
    IF dId == 19 THEN
      zz := (dX - pr1)/pr2;
      IF zz > 1 THEN RETURN st_atan(1/zz)/3.14159265358979; END;
      RETURN 0.5 - st_atan(zz)/3.14159265358979;
    END;
    // 20 Frechet: 1 - e^(-u), u = z^(-alpha); the series for small u,
    //    as the Gumbel
    IF dId == 20 THEN
      IF dX <= pr3 THEN RETURN 1; END;
      sfU := ((dX - pr3)/pr2)^(-pr1);
      IF sfU < 0.0001 THEN RETURN sfU - sfU^2/2 + sfU^3/6; END;
      RETURN 1 - EXP(-sfU);
    END;
    // 21 Gompertz: e^(-eta (e^(b x) - 1))
    IF dId == 21 THEN
      IF dX <= 0 THEN RETURN 1; END;
      IF pr2*dX > 700 THEN RETURN 0; END;
      RETURN EXP(-pr1*(EXP(pr2*dX) - 1));
    END;
    // 22 Maxwell-Boltzmann: 2 Phi(-z) + sqrt(2/pi) z e^(-z^2/2), z = x/a:
    //    two positive parts, so nothing cancels
    IF dId == 22 THEN
      IF dX <= 0 THEN RETURN 1; END;
      zz := dX/pr1;
      RETURN 2*NORMALD_CDF(0, 1, -zz) + 0.797884560802865*zz*EXP(-zz^2/2);
    END;
    // 23 Erlang(k, lambda): Q(k, lambda x), our own st_gammaQ up to
    //    k = 5000; beyond, the same number from the Poisson count: the
    //    k-th call comes after time x exactly when FEWER than k calls
    //    arrive by x, POISSON_CDF(lambda x, k - 1) -- a lower tail
    IF dId == 23 THEN
      IF dX <= 0 THEN RETURN 1; END;
      IF pr1 <= 5000 THEN RETURN st_gammaQ(pr1, pr2*dX); END;
      RETURN POISSON_CDF(pr2*dX, pr1 - 1);
    END;
    // 24 Gamma(k, theta): Q(k, x/theta) (the form keeps k at most 5000)
    IF dId == 24 THEN
      IF dX <= 0 THEN RETURN 1; END;
      RETURN st_gammaQ(pr1, dX/pr2);
    END;
    // 25 Beta(alpha, beta): 1 - X is Beta(beta, alpha), so P(X > x) =
    //    P(1 - X < 1 - x) = I_(1-x)(beta, alpha), a lower tail
    IF dId == 25 THEN RETURN st_betaI(pr2, pr1, 1 - dX); END;
    // 26 Negative binomial(r, p): more than k failures before the r-th
    //    success exactly when the first k + r trials hold FEWER than r
    //    successes: BINOMIAL_CDF(k + r, p, r - 1)
    IF dId == 26 THEN
      IF dX < 0 THEN RETURN 1; END;
      RETURN BINOMIAL_CDF(dX + pr1, pr2, pr1 - 1);
    END;
    // 27 Hypergeometric: add P(X = j) from k + 1 to the largest possible
    //    j (the CDF adds from the smallest); the loop's start is never
    //    past its end (the two IFs), so no empty FOR is relied on
    IF dId == 27 THEN
      sfLo := MAX(0, pr3 - (pr1 - pr2)); sfHi := MIN(pr3, pr2);
      IF dX >= sfHi THEN RETURN 0; END;
      IF dX < sfLo THEN RETURN 1; END;
      zz := 0;
      FOR hk FROM dX + 1 TO sfHi DO
        zz := zz + st_hyp(pr1, pr2, pr3, hk);
      END;
      RETURN MIN(1, zz);
    END;
  END;
  RETURN 0;
END;

// ---- a quantile by bisection, for the CDFs with no inverse ----
// First widen [qLo, qHi] until it holds the answer (doubling the
// step), then halve it 34 times: 2^-34 of the bracket is far below
// what a pixel can show.
st_qbis(dId, dPl, qP, qLo, qHi)
BEGIN
  LOCAL qa, qb, qm, qw, qi;
  qa := qLo; qb := qHi; qw := qb - qa;
  qi := 0;
  WHILE st_cdf(dId, dPl, qb) < qP AND qi < 200 DO
    qa := qb; qb := qb + qw; qw := 2*qw; qi := qi + 1;
  END;
  qi := 0; qw := qb - qa;
  WHILE st_cdf(dId, dPl, qa) > qP AND qi < 200 DO
    qb := qa; qa := qa - qw; qw := 2*qw; qi := qi + 1;
  END;
  FOR qi FROM 1 TO 34 DO
    qm := (qa + qb)/2;
    IF st_cdf(dId, dPl, qm) < qP THEN qa := qm; ELSE qb := qm; END;
  END;
  RETURN (qa + qb)/2;
END;

// ---- a discrete quantile: the smallest whole k in [qLo, qHi] with
// P(X <= k) >= qP, by halving the range of whole numbers ----
st_dq(dId, dPl, qP, qLo, qHi)
BEGIN
  LOCAL qa, qb, qm;
  qa := qLo; qb := qHi;
  WHILE qb > qa DO
    qm := IP((qa + qb)/2);
    IF st_cdf(dId, dPl, qm) >= qP THEN qb := qm; ELSE qa := qm + 1; END;
  END;
  RETURN qa;
END;

// ---- the quantile: the x with F(x) = qP (used to fit the window) ----
// Closed forms where they exist; verified built-in inverses where the
// df is whole; bisection on the CDF otherwise.
st_icdf(dId, dPl, qP)
BEGIN
  LOCAL pr1, pr2, pr3, qF;
  pr1 := dPl(1); pr2 := dPl(2); pr3 := dPl(3);
  CASE
    IF dId == 1 THEN RETURN NORMALD_ICDF(pr1, pr2, qP); END;
    IF dId == 2 THEN RETURN pr1 + qP*(pr2 - pr1); END;
    IF dId == 3 THEN RETURN -LN(1 - qP)/pr1; END;
    IF dId == 4 THEN RETURN pr2*(-LN(1 - qP))^(1/pr1); END;
    IF dId == 5 THEN RETURN EXP(pr1 + pr2*NORMALD_ICDF(0, 1, qP)); END;
    IF dId == 6 THEN
      IF qP < 0.5 THEN RETURN pr1 + pr2*LN(2*qP); END;
      RETURN pr1 - pr2*LN(2*(1 - qP));
    END;
    IF dId == 7 THEN RETURN st_qbis(7, dPl, qP, -1, 1); END;
    // chi-square: CHISQUARE_ICDF takes a whole df only (GfxProbe T11);
    // otherwise bisection on st_cdf, which (v2.7b) then uses st_gammaP
    IF dId == 8 THEN
      IF FP(pr1) == 0 THEN RETURN CHISQUARE_ICDF(pr1, qP); END;
      RETURN st_qbis(8, dPl, qP, 0, pr1 + 1);
    END;
    IF dId == 9 THEN RETURN st_qbis(9, dPl, qP, 0, 2); END;
    IF dId == 10 THEN RETURN pr3*((1 - qP)^(-1/pr2) - 1)^(1/pr1); END;
    IF dId == 11 THEN RETURN pr1 - pr2*LN(-LN(qP)); END;
    IF dId == 12 THEN RETURN pr1*(1 - qP)^(-1/pr2); END;
    IF dId == 13 THEN RETURN pr1*(-2*LN(1 - qP))^0.5; END;
    IF dId == 14 THEN RETURN pr1 + pr2*LN(qP/(1 - qP)); END;
    IF dId == 15 THEN RETURN st_dq(15, dPl, qP, 0, pr1); END;
    IF dId == 16 THEN RETURN st_dq(16, dPl, qP, 0, IP(pr1 + 12*pr1^0.5 + 12)); END;
    IF dId == 17 THEN RETURN MAX(1, CEILING(LN(1 - qP)/LN(1 - pr1) - 0.000000001)); END;
    IF dId == 18 THEN
      qF := (pr2 - pr1)/(pr3 - pr1);
      IF qP < qF THEN RETURN pr1 + (qP*(pr3 - pr1)*(pr2 - pr1))^0.5; END;
      RETURN pr3 - ((1 - qP)*(pr3 - pr1)*(pr3 - pr2))^0.5;
    END;
    IF dId == 19 THEN RETURN pr1 + pr2*st_tanr(3.14159265358979*(qP - 0.5)); END;
    IF dId == 20 THEN RETURN pr3 + pr2*(-LN(qP))^(-1/pr1); END;
    IF dId == 21 THEN RETURN LN(1 - LN(1 - qP)/pr1)/pr2; END;
    // Maxwell: (X/a)^2 is chi-square with 3 df -- a whole df, verified
    IF dId == 22 THEN RETURN pr1*CHISQUARE_ICDF(3, qP)^0.5; END;
    // Erlang: 2 lambda X is chi-square with 2k df -- whole, verified
    IF dId == 23 THEN RETURN CHISQUARE_ICDF(2*pr1, qP)/(2*pr2); END;
    IF dId == 24 THEN RETURN st_qbis(24, dPl, qP, 0, pr1*pr2 + pr2); END;
    IF dId == 25 THEN RETURN st_qbis(25, dPl, qP, 0, 1); END;
    IF dId == 26 THEN RETURN st_dq(26, dPl, qP, 0, IP(pr1*(1 - pr2)/pr2 + 15*(pr1*(1 - pr2))^0.5/pr2 + 15)); END;
    IF dId == 27 THEN RETURN st_dq(27, dPl, qP, MAX(0, pr3 - (pr1 - pr2)), MIN(pr3, pr2)); END;
  END;
  RETURN 0;
END;

// ---- mean and variance {mean, variance} ----
// Either can be a STRING ("undef ...") when the moment does not exist --
// itself a lesson: a Cauchy has no mean at all, a Pareto with alpha <= 1
// none either.
st_Mom(dId, dPl)
BEGIN
  LOCAL pr1, pr2, pr3, mM, mV, mE2;
  pr1 := dPl(1); pr2 := dPl(2); pr3 := dPl(3);
  mM := 0; mV := 0;
  CASE
    IF dId == 1 THEN mM := pr1; mV := pr2^2; END;
    IF dId == 2 THEN mM := (pr1 + pr2)/2; mV := (pr2 - pr1)^2/12; END;
    IF dId == 3 THEN mM := 1/pr1; mV := 1/pr1^2; END;
    // Weibull: lambda Gamma(1 + 1/k), lambda^2 Gamma(1 + 2/k) - mean^2
    IF dId == 4 THEN
      mM := pr2*EXP(st_lngamma(1 + 1/pr1));
      mV := MAX(0, pr2^2*EXP(st_lngamma(1 + 2/pr1)) - mM^2);
    END;
    IF dId == 5 THEN mM := EXP(pr1 + pr2^2/2); mV := (EXP(pr2^2) - 1)*EXP(2*pr1 + pr2^2); END;
    IF dId == 6 THEN mM := pr1; mV := 2*pr2^2; END;
    IF dId == 7 THEN
      IF pr1 > 1 THEN mM := 0; ELSE mM := "undef (df<=1)"; END;
      IF pr1 > 2 THEN mV := pr1/(pr1 - 2); ELSE mV := "undef (df<=2)"; END;
    END;
    IF dId == 8 THEN mM := pr1; mV := 2*pr1; END;
    IF dId == 9 THEN
      IF pr2 > 2 THEN mM := pr2/(pr2 - 2); ELSE mM := "undef (df2<=2)"; END;
      IF pr2 > 4 THEN
        mV := 2*pr2^2*(pr1 + pr2 - 2)/(pr1*(pr2 - 2)^2*(pr2 - 4));
      ELSE
        mV := "undef (df2<=4)";
      END;
    END;
    // Burr XII: the v2.5 Lanczos moments; the mean needs K C > 1,
    // the variance K C > 2
    IF dId == 10 THEN
      IF pr2*pr1 > 1 THEN
        mM := pr3*EXP(st_lngamma(1 + 1/pr1) + st_lngamma(pr2 - 1/pr1) - st_lngamma(pr2));
      ELSE
        mM := "undef (K*C<=1)";
      END;
      IF pr2*pr1 > 2 THEN
        mE2 := pr3^2*EXP(st_lngamma(1 + 2/pr1) + st_lngamma(pr2 - 2/pr1) - st_lngamma(pr2));
        mV := MAX(0, mE2 - mM^2);
      ELSE
        mV := "undef (K*C<=2)";
      END;
    END;
    // Gumbel: mu + Euler's gamma beta, (pi^2/6) beta^2
    IF dId == 11 THEN mM := pr1 + 0.577215664901533*pr2; mV := 1.64493406684823*pr2^2; END;
    IF dId == 12 THEN
      IF pr2 > 1 THEN mM := pr2*pr1/(pr2 - 1); ELSE mM := "undef (α<=1)"; END;
      IF pr2 > 2 THEN mV := pr1^2*pr2/((pr2 - 1)^2*(pr2 - 2)); ELSE mV := "undef (α<=2)"; END;
    END;
    // Rayleigh: sigma sqrt(pi/2), (4 - pi)/2 sigma^2
    IF dId == 13 THEN mM := 1.2533141373155*pr1; mV := 0.429203673205103*pr1^2; END;
    // Logistic: mu, (pi^2/3) s^2
    IF dId == 14 THEN mM := pr1; mV := 3.28986813369645*pr2^2; END;
    IF dId == 15 THEN mM := pr1*pr2; mV := pr1*pr2*(1 - pr2); END;
    IF dId == 16 THEN mM := pr1; mV := pr1; END;
    IF dId == 17 THEN mM := 1/pr1; mV := (1 - pr1)/pr1^2; END;
    IF dId == 18 THEN
      mM := (pr1 + pr2 + pr3)/3;
      mV := (pr1^2 + pr2^2 + pr3^2 - pr1*pr2 - pr1*pr3 - pr2*pr3)/18;
    END;
    IF dId == 19 THEN mM := "undef (no mean)"; mV := "undef"; END;
    // Frechet: m + s Gamma(1 - 1/alpha) for alpha > 1;
    // s^2 (Gamma(1 - 2/alpha) - Gamma(1 - 1/alpha)^2) for alpha > 2
    IF dId == 20 THEN
      IF pr1 > 1 THEN mM := pr3 + pr2*EXP(st_lngamma(1 - 1/pr1)); ELSE mM := "undef (α<=1)"; END;
      IF pr1 > 2 THEN
        mV := MAX(0, pr2^2*(EXP(st_lngamma(1 - 2/pr1)) - EXP(2*st_lngamma(1 - 1/pr1))));
      ELSE
        mV := "undef (α<=2)";
      END;
    END;
    IF dId == 21 THEN mE2 := st_gompMom(pr1, pr2); mM := mE2(1); mV := mE2(2); END;
    // Maxwell: 2 a sqrt(2/pi), a^2 (3 pi - 8)/pi
    IF dId == 22 THEN mM := 1.59576912160573*pr1; mV := 0.453520910529675*pr1^2; END;
    IF dId == 23 THEN mM := pr1/pr2; mV := pr1/pr2^2; END;
    IF dId == 24 THEN mM := pr1*pr2; mV := pr1*pr2^2; END;
    IF dId == 25 THEN mM := pr1/(pr1 + pr2); mV := pr1*pr2/((pr1 + pr2)^2*(pr1 + pr2 + 1)); END;
    IF dId == 26 THEN mM := pr1*(1 - pr2)/pr2; mV := pr1*(1 - pr2)/pr2^2; END;
    // Hypergeometric: n K/N, n (K/N)(1 - K/N)(N - n)/(N - 1)
    IF dId == 27 THEN
      mM := pr3*pr2/pr1;
      IF pr1 > 1 THEN mV := pr3*(pr2/pr1)*(1 - pr2/pr1)*(pr1 - pr3)/(pr1 - 1); ELSE mV := 0; END;
    END;
  END;
  RETURN {mM, mV};
END;

// The parameters as one line: "μ = 0, σ = 1".
st_ParStr(dId, dPl)
BEGIN
  LOCAL psA, psB, psC;
  // v2.7d: each value carries its own " = " (st_Eq), so a sigma of 0.00005
  // reads "σ < 0.0001", not "σ = < 0.0001"
  psA := st_Eq(st_fmtx(dPl(1))); psB := st_Eq(st_fmtx(dPl(2))); psC := st_Eq(st_fmtx(dPl(3)));
  CASE
    IF dId == 1 THEN RETURN "μ" + psA + ", σ" + psB; END;
    IF dId == 2 THEN RETURN "a" + psA + ", b" + psB; END;
    IF dId == 3 THEN RETURN "rate λ" + psA; END;
    IF dId == 4 THEN RETURN "shape k" + psA + ", scale λ" + psB; END;
    IF dId == 5 THEN RETURN "μ" + psA + ", σ" + psB + " (of ln X)"; END;
    IF dId == 6 THEN RETURN "μ" + psA + ", scale b" + psB; END;
    IF dId == 7 OR dId == 8 THEN RETURN "df" + psA; END;
    IF dId == 9 THEN RETURN "df1" + psA + ", df2" + psB; END;
    IF dId == 10 THEN RETURN "C" + psA + ", K" + psB + ", S" + psC; END;
    IF dId == 11 THEN RETURN "μ" + psA + ", β" + psB; END;
    IF dId == 12 THEN RETURN "x_m" + psA + ", α" + psB; END;
    IF dId == 13 THEN RETURN "σ" + psA; END;
    IF dId == 14 THEN RETURN "μ" + psA + ", s" + psB; END;
    IF dId == 15 THEN RETURN "n" + psA + ", p" + psB; END;
    IF dId == 16 THEN RETURN "λ" + psA; END;
    IF dId == 17 THEN RETURN "p" + psA; END;
    IF dId == 18 THEN RETURN "min" + psA + ", mode" + psB + ", max" + psC; END;
    IF dId == 19 THEN RETURN "x0" + psA + ", scale" + psB; END;
    IF dId == 20 THEN RETURN "α" + psA + ", s" + psB + ", m" + psC; END;
    IF dId == 21 THEN RETURN "eta" + psA + ", b" + psB; END;
    IF dId == 22 THEN RETURN "a" + psA; END;
    IF dId == 23 THEN RETURN "k" + psA + ", λ" + psB; END;
    IF dId == 24 THEN RETURN "k" + psA + ", theta" + psB; END;
    IF dId == 25 THEN RETURN "α" + psA + ", β" + psB; END;
    IF dId == 26 THEN RETURN "r" + psA + ", p" + psB; END;
    IF dId == 27 THEN RETURN "N" + psA + ", K" + psB + ", n" + psC; END;
  END;
  RETURN "";
END;

// The default a and b for Between -- the v2.5 rules, stated in each
// bound's help line (st_ABHelp).
st_DefAB(dId, dPl)
BEGIN
  LOCAL pr1, pr2, pr3, dbM, dbMo, dbSd;
  pr1 := dPl(1); pr2 := dPl(2); pr3 := dPl(3);
  dbMo := st_Mom(dId, dPl);
  dbM := dbMo(1);
  dbSd := 0;
  IF TYPE(dbMo(2)) == 0 THEN dbSd := MAX(0, dbMo(2))^0.5; END;
  CASE
    IF dId == 1 OR dId == 6 OR dId == 14 OR dId == 19 THEN RETURN {pr1 - pr2, pr1 + pr2}; END;
    IF dId == 2 THEN RETURN {pr1, pr2}; END;
    IF dId == 3 OR dId == 4 OR dId == 13 OR dId == 22 OR dId == 23 OR dId == 24 THEN RETURN {0, dbM}; END;
    IF dId == 5 THEN RETURN {0.5, dbM}; END;
    IF dId == 7 THEN RETURN {-2, 2}; END;
    IF dId == 8 THEN RETURN {pr1/2, 2*pr1}; END;
    IF dId == 9 OR dId == 10 THEN RETURN {0.5, 2}; END;
    IF dId == 11 THEN RETURN {pr1 - pr2, pr1 + 2*pr2}; END;
    IF dId == 12 THEN RETURN {pr1, 2*pr1}; END;
    IF dId == 15 THEN RETURN {IP(dbM - dbSd), IP(dbM + dbSd)}; END;
    IF dId == 16 THEN RETURN {MAX(0, IP(pr1 - 2*pr1^0.5)), IP(pr1 + 2*pr1^0.5)}; END;
    IF dId == 17 THEN RETURN {1, IP(dbM + 2*dbSd)}; END;
    IF dId == 18 THEN RETURN {pr1, pr2}; END;
    IF dId == 20 THEN RETURN {pr3, pr3 + 2*pr2}; END;
    IF dId == 21 THEN RETURN {MAX(0, dbM - dbSd), dbM + dbSd}; END;
    IF dId == 25 THEN RETURN {MAX(0, dbM - dbSd), MIN(1, dbM + dbSd)}; END;
    IF dId == 26 THEN RETURN {0, IP(dbM)}; END;
    IF dId == 27 THEN RETURN {MAX(0, pr3 - (pr1 - pr2)), IP(dbM)}; END;
  END;
  RETURN {0, 1};
END;

// The help line under each bound on the Between form: the rule that
// st_DefAB uses (v2.5's wording for the first seventeen).
st_ABHelp(dId)
BEGIN
  CASE
    IF dId == 1 THEN RETURN {"Lower bound (default = μ - σ)", "Upper bound (default = μ + σ)"}; END;
    IF dId == 2 THEN RETURN {"Lower bound (default = a)", "Upper bound (default = b)"}; END;
    IF dId == 3 THEN RETURN {"Lower >= 0 (default = 0)", "Upper > a (default = μ = 1/λ)"}; END;
    IF dId == 4 OR dId == 13 THEN RETURN {"Lower >= 0 (default = 0)", "Upper > a (default = μ)"}; END;
    IF dId == 5 THEN RETURN {"Lower > 0 (default = 0.5)", "Upper > a (default = mean)"}; END;
    IF dId == 6 THEN RETURN {"Lower (default = μ - b)", "Upper > a (default = μ + b)"}; END;
    IF dId == 7 THEN RETURN {"Lower bound (default = -2)", "Upper > a (default = 2)"}; END;
    IF dId == 8 THEN RETURN {"Lower >= 0 (default = df/2)", "Upper > a (default = 2 df)"}; END;
    IF dId == 9 THEN RETURN {"Lower >= 0 (default = 0.5)", "Upper > a (default = 2)"}; END;
    IF dId == 10 THEN RETURN {"Lower > 0 (default = 0.5)", "Upper > a (default = 2)"}; END;
    IF dId == 11 THEN RETURN {"Lower (default = μ - β)", "Upper > a (default = μ + 2β)"}; END;
    IF dId == 12 THEN RETURN {"Lower >= x_m (default = x_m)", "Upper > a (default = 2 x_m)"}; END;
    IF dId == 14 THEN RETURN {"Lower (default = μ - s)", "Upper > a (default = μ + s)"}; END;
    IF dId == 15 THEN RETURN {"Lower int (default = μ - σ)", "Upper int (default = μ + σ)"}; END;
    IF dId == 16 THEN RETURN {"Lower int >= 0 (default = λ - 2σ)", "Upper int (default = λ + 2σ)"}; END;
    IF dId == 17 THEN RETURN {"Lower int >= 1 (default = 1)", "Upper int (default = μ + 2σ)"}; END;
    IF dId == 18 THEN RETURN {"Lower (default = min)", "Upper > a (default = mode)"}; END;
    IF dId == 19 THEN RETURN {"Lower (default = x0 - scale)", "Upper (default = x0 + scale)"}; END;
    IF dId == 20 THEN RETURN {"Lower >= m (default = m)", "Upper > a (default = m + 2s)"}; END;
    IF dId == 21 THEN RETURN {"Lower >= 0 (default = mean - sd)", "Upper (default = mean + sd)"}; END;
    IF dId == 22 OR dId == 23 OR dId == 24 THEN RETURN {"Lower >= 0 (default = 0)", "Upper > a (default = mean)"}; END;
    IF dId == 25 THEN RETURN {"Lower >= 0 (default = mean - sd)", "Upper <= 1 (default = mean + sd)"}; END;
    IF dId == 26 THEN RETURN {"Lower int >= 0 (default = 0)", "Upper int (default = mean)"}; END;
    IF dId == 27 THEN RETURN {"Lower int (default = smallest k)", "Upper int (default = mean)"}; END;
  END;
  RETURN {"Lower bound", "Upper bound"};
END;


// ============================================================
// 5. THE FORMS
// Every distribution follows the same v2.5 order:
//   1. ask which probability (lower tail, upper tail, between);
//   2. INPUT the parameters -- plus x for the two tails, nothing
//      else for Between, which asks a and b on its own form with
//      defaults computed from the parameters just typed;
//   3. check the inputs (a one-line MSGBOX and back if wrong);
//   4. draw the plot.
// Each st_ask_* takes the tail (aTl), the values to show (aCur =
// {p1, p2, p3, x, a, b}: the defaults, the gallery's values, or the
// plot's own when Edit re-opens the form) and aFull: 1 = ask the
// Between bounds too, 0 = keep to the parameter form (Edit).  It
// returns {p1, p2, p3, x, a, b}, or 0 when cancelled or wrong.
// The INPUT title carries the distribution's typical use.
// ============================================================

// Returns 1=lower, 2=upper, 3=between, 0=cancel
st_get_tail_type()
BEGIN
  LOCAL idx;
  idx := 0;
  CHOOSE(idx, "Probability Type",
    "Lower tail: P(X <= a)", "Upper tail: P(X > a)",
    "Between: P(a < X < b)", "Cancel");
  IF idx == 0 OR idx == 4 THEN RETURN 0; END;
  RETURN idx;
END;

// The last step of every form: the Between bounds.  For the two tails
// a and b are simply set by the rule (they wait in case the a<X<b soft
// key is pressed later).  For Between, the "Between Bounds" form opens
// with the rule's values (aFull = 1), or they are taken as they stand
// (aFull = 0, from Edit).  A discrete bound is rounded down (with a
// box) and kept inside the support, as v2.5 did; bounds typed the
// wrong way round are swapped.
st_AskAB(dId, abPl, abX, abTl, abFull)
BEGIN
  LOCAL a_val, b_val, abD, abH, abHa, abHb;
  LOCAL abSp, abT, okI;
  abD := st_DefAB(dId, abPl);
  a_val := abD(1); b_val := abD(2);
  IF abTl == 3 AND abFull == 1 THEN
    abH := st_ABHelp(dId);
    abHa := abH(1); abHb := abH(2);
    IF st_Disc(dId) == 1 THEN
      okI := INPUT({a_val, b_val}, "Between Bounds", {"Lower k:", "Upper k:"},
        {abHa, abHb}, {abD(1), abD(2)});
      IF okI == 0 THEN RETURN 0; END;
      a_val := st_Int(a_val, "lower k");
      b_val := st_Int(b_val, "upper k");
    ELSE
      okI := INPUT({a_val, b_val}, "Between Bounds", {"Lower (a):", "Upper (b):"},
        {abHa, abHb}, {abD(1), abD(2)});
      IF okI == 0 THEN RETURN 0; END;
    END;
  END;
  IF b_val < a_val THEN abT := a_val; a_val := b_val; b_val := abT; END;
  IF st_Disc(dId) == 1 THEN
    abSp := st_Sup(dId, abPl);
    a_val := MAX(abSp(1), a_val); b_val := MIN(abSp(2), b_val);
  END;
  RETURN CONCAT(abPl, {abX, a_val, b_val});
END;

// ─── 1 NORMAL -- typical use: test scores, measurement errors ───
st_ask_normal(aTl, aCur, aFull)
BEGIN
  LOCAL mu, sg, x_val, okI;
  mu := aCur(1); sg := aCur(2); x_val := aCur(4);
  IF aTl == 3 THEN
    okI := INPUT({mu, sg}, "Normal: test scores, errors",
      {"Mean (μ):", "Std Dev (σ):"},
      {"Location", "Scale > 0"}, {0, 1});
    IF okI == 0 THEN RETURN 0; END;
  ELSE
    okI := INPUT({mu, sg, x_val}, "Normal: test scores, errors",
      {"Mean (μ):", "Std Dev (σ):", "X value:"},
      {"Location", "Scale > 0", "Value to evaluate"},
      {0, 1, 0});
    IF okI == 0 THEN RETURN 0; END;
  END;
  IF sg <= 0 THEN MSGBOX("Error: Std Dev must be > 0"); RETURN 0; END;
  RETURN st_AskAB(1, {mu, sg, 0}, x_val, aTl, aFull);
END;

// ─── 2 UNIFORM -- typical use: time of a crime within a shift ───
st_ask_uniform(aTl, aCur, aFull)
BEGIN
  LOCAL a_param, b_param, x_val, okI;
  a_param := aCur(1); b_param := aCur(2); x_val := aCur(4);
  IF aTl == 3 THEN
    okI := INPUT({a_param, b_param}, "Uniform: time within a shift",
      {"Min (a):", "Max (b):"},
      {"Lower bound", "Upper bound > a"}, {0, 1});
    IF okI == 0 THEN RETURN 0; END;
  ELSE
    okI := INPUT({a_param, b_param, x_val}, "Uniform: time within a shift",
      {"Min (a):", "Max (b):", "X value:"},
      {"Lower bound", "Upper bound > a", "Value to evaluate"},
      {0, 1, 0.5});
    IF okI == 0 THEN RETURN 0; END;
  END;
  IF b_param <= a_param THEN MSGBOX("Error: b must be > a"); RETURN 0; END;
  RETURN st_AskAB(2, {a_param, b_param, 0}, x_val, aTl, aFull);
END;

// ─── 3 EXPONENTIAL -- typical use: time to the next incident ───
st_ask_exponential(aTl, aCur, aFull)
BEGIN
  LOCAL lambda, x_val, okI;
  lambda := aCur(1); x_val := aCur(4);
  IF aTl == 3 THEN
    okI := INPUT({lambda}, "Exponential: time to incident",
      {"Rate (λ):"},
      {"Rate > 0"}, {1});
    IF okI == 0 THEN RETURN 0; END;
  ELSE
    okI := INPUT({lambda, x_val}, "Exponential: time to incident",
      {"Rate (λ):", "X value:"},
      {"Rate > 0", "Value >= 0"}, {1, 0});
    IF okI == 0 THEN RETURN 0; END;
    IF x_val < 0 THEN MSGBOX("Error: x must be >= 0"); RETURN 0; END;
  END;
  IF lambda <= 0 THEN MSGBOX("Error: λ must be > 0"); RETURN 0; END;
  RETURN st_AskAB(3, {lambda, 0, 0}, x_val, aTl, aFull);
END;

// ─── 4 WEIBULL (shape k, scale) -- typical use: time to re-offence ───
// v2.5: the scale is stS (S is a Home global).  Weibull is on both the
// Common and the Advanced menu; that is deliberate.
st_ask_weibull(aTl, aCur, aFull)
BEGIN
  LOCAL ksh, stS, x_val, okI;
  ksh := aCur(1); stS := aCur(2); x_val := aCur(4);
  IF aTl == 3 THEN
    okI := INPUT({ksh, stS}, "Weibull: time to re-offence",
      {"Shape k:", "Scale λ:"},
      {"k > 0", "λ > 0"}, {2, 1});
    IF okI == 0 THEN RETURN 0; END;
  ELSE
    okI := INPUT({ksh, stS, x_val}, "Weibull: time to re-offence",
      {"Shape k:", "Scale λ:", "X value:"},
      {"k > 0", "λ > 0", "x >= 0"}, {2, 1, 1});
    IF okI == 0 THEN RETURN 0; END;
    IF x_val < 0 THEN MSGBOX("Error: x must be >= 0"); RETURN 0; END;
  END;
  IF ksh <= 0 THEN MSGBOX("Error: k must be > 0"); RETURN 0; END;
  IF stS <= 0 THEN MSGBOX("Error: λ must be > 0"); RETURN 0; END;
  RETURN st_AskAB(4, {ksh, stS, 0}, x_val, aTl, aFull);
END;

// ─── 5 LOGNORMAL (μ, σ of the underlying normal) -- sentence lengths ───
st_ask_lognormal(aTl, aCur, aFull)
BEGIN
  LOCAL mu, sg, x_val, okI;
  mu := aCur(1); sg := aCur(2); x_val := aCur(4);
  IF aTl == 3 THEN
    okI := INPUT({mu, sg}, "LogNormal: sentence lengths",
      {"μ (log-mean):", "σ (log-sd):"},
      {"μ real", "σ > 0"}, {0, 1});
    IF okI == 0 THEN RETURN 0; END;
  ELSE
    okI := INPUT({mu, sg, x_val}, "LogNormal: sentence lengths",
      {"μ (log-mean):", "σ (log-sd):", "X value:"},
      {"μ real", "σ > 0", "x > 0"}, {0, 1, 1});
    IF okI == 0 THEN RETURN 0; END;
    IF x_val <= 0 THEN MSGBOX("Error: x must be > 0"); RETURN 0; END;
  END;
  IF sg <= 0 THEN MSGBOX("Error: σ must be > 0"); RETURN 0; END;
  RETURN st_AskAB(5, {mu, sg, 0}, x_val, aTl, aFull);
END;

// ─── 6 LAPLACE (location μ, scale b) -- age-estimate errors ───
st_ask_laplace(aTl, aCur, aFull)
BEGIN
  LOCAL mu, bb, x_val, okI;
  mu := aCur(1); bb := aCur(2); x_val := aCur(4);
  IF aTl == 3 THEN
    okI := INPUT({mu, bb}, "Laplace: age-estimate errors",
      {"Location μ:", "Scale b:"},
      {"μ real", "b > 0"}, {0, 1});
    IF okI == 0 THEN RETURN 0; END;
  ELSE
    okI := INPUT({mu, bb, x_val}, "Laplace: age-estimate errors",
      {"Location μ:", "Scale b:", "X value:"},
      {"μ real", "b > 0", "Value"}, {0, 1, 0});
    IF okI == 0 THEN RETURN 0; END;
  END;
  IF bb <= 0 THEN MSGBOX("Error: b must be > 0"); RETURN 0; END;
  RETURN st_AskAB(6, {mu, bb, 0}, x_val, aTl, aFull);
END;

// ─── 7 STUDENT'S t -- typical use: small-sample means ───
st_ask_student(aTl, aCur, aFull)
BEGIN
  LOCAL df, x_val, okI;
  df := aCur(1); x_val := aCur(4);
  IF aTl == 3 THEN
    okI := INPUT({df}, "Student t: small-sample means",
      {"df:"},
      {"df > 0"}, {10});
    IF okI == 0 THEN RETURN 0; END;
  ELSE
    okI := INPUT({df, x_val}, "Student t: small-sample means",
      {"df:", "X value:"},
      {"df > 0", "Value to evaluate"}, {10, 0});
    IF okI == 0 THEN RETURN 0; END;
  END;
  IF df <= 0 THEN MSGBOX("Error: df must be > 0"); RETURN 0; END;
  RETURN st_AskAB(7, {df, 0, 0}, x_val, aTl, aFull);
END;

// ─── 8 CHI-SQUARE -- typical use: goodness of fit, crosstabs ───
st_ask_chisquare(aTl, aCur, aFull)
BEGIN
  LOCAL df, x_val, okI;
  df := aCur(1); x_val := aCur(4);
  IF aTl == 3 THEN
    // v2.7b: "df:" (the emulator cut "Degrees of freedom:" to "Degrees
    // of fr..."); any df > 0, whole or not
    okI := INPUT({df}, "Chi-square: goodness of fit",
      {"df:"},
      {"0 < df <= 5000 (need not be whole)"}, {10});
    IF okI == 0 THEN RETURN 0; END;
  ELSE
    okI := INPUT({df, x_val}, "Chi-square: goodness of fit",
      {"df:", "X value:"},
      {"0 < df <= 5000 (need not be whole)", "Value >= 0"}, {10, 5});
    IF okI == 0 THEN RETURN 0; END;
    IF x_val < 0 THEN MSGBOX("Error: x must be >= 0"); RETURN 0; END;
  END;
  IF df <= 0 THEN MSGBOX("Error: df must be > 0"); RETURN 0; END;
  // v2.7e (B1): the program's own incomplete gamma (st_gammaP, and
  // st_gammaQ for the right tail) works with the shape df/2 and stops
  // after 500 passes; a shape of 5000 needs about 490, and from about
  // 10000 on the answer is quietly wrong (Gamma(20000, 1) at its mean
  // gave 0.5007, not 0.5009).  So df is kept at most 5000, the Gamma's
  // bound for k, with room to spare -- no classroom table goes near it.
  IF df > 5000 THEN MSGBOX("Error: df must be <= 5000 here" + CHAR(10) + "(the incomplete gamma stops after 500 steps)"); RETURN 0; END;
  RETURN st_AskAB(8, {df, 0, 0}, x_val, aTl, aFull);
END;

// ─── 9 FISHER F -- typical use: ANOVA, a ratio of two variances ───
st_ask_fisher(aTl, aCur, aFull)
BEGIN
  LOCAL df1, df2, x_val, okI;
  df1 := aCur(1); df2 := aCur(2); x_val := aCur(4);
  IF aTl == 3 THEN
    // v2.7b: short labels ("Denominator df:" would be cut); the help
    // line says which df is which
    okI := INPUT({df1, df2}, "Fisher F: ANOVA variances",
      {"df1:", "df2:"},
      {"Numerator df > 0", "Denominator df > 0"}, {10, 15});
    IF okI == 0 THEN RETURN 0; END;
  ELSE
    okI := INPUT({df1, df2, x_val}, "Fisher F: ANOVA variances",
      {"df1:", "df2:", "X value:"},
      {"Numerator df > 0", "Denominator df > 0", "Value >= 0"}, {10, 15, 1});
    IF okI == 0 THEN RETURN 0; END;
    IF x_val < 0 THEN MSGBOX("Error: x must be >= 0"); RETURN 0; END;
  END;
  IF df1 <= 0 THEN MSGBOX("Error: df1 must be > 0"); RETURN 0; END;
  IF df2 <= 0 THEN MSGBOX("Error: df2 must be > 0"); RETURN 0; END;
  RETURN st_AskAB(9, {df1, df2, 0}, x_val, aTl, aFull);
END;

// ─── 10 BURR TYPE XII -- typical use: time to re-arrest ───
// v2.5: the shapes and scale are stC, stK, stS -- C, K and S are
// Home's built-in global reals.
st_ask_burr(aTl, aCur, aFull)
BEGIN
  LOCAL stC, stK, stS, x_val, okI;
  stC := aCur(1); stK := aCur(2); stS := aCur(3); x_val := aCur(4);
  IF aTl == 3 THEN
    okI := INPUT({stC, stK, stS}, "Burr XII: time to re-arrest",
      {"Shape C:", "Shape K:", "Scale S:"},
      {"C > 0", "K > 0", "S > 0"}, {2, 2, 1});
    IF okI == 0 THEN RETURN 0; END;
  ELSE
    okI := INPUT({stC, stK, stS, x_val}, "Burr XII: time to re-arrest",
      {"Shape C:", "Shape K:", "Scale S:", "X value:"},
      {"C > 0", "K > 0", "S > 0", "x > 0"}, {2, 2, 1, 1});
    IF okI == 0 THEN RETURN 0; END;
    IF x_val <= 0 THEN MSGBOX("Error: x must be > 0"); RETURN 0; END;
  END;
  IF stC <= 0 THEN MSGBOX("Error: C must be > 0"); RETURN 0; END;
  IF stK <= 0 THEN MSGBOX("Error: K must be > 0"); RETURN 0; END;
  IF stS <= 0 THEN MSGBOX("Error: S must be > 0"); RETURN 0; END;
  RETURN st_AskAB(10, {stC, stK, stS}, x_val, aTl, aFull);
END;

// ─── 11 GUMBEL (max) -- typical use: the worst day of the year ───
st_ask_gumbel(aTl, aCur, aFull)
BEGIN
  LOCAL mu, be, x_val, okI;
  mu := aCur(1); be := aCur(2); x_val := aCur(4);
  IF aTl == 3 THEN
    okI := INPUT({mu, be}, "Gumbel (max): worst day",
      {"Location μ:", "Scale β:"},
      {"μ real", "β > 0"}, {0, 1});
    IF okI == 0 THEN RETURN 0; END;
  ELSE
    okI := INPUT({mu, be, x_val}, "Gumbel (max): worst day",
      {"Location μ:", "Scale β:", "X value:"},
      {"μ real", "β > 0", "Value"}, {0, 1, 0});
    IF okI == 0 THEN RETURN 0; END;
  END;
  IF be <= 0 THEN MSGBOX("Error: β must be > 0"); RETURN 0; END;
  RETURN st_AskAB(11, {mu, be, 0}, x_val, aTl, aFull);
END;

// ─── 12 PARETO -- typical use: a few offenders commit most crime ───
st_ask_pareto(aTl, aCur, aFull)
BEGIN
  LOCAL xm, al, x_val, okI;
  xm := aCur(1); al := aCur(2); x_val := aCur(4);
  IF aTl == 3 THEN
    okI := INPUT({xm, al}, "Pareto: offending (80/20)",
      {"Scale x_m:", "Shape α:"},
      {"x_m > 0", "α > 0"}, {1, 3});
    IF okI == 0 THEN RETURN 0; END;
  ELSE
    okI := INPUT({xm, al, x_val}, "Pareto: offending (80/20)",
      {"Scale x_m:", "Shape α:", "X value:"},
      {"x_m > 0", "α > 0", "x >= x_m"}, {1, 3, 2});
    IF okI == 0 THEN RETURN 0; END;
  END;
  IF xm <= 0 THEN MSGBOX("Error: x_m must be > 0"); RETURN 0; END;
  IF al <= 0 THEN MSGBOX("Error: α must be > 0"); RETURN 0; END;
  IF aTl <> 3 AND x_val < xm THEN MSGBOX("Error: x must be >= x_m"); RETURN 0; END;
  RETURN st_AskAB(12, {xm, al, 0}, x_val, aTl, aFull);
END;

// ─── 13 RAYLEIGH -- typical use: a shot's distance from the aim point ───
st_ask_rayleigh(aTl, aCur, aFull)
BEGIN
  LOCAL sg, x_val, okI;
  sg := aCur(1); x_val := aCur(4);
  IF aTl == 3 THEN
    okI := INPUT({sg}, "Rayleigh: shot miss distance",
      {"Scale σ:"}, {"σ > 0"}, {1.5});
    IF okI == 0 THEN RETURN 0; END;
  ELSE
    okI := INPUT({sg, x_val}, "Rayleigh: shot miss distance",
      {"Scale σ:", "X value:"}, {"σ > 0", "x >= 0"}, {1.5, 2});
    IF okI == 0 THEN RETURN 0; END;
    IF x_val < 0 THEN MSGBOX("Error: x must be >= 0"); RETURN 0; END;
  END;
  IF sg <= 0 THEN MSGBOX("Error: σ must be > 0"); RETURN 0; END;
  RETURN st_AskAB(13, {sg, 0, 0}, x_val, aTl, aFull);
END;

// ─── 14 LOGISTIC -- typical use: risk-assessment scores ───
st_ask_logistic(aTl, aCur, aFull)
BEGIN
  LOCAL mu, sc, x_val, okI;
  mu := aCur(1); sc := aCur(2); x_val := aCur(4);
  IF aTl == 3 THEN
    okI := INPUT({mu, sc}, "Logistic: risk scores",
      {"Location μ:", "Scale s:"},
      {"μ real", "s > 0"}, {0, 1});
    IF okI == 0 THEN RETURN 0; END;
  ELSE
    okI := INPUT({mu, sc, x_val}, "Logistic: risk scores",
      {"Location μ:", "Scale s:", "X value:"},
      {"μ real", "s > 0", "Value"}, {0, 1, 0});
    IF okI == 0 THEN RETURN 0; END;
  END;
  IF sc <= 0 THEN MSGBOX("Error: s must be > 0"); RETURN 0; END;
  RETURN st_AskAB(14, {mu, sc, 0}, x_val, aTl, aFull);
END;

// ─── 15 BINOMIAL -- typical use: convictions in n cases ───
// v2.5: n and p are stN and stP.  A fraction typed for n or k is
// rounded down, and a one-line box says so.
st_ask_binomial(aTl, aCur, aFull)
BEGIN
  LOCAL stN, stP, k_val, okI;
  stN := aCur(1); stP := aCur(2); k_val := aCur(4);
  IF aTl == 3 THEN
    okI := INPUT({stN, stP}, "Binomial: convictions in n",
      {"n (trials):", "p (success):"},
      {"Trials", "0 < p < 1"}, {20, 0.5});
    IF okI == 0 THEN RETURN 0; END;
  ELSE
    okI := INPUT({stN, stP, k_val}, "Binomial: convictions in n",
      {"n (trials):", "p (success):", "k (successes):"},
      {"Trials", "0 < p < 1", "Successes"}, {20, 0.5, 10});
    IF okI == 0 THEN RETURN 0; END;
  END;
  stN := st_Int(stN, "n");
  IF aTl <> 3 THEN k_val := st_Int(k_val, "k"); END;
  IF stN < 1 THEN MSGBOX("Error: n must be >= 1"); RETURN 0; END;
  IF stP <= 0 OR stP >= 1 THEN MSGBOX("Error: 0 < p < 1"); RETURN 0; END;
  IF aTl <> 3 AND (k_val < 0 OR k_val > stN) THEN MSGBOX("Error: 0 <= k <= n"); RETURN 0; END;
  RETURN st_AskAB(15, {stN, stP, 0}, k_val, aTl, aFull);
END;

// ─── 16 POISSON -- typical use: calls per shift ───
st_ask_poisson(aTl, aCur, aFull)
BEGIN
  LOCAL lambda, k_val, okI;
  lambda := aCur(1); k_val := aCur(4);
  IF aTl == 3 THEN
    okI := INPUT({lambda}, "Poisson: calls per shift",
      {"λ (rate):"},
      {"Rate > 0"}, {5});
    IF okI == 0 THEN RETURN 0; END;
  ELSE
    okI := INPUT({lambda, k_val}, "Poisson: calls per shift",
      {"λ (rate):", "k (events):"},
      {"Rate > 0", "Events"}, {5, 3});
    IF okI == 0 THEN RETURN 0; END;
    k_val := st_Int(k_val, "k");
    IF k_val < 0 THEN MSGBOX("Error: k must be >= 0"); RETURN 0; END;
  END;
  IF lambda <= 0 THEN MSGBOX("Error: λ must be > 0"); RETURN 0; END;
  RETURN st_AskAB(16, {lambda, 0, 0}, k_val, aTl, aFull);
END;

// ─── 17 GEOMETRIC (k = the trial of the first success, k >= 1) ───
//     typical use: searches until the first hit.  v2.5: p is stP.
st_ask_geometric(aTl, aCur, aFull)
BEGIN
  LOCAL stP, k_val, okI;
  stP := aCur(1); k_val := aCur(4);
  IF aTl == 3 THEN
    okI := INPUT({stP}, "Geometric: searches to 1st hit",
      {"p (success):"},
      {"0 < p < 1"}, {0.3});
    IF okI == 0 THEN RETURN 0; END;
  ELSE
    okI := INPUT({stP, k_val}, "Geometric: searches to 1st hit",
      {"p (success):", "k (trial):"},
      {"0 < p < 1", "Trial of 1st success (>=1)"}, {0.3, 3});
    IF okI == 0 THEN RETURN 0; END;
    k_val := st_Int(k_val, "k");
    IF k_val < 1 THEN MSGBOX("Error: k must be >= 1"); RETURN 0; END;
  END;
  IF stP <= 0 OR stP >= 1 THEN MSGBOX("Error: 0 < p < 1"); RETURN 0; END;
  RETURN st_AskAB(17, {stP, 0, 0}, k_val, aTl, aFull);
END;

// ─── 18 TRIANGULAR (min, mode, max) -- v2.6, was a stub ───
//     typical use: an expert's rough guess, e.g. hours since death
//     "at least 2, most likely 5, at most 12".
st_ask_triangular(aTl, aCur, aFull)
BEGIN
  LOCAL trA, trC, trB, x_val, okI;
  trA := aCur(1); trC := aCur(2); trB := aCur(3); x_val := aCur(4);
  IF aTl == 3 THEN
    okI := INPUT({trA, trC, trB}, "Triangular: min, likeliest, max",
      {"Min (a):", "Mode (c):", "Max (b):"},
      {"Smallest possible value", "Most likely value", "Largest possible value"},
      {2, 5, 12});
    IF okI == 0 THEN RETURN 0; END;
  ELSE
    okI := INPUT({trA, trC, trB, x_val}, "Triangular: min, likeliest, max",
      {"Min (a):", "Mode (c):", "Max (b):", "X value:"},
      {"Smallest possible value", "Most likely value", "Largest possible value", "Value to evaluate"},
      {2, 5, 12, 5});
    IF okI == 0 THEN RETURN 0; END;
  END;
  IF trB <= trA OR trC < trA OR trC > trB THEN
    MSGBOX("Error: need min <= mode <= max, min < max"); RETURN 0;
  END;
  RETURN st_AskAB(18, {trA, trC, trB}, x_val, aTl, aFull);
END;

// ─── 19 CAUCHY (location x0, scale) -- v2.6, was a stub ───
//     typical use: the ratio of two measurement errors.  It has NO
//     mean: the average of Cauchy values never settles down.
st_ask_cauchy(aTl, aCur, aFull)
BEGIN
  LOCAL cLoc, cScl, x_val, okI;
  cLoc := aCur(1); cScl := aCur(2); x_val := aCur(4);
  IF aTl == 3 THEN
    okI := INPUT({cLoc, cScl}, "Cauchy: ratio of errors",
      {"Location (x0):", "Scale:"},
      {"The median", "Scale > 0 (half the IQR)"}, {0, 1});
    IF okI == 0 THEN RETURN 0; END;
  ELSE
    okI := INPUT({cLoc, cScl, x_val}, "Cauchy: ratio of errors",
      {"Location (x0):", "Scale:", "X value:"},
      {"The median", "Scale > 0 (half the IQR)", "Value to evaluate"},
      {0, 1, 1});
    IF okI == 0 THEN RETURN 0; END;
  END;
  IF cScl <= 0 THEN MSGBOX("Error: scale must be > 0"); RETURN 0; END;
  RETURN st_AskAB(19, {cLoc, cScl, 0}, x_val, aTl, aFull);
END;

// ─── 20 FRECHET (shape α, scale s, location m) -- v2.6, was a stub ───
//     typical use: the largest of many losses (heavy right tail).
st_ask_frechet(aTl, aCur, aFull)
BEGIN
  LOCAL frA, frS, frM, x_val, okI;
  frA := aCur(1); frS := aCur(2); frM := aCur(3); x_val := aCur(4);
  IF aTl == 3 THEN
    okI := INPUT({frA, frS, frM}, "Frechet: largest losses",
      {"Shape (α):", "Scale (s):", "Location (m):"},
      {"α > 0 (mean needs α > 1)", "s > 0", "Smallest possible value"},
      {3, 1, 0});
    IF okI == 0 THEN RETURN 0; END;
  ELSE
    okI := INPUT({frA, frS, frM, x_val}, "Frechet: largest losses",
      {"Shape (α):", "Scale (s):", "Location (m):", "X value:"},
      {"α > 0 (mean needs α > 1)", "s > 0", "Smallest possible value", "Value to evaluate"},
      {3, 1, 0, 1});
    IF okI == 0 THEN RETURN 0; END;
  END;
  IF frA <= 0 THEN MSGBOX("Error: α must be > 0"); RETURN 0; END;
  IF frS <= 0 THEN MSGBOX("Error: s must be > 0"); RETURN 0; END;
  RETURN st_AskAB(20, {frA, frS, frM}, x_val, aTl, aFull);
END;

// ─── 21 GOMPERTZ (shape eta, rate b) -- v2.6, was a stub ───
//     typical use: years lived past age 20 (the Gompertz law: the
//     death rate b eta e^(b x) doubles every ln(2)/b years; b = 0.085
//     is about 8 years, close to human adult mortality).
st_ask_gompertz(aTl, aCur, aFull)
BEGIN
  LOCAL gmEt, gmBr, x_val, okI;
  gmEt := aCur(1); gmBr := aCur(2); x_val := aCur(4);
  IF aTl == 3 THEN
    okI := INPUT({gmEt, gmBr}, "Gompertz: years lived past 20",
      {"Shape (eta):", "Rate (b):"},
      {"eta > 0 (smaller = longer lives)", "b > 0 (doubling time ln2/b)"},
      {0.01, 0.085});
    IF okI == 0 THEN RETURN 0; END;
  ELSE
    okI := INPUT({gmEt, gmBr, x_val}, "Gompertz: years lived past 20",
      {"Shape (eta):", "Rate (b):", "X value:"},
      {"eta > 0 (smaller = longer lives)", "b > 0 (doubling time ln2/b)", "x >= 0"},
      {0.01, 0.085, 50});
    IF okI == 0 THEN RETURN 0; END;
    IF x_val < 0 THEN MSGBOX("Error: x must be >= 0"); RETURN 0; END;
  END;
  IF gmEt <= 0 THEN MSGBOX("Error: eta must be > 0"); RETURN 0; END;
  IF gmBr <= 0 THEN MSGBOX("Error: b must be > 0"); RETURN 0; END;
  RETURN st_AskAB(21, {gmEt, gmBr, 0}, x_val, aTl, aFull);
END;

// ─── 22 MAXWELL-BOLTZMANN (scale a) -- v2.6, was a stub ───
//     typical use: a 3-D miss distance (the length of a vector of three
//     independent normal errors), or a gas molecule's speed.
st_ask_maxwell(aTl, aCur, aFull)
BEGIN
  LOCAL mxA, x_val, okI;
  mxA := aCur(1); x_val := aCur(4);
  IF aTl == 3 THEN
    okI := INPUT({mxA}, "Maxwell: 3-D miss distance",
      {"Scale (a):"}, {"a > 0 (the sd of each axis)"}, {1});
    IF okI == 0 THEN RETURN 0; END;
  ELSE
    okI := INPUT({mxA, x_val}, "Maxwell: 3-D miss distance",
      {"Scale (a):", "X value:"}, {"a > 0 (the sd of each axis)", "x >= 0"},
      {1, 1.5});
    IF okI == 0 THEN RETURN 0; END;
    IF x_val < 0 THEN MSGBOX("Error: x must be >= 0"); RETURN 0; END;
  END;
  IF mxA <= 0 THEN MSGBOX("Error: a must be > 0"); RETURN 0; END;
  RETURN st_AskAB(22, {mxA, 0, 0}, x_val, aTl, aFull);
END;

// ─── 23 ERLANG (whole k, rate λ) -- v2.6, was a stub ───
//     typical use: time until the k-th call (k exponential waits).
st_ask_erlang(aTl, aCur, aFull)
BEGIN
  LOCAL erK, erL, x_val, okI;
  erK := aCur(1); erL := aCur(2); x_val := aCur(4);
  IF aTl == 3 THEN
    okI := INPUT({erK, erL}, "Erlang: time to the k-th call",
      {"Shape k:", "Rate (λ):"},
      {"Whole number >= 1", "λ > 0"}, {3, 1});
    IF okI == 0 THEN RETURN 0; END;
  ELSE
    okI := INPUT({erK, erL, x_val}, "Erlang: time to the k-th call",
      {"Shape k:", "Rate (λ):", "X value:"},
      {"Whole number >= 1", "λ > 0", "x >= 0"}, {3, 1, 2});
    IF okI == 0 THEN RETURN 0; END;
    IF x_val < 0 THEN MSGBOX("Error: x must be >= 0"); RETURN 0; END;
  END;
  erK := st_Int(erK, "k");
  IF erK < 1 THEN MSGBOX("Error: k must be >= 1"); RETURN 0; END;
  IF erL <= 0 THEN MSGBOX("Error: λ must be > 0"); RETURN 0; END;
  RETURN st_AskAB(23, {erK, erL, 0}, x_val, aTl, aFull);
END;

// ─── 24 GAMMA (shape k, scale theta) -- v2.6, was a stub ───
//     typical use: case processing time.  Any k > 0 (k need not be a
//     whole number: the CDF is our own incomplete gamma).
st_ask_gamma(aTl, aCur, aFull)
BEGIN
  LOCAL gaK, gaT, x_val, okI;
  gaK := aCur(1); gaT := aCur(2); x_val := aCur(4);
  IF aTl == 3 THEN
    okI := INPUT({gaK, gaT}, "Gamma: case processing time",
      {"Shape (k):", "Scale (theta):"},
      {"0 < k <= 5000 (any, not just whole)", "theta > 0"}, {2.5, 2});
    IF okI == 0 THEN RETURN 0; END;
  ELSE
    okI := INPUT({gaK, gaT, x_val}, "Gamma: case processing time",
      {"Shape (k):", "Scale (theta):", "X value:"},
      {"0 < k <= 5000 (any, not just whole)", "theta > 0", "x >= 0"}, {2.5, 2, 4});
    IF okI == 0 THEN RETURN 0; END;
    IF x_val < 0 THEN MSGBOX("Error: x must be >= 0"); RETURN 0; END;
  END;
  IF gaK <= 0 THEN MSGBOX("Error: k must be > 0"); RETURN 0; END;
  // v2.7e (B1): k at most 5000, as the chi-square's df: the program's
  // own incomplete gamma stops after 500 passes, and from about
  // k = 10000 that is too few (Gamma(20000, 1) at its mean: 0.5007,
  // where 0.5009 is right)
  IF gaK > 5000 THEN MSGBOX("Error: k must be <= 5000 here" + CHAR(10) + "(the incomplete gamma stops after 500 steps)"); RETURN 0; END;
  IF gaT <= 0 THEN MSGBOX("Error: theta must be > 0"); RETURN 0; END;
  RETURN st_AskAB(24, {gaK, gaT, 0}, x_val, aTl, aFull);
END;

// ─── 25 BETA (shapes α, β) -- v2.6, was a stub ───
//     typical use: proportions, such as clearance rates.
st_ask_beta(aTl, aCur, aFull)
BEGIN
  LOCAL btA, btB, x_val, okI;
  btA := aCur(1); btB := aCur(2); x_val := aCur(4);
  IF aTl == 3 THEN
    okI := INPUT({btA, btB}, "Beta: clearance rates",
      {"Shape (α):", "Shape (β):"},
      {"α > 0", "β > 0"}, {2, 5});
    IF okI == 0 THEN RETURN 0; END;
  ELSE
    okI := INPUT({btA, btB, x_val}, "Beta: clearance rates",
      {"Shape (α):", "Shape (β):", "X value:"},
      {"α > 0", "β > 0", "0 <= x <= 1"}, {2, 5, 0.3});
    IF okI == 0 THEN RETURN 0; END;
    IF x_val < 0 OR x_val > 1 THEN MSGBOX("Error: 0 <= x <= 1"); RETURN 0; END;
  END;
  IF btA <= 0 OR btB <= 0 THEN MSGBOX("Error: α and β must be > 0"); RETURN 0; END;
  RETURN st_AskAB(25, {btA, btB, 0}, x_val, aTl, aFull);
END;

// ─── 26 NEGATIVE BINOMIAL (r, p; k = failures before the r-th success) ───
//     v2.6, was a stub.  Typical use: stops before the r-th arrest.
st_ask_negbinom(aTl, aCur, aFull)
BEGIN
  LOCAL nbR, nbP, k_val, okI;
  nbR := aCur(1); nbP := aCur(2); k_val := aCur(4);
  IF aTl == 3 THEN
    okI := INPUT({nbR, nbP}, "Neg. binomial: stops, r arrests",
      {"r (successes):", "p (success):"},
      {"Whole number >= 1", "0 < p < 1"}, {3, 0.4});
    IF okI == 0 THEN RETURN 0; END;
  ELSE
    okI := INPUT({nbR, nbP, k_val}, "Neg. binomial: stops, r arrests",
      {"r (successes):", "p (success):", "k (failures):"},
      {"Whole number >= 1", "0 < p < 1", "Failures before the r-th success"},
      {3, 0.4, 4});
    IF okI == 0 THEN RETURN 0; END;
    k_val := st_Int(k_val, "k");
    IF k_val < 0 THEN MSGBOX("Error: k must be >= 0"); RETURN 0; END;
  END;
  nbR := st_Int(nbR, "r");
  IF nbR < 1 THEN MSGBOX("Error: r must be >= 1"); RETURN 0; END;
  IF nbP <= 0 OR nbP >= 1 THEN MSGBOX("Error: 0 < p < 1"); RETURN 0; END;
  RETURN st_AskAB(26, {nbR, nbP, 0}, k_val, aTl, aFull);
END;

// ─── 27 HYPERGEOMETRIC (N in all, K of them "positive", n drawn) ───
//     v2.6, was a stub.  Typical use: sampling seized packets without
//     replacement -- if only 40 of 50 packets hold the drug, how
//     likely is it that all 6 tested are positive?  P(X > 5) = 0.2415.
st_ask_hypergeom(aTl, aCur, aFull)
BEGIN
  LOCAL hgN, hgK, hgM, k_val, okI;
  hgN := aCur(1); hgK := aCur(2); hgM := aCur(3); k_val := aCur(4);
  IF aTl == 3 THEN
    // v2.7b: shorter labels (the long ones would be cut); the help
    // lines say what each number is
    okI := INPUT({hgN, hgK, hgM}, "Hypergeometric: seized packets",
      {"N (pop.):", "K (pos. in N):", "n (drawn):"},
      {"Population size, 1 to 1000", "Positives in N, 0 to N", "Drawn, 1 to N"},
      {50, 40, 6});
    IF okI == 0 THEN RETURN 0; END;
  ELSE
    okI := INPUT({hgN, hgK, hgM, k_val}, "Hypergeometric: seized packets",
      {"N (pop.):", "K (pos. in N):", "n (drawn):", "k (pos. in n):"},
      {"Population size, 1 to 1000", "Positives in N, 0 to N", "Drawn, 1 to N", "Positives drawn, 0 <= k <= n"},
      {50, 40, 6, 5});
    IF okI == 0 THEN RETURN 0; END;
    k_val := st_Int(k_val, "k");
  END;
  hgN := st_Int(hgN, "N"); hgK := st_Int(hgK, "K"); hgM := st_Int(hgM, "n");
  IF hgN < 1 OR hgN > 1000 THEN MSGBOX("Error: N from 1 to 1000"); RETURN 0; END;
  IF hgK < 0 OR hgK > hgN THEN MSGBOX("Error: 0 <= K <= N"); RETURN 0; END;
  IF hgM < 1 OR hgM > hgN THEN MSGBOX("Error: 1 <= n <= N"); RETURN 0; END;
  IF aTl <> 3 AND (k_val < 0 OR k_val > hgM) THEN MSGBOX("Error: 0 <= k <= n"); RETURN 0; END;
  RETURN st_AskAB(27, {hgN, hgK, hgM}, k_val, aTl, aFull);
END;

// ---- the form for distribution dId: one place that knows them all ----
st_Ask(dId, aTl, aCur, aFull)
BEGIN
  CASE
    IF dId == 1 THEN RETURN st_ask_normal(aTl, aCur, aFull); END;
    IF dId == 2 THEN RETURN st_ask_uniform(aTl, aCur, aFull); END;
    IF dId == 3 THEN RETURN st_ask_exponential(aTl, aCur, aFull); END;
    IF dId == 4 THEN RETURN st_ask_weibull(aTl, aCur, aFull); END;
    IF dId == 5 THEN RETURN st_ask_lognormal(aTl, aCur, aFull); END;
    IF dId == 6 THEN RETURN st_ask_laplace(aTl, aCur, aFull); END;
    IF dId == 7 THEN RETURN st_ask_student(aTl, aCur, aFull); END;
    IF dId == 8 THEN RETURN st_ask_chisquare(aTl, aCur, aFull); END;
    IF dId == 9 THEN RETURN st_ask_fisher(aTl, aCur, aFull); END;
    IF dId == 10 THEN RETURN st_ask_burr(aTl, aCur, aFull); END;
    IF dId == 11 THEN RETURN st_ask_gumbel(aTl, aCur, aFull); END;
    IF dId == 12 THEN RETURN st_ask_pareto(aTl, aCur, aFull); END;
    IF dId == 13 THEN RETURN st_ask_rayleigh(aTl, aCur, aFull); END;
    IF dId == 14 THEN RETURN st_ask_logistic(aTl, aCur, aFull); END;
    IF dId == 15 THEN RETURN st_ask_binomial(aTl, aCur, aFull); END;
    IF dId == 16 THEN RETURN st_ask_poisson(aTl, aCur, aFull); END;
    IF dId == 17 THEN RETURN st_ask_geometric(aTl, aCur, aFull); END;
    IF dId == 18 THEN RETURN st_ask_triangular(aTl, aCur, aFull); END;
    IF dId == 19 THEN RETURN st_ask_cauchy(aTl, aCur, aFull); END;
    IF dId == 20 THEN RETURN st_ask_frechet(aTl, aCur, aFull); END;
    IF dId == 21 THEN RETURN st_ask_gompertz(aTl, aCur, aFull); END;
    IF dId == 22 THEN RETURN st_ask_maxwell(aTl, aCur, aFull); END;
    IF dId == 23 THEN RETURN st_ask_erlang(aTl, aCur, aFull); END;
    IF dId == 24 THEN RETURN st_ask_gamma(aTl, aCur, aFull); END;
    IF dId == 25 THEN RETURN st_ask_beta(aTl, aCur, aFull); END;
    IF dId == 26 THEN RETURN st_ask_negbinom(aTl, aCur, aFull); END;
    IF dId == 27 THEN RETURN st_ask_hypergeom(aTl, aCur, aFull); END;
  END;
  RETURN 0;
END;


// ============================================================
// 6. HELP
// ============================================================

// v2.5: a drawn page; v2.6 adds every gesture and soft key (two pages).
st_show_help()
BEGIN
  st_Page("Statistical Distributions v2.7e", {
    "Probabilities for 27 distributions.",
    "Gallery: tap a family (or arrows),",
    "Menus: the tail first, then the form.",
    "",
    "ON THE PLOT",
    "Tap or drag in the plot = move x (k);",
    "  on lift x snaps to the grid.",
    "Soft keys, page 1:",
    "P<=x, P>x = which tail is shaded.",
    "a<X<b: tap a, then tap b; after that",
    "  a drag moves the nearer line.",
    "Zoom-, Zoom+ = about x.",
    "More (slot 6) = the next page:",
    "Edit = the form (in a<X<b it resets",
    "  a and b), Numbers = the numbers",
    "  page, CDF = the S-curve, Approx /",
    "  N(0,1) = the overlay (else Refit),",
    "  Menu = back to the menu.",
    "More again = page 3: Formula = the",
    "  CAS works out E[X], E[X^2], Var",
    "  (and F(x) for Exponential, Weibull,",
    "  Pareto; an exact P for discrete),",
    "  each checked against the program;",
    "  Help = these pages.",
    "Left/Right = x by one grid step.",
    "Up = upper tail, Down = lower tail;",
    "  in a<X<b, Up = move b, Down = a.",
    "+/- = zoom, 0 = refit the window.",
    "Ent = options: numbers page, Formula",
    "  page, CDF, normal approx., N(0,1)",
    "  behind t.",
    "Esc = back to the menu.",
    "",
    "THE FORMULA PAGE (the CAS)",
    "exact (CAS) = the CAS's answers,",
    "  each one checked; decimal = the",
    "  CAS was not used. 'No closed form",
    "  shown' = the CAS could not do it",
    "  (or its answer failed the check).",
    "",
    "COLOURS: blue = density or pmf,",
    "  pale blue = the probability,",
    "  red = x, green dots = the mean,",
    "  magenta = P and the CDF.",
    "",
    "R. Metcalf, UT Arlington CRCJ 3340"});
END;


// ============================================================
// 7. THE PLOT
// Window: stLo..stHi across the plot box stPL..stPR.  Continuous:
// pixel column c (0..stW) is x = stLo + c (stHi - stLo)/stW, and the
// caches hold entry c + 1.  Discrete: k = stLo..stHi, one bar each,
// stW/(stHi - stLo + 1) pixels wide, cache entry k - stLo + 1.
// ============================================================

// The state as one list, in the shape the forms take.
st_Cur()
BEGIN
  RETURN CONCAT(stPar, {stXv, stAv, stBv});
END;

// The pixel column of value pxV (a discrete k: the centre of its bar).
st_Px(pxV)
BEGIN
  IF stDisc == 1 THEN
    RETURN stPL + IP((pxV - stLo + 0.5)*stW/(stHi - stLo + 1));
  END;
  RETURN stPL + IP((pxV - stLo)/(stHi - stLo)*stW + 0.5);
END;

// The value under pixel column vaPx (a discrete k: the bar it is in).
st_Val(vaPx)
BEGIN
  LOCAL vaK;
  IF stDisc == 1 THEN
    vaK := stLo + IP((vaPx - stPL)*(stHi - stLo + 1)/stW);
    RETURN MIN(stHi, MAX(stLo, vaK));
  END;
  RETURN stLo + (vaPx - stPL)/stW*(stHi - stLo);
END;

// The cache entry for value coV.
st_Col(coV)
BEGIN
  IF stDisc == 1 THEN RETURN MIN(stHi - stLo + 1, MAX(1, coV - stLo + 1)); END;
  RETURN MIN(stW + 1, MAX(1, IP((coV - stLo)/(stHi - stLo)*stW + 0.5) + 1));
END;

// F(v): from the cache while a finger drags (the finger always sits on
// a cached column, so the value is exact), and for a discrete k on the
// plot; otherwise straight from st_cdf.
st_CdfAt(cfV)
BEGIN
  IF stDisc == 1 THEN
    IF cfV >= stLo AND cfV <= stHi THEN RETURN stCd(cfV - stLo + 1); END;
    RETURN st_cdf(stId, stPar, cfV);
  END;
  IF stDrg == 1 AND cfV >= stLo AND cfV <= stHi THEN RETURN stCd(st_Col(cfV)); END;
  RETURN st_cdf(stId, stPar, cfV);
END;

// P(X > v), the right tail, the same way (v2.7e): from the stSv cache
// while a finger drags and for a discrete k on the plot; otherwise
// straight from st_sf -- never as 1 - F(v), which loses the digits of a
// far right tail on the Prime's 12-digit reals.
st_SfAt(sfV)
BEGIN
  IF stDisc == 1 THEN
    IF sfV >= stLo AND sfV <= stHi THEN RETURN stSv(sfV - stLo + 1); END;
    RETURN st_sf(stId, stPar, sfV);
  END;
  IF stDrg == 1 AND sfV >= stLo AND sfV <= stHi THEN RETURN stSv(st_Col(sfV)); END;
  RETURN st_sf(stId, stPar, sfV);
END;

// f(v) the same way.
st_PdfAt(pfV)
BEGIN
  IF stDisc == 1 THEN
    IF pfV >= stLo AND pfV <= stHi THEN RETURN stPd(pfV - stLo + 1); END;
    RETURN st_pdf(stId, stPar, pfV);
  END;
  IF stDrg == 1 AND pfV >= stLo AND pfV <= stHi THEN RETURN stPd(st_Col(pfV)); END;
  RETURN st_pdf(stId, stPar, pfV);
END;

// The probability asked for.  Discrete Between counts both ends:
// P(a <= X <= b) = F(b) - F(a - 1).
// v2.7e: the right tail is st_SfAt, never 1 - F(x) (B3, B4, D9).  And a
// Between that lies in the upper half (F(a) above 1/2) is the difference
// of the two RIGHT tails, S(a) - S(b): there both are small and keep
// their digits, where F(b) - F(a) would subtract two numbers that both
// start 0.99999... -- Normal(0, 1), P(8 < X < 9) came out exactly 0,
// while P(-9 < X < -8) read "< 0.000001".
st_Prob()
BEGIN
  IF stTl == 1 THEN RETURN st_CdfAt(stXv); END;
  IF stTl == 2 THEN RETURN st_SfAt(stXv); END;
  IF stDisc == 1 THEN
    IF st_CdfAt(stAv - 1) > 0.5 THEN RETURN st_SfAt(stAv - 1) - st_SfAt(stBv); END;
    RETURN st_CdfAt(stBv) - st_CdfAt(stAv - 1);
  END;
  IF st_CdfAt(stAv) > 0.5 THEN RETURN st_SfAt(stAv) - st_SfAt(stBv); END;
  RETURN st_CdfAt(stBv) - st_CdfAt(stAv);
END;

// The x grid: ticks about six to the window, arrow steps about sixty.
st_Grid()
BEGIN
  IF stDisc == 1 THEN
    stGs := 1;
    stTk := MAX(1, st_Nice((stHi - stLo + 1)/8));
  ELSE
    stTk := st_Nice((stHi - stLo)/6);
    stGs := st_Nice((stHi - stLo)/60);
  END;
END;

// How much probability lies left and right of the window -- for the
// "beyond" notes, so a cut tail is never hidden in silence.
// v2.7e: the right side is st_sf, the right tail itself, not 1 - F.
st_CutCalc()
BEGIN
  IF stDisc == 1 THEN
    stCut := {st_cdf(stId, stPar, stLo - 1), st_sf(stId, stPar, stHi)};
  ELSE
    stCut := {st_cdf(stId, stPar, stLo), st_sf(stId, stPar, stHi)};
  END;
END;

// Keep a discrete x, a, b inside the support (after an Edit).
st_Clamp()
BEGIN
  LOCAL clS;
  IF stDisc == 1 THEN
    clS := st_Sup(stId, stPar);
    stXv := MIN(clS(2), MAX(clS(1), stXv));
    stAv := MIN(clS(2), MAX(clS(1), stAv));
    stBv := MIN(clS(2), MAX(clS(1), stBv));
  END;
END;

// ---- fit the window ----
// Continuous: from the 0.1 % to the 99.9 % quantile.  A heavy tail
// (Cauchy, Pareto, t with 1 or 2 df, LogNormal ...) would squeeze the
// whole curve into a few pixels, so each side is cut at 8
// interquartile ranges beyond the quartile, and the plot then says
// how much lies beyond.  A finite end of the support close to the
// window becomes the window's edge (an Exponential starts at 0); a
// support finite at both ends (Uniform, Triangular, Beta) is shown
// whole with a 5 % margin.  Discrete: the whole support when it has
// at most 41 values, else the 0.1 % .. 99.9 % k's, at least nine
// bars.  fiStr = 1 stretches the window to show x (or a and b), as
// Norm3 keeps x0 visible.
st_Fit(fiStr)
BEGIN
  LOCAL fiS, fiQ1, fiQ9, fiQa, fiQb, fiIq;
  LOCAL fiLo, fiHi, fiW, fiN, fiMo;
  stZm := 0;
  fiS := st_Sup(stId, stPar);
  IF stDisc == 1 THEN
    IF fiS(2) - fiS(1) <= 40 THEN
      fiLo := fiS(1); fiHi := fiS(2);
    ELSE
      fiLo := st_icdf(stId, stPar, 0.001);
      fiHi := st_icdf(stId, stPar, 0.999);
      fiN := 0;
      WHILE fiHi - fiLo < 8 AND fiN < 20 DO
        IF fiHi < fiS(2) THEN fiHi := fiHi + 1; END;
        IF fiHi - fiLo < 8 AND fiLo > fiS(1) THEN fiLo := fiLo - 1; END;
        fiN := fiN + 1;
      END;
    END;
    IF fiStr == 1 THEN
      IF stTl == 3 THEN
        fiLo := MIN(fiLo, stAv); fiHi := MAX(fiHi, stBv);
      ELSE
        fiLo := MIN(fiLo, stXv); fiHi := MAX(fiHi, stXv);
      END;
    END;
    // never more than 251 bars (a bar would be thinner than a pixel)
    IF fiHi - fiLo > 250 THEN
      fiMo := st_Mom(stId, stPar);
      fiLo := MAX(fiLo, IP(fiMo(1)) - 125);
      fiHi := fiLo + 250;
    END;
  ELSE
    IF fiS(1) > -1E299 AND fiS(2) < 1E299 THEN
      fiW := fiS(2) - fiS(1);
      fiLo := fiS(1) - 0.05*fiW; fiHi := fiS(2) + 0.05*fiW;
    ELSE
      fiQ1 := st_icdf(stId, stPar, 0.001);
      fiQ9 := st_icdf(stId, stPar, 0.999);
      fiQa := st_icdf(stId, stPar, 0.25);
      fiQb := st_icdf(stId, stPar, 0.75);
      fiIq := fiQb - fiQa;
      fiLo := MAX(fiQ1, fiQa - 8*fiIq);
      fiHi := MIN(fiQ9, fiQb + 8*fiIq);
      IF fiS(1) > -1E299 AND fiLo - fiS(1) < 0.1*(fiHi - fiLo) THEN fiLo := fiS(1); END;
    END;
    IF fiStr == 1 THEN
      fiW := fiHi - fiLo;
      IF stTl == 3 THEN
        IF stAv < fiLo THEN fiLo := stAv - 0.05*fiW; END;
        IF stBv > fiHi THEN fiHi := stBv + 0.05*fiW; END;
      ELSE
        IF stXv < fiLo THEN fiLo := stXv - 0.05*fiW; END;
        IF stXv > fiHi THEN fiHi := stXv + 0.05*fiW; END;
      END;
    END;
    IF fiHi <= fiLo THEN fiHi := fiLo + 1; END;
  END;
  stLo := fiLo; stHi := fiHi;
  st_Grid();
  st_CutCalc();
END;

// ---- build the caches: once per parameter or window change ----
// Every frame -- every sample of a drag -- then draws from these lists
// and evaluates nothing (Lessons Learned s6: cache the density).
st_Cache()
BEGIN
  LOCAL caI, caN, caX, caV, caMx, caMo;
  LOCAL caSd, caLo, caHi, caJ;
  caMo := st_Mom(stId, stPar);
  stMu := caMo(1); stVr := caMo(2);
  IF stDisc == 1 THEN caN := stHi - stLo + 1; ELSE caN := stW + 1; END;
  stPd := MAKELIST(0, K, 1, caN, 1);
  stCd := MAKELIST(0, K, 1, caN, 1);
  stPy := MAKELIST(0, K, 1, caN, 1);
  // v2.7e: the right-tail cache, with a LOCAL counter (caJ), so Home's
  // K is not written (GfxProbe T19: a LOCAL counter works)
  stSv := MAKELIST(0, caJ, 1, caN, 1);
  IF stDisc == 1 THEN
    // a running sum of the pmf, started from the exact CDF just left of
    // the window, gives the CDF at every bar
    caV := st_cdf(stId, stPar, stLo - 1);
    FOR caI FROM 1 TO caN DO
      stPd(caI) := st_pdf(stId, stPar, stLo + caI - 1);
      caV := caV + stPd(caI);
      stCd(caI) := MIN(1, caV);
    END;
    // v2.7e: the mirror image for the right tail: start from the exact
    // P(X > last bar) and add the pmf from the right, bar by bar --
    // P(X > k) = P(X > k + 1) + P(X = k + 1).  Only positive numbers are
    // added, so a tiny tail keeps its digits (1 - F(k) did not).  A
    // WHILE counts down; it does nothing when there is one bar.
    caV := st_sf(stId, stPar, stHi);
    stSv(caN) := MIN(1, caV);
    caI := caN - 1;
    WHILE caI >= 1 DO
      caV := caV + stPd(caI + 1);
      stSv(caI) := MIN(1, caV);
      caI := caI - 1;
    END;
  ELSE
    FOR caI FROM 1 TO caN DO
      caX := stLo + (caI - 1)*(stHi - stLo)/stW;
      stPd(caI) := st_pdf(stId, stPar, caX);
      stCd(caI) := st_cdf(stId, stPar, caX);
      // v2.7e: the right tail.  Where F(x) is at most 1/2, 1 - F(x) is
      // exact to the last digit (nothing near 1 is subtracted); above
      // 1/2 it is st_sf, the family's own right tail.
      stSv(caI) := 1 - stCd(caI);
      IF stCd(caI) > 0.5 THEN stSv(caI) := st_sf(stId, stPar, caX); END;
    END;
  END;
  // The top of the y axis: the tallest value -- ignoring the two
  // columns at each end of a curve, where a density that rises without
  // limit (chi-square df = 1, Gamma k < 1) would flatten everything.
  caMx := 0;
  caLo := 1; caHi := caN;
  IF stDisc == 0 THEN caLo := 3; caHi := caN - 2; END;
  FOR caI FROM caLo TO caHi DO
    IF stPd(caI) > caMx THEN caMx := stPd(caI); END;
  END;
  // an overlay needs room too: check BOTH curves (Lessons Learned s6)
  IF stOvT == 1 AND stId == 7 THEN caMx := MAX(caMx, 0.398942280401433); END;
  IF stOvN == 1 AND (stId == 15 OR stId == 16) THEN
    caMx := MAX(caMx, 0.398942280401433/stVr^0.5);
  END;
  IF caMx <= 0 THEN caMx := 1; END;
  stYmx := caMx*1.12;
  // the pixel heights, clipped at the top of the box
  stClip := 0;
  FOR caI FROM 1 TO caN DO
    caV := stPd(caI)/stYmx;
    IF caV > 1 THEN caV := 1; stClip := 1; END;
    stPy(caI) := stPB - IP(caV*stH + 0.5);
  END;
  // the overlay curve, one pixel height per column
  stNy := MAKELIST(0, K, 1, stW + 1, 1);
  IF stOvT == 1 AND stId == 7 THEN
    FOR caI FROM 1 TO stW + 1 DO
      caX := stLo + (caI - 1)*(stHi - stLo)/stW;
      stNy(caI) := stPB - IP(MIN(NORMALD(0, 1, caX)/stYmx, 1)*stH + 0.5);
    END;
  END;
  IF stOvN == 1 AND (stId == 15 OR stId == 16) THEN
    // N(mean, sd) across the bars: column c sits at k = stLo - 1/2 + ...
    caSd := stVr^0.5;
    FOR caI FROM 1 TO stW + 1 DO
      caX := stLo - 0.5 + (caI - 1)*(stHi - stLo + 1)/stW;
      stNy(caI) := stPB - IP(MIN(NORMALD(stMu, caSd, caX)/stYmx, 1)*stH + 0.5);
    END;
  END;
END;

// A red line for x (or a, or b) with its letter; two pixels wide for
// the bound the arrows move.  Nothing is drawn off the plot.
st_Mark(mkV, mkL, mkB)
BEGIN
  LOCAL mkPx;
  mkPx := st_Px(mkV);
  IF mkPx < stPL OR mkPx > stPR THEN RETURN 0; END;
  LINE_P(G1, mkPx, stPT, mkPx, stPB, stRED);
  IF mkB == 1 THEN LINE_P(G1, mkPx + 1, stPT, mkPx + 1, stPB, stRED); END;
  TEXTOUT_P(mkL, G1, mkPx + 3, stPT + 12, 1, stRED, 12);
  RETURN 1;
END;

// ---- a continuous curve: shading, overlays and the density ----
st_DrawCurve()
BEGIN
  LOCAL dcI, dcC1, dcC2, dcY, dcY2;
  // The cache entries to shade: P(X <= x) left of the red line,
  // P(X > x) right of it, a < X < b between the two lines.  Shading
  // goes pixel column by pixel column (Lessons Learned s6: no gaps).
  IF stTl == 1 THEN dcC1 := 1; dcC2 := st_Px(stXv) - stPL + 1; END;
  IF stTl == 2 THEN dcC1 := st_Px(stXv) - stPL + 2; dcC2 := stW + 1; END;
  IF stTl == 3 THEN dcC1 := st_Px(stAv) - stPL + 1; dcC2 := st_Px(stBv) - stPL + 1; END;
  dcC1 := MAX(1, dcC1); dcC2 := MIN(stW + 1, dcC2);
  FOR dcI FROM dcC1 TO dcC2 DO
    dcY := stPy(dcI);
    IF dcY < stPB - 1 THEN LINE_P(G1, stPL + dcI - 1, stPB - 1, stPL + dcI - 1, dcY, stSHD); END;
  END;
  // N(0,1) behind Student t, in pale grey: t has heavier tails
  IF stOvT == 1 AND stId == 7 THEN
    FOR dcI FROM 2 TO stW + 1 DO
      LINE_P(G1, stPL + dcI - 2, stNy(dcI - 1), stPL + dcI - 1, stNy(dcI), stPALE);
    END;
  END;
  // the density itself, from the cache
  FOR dcI FROM 2 TO stW + 1 DO
    LINE_P(G1, stPL + dcI - 2, stPy(dcI - 1), stPL + dcI - 1, stPy(dcI), stBLU);
  END;
  // the CDF as a magenta S-curve, 0 at the bottom of the box, 1 at the top
  IF stOvC == 1 THEN
    dcY := stPB - IP(stCd(1)*stH + 0.5);
    FOR dcI FROM 2 TO stW + 1 DO
      dcY2 := stPB - IP(stCd(dcI)*stH + 0.5);
      LINE_P(G1, stPL + dcI - 2, dcY, stPL + dcI - 1, dcY2, stMAG);
      dcY := dcY2;
    END;
  END;
END;

// ---- discrete bars: counted bars filled, the chosen k red ----
// The chosen bar is filled red when it is counted in the probability
// and only outlined in red when it is not -- so P(X > 3) shows bar 3
// red but EMPTY: "more than 3" starts at 4.
st_DrawBars()
BEGIN
  LOCAL dbI, dbK, dbX1, dbX2, dbY, dbIn;
  LOCAL dbBw, dbGp, dbCh, dbE, dbF, dbY0;
  LOCAL dbN, dbL;
  dbN := stHi - stLo + 1;
  dbBw := stW/dbN;
  dbGp := 0;
  IF dbBw >= 4 THEN dbGp := 1; END;        // a gap between wide bars
  FOR dbI FROM 1 TO dbN DO
    dbK := stLo + dbI - 1;
    dbX1 := stPL + IP((dbI - 1)*dbBw) + dbGp;
    dbX2 := stPL + IP(dbI*dbBw) - 1 - dbGp;
    IF dbX2 < dbX1 THEN dbX2 := dbX1; END;
    // is this bar counted?  is it a chosen k (x, or a or b)?
    dbIn := 0; dbCh := 0; dbL := "k";
    IF stTl == 1 AND dbK <= stXv THEN dbIn := 1; END;
    IF stTl == 2 AND dbK > stXv THEN dbIn := 1; END;
    IF stTl == 3 AND dbK >= stAv AND dbK <= stBv THEN dbIn := 1; END;
    IF stTl <> 3 AND dbK == stXv THEN dbCh := 1; END;
    IF stTl == 3 AND dbK == stAv THEN dbCh := 1; dbL := "a"; END;
    IF stTl == 3 AND dbK == stBv THEN dbCh := 1; dbL := "b"; END;
    dbE := stBLU; dbF := #FFFFFFh;
    IF dbIn == 1 THEN dbF := stSHD; END;
    IF dbCh == 1 THEN
      dbE := stRED;
      IF dbIn == 1 THEN dbF := stRED; END;
    END;
    dbY := stPy(dbI);
    IF dbY < stPB - 1 THEN RECT_P(G1, dbX1, dbY, dbX2, stPB - 1, dbE, dbF); END;
    IF dbCh == 1 THEN
      // the letter over the chosen bar, and a red tick under the axis
      // in case the bar is too short to see
      TEXTOUT_P(dbL, G1, IP((dbX1 + dbX2)/2) - 2, MAX(stPT, dbY - 11), 1, stRED, 10);
      LINE_P(G1, IP((dbX1 + dbX2)/2), stPB + 1, IP((dbX1 + dbX2)/2), stPB + 3, stRED);
    END;
  END;
  // the normal approximation, magenta, over the bars
  IF stOvN == 1 AND (stId == 15 OR stId == 16) THEN
    FOR dbI FROM 2 TO stW + 1 DO
      LINE_P(G1, stPL + dbI - 2, stNy(dbI - 1), stPL + dbI - 1, stNy(dbI), stMAG);
    END;
  END;
  // the CDF of a discrete distribution is a staircase
  IF stOvC == 1 THEN
    dbY0 := stPB - IP((stCd(1) - stPd(1))*stH + 0.5);
    FOR dbI FROM 1 TO dbN DO
      dbX1 := stPL + IP((dbI - 1)*dbBw);
      dbX2 := stPL + IP(dbI*dbBw);
      dbY := stPB - IP(stCd(dbI)*stH + 0.5);
      LINE_P(G1, dbX1, dbY0, dbX1, dbY, stMAG);
      LINE_P(G1, dbX1, dbY, dbX2, dbY, stMAG);
      dbY0 := dbY;
    END;
  END;
END;

// ---- the orange panel under the plot ----
st_Panel()
BEGIN
  LOCAL pnS, pnP, pnSd, pnN, pnC;
  LOCAL pnA, pnZ, pnPoor, pnX;
  RECT_P(G1, 0, 152, 319, 217, stORN, stLORN);
  // line 1: the parameters
  TEXTOUT_P(st_ParStr(stId, stPar), G1, 6, 155, 1, #000000h, 308);
  // line 2: mean and sd, green like the dotted mean line
  IF TYPE(stVr) == 0 THEN pnSd := st_fmtx(MAX(0, stVr)^0.5); ELSE pnSd := "undef"; END;
  // v2.7d: st_Eq -- the relation without "=" ("P(X <= -6) < 0.000001")
  TEXTOUT_P("mean" + st_Eq(st_fmtx(stMu)) + ", sd" + st_Eq(pnSd), G1, 6, 167, 1, stGRN, 308);
  // line 3: the height of the curve -- a density is NOT a probability
  // (it can be bigger than 1); a bar's height IS one
  IF stDisc == 1 THEN
    IF stTl == 3 THEN
      // v2.7d: st_Eq -- the relation without "=" ("P(X <= -6) < 0.000001")
      pnS := "P(X = " + st_fmtx(stAv) + ")" + st_Eq(st_fmtp(st_PdfAt(stAv), 4)) + ", P(X = " + st_fmtx(stBv) + ")" + st_Eq(st_fmtp(st_PdfAt(stBv), 4));
    ELSE
      // v2.7d: st_Eq -- the relation without "=" ("P(X <= -6) < 0.000001")
      pnS := "P(X = " + st_fmtx(stXv) + ")" + st_Eq(st_fmtp(st_PdfAt(stXv), 6)) + ": the red bar";
    END;
  ELSE
    IF stTl == 3 THEN
      // v2.7d: st_Eq -- the relation without "=" ("P(X <= -6) < 0.000001")
      pnS := "f(a)" + st_Eq(st_fmtp(st_PdfAt(stAv), 4)) + ", f(b)" + st_Eq(st_fmtp(st_PdfAt(stBv), 4)) + ": densities";
    ELSE
      // v2.7d: st_Eq -- the relation without "=" ("P(X <= -6) < 0.000001")
      pnS := "f(x)" + st_Eq(st_fmtp(st_PdfAt(stXv), 4)) + ": a density, not a probability";
    END;
  END;
  TEXTOUT_P(pnS, G1, 6, 179, 1, #000000h, 308);
  // line 4: THE answer, magenta, in the bigger font
  pnP := st_Prob();
  // v2.7d: st_Eq -- the relation without "=" ("P(X <= -6) < 0.000001")
  // v2.7e: x, a and b through st_XStr -- a value, never "< 0.0001" (B7)
  IF stTl == 1 THEN pnS := "P(X <= " + st_XStr(stXv, st_fmtx(stXv)) + ")" + st_Eq(st_fmtp(pnP, 6)); END;
  IF stTl == 2 THEN
    IF stDisc == 1 THEN
      // v2.7d: st_Eq -- the relation without "=" ("P(X <= -6) < 0.000001")
      pnS := "P(X > " + st_fmtx(stXv) + ") = P(X >= " + st_fmtx(stXv + 1) + ")" + st_Eq(st_fmtp(pnP, 6));
    ELSE
      // v2.7d: st_Eq -- the relation without "=" ("P(X <= -6) < 0.000001")
      pnS := "P(X > " + st_XStr(stXv, st_fmtx(stXv)) + ")" + st_Eq(st_fmtp(pnP, 6));
    END;
  END;
  IF stTl == 3 THEN
    IF stDisc == 1 THEN
      // v2.7d: st_Eq -- the relation without "=" ("P(X <= -6) < 0.000001")
      pnS := "P(" + st_fmtx(stAv) + " <= X <= " + st_fmtx(stBv) + ")" + st_Eq(st_fmtp(pnP, 6));
    ELSE
      // v2.7d: st_Eq -- the relation without "=" ("P(X <= -6) < 0.000001")
      pnS := "P(" + st_XStr(stAv, st_fmtx(stAv)) + " < X < " + st_XStr(stBv, st_fmtx(stBv)) + ")" + st_Eq(st_fmtp(pnP, 6));
    END;
  END;
  TEXTOUT_P(pnS, G1, 6, 191, 2, stMAG, 308);
  stLast := pnS;                       // returned to Home on exit
  // line 5: one note, the most useful one
  pnN := ""; pnC := stFOOT;
  IF stTl == 3 AND stPh == 1 THEN pnN := "tap the plot to set a (the lower bound)"; pnC := stRED; END;
  IF stTl == 3 AND stPh == 2 THEN pnN := "tap the plot to set b (the upper bound)"; pnC := stRED; END;
  IF DIM(pnN) == 0 AND stOvN == 1 AND (stId == 15 OR stId == 16) THEN
    // the area under the normal curve, continuity corrected (+/- 0.5)
    pnZ := stVr^0.5;
    IF stTl == 1 THEN pnA := NORMALD_CDF(stMu, pnZ, stXv + 0.5); END;
    // v2.7d: the right tail as the LOWER tail at 2 mean - (x + 0.5), the same
    // number by symmetry: 1 - NORMALD_CDF is exactly 0 for a far tail
    IF stTl == 2 THEN pnA := NORMALD_CDF(stMu, pnZ, 2*stMu - (stXv + 0.5)); END;
    IF stTl == 3 THEN pnA := NORMALD_CDF(stMu, pnZ, stBv + 0.5) - NORMALD_CDF(stMu, pnZ, stAv - 0.5); END;
    // the rule of thumb: np and n(1-p) at least 10; lambda at least 10
    pnPoor := 0;
    IF stId == 15 AND (stPar(1)*stPar(2) < 10 OR stPar(1)*(1 - stPar(2)) < 10) THEN pnPoor := 1; END;
    IF stId == 16 AND stPar(1) < 10 THEN pnPoor := 1; END;
    IF pnPoor == 1 THEN
      IF stId == 15 THEN pnX := "np or n(1-p) < 10"; ELSE pnX := "λ < 10"; END;
      // v2.7e: 46 characters, not 58 -- the old line ran past the
      // screen's edge at font 1 (B2); " = " through st_Eq as below
      pnN := "approx." + st_Eq(st_fmtp(pnA, 4)) + ": poor here, " + pnX;
      pnC := stRED;
    ELSE
      // v2.7d: st_Eq -- the relation without "=" ("P(X <= -6) < 0.000001")
      pnN := "normal approx." + st_Eq(st_fmtp(pnA, 4)) + " (continuity corr.)";
      pnC := stMAG;
    END;
  END;
  IF DIM(pnN) == 0 AND stOvT == 1 AND stId == 7 THEN
    pnN := "grey = N(0,1): t has heavier tails"; pnC := stGRY;
  END;
  IF DIM(pnN) == 0 AND stClip == 1 THEN pnN := "curve clipped at the top of the plot"; pnC := stGRY; END;
  IF DIM(pnN) == 0 AND stZm == 1 THEN pnN := "zoomed, 0 = refit the window"; END;
  IF DIM(pnN) == 0 THEN pnN := "Ent = options, +/- = zoom, 0 = refit, Esc = menu"; END;
  TEXTOUT_P(pnN, G1, 6, 205, 1, pnC, 308);
END;

// ---- compose one frame in G1, blit it once, then the soft keys ----
st_Draw()
BEGIN
  LOCAL dwPx, dwV, dwN, dwS, dwD, dwY;
  LOCAL dwYs, dwF;
  DIMGROB_P(G1, 320, 240, #FFFFFFh);
  // title: the name in blue, the typical use in grey
  dwS := st_Name(stId);
  TEXTOUT_P(dwS, G1, 4, 1, 2, stBLU, 200);
  TEXTOUT_P(st_Use(stId), G1, 12 + 8*DIM(dwS), 4, 1, stGRY, 300 - 8*DIM(dwS));
  // the x grid and tick labels at nice values (below the axis)
  dwD := MAX(0, -FLOOR(LOG(stTk)));
  dwN := CEILING(stLo/stTk - 0.000001);
  dwV := dwN*stTk;
  WHILE dwV <= stHi + stTk/1000 DO
    dwPx := st_Px(dwV);
    IF dwPx >= stPL AND dwPx <= stPR THEN
      LINE_P(G1, dwPx, stPT, dwPx, stPB, stGRID);
      LINE_P(G1, dwPx, stPB, dwPx, stPB + 3, #000000h);
      IF ABS(dwV) < stTk/1000 THEN dwV := 0; END;   // no "-0.0"
      dwS := STRING(dwV, 2, dwD);
      TEXTOUT_P(dwS, G1, dwPx - 3*DIM(dwS), stPB + 4, 1, stGRY, 60);
    END;
    dwN := dwN + 1;
    dwV := dwN*stTk;
  END;
  // the y grid and labels (a density's scale, or a pmf's)
  dwYs := st_Nice(stYmx/4);
  dwD := MAX(0, -FLOOR(LOG(dwYs)));
  dwV := dwYs;
  WHILE dwV < stYmx DO
    dwY := stPB - IP(dwV/stYmx*stH + 0.5);
    LINE_P(G1, stPL + 1, dwY, stPR, dwY, stGRID);
    LINE_P(G1, stPL - 3, dwY, stPL, dwY, #000000h);
    IF dwD <= 4 THEN TEXTOUT_P(STRING(dwV, 2, dwD), G1, 1, dwY - 4, 1, stGRY, stPL - 4); END;
    dwV := dwV + dwYs;
  END;
  // the distribution, with its shading
  IF stDisc == 1 THEN st_DrawBars(); ELSE st_DrawCurve(); END;
  // the mean: a green dotted line (PIXON_P every 4th pixel -- LINE_P
  // has no dash style), only when the mean exists
  IF TYPE(stMu) == 0 THEN
    dwPx := st_Px(stMu);
    IF dwPx > stPL AND dwPx < stPR THEN
      FOR dwY FROM stPT TO stPB STEP 4 DO PIXON_P(G1, dwPx, dwY, stGRN); END;
    END;
  END;
  // x, or a and b: red lines (a discrete k is the red bar itself)
  IF stDisc == 0 THEN
    IF stTl == 3 THEN
      st_Mark(stAv, "a", stAct == 1);
      st_Mark(stBv, "b", stAct == 2);
    ELSE
      st_Mark(stXv, "x", 0);
    END;
  END;
  // the CDF's own axis on the right, 0 to 1, and F(x) read across to it
  IF stOvC == 1 THEN
    LINE_P(G1, stPR, stPT, stPR, stPB, stMAG);
    TEXTOUT_P("1", G1, stPR - 8, stPT + 1, 1, stMAG, 10);
    TEXTOUT_P(".5", G1, stPR - 12, stPB - IP(stH/2) - 5, 1, stMAG, 14);
    IF stTl <> 3 THEN
      dwY := stPB - IP(st_CdfAt(stXv)*stH + 0.5);
      dwPx := MAX(stPL, MIN(stPR, st_Px(stXv)));
      WHILE dwPx < stPR DO
        PIXON_P(G1, dwPx, dwY, stMAG);
        dwPx := dwPx + 4;
      END;
    END;
  END;
  // the axes
  LINE_P(G1, stPL, stPT, stPL, stPB, #000000h);
  LINE_P(G1, stPL, stPB, stPR, stPB, #000000h);
  // say when part of the distribution lies off the plot (a cut tail)
  IF stCut(2) >= 0.0015 THEN
    dwF := STRING(100*stCut(2), 2, 1) + "% beyond >>";
    TEXTOUT_P(dwF, G1, stPR - 4 - 6*DIM(dwF), stPT + 2, 1, stGRY, 100);
  END;
  IF stCut(1) >= 0.0015 THEN
    TEXTOUT_P("<< " + STRING(100*stCut(1), 2, 1) + "% beyond", G1, stPL + 3, stPT + 2, 1, stGRY, 100);
  END;
  st_Panel();
  BLIT_P(G0, 0, 0, G1);                  // one blit, then the soft keys
  // Three pages of soft keys (the house "More" convention, 16 Oct 2026):
  // "More" is always slot 6 and turns the page.  Page 1 = the tail and
  // the zoom; page 2 = everything the Enter menu offers, for a student
  // who is only using the touch screen; page 3 (v2.7) = the Formula
  // page and Help, with Menu in the same slot as on page 2.
  IF stSKP == 1 THEN
    IF stDisc == 1 THEN
      DRAWMENU({"P<=k", "P>k", "a<X<b", "Zoom−", "Zoom+", "More"});
    ELSE
      DRAWMENU({"P<=x", "P>x", "a<X<b", "Zoom−", "Zoom+", "More"});
    END;
  ELSE
    IF stSKP == 2 THEN
      IF stId == 15 OR stId == 16 THEN
        DRAWMENU({"Edit", "Numbers", "CDF", "Approx", "Menu", "More"});
      ELSE
        IF stId == 7 THEN
          DRAWMENU({"Edit", "Numbers", "CDF", "N(0,1)", "Menu", "More"});
        ELSE
          DRAWMENU({"Edit", "Numbers", "CDF", "Refit", "Menu", "More"});
        END;
      END;
    ELSE
      DRAWMENU({"Formula", "Help", "", "", "Menu", "More"});
    END;
  END;
END;

// Put value smV where the current marker is: x, or the bound in hand.
st_SetMark(smV)
BEGIN
  IF stTl == 3 THEN
    IF stAct == 1 THEN stAv := smV; ELSE stBv := smV; END;
  ELSE
    stXv := smV;
  END;
END;

// Keep a <= b; if they cross, the bound in hand swaps letters with them.
st_Order()
BEGIN
  LOCAL orT;
  IF stTl == 3 AND stAv > stBv THEN
    orT := stAv; stAv := stBv; stBv := orT;
    stAct := 3 - stAct;
  END;
END;

// The soft keys P<=x (1), P>x (2) and a<X<b (3) change the tail in
// place.  a<X<b starts the two-tap set-up; a and b begin at the
// default rule unless both already sit on the plot in order.
st_SetTail(teT)
BEGIN
  LOCAL teD;
  IF teT == 3 THEN
    IF stAv < stLo OR stAv > stHi OR stBv < stLo OR stBv > stHi OR stAv >= stBv THEN
      teD := st_DefAB(stId, stPar);
      stAv := teD(1); stBv := teD(2);
      st_Clamp();
    END;
    stPh := 1; stAct := 1;
  ELSE
    stPh := 0;
  END;
  stTl := teT;
END;

// Up (2) and Down (12): the upper or lower tail; in a<X<b, which bound
// the arrows move (Up = b, Down = a).
st_UpDn(udK)
BEGIN
  IF stTl == 3 THEN
    stPh := 0;
    IF udK == 2 THEN stAct := 2; ELSE stAct := 1; END;
  ELSE
    IF udK == 2 THEN st_SetTail(2); ELSE st_SetTail(1); END;
  END;
END;

// When x walks off the plot, slide the window after it.
st_Pan(paV)
BEGIN
  LOCAL paW, paS;
  IF paV >= stLo AND paV <= stHi THEN RETURN 0; END;
  paW := stHi - stLo;
  IF stDisc == 1 THEN
    IF paV < stLo THEN stLo := paV - IP(paW/4); ELSE stLo := paV - IP(3*paW/4); END;
    paS := st_Sup(stId, stPar);
    stLo := MAX(paS(1), stLo);
  ELSE
    IF paV < stLo THEN stLo := paV - paW/4; ELSE stLo := paV - 3*paW/4; END;
  END;
  stHi := stLo + paW;
  stZm := 1;
  st_Cache(); st_CutCalc();
  RETURN 1;
END;

// Left (-1) and Right (+1): x (or the bound in hand) to the next grid
// line.  Each press lands ON the grid (NormExplorer's idiom), so a
// typed 1.96 steps to 2.0, then 2.1 -- never 2.06, 2.16.
st_Step(spD)
BEGIN
  LOCAL spV, spS;
  IF stTl == 3 THEN
    IF stAct == 1 THEN spV := stAv; ELSE spV := stBv; END;
  ELSE
    spV := stXv;
  END;
  IF stDisc == 1 THEN
    spS := st_Sup(stId, stPar);
    spV := MIN(spS(2), MAX(spS(1), spV + spD));
  ELSE
    IF spD > 0 THEN
      spV := (FLOOR(spV/stGs + 0.000001) + 1)*stGs;
    ELSE
      spV := (CEILING(spV/stGs - 0.000001) - 1)*stGs;
    END;
    spV := ROUND(spV, 10);
  END;
  st_SetMark(spV);
  st_Order();
  st_Pan(spV);
END;

// Zoom about x (or the middle of a..b): zmF = 0.5 halves the window
// (Zoom+), 2 doubles it (Zoom-).  The caches are rebuilt.
st_Zoom(zmF)
BEGIN
  LOCAL zmC, zmH, zmS, zmN;
  IF stTl == 3 THEN zmC := (stAv + stBv)/2; ELSE zmC := stXv; END;
  IF stDisc == 1 THEN
    zmS := st_Sup(stId, stPar);
    zmN := MAX(5, MIN(251, ROUND((stHi - stLo + 1)*zmF, 0)));   // bars on screen
    stLo := MAX(zmS(1), IP(zmC - zmN/2));
    stHi := stLo + zmN - 1;
    IF stHi > zmS(2) THEN stHi := zmS(2); stLo := MAX(zmS(1), stHi - zmN + 1); END;
  ELSE
    zmH := (stHi - stLo)*zmF/2;
    IF zmH < 0.000000001*(1 + ABS(zmC)) OR zmH > 1E12 THEN RETURN 0; END;
    stLo := zmC - zmH; stHi := zmC + zmH;
  END;
  st_Grid();
  stZm := 1;
  st_Cache(); st_CutCalc();
  RETURN 1;
END;

// Snap a value onto the x grid (discrete values are whole already).
st_Snap(snV)
BEGIN
  IF stDisc == 1 THEN RETURN snV; END;
  RETURN ROUND(ROUND(snV/stGs, 0)*stGs, 10);
END;

// ---- a finger on the plot ----
// The first sample is the press.  In a<X<b it decides which line moves:
// the pending one during the two-tap set-up, else the nearer one.
// While the finger stays down, every new pixel column moves the marker
// and redraws from the caches (stDrg = 1: nothing is evaluated).  On
// lift the value snaps to the grid and the readout is recomputed
// exactly.  The window stays frozen during the drag (Lessons Learned
// s6).  Esc, or a reading off the screen, also ends the drag.
st_Drag(dgP0)
BEGIN
  LOCAL dgM, dgM1, dgK, dgV, dgLast, dgGo;
  LOCAL dgPx;
  dgPx := dgP0;
  stDrg := 1;
  IF stTl == 3 THEN
    IF stPh == 1 THEN stAct := 1; END;
    IF stPh == 2 THEN stAct := 2; END;
    IF stPh == 0 THEN
      IF ABS(dgPx - st_Px(stAv)) <= ABS(dgPx - st_Px(stBv)) THEN stAct := 1; ELSE stAct := 2; END;
    END;
  END;
  dgPx := MIN(stPR, MAX(stPL, dgPx));
  dgLast := -999; dgGo := 1; dgV := 0;
  WHILE dgGo == 1 DO
    IF dgPx <> dgLast THEN
      dgV := st_Val(dgPx);
      st_SetMark(dgV);
      // the first a<X<b tap places a AND b together (P = 0); the second
      // stretches the interval.  The bounds stay in order on every sample,
      // so a crossing never shows a negative probability.
      IF stTl == 3 AND stPh == 1 THEN stBv := dgV; END;
      st_Order();
      st_Draw();
      dgLast := dgPx;
    END;
    dgM := MOUSE;
    dgM1 := dgM(1);
    dgK := GETKEY;
    IF SIZE(dgM1) == 0 OR dgK == 4 THEN
      dgGo := 0;                                  // the finger lifted
    ELSE
      IF dgM1(1) >= 0 AND dgM1(1) <= 319 AND dgM1(2) >= 0 AND dgM1(2) <= 239 THEN
        dgPx := MIN(stPR, MAX(stPL, dgM1(1)));
      ELSE
        dgGo := 0;                                // a phantom reading
      END;
    END;
  END;
  stDrg := 0;
  st_SetMark(st_Snap(dgV));
  IF stTl == 3 AND stPh == 1 THEN stBv := stAv; END;
  IF stTl == 3 THEN
    IF stPh == 1 THEN
      stPh := 2; stAct := 2;                      // now tap for b
    ELSE
      IF stPh == 2 THEN stPh := 0; END;           // both set
    END;
    st_Order();
  END;
  st_Drain();
END;

// ---- the numbers page: v2.5's result text, built from the plot ----
// It follows whatever is on the plot now -- a dragged x, a new tail,
// edited parameters -- so nothing on the v2.5 pages is lost.
st_Numbers()
BEGIN
  LOCAL nuL, nuM, nuSd, nuF, nuFa, nuFb;
  LOCAL nuX, nuP, nuS;
  nuL := {"Typical use: " + st_Use(stId), st_ParStr(stId, stPar)};
  nuM := st_fmtv(stMu, 4);
  IF TYPE(stVr) == 0 THEN
    nuSd := st_fmtp(MAX(0, stVr)^0.5, 4);
    // v2.7d: st_Eq -- the relation without "=" ("P(X <= -6) < 0.000001")
    nuL := CONCAT(nuL, {"mean" + st_Eq(nuM) + ",  sd" + st_Eq(nuSd), "variance" + st_Eq(st_fmtp(stVr, 4))});
  ELSE
    // v2.7d: st_Eq -- the relation without "=" ("P(X <= -6) < 0.000001")
    nuL := CONCAT(nuL, {"mean" + st_Eq(nuM), "variance = " + stVr});
  END;
  nuL := CONCAT(nuL, {""});
  nuP := st_Prob();
  IF stTl == 3 THEN
    IF stDisc == 1 THEN
      nuFa := st_cdf(stId, stPar, stAv - 1); nuFb := st_cdf(stId, stPar, stBv);
      // v2.7d: st_Eq -- the relation without "=" ("P(X <= -6) < 0.000001")
      nuL := CONCAT(nuL, {"CDF(b) = P(X <= " + STRING(stBv, 2, 0) + ")" + st_Eq(st_fmtp(nuFb, 6)),
        "CDF(a-1) = P(X <= " + STRING(stAv - 1, 2, 0) + ")" + st_Eq(st_fmtp(nuFa, 6)),
        "P = CDF(b) - CDF(a-1)", "",
        "P(" + STRING(stAv, 2, 0) + " <= X <= " + STRING(stBv, 2, 0) + ")" + st_Eq(st_fmtp(nuP, 6))});
    ELSE
      nuFa := st_cdf(stId, stPar, stAv); nuFb := st_cdf(stId, stPar, stBv);
      // v2.7d: st_Eq -- the relation without "=" ("P(X <= -6) < 0.000001")
      // v2.7e: a and b through st_XStr -- a value, never "< 0.0001" (B7)
      nuL := CONCAT(nuL, {"CDF(b)" + st_Eq(st_fmtp(nuFb, 6)) + " at b" + st_Eq(st_XStr(stBv, st_fmtp(stBv, 4))),
        "CDF(a)" + st_Eq(st_fmtp(nuFa, 6)) + " at a" + st_Eq(st_XStr(stAv, st_fmtp(stAv, 4))),
        "P = CDF(b) - CDF(a)", "",
        "P(" + st_XStr(stAv, st_fmtp(stAv, 4)) + " < X < " + st_XStr(stBv, st_fmtp(stBv, 4)) + ")" + st_Eq(st_fmtp(nuP, 6))});
    END;
  ELSE
    nuF := st_cdf(stId, stPar, stXv);
    IF stDisc == 1 THEN
      nuX := STRING(stXv, 2, 0);
      nuL := CONCAT(nuL, {"At k = " + nuX + ":",
        "  PMF: " + st_fmtp(st_pdf(stId, stPar, stXv), 6),
        "  CDF: " + st_fmtp(nuF, 6), ""});
      IF stTl == 1 THEN
        // v2.7d: st_Eq -- the relation without "=" ("P(X <= -6) < 0.000001")
        nuS := "P(X <= " + nuX + ")" + st_Eq(st_fmtp(nuP, 6));
      ELSE
        // v2.7d: st_Eq -- the relation without "=" ("P(X <= -6) < 0.000001")
        nuS := "P(X > " + nuX + ") = P(X >= " + STRING(stXv + 1, 2, 0) + ")" + st_Eq(st_fmtp(nuP, 6));
      END;
    ELSE
      // v2.7e: x through st_XStr -- a value, never "< 0.0001" (B7)
      nuX := st_XStr(stXv, st_fmtp(stXv, 4));
      // v2.7d: st_Eq -- the relation without "=" ("P(X <= -6) < 0.000001")
      nuL := CONCAT(nuL, {"At x" + st_Eq(nuX) + ":",
        "  PDF: " + st_fmtp(st_pdf(stId, stPar, stXv), 6) + " (a density)",
        "  CDF: " + st_fmtp(nuF, 6), ""});
      IF stTl == 1 THEN
        // v2.7d: st_Eq -- the relation without "=" ("P(X <= -6) < 0.000001")
        nuS := "P(X <= " + nuX + ")" + st_Eq(st_fmtp(nuP, 6));
      ELSE
        // v2.7d: st_Eq -- the relation without "=" ("P(X <= -6) < 0.000001")
        nuS := "P(X > " + nuX + ")" + st_Eq(st_fmtp(nuP, 6));
      END;
    END;
    nuL := CONCAT(nuL, {nuS});
  END;
  st_Page(st_Name(stId) + " distribution", nuL);
END;

// Edit: the parameter form again, with the plot's values in it.  In
// a<X<b, a and b go back to their default rule (the old ones may not
// fit the new parameters).
st_Edit()
BEGIN
  LOCAL edS;
  edS := st_Ask(stId, stTl, st_Cur(), 0);
  IF TYPE(edS) == 6 THEN
    stPar := {edS(1), edS(2), edS(3)};
    stXv := edS(4);
    IF stTl == 3 THEN stAv := edS(5); stBv := edS(6); stPh := 0; END;
    st_Clamp();
    st_Fit(1); st_Cache();
  END;
  RETURN 0;
END;

// ============================================================
// 7b. THE FORMULA PAGE (v2.7): the CAS does the calculus
// "Where does E[X] come from?", answered for the distribution on the
// plot with the student's own numbers.  Continuous: E[X], E[X^2]
// and Var(X) as integrals (and the CDF for the Exponential, Weibull
// and Pareto); discrete: the probability on the plot as an exact
// sum.  Built when the page opens (st_FmBuild), kept until the plot
// changes (st_FmSame), shown by st_FmPage.  Each CAS answer passes
// st_FmAsk's check against the program's own number first.
// ============================================================

// The names the parameters go by on the page, for a symbolic answer.
st_FmNm(nmI)
BEGIN
  CASE
    IF nmI == 1 OR nmI == 5 THEN RETURN {"μ", "σ", ""}; END;
    IF nmI == 2 THEN RETURN {"a", "b", ""}; END;
    IF nmI == 3 OR nmI == 16 THEN RETURN {"λ", "", ""}; END;
    IF nmI == 4 THEN RETURN {"k", "λ", ""}; END;
    IF nmI == 6 THEN RETURN {"μ", "b", ""}; END;
    IF nmI == 7 OR nmI == 8 THEN RETURN {"df", "", ""}; END;
    IF nmI == 9 THEN RETURN {"df1", "df2", ""}; END;
    IF nmI == 10 THEN RETURN {"C", "K", "S"}; END;
    IF nmI == 11 THEN RETURN {"μ", "β", ""}; END;
    IF nmI == 12 THEN RETURN {"x_m", "α", ""}; END;
    IF nmI == 13 THEN RETURN {"σ", "", ""}; END;
    IF nmI == 14 THEN RETURN {"μ", "s", ""}; END;
    IF nmI == 15 THEN RETURN {"n", "p", ""}; END;
    IF nmI == 17 THEN RETURN {"p", "", ""}; END;
    IF nmI == 18 THEN RETURN {"min", "mode", "max"}; END;
    IF nmI == 19 THEN RETURN {"x0", "scale", ""}; END;
    IF nmI == 20 THEN RETURN {"α", "s", "m"}; END;
    IF nmI == 21 THEN RETURN {"eta", "b", ""}; END;
    IF nmI == 22 THEN RETURN {"a", "", ""}; END;
    IF nmI == 23 THEN RETURN {"k", "λ", ""}; END;
    IF nmI == 24 THEN RETURN {"k", "theta", ""}; END;
    IF nmI == 25 THEN RETURN {"α", "β", ""}; END;
    IF nmI == 26 THEN RETURN {"r", "p", ""}; END;
    IF nmI == 27 THEN RETURN {"N", "K", "n"}; END;
  END;
  RETURN {"", "", ""};
END;

// Which parameters must be positive (scales and shapes): a symbolic
// one is handed to the CAS with that assumption ("| stqa>0").
st_FmPos(poI)
BEGIN
  CASE
    IF poI == 1 OR poI == 5 OR poI == 6 OR poI == 11 OR poI == 14 OR poI == 19 THEN RETURN {0, 1, 0}; END;
    IF poI == 3 OR poI == 7 OR poI == 8 OR poI == 13 OR poI == 16 OR poI == 22 THEN RETURN {1, 0, 0}; END;
    IF poI == 4 OR poI == 9 OR poI == 12 OR poI == 21 OR poI == 23 OR poI == 24 THEN RETURN {1, 1, 0}; END;
    IF poI == 25 THEN RETURN {1, 1, 0}; END;
    IF poI == 10 THEN RETURN {1, 1, 1}; END;
    IF poI == 20 THEN RETURN {1, 1, 0}; END;
    IF poI == 15 OR poI == 26 THEN RETURN {0, 1, 0}; END;
    IF poI == 17 THEN RETURN {1, 0, 0}; END;
  END;
  RETURN {0, 0, 0};
END;

// The density (or pmf) in words, as the page's second line.
st_FmF(ffI)
BEGIN
  CASE
    IF ffI == 1 THEN RETURN "f(x) = e^(-(x-μ)^2/(2σ^2)) / (σ sqrt(2 pi))"; END;
    IF ffI == 2 THEN RETURN "f(x) = 1/(b - a) for a <= x <= b"; END;
    IF ffI == 3 THEN RETURN "f(x) = λ e^(-λx) for x >= 0"; END;
    IF ffI == 4 THEN RETURN "f(x) = (k/λ) (x/λ)^(k-1) e^(-(x/λ)^k), x >= 0"; END;
    IF ffI == 5 THEN RETURN "f(x) = e^(-(ln x - μ)^2/(2σ^2)) / (x σ sqrt(2 pi))"; END;
    IF ffI == 6 THEN RETURN "f(x) = e^(-|x - μ|/b) / (2b)"; END;
    IF ffI == 7 THEN RETURN "f(x) = Gamma((df+1)/2) / (sqrt(df pi) Gamma(df/2)) (1 + x^2/df)^(-(df+1)/2)"; END;
    IF ffI == 8 THEN RETURN "f(x) = x^(df/2-1) e^(-x/2) / (2^(df/2) Gamma(df/2))"; END;
    IF ffI == 9 THEN RETURN "f(x) = sqrt((d1 x)^d1 d2^d2 / (d1 x + d2)^(d1+d2)) / (x B(d1/2, d2/2))"; END;
    IF ffI == 10 THEN RETURN "f(x) = K C (x/S)^C / (x (1 + (x/S)^C)^(K+1)), x > 0"; END;
    IF ffI == 11 THEN RETURN "f(x) = e^(-z - e^(-z)) / β, z = (x - μ)/β"; END;
    IF ffI == 12 THEN RETURN "f(x) = α x_m^α / x^(α+1) for x >= x_m"; END;
    IF ffI == 13 THEN RETURN "f(x) = (x/σ^2) e^(-x^2/(2σ^2)), x >= 0"; END;
    IF ffI == 14 THEN RETURN "f(x) = e^(-z) / (s (1 + e^(-z))^2), z = (x - μ)/s"; END;
    IF ffI == 15 THEN RETURN "P(X = k) = C(n, k) p^k (1 - p)^(n-k)"; END;
    IF ffI == 16 THEN RETURN "P(X = k) = e^(-λ) λ^k / k!"; END;
    IF ffI == 17 THEN RETURN "P(X = k) = p (1 - p)^(k-1), k = 1, 2, 3 ..."; END;
    IF ffI == 18 THEN RETURN "f(x) = 2(x - a)/((b - a)(c - a)) up to the mode c, then 2(b - x)/((b - a)(b - c))"; END;
    IF ffI == 19 THEN RETURN "f(x) = 1 / (pi scale (1 + z^2)), z = (x - x0)/scale"; END;
    IF ffI == 20 THEN RETURN "f(x) = (α/s) z^(-1-α) e^(-z^(-α)), z = (x - m)/s"; END;
    IF ffI == 21 THEN RETURN "f(x) = b eta e^(eta + b x - eta e^(b x)), x >= 0"; END;
    IF ffI == 22 THEN RETURN "f(x) = sqrt(2/pi) x^2 e^(-x^2/(2a^2)) / a^3"; END;
    IF ffI == 23 THEN RETURN "f(x) = λ^k x^(k-1) e^(-λx) / (k-1)!"; END;
    IF ffI == 24 THEN RETURN "f(x) = x^(k-1) e^(-x/theta) / (Gamma(k) theta^k)"; END;
    IF ffI == 25 THEN RETURN "f(x) = x^(α-1) (1 - x)^(β-1) / B(α, β), 0 < x < 1"; END;
    IF ffI == 26 THEN RETURN "P(X = k) = C(k+r-1, k) p^r (1 - p)^k"; END;
    IF ffI == 27 THEN RETURN "P(X = k) = C(K, k) C(N-K, n-k) / C(N, n)"; END;
  END;
  RETURN "";
END;

// The density as CAS text in the variable fdV, with the parameter
// texts fdP = {p1, p2, p3} (exact ratios, or symbols), and the limits
// it lives between: {f, lo, hi}, or {f1, lo1, hi1, f2, lo2, hi2} for
// the two-piece Laplace and Triangular.  {} = no template (discrete).
// These are the templates checked in desktop Giac, family by family
// (verify_StatDist_v2.7.py).
st_FmDens(fdI, fdP, fdV)
BEGIN
  LOCAL fdA, fdB, fdC, fdX, fdR;
  // ("" + marks them as text: they are pieces of CAS text, never numbers)
  fdA := "" + fdP(1); fdB := "" + fdP(2); fdC := "" + fdP(3); fdX := "" + fdV;
  fdR := {};
  CASE
    IF fdI == 1 THEN fdR := {"exp(-(" + fdX + "-" + fdA + ")^2/(2*" + fdB + "^2))/(" + fdB + "*sqrt(2*pi))", "-inf", "inf"}; END;
    IF fdI == 2 THEN fdR := {"1/(" + fdB + "-" + fdA + ")", fdA, fdB}; END;
    IF fdI == 3 THEN fdR := {fdA + "*exp(-" + fdA + "*" + fdX + ")", "0", "inf"}; END;
    IF fdI == 4 THEN fdR := {"(" + fdA + "/" + fdB + ")*(" + fdX + "/" + fdB + ")^(" + fdA + "-1)*exp(-(" + fdX + "/" + fdB + ")^" + fdA + ")", "0", "inf"}; END;
    IF fdI == 5 THEN fdR := {"exp(-(ln(" + fdX + ")-" + fdA + ")^2/(2*" + fdB + "^2))/(" + fdX + "*" + fdB + "*sqrt(2*pi))", "0", "inf"}; END;
    IF fdI == 6 THEN
      // two halves, so no absolute value reaches the CAS
      fdR := {"exp((" + fdX + "-" + fdA + ")/" + fdB + ")/(2*" + fdB + ")", "-inf", fdA,
        "exp(-(" + fdX + "-" + fdA + ")/" + fdB + ")/(2*" + fdB + ")", fdA, "inf"};
    END;
    IF fdI == 7 THEN fdR := {"Gamma((" + fdA + "+1)/2)/(sqrt(" + fdA + "*pi)*Gamma(" + fdA + "/2))*(1+" + fdX + "^2/" + fdA + ")^(-(" + fdA + "+1)/2)", "-inf", "inf"}; END;
    IF fdI == 8 THEN fdR := {fdX + "^(" + fdA + "/2-1)*exp(-" + fdX + "/2)/(2^(" + fdA + "/2)*Gamma(" + fdA + "/2))", "0", "inf"}; END;
    IF fdI == 9 THEN fdR := {"sqrt((" + fdA + "*" + fdX + ")^" + fdA + "*" + fdB + "^" + fdB + "/(" + fdA + "*" + fdX + "+" + fdB + ")^(" + fdA + "+" + fdB + "))/(" + fdX + "*Beta(" + fdA + "/2," + fdB + "/2))", "0", "inf"}; END;
    IF fdI == 10 THEN fdR := {fdB + "*" + fdA + "*(" + fdX + "/" + fdC + ")^" + fdA + "/(" + fdX + "*(1+(" + fdX + "/" + fdC + ")^" + fdA + ")^(" + fdB + "+1))", "0", "inf"}; END;
    IF fdI == 11 THEN fdR := {"exp(-(" + fdX + "-" + fdA + ")/" + fdB + "-exp(-(" + fdX + "-" + fdA + ")/" + fdB + "))/" + fdB, "-inf", "inf"}; END;
    IF fdI == 12 THEN fdR := {fdB + "*" + fdA + "^" + fdB + "/" + fdX + "^(" + fdB + "+1)", fdA, "inf"}; END;
    IF fdI == 13 THEN fdR := {fdX + "/" + fdA + "^2*exp(-" + fdX + "^2/(2*" + fdA + "^2))", "0", "inf"}; END;
    IF fdI == 14 THEN fdR := {"exp(-(" + fdX + "-" + fdA + ")/" + fdB + ")/(" + fdB + "*(1+exp(-(" + fdX + "-" + fdA + ")/" + fdB + "))^2)", "-inf", "inf"}; END;
    IF fdI == 18 THEN
      // min a = p1, mode c = p2, max b = p3; a side of zero width
      // (mode at an end) is left out
      IF stPar(2) > stPar(1) THEN
        fdR := {"2*(" + fdX + "-" + fdA + ")/((" + fdC + "-" + fdA + ")*(" + fdB + "-" + fdA + "))", fdA, fdB};
      END;
      IF stPar(3) > stPar(2) THEN
        fdR := CONCAT(fdR, {"2*(" + fdC + "-" + fdX + ")/((" + fdC + "-" + fdA + ")*(" + fdC + "-" + fdB + "))", fdB, fdC});
      END;
    END;
    IF fdI == 19 THEN fdR := {"1/(pi*" + fdB + "*(1+((" + fdX + "-" + fdA + ")/" + fdB + ")^2))", "-inf", "inf"}; END;
    IF fdI == 20 THEN fdR := {"(" + fdA + "/" + fdB + ")*((" + fdX + "-" + fdC + ")/" + fdB + ")^(-1-" + fdA + ")*exp(-((" + fdX + "-" + fdC + ")/" + fdB + ")^(-" + fdA + "))", fdC, "inf"}; END;
    IF fdI == 21 THEN fdR := {fdB + "*" + fdA + "*exp(" + fdA + "+" + fdB + "*" + fdX + "-" + fdA + "*exp(" + fdB + "*" + fdX + "))", "0", "inf"}; END;
    IF fdI == 22 THEN fdR := {"sqrt(2/pi)*" + fdX + "^2*exp(-" + fdX + "^2/(2*" + fdA + "^2))/" + fdA + "^3", "0", "inf"}; END;
    IF fdI == 23 THEN fdR := {fdB + "^" + fdA + "*" + fdX + "^(" + fdA + "-1)*exp(-" + fdB + "*" + fdX + ")/(" + fdA + "-1)!", "0", "inf"}; END;
    IF fdI == 24 THEN fdR := {fdX + "^(" + fdA + "-1)*exp(-" + fdX + "/" + fdB + ")/(Gamma(" + fdA + ")*" + fdB + "^" + fdA + ")", "0", "inf"}; END;
    IF fdI == 25 THEN fdR := {fdX + "^(" + fdA + "-1)*(1-" + fdX + ")^(" + fdB + "-1)/Beta(" + fdA + "," + fdB + ")", "0", "1"}; END;
  END;
  RETURN fdR;
END;

// The pmf as CAS text, for k given as the text fkK (the variable
// stqk, or a whole number in brackets).  Whole-number parameters
// (n, r, N, K) are written by st_IStr; p and lambda come from fkP.
st_FmPmf(fkI, fkP, fkK)
BEGIN
  LOCAL fkA, fkB, fkN, fkR, fkKt;
  // ("" + marks them as text: pieces of CAS text, never numbers)
  fkA := "" + fkP(1); fkB := "" + fkP(2); fkKt := "" + fkK;
  fkN := "" + st_IStr(stPar(1));
  fkR := "";
  CASE
    IF fkI == 15 THEN fkR := "comb(" + fkN + "," + fkKt + ")*" + fkB + "^" + fkKt + "*(1-" + fkB + ")^(" + fkN + "-" + fkKt + ")"; END;
    IF fkI == 16 THEN fkR := "exp(-" + fkA + ")*" + fkA + "^" + fkKt + "/" + fkKt + "!"; END;
    IF fkI == 17 THEN fkR := fkA + "*(1-" + fkA + ")^(" + fkKt + "-1)"; END;
    IF fkI == 26 THEN fkR := "comb(" + fkKt + "+" + fkN + "-1," + fkKt + ")*" + fkB + "^" + fkN + "*(1-" + fkB + ")^" + fkKt; END;
    IF fkI == 27 THEN
      fkR := "comb(" + st_IStr(stPar(2)) + "," + fkKt + ")*comb(" + fkN + "-" + st_IStr(stPar(2)) + "," + st_IStr(stPar(3)) + "-" + fkKt + ")/comb(" + fkN + "," + st_IStr(stPar(3)) + ")";
    END;
  END;
  RETURN fkR;
END;

// The parameters as CAS text: an exact ratio for each simple one, a
// symbol (stqa, stqb, stqc) for one that is not.  Sets stFmSy and
// stFmR (the symbols' values as 12-digit ratios, for the checks).
st_FmPars()
BEGIN
  LOCAL fpI, fpR, fpT, fpV, fpN;
  fpT := {"", "", ""}; stFmSy := {0, 0, 0}; stFmR := {"(0)", "(0)", "(0)"};
  fpN := {"stqa", "stqb", "stqc"};
  FOR fpI FROM 1 TO 3 DO
    fpV := stPar(fpI);
    fpR := st_Rat(fpV);
    IF SIZE(fpR) == 2 THEN
      fpT(fpI) := st_RTxt(fpR);
    ELSE
      fpT(fpI) := "(" + fpN(fpI) + ")";
      stFmSy(fpI) := 1;
      stFmR(fpI) := st_LTxt(fpV);
    END;
  END;
  stFmPT := fpT;
  RETURN fpT;
END;

// The assumptions for the symbolic parameters that must be positive,
// as " | stqa>0 and stqb>0", or "" when there are none.
st_FmAsm()
BEGIN
  LOCAL amS, amI, amP, amN;
  amS := ""; amP := st_FmPos(stId); amN := {"stqa", "stqb", "stqc"};
  FOR amI FROM 1 TO 3 DO
    IF stFmSy(amI) == 1 AND amP(amI) == 1 THEN
      IF DIM(amS) > 0 THEN amS := amS + " and "; END;
      amS := amS + amN(amI) + ">0";
    END;
  END;
  IF DIM(amS) > 0 THEN amS := " | " + amS; END;
  RETURN amS;
END;

// The integral of v^pw f(v) over every piece of fiD, as CAS text
// (pw = 0 for the CDF, 1 for E[X], 2 for E[X^2]).
st_FmInt(fiD, fiPw, fiV)
BEGIN
  LOCAL fiS, fiI, fiM, fiX;
  fiS := ""; fiM := ""; fiX := "" + fiV;    // the variable's name, as text
  IF fiPw == 1 THEN fiM := fiX + "*"; END;
  IF fiPw == 2 THEN fiM := fiX + "^2*"; END;
  FOR fiI FROM 1 TO SIZE(fiD)/3 DO
    IF fiI > 1 THEN fiS := fiS + "+"; END;
    fiS := fiS + "int(" + fiM + "(" + fiD(3*fiI - 2) + ")," + fiX + "," + fiD(3*fiI - 1) + "," + fiD(3*fiI) + ")";
  END;
  RETURN fiS;
END;

// A CAS answer fit to be shown: not empty, not an integral or sum
// the CAS left undone, not a conditional or piecewise answer.
st_FmGood(fgT)
BEGIN
  IF DIM(fgT) == 0 THEN RETURN 0; END;
  IF INSTRING(fgT, "integrate") > 0 OR INSTRING(fgT, "int(") > 0 OR INSTRING(fgT, "sum(") > 0 THEN RETURN 0; END;
  IF INSTRING(fgT, "seq(") > 0 OR INSTRING(fgT, "limit") > 0 OR INSTRING(fgT, "?") > 0 THEN RETURN 0; END;
  IF INSTRING(fgT, "when") > 0 OR INSTRING(fgT, "piecewise") > 0 OR INSTRING(fgT, "infinity") > 0 THEN RETURN 0; END;
  IF INSTRING(fgT, CHAR(8734)) > 0 THEN RETURN 0; END;     // the infinity glyph
  // v2.7a: an undone integral or sum written with its sign (the
  // integral sign, the summation sign, capital sigma) -- not measured
  // on the Prime, refused in case
  IF INSTRING(fgT, CHAR(8747)) > 0 OR INSTRING(fgT, CHAR(8721)) > 0 OR INSTRING(fgT, CHAR(931)) > 0 THEN RETURN 0; END;
  RETURN 1;
END;

// The answer at the user's own values: each symbolic parameter (and x,
// in a CDF) replaced by its 12-digit ratio, for the check.
st_FmSub(fsT)
BEGIN
  LOCAL fsS, fsI, fsN;
  fsS := fsT; fsN := {"stqa", "stqb", "stqc"};
  FOR fsI FROM 1 TO 3 DO
    IF stFmSy(fsI) == 1 THEN fsS := "subst(" + fsS + "," + fsN(fsI) + "," + stFmR(fsI) + ")"; END;
  END;
  IF INSTRING(fsS, "stqx") > 0 THEN fsS := "subst(" + fsS + ",stqx," + stFmXR + ")"; END;
  RETURN fsS;
END;

// Ask the CAS faT, and keep its answer only if the CAS's own decimal
// of it agrees with the program's value faV to 6 significant digits.
// Returns the answer, or "" (the page then says "no closed form
// shown").  This is the guard: a CAS slip (desktop Giac gives the
// Logistic's E[X^2] as -pi^2/3) never reaches the screen.
st_FmAsk(faT, faV)
BEGIN
  LOCAL faR, faQ;
  faR := st_Cas(faT);
  IF st_FmGood(faR) == 0 THEN RETURN ""; END;
  faQ := st_CasNum(st_FmSub(faR), st_Scale(faV));
  IF faQ(2) == 0 THEN RETURN ""; END;
  IF ABS(faQ(1) - faV) > 0.000001*MAX(ABS(faV), 0.000000001) THEN RETURN ""; END;
  RETURN faR;
END;

// Did the CAS find the integral infinite?  v2.7a: st_Cas refuses any
// answer holding the infinity sign (the Prime prints infinity as
// CHAR(8734), CASProbe beta 1b T41), so the question the CAS is asked
// is exp(-(integral)), and its answer is "0" exactly when the integral
// is plus infinity.  A finite integral v gives e^(-v), never 0; minus
// infinity gives an infinity, which st_Cas refuses; an integral the
// CAS cannot do stays unevaluated.  So only "0" means infinite.
st_FmInf(fnT)
BEGIN
  IF fnT == "0" THEN RETURN 1; END;
  RETURN 0;
END;

// CAS text made readable: * as the times sign (^ and / are kept), the
// tagged CAS names as plain letters, the parameters by their names.
st_FmShow(fwT)
BEGIN
  LOCAL fwS, fwN;
  fwN := st_FmNm(stId);
  fwS := st_Repl(fwT, "*", CHAR(215));
  fwS := st_Repl(fwS, "+infinity", "infinity");
  fwS := st_Repl(fwS, "stqx", "x");
  fwS := st_Repl(fwS, "stqt", "t");
  fwS := st_Repl(fwS, "stqk", "k");
  fwS := st_Repl(fwS, "stqa", fwN(1));
  fwS := st_Repl(fwS, "stqb", fwN(2));
  fwS := st_Repl(fwS, "stqc", fwN(3));
  RETURN fwS;
END;

// A limit of integration as words: "infinity", or the parameter it
// is, as the user typed it (or by its name, when it is a symbol).
st_FmLimS(flT)
BEGIN
  LOCAL flI, flN;
  IF flT == "-inf" THEN RETURN "-infinity"; END;
  IF flT == "inf" THEN RETURN "infinity"; END;
  flN := st_FmNm(stId);
  FOR flI FROM 1 TO 3 DO
    IF flT == stFmPT(flI) THEN
      IF stFmSy(flI) == 1 THEN RETURN flN(flI); END;
      RETURN st_fmtx(stPar(flI));
    END;
  END;
  RETURN st_FmShow(flT);
END;

// One more line (and its colour) for the page, and a long one wrapped
// to the screen (font 1 is about 6.5 px a character: 46 fit; the
// lines after the first are indented).
st_FmAdd(faL, faC)
BEGIN
  stFmL := CONCAT(stFmL, {faL});
  stFmC := CONCAT(stFmC, {faC});
END;

st_FmWrap(fwL, fwC)
BEGIN
  LOCAL fwW, fwI;
  fwW := st_Wrap(fwL, 46);
  FOR fwI FROM 1 TO SIZE(fwW) DO
    IF fwI == 1 THEN st_FmAdd(fwW(fwI), fwC); ELSE st_FmAdd("   " + fwW(fwI), fwC); END;
  END;
END;

// A checked CAS answer and the program's matching decimal.  An answer
// that will not fit is described, never cut: a fraction (only digits,
// / and - in it) longer than one line becomes "an exact fraction of
// 77 digits"; any other form longer than three lines, "an exact form
// of 300 characters".  The decimal below it is the program's.
st_FmRes(frT, frV)
BEGIN
  LOCAL frS, frI, frF;
  frS := st_FmShow(frT);
  frF := 1;
  FOR frI FROM 1 TO DIM(frS) DO
    IF INSTRING("0123456789/-", MID(frS, frI, 1)) == 0 THEN frF := 0; END;
  END;
  IF frF AND DIM(frS) > 42 THEN
    frS := "an exact fraction of " + st_IStr(st_Digits(frS)) + " digits";
  ELSE
    IF DIM(frS) > 130 THEN frS := "an exact form of " + st_IStr(DIM(frS)) + " characters"; END;
  END;
  st_FmWrap("  = " + frS, stMAG);
  // v2.7d: st_Eq -- the relation without "=" ("P(X <= -6) < 0.000001")
  st_FmAdd(" " + st_Eq(st_fmtv(frV, 6)) + ", the program's value: agree", stGRN);
  stCasOK := 1;
END;

// No CAS answer to show: say so, with the program's own decimal.
st_FmNone(fnV)
BEGIN
  st_FmAdd("  no closed form shown", stGRY);
  st_FmAdd("  the program's value: " + st_fmtv(fnV, 6), stGRN);
END;

// ---- the continuous families ----
st_FmCont(fcP)
BEGIN
  LOCAL fcD, fcN, fcM, fcE, fcSf;
  LOCAL fcLim, fcV, fcX, fcH, fcR;
  fcSf := st_FmAsm();
  fcD := st_FmDens(stId, fcP, "stqx");
  fcN := IP(SIZE(fcD)/3);
  IF fcN < 1 THEN st_FmNone(stMu); RETURN 0; END;
  fcLim := "x from " + st_FmLimS(fcD(2)) + " to " + st_FmLimS(fcD(3*fcN)) + ":";
  // ---- E[X] ----
  st_FmWrap("Mean: E[X] = integral of x f(x) dx, " + fcLim, #000000h);
  fcM := "";
  IF TYPE(stMu) == 0 THEN
    fcM := st_FmAsk("simplify(" + st_FmInt(fcD, 1, "stqx") + ")" + fcSf, stMu);
    IF DIM(fcM) > 0 THEN st_FmRes(fcM, stMu); ELSE st_FmNone(stMu); END;
  ELSE
    // no mean: the CAS should find the integral infinite.  For the
    // symmetric t and Cauchy, the half from the centre up is enough.
    fcH := fcD;
    IF stId == 7 THEN fcH(2) := "0"; END;
    IF stId == 19 THEN fcH(2) := fcP(1); END;
    // v2.7a: asked as exp(-(integral)), which is 0 when it is infinite
    fcR := st_Cas("exp(-(simplify(" + st_FmInt(fcH, 1, "stqx") + ")" + fcSf + "))");
    IF st_FmInf(fcR) THEN
      IF stId == 7 OR stId == 19 THEN
        st_FmAdd("  from the centre up alone: infinity (CAS).", stMAG);
      ELSE
        st_FmAdd("  = infinity (CAS).", stMAG);
      END;
      st_FmAdd("  So there is NO mean (the program: undef).", stGRN);
      stCasOK := 1;
    ELSE
      st_FmAdd("  no closed form shown (the program: undef)", stGRY);
    END;
  END;
  // ---- E[X^2] ----
  st_FmAdd("E[X^2] = integral of x^2 f(x) dx, same limits:", #000000h);
  fcE := "";
  IF TYPE(stMu) == 0 AND TYPE(stVr) == 0 THEN
    fcV := stVr + stMu^2;
    fcE := st_FmAsk("simplify(" + st_FmInt(fcD, 2, "stqx") + ")" + fcSf, fcV);
    IF DIM(fcE) > 0 THEN st_FmRes(fcE, fcV); ELSE st_FmNone(fcV); END;
  ELSE
    // v2.7a: asked as exp(-(integral)), which is 0 when it is infinite
    fcR := st_Cas("exp(-(simplify(" + st_FmInt(fcD, 2, "stqx") + ")" + fcSf + "))");
    IF st_FmInf(fcR) THEN
      st_FmAdd("  = infinity (CAS): there is NO variance", stMAG);
      st_FmAdd("  (the program: undef).", stGRN);
      stCasOK := 1;
    ELSE
      st_FmAdd("  no closed form shown (the program: undef)", stGRY);
    END;
  END;
  // ---- Var(X), from the two answers above ----
  IF DIM(fcM) > 0 AND DIM(fcE) > 0 THEN
    st_FmAdd("Var(X) = E[X^2] - E[X]^2:", #000000h);
    fcR := st_FmAsk("simplify((" + fcE + ")-(" + fcM + ")^2)" + fcSf, stVr);
    IF DIM(fcR) > 0 THEN st_FmRes(fcR, stVr); ELSE st_FmNone(stVr); END;
  END;
  // ---- the CDF, for the Exponential, Weibull and Pareto ----
  IF stId == 3 OR stId == 4 OR stId == 12 THEN
    fcX := stXv;
    IF stTl == 3 THEN fcX := stBv; END;
    fcD := st_FmDens(stId, fcP, "stqt");
    st_FmWrap("CDF: F(x) = integral of f(t) dt, t = " + st_FmLimS(fcD(2)) + " to x:", #000000h);
    fcR := st_FmAsk("simplify(int(" + fcD(1) + ",stqt," + fcD(2) + ",stqx))" + fcSf, st_cdf(stId, stPar, fcX));
    IF DIM(fcR) > 0 THEN
      st_FmWrap("  F(x) = " + st_FmShow(fcR), stMAG);
      // v2.7d: st_Eq -- the relation without "=" ("P(X <= -6) < 0.000001")
      // v2.7e: x through st_XStr -- a value, never "< 0.0001" (B7)
      st_FmAdd("  at x" + st_Eq(st_XStr(fcX, st_fmtx(fcX))) + ": " + st_fmtp(st_cdf(stId, stPar, fcX), 6) + ", the program's: agree", stGRN);
      stCasOK := 1;
    ELSE
      st_FmNone(st_cdf(stId, stPar, fcX));
    END;
  END;
  RETURN 1;
END;

// ---- the discrete families ----
st_FmDisc(ddP)
BEGIN
  LOCAL ddLo, ddA, ddB, ddPm, ddR;
  LOCAL ddBig, ddHi, ddS, ddV;
  // the smallest k the distribution can take
  ddLo := 0;
  IF stId == 17 THEN ddLo := 1; END;
  IF stId == 27 THEN ddLo := MAX(0, stPar(3) - (stPar(1) - stPar(2))); END;
  // the k's the probability on the plot adds up
  IF stTl == 3 THEN ddA := stAv; ddB := stBv; ELSE ddA := ddLo; ddB := stXv; END;
  // an exact sum is only asked for at moderate sizes (the fractions
  // grow with every term, and the G1's CAS speed is unknown)
  ddBig := 0;
  IF stId == 15 AND stPar(1) > 200 THEN ddBig := 1; END;
  IF stId == 16 AND ddB > 150 THEN ddBig := 1; END;
  IF stId == 17 AND ddB > 500 THEN ddBig := 1; END;
  IF stId == 26 AND (ddB > 200 OR stPar(1) > 200) THEN ddBig := 1; END;
  ddV := st_Prob();
  IF ddBig THEN
    st_FmAdd("These numbers are too large for an exact sum", stGRY);
    st_FmAdd("here (n or k over 150-500).", stGRY);
    st_FmAdd("  the program's value: " + st_fmtp(ddV, 6), stGRN);
    RETURN 0;
  END;
  ddPm := st_FmPmf(stId, ddP, "stqk");
  // ---- the probability on the plot, as a sum ----
  IF stTl == 1 THEN st_FmAdd("P(X <= " + st_IStr(stXv) + ") = the sum of P(X = k),", #000000h); END;
  IF stTl == 2 THEN st_FmAdd("P(X > " + st_IStr(stXv) + ") = 1 - the sum of P(X = k),", #000000h); END;
  IF stTl == 3 THEN st_FmAdd("P(" + st_IStr(stAv) + " <= X <= " + st_IStr(stBv) + ") = the sum of P(X = k),", #000000h); END;
  st_FmAdd("  k from " + st_IStr(ddA) + " to " + st_IStr(ddB) + ":", #000000h);
  IF ddB < ddA THEN
    // no k at all in the range: nothing to add up
    IF stTl == 2 THEN st_FmAdd("  = 1 (no k to take away)", stMAG); ELSE st_FmAdd("  = 0 (no k in the range)", stMAG); END;
  ELSE
    ddS := "sum(seq(" + ddPm + ",stqk," + st_IStr(ddA) + "," + st_IStr(ddB) + "))";
    IF stTl == 2 THEN ddS := "1-" + ddS; END;
    ddR := st_FmAsk("simplify(" + ddS + ")" + st_FmAsm(), ddV);
    IF DIM(ddR) > 0 THEN st_FmRes(ddR, ddV); ELSE st_FmNone(ddV); END;
  END;
  // ---- the red bar, P(X = k), for the two tails ----
  IF stTl <> 3 THEN
    st_FmAdd("P(X = " + st_IStr(stXv) + "), the red bar:", #000000h);
    ddV := st_pdf(stId, stPar, stXv);
    ddR := st_FmAsk("simplify(" + st_FmPmf(stId, ddP, "(" + st_IStr(stXv) + ")") + ")" + st_FmAsm(), ddV);
    IF DIM(ddR) > 0 THEN st_FmRes(ddR, ddV); ELSE st_FmNone(ddV); END;
  END;
  // ---- the mean as a finite sum: Binomial and Hypergeometric ----
  IF (stId == 15 OR stId == 27) AND TYPE(stMu) == 0 THEN
    ddHi := stPar(1);
    IF stId == 27 THEN ddHi := MIN(stPar(3), stPar(2)); END;
    st_FmAdd("E[X] = the sum of k P(X = k), k = " + st_IStr(ddLo) + " to " + st_IStr(ddHi) + ":", #000000h);
    ddR := st_FmAsk("simplify(sum(seq(stqk*(" + ddPm + "),stqk," + st_IStr(ddLo) + "," + st_IStr(ddHi) + ")))" + st_FmAsm(), stMu);
    IF DIM(ddR) > 0 THEN st_FmRes(ddR, stMu); ELSE st_FmNone(stMu); END;
  END;
  RETURN 1;
END;

// The plot's state, as the Formula page's key: the page is rebuilt
// (and the CAS asked) only when this changes.
st_FmKeyNow()
BEGIN
  RETURN {stId, stPar(1), stPar(2), stPar(3), stTl, stXv, stAv, stBv};
END;

st_FmSame()
BEGIN
  LOCAL fkN, fkI;
  fkN := st_FmKeyNow();
  IF SIZE(stFmK) <> SIZE(fkN) THEN RETURN 0; END;
  FOR fkI FROM 1 TO SIZE(fkN) DO
    IF stFmK(fkI) <> fkN(fkI) THEN RETURN 0; END;
  END;
  RETURN 1;
END;

// Build the page: ask the CAS everything it needs, once.
st_FmBuild()
BEGIN
  LOCAL fbP, fbI, fbN, fbX;
  // the CAS's speed on a G1 is not known yet: never a blank wait
  DIMGROB_P(G1, 320, 240, #FFFFFFh);
  TEXTOUT_P("Asking the CAS ...", G1, 40, 96, 2, stBLU);
  TEXTOUT_P("exact arithmetic: it can take a few", G1, 40, 122, 1, stGRY);
  TEXTOUT_P("seconds, once (the answer is kept)", G1, 40, 134, 1, stGRY);
  BLIT_P(G0, 0, 0, G1);
  stFmL := {}; stFmC := {}; stCasOK := 0;
  st_FmAdd(st_ParStr(stId, stPar), #000000h);
  st_FmWrap(st_FmF(stId), stBLU);
  fbP := st_FmPars();
  fbX := stXv;
  IF stTl == 3 THEN fbX := stBv; END;
  stFmXR := st_LTxt(fbX);
  // a parameter that is not a simple ratio stays a symbol: say so
  fbN := st_FmNm(stId);
  FOR fbI FROM 1 TO 3 DO
    IF stFmSy(fbI) == 1 AND DIM(fbN(fbI)) > 0 THEN
      st_FmWrap(fbN(fbI) + " = " + STRING(stPar(fbI), 2, 12) + " is not a simple ratio, so the CAS keeps it as the symbol " + fbN(fbI) + ":", stGRY);
    END;
  END;
  IF stDisc == 1 THEN st_FmDisc(fbP); ELSE st_FmCont(fbP); END;
  stFmOK := stCasOK;
  stFmK := st_FmKeyNow();
END;

// The page: coloured lines in font 1, 12 px apart, 16 to a screen
// (more turn over), the house footer and the badge in dodger blue.
st_FmPage(fgT, fgL, fgC, fgOk)
BEGIN
  LOCAL fgN, fgF, fgLa, fgI, fgY, fgB;
  fgN := SIZE(fgL); fgF := 1;
  fgB := "decimal";
  IF fgOk THEN fgB := "exact (CAS)"; END;
  REPEAT
    fgLa := MIN(fgN, fgF + 15);
    DIMGROB_P(G1, 320, 240, #FFFFFFh);
    RECT_P(G1, 0, 0, 319, 21, #FF4500h, #FF4500h);   // orange title band
    TEXTOUT_P(fgT, G1, 6, 4, 2, #FFFFFFh, 308);
    fgY := 26;
    FOR fgI FROM fgF TO fgLa DO
      TEXTOUT_P(fgL(fgI), G1, 6, fgY, 1, fgC(fgI), 308);
      fgY := fgY + 12;
    END;
    IF fgLa < fgN THEN
      TEXTOUT_P("any key or tap = next page", G1, 6, 224, 1, stFOOT, 220);
    ELSE
      TEXTOUT_P("any key or tap = back", G1, 6, 224, 1, stFOOT, 220);
    END;
    TEXTOUT_P(fgB, G1, 314 - 7*DIM(fgB), 224, 1, stFOOT, 90);
    BLIT_P(G0, 0, 0, G1);
    st_WaitAny();                                     // drain, wait, drain
    fgF := fgLa + 1;
  UNTIL fgF > fgN;
END;

// THE FORMULA PAGE: built (with the CAS) only if the plot changed since
// it was last built; otherwise the kept page is shown at once.
st_Formula()
BEGIN
  IF st_FmSame() == 0 THEN st_FmBuild(); END;
  stCasOK := stFmOK;
  st_FmPage(st_Name(stId) + ": the formulas", stFmL, stFmC, stFmOK);
END;

// Enter: the options.  The overlay item depends on the distribution,
// so there are three menus; opMap turns the item into an action.
// Returns 0 to leave the plot, 1 to stay.
st_Options()
BEGIN
  LOCAL opI, opMap, opA;
  opI := 0;
  // v2.7: "Formula page (CAS)" is item 2 in all three menus (action 10)
  IF stId == 15 OR stId == 16 THEN
    CHOOSE(opI, "Plot options", {"Numbers page", "Formula page (CAS)",
      "CDF curve on/off", "Normal approx. on/off", "Between a and b",
      "Edit parameters", "Refit the window", "Help", "Back to the menu"});
    opMap := {1, 10, 2, 3, 4, 5, 6, 7, 8};
  ELSE
    IF stId == 7 THEN
      CHOOSE(opI, "Plot options", {"Numbers page", "Formula page (CAS)",
        "CDF curve on/off", "N(0,1) behind t on/off", "Between a and b",
        "Edit parameters", "Refit the window", "Help", "Back to the menu"});
      opMap := {1, 10, 2, 9, 4, 5, 6, 7, 8};
    ELSE
      CHOOSE(opI, "Plot options", {"Numbers page", "Formula page (CAS)",
        "CDF curve on/off", "Between a and b", "Edit parameters",
        "Refit the window", "Help", "Back to the menu"});
      opMap := {1, 10, 2, 4, 5, 6, 7, 8};
    END;
  END;
  IF opI == 0 THEN RETURN 1; END;
  opA := opMap(opI);
  CASE
    IF opA == 1 THEN st_Numbers(); END;
    IF opA == 10 THEN st_Formula(); END;
    IF opA == 2 THEN stOvC := 1 - stOvC; END;
    IF opA == 3 THEN stOvN := 1 - stOvN; st_Cache(); END;
    IF opA == 9 THEN stOvT := 1 - stOvT; st_Cache(); END;
    IF opA == 4 THEN st_SetTail(3); END;
    IF opA == 5 THEN st_Edit(); END;
    IF opA == 6 THEN st_Fit(1); st_Cache(); END;
    IF opA == 7 THEN st_show_help(); END;
    IF opA == 8 THEN RETURN 0; END;
  END;
  RETURN 1;
END;

// ---- the plot screen ----
// Draw, then wait for a key or a touch (both polled together, Lessons
// Learned s6).  A touch at y >= 218 is a soft key, slot MIN(6, IP(x/
// 53.4) + 1); a touch in or just around the plot box moves x; keys are
// the shortcuts, so everything works without touching the screen.
st_Plot()
BEGIN
  LOCAL plK, plM, plM1, plGo, plSl, plK2;
  stDisc := st_Disc(stId);
  stPh := 0; stAct := 2; stDrg := 0; stSKP := 1;
  st_Clamp();
  st_Fit(1);
  st_Cache();
  plGo := 1;
  WHILE plGo == 1 DO
    st_Draw();
    REPEAT plK := GETKEY; UNTIL plK == -1;      // drain: a key left by a form
    st_Drain();                                 // ... or a finger still down
    REPEAT
      plK := GETKEY;
      plM := MOUSE;
      plM1 := plM(1);
    UNTIL plK > -1 OR SIZE(plM1) > 0;
    IF SIZE(plM1) > 0 THEN
      IF plM1(2) >= 218 THEN
        st_Drain();
        plSl := MIN(6, IP(plM1(1)/53.4) + 1);
        IF plSl == 6 THEN
          stSKP := stSKP + 1;                     // More: the next page (v2.7: of three)
          IF stSKP > 3 THEN stSKP := 1; END;
        ELSE
          IF stSKP == 1 THEN
            IF plSl <= 3 THEN st_SetTail(plSl); END;
            IF plSl == 4 THEN st_Zoom(2); END;    // Zoom-: twice as wide
            IF plSl == 5 THEN st_Zoom(0.5); END;  // Zoom+: half as wide
          ELSE
            IF stSKP == 2 THEN
              IF plSl == 1 THEN st_Edit(); END;
              IF plSl == 2 THEN st_Numbers(); END;
              IF plSl == 3 THEN stOvC := 1 - stOvC; END;
              IF plSl == 4 THEN
                // the overlay this distribution has, else a refit
                CASE
                  IF stId == 15 OR stId == 16 THEN stOvN := 1 - stOvN; st_Cache(); END;
                  IF stId == 7 THEN stOvT := 1 - stOvT; st_Cache(); END;
                  DEFAULT st_Fit(1); st_Cache();
                END;
              END;
              IF plSl == 5 THEN plGo := 0; END;   // Menu: leave the plot
            ELSE
              // v2.7, page 3: the Formula page (the CAS), Help, Menu
              IF plSl == 1 THEN st_Formula(); END;
              IF plSl == 2 THEN st_show_help(); END;
              IF plSl == 5 THEN plGo := 0; END;   // Menu: leave the plot
            END;
          END;
        END;
      ELSE
        IF plM1(1) >= stPL - 6 AND plM1(1) <= stPR + 6 AND plM1(2) >= stPT AND plM1(2) <= stPB + 12 THEN
          st_Drag(plM1(1));
        ELSE
          st_Drain();                             // a tap elsewhere does nothing
        END;
      END;
    ELSE
      REPEAT plK2 := GETKEY; UNTIL plK2 == -1;     // drain this key
      CASE
        IF plK == 4 THEN plGo := 0; END;          // Esc: back to the menu
        IF plK == 7 THEN st_Step(-1); END;        // Left
        IF plK == 8 THEN st_Step(1); END;         // Right
        IF plK == 2 OR plK == 12 THEN st_UpDn(plK); END;
        IF plK == 50 THEN st_Zoom(0.5); END;      // +
        IF plK == 45 THEN st_Zoom(2); END;        // -
        IF plK == 42 THEN st_Fit(1); st_Cache(); END;   // 0: refit
        IF plK == 30 THEN plGo := st_Options(); END;    // Enter
      END;
    END;
  END;
END;

// ---- one distribution, start to finish ----
// The tail first, then the form, then the plot (v2.5's order).  rCur
// is what the form shows: the defaults, or the gallery's values.
// rOvT = 1 opens Student t with N(0,1) already behind it.
st_Run(dId, rCur, rOvT)
BEGIN
  LOCAL rTl, rSt;
  rTl := st_get_tail_type();
  IF rTl == 0 THEN RETURN 0; END;
  rSt := st_Ask(dId, rTl, rCur, 1);
  IF TYPE(rSt) <> 6 THEN RETURN 0; END;
  stId := dId; stTl := rTl;
  stPar := {rSt(1), rSt(2), rSt(3)};
  stXv := rSt(4); stAv := rSt(5); stBv := rSt(6);
  stOvN := 0; stOvC := 0; stOvT := rOvT;
  st_Plot();
  RETURN 1;
END;


// ============================================================
// 8. THE GALLERY
// Six thumbnails, computed once into one flat list stGH (pixel
// heights 0..64): curves 1-4 at 97 columns each (entries 1..388),
// N(0,1) for behind t (389..485), Binomial(10, 0.3) bars (486..496)
// and Poisson(3) bars (497..507).  Redrawing a frame then costs no
// maths at all -- only the selection frame moves.
// ============================================================
st_GalBuild()
BEGIN
  LOCAL gbT, gbC, gbX, gbV, gbM;
  stGH := MAKELIST(0, K, 1, 507, 1);
  FOR gbT FROM 1 TO 5 DO
    FOR gbC FROM 0 TO 96 DO
      gbV := 0; gbM := 0.43;
      CASE
        IF gbT == 1 THEN gbV := NORMALD(0, 1, -3.5 + gbC*7/96); END;
        IF gbT == 2 THEN gbV := EXP(-gbC*5/96); gbM := 1.08; END;
        IF gbT == 3 THEN
          gbX := gbC*14/96;
          IF gbX > 0 THEN gbV := CHISQUARE(4, gbX); END;
          gbM := 0.2;
        END;
        IF gbT == 4 THEN gbV := STUDENT(3, -4.5 + gbC*9/96); END;
        DEFAULT gbV := NORMALD(0, 1, -4.5 + gbC*9/96);
      END;
      stGH((gbT - 1)*97 + gbC + 1) := IP(MIN(gbV/gbM, 1)*64 + 0.5);
    END;
  END;
  FOR gbC FROM 0 TO 10 DO
    stGH(486 + gbC) := IP(BINOMIAL(10, 0.3, gbC)/0.29*64 + 0.5);
    stGH(497 + gbC) := IP(POISSON(3, gbC)/0.29*64 + 0.5);
  END;
END;

// Compose the gallery from stGH, blit once, then the soft keys.  The
// tap boxes stGX, stGY are recorded while drawing (ChiTable's idiom),
// so the drawing and the touch code cannot disagree.
st_GalDraw()
BEGIN
  LOCAL gdI, gdC, gdR, gdX, gdY, gdB;
  LOCAL gdS, gdH, gdP, gdK, gdLb, gdF;
  LOCAL gdMn;
  gdLb := {"Normal(0, 1)", "Exponential(1)", "Chi-square(4)",
    "t(3) and N(0,1)", "Binomial(10,0.3)", "Poisson(3)"};
  gdMn := {49, 20, 28, 49, 32, 32};      // pixel offset of each mean
  DIMGROB_P(G1, 320, 240, #FFFFFFh);
  RECT_P(G1, 0, 0, 319, 21, stORN, stORN);
  TEXTOUT_P("Distribution families", G1, 6, 4, 2, #FFFFFFh, 190);
  TEXTOUT_P("tap one to open it", G1, 204, 6, 1, #FFFFFFh, 112);
  FOR gdI FROM 1 TO 6 DO
    gdC := gdI - 1 - 3*IP((gdI - 1)/3);   // column 0..2
    gdR := IP((gdI - 1)/3);               // row 0..1
    gdX := 6 + gdC*104; gdY := 26 + gdR*96;
    stGX(gdI) := gdX; stGY(gdI) := gdY;
    gdF := stGRY;
    IF gdI == stGSel THEN
      RECT_P(G1, gdX - 2, gdY - 2, gdX + 100, gdY + 92, stRED, stRED);
      gdF := stRED;
    END;
    RECT_P(G1, gdX, gdY, gdX + 98, gdY + 90, gdF, #FFFFFFh);
    TEXTOUT_P(gdLb(gdI), G1, gdX + 4, gdY + 3, 1, #000000h, 92);
    gdB := gdY + 84;                      // the thumbnail's baseline
    IF gdI <= 4 THEN
      gdS := (gdI - 1)*97;
      FOR gdK FROM 0 TO 96 DO             // shade under the curve
        gdH := stGH(gdS + gdK + 1);
        IF gdH > 0 THEN LINE_P(G1, gdX + 1 + gdK, gdB - 1, gdX + 1 + gdK, gdB - gdH, stSHD); END;
      END;
      IF gdI == 4 THEN                    // N(0,1) in pale grey behind t
        FOR gdK FROM 1 TO 96 DO
          LINE_P(G1, gdX + gdK, gdB - stGH(388 + gdK), gdX + gdK + 1, gdB - stGH(389 + gdK), stPALE);
        END;
      END;
      FOR gdK FROM 1 TO 96 DO             // the curve
        LINE_P(G1, gdX + gdK, gdB - stGH(gdS + gdK), gdX + gdK + 1, gdB - stGH(gdS + gdK + 1), stBLU);
      END;
    ELSE
      gdS := 485 + (gdI - 5)*11;          // eleven bars, k = 0..10
      FOR gdK FROM 0 TO 10 DO
        gdH := stGH(gdS + gdK + 1);
        gdP := gdX + 5 + gdK*8;
        IF gdH > 0 THEN RECT_P(G1, gdP, gdB - gdH, gdP + 6, gdB - 1, stBLU, stSHD); END;
      END;
    END;
    LINE_P(G1, gdX + 1, gdB, gdX + 97, gdB, stGRY);
    gdP := gdX + gdMn(gdI);               // the mean, green dots
    FOR gdH FROM gdY + 16 TO gdB STEP 4 DO PIXON_P(G1, gdP, gdH, stGRN); END;
  END;
  BLIT_P(G0, 0, 0, G1);
  DRAWMENU({"Open", "", "", "", "Help", "Back"});
END;

// Open thumbnail goI: its distribution's form, with the thumbnail's
// parameters as the defaults ({p1, p2, p3, x, a, b}).
st_GalOpen(goI)
BEGIN
  IF goI == 1 THEN st_Run(1, {0, 1, 0, 1, -1, 1}, 0); END;
  IF goI == 2 THEN st_Run(3, {1, 0, 0, 1, 0, 1}, 0); END;
  IF goI == 3 THEN st_Run(8, {4, 0, 0, 5, 2, 8}, 0); END;
  IF goI == 4 THEN st_Run(7, {3, 0, 0, 2, -2, 2}, 1); END;
  IF goI == 5 THEN st_Run(15, {10, 0.3, 0, 3, 2, 4}, 0); END;
  IF goI == 6 THEN st_Run(16, {3, 0, 0, 3, 1, 5}, 0); END;
END;

// The gallery screen: tap a thumbnail (or move the red frame with the
// arrows and press Enter, or the Open soft key); Back or Esc returns.
st_Gallery()
BEGIN
  LOCAL glK, glM, glM1, glGo, glI, glHit;
  LOCAL glSl, glK2;
  st_GalBuild();
  glGo := 1;
  WHILE glGo == 1 DO
    st_GalDraw();
    REPEAT glK := GETKEY; UNTIL glK == -1;
    st_Drain();
    REPEAT
      glK := GETKEY;
      glM := MOUSE;
      glM1 := glM(1);
    UNTIL glK > -1 OR SIZE(glM1) > 0;
    glHit := 0;
    IF SIZE(glM1) > 0 THEN
      st_Drain();
      IF glM1(2) >= 218 THEN
        glSl := MIN(6, IP(glM1(1)/53.4) + 1);
        IF glSl == 1 THEN glHit := stGSel; END;
        IF glSl == 5 THEN
          st_Page("Gallery", {"Six distribution families, drawn",
            "with fixed parameters.", "",
            "Tap one to open it: the tail menu",
            "and the form come up with these",
            "parameters filled in.", "",
            "Arrows = move the red frame,",
            "Ent or Open = open it,",
            "Esc or Back = the main menu."});
        END;
        IF glSl == 6 THEN glGo := 0; END;
      ELSE
        FOR glI FROM 1 TO 6 DO
          IF glM1(1) >= stGX(glI) AND glM1(1) <= stGX(glI) + 98 AND glM1(2) >= stGY(glI) AND glM1(2) <= stGY(glI) + 90 THEN
            glHit := glI;
          END;
        END;
      END;
    ELSE
      REPEAT glK2 := GETKEY; UNTIL glK2 == -1;
      IF glK == 7 THEN stGSel := MAX(1, stGSel - 1); END;
      IF glK == 8 THEN stGSel := MIN(6, stGSel + 1); END;
      IF glK == 2 AND stGSel > 3 THEN stGSel := stGSel - 3; END;
      IF glK == 12 AND stGSel <= 3 THEN stGSel := stGSel + 3; END;
      IF glK == 30 THEN glHit := stGSel; END;
      IF glK == 4 THEN glGo := 0; END;
    END;
    IF glHit > 0 THEN
      stGSel := glHit;
      st_GalOpen(glHit);
    END;
  END;
END;


// ============================================================
// 9. THE MENUS
// Each item calls st_Run with the distribution's number and its
// default form values.
// ============================================================
st_common_continuous()
BEGIN
  LOCAL idx;
  idx := 0;
  CHOOSE(idx, "Common Continuous",
    "Normal", "Uniform", "Exponential", "Triangular",
    "Weibull", "LogNormal", "Cauchy", "Laplace", "Back");
  IF idx == 0 THEN RETURN; END;
  CASE
    IF idx == 1 THEN st_Run(1, st_Def(1), 0); END;
    IF idx == 2 THEN st_Run(2, st_Def(2), 0); END;
    IF idx == 3 THEN st_Run(3, st_Def(3), 0); END;
    IF idx == 4 THEN st_Run(18, st_Def(18), 0); END;
    IF idx == 5 THEN st_Run(4, st_Def(4), 0); END;
    IF idx == 6 THEN st_Run(5, st_Def(5), 0); END;
    IF idx == 7 THEN st_Run(19, st_Def(19), 0); END;
    IF idx == 8 THEN st_Run(6, st_Def(6), 0); END;
  END;
END;

st_advanced_continuous()
BEGIN
  LOCAL idx;
  idx := 0;
  // v2.7e (B5): "Chi-square" and "Maxwell-Boltzmann", spelt as on the
  // plot and the Numbers page (st_Name).  No comment goes inside the
  // list's brackets: a // there is a syntax error on the Prime.
  CHOOSE(idx, "Advanced Continuous",
    {"Student t", "Chi-square", "Fisher F", "Gamma", "Beta",
    "Burr Type XII", "Weibull", "Gumbel", "Frechet", "Pareto",
    "Rayleigh", "Gompertz", "Erlang", "Maxwell-Boltzmann", "Logistic", "Back"});
  IF idx == 0 THEN RETURN; END;
  CASE
    IF idx == 1 THEN st_Run(7, st_Def(7), 0); END;
    IF idx == 2 THEN st_Run(8, st_Def(8), 0); END;
    IF idx == 3 THEN st_Run(9, st_Def(9), 0); END;
    IF idx == 4 THEN st_Run(24, st_Def(24), 0); END;
    IF idx == 5 THEN st_Run(25, st_Def(25), 0); END;
    IF idx == 6 THEN st_Run(10, st_Def(10), 0); END;
    IF idx == 7 THEN st_Run(4, st_Def(4), 0); END;
    IF idx == 8 THEN st_Run(11, st_Def(11), 0); END;
    IF idx == 9 THEN st_Run(20, st_Def(20), 0); END;
    IF idx == 10 THEN st_Run(12, st_Def(12), 0); END;
    IF idx == 11 THEN st_Run(13, st_Def(13), 0); END;
    IF idx == 12 THEN st_Run(21, st_Def(21), 0); END;
    IF idx == 13 THEN st_Run(23, st_Def(23), 0); END;
    IF idx == 14 THEN st_Run(22, st_Def(22), 0); END;
    IF idx == 15 THEN st_Run(14, st_Def(14), 0); END;
  END;
END;

st_discrete_distributions()
BEGIN
  LOCAL idx;
  idx := 0;
  CHOOSE(idx, "Discrete Distributions",
    "Binomial", "Poisson", "Geometric",
    "Negative Binomial", "Hypergeometric", "Back");
  IF idx == 0 THEN RETURN; END;
  CASE
    IF idx == 1 THEN st_Run(15, st_Def(15), 0); END;
    IF idx == 2 THEN st_Run(16, st_Def(16), 0); END;
    IF idx == 3 THEN st_Run(17, st_Def(17), 0); END;
    IF idx == 4 THEN st_Run(26, st_Def(26), 0); END;
    IF idx == 5 THEN st_Run(27, st_Def(27), 0); END;
  END;
END;

// ═════════════════════════════════════════════════════════════════
// The entry point.  v2.6: Gallery is the first item; the function is
// at the end of the file so that every routine it reaches is defined
// before it.  Returns the last probability shown, for Home's history.
EXPORT StatDist()
BEGIN
  LOCAL idx;
  st_Init();
  REPEAT
    idx := 0;
    CHOOSE(idx, "Statistical Distributions",
      "Gallery", "Common Continuous", "Advanced Continuous",
      "Discrete Distributions", "About/Help", "Exit");
    IF idx == 0 OR idx == 6 THEN RETURN "StatDist: " + stLast; END;
    CASE
      IF idx == 1 THEN st_Gallery(); END;
      IF idx == 2 THEN st_common_continuous(); END;
      IF idx == 3 THEN st_advanced_continuous(); END;
      IF idx == 4 THEN st_discrete_distributions(); END;
      IF idx == 5 THEN st_show_help(); END;
    END;
  UNTIL 0;
END;
