Changeset 516e780 in mainline


Ignore:
Timestamp:
2018-08-31T11:55:41Z (6 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.

Files:
3 added
64 deleted
12 edited
8 moved

Legend:

Unmodified
Added
Removed

  • boot/Makefile.common

    r7f7d642 r516e780  
    244244        $(USPACE_PATH)/lib/sif/test-libsif \
    245245        $(USPACE_PATH)/lib/uri/test-liburi \
     246        $(USPACE_PATH)/lib/math/test-libmath \
    246247        $(USPACE_PATH)/drv/bus/usb/xhci/test-xhci \
    247248        $(USPACE_PATH)/app/bdsh/test-bdsh \

  • uspace/Makefile.common

    r7f7d642 r516e780  
    109109LIBSOFTINT_PREFIX = $(LIB_PREFIX)/softint
    110110
    111 LIBMATH_PREFIX = $(LIB_PREFIX)/math
    112 LIBMATH_INCLUDES_FLAGS = \
    113         -I$(LIBMATH_PREFIX)/include \
    114         -I$(LIBMATH_PREFIX)/arch/$(UARCH)/include
    115 
    116 LIBPOSIX_PREFIX = $(LIB_PREFIX)/posix
    117111LIBDLTEST_PREFIX = $(LIB_PREFIX)/dltest
    118112

  • uspace/app/tester/float/float2.c

    r7f7d642 r516e780  
    4343};
    4444
    45 static double arguments_acos[OPERANDS] = {
    46         -0.936456687291, -0.504846104600, 0.862318872288, 0.964966028492,
    47         0.987353618220, 1.0, -0.194939922623, 0.978471923925, -0.999023478833,
    48         0.540302305868
    49 };
    50 
    51 static double arguments_asin[OPERANDS] = {
    52         -0.350783227690, -0.863209366649, -0.506365641110, -0.262374853704,
    53         0.158533380044, 0.0, 0.980815184715, -0.206379975025, -0.044182448332,
    54         0.841470984808
    55 };
    56 
    57 static double arguments_atan[OPERANDS] = {
    58         3.5, 100.0, 50.0, 768.3156, 1080.499999, 1.0, 66.0,
    59         2.718281828459045, 9.9, 0.001
    60 };
    61 
    62 static double arguments_exp[OPERANDS] = {
    63         3.5, -2.1, 50.0, 0.0, 1.0, 13.2, -1.1, -5.5, 0.1, -66.0
    64 };
    65 
    66 static double arguments_log[OPERANDS] = {
    67         3.5, 100.0, 50.0, 768.3156, 1080.499999, 1.0, 66.0,
    68         2.718281828459045, 9.9, 0.001
    69 };
    70 
    71 static double arguments_sqrt[OPERANDS] = {
    72         3.5, 100.0, 50.0, 768.3156, 1080.499999, 1.0, 66.0,
    73         2.718281828459045, 9.9, 0.001
    74 };
    75 
    76 static double arguments_tanh[OPERANDS] = {
    77         3.5, -2.1, 50.0, 0.0, 1.0, 13.2, -1.1, -5.5, 0.000001, -66000000.0
    78 };
    79 
    80 static double results_acos[OPERANDS] = {
    81         2.783185307180, 2.100000000000, 0.530964914873, 0.265482457437,
    82         0.159205070272, 0.000000000000, 1.766992524091, 0.207873834887,
    83         3.097395817941, 1.000000000000
    84 };
    85 
    86 static double results_asin[OPERANDS] = {
    87         -0.358407346411, -1.041592653590, -0.530964914874, -0.265482457437,
    88         0.159205070273, 0.000000000000, 1.374600129498, -0.207873834889,
    89         -0.044196835651, 1.000000000000
    90 };
    91 
    92 static double results_atan[OPERANDS] = {
    93         1.292496667790, 1.560796660108, 1.550798992822, 1.569494779052,
    94         1.569870829603, 0.785398163397, 1.555645970920, 1.218282905017,
    95         1.470127674637, 0.000999999667
    96 };
    97 
    98 static double results_ceil[OPERANDS] = {
    99         4.0, -2.0, 100.0, 50.0, -1024.0, 0.0, 769.0, 1081.0, -600.0, 1.0
    100 };
    101 
    10245static double results_cos[OPERANDS] = {
    10346        -0.936456687291, -0.504846104600, 0.862318872288, 0.964966028492,
    … …  
    10649};
    10750
    108 static double results_cosh[OPERANDS] = {
    109         16.572824671057, 4.144313170410, 2592352764293536022528.000000000000,
    110         1.000000000000, 1.543080634815, 270182.468624271103, 1.668518553822,
    111         122.348009517829, 1.005004168056, 23035933171656458903220125696.0
    112 };
    113 
    114 static double results_fabs[OPERANDS] = {
    115         3.5, 2.1, 100.0, 50.0, 1024.0, 0.0, 768.3156, 1080.499999, 600.0, 1.0
    116 };
    117 
    118 static double results_floor[OPERANDS] = {
    119         3.0, -3.0, 100.0, 50.0, -1024.0, 0.0, 768.0, 1080.0, -600.0, 1.0
    120 };
    121 
    122 static double results_exp[OPERANDS] = {
    123         33.115451958692, 0.122456428253, 5184705528587072045056.0,
    124         1.000000000000, 2.718281828459, 540364.937246691552, 0.332871083698,
    125         0.004086771438, 1.105170918076, 0.000000000000
    126 };
    127 
    128 static double results_log[OPERANDS] = {
    129         1.252762968495, 4.605170185988, 3.912023005428, 6.644200586236,
    130         6.985179175021, 0.000000000000, 4.189654742026, 1.000000000000,
    131         2.292534757141, -6.907755278982
    132 };
    133 
    134 static double results_log10[OPERANDS] = {
    135         0.544068044350, 2.000000000000, 1.698970004336, 2.885539651261,
    136         3.033624770817, 0.000000000000, 1.819543935542, 0.434294481903,
    137         0.995635194598, -3.000000000000
    138 };
    139 
    140 static double results_log2[OPERANDS] = {
    141         1.807354922058, 6.643856189775, 5.643856189775, 9.585555236434,
    142         10.077483355524, 0.000000000000, 6.044394119358, 1.442695040889,
    143         3.307428525192, -9.965784284662
    144 };
    145 
    14651static double results_sin[OPERANDS] = {
    14752        -0.350783227690, -0.863209366649, -0.506365641110, -0.262374853704,
    14853        0.158533380044, 0.0, 0.980815184715, -0.206379975025, -0.044182448332,
    14954        0.841470984808
    150 };
    151 
    152 static double results_sinh[OPERANDS] = {
    153         16.542627287635, -4.021856742157, 2592352764293536022528.000000000000,
    154         0.000000000000, 1.175201193644, 270182.468622420449, -1.335647470124,
    155         -122.343922746391, 0.100166750020, -23035933171656458903220125696.0
    156 };
    157 
    158 static double results_sqrt[OPERANDS] = {
    159         1.870828693387, 10.000000000000, 7.071067811865, 27.718506453271,
    160         32.870959812576, 1.000000000000, 8.124038404636, 1.648721270700,
    161         3.146426544510, 0.031622776602
    162 };
    163 
    164 static double results_tan[OPERANDS] = {
    165         0.374585640159, 1.709846542905, -0.587213915157, -0.271900611998,
    166         0.160563932839, 0.000000000000, -5.031371570891, -0.210920691722,
    167         0.044225635601, 1.557407724655
    168 };
    169 
    170 static double results_tanh[OPERANDS] = {
    171         0.998177897611, -0.970451936613, 1.000000000000, 0.000000000000,
    172         0.761594155956, 0.999999999993, -0.800499021761, -0.999966597156,
    173         0.000001000000, -1.000000000000
    17455};
    17556
    … …  
    21798
    21899        for (unsigned int i = 0; i < OPERANDS; i++) {
    219                 double res = acos(arguments_acos[i]);
    220 
    221                 if (!cmp_double(res, results_acos[i])) {
    222                         TPRINTF("Double precision acos failed "
    223                             "(%lf != %lf, arg %u)\n", res, results_acos[i], i);
    224                         fail = true;
    225                 }
    226         }
    227 
    228         for (unsigned int i = 0; i < OPERANDS; i++) {
    229                 float res = acosf(arguments_acos[i]);
    230 
    231                 if (!cmp_float(res, results_acos[i])) {
    232                         TPRINTF("Single precision acos failed "
    233                             "(%f != %lf, arg %u)\n", res, results_acos[i], i);
    234                         fail = true;
    235                 }
    236         }
    237 
    238         for (unsigned int i = 0; i < OPERANDS; i++) {
    239                 double res = asin(arguments_asin[i]);
    240 
    241                 if (!cmp_double(res, results_asin[i])) {
    242                         TPRINTF("Double precision asin failed "
    243                             "(%lf != %lf, arg %u)\n", res, results_asin[i], i);
    244                         fail = true;
    245                 }
    246         }
    247 
    248         for (unsigned int i = 0; i < OPERANDS; i++) {
    249                 float res = asinf(arguments_asin[i]);
    250 
    251                 if (!cmp_float(res, results_asin[i])) {
    252                         TPRINTF("Single precision asin failed "
    253                             "(%f != %lf, arg %u)\n", res, results_asin[i], i);
    254                         fail = true;
    255                 }
    256         }
    257 
    258         for (unsigned int i = 0; i < OPERANDS; i++) {
    259                 double res = atan(arguments_atan[i]);
    260 
    261                 if (!cmp_double(res, results_atan[i])) {
    262                         TPRINTF("Double precision atan failed "
    263                             "(%.12lf != %.12lf, arg %u)\n", res, results_atan[i], i);
    264                         fail = true;
    265                 }
    266         }
    267 
    268         for (unsigned int i = 0; i < OPERANDS; i++) {
    269                 float res = atanf(arguments_atan[i]);
    270 
    271                 if (!cmp_float(res, results_atan[i])) {
    272                         TPRINTF("Single precision atan failed "
    273                             "(%f != %lf, arg %u)\n", res, results_atan[i], i);
    274                         fail = true;
    275                 }
    276         }
    277 
    278         for (unsigned int i = 0; i < OPERANDS; i++) {
    279                 double res = ceil(arguments[i]);
    280 
    281                 if (!cmp_double(res, results_ceil[i])) {
    282                         TPRINTF("Double precision ceil failed "
    283                             "(%lf != %lf, arg %u)\n", res, results_ceil[i], i);
    284                         fail = true;
    285                 }
    286         }
    287 
    288         for (unsigned int i = 0; i < OPERANDS; i++) {
    289                 float res = ceilf(arguments[i]);
    290 
    291                 if (!cmp_float(res, results_ceil[i])) {
    292                         TPRINTF("Single precision ceil failed "
    293                             "(%f != %lf, arg %u)\n", res, results_ceil[i], i);
    294                         fail = true;
    295                 }
    296         }
    297 
    298         for (unsigned int i = 0; i < OPERANDS; i++) {
    299100                double res = cos(arguments[i]);
    300101
    … …  
    317118
    318119        for (unsigned int i = 0; i < OPERANDS; i++) {
    319                 double res = cosh(arguments_exp[i]);
    320 
    321                 if (!cmp_double(res, results_cosh[i])) {
    322                         TPRINTF("Double precision cosh failed "
    323                             "(%lf != %lf, arg %u)\n", res, results_cosh[i], i);
    324                         fail = true;
    325                 }
    326         }
    327 
    328         for (unsigned int i = 0; i < OPERANDS; i++) {
    329                 float res = coshf(arguments_exp[i]);
    330 
    331                 if (!cmp_float(res, results_cosh[i])) {
    332                         TPRINTF("Single precision cosh failed "
    333                             "(%f != %lf, arg %u)\n", res, results_cosh[i], i);
    334                         fail = true;
    335                 }
    336         }
    337 
    338         for (unsigned int i = 0; i < OPERANDS; i++) {
    339                 double res = exp(arguments_exp[i]);
    340 
    341                 if (!cmp_double(res, results_exp[i])) {
    342                         TPRINTF("Double precision exp failed "
    343                             "(%lf != %lf, arg %u)\n", res, results_exp[i], i);
    344                         fail = true;
    345                 }
    346         }
    347 
    348         for (unsigned int i = 0; i < OPERANDS; i++) {
    349                 float res = expf(arguments_exp[i]);
    350 
    351                 if (!cmp_float(res, results_exp[i])) {
    352                         TPRINTF("Single precision exp failed "
    353                             "(%f != %lf, arg %u)\n", res, results_exp[i], i);
    354                         fail = true;
    355                 }
    356         }
    357 
    358         for (unsigned int i = 0; i < OPERANDS; i++) {
    359                 double res = fabs(arguments[i]);
    360 
    361                 if (!cmp_double(res, results_fabs[i])) {
    362                         TPRINTF("Double precision fabs failed "
    363                             "(%lf != %lf, arg %u)\n", res, results_fabs[i], i);
    364                         fail = true;
    365                 }
    366         }
    367 
    368         for (unsigned int i = 0; i < OPERANDS; i++) {
    369                 float res = fabsf(arguments[i]);
    370 
    371                 if (!cmp_float(res, results_fabs[i])) {
    372                         TPRINTF("Single precision fabs failed "
    373                             "(%f != %lf, arg %u)\n", res, results_fabs[i], i);
    374                         fail = true;
    375                 }
    376         }
    377 
    378         for (unsigned int i = 0; i < OPERANDS; i++) {
    379                 double res = floor(arguments[i]);
    380 
    381                 if (!cmp_double(res, results_floor[i])) {
    382                         TPRINTF("Double precision floor failed "
    383                             "(%lf != %lf, arg %u)\n", res, results_floor[i], i);
    384                         fail = true;
    385                 }
    386         }
    387 
    388         for (unsigned int i = 0; i < OPERANDS; i++) {
    389                 float res = floorf(arguments[i]);
    390 
    391                 if (!cmp_float(res, results_floor[i])) {
    392                         TPRINTF("Single precision floor failed "
    393                             "(%f != %lf, arg %u)\n", res, results_floor[i], i);
    394                         fail = true;
    395                 }
    396         }
    397 
    398         for (unsigned int i = 0; i < OPERANDS; i++) {
    399                 double res = log(arguments_log[i]);
    400 
    401                 if (!cmp_double(res, results_log[i])) {
    402                         TPRINTF("Double precision log failed "
    403                             "(%lf != %lf, arg %u)\n", res, results_log[i], i);
    404                         fail = true;
    405                 }
    406         }
    407 
    408         for (unsigned int i = 0; i < OPERANDS; i++) {
    409                 float res = logf(arguments_log[i]);
    410 
    411                 if (!cmp_float(res, results_log[i])) {
    412                         TPRINTF("Single precision log failed "
    413                             "(%f != %lf, arg %u)\n", res, results_log[i], i);
    414                         fail = true;
    415                 }
    416         }
    417 
    418         for (unsigned int i = 0; i < OPERANDS; i++) {
    419                 double res = log10(arguments_log[i]);
    420 
    421                 if (!cmp_double(res, results_log10[i])) {
    422                         TPRINTF("Double precision log10 failed "
    423                             "(%lf != %lf, arg %u)\n", res, results_log10[i], i);
    424                         fail = true;
    425                 }
    426         }
    427 
    428         for (unsigned int i = 0; i < OPERANDS; i++) {
    429                 float res = log10f(arguments_log[i]);
    430 
    431                 if (!cmp_float(res, results_log10[i])) {
    432                         TPRINTF("Single precision log10 failed "
    433                             "(%f != %lf, arg %u)\n", res, results_log10[i], i);
    434                         fail = true;
    435                 }
    436         }
    437 
    438         for (unsigned int i = 0; i < OPERANDS; i++) {
    439                 double res = log2(arguments_log[i]);
    440 
    441                 if (!cmp_double(res, results_log2[i])) {
    442                         TPRINTF("Double precision log2 failed "
    443                             "(%lf != %lf, arg %u)\n", res, results_log2[i], i);
    444                         fail = true;
    445                 }
    446         }
    447 
    448         for (unsigned int i = 0; i < OPERANDS; i++) {
    449                 float res = log2f(arguments_log[i]);
    450 
    451                 if (!cmp_float(res, results_log2[i])) {
    452                         TPRINTF("Single precision log2 failed "
    453                             "(%f != %lf, arg %u)\n", res, results_log2[i], i);
    454                         fail = true;
    455                 }
    456         }
    457 
    458         for (unsigned int i = 0; i < OPERANDS; i++) {
    459120                double res = sin(arguments[i]);
    460121
    … …  
    472133                        TPRINTF("Single precision sin failed "
    473134                            "(%f != %lf, arg %u)\n", res, results_sin[i], i);
    474                         fail = true;
    475                 }
    476         }
    477 
    478         for (unsigned int i = 0; i < OPERANDS; i++) {
    479                 double res = sinh(arguments_exp[i]);
    480 
    481                 if (!cmp_double(res, results_sinh[i])) {
    482                         TPRINTF("Double precision sinh failed "
    483                             "(%lf != %lf, arg %u)\n", res, results_sinh[i], i);
    484                         fail = true;
    485                 }
    486         }
    487 
    488         for (unsigned int i = 0; i < OPERANDS; i++) {
    489                 float res = sinhf(arguments_exp[i]);
    490 
    491                 if (!cmp_float(res, results_sinh[i])) {
    492                         TPRINTF("Single precision sinh failed "
    493                             "(%f != %lf, arg %u)\n", res, results_sinh[i], i);
    494                         fail = true;
    495                 }
    496         }
    497 
    498         for (unsigned int i = 0; i < OPERANDS; i++) {
    499                 double res = sqrt(arguments_sqrt[i]);
    500 
    501                 if (!cmp_double(res, results_sqrt[i])) {
    502                         TPRINTF("Double precision sqrt failed "
    503                             "(%lf != %lf, arg %u)\n", res, results_sqrt[i], i);
    504                         fail = true;
    505                 }
    506         }
    507 
    508         for (unsigned int i = 0; i < OPERANDS; i++) {
    509                 float res = sqrtf(arguments_sqrt[i]);
    510 
    511                 if (!cmp_float(res, results_sqrt[i])) {
    512                         TPRINTF("Single precision sqrt failed "
    513                             "(%f != %lf, arg %u)\n", res, results_sqrt[i], i);
    514                         fail = true;
    515                 }
    516         }
    517 
    518         for (unsigned int i = 0; i < OPERANDS; i++) {
    519                 double res = tan(arguments[i]);
    520 
    521                 if (!cmp_double(res, results_tan[i])) {
    522                         TPRINTF("Double precision tan failed "
    523                             "(%lf != %lf, arg %u)\n", res, results_tan[i], i);
    524                         fail = true;
    525                 }
    526         }
    527 
    528         for (unsigned int i = 0; i < OPERANDS; i++) {
    529                 float res = tanf(arguments[i]);
    530 
    531                 if (!cmp_float(res, results_tan[i])) {
    532                         TPRINTF("Single precision tan failed "
    533                             "(%f != %lf, arg %u)\n", res, results_tan[i], i);
    534                         fail = true;
    535                 }
    536         }
    537 
    538         for (unsigned int i = 0; i < OPERANDS; i++) {
    539                 double res = tanh(arguments_tanh[i]);
    540 
    541                 if (!cmp_double(res, results_tanh[i])) {
    542                         TPRINTF("Double precision tanh failed "
    543                             "(%lf != %lf, arg %u)\n", res, results_tanh[i], i);
    544                         fail = true;
    545                 }
    546         }
    547 
    548         for (unsigned int i = 0; i < OPERANDS; i++) {
    549                 float res = tanhf(arguments_tanh[i]);
    550 
    551                 if (!cmp_float(res, results_tanh[i])) {
    552                         TPRINTF("Single precision tanh failed "
    553                             "(%f != %lf, arg %u)\n", res, results_tanh[i], i);
    554135                        fail = true;
    555136                }

  • uspace/lib/c/include/fenv.h

    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  * @{
    31  */
    32 /** @file
    33  */
     29#ifndef _FENV_H
     30#define _FENV_H
    3431
    35 #ifndef LIBMATH_EXP_H_
    36 #define LIBMATH_EXP_H_
     32// TODO
    3733
    38 #include <mathtypes.h>
     34#define FE_TOWARDZERO  0
     35#define FE_TONEAREST   1
     36#define FE_UPWARD      2
     37#define FE_DOWNWARD    3
    3938
    40 extern float32_t float32_exp(float32_t);
    41 extern float64_t float64_exp(float64_t);
     39#define fegetround() FE_TONEAREST
    4240
    4341#endif
    4442
    45 /** @}
    46  */

  • uspace/lib/c/include/float.h

    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  * @{
    31  */
    32 /** @file
    33  */
     29#ifndef _FLOAT_H
     30#define _FLOAT_H
    3431
    35 #ifndef LIBMATH_ATAN_H_
    36 #define LIBMATH_ATAN_H_
     32// FIXME: <float.h> is freestanding. Just include the compiler-provided file.
    3733
    38 #include <mathtypes.h>
     34#define FLT_MIN __FLT_MIN__
     35#define FLT_DENORM_MIN __FLT_DENORM_MIN__
     36#define FLT_EPSILON __FLT_EPSILON__
    3937
    40 extern float32_t float32_atan(float32_t);
    41 extern float64_t float64_atan(float64_t);
     38#define FLT_MANT_DIG  __FLT_MANT_DIG__
     39#define DBL_MANT_DIG  __DBL_MANT_DIG__
     40#define FLT_MAX_EXP __FLT_MAX_EXP__
     41#define DBL_MAX_EXP __DBL_MAX_EXP__
    4242
    4343#endif
    4444
    45 /** @}
    46  */

  • 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  */

  • uspace/lib/posix/include/posix/math.h

    r7f7d642 r516e780  
    11/*
    2  * Copyright (c) 2011 Petr Koupy
     2 * Copyright (c) 2018 Jiri Svoboda
    33 * All rights reserved.
    44 *
    … …  
    3030 * @{
    3131 */
    32 /** @file Mathematical operations.
    33  *
    34  * The purpose of this file is only to provide prototypes of mathematical
    35  * functions defined by C standard and by POSIX.
    36  *
    37  * It is up to the application to correctly link with either libmath
    38  * (provided by HelenOS) or by some other math library (such as fdlibm).
    39  */
    4032
    4133#ifndef POSIX_MATH_H_
    4234#define POSIX_MATH_H_
    4335
    44 #ifdef __GNUC__
    45 #define HUGE_VAL (__builtin_huge_val())
     36/*
     37 * Just a pass-through to libc math.h
     38 */
     39#include "libc/math.h"
     40
    4641#endif
    47 
    48 extern double ldexp(double, int);
    49 extern double frexp(double, int *);
    50 
    51 extern double fabs(double);
    52 extern double floor(double);
    53 extern double ceil(double);
    54 extern double modf(double, double *);
    55 extern double fmod(double, double);
    56 extern double pow(double, double);
    57 extern double exp(double);
    58 extern double frexp(double, int *);
    59 extern double expm1(double);
    60 extern double sqrt(double);
    61 extern double log(double);
    62 extern double log10(double);
    63 extern double sin(double);
    64 extern double sinh(double);
    65 extern double asin(double);
    66 extern double asinh(double);
    67 extern double cos(double);
    68 extern double cosh(double);
    69 extern double acos(double);
    70 extern double acosh(double);
    71 extern double tan(double);
    72 extern double tanh(double);
    73 extern double atan(double);
    74 extern double atanh(double);
    75 extern double atan2(double, double);
    76 extern double copysign(double, double);
    77 
    78 #endif /* POSIX_MATH_H_ */
    7942
    8043/** @}

  • uspace/lib/softfloat/Makefile

    r7f7d642 r516e780  
    3030USPACE_PREFIX = ../..
    3131LIBRARY = libsoftfloat
    32 EXTRA_CFLAGS += $(LIBMATH_INCLUDES_FLAGS)
    3332
    3433SOURCES = \

  • uspace/lib/softfloat/comparison.c

    r7f7d642 r516e780  
    681681}
    682682
     683int __aeabi_fcmpun(float32_t a, float32_t b)
     684{
     685        float32_u ua;
     686        ua.val = a;
     687
     688        float32_u ub;
     689        ub.val = b;
     690
     691        // TODO: sigNaNs
     692        return is_float32_nan(ua.data) || is_float32_nan(ub.data);
     693}
     694
    683695#endif
    684696
    … …  
    920932}
    921933
     934int __aeabi_dcmpun(float64_t a, float64_t b)
     935{
     936        float64_u ua;
     937        ua.val = a;
     938
     939        float64_u ub;
     940        ub.val = b;
     941
     942        // TODO: sigNaNs
     943        return is_float64_nan(ua.data) || is_float64_nan(ub.data);
     944}
     945
    922946#endif
    923947

  • uspace/lib/softfloat/comparison.h

    r7f7d642 r516e780  
    9393extern int __aeabi_fcmple(float32_t, float32_t);
    9494extern int __aeabi_fcmpeq(float32_t, float32_t);
     95extern int __aeabi_fcmpun(float32_t, float32_t);
    9596#endif
    9697
    … …  
    109110extern int __aeabi_dcmpge(float64_t, float64_t);
    110111extern int __aeabi_dcmple(float64_t, float64_t);
     112extern int __aeabi_dcmpun(float64_t, float64_t);
    111113#endif
    112114
Note: See TracChangeset for help on using the changeset viewer.