Changeset 516e780 in mainline for uspace/lib/math


Ignore:
Timestamp:
2018-08-31T11:55:41Z (7 years ago)
Author:
GitHub <noreply@…>
Branches:
lfn, master, serial, ticket/834-toolchain-update, topic/msim-upgrade, topic/simplify-dev-export
Children:
fa86fff
Parents:
7f7d642
git-author:
Jiří Zárevúcky <zarevucky.jiri@…> (2018-08-31 11:55:41)
git-committer:
GitHub <noreply@…> (2018-08-31 11:55:41)
Message:

Strip down libmath. (#45)

libmath is mostly unused (except for trunc(), sin() and cos()), and most functions in it are either very imprecise or downright broken. Additionally, it is implemented in manner that conflicts with C standard. Instead of trying to fix all the shortcomings while maintaining unused functionality, I'm opting to simply remove most of it and only keep the parts that are currently necessary.

Later readdition of the removed functions is possible, but there needs to be a reliable way to evaluate their quality first.

Location:
uspace/lib/math
Files:
2 added
64 deleted
5 edited
5 moved

Legend:

Unmodified
Added
Removed

  • uspace/lib/math/Makefile

    r7f7d642 r516e780  
    3030ROOT_PATH = $(USPACE_PREFIX)/..
    3131
    32 CONFIG_MAKEFILE = $(ROOT_PATH)/Makefile.config
    33 
    3432LIBRARY = libmath
    3533SOVERSION = 0.0
    3634
    37 EXTRA_CFLAGS += -Iarch/$(UARCH)/include
    38 
    39 -include $(CONFIG_MAKEFILE)
    40 -include arch/$(UARCH)/Makefile.inc
    41 
    42 GENERIC_SOURCES = \
    43         generic/acos.c \
    44         generic/asin.c \
    45         generic/atan.c \
    46         generic/atan2.c \
    47         generic/ceil.c \
    48         generic/cosh.c \
    49         generic/exp.c \
     35SOURCES = \
     36        generic/__fcompare.c \
     37        generic/__fpclassify.c \
     38        generic/__signbit.c \
    5039        generic/fabs.c \
    51         generic/floor.c \
    5240        generic/fmod.c \
    53         generic/frexp.c \
    54         generic/ldexp.c \
    55         generic/log.c \
    56         generic/log10.c \
    57         generic/log2.c \
    58         generic/modf.c \
    59         generic/pow.c \
    60         generic/sinh.c \
    61         generic/sqrt.c \
    62         generic/tan.c \
    63         generic/tanh.c \
     41        generic/nearbyint.c \
     42        generic/round.c \
    6443        generic/trig.c \
    6544        generic/trunc.c
    6645
    67 SOURCES = \
    68         $(GENERIC_SOURCES) \
    69         $(ARCH_SOURCES)
     46TEST_SOURCES = \
     47        test/rounding.c \
     48        test/main.c
    7049
    7150include $(USPACE_PREFIX)/Makefile.common

  • uspace/lib/math/generic/__fcompare.c

    r7f7d642 r516e780  
    11/*
    2  * Copyright (c) 2015 Jiri Svoboda
     2 * Copyright (c) 2018 CZ.NIC, z.s.p.o.
    33 * All rights reserved.
    44 *
    … …  
    2727 */
    2828
    29 /** @addtogroup libmath
    30  * @{
     29#include <math.h>
     30#include <stdarg.h>
     31
     32/**
     33 * Fallback symbol used when code including <math.h> is compiled with something
     34 * other than GCC or Clang. The function itself must be built with GCC or Clang.
    3135 */
    32 /** @file
    33  */
     36int __fcompare(size_t sz1, size_t sz2, ...)
     37{
     38        va_list ap;
     39        va_start(ap, sz2);
    3440
    35 #include <ceil.h>
    36 #include <math.h>
    37 #include <mathtypes.h>
     41        long double val1;
     42        long double val2;
    3843
    39 /** Ceiling (round towards positive infinity, 32-bit floating point)
    40  *
    41  * @param val Floating point number.
    42  *
    43  * @return Number rounded towards positive infinity.
    44  */
    45 float32_t float32_ceil(float32_t val)
    46 {
    47         float32_u t;
    48         float32_u v;
    49         float32_u r;
    50 
    51         v.val = val;
    52         t.val = trunc_f32(val);
    53 
    54         if (v.data.parts.sign == 1 || val == t.val) {
    55                 r = t;
    56         } else {
    57                 r.val = t.val + 1.0;
     44        switch (sz1) {
     45        case 4:
     46                val1 = (long double) va_arg(ap, double);
     47                break;
     48        case 8:
     49                val1 = (long double) va_arg(ap, double);
     50                break;
     51        default:
     52                val1 = va_arg(ap, long double);
     53                break;
    5854        }
    5955
    60         return r.val;
     56        switch (sz2) {
     57        case 4:
     58                val2 = (long double) va_arg(ap, double);
     59                break;
     60        case 8:
     61                val2 = (long double) va_arg(ap, double);
     62                break;
     63        default:
     64                val2 = va_arg(ap, long double);
     65                break;
     66        }
     67
     68        va_end(ap);
     69
     70        if (isgreaterequal(val1, val2)) {
     71                if (isgreater(val1, val2))
     72                        return __FCOMPARE_GREATER;
     73                else
     74                        return __FCOMPARE_EQUAL;
     75        } else {
     76                if (isless(val1, val2))
     77                        return __FCOMPARE_LESS;
     78                else
     79                        return 0;
     80        }
    6181}
    6282
    63 /** Ceiling (round towards positive infinity, 64-bit floating point)
    64  *
    65  * @param val Floating point number.
    66  *
    67  * @return Number rounded towards positive infinity.
    68  */
    69 float64_t float64_ceil(float64_t val)
    70 {
    71         float64_u t;
    72         float64_u v;
    73         float64_u r;
    74 
    75         v.val = val;
    76         t.val = trunc_f64(val);
    77 
    78         if (v.data.parts.sign == 1 || val == t.val) {
    79                 r = t;
    80         } else {
    81                 r.val = t.val + 1.0;
    82         }
    83 
    84         return r.val;
    85 }
    86 
    87 /** @}
    88  */

  • uspace/lib/math/generic/__fpclassify.c

    r7f7d642 r516e780  
    11/*
    2  * Copyright (c) 2015 Jiri Svoboda
     2 * Copyright (c) 2018 CZ.NIC, z.s.p.o.
    33 * All rights reserved.
    44 *
    … …  
    2727 */
    2828
    29 /** @addtogroup libmath
    30  * @{
     29#include <math.h>
     30#include <stdarg.h>
     31
     32/**
     33 * Fallback symbol used when code including <math.h> is compiled with something
     34 * other than GCC or Clang. The function itself must be built with GCC or Clang.
    3135 */
    32 /** @file
    33  */
     36int __fpclassify(size_t sz, ...)
     37{
     38        va_list ap;
     39        va_start(ap, sz);
    3440
    35 #include <math.h>
    36 #include <sqrt.h>
     41        int result;
    3742
    38 /** Single precision square root
    39  *
    40  * Compute square root.
    41  *
    42  * @param val Value
    43  *
    44  * @return Square root.
    45  *
    46  */
    47 float32_t float32_sqrt(float32_t val)
    48 {
    49         return pow_f32(val, 0.5);
     43        switch (sz) {
     44        case 4:
     45                result = fpclassify((float) va_arg(ap, double));
     46                break;
     47        case 8:
     48                result = fpclassify(va_arg(ap, double));
     49                break;
     50        default:
     51                result = fpclassify(va_arg(ap, long double));
     52                break;
     53        }
     54
     55        va_end(ap);
     56        return result;
    5057}
    5158
    52 /** Double precision square root
    53  *
    54  * Compute squre root.
    55  *
    56  * @param val Value
    57  *
    58  * @return Square root.
    59  *
    60  */
    61 float64_t float64_sqrt(float64_t val)
    62 {
    63         return pow_f64(val, 0.5);
    64 }
    65 
    66 /** @}
    67  */

  • uspace/lib/math/generic/__signbit.c

    r7f7d642 r516e780  
    11/*
    2  * Copyright (c) 2015 Jiri Svoboda
     2 * Copyright (c) 2018 CZ.NIC, z.s.p.o.
    33 * All rights reserved.
    44 *
    … …  
    2727 */
    2828
    29 /** @addtogroup libmath
    30  * @{
     29#include <math.h>
     30#include <stdarg.h>
     31
     32/**
     33 * Fallback symbol used when code including <math.h> is compiled with something
     34 * other than GCC or Clang. The function itself must be built with GCC or Clang.
    3135 */
    32 /** @file
    33  */
     36int __signbit(size_t sz, ...)
     37{
     38        va_list ap;
     39        va_start(ap, sz);
    3440
    35 #include <math.h>
    36 #include <tan.h>
     41        int result;
    3742
    38 /** Tangent (32-bit floating point)
    39  *
    40  * Compute tangent value.
    41  *
    42  * @param arg Tangent argument.
    43  *
    44  * @return Tangent value.
    45  *
    46  */
    47 float32_t float32_tan(float32_t arg)
    48 {
    49         return sin_f32(arg) / cos_f32(arg);
     43        switch (sz) {
     44        case 4:
     45                result = signbit(va_arg(ap, double));
     46                break;
     47        case 8:
     48                result = signbit(va_arg(ap, double));
     49                break;
     50        default:
     51                result = signbit(va_arg(ap, long double));
     52                break;
     53        }
     54
     55        va_end(ap);
     56        return result;
    5057}
    5158
    52 /** Sine (64-bit floating point)
    53  *
    54  * Compute sine value.
    55  *
    56  * @param arg Sine argument.
    57  *
    58  * @return Sine value.
    59  *
    60  */
    61 float64_t float64_tan(float64_t arg)
    62 {
    63         return sin_f64(arg) / cos_f64(arg);
    64 }
    65 
    66 
    67 /** @}
    68  */

  • uspace/lib/math/generic/fabs.c

    r7f7d642 r516e780  
    11/*
    2  * Copyright (c) 2015 Jiri Svoboda
     2 * Copyright (c) 2014 Martin Decky
    33 * All rights reserved.
    44 *
    … …  
    3434
    3535#include <math.h>
    36 #include <fabs.h>
    3736
    38 /** Absolute value (32-bit floating point)
    39  *
    40  * Compute absolute value.
    41  *
    42  * @param arg Argument.
    43  *
    44  * @return Absolute value.
    45  *
    46  */
    47 float32_t float32_fabs(float32_t arg)
     37float fabsf(float val)
    4838{
    49         if (arg < 0.0)
    50                 return -arg;
    51         else
    52                 return arg;
     39        return copysignf(val, 1.0f);
    5340}
    5441
    55 /** Absolute value (64-bit floating point)
    56  *
    57  * Compute absolute value.
    58  *
    59  * @param arg Argument.
    60  *
    61  * @return Absolute value.
    62  *
    63  */
    64 float64_t float64_fabs(float64_t arg)
     42double fabs(double val)
    6543{
    66         if (arg < 0.0)
    67                 return -arg;
    68         else
    69                 return arg;
     44        return copysign(val, 1.0);
    7045}
    7146

  • uspace/lib/math/generic/fmod.c

    r7f7d642 r516e780  
    3333 */
    3434
    35 #include <fmod.h>
    3635#include <math.h>
    3736
    … …  
    5251 *
    5352 */
    54 float32_t float32_fmod(float32_t dividend, float32_t divisor)
     53float fmodf(float dividend, float divisor)
    5554{
    5655        // FIXME: replace with exact arithmetics
    5756
    58         float32_t quotient = trunc_f32(dividend / divisor);
     57        float quotient = truncf(dividend / divisor);
    5958
    6059        return (dividend - quotient * divisor);
    … …  
    7776 *
    7877 */
    79 float64_t float64_fmod(float64_t dividend, float64_t divisor)
     78double fmod(double dividend, double divisor)
    8079{
    8180        // FIXME: replace with exact arithmetics
    8281
    83         float64_t quotient = trunc_f64(dividend / divisor);
     82        double quotient = trunc(dividend / divisor);
    8483
    8584        return (dividend - quotient * divisor);

  • uspace/lib/math/generic/round.c

    r7f7d642 r516e780  
    3535
    3636#include <math.h>
    37 #include <pow.h>
     37#include <float.h>
     38#include <stdint.h>
    3839
    39 /** Single precision power
    40  *
    41  * Compute power value.
    42  *
    43  * @param x Base
    44  * @param y Exponent
    45  *
    46  * @return Cosine value.
    47  *
     40/**
     41 * Rounds its argument to the nearest integer value in floating-point format,
     42 * rounding halfway cases away from zero, regardless of the current rounding
     43 * direction.
    4844 */
    49 float32_t float32_pow(float32_t x, float32_t y)
     45float roundf(float val)
    5046{
    51         /* x^y = (e ^ log(x))^y = e ^ (log(x) * y) */
    52         return exp_f32(log_f32(x) * y);
     47        const int exp_bias = FLT_MAX_EXP - 1;
     48        const int mant_bits = FLT_MANT_DIG - 1;
     49
     50        union {
     51                float f;
     52                uint32_t i;
     53        } u = { .f = fabsf(val) };
     54
     55        int exp = (u.i >> mant_bits) - exp_bias;
     56
     57        /* If value is less than 0.5, return zero with appropriate sign. */
     58        if (exp < -1)
     59                return copysignf(0.0f, val);
     60
     61        /* If exponent is exactly mant_bits, adding 0.5 could change the result. */
     62        if (exp >= mant_bits)
     63                return val;
     64
     65        /* Use trunc with adjusted value to do the rounding. */
     66        return copysignf(truncf(u.f + 0.5), val);
    5367}
    5468
    55 /** Double precision power
    56  *
    57  * Compute power value.
    58  *
    59  * @param x Base
    60  * @param y Exponent
    61  *
    62  * @return Cosine value.
    63  *
     69/**
     70 * Rounds its argument to the nearest integer value in floating-point format,
     71 * rounding halfway cases away from zero, regardless of the current rounding
     72 * direction.
    6473 */
    65 float64_t float64_pow(float64_t x, float64_t y)
     74double round(double val)
    6675{
    67         /* x^y = (e ^ log(x))^y = e ^ (log(x) * y) */
    68         return exp_f64(log_f64(x) * y);
     76        const int exp_bias = DBL_MAX_EXP - 1;
     77        const int mant_bits = DBL_MANT_DIG - 1;
     78
     79        union {
     80                double f;
     81                uint64_t i;
     82        } u = { .f = fabs(val) };
     83
     84        int exp = ((int)(u.i >> mant_bits)) - exp_bias;
     85
     86        /* If value is less than 0.5, return zero with appropriate sign. */
     87        if (exp < -1)
     88                return copysign(0.0, val);
     89
     90        /* If exponent is exactly mant_bits, adding 0.5 could change the result. */
     91        if (exp >= mant_bits)
     92                return val;
     93
     94        /* Use trunc with adjusted value to do the rounding. */
     95        return copysign(trunc(u.f + 0.5), val);
    6996}
    7097

  • uspace/lib/math/generic/trig.c

    r7f7d642 r516e780  
    3535
    3636#include <math.h>
    37 #include <trig.h>
    3837
    3938#define TAYLOR_DEGREE_32 13
    … …  
    4140
    4241/** Precomputed values for factorial (starting from 1!) */
    43 static float64_t factorials[TAYLOR_DEGREE_64] = {
     42static double factorials[TAYLOR_DEGREE_64] = {
    4443        1, 2, 6, 24, 120, 720, 5040, 40320, 362880, 3628800, 39916800,
    4544        479001600, 6227020800.0L, 87178291200.0L, 1307674368000.0L,
    … …  
    6059 *
    6160 */
    62 static float32_t taylor_sin_32(float32_t arg)
    63 {
    64         float32_t ret = 0;
    65         float32_t nom = 1;
     61static float taylor_sin_32(float arg)
     62{
     63        float ret = 0;
     64        float nom = 1;
    6665
    6766        for (unsigned int i = 0; i < TAYLOR_DEGREE_32; i++) {
    … …  
    8988 *
    9089 */
    91 static float64_t taylor_sin_64(float64_t arg)
    92 {
    93         float64_t ret = 0;
    94         float64_t nom = 1;
     90static double taylor_sin_64(double arg)
     91{
     92        double ret = 0;
     93        double nom = 1;
    9594
    9695        for (unsigned int i = 0; i < TAYLOR_DEGREE_64; i++) {
    … …  
    118117 *
    119118 */
    120 static float32_t taylor_cos_32(float32_t arg)
    121 {
    122         float32_t ret = 1;
    123         float32_t nom = 1;
     119static float taylor_cos_32(float arg)
     120{
     121        float ret = 1;
     122        float nom = 1;
    124123
    125124        for (unsigned int i = 0; i < TAYLOR_DEGREE_32; i++) {
    … …  
    147146 *
    148147 */
    149 static float64_t taylor_cos_64(float64_t arg)
    150 {
    151         float64_t ret = 1;
    152         float64_t nom = 1;
     148static double taylor_cos_64(double arg)
     149{
     150        double ret = 1;
     151        double nom = 1;
    153152
    154153        for (unsigned int i = 0; i < TAYLOR_DEGREE_64; i++) {
    … …  
    176175 *
    177176 */
    178 static float32_t base_sin_32(float32_t arg)
     177static float base_sin_32(float arg)
    179178{
    180179        unsigned int period = arg / (M_PI / 4);
    … …  
    209208 *
    210209 */
    211 static float64_t base_sin_64(float64_t arg)
     210static double base_sin_64(double arg)
    212211{
    213212        unsigned int period = arg / (M_PI / 4);
    … …  
    242241 *
    243242 */
    244 static float32_t base_cos_32(float32_t arg)
     243static float base_cos_32(float arg)
    245244{
    246245        unsigned int period = arg / (M_PI / 4);
    … …  
    275274 *
    276275 */
    277 static float64_t base_cos_64(float64_t arg)
     276static double base_cos_64(double arg)
    278277{
    279278        unsigned int period = arg / (M_PI / 4);
    … …  
    305304 *
    306305 */
    307 float32_t float32_sin(float32_t arg)
    308 {
    309         float32_t base_arg = fmod_f32(arg, 2 * M_PI);
     306float sinf(float arg)
     307{
     308        float base_arg = fmodf(arg, 2 * M_PI);
    310309
    311310        if (base_arg < 0)
    … …  
    324323 *
    325324 */
    326 float64_t float64_sin(float64_t arg)
    327 {
    328         float64_t base_arg = fmod_f64(arg, 2 * M_PI);
     325double sin(double arg)
     326{
     327        double base_arg = fmod(arg, 2 * M_PI);
    329328
    330329        if (base_arg < 0)
    … …  
    343342 *
    344343 */
    345 float32_t float32_cos(float32_t arg)
    346 {
    347         float32_t base_arg = fmod_f32(arg, 2 * M_PI);
     344float cosf(float arg)
     345{
     346        float base_arg = fmodf(arg, 2 * M_PI);
    348347
    349348        if (base_arg < 0)
    … …  
    362361 *
    363362 */
    364 float64_t float64_cos(float64_t arg)
    365 {
    366         float64_t base_arg = fmod_f64(arg, 2 * M_PI);
     363double cos(double arg)
     364{
     365        double base_arg = fmod(arg, 2 * M_PI);
    367366
    368367        if (base_arg < 0)

  • uspace/lib/math/generic/trunc.c

    r7f7d642 r516e780  
    3434 */
    3535
    36 #include <mathtypes.h>
    37 #include <trunc.h>
     36#include <math.h>
     37#include <float.h>
     38#include <stdint.h>
    3839
    3940/** Truncate fractional part (round towards zero)
    … …  
    5253 *
    5354 */
    54 float32_t float32_trunc(float32_t val)
     55float truncf(float val)
    5556{
    56         float32_u v;
    57         int32_t exp;
     57        const int exp_bias = FLT_MAX_EXP - 1;
     58        const int mant_bits = FLT_MANT_DIG - 1;
     59        const uint32_t mant_mask = (UINT32_C(1) << mant_bits) - 1;
    5860
    59         v.val = val;
    60         exp = v.data.parts.exp - FLOAT32_BIAS;
     61        union {
     62                float f;
     63                uint32_t i;
     64        } u = { .f = fabsf(val) };
    6165
    62         if (exp < 0) {
    63                 /* -1 < val < 1 => result is +0 or -0 */
    64                 v.data.parts.exp = 0;
    65                 v.data.parts.fraction = 0;
    66         } else if (exp >= FLOAT32_FRACTION_SIZE) {
    67                 if (exp == 1024) {
    68                         /* val is +inf, -inf or NaN => trigger an exception */
    69                         // FIXME TODO
    70                 }
     66        int exp = (u.i >> mant_bits) - exp_bias;
    7167
    72                 /* All bits in val are relevant for the result */
    73         } else {
    74                 /* Truncate irrelevant fraction bits */
    75                 v.data.parts.fraction &= ~(UINT32_C(0x007fffff) >> exp);
    76         }
     68        /* If value is less than one, return zero with appropriate sign. */
     69        if (exp < 0)
     70                return copysignf(0.0f, val);
    7771
    78         return v.val;
     72        if (exp >= mant_bits)
     73                return val;
     74
     75        /* Truncate irrelevant fraction bits */
     76        u.i &= ~(mant_mask >> exp);
     77        return copysignf(u.f, val);
    7978}
    8079
    … …  
    9493 *
    9594 */
    96 float64_t float64_trunc(float64_t val)
     95double trunc(double val)
    9796{
    98         float64_u v;
    99         int32_t exp;
     97        const int exp_bias = DBL_MAX_EXP - 1;
     98        const int mant_bits = DBL_MANT_DIG - 1;
     99        const uint64_t mant_mask = (UINT64_C(1) << mant_bits) - 1;
    100100
    101         v.val = val;
    102         exp = v.data.parts.exp - FLOAT64_BIAS;
     101        union {
     102                double f;
     103                uint64_t i;
     104        } u = { .f = fabs(val) };
    103105
    104         if (exp < 0) {
    105                 /* -1 < val < 1 => result is +0 or -0 */
    106                 v.data.parts.exp = 0;
    107                 v.data.parts.fraction = 0;
    108         } else if (exp >= FLOAT64_FRACTION_SIZE) {
    109                 if (exp == 1024) {
    110                         /* val is +inf, -inf or NaN => trigger an exception */
    111                         // FIXME TODO
    112                 }
     106        int exp = ((int)(u.i >> mant_bits)) - exp_bias;
    113107
    114                 /* All bits in val are relevant for the result */
    115         } else {
    116                 /* Truncate irrelevant fraction bits */
    117                 v.data.parts.fraction &= ~(UINT64_C(0x000fffffffffffff) >> exp);
    118         }
     108        /* If value is less than one, return zero with appropriate sign. */
     109        if (exp < 0)
     110                return copysign(0.0, val);
    119111
    120         return v.val;
     112        if (exp >= mant_bits)
     113                return val;
     114
     115        /* Truncate irrelevant fraction bits */
     116        u.i &= ~(mant_mask >> exp);
     117        return copysign(u.f, val);
    121118}
    122119

  • uspace/lib/math/test/main.c

    r7f7d642 r516e780  
    11/*
    2  * Copyright (c) 2015 Jiri Svoboda
     2 * Copyright (c) 2014 Vojtech Horky
    33 * All rights reserved.
    44 *
    … …  
    2727 */
    2828
    29 /** @addtogroup libmath
    30  * @{
    31  */
    32 /** @file
    33  */
     29#include <stdio.h>
     30#include <pcut/pcut.h>
    3431
    35 #ifndef LIBMATH_LOG_H_
    36 #define LIBMATH_LOG_H_
     32PCUT_INIT;
    3733
    38 #include <mathtypes.h>
     34PCUT_IMPORT(rounding);
    3935
    40 extern float32_t float32_log(float32_t);
    41 extern float64_t float64_log(float64_t);
     36PCUT_MAIN();
    4237
    43 #endif
    44 
    45 /** @}
    46  */
Note: See TracChangeset for help on using the changeset viewer.