[FFmpeg-devel] [PATCHv2] configure+libm.h: add hypot emulation

[Michael Niedermayer]

On Sun, Nov 15, 2015 at 11:01:58AM -0500, Ganesh Ajjanagadde wrote:
> On Sun, Nov 15, 2015 at 10:56 AM, Michael Niedermayer
> <michael at niedermayer.cc> wrote:
> > On Sun, Nov 15, 2015 at 10:03:37AM -0500, Ganesh Ajjanagadde wrote:
> >> It is known that the naive sqrt(x*x + y*y) approach for computing the
> >> hypotenuse suffers from overflow and accuracy issues, see e.g
> >> http://www.johndcook.com/blog/2010/06/02/whats-so-hard-about-finding-a-hypotenuse/.
> >> This adds hypot support to FFmpeg, a C99 function.
> >>
> >> On platforms without hypot, this patch does a reaonable workaround, that
> >> although not as accurate as GNU libm, is readable and does not suffer
> >> from the overflow issue. Improvements can be made separately.
> >>
> >> Signed-off-by: Ganesh Ajjanagadde <gajjanagadde at gmail.com>
> >> ---
> >>  configure        |  2 ++
> >>  libavutil/libm.h | 23 +++++++++++++++++++++++
> >>  2 files changed, 25 insertions(+)
> >>
> >> diff --git a/configure b/configure
> >> index d518b21..45df724 100755
> >> --- a/configure
> >> +++ b/configure
> >> @@ -1774,6 +1774,7 @@ MATH_FUNCS="
> >>      exp2
> >>      exp2f
> >>      expf
> >> +    hypot
> >>      isinf
> >>      isnan
> >>      ldexpf
> >> @@ -5309,6 +5310,7 @@ disabled crystalhd || check_lib libcrystalhd/libcrystalhd_if.h DtsCrystalHDVersi
> >>
> >>  atan2f_args=2
> >>  copysign_args=2
> >> +hypot_args=2
> >>  ldexpf_args=2
> >>  powf_args=2
> >>
> >> diff --git a/libavutil/libm.h b/libavutil/libm.h
> >> index 6c17b28..f7a2b41 100644
> >> --- a/libavutil/libm.h
> >> +++ b/libavutil/libm.h
> >> @@ -102,6 +102,29 @@ static av_always_inline av_const int isnan(float x)
> >>  }
> >>  #endif /* HAVE_ISNAN */
> >>
> >> +#if !HAVE_HYPOT
> >> +#undef hypot
> >> +static inline av_const double hypot(double x, double y)
> >> +{
> >> +    double ret, temp;
> >> +    x = fabs(x);
> >> +    y = fabs(y);
> >> +
> >> +    if (isinf(x) || isinf(y))
> >> +        return av_int2double(0x7ff0000000000000);
> >
> > if either is NaN the result should be NaN i think
> > return x+y
> > might achive this
> 
> No, quoting the man page/standard:
>        If x or y is an infinity, positive infinity is returned.
> 
>        If x or y is a NaN, and the other argument is not an infinity,
> a NaN is returned.

indeed, the spec says thats how it should be.

this is not what i expected though and renders the function
problematic in practice IMHO.
For example a big use of NaN is to detect errors.
One does a big complicated computation and if at the end the result is
NaN or +-Inf then one knows there was something wrong in it.
if NaN is infective then any computation that returns it can reliably
be detected. These exceptions in C like hypot() break this.
1/hypot(x,y) should be NaN if either x or y was NaN

also mathematically its wrong to ignore a NaN argument
consider
hypot(sqrt(-x), sqrt(x)) for x->infinite

of course theres nothing we can or should do about hypot() its defined
in C as it is defined but its something one should be aware of if
one expects that NaNs can be used as a reliable means to detect
NaNs from intermediate steps of a complicated calculation

[...]
-- 
Michael     GnuPG fingerprint: 9FF2128B147EF6730BADF133611EC787040B0FAB

Why not whip the teacher when the pupil misbehaves? -- Diogenes of Sinope
-------------- next part --------------
A non-text attachment was scrubbed...
Name: signature.asc
Type: application/pgp-signature
Size: 181 bytes
Desc: Digital signature
URL: <http://ffmpeg.org/pipermail/ffmpeg-devel/attachments/20151115/957241de/attachment.sig>