Changeset 516e780 in mainline for uspace/lib/math
- Timestamp:
- 2018-08-31T11:55:41Z (7 years ago)
- 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)
- Location:
- uspace/lib/math
- Files:
-
- 2 added
- 64 deleted
- 5 edited
- 5 moved
Legend:
- Unmodified
- Added
- Removed
-
uspace/lib/math/Makefile
r7f7d642 r516e780 30 30 ROOT_PATH = $(USPACE_PREFIX)/.. 31 31 32 CONFIG_MAKEFILE = $(ROOT_PATH)/Makefile.config33 34 32 LIBRARY = libmath 35 33 SOVERSION = 0.0 36 34 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 \ 35 SOURCES = \ 36 generic/__fcompare.c \ 37 generic/__fpclassify.c \ 38 generic/__signbit.c \ 50 39 generic/fabs.c \ 51 generic/floor.c \52 40 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 \ 64 43 generic/trig.c \ 65 44 generic/trunc.c 66 45 67 SOURCES = \68 $(GENERIC_SOURCES)\69 $(ARCH_SOURCES)46 TEST_SOURCES = \ 47 test/rounding.c \ 48 test/main.c 70 49 71 50 include $(USPACE_PREFIX)/Makefile.common -
uspace/lib/math/generic/__fcompare.c
r7f7d642 r516e780 1 1 /* 2 * Copyright (c) 201 5 Jiri Svoboda2 * Copyright (c) 2018 CZ.NIC, z.s.p.o. 3 3 * All rights reserved. 4 4 * … … 27 27 */ 28 28 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. 31 35 */ 32 /** @file 33 */ 36 int __fcompare(size_t sz1, size_t sz2, ...) 37 { 38 va_list ap; 39 va_start(ap, sz2); 34 40 35 #include <ceil.h> 36 #include <math.h> 37 #include <mathtypes.h> 41 long double val1; 42 long double val2; 38 43 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; 58 54 } 59 55 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 } 61 81 } 62 82 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 1 1 /* 2 * Copyright (c) 201 5 Jiri Svoboda2 * Copyright (c) 2018 CZ.NIC, z.s.p.o. 3 3 * All rights reserved. 4 4 * … … 27 27 */ 28 28 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. 31 35 */ 32 /** @file 33 */ 36 int __fpclassify(size_t sz, ...) 37 { 38 va_list ap; 39 va_start(ap, sz); 34 40 35 #include <math.h> 36 #include <sqrt.h> 41 int result; 37 42 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; 50 57 } 51 58 52 /** Double precision square root53 *54 * Compute squre root.55 *56 * @param val Value57 *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 1 1 /* 2 * Copyright (c) 201 5 Jiri Svoboda2 * Copyright (c) 2018 CZ.NIC, z.s.p.o. 3 3 * All rights reserved. 4 4 * … … 27 27 */ 28 28 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. 31 35 */ 32 /** @file 33 */ 36 int __signbit(size_t sz, ...) 37 { 38 va_list ap; 39 va_start(ap, sz); 34 40 35 #include <math.h> 36 #include <tan.h> 41 int result; 37 42 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; 50 57 } 51 58 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 1 1 /* 2 * Copyright (c) 201 5 Jiri Svoboda2 * Copyright (c) 2014 Martin Decky 3 3 * All rights reserved. 4 4 * … … 34 34 35 35 #include <math.h> 36 #include <fabs.h>37 36 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) 37 float fabsf(float val) 48 38 { 49 if (arg < 0.0) 50 return -arg; 51 else 52 return arg; 39 return copysignf(val, 1.0f); 53 40 } 54 41 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) 42 double fabs(double val) 65 43 { 66 if (arg < 0.0) 67 return -arg; 68 else 69 return arg; 44 return copysign(val, 1.0); 70 45 } 71 46 -
uspace/lib/math/generic/fmod.c
r7f7d642 r516e780 33 33 */ 34 34 35 #include <fmod.h>36 35 #include <math.h> 37 36 … … 52 51 * 53 52 */ 54 float 32_t float32_fmod(float32_t dividend, float32_t divisor)53 float fmodf(float dividend, float divisor) 55 54 { 56 55 // FIXME: replace with exact arithmetics 57 56 58 float 32_t quotient = trunc_f32(dividend / divisor);57 float quotient = truncf(dividend / divisor); 59 58 60 59 return (dividend - quotient * divisor); … … 77 76 * 78 77 */ 79 float64_t float64_fmod(float64_t dividend, float64_tdivisor)78 double fmod(double dividend, double divisor) 80 79 { 81 80 // FIXME: replace with exact arithmetics 82 81 83 float64_t quotient = trunc_f64(dividend / divisor);82 double quotient = trunc(dividend / divisor); 84 83 85 84 return (dividend - quotient * divisor); -
uspace/lib/math/generic/round.c
r7f7d642 r516e780 35 35 36 36 #include <math.h> 37 #include <pow.h> 37 #include <float.h> 38 #include <stdint.h> 38 39 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. 48 44 */ 49 float 32_t float32_pow(float32_t x, float32_t y)45 float roundf(float val) 50 46 { 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); 53 67 } 54 68 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. 64 73 */ 65 float64_t float64_pow(float64_t x, float64_t y)74 double round(double val) 66 75 { 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); 69 96 } 70 97 -
uspace/lib/math/generic/trig.c
r7f7d642 r516e780 35 35 36 36 #include <math.h> 37 #include <trig.h>38 37 39 38 #define TAYLOR_DEGREE_32 13 … … 41 40 42 41 /** Precomputed values for factorial (starting from 1!) */ 43 static float64_tfactorials[TAYLOR_DEGREE_64] = {42 static double factorials[TAYLOR_DEGREE_64] = { 44 43 1, 2, 6, 24, 120, 720, 5040, 40320, 362880, 3628800, 39916800, 45 44 479001600, 6227020800.0L, 87178291200.0L, 1307674368000.0L, … … 60 59 * 61 60 */ 62 static float 32_t taylor_sin_32(float32_t arg)63 { 64 float 32_tret = 0;65 float 32_tnom = 1;61 static float taylor_sin_32(float arg) 62 { 63 float ret = 0; 64 float nom = 1; 66 65 67 66 for (unsigned int i = 0; i < TAYLOR_DEGREE_32; i++) { … … 89 88 * 90 89 */ 91 static float64_t taylor_sin_64(float64_targ)92 { 93 float64_tret = 0;94 float64_tnom = 1;90 static double taylor_sin_64(double arg) 91 { 92 double ret = 0; 93 double nom = 1; 95 94 96 95 for (unsigned int i = 0; i < TAYLOR_DEGREE_64; i++) { … … 118 117 * 119 118 */ 120 static float 32_t taylor_cos_32(float32_t arg)121 { 122 float 32_tret = 1;123 float 32_tnom = 1;119 static float taylor_cos_32(float arg) 120 { 121 float ret = 1; 122 float nom = 1; 124 123 125 124 for (unsigned int i = 0; i < TAYLOR_DEGREE_32; i++) { … … 147 146 * 148 147 */ 149 static float64_t taylor_cos_64(float64_targ)150 { 151 float64_tret = 1;152 float64_tnom = 1;148 static double taylor_cos_64(double arg) 149 { 150 double ret = 1; 151 double nom = 1; 153 152 154 153 for (unsigned int i = 0; i < TAYLOR_DEGREE_64; i++) { … … 176 175 * 177 176 */ 178 static float 32_t base_sin_32(float32_t arg)177 static float base_sin_32(float arg) 179 178 { 180 179 unsigned int period = arg / (M_PI / 4); … … 209 208 * 210 209 */ 211 static float64_t base_sin_64(float64_targ)210 static double base_sin_64(double arg) 212 211 { 213 212 unsigned int period = arg / (M_PI / 4); … … 242 241 * 243 242 */ 244 static float 32_t base_cos_32(float32_t arg)243 static float base_cos_32(float arg) 245 244 { 246 245 unsigned int period = arg / (M_PI / 4); … … 275 274 * 276 275 */ 277 static float64_t base_cos_64(float64_targ)276 static double base_cos_64(double arg) 278 277 { 279 278 unsigned int period = arg / (M_PI / 4); … … 305 304 * 306 305 */ 307 float 32_t float32_sin(float32_t arg)308 { 309 float 32_t base_arg = fmod_f32(arg, 2 * M_PI);306 float sinf(float arg) 307 { 308 float base_arg = fmodf(arg, 2 * M_PI); 310 309 311 310 if (base_arg < 0) … … 324 323 * 325 324 */ 326 float64_t float64_sin(float64_targ)327 { 328 float64_t base_arg = fmod_f64(arg, 2 * M_PI);325 double sin(double arg) 326 { 327 double base_arg = fmod(arg, 2 * M_PI); 329 328 330 329 if (base_arg < 0) … … 343 342 * 344 343 */ 345 float 32_t float32_cos(float32_t arg)346 { 347 float 32_t base_arg = fmod_f32(arg, 2 * M_PI);344 float cosf(float arg) 345 { 346 float base_arg = fmodf(arg, 2 * M_PI); 348 347 349 348 if (base_arg < 0) … … 362 361 * 363 362 */ 364 float64_t float64_cos(float64_targ)365 { 366 float64_t base_arg = fmod_f64(arg, 2 * M_PI);363 double cos(double arg) 364 { 365 double base_arg = fmod(arg, 2 * M_PI); 367 366 368 367 if (base_arg < 0) -
uspace/lib/math/generic/trunc.c
r7f7d642 r516e780 34 34 */ 35 35 36 #include <mathtypes.h> 37 #include <trunc.h> 36 #include <math.h> 37 #include <float.h> 38 #include <stdint.h> 38 39 39 40 /** Truncate fractional part (round towards zero) … … 52 53 * 53 54 */ 54 float 32_t float32_trunc(float32_t val)55 float truncf(float val) 55 56 { 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; 58 60 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) }; 61 65 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; 71 67 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); 77 71 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); 79 78 } 80 79 … … 94 93 * 95 94 */ 96 float64_t float64_trunc(float64_tval)95 double trunc(double val) 97 96 { 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; 100 100 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) }; 103 105 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; 113 107 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); 119 111 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); 121 118 } 122 119 -
uspace/lib/math/test/main.c
r7f7d642 r516e780 1 1 /* 2 * Copyright (c) 201 5 Jiri Svoboda2 * Copyright (c) 2014 Vojtech Horky 3 3 * All rights reserved. 4 4 * … … 27 27 */ 28 28 29 /** @addtogroup libmath 30 * @{ 31 */ 32 /** @file 33 */ 29 #include <stdio.h> 30 #include <pcut/pcut.h> 34 31 35 #ifndef LIBMATH_LOG_H_ 36 #define LIBMATH_LOG_H_ 32 PCUT_INIT; 37 33 38 #include <mathtypes.h> 34 PCUT_IMPORT(rounding); 39 35 40 extern float32_t float32_log(float32_t); 41 extern float64_t float64_log(float64_t); 36 PCUT_MAIN(); 42 37 43 #endif44 45 /** @}46 */
Note:
See TracChangeset
for help on using the changeset viewer.