Ignore:
File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r9adb61d rc0c38c7c  
    11/*
    2  * Copyright (c) 2015 Jiri Svoboda
    32 * Copyright (c) 2014 Martin Decky
    43 * All rights reserved.
    … …  
    3736#include <trig.h>
    3837
    39 #define TAYLOR_DEGREE_32 13
    40 #define TAYLOR_DEGREE_64 21
     38#define TAYLOR_DEGREE  13
    4139
    4240/** Precomputed values for factorial (starting from 1!) */
    43 static float64_t factorials[TAYLOR_DEGREE_64] = {
     41static float64_t factorials[TAYLOR_DEGREE] = {
    4442        1, 2, 6, 24, 120, 720, 5040, 40320, 362880, 3628800, 39916800,
    45         479001600, 6227020800.0L, 87178291200.0L, 1307674368000.0L,
    46         20922789888000.0L, 355687428096000.0L, 6402373705728000.0L,
    47         121645100408832000.0L, 2432902008176640000.0L, 51090942171709440000.0L
     43        479001600, 6227020800
    4844};
    4945
    50 /** Sine approximation by Taylor series (32-bit floating point)
     46/** Sine approximation by Taylor series
    5147 *
    5248 * Compute the approximation of sine by a Taylor
    … …  
    6056 *
    6157 */
    62 static 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++) {
     58static float64_t taylor_sin(float64_t arg)
     59{
     60        float64_t ret = 0;
     61        float64_t nom = 1;
     62       
     63        for (unsigned int i = 0; i < TAYLOR_DEGREE; i++) {
    6864                nom *= arg;
    6965               
    … …  
    7773}
    7874
    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  */
    91 static float64_t taylor_sin_64(float64_t arg)
    92 {
    93         float64_t ret = 0;
    94         float64_t nom = 1;
    95        
    96         for (unsigned int i = 0; i < TAYLOR_DEGREE_64; i++) {
    97                 nom *= arg;
    98                
    99                 if ((i % 4) == 0)
    100                         ret += nom / factorials[i];
    101                 else if ((i % 4) == 2)
    102                         ret -= nom / factorials[i];
    103         }
    104        
    105         return ret;
    106 }
    107 
    108 /** Cosine approximation by Taylor series (32-bit floating point)
     75/** Cosine approximation by Taylor series
    10976 *
    11077 * Compute the approximation of cosine by a Taylor
    … …  
    11885 *
    11986 */
    120 static 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++) {
     87static float64_t taylor_cos(float64_t arg)
     88{
     89        float64_t ret = 1;
     90        float64_t nom = 1;
     91       
     92        for (unsigned int i = 0; i < TAYLOR_DEGREE; i++) {
    12693                nom *= arg;
    12794               
    … …  
    135102}
    136103
    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  */
    149 static float64_t taylor_cos_64(float64_t arg)
    150 {
    151         float64_t ret = 1;
    152         float64_t nom = 1;
    153        
    154         for (unsigned int i = 0; i < TAYLOR_DEGREE_64; i++) {
    155                 nom *= arg;
    156                
    157                 if ((i % 4) == 1)
    158                         ret -= nom / factorials[i];
    159                 else if ((i % 4) == 3)
    160                         ret += nom / factorials[i];
    161         }
    162        
    163         return ret;
    164 }
    165 
    166 /** Sine value for values within base period (32-bit floating point)
     104/** Sine value for values within base period
    167105 *
    168106 * Compute the value of sine for arguments within
    … …  
    176114 *
    177115 */
    178 static float32_t base_sin_32(float32_t arg)
     116static float64_t base_sin(float64_t arg)
    179117{
    180118        unsigned int period = arg / (M_PI / 4);
    … …  
    182120        switch (period) {
    183121        case 0:
    184                 return taylor_sin_32(arg);
     122                return taylor_sin(arg);
    185123        case 1:
    186124        case 2:
    187                 return taylor_cos_32(arg - M_PI / 2);
     125                return taylor_cos(arg - M_PI / 2);
    188126        case 3:
    189127        case 4:
    190                 return -taylor_sin_32(arg - M_PI);
     128                return -taylor_sin(arg - M_PI);
    191129        case 5:
    192130        case 6:
    193                 return -taylor_cos_32(arg - 3 * M_PI / 2);
     131                return -taylor_cos(arg - 3 * M_PI / 2);
    194132        default:
    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  */
    211 static 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)
     133                return taylor_sin(arg - 2 * M_PI);
     134        }
     135}
     136
     137/** Cosine value for values within base period
    233138 *
    234139 * Compute the value of cosine for arguments within
    … …  
    242147 *
    243148 */
    244 static float32_t base_cos_32(float32_t arg)
     149static float64_t base_cos(float64_t arg)
    245150{
    246151        unsigned int period = arg / (M_PI / 4);
    … …  
    248153        switch (period) {
    249154        case 0:
    250                 return taylor_cos_32(arg);
     155                return taylor_cos(arg);
    251156        case 1:
    252157        case 2:
    253                 return -taylor_sin_32(arg - M_PI / 2);
     158                return -taylor_sin(arg - M_PI / 2);
    254159        case 3:
    255160        case 4:
    256                 return -taylor_cos_32(arg - M_PI);
     161                return -taylor_cos(arg - M_PI);
    257162        case 5:
    258163        case 6:
    259                 return taylor_sin_32(arg - 3 * M_PI / 2);
     164                return taylor_sin(arg - 3 * M_PI / 2);
    260165        default:
    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.
     166                return taylor_cos(arg - 2 * M_PI);
     167        }
     168}
     169
     170/** Double precision sine
     171 *
     172 * Compute sine value.
     173 *
     174 * @param arg Sine argument.
     175 *
     176 * @return Sine value.
     177 *
     178 */
     179float64_t float64_sin(float64_t arg)
     180{
     181        float64_t base_arg = fmod(arg, 2 * M_PI);
     182       
     183        if (base_arg < 0)
     184                return -base_sin(-base_arg);
     185       
     186        return base_sin(base_arg);
     187}
     188
     189/** Double precision cosine
     190 *
     191 * Compute cosine value.
    271192 *
    272193 * @param arg Cosine argument.
    … …  
    275196 *
    276197 */
    277 static 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)
    299  *
    300  * Compute sine value.
    301  *
    302  * @param arg Sine argument.
    303  *
    304  * @return Sine value.
    305  *
    306  */
    307 float32_t float32_sin(float32_t arg)
    308 {
    309         float32_t base_arg = fmod_f32(arg, 2 * M_PI);
     198float64_t float64_cos(float64_t arg)
     199{
     200        float64_t base_arg = fmod(arg, 2 * M_PI);
    310201       
    311202        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  */
    326 float64_t float64_sin(float64_t arg)
    327 {
    328         float64_t base_arg = fmod_f64(arg, 2 * M_PI);
    329        
    330         if (base_arg < 0)
    331                 return -base_sin_64(-base_arg);
    332        
    333         return base_sin_64(base_arg);
    334 }
    335 
    336 /** Cosine (32-bit floating point)
    337  *
    338  * Compute cosine value.
    339  *
    340  * @param arg Cosine argument.
    341  *
    342  * @return Cosine value.
    343  *
    344  */
    345 float32_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  */
    364 float64_t float64_cos(float64_t arg)
    365 {
    366         float64_t base_arg = fmod_f64(arg, 2 * M_PI);
    367        
    368         if (base_arg < 0)
    369                 return base_cos_64(-base_arg);
    370        
    371         return base_cos_64(base_arg);
     203                return base_cos(-base_arg);
     204       
     205        return base_cos(base_arg);
    372206}
    373207
Note: See TracChangeset for help on using the changeset viewer.