Maybe a bug in addplot3

[XiangyunHuang]

The minimal example

\documentclass{standalone}
\usepackage{tikz}
\usepackage{pgfplots}
\pgfplotsset{width=13cm,compat=1.18}
\begin{document}

\begin{tikzpicture}
\begin{axis}
[ 
  view/h=-30,
  colormap/viridis,
  samples=20,
  axis equal image,
  z buffer=sort,
  opacity=0.7,
]
\addplot3 [surf,shader=interp,domain=-1:0] (
 {x},
 {y},
 {x*y + sqrt(x^2*y^2 - x^2 - y^2 + 1)}
);
\end{axis}
\end{tikzpicture}

\end{document}

The output is

[XiangyunHuang]

or

\documentclass{standalone}
\usepackage{tikz}
\usepackage{pgfplots}
\pgfplotsset{width=13cm,compat=1.18}
\begin{document}

\begin{tikzpicture}
\begin{axis}
[ 
  view/h=-30,
  colormap/viridis,
  samples=20,
  axis equal image,
  % z buffer=sort,
  opacity=0.7,
]
\addplot3 [surf,shader=interp,domain=-1:1] (
 {x},
 {y},
 % {x*y}
 {x*y + sqrt(x^2*y^2 - x^2 - y^2 + 1)}
);
\end{axis}
\end{tikzpicture}

\end{document}

[muzimuzhi]

Compile the first example and I find

NOTE: coordinate (2Y1.0e0],2Y2.6318e-1],3Y0.0e0]) has been dropped because it i
s unbounded (in z). (see also unbounded coords=jump). 

in log. Here the 3Y0.0e0 stands for NaN (not a number) in the underlying fpu library of pgf/tikz. See https://tikz.dev/library-fpu#pgf.back/pgfmathfloatparsenumber.

The NaN disappears if I rewrite the function for z-coordinate from x*y + sqrt(x^2*y^2 - x^2 - y^2 + 1) to x*y + sqrt((x^2-1)*(y^2-1)), a mathematically equivalent form. Thus I infer it's a fpu problem.

Full example: z-coord function in two forms

\documentclass{standalone}
\usepackage{pgfplots}
\pgfplotsset{width=13cm,compat=1.18}
\begin{document}

\begin{tikzpicture}
\begin{axis}
[ 
  view/h=-30,
  colormap/viridis,
  samples=20,
  axis equal image,
  z buffer=sort,
  opacity=0.7,
]
\addplot3 [surf,shader=interp,domain=-1:0] (
  {x},
  {y},
  {x*y + sqrt(x^2*y^2 - x^2 - y^2 + 1)}
);
\end{axis}
\end{tikzpicture}

\begin{tikzpicture}
\begin{axis}
[ 
  view/h=-30,
  colormap/viridis,
  samples=20,
  axis equal image,
  z buffer=sort,
  opacity=0.7,
]
\addplot3 [surf,shader=interp,domain=-1:0] (
  {x},
  {y},
  {x*y + sqrt((x^2-1)*(y^2-1))}
);
\end{axis}
\end{tikzpicture}

\end{document}

[muzimuzhi]

For x = -1.0, y = -0.26138, x^2*y^2 - x^2 - y^2 + 1 should evaluate to zero, but fpu thinks it evaluates to -1.0e-7, an almost-zero negative number, but still negative. In floating-number arithmetic, square root of a negative is NaN and for all x, x + NaN = NaN, therefore the final result is NaN.

l3kernel's l3fp works well in this case.

Example 2: Comparison of step-by-step results

\documentclass{article}
\usepackage{pgf}
\usepgflibrary{fpu}

\leftskip=1in
\parindent=0pt

\def\tests{%
  \testfpu{           x^2}%
  \testfpu{               y^2}%
  \testfpu{           x^2*y^2}%
  \testfpu{           x^2*y^2 - x^2}%
  \testfpu{           x^2*y^2 - x^2 - y^2}%
  \testfpu{           x^2*y^2 - x^2 - y^2 + 1}%
  \testfpu{      sqrt(x^2*y^2 - x^2 - y^2 + 1)}%
  \testfpu{x*y + sqrt(x^2*y^2 - x^2 - y^2 + 1)}%
  \testfpu{x*y + sqrt((x^2-1)*(y^2-1))}%
}

\begin{document}

% https://github.com/pgf-tikz/pgfplots/issues/461
% NOTE: coordinate (2Y1.0e0],2Y2.6318e-1],3Y0.0e0]) has been dropped because it
% is unbounded (in z). (see also unbounded coords=jump).

\wlog{x = -1.0, y = -0.26138}

\pgfmathdeclarefunction{x}{0}{\def\pgfmathresult{-1.0}}
\pgfmathdeclarefunction{y}{0}{\def\pgfmathresult{-0.26138}}

\def\testfpu#1{%
  \pgfmathparse{#1}%
  \leavevmode
  % \wlog{\detokenize{#1} = \pgfmathresult}%
  \texttt{\llap{\detokenize{#1}} = \pgfmathresult}\par
}

\subsection*{/pgf/fpu=false}
\pgfset{fpu=false}
\tests

\subsection*{/pgf/fpu=true}
\pgfset{fpu=true, fpu/output format=fixed}
\tests

\subsection*{l3fp in l3kernel}
\ExplSyntaxOn
\fp_new_variable:n { x }
\fp_new_variable:n { y }
\fp_set_variable:nn { x } { -1.0 }
\fp_set_variable:nn { y } { -0.26138 }

\cs_set_protected:Npn \testfpu #1
  {
    \mode_leave_vertical:
    \texttt{\llap{\tl_to_str:n {#1}} ~=~ \fp_eval:n {#1}} \par
    % \iow_log:e { \tl_to_str:n {#1} ~=~ \fp_eval:n {#1} }
  }

\tests
\ExplSyntaxOff
\end{document}

[muzimuzhi]

l3kernel's l3fp works well in this case.

This only means l3fp gives more accurate results than pgf's fpu in some cases.

But in general there're essentially cases that the minor rounding errors of floating-point arithmetic will cause unexpected wrong results. So it's basically user's responsibility to guard unbounded results, for example use

x*y + sqrt(max(0, x^2*y^2 - x^2 - y^2 + 1))

as the z-coordinate plot function in original examples.

[XiangyunHuang]

Thanks for your help, I will close the issue.