[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Re: [PATCH 2/2] libm: pow{f,l}(), cbrt{,f,l}() und hypot{,f,l}()

Am 12.04.2021 um 23:17 hat Hanna Reitz geschrieben:
> + float- und long-double-Versionen von pow(), alle Versionen von cbrt(),
>   und alle Versionen von hypot().  Alles von µxoµcota kopiert, mit
>   geänderter Lizenz (von MIT auf BSD-2).
> 
> Signed-off-by: Hanna Reitz <hanna.reitz@xxxxxxxxxxx>
> ---
>  src/modules/include/math.h           | 10 +++++++
>  src/modules/lib/stdlibc/math/cbrt.c  | 44 ++++++++++++++++++++++++++++
>  src/modules/lib/stdlibc/math/hypot.c | 44 ++++++++++++++++++++++++++++
>  src/modules/lib/stdlibc/math/pow.c   | 31 +++++++++++++++++++-
>  4 files changed, 128 insertions(+), 1 deletion(-)
>  create mode 100644 src/modules/lib/stdlibc/math/cbrt.c
>  create mode 100644 src/modules/lib/stdlibc/math/hypot.c
> 
> diff --git a/src/modules/include/math.h b/src/modules/include/math.h
> index 8857b418..eef5f2b2 100644
> --- a/src/modules/include/math.h
> +++ b/src/modules/include/math.h
> @@ -141,6 +141,8 @@ long double nearbyintl(long double x);
>  
>  double      rint(double x);
>  double      pow(double x, double y);
> +float       powf(float x, float y);
> +long double powl(long double x, long double y);
>  
>  long        lround(double x);
>  long        lroundf(float x);
> @@ -168,6 +170,14 @@ float       modff(float x, float* i);
>  double      modf(double x, double* i);
>  long double modfl(long double x, long double* i);
>  
> +double      cbrt(double x);
> +float       cbrtf(float x);
> +long double cbrtl(long double x);
> +
> +double      hypot(double x, double y);
> +float       hypotf(float x, float y);
> +long double hypotl(long double x, long double y);
> +
>  #ifdef __cplusplus
>  }; // extern "C"
>  #endif
> diff --git a/src/modules/lib/stdlibc/math/cbrt.c b/src/modules/lib/stdlibc/math/cbrt.c
> new file mode 100644
> index 00000000..09e32977
> --- /dev/null
> +++ b/src/modules/lib/stdlibc/math/cbrt.c
> @@ -0,0 +1,44 @@
> +/*
> + * Copyright (c) 2009 The tyndur Project. All rights reserved.
> + *
> + * This code is derived from software contributed to the tyndur Project
> + * by Hanna Reitz.
> + *
> + * Redistribution and use in source and binary forms, with or without
> + * modification, are permitted provided that the following conditions
> + * are met:
> + * 1. Redistributions of source code must retain the above copyright
> + *    notice, this list of conditions and the following disclaimer.
> + * 2. Redistributions in binary form must reproduce the above copyright
> + *    notice, this list of conditions and the following disclaimer in the
> + *    documentation and/or other materials provided with the distribution.
> + *
> + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
> + * ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
> + * TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
> + * PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE COPYRIGHT HOLDERS OR
> + * CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
> + * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
> + * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS;
> + * OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,
> + * WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
> + * OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
> + * ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
> + */
> +
> +#include <math.h>
> +
> +double cbrt(double x)
> +{
> +    return pow(x, 1. / 3.);
> +}
> +
> +float cbrtf(float x)
> +{
> +    return powf(x, 1.f / 3.f);
> +}
> +
> +long double cbrtl(long double x)
> +{
> +    return powl(x, 1.L / 3.L);
> +}

RATIONALE
       For some applications, a true cube root function, which returns
       negative results for negative arguments, is more appropriate than
       pow(x, 1.0/3.0), which returns a NaN for x less than 0.

Hm, okay, überlesen wir das mal ganz schnell... ;-)

Also, nicht dass ich wüsste, was unser pow() da tut.

> diff --git a/src/modules/lib/stdlibc/math/hypot.c b/src/modules/lib/stdlibc/math/hypot.c
> new file mode 100644
> index 00000000..a78f8798
> --- /dev/null
> +++ b/src/modules/lib/stdlibc/math/hypot.c
> @@ -0,0 +1,44 @@
> +/*
> + * Copyright (c) 2009 The tyndur Project. All rights reserved.
> + *
> + * This code is derived from software contributed to the tyndur Project
> + * by Hanna Reitz.
> + *
> + * Redistribution and use in source and binary forms, with or without
> + * modification, are permitted provided that the following conditions
> + * are met:
> + * 1. Redistributions of source code must retain the above copyright
> + *    notice, this list of conditions and the following disclaimer.
> + * 2. Redistributions in binary form must reproduce the above copyright
> + *    notice, this list of conditions and the following disclaimer in the
> + *    documentation and/or other materials provided with the distribution.
> + *
> + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
> + * ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
> + * TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
> + * PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE COPYRIGHT HOLDERS OR
> + * CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
> + * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
> + * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS;
> + * OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,
> + * WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
> + * OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
> + * ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
> + */
> +
> +#include <math.h>
> +
> +double hypot(double x, double y)
> +{
> +    return sqrt(x * x + y * y);
> +}
> +
> +float hypotf(float x, float y)
> +{
> +    return sqrtf(x * x + y * y);
> +}
> +
> +long double hypotl(long double x, long double y)
> +{
> +    return sqrtl(x * x + y * y);
> +}

APPLICATION USAGE
       [...]
       These functions take precautions against overflow during intermediate steps of the computation.

Ich weiß auf Anhieb nicht, welche Precautions man da taken soll, vor
allem bei der long-double-Version. Bei den anderen könnte man vielleicht
auch einfach in long double rechnen. Aber ich nehme mal an, diese
Implementierung ist zumindest besser als gar keine.

> diff --git a/src/modules/lib/stdlibc/math/pow.c b/src/modules/lib/stdlibc/math/pow.c
> index b1f791ae..96856ffb 100644
> --- a/src/modules/lib/stdlibc/math/pow.c
> +++ b/src/modules/lib/stdlibc/math/pow.c
> @@ -28,8 +28,37 @@
>  
>  #include <math.h>
>  
> +long double powl(long double x, long double y)
> +{
> +    long double res;
> +
> +    __asm__ __volatile__ (
> +        // ST(0) == x; ST(1) == y
> +        "fld1;"
> +        // ST(0) == 1; ST(1) == x; ST(2) == y
> +        "fxch;"
> +        // ST(0) == x; ST(1) == 1; ST(2) == y
> +        "fyl2x;"
> +        // ST(0) == log2(x); ST(1) == y
> +        "fmulp;"
> +        // ST(0) == y * log2(x)

Wieso packen wir denn aufwendig eine 1 nach ST(1), damit die
Multiplikation von fyl2x nichts macht, wenn wir direkt danach manuell
multiplizieren?

Würde nicht das fyl2x allein ohne den ganzen Rest genau das tun, was man
haben will?

> +        "f2xm1;"
> +        // ST(0) == 2^(y * log2(x)) - 1 == x^y - 1
> +        "fld1;"
> +        // ST(0) == 1; ST(1) == x^y - 1
> +        "faddp"
> +        // ST(0) == x^y
> +        : "=t"(res) : "0"(x), "u"(y));
> +
> +    return res;
> +}

Kevin