[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 02.05.21 um 11:08 hat Kevin Wolf geschrieben:
> 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... ;-)

Oh, äh, hm.

Wir könnten auch einfach copysign(pow(fabs(x), 1. / 3.), x).
Oh, wir haben kein copysign(). Und ich hab gar keins zum Kopieren. (Und
auch ein bisschen schade, dass unser fabs() nicht einfach "fabs" ist.)

Wie macht man copysign() überhaupt? Anscheinend gibts keine Instruktion
dafür, das wär auch zu schön gewesen. Dann könnte es wohl das billige (y
< 0. ? -fabs(x) : fabs(x)) werden; oder richtig wäre vermutlich "fstp;
fstp", dann tatsächliches Kopieren des Bits, und "fld".

> 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.

Irgendwie ist in long double rechnen, wenn man die float-Funktion
aufruft, doch auch doof.

Besser weiß ichs leider auch nicht.

>> 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?

Tja, gute Frage. O:)

Ja, einfach fyl2x; f2xm1; fld1; faddp sollten tun. (Hab ich vielleicht
einfach mit einem reinen log() angefangen?)

Jedenfalls erklärt das, was unser pow() bei negativem x tut, nämlich
eine floating-point invalid-operation auslösen...

Hanna

>> +        "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