Changeset 992ffa6 in mainline


Ignore:
Timestamp:
2015-09-04T06:40:20Z (9 years ago)
Author:
Jiri Svoboda <jiri@…>
Branches:
lfn, master, serial, ticket/834-toolchain-update, topic/msim-upgrade, topic/simplify-dev-export
Children:
01cdd5a
Parents:
bae1e1f
Message:

Add exp(f), log(f), pow(f). Improve precision of sin, cos.

Location:
uspace
Files:
6 added
9 edited

Legend:

Unmodified
Added
Removed

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

    rbae1e1f r992ffa6  
    11/*
    22 * Copyright (c) 2014 Martin Decky
     3 * Copyright (c) 2015 Jiri Svoboda
    34 * All rights reserved.
    45 *
    … …  
    2728 */
    2829
     30#include <stdbool.h>
    2931#include <stdio.h>
    3032#include <stdlib.h>
    … …  
    3234#include "../tester.h"
    3335
    34 #define OPERANDS   10
    35 #define PRECISION  100000000
     36#define OPERANDS         10
     37#define PRECISIONF    10000
     38#define PRECISION 100000000
    3639
    3740static double arguments[OPERANDS] = {
    3841        3.5, -2.1, 100.0, 50.0, -1024.0, 0.0, 768.3156, 1080.499999, -600.0, 1.0
    3942};
     43
     44static double arguments_exp[OPERANDS] = {
     45        3.5, -2.1, 50.0, 0.0, 1.0, 13.2, -1.1, -5.5, 0.1, -66.0
     46};
     47
     48static double arguments_log[OPERANDS] = {
     49        3.5, 100.0, 50.0, 768.3156, 1080.499999, 1.0, 66.0,
     50        2.718281828459045, 9.9, 0.001
     51};
     52
    4053
    4154static double results_ceil[OPERANDS] = {
    … …  
    6376};
    6477
     78static double results_log[OPERANDS] = {
     79        1.252762968495, 4.605170185988, 3.912023005428, 6.644200586236,
     80        6.985179175021, 0.000000000000, 4.189654742026, 1.000000000000,
     81        2.292534757141, -6.907755278982
     82};
     83
     84static double results_exp[OPERANDS] = {
     85        33.115451958692, 0.122456428253, 5184705528587072045056.0,
     86        1.000000000000, 2.718281828459, 540364.937246691552, 0.332871083698,
     87        0.004086771438, 1.105170918076, 0.000000000000
     88};
     89
     90static bool cmp_float(float a, float b)
     91{
     92        float r;
     93
     94        /* XXX Need fabsf() */
     95        if (b < 1.0 / PRECISIONF && b > -1.0 / PRECISIONF)
     96                r = a;
     97        else
     98                r = a / b - 1.0;
     99
     100        /* XXX Need fabsf() */
     101        if (r < 0.0)
     102                r = -r;
     103
     104        return r < 1.0 / PRECISIONF;
     105}
     106
     107static bool cmp_double(double a, double b)
     108{
     109        double r;
     110
     111        /* XXX Need fabs() */
     112        if (b < 1.0 / PRECISION && b > -1.0 / PRECISION)
     113                r = a;
     114        else
     115                r = a / b - 1.0;
     116
     117        /* XXX Need fabs() */
     118        if (r < 0.0)
     119                r = -r;
     120
     121        return r < 1.0 / PRECISION;
     122}
     123
    65124const char *test_float2(void)
    66125{
    67126        bool fail = false;
    68        
     127
    69128        for (unsigned int i = 0; i < OPERANDS; i++) {
    70129                double res = floor(arguments[i]);
    … …  
    73132               
    74133                if (res_int != corr_int) {
    75                         TPRINTF("Double floor failed (%" PRId64 " != %" PRId64
    76                             ", arg %u)\n", res_int, corr_int, i);
     134                        TPRINTF("Double precision floor failed (%" PRId64
     135                            " != %" PRId64 ", arg %u)\n", res_int, corr_int, i);
    77136                        fail = true;
    78137                }
    … …  
    85144               
    86145                if (res_int != corr_int) {
    87                         TPRINTF("Double ceil failed (%" PRId64 " != %" PRId64
    88                             ", arg %u)\n", res_int, corr_int, i);
     146                        TPRINTF("Double precision ceil failed (%" PRId64
     147                            " != %" PRId64 ", arg %u)\n", res_int, corr_int, i);
    89148                        fail = true;
    90149                }
    … …  
    97156               
    98157                if (res_int != corr_int) {
    99                         TPRINTF("Double truncation failed (%" PRId64 " != %" PRId64
    100                             ", arg %u)\n", res_int, corr_int, i);
     158                        TPRINTF("Double precisiontruncation failed (%" PRId64
     159                            " != %" PRId64 ", arg %u)\n", res_int, corr_int, i);
    101160                        fail = true;
    102161                }
    … …  
    109168               
    110169                if (res_int != corr_int) {
    111                         TPRINTF("Double sine failed (%" PRId64 " != %" PRId64
    112                             ", arg %u)\n", res_int, corr_int, i);
     170                        TPRINTF("Double precision sine failed (%" PRId64
     171                            " != %" PRId64 ", arg %u)\n", res_int, corr_int, i);
    113172                        fail = true;
    114173                }
    … …  
    121180               
    122181                if (res_int != corr_int) {
    123                         TPRINTF("Double cosine failed (%" PRId64 " != %" PRId64
    124                             ", arg %u)\n", res_int, corr_int, i);
     182                        TPRINTF("Double precision cosine failed (%" PRId64
     183                            " != %" PRId64 ", arg %u)\n", res_int, corr_int, i);
     184                        fail = true;
     185                }
     186        }
     187       
     188        for (unsigned int i = 0; i < OPERANDS; i++) {
     189                float res = logf(arguments_log[i]);
     190               
     191                if (!cmp_float(res, results_log[i])) {
     192                        TPRINTF("Single precision logarithm failed "
     193                            "(%lf != %lf, arg %u)\n", res, results_log[i], i);
     194                        fail = true;
     195                }
     196        }
     197       
     198        for (unsigned int i = 0; i < OPERANDS; i++) {
     199                double res = log(arguments_log[i]);
     200               
     201                if (!cmp_double(res, results_log[i])) {
     202                        TPRINTF("Double precision logarithm failed "
     203                            "(%lf != %lf, arg %u)\n", res, results_log[i], i);
     204                        fail = true;
     205                }
     206        }
     207       
     208        for (unsigned int i = 0; i < OPERANDS; i++) {
     209                float res = exp(arguments_exp[i]);
     210               
     211                if (!cmp_float(res, results_exp[i])) {
     212                        TPRINTF("Single precision exponential failed "
     213                            "(%lf != %lf, arg %u)\n", res, results_exp[i], i);
     214                        fail = true;
     215                }
     216        }
     217       
     218        for (unsigned int i = 0; i < OPERANDS; i++) {
     219                double res = exp(arguments_exp[i]);
     220               
     221                if (!cmp_double(res, results_exp[i])) {
     222                        TPRINTF("Double precision exponential failed "
     223                            "(%lf != %lf, arg %u)\n", res, results_exp[i], i);
    125224                        fail = true;
    126225                }

  • uspace/lib/math/Makefile

    rbae1e1f r992ffa6  
    4242GENERIC_SOURCES = \
    4343        generic/ceil.c \
     44        generic/exp.c \
    4445        generic/floor.c \
     46        generic/log.c \
    4547        generic/trig.c \
     48        generic/pow.c \
    4649        generic/mod.c \
    4750        generic/trunc.c

  • uspace/lib/math/arch/amd64/include/libarch/math.h

    rbae1e1f r992ffa6  
    3737
    3838#include <ceil.h>
     39#include <exp.h>
    3940#include <floor.h>
     41#include <log.h>
    4042#include <mathtypes.h>
    4143#include <mod.h>
     44#include <pow.h>
    4245
    4346static inline float64_t fmod(float64_t dividend, float64_t divisor)
    … …  
    6063}
    6164
     65static inline float32_t expf(float32_t val)
     66{
     67        return float32_exp(val);
     68}
     69
     70static inline float64_t exp(float64_t val)
     71{
     72        return float64_exp(val);
     73}
     74
    6275static inline float64_t floor(float64_t val)
    6376{
    … …  
    7184}
    7285
     86static inline float32_t logf(float32_t val)
     87{
     88        return float32_log(val);
     89}
     90
     91static inline float64_t log(float64_t val)
     92{
     93        return float64_log(val);
     94}
     95
     96static inline float32_t powf(float32_t x, float32_t y)
     97{
     98        return float32_pow(x, y);
     99}
     100
     101static inline float64_t pow(float64_t x, float64_t y)
     102{
     103        return float64_pow(x, y);
     104}
     105
    73106extern float64_t trunc(float64_t);
    74107

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

    rbae1e1f r992ffa6  
    4848        float64_u v;
    4949        float64_u r;
    50 
     50       
    5151        v.data = val;
    5252        t.data = trunc_float64(val);
    53 
     53       
    5454        if (val.parts.sign == 1 || v.val == t.val) {
    5555                r = t;
    … …  
    5757                r.val = t.val + 1.0;
    5858        }
    59 
     59       
    6060        return r.data;
    6161}

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

    rbae1e1f r992ffa6  
    4949        float64_u v;
    5050        float64_u r;
    51 
     51       
    5252        v.data = val;
    5353        t.data = trunc_float64(val);
    54 
     54       
    5555        if (val.parts.sign == 0 || v.val == t.val) {
    5656                r = t;
    … …  
    5858                r.val = t.val - 1.0;
    5959        }
    60 
     60       
    6161        return r.data;
    6262}

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

    rbae1e1f r992ffa6  
    3636#include <trig.h>
    3737
    38 #define TAYLOR_DEGREE  13
     38#define TAYLOR_DEGREE_32 13
     39#define TAYLOR_DEGREE_64 21
    3940
    4041/** Precomputed values for factorial (starting from 1!) */
    41 static float64_t factorials[TAYLOR_DEGREE] = {
     42static float64_t factorials[TAYLOR_DEGREE_64] = {
    4243        1, 2, 6, 24, 120, 720, 5040, 40320, 362880, 3628800, 39916800,
    43         479001600, 6227020800
     44        479001600, 6227020800.0L, 87178291200.0L, 1307674368000.0L,
     45        20922789888000.0L, 355687428096000.0L, 6402373705728000.0L,
     46        121645100408832000.0L, 2432902008176640000.0L, 51090942171709440000.0L
    4447};
    4548
    … …  
    6164        float64_t nom = 1;
    6265       
    63         for (unsigned int i = 0; i < TAYLOR_DEGREE; i++) {
     66        for (unsigned int i = 0; i < TAYLOR_DEGREE_64; i++) {
    6467                nom *= arg;
    6568               
    … …  
    9093        float64_t nom = 1;
    9194       
    92         for (unsigned int i = 0; i < TAYLOR_DEGREE; i++) {
     95        for (unsigned int i = 0; i < TAYLOR_DEGREE_64; i++) {
    9396                nom *= arg;
    9497               

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

    rbae1e1f r992ffa6  
    3535#include <mathtypes.h>
    3636#include <trunc.h>
     37
     38/** Truncate fractional part (round towards zero)
     39 *
     40 * Truncate the fractional part of IEEE 754 single
     41 * precision floating point number by zeroing fraction
     42 * bits, effectively rounding the number towards zero
     43 * to the nearest whole number.
     44 *
     45 * If the argument is infinity or NaN, an exception
     46 * should be indicated. This is not implemented yet.
     47 *
     48 * @param val Floating point number.
     49 *
     50 * @return Number rounded towards zero.
     51 *
     52 */
     53float32 trunc_float32(float32 val)
     54{
     55        int32_t exp = val.parts.exp - FLOAT32_BIAS;
     56       
     57        if (exp < 0) {
     58                /* -1 < val < 1 => result is +0 or -0 */
     59                val.parts.exp = 0;
     60                val.parts.fraction = 0;
     61        } else if (exp >= FLOAT32_FRACTION_SIZE) {
     62                if (exp == 1024) {
     63                        /* val is +inf, -inf or NaN => trigger an exception */
     64                        // FIXME TODO
     65                }
     66               
     67                /* All bits in val are relevant for the result */
     68        } else {
     69                /* Truncate irrelevant fraction bits */
     70                val.parts.fraction &= ~(UINT32_C(0x007fffff) >> exp);
     71        }
     72       
     73        return val;
     74}
    3775
    3876/** Truncate fractional part (round towards zero)

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

    rbae1e1f r992ffa6  
    3838#include <libarch/math.h>
    3939
    40 #define M_PI  3.14159265358979323846
     40#define M_LN2 0.69314718055994530942
     41#define M_LOG2E 1.4426950408889634074
     42#define M_PI 3.14159265358979323846
    4143
    4244#endif

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

    rbae1e1f r992ffa6  
    3838#include <mathtypes.h>
    3939
     40extern float32 trunc_float32(float32);
    4041extern float64 trunc_float64(float64);
    4142
Note: See TracChangeset for help on using the changeset viewer.