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

[Ganesh Ajjanagadde]

On Sun, Nov 15, 2015 at 11:34 AM, Michael Niedermayer
<michael at niedermayer.cc> wrote:
> 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

Yes, I was extremely surprised myself, and you are right that it
defeats the NaN's purpose. Some day I may dig up the committee's
rationale for this, am curious. I doubt it was oversight, since they
are usually very careful about such things, and the defined behavior
is very specific suggesting deep thought.

Anyway, I do not take this as a formal ack yet. Hopefully we don't run
into Carl's weird debian thing that forced disabling fmax, fmin
emulation.

>
> [...]
> --
> Michael     GnuPG fingerprint: 9FF2128B147EF6730BADF133611EC787040B0FAB
>
> Why not whip the teacher when the pupil misbehaves? -- Diogenes of Sinope
>
> _______________________________________________
> ffmpeg-devel mailing list
> ffmpeg-devel at ffmpeg.org
> http://ffmpeg.org/mailman/listinfo/ffmpeg-devel
>