bug#58929: 29.0.50; Calc: finding roots utpn doesn't work as advertised

[Matt Armstrong <matt@rfc20.org>]

Quick facts:
- all of these examples work with emacs -Q
- reproduces in Emacs 29, 27, 25

I think I have found that the calc package's "root" finding command does
not work as documented with calc's probability distribution functions.
Or, at least, not with "utpn", the probability distribution function for
the normal (gaussian) distribution.  I have found similar issues with
"utpt" (the Student's t-distribution), but will focus on "utpn" here.

The Emacs Calc package info page for Probability Distribution Functions
(https://www.gnu.org/software/emacs/manual/html_node/calc/Probability-Distribution-Functions.html)
has this to say about the normal distribution function:

> The ‘utpn(x,m,s)’ function uses a normal (Gaussian) distribution with
> mean ‘m’ and standard deviation ‘s’. It is the probability that such a
> normal-distributed random variable would exceed ‘x’.

I am interested in solving `utpn(x,m,s) = p` for `p` given a fixed `m`
and `s`.  The same info page has this to say about that:

> While Calc does not provide inverses of the probability distribution
> functions, the a R command can be used to solve for the inverse. Since
> the distribution functions are monotonic, a R is guaranteed to be able
> to find a solution given any initial guess.

SUPPOSITION: the above is true, especially this: "a R is guaranteed to be
able to find a solution given any initial guess."

The Root Finding info page
(https://www.gnu.org/software/emacs/manual/html_node/calc/Root-Finding.html)
has this to say about the a R command:

> The a R (calc-find-root) [root] command finds a numerical solution (or
> root) of an equation. (This command treats inequalities the same as
> equations. If the input is any other kind of formula, it is interpreted
> as an equation of the form ‘X = 0’.)
>
> The a R command requires an initial guess on the top of the stack, and a
> formula in the second-to-top position. It prompts for a solution
> variable, which must appear in the formula. All other variables that
> appear in the formula must have assigned values, i.e., when a value is
> assigned to the solution variable and the formula is evaluated with =, it
> should evaluate to a number. Any assigned value for the solution variable
> itself is ignored and unaffected by this command.
>
> When the command completes, the initial guess is replaced on the stack by
> a vector of two numbers: The value of the solution variable that solves
> the equation, and the difference between the lefthand and righthand sides
> of the equation at that value. Ordinarily, the second number will be zero
> or very nearly zero. (Note that Calc uses a slightly higher precision
> while finding the root, and thus the second number may be slightly
> different from the value you would compute from the equation yourself.)

WHAT I HAVE TRIED

I have tried entering the formula `utpn(x, 10, 1) = 0.5` and invoking `a
R` with an initial guess of `10`.  I then solve for `x`, expecting
`x = 10`.

Note that the correct answer of `x = 10` is "obviously" correct.  I am
passing `m = 10` for the "mean" arg to `utpn`, and by definition a
random value from the normal distribution has an `0.5` chance of being
larger than the mean.

To do this I type "C-x * *" to enter the *Calculator* buffer and then
type:

    ` utpn(x, 10, 1) = 0.5 RET
    10 RET                      <-- initial guess
    a R RET                     <-- invoke "root" function
    x RET                       <-- answer to prompt; solve for "x"

and I end up with this in the `*Calculator*` buffer:

    2:  utpn(x, 10, 1) = 0.5
    1:  [10., 0.]

This is correct.  The 'a R' command has returned a vector whose first
element is `x` and the second is an error delta (not important here).

Alternatively, evaluate this lisp to do the same thing:

    ELISP> (calc-eval '("root(utpn(x, 10, 1) = 0.5, x, 10)"))
    "[10., 0.]"

While correct, here I have *given* the answer to calc in the initial
guess (the root function's third arg).  The problem is that the function
doesn't do well when given other guesses.

```
` utpn(x,10,1)=0.5 RET
1.0 RET                  <---- a different guess
a R RET
x RET
```

results in this:

```
2:  utpn(x, 10, 1) = 0.5
1:  root(utpn(x, 10, 1) = 0.5, x, 1)
```

With this in the minibuffer:

> *Newton's method failed to converge: 4.8639203627839e17*

I'm finding that if the initial guess is close enough to the right answer,
in this case `10`, then `a R` works, otherwise not.  E.g. guessing `8.6`
fails but `8.7` works.

This behavior makes `a R` much less useful.  Also, documentation
suggests that any guess is expected to work (re-quoting from above with
my emphasis):

> Since the distribution functions are monotonic, 'a R' is guaranteed to
> be able to find a solution *given any initial guess*.

Have I found a bug here, or am I doing something wrong?

ADDENDUM

This shows the function failing to find a root for guesses not all that
far away from the actual answer:

ELISP> (cl-loop for x from 8.0 to 12.0 by 0.1 collect (list (format "%.1f" x) (calc-eval '("root(utpn(x,10,1)=0.5,x,$)") nil (prin1-to-string x))))
(("8.0" "root(utpn(x, 10, 1) = 0.5, x, 8.)")
 ("8.1" "root(utpn(x, 10, 1) = 0.5, x, 8.1)")
 ("8.2" "root(utpn(x, 10, 1) = 0.5, x, 8.2)")
 ("8.3" "root(utpn(x, 10, 1) = 0.5, x, 8.3)")
 ("8.4" "root(utpn(x, 10, 1) = 0.5, x, 8.4)")
 ("8.5" "root(utpn(x, 10, 1) = 0.5, x, 8.5)")
 ("8.6" "root(utpn(x, 10, 1) = 0.5, x, 8.6)")
 ("8.7" "[10., 0.]")
 ("8.8" "[10., 0.]")
 ("8.9" "[10., 0.]")
 ("9.0" "[10., 0.]")
 ("9.1" "[10., 0.]")
 ("9.2" "[10., 0.]")
 ("9.3" "[10., 0.]")
 ("9.4" "[10., 0.]")
 ("9.5" "[10., 0.]")
 ("9.6" "[10., 0.]")
 ("9.7" "[10., 0.]")
 ("9.8" "[10., 0.]")
 ("9.9" "[10., 0.]")
 ("10.0" "[10., 0.]")
 ("10.1" "[10., 0.]")
 ("10.2" "[10., 0.]")
 ("10.3" "[10., 0.]")
 ("10.4" "[10., 0.]")
 ("10.5" "[10., 0.]")
 ("10.6" "[10., 0.]")
 ("10.7" "[10., 0.]")
 ("10.8" "[10., 0.]")
 ("10.9" "[10., 0.]")
 ("11.0" "[10., 0.]")
 ("11.1" "[10., 0.]")
 ("11.2" "[10., 0.]")
 ("11.3" "[10., 0.]")
 ("11.4" "root(utpn(x, 10, 1) = 0.5, x, 11.4)")
 ("11.5" "root(utpn(x, 10, 1) = 0.5, x, 11.5)")
 ("11.6" "root(utpn(x, 10, 1) = 0.5, x, 11.6)")
 ("11.7" "root(utpn(x, 10, 1) = 0.5, x, 11.7)")
 ("11.8" "root(utpn(x, 10, 1) = 0.5, x, 11.8)")
 ("11.9" "root(utpn(x, 10, 1) = 0.5, x, 11.9)")
 ("12.0" "root(utpn(x, 10, 1) = 0.5, x, 12.)"))



In GNU Emacs 29.0.50 (build 3, x86_64-pc-linux-gnu, GTK+ Version
 3.24.34, cairo version 1.16.0) of 2022-10-30 built on naz
Repository revision: 73953b23aad5af80275838cd3b1c88e045e1856e
Repository branch: master
Windowing system distributor 'The X.Org Foundation', version 11.0.12201003
System Description: Debian GNU/Linux bookworm/sid

Configured using:
 'configure 'CFLAGS=-O2 -g3''

Configured features:
ACL CAIRO DBUS FREETYPE GIF GLIB GMP GNUTLS GPM GSETTINGS HARFBUZZ JPEG
JSON LCMS2 LIBOTF LIBSELINUX LIBSYSTEMD LIBXML2 M17N_FLT MODULES NOTIFY
INOTIFY PDUMPER PNG RSVG SECCOMP SOUND SQLITE3 THREADS TIFF
TOOLKIT_SCROLL_BARS X11 XDBE XIM XINPUT2 XPM GTK3 ZLIB

Important settings:
  value of $LANG: en_US.UTF-8
  value of $XMODIFIERS: @im=ibus
  locale-coding-system: utf-8-unix

Major mode: Text

Minor modes in effect:
  msb-mode: t
  display-time-mode: t
  flyspell-mode: t
  global-tab-line-mode: t
  tab-line-mode: t
  envrc-global-mode: t
  envrc-mode: t
  shell-dirtrack-mode: t
  auto-insert-mode: t
  keyfreq-autosave-mode: t
  keyfreq-mode: t
  savehist-mode: t
  icomplete-vertical-mode: t
  icomplete-mode: t
  editorconfig-mode: t
  which-key-mode: t
  electric-pair-mode: t
  override-global-mode: t
  tooltip-mode: t
  global-eldoc-mode: t
  show-paren-mode: t
  electric-indent-mode: t
  mouse-wheel-mode: t
  tab-bar-mode: t
  menu-bar-mode: t
  file-name-shadow-mode: t
  global-font-lock-mode: t
  font-lock-mode: t
  blink-cursor-mode: t
  column-number-mode: t
  line-number-mode: t
  auto-fill-function: do-auto-fill
  transient-mark-mode: t
  auto-composition-mode: t
  auto-encryption-mode: t
  auto-compression-mode: t
  temp-buffer-resize-mode: t
  abbrev-mode: t

Load-path shadows:
~/env/elisp/ol-notmuch hides /home/matt/.config/emacs/elpa/ol-notmuch-20220428.1337/ol-notmuch
/home/matt/.config/emacs/elpa/transient-20220918.2101/transient hides /home/matt/git/e/daily-driver/lisp/transient
/home/matt/.config/emacs/elpa/eglot-20221005.1955/eglot hides /home/matt/git/e/daily-driver/lisp/progmodes/eglot

Features:
(shadow emacsbug ielm calc-yank shortdoc calc-undo edebug dired-aux
calcalg3 imenu cl-print calc-stuff calc-funcs dabbrev calcalg2 calc-poly
calc-vec calc-mtx calc-lang calc-cplx calccomp calc-menu calc-frac
calc-comb calc-aent calc-math calc-arith calc-alg calc-bin calc-tests
calc-misc calc-prog calc-forms calc-units calc-ext calc calc-loaddefs
rect calc-macs ert pp ewoc debug backtrace mule-util magit-base crm
misearch multi-isearch copyright help-fns radix-tree pcase smerge-mode
diff qp sort smiley gnus-cite mail-extr textsec uni-scripts idna-mapping
ucs-normalize uni-confusable textsec-check gnus-async gnus-bcklg
gnus-agent gnus-srvr gnus-score score-mode nnvirtual nntp gnus-ml
gnus-msg nndoc gnus-cache gnus-dup mm-archive network-stream url-cache
display-line-numbers debbugs-gnu add-log debbugs-compat debbugs
soap-client url-http url-auth url-gw nsm rng-xsd rng-dt rng-util
xsd-regexp debbugs-browse vc-git diff-mode vc-dispatcher bug-reference
editorconfig-core editorconfig-core-handle editorconfig-fnmatch protbuf
msb time flyspell ispell org-element avl-tree generator ol-w3m ol-rmail
ol-mhe ol-irc ol-info org-habit org-agenda org-refile ol-gnus nnselect
gnus-art mm-uu mml2015 mm-view mml-smime smime gnutls dig gnus-sum
gnus-group gnus-undo gnus-start gnus-dbus dbus gnus-cloud nnimap nnmail
mail-source utf7 nnoo parse-time gnus-spec gnus-int gnus-range message
sendmail yank-media rfc822 mml mml-sec epa derived epg rfc6068
epg-config mm-decode mm-bodies mm-encode mail-parse rfc2231 rfc2047
rfc2045 ietf-drums mailabbrev gmm-utils mailheader gnus-win ol-eww eww
xdg url-queue shr pixel-fill kinsoku url-file svg xml dom puny mm-url
gnus nnheader gnus-util text-property-search mail-utils range wid-edit
mm-util mail-prsvr ol-doi org-link-doi ol-docview doc-view filenotify
jka-compr image-mode exif dired dired-loaddefs ol-bibtex ol-bbdb
tab-line server envrc inheritenv web-mode disp-table nix-mode ffap
thingatpt smie nix-repl nix-shell nix-store magit-section dash compat-27
compat-26 nix-instantiate nix-shebang nix-format nix dirtrack ob-shell
shell ob-ruby ob-python python compat compat-macs ob-dot org-protocol
org ob ob-tangle ob-ref ob-lob ob-table ob-exp org-macro org-footnote
org-src ob-comint org-pcomplete pcomplete comint ansi-osc ansi-color
ring org-list org-faces org-entities noutline outline org-version
ob-emacs-lisp ob-core ob-eval org-table oc-basic bibtex iso8601
time-date org-keys oc org-loaddefs find-func cal-menu calendar
cal-loaddefs finder-inf ol-notmuch ol rx org-compat org-macs format-spec
skeleton autoinsert advice keyfreq project edmacro kmacro warnings icons
savehist icomplete editorconfig which-key package browse-url url
url-proxy url-privacy url-expand url-methods url-history url-cookie
generate-lisp-file url-domsuf url-util mailcap url-handlers url-parse
auth-source eieio eieio-core password-cache json subr-x map byte-opt
url-vars cl-extra help-mode cl-macs gv cl-seq elec-pair use-package
use-package-ensure use-package-delight use-package-diminish
use-package-bind-key bind-key easy-mmode use-package-core cl-loaddefs
cl-lib bytecomp byte-compile info bazel-autoloads
clang-format+-autoloads clang-format-autoloads cmake-mode-autoloads
d-mode-autoloads debbugs-autoloads editorconfig-autoloads
eglot-autoloads elpy-autoloads company-autoloads envrc-autoloads
exec-path-from-shell-autoloads flymake-ruby-autoloads
flymake-easy-autoloads flymake-yamllint-autoloads go-mode-autoloads
google-c-style-autoloads graphviz-dot-mode-autoloads
highlight-indentation-autoloads inheritenv-autoloads magit-autoloads
git-commit-autoloads markdown-mode-autoloads meson-mode-autoloads
nix-mode-autoloads magit-section-autoloads dash-autoloads
nixpkgs-fmt-autoloads ol-notmuch-autoloads notmuch-autoloads
orderless-autoloads org-drill-autoloads ox-hugo-autoloads
persist-autoloads pylint-autoloads pyvenv-autoloads s-autoloads
shfmt-autoloads reformatter-autoloads tomelr-autoloads
transient-autoloads use-package-autoloads bind-key-autoloads
vertico-autoloads web-mode-autoloads which-key-autoloads
with-editor-autoloads compat-autoloads yaml-mode-autoloads
yasnippet-autoloads rmc iso-transl tooltip cconv eldoc paren electric
uniquify ediff-hook vc-hooks lisp-float-type elisp-mode mwheel
term/x-win x-win term/common-win x-dnd tool-bar dnd fontset image
regexp-opt fringe tabulated-list replace newcomment text-mode lisp-mode
prog-mode register page tab-bar menu-bar rfn-eshadow isearch easymenu
timer select scroll-bar mouse jit-lock font-lock syntax font-core
term/tty-colors frame minibuffer nadvice seq simple cl-generic
indonesian philippine cham georgian utf-8-lang misc-lang vietnamese
tibetan thai tai-viet lao korean japanese eucjp-ms cp51932 hebrew greek
romanian slovak czech european ethiopic indian cyrillic chinese
composite emoji-zwj charscript charprop case-table epa-hook
jka-cmpr-hook help abbrev obarray oclosure cl-preloaded button loaddefs
theme-loaddefs faces cus-face macroexp files window text-properties
overlay sha1 md5 base64 format env code-pages mule custom widget keymap
hashtable-print-readable backquote threads dbusbind inotify lcms2
dynamic-setting system-font-setting font-render-setting cairo
move-toolbar gtk x-toolkit xinput2 x multi-tty make-network-process
emacs)

Memory information:
((conses 16 977488 266858)
 (symbols 48 43332 287)
 (strings 32 225511 12569)
 (string-bytes 1 5996788)
 (vectors 16 103426)
 (vector-slots 8 1994918 188498)
 (floats 8 517 940)
 (intervals 56 9190 882)
 (buffers 984 38))