Changeset 9adb61d in mainline for uspace/lib/math/generic/trig.c


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

File:
1 edited

Legend:

Unmodified
Added
Removed

  • 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
Note: See TracChangeset for help on using the changeset viewer.