Changeset 9adb61d in mainline for uspace/lib/math/generic


Ignore:
Timestamp:
2015-09-05T11:50:00Z (10 years ago)
Author:
Jiri Svoboda <jiri@…>
Branches:
lfn, master, serial, ticket/834-toolchain-update, topic/msim-upgrade, topic/simplify-dev-export
Children:
996dc042, ba8eecf
Parents:
e6f5766
Message:

Add single-precision variant for all functions. Allow generic implementations to call other functions while selecting the number of bits of precision, but not the implementation (generic or arch-specific).

Location:
uspace/lib/math/generic
Files:
8 edited

Legend:

Unmodified
Added
Removed

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

    re6f5766 r9adb61d  
    3434
    3535#include <ceil.h>
     36#include <math.h>
    3637#include <mathtypes.h>
    37 #include <trunc.h>
    3838
    39 /** Ceiling (round towards positive infinity)
     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 */
     45float32_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;
     58        }
     59       
     60        return r.val;
     61}
     62
     63/** Ceiling (round towards positive infinity, 64-bit floating point)
    4064 *
    4165 * @param val Floating point number.
    … …  
    5074       
    5175        v.val = val;
    52         t.val = float64_trunc(val);
     76        t.val = trunc_f64(val);
    5377       
    5478        if (v.data.parts.sign == 1 || val == t.val) {

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

    re6f5766 r9adb61d  
    4949};
    5050
    51 /** Exponential approximation by Taylor series (32-))
     51/** Exponential approximation by Taylor series (32-bit floating point)
    5252 *
    5353 * Compute the approximation of exponential by a Taylor
    … …  
    6161 *
    6262 */
    63 static float64_t taylor_exp_32(float64_t arg)
     63static float32_t taylor_exp_32(float32_t arg)
    6464{
    65         float64_t ret = 1;
    66         float64_t nom = 1;
     65        float32_t ret = 1;
     66        float32_t nom = 1;
    6767       
    6868        for (unsigned int i = 0; i < TAYLOR_DEGREE_32; i++) {
    … …  
    7474}
    7575
    76 /** Exponential approximation by Taylor series
     76/** Exponential approximation by Taylor series (64-bit floating point)
    7777 *
    7878 * Compute the approximation of exponential by a Taylor
    … …  
    9999}
    100100
    101 /** Single precision exponential
     101/** Exponential (32-bit floating point)
    102102 *
    103103 * Compute exponential value.
    … …  
    121121         */
    122122
    123         i = float32_trunc(arg * M_LOG2E);
     123        i = trunc_f32(arg * M_LOG2E);
    124124        f = arg * M_LOG2E - i;
    125125
    … …  
    129129}
    130130
    131 /** Double precision exponential
     131/** Exponential (64-bit floating point)
    132132 *
    133133 * Compute exponential value.
    … …  
    151151         */
    152152
    153         i = float64_trunc(arg * M_LOG2E);
     153        i = trunc_f64(arg * M_LOG2E);
    154154        f = arg * M_LOG2E - i;
    155155

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

    re6f5766 r9adb61d  
    3434
    3535#include <floor.h>
     36#include <math.h>
    3637#include <mathtypes.h>
    37 #include <trunc.h>
    3838
    39 /** Ceiling (round towards negative infinity)
     39/** Ceiling (round towards negative infinity, 32-bit floating point)
     40 *
     41 * @param val Floating point number.
     42 *
     43 * @return Number rounded towards negative infinity.
     44 */
     45
     46float32_t float32_floor(float32_t val)
     47{
     48        float32_t t;
     49        float32_u v;
     50       
     51        v.val = val;
     52        t = trunc_f32(val);
     53       
     54        if (v.data.parts.sign == 0 || val == t)
     55                return t;
     56        else
     57                return t - 1.0;
     58}
     59
     60/** Ceiling (round towards negative infinity, 64-bit floating point)
    4061 *
    4162 * @param val Floating point number.
    … …  
    5071       
    5172        v.val = val;
    52         t = float64_trunc(val);
     73        t = trunc_f64(val);
    5374       
    5475        if (v.data.parts.sign == 0 || val == t)

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

    re6f5766 r9adb61d  
    4040#define TAYLOR_DEGREE_64  63
    4141
    42 /** Single precision log(1 - arg) approximation by Taylor series
     42/** log(1 - arg) approximation by Taylor series (32-bit floating point)
    4343 *
    4444 * Compute the approximation of log(1 - arg) by a Taylor
    … …  
    6868}
    6969
    70 /** Double precision log(1 - arg) approximation by Taylor series
     70/** log(1 - arg) approximation by Taylor series (64-bit floating point)
    7171 *
    7272 * Compute the approximation of log(1 - arg) by a Taylor
    … …  
    9696}
    9797
    98 /** Single precision logarithm
     98/** Natural logarithm (32-bit floating point)
    9999 *
    100100 * @param arg Argument.
    … …  
    126126}
    127127
    128 /** Double precision logarithm
     128/** Natural logarithm (64-bit floating point)
    129129 *
    130130 * @param arg Argument.

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

    re6f5766 r9adb61d  
    3636#include <mod.h>
    3737
    38 /** Double precision modulo
     38/** Remainder function (32-bit floating point)
     39 *
     40 * Calculate the modulo of dividend by divisor.
     41 *
     42 * This is a very basic implementation that uses
     43 * division and multiplication (instead of exact
     44 * arithmetics). Thus the result might be very
     45 * imprecise (depending on the magnitude of the
     46 * arguments).
     47 *
     48 * @param dividend Dividend.
     49 * @param divisor  Divisor.
     50 *
     51 * @return Modulo.
     52 *
     53 */
     54float32_t float32_mod(float32_t dividend, float32_t divisor)
     55{
     56        // FIXME: replace with exact arithmetics
     57       
     58        float32_t quotient = trunc_f32(dividend / divisor);
     59       
     60        return (dividend - quotient * divisor);
     61}
     62
     63/** Remainder function (64-bit floating point)
    3964 *
    4065 * Calculate the modulo of dividend by divisor.
    … …  
    5681        // FIXME: replace with exact arithmetics
    5782       
    58         float64_t quotient = trunc(dividend / divisor);
     83        float64_t quotient = trunc_f64(dividend / divisor);
    5984       
    6085        return (dividend - quotient * divisor);

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

    re6f5766 r9adb61d  
    5252{
    5353        /* x^y = (e ^ log(x))^y = e ^ (log(x) * y) */
    54         return float32_exp(float32_log(x) * y);
     54        return exp_f32(log_f32(x) * y);
    5555}
    5656
    … …  
    6868{
    6969        /* x^y = (e ^ log(x))^y = e ^ (log(x) * y) */
    70         return float64_exp(float64_log(x) * y);
     70        return exp_f64(log_f64(x) * y);
    7171}
    7272

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

    re6f5766 r9adb61d  
    11/*
     2 * Copyright (c) 2015 Jiri Svoboda
    23 * Copyright (c) 2014 Martin Decky
    34 * All rights reserved.
    … …  
    4748};
    4849
    49 /** Sine approximation by Taylor series
     50/** Sine approximation by Taylor series (32-bit floating point)
    5051 *
    5152 * Compute the approximation of sine by a Taylor
    … …  
    5960 *
    6061 */
    61 static float64_t taylor_sin(float64_t arg)
     62static float32_t taylor_sin_32(float32_t arg)
     63{
     64        float32_t ret = 0;
     65        float32_t nom = 1;
     66       
     67        for (unsigned int i = 0; i < TAYLOR_DEGREE_32; i++) {
     68                nom *= arg;
     69               
     70                if ((i % 4) == 0)
     71                        ret += nom / factorials[i];
     72                else if ((i % 4) == 2)
     73                        ret -= nom / factorials[i];
     74        }
     75       
     76        return ret;
     77}
     78
     79/** Sine approximation by Taylor series (64-bit floating point)
     80 *
     81 * Compute the approximation of sine by a Taylor
     82 * series (using the first TAYLOR_DEGREE terms).
     83 * The approximation is reasonably accurate for
     84 * arguments within the interval [-pi/4, pi/4].
     85 *
     86 * @param arg Sine argument.
     87 *
     88 * @return Sine value approximation.
     89 *
     90 */
     91static float64_t taylor_sin_64(float64_t arg)
    6292{
    6393        float64_t ret = 0;
    … …  
    76106}
    77107
    78 /** Cosine approximation by Taylor series
     108/** Cosine approximation by Taylor series (32-bit floating point)
    79109 *
    80110 * Compute the approximation of cosine by a Taylor
    … …  
    88118 *
    89119 */
    90 static float64_t taylor_cos(float64_t arg)
     120static float32_t taylor_cos_32(float32_t arg)
     121{
     122        float32_t ret = 1;
     123        float32_t nom = 1;
     124       
     125        for (unsigned int i = 0; i < TAYLOR_DEGREE_32; i++) {
     126                nom *= arg;
     127               
     128                if ((i % 4) == 1)
     129                        ret -= nom / factorials[i];
     130                else if ((i % 4) == 3)
     131                        ret += nom / factorials[i];
     132        }
     133       
     134        return ret;
     135}
     136
     137/** Cosine approximation by Taylor series (64-bit floating point)
     138 *
     139 * Compute the approximation of cosine by a Taylor
     140 * series (using the first TAYLOR_DEGREE terms).
     141 * The approximation is reasonably accurate for
     142 * arguments within the interval [-pi/4, pi/4].
     143 *
     144 * @param arg Cosine argument.
     145 *
     146 * @return Cosine value approximation.
     147 *
     148 */
     149static float64_t taylor_cos_64(float64_t arg)
    91150{
    92151        float64_t ret = 1;
    … …  
    105164}
    106165
    107 /** Sine value for values within base period
     166/** Sine value for values within base period (32-bit floating point)
    108167 *
    109168 * Compute the value of sine for arguments within
    … …  
    117176 *
    118177 */
    119 static float64_t base_sin(float64_t arg)
     178static float32_t base_sin_32(float32_t arg)
    120179{
    121180        unsigned int period = arg / (M_PI / 4);
    … …  
    123182        switch (period) {
    124183        case 0:
    125                 return taylor_sin(arg);
     184                return taylor_sin_32(arg);
    126185        case 1:
    127186        case 2:
    128                 return taylor_cos(arg - M_PI / 2);
     187                return taylor_cos_32(arg - M_PI / 2);
    129188        case 3:
    130189        case 4:
    131                 return -taylor_sin(arg - M_PI);
     190                return -taylor_sin_32(arg - M_PI);
    132191        case 5:
    133192        case 6:
    134                 return -taylor_cos(arg - 3 * M_PI / 2);
     193                return -taylor_cos_32(arg - 3 * M_PI / 2);
    135194        default:
    136                 return taylor_sin(arg - 2 * M_PI);
    137         }
    138 }
    139 
    140 /** Cosine value for values within base period
     195                return taylor_sin_32(arg - 2 * M_PI);
     196        }
     197}
     198
     199/** Sine value for values within base period (64-bit floating point)
     200 *
     201 * Compute the value of sine for arguments within
     202 * the base period [0, 2pi]. For arguments outside
     203 * the base period the returned values can be
     204 * very inaccurate or even completely wrong.
     205 *
     206 * @param arg Sine argument.
     207 *
     208 * @return Sine value.
     209 *
     210 */
     211static float64_t base_sin_64(float64_t arg)
     212{
     213        unsigned int period = arg / (M_PI / 4);
     214       
     215        switch (period) {
     216        case 0:
     217                return taylor_sin_64(arg);
     218        case 1:
     219        case 2:
     220                return taylor_cos_64(arg - M_PI / 2);
     221        case 3:
     222        case 4:
     223                return -taylor_sin_64(arg - M_PI);
     224        case 5:
     225        case 6:
     226                return -taylor_cos_64(arg - 3 * M_PI / 2);
     227        default:
     228                return taylor_sin_64(arg - 2 * M_PI);
     229        }
     230}
     231
     232/** Cosine value for values within base period (32-bit floating point)
    141233 *
    142234 * Compute the value of cosine for arguments within
    … …  
    150242 *
    151243 */
    152 static float64_t base_cos(float64_t arg)
     244static float32_t base_cos_32(float32_t arg)
    153245{
    154246        unsigned int period = arg / (M_PI / 4);
    … …  
    156248        switch (period) {
    157249        case 0:
    158                 return taylor_cos(arg);
     250                return taylor_cos_32(arg);
    159251        case 1:
    160252        case 2:
    161                 return -taylor_sin(arg - M_PI / 2);
     253                return -taylor_sin_32(arg - M_PI / 2);
    162254        case 3:
    163255        case 4:
    164                 return -taylor_cos(arg - M_PI);
     256                return -taylor_cos_32(arg - M_PI);
    165257        case 5:
    166258        case 6:
    167                 return taylor_sin(arg - 3 * M_PI / 2);
     259                return taylor_sin_32(arg - 3 * M_PI / 2);
    168260        default:
    169                 return taylor_cos(arg - 2 * M_PI);
    170         }
    171 }
    172 
    173 /** Double precision sine
     261                return taylor_cos_32(arg - 2 * M_PI);
     262        }
     263}
     264
     265/** Cosine value for values within base period (64-bit floating point)
     266 *
     267 * Compute the value of cosine for arguments within
     268 * the base period [0, 2pi]. For arguments outside
     269 * the base period the returned values can be
     270 * very inaccurate or even completely wrong.
     271 *
     272 * @param arg Cosine argument.
     273 *
     274 * @return Cosine value.
     275 *
     276 */
     277static float64_t base_cos_64(float64_t arg)
     278{
     279        unsigned int period = arg / (M_PI / 4);
     280       
     281        switch (period) {
     282        case 0:
     283                return taylor_cos_64(arg);
     284        case 1:
     285        case 2:
     286                return -taylor_sin_64(arg - M_PI / 2);
     287        case 3:
     288        case 4:
     289                return -taylor_cos_64(arg - M_PI);
     290        case 5:
     291        case 6:
     292                return taylor_sin_64(arg - 3 * M_PI / 2);
     293        default:
     294                return taylor_cos_64(arg - 2 * M_PI);
     295        }
     296}
     297
     298/** Sine (32-bit floating point)
    174299 *
    175300 * Compute sine value.
    … …  
    180305 *
    181306 */
     307float32_t float32_sin(float32_t arg)
     308{
     309        float32_t base_arg = fmod_f32(arg, 2 * M_PI);
     310       
     311        if (base_arg < 0)
     312                return -base_sin_32(-base_arg);
     313       
     314        return base_sin_32(base_arg);
     315}
     316
     317/** Sine (64-bit floating point)
     318 *
     319 * Compute sine value.
     320 *
     321 * @param arg Sine argument.
     322 *
     323 * @return Sine value.
     324 *
     325 */
    182326float64_t float64_sin(float64_t arg)
    183327{
    184         float64_t base_arg = fmod(arg, 2 * M_PI);
     328        float64_t base_arg = fmod_f64(arg, 2 * M_PI);
    185329       
    186330        if (base_arg < 0)
    187                 return -base_sin(-base_arg);
    188        
    189         return base_sin(base_arg);
    190 }
    191 
    192 /** Double precision cosine
     331                return -base_sin_64(-base_arg);
     332       
     333        return base_sin_64(base_arg);
     334}
     335
     336/** Cosine (32-bit floating point)
    193337 *
    194338 * Compute cosine value.
    … …  
    199343 *
    200344 */
     345float32_t float32_cos(float32_t arg)
     346{
     347        float32_t base_arg = fmod_f32(arg, 2 * M_PI);
     348       
     349        if (base_arg < 0)
     350                return base_cos_32(-base_arg);
     351       
     352        return base_cos_32(base_arg);
     353}
     354
     355/** Cosine (64-bit floating point)
     356 *
     357 * Compute cosine value.
     358 *
     359 * @param arg Cosine argument.
     360 *
     361 * @return Cosine value.
     362 *
     363 */
    201364float64_t float64_cos(float64_t arg)
    202365{
    203         float64_t base_arg = fmod(arg, 2 * M_PI);
     366        float64_t base_arg = fmod_f64(arg, 2 * M_PI);
    204367       
    205368        if (base_arg < 0)
    206                 return base_cos(-base_arg);
    207        
    208         return base_cos(base_arg);
     369                return base_cos_64(-base_arg);
     370       
     371        return base_cos_64(base_arg);
    209372}
    210373

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

    re6f5766 r9adb61d  
    11/*
     2 * Copyright (c) 2015 Jiri Svoboda
    23 * Copyright (c) 2014 Martin Decky
    34 * All rights reserved.
Note: See TracChangeset for help on using the changeset viewer.