Index: uspace/app/tester/float/float2.c =================================================================== --- uspace/app/tester/float/float2.c (revision bae1e1f6864fcff51386ca7dda92f703816ffeef) +++ uspace/app/tester/float/float2.c (revision 992ffa61925d507ddd9ca906e90fb41d1ddbaacb) @@ -1,4 +1,5 @@ /* * Copyright (c) 2014 Martin Decky + * Copyright (c) 2015 Jiri Svoboda * All rights reserved. * @@ -27,4 +28,5 @@ */ +#include #include #include @@ -32,10 +34,21 @@ #include "../tester.h" -#define OPERANDS 10 -#define PRECISION 100000000 +#define OPERANDS 10 +#define PRECISIONF 10000 +#define PRECISION 100000000 static double arguments[OPERANDS] = { 3.5, -2.1, 100.0, 50.0, -1024.0, 0.0, 768.3156, 1080.499999, -600.0, 1.0 }; + +static double arguments_exp[OPERANDS] = { + 3.5, -2.1, 50.0, 0.0, 1.0, 13.2, -1.1, -5.5, 0.1, -66.0 +}; + +static double arguments_log[OPERANDS] = { + 3.5, 100.0, 50.0, 768.3156, 1080.499999, 1.0, 66.0, + 2.718281828459045, 9.9, 0.001 +}; + static double results_ceil[OPERANDS] = { @@ -63,8 +76,54 @@ }; +static double results_log[OPERANDS] = { + 1.252762968495, 4.605170185988, 3.912023005428, 6.644200586236, + 6.985179175021, 0.000000000000, 4.189654742026, 1.000000000000, + 2.292534757141, -6.907755278982 +}; + +static double results_exp[OPERANDS] = { + 33.115451958692, 0.122456428253, 5184705528587072045056.0, + 1.000000000000, 2.718281828459, 540364.937246691552, 0.332871083698, + 0.004086771438, 1.105170918076, 0.000000000000 +}; + +static bool cmp_float(float a, float b) +{ + float r; + + /* XXX Need fabsf() */ + if (b < 1.0 / PRECISIONF && b > -1.0 / PRECISIONF) + r = a; + else + r = a / b - 1.0; + + /* XXX Need fabsf() */ + if (r < 0.0) + r = -r; + + return r < 1.0 / PRECISIONF; +} + +static bool cmp_double(double a, double b) +{ + double r; + + /* XXX Need fabs() */ + if (b < 1.0 / PRECISION && b > -1.0 / PRECISION) + r = a; + else + r = a / b - 1.0; + + /* XXX Need fabs() */ + if (r < 0.0) + r = -r; + + return r < 1.0 / PRECISION; +} + const char *test_float2(void) { bool fail = false; - + for (unsigned int i = 0; i < OPERANDS; i++) { double res = floor(arguments[i]); @@ -73,6 +132,6 @@ if (res_int != corr_int) { - TPRINTF("Double floor failed (%" PRId64 " != %" PRId64 - ", arg %u)\n", res_int, corr_int, i); + TPRINTF("Double precision floor failed (%" PRId64 + " != %" PRId64 ", arg %u)\n", res_int, corr_int, i); fail = true; } @@ -85,6 +144,6 @@ if (res_int != corr_int) { - TPRINTF("Double ceil failed (%" PRId64 " != %" PRId64 - ", arg %u)\n", res_int, corr_int, i); + TPRINTF("Double precision ceil failed (%" PRId64 + " != %" PRId64 ", arg %u)\n", res_int, corr_int, i); fail = true; } @@ -97,6 +156,6 @@ if (res_int != corr_int) { - TPRINTF("Double truncation failed (%" PRId64 " != %" PRId64 - ", arg %u)\n", res_int, corr_int, i); + TPRINTF("Double precisiontruncation failed (%" PRId64 + " != %" PRId64 ", arg %u)\n", res_int, corr_int, i); fail = true; } @@ -109,6 +168,6 @@ if (res_int != corr_int) { - TPRINTF("Double sine failed (%" PRId64 " != %" PRId64 - ", arg %u)\n", res_int, corr_int, i); + TPRINTF("Double precision sine failed (%" PRId64 + " != %" PRId64 ", arg %u)\n", res_int, corr_int, i); fail = true; } @@ -121,6 +180,46 @@ if (res_int != corr_int) { - TPRINTF("Double cosine failed (%" PRId64 " != %" PRId64 - ", arg %u)\n", res_int, corr_int, i); + TPRINTF("Double precision cosine failed (%" PRId64 + " != %" PRId64 ", arg %u)\n", res_int, corr_int, i); + fail = true; + } + } + + for (unsigned int i = 0; i < OPERANDS; i++) { + float res = logf(arguments_log[i]); + + if (!cmp_float(res, results_log[i])) { + TPRINTF("Single precision logarithm failed " + "(%lf != %lf, arg %u)\n", res, results_log[i], i); + fail = true; + } + } + + for (unsigned int i = 0; i < OPERANDS; i++) { + double res = log(arguments_log[i]); + + if (!cmp_double(res, results_log[i])) { + TPRINTF("Double precision logarithm failed " + "(%lf != %lf, arg %u)\n", res, results_log[i], i); + fail = true; + } + } + + for (unsigned int i = 0; i < OPERANDS; i++) { + float res = exp(arguments_exp[i]); + + if (!cmp_float(res, results_exp[i])) { + TPRINTF("Single precision exponential failed " + "(%lf != %lf, arg %u)\n", res, results_exp[i], i); + fail = true; + } + } + + for (unsigned int i = 0; i < OPERANDS; i++) { + double res = exp(arguments_exp[i]); + + if (!cmp_double(res, results_exp[i])) { + TPRINTF("Double precision exponential failed " + "(%lf != %lf, arg %u)\n", res, results_exp[i], i); fail = true; } Index: uspace/lib/math/Makefile =================================================================== --- uspace/lib/math/Makefile (revision bae1e1f6864fcff51386ca7dda92f703816ffeef) +++ uspace/lib/math/Makefile (revision 992ffa61925d507ddd9ca906e90fb41d1ddbaacb) @@ -42,6 +42,9 @@ GENERIC_SOURCES = \ generic/ceil.c \ + generic/exp.c \ generic/floor.c \ + generic/log.c \ generic/trig.c \ + generic/pow.c \ generic/mod.c \ generic/trunc.c Index: uspace/lib/math/arch/amd64/include/libarch/math.h =================================================================== --- uspace/lib/math/arch/amd64/include/libarch/math.h (revision bae1e1f6864fcff51386ca7dda92f703816ffeef) +++ uspace/lib/math/arch/amd64/include/libarch/math.h (revision 992ffa61925d507ddd9ca906e90fb41d1ddbaacb) @@ -37,7 +37,10 @@ #include +#include #include +#include #include #include +#include static inline float64_t fmod(float64_t dividend, float64_t divisor) @@ -60,4 +63,14 @@ } +static inline float32_t expf(float32_t val) +{ + return float32_exp(val); +} + +static inline float64_t exp(float64_t val) +{ + return float64_exp(val); +} + static inline float64_t floor(float64_t val) { @@ -71,4 +84,24 @@ } +static inline float32_t logf(float32_t val) +{ + return float32_log(val); +} + +static inline float64_t log(float64_t val) +{ + return float64_log(val); +} + +static inline float32_t powf(float32_t x, float32_t y) +{ + return float32_pow(x, y); +} + +static inline float64_t pow(float64_t x, float64_t y) +{ + return float64_pow(x, y); +} + extern float64_t trunc(float64_t); Index: uspace/lib/math/generic/ceil.c =================================================================== --- uspace/lib/math/generic/ceil.c (revision bae1e1f6864fcff51386ca7dda92f703816ffeef) +++ uspace/lib/math/generic/ceil.c (revision 992ffa61925d507ddd9ca906e90fb41d1ddbaacb) @@ -48,8 +48,8 @@ float64_u v; float64_u r; - + v.data = val; t.data = trunc_float64(val); - + if (val.parts.sign == 1 || v.val == t.val) { r = t; @@ -57,5 +57,5 @@ r.val = t.val + 1.0; } - + return r.data; } Index: uspace/lib/math/generic/exp.c =================================================================== --- uspace/lib/math/generic/exp.c (revision 992ffa61925d507ddd9ca906e90fb41d1ddbaacb) +++ uspace/lib/math/generic/exp.c (revision 992ffa61925d507ddd9ca906e90fb41d1ddbaacb) @@ -0,0 +1,166 @@ +/* + * Copyright (c) 2015 Jiri Svoboda + * Copyright (c) 2014 Martin Decky + * All rights reserved. + * + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions + * are met: + * + * - Redistributions of source code must retain the above copyright + * notice, this list of conditions and the following disclaimer. + * - Redistributions in binary form must reproduce the above copyright + * notice, this list of conditions and the following disclaimer in the + * documentation and/or other materials provided with the distribution. + * - The name of the author may not be used to endorse or promote products + * derived from this software without specific prior written permission. + * + * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR + * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES + * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. + * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT, + * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT + * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, + * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY + * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT + * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF + * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + */ + +/** @addtogroup libmath + * @{ + */ +/** @file + */ + +#include +#include +#include + +#define TAYLOR_DEGREE_32 13 +#define TAYLOR_DEGREE_64 21 + +/** Precomputed values for factorial (starting from 1!) */ +static float64_t factorials[TAYLOR_DEGREE_64] = { + 1, 2, 6, 24, 120, 720, 5040, 40320, 362880, 3628800, 39916800, + 479001600, 6227020800.0L, 87178291200.0L, 1307674368000.0L, + 20922789888000.0L, 355687428096000.0L, 6402373705728000.0L, + 121645100408832000.0L, 2432902008176640000.0L, 51090942171709440000.0L +}; + +/** Exponential approximation by Taylor series (32-)) + * + * Compute the approximation of exponential by a Taylor + * series (using the first TAYLOR_DEGREE terms). + * The approximation is reasonably accurate for + * arguments within the interval XXXX. + * + * @param arg Argument. + * + * @return Exponential value approximation. + * + */ +static float64_t taylor_exp_32(float64_t arg) +{ + float64_t ret = 1; + float64_t nom = 1; + + for (unsigned int i = 0; i < TAYLOR_DEGREE_32; i++) { + nom *= arg; + ret += nom / factorials[i]; + } + + return ret; +} + +/** Exponential approximation by Taylor series + * + * Compute the approximation of exponential by a Taylor + * series (using the first TAYLOR_DEGREE terms). + * The approximation is reasonably accurate for + * arguments within the interval XXXX. + * + * @param arg Argument. + * + * @return Exponential value approximation. + * + */ +static float64_t taylor_exp_64(float64_t arg) +{ + float64_t ret = 1; + float64_t nom = 1; + + for (unsigned int i = 0; i < TAYLOR_DEGREE_64; i++) { + nom *= arg; + ret += nom / factorials[i]; + } + + return ret; +} + +/** Single precision exponential + * + * Compute exponential value. + * + * @param arg Exponential argument. + * + * @return Exponential value. + * + */ +float32_t float32_exp(float32_t arg) +{ + float32_t f; + float32_u i; + float32_u m; + float32_u r; + + /* + * e^a = (2 ^ log2(e))^a = 2 ^ (log2(e) * a) + * log2(e) * a = i + f | f in [0, 1] + * e ^ a = 2 ^ (i + f) = 2^f * 2^i = (e ^ log(2))^f * 2^i = + * e^(log(2)*f) * 2^i + */ + + m.val = arg * M_LOG2E; + i.data = trunc_float32(m.data); + f = arg * M_LOG2E - i.val; + + r.val = taylor_exp_32(M_LN2 * f); + r.data.parts.exp += i.val; + return r.val; +} + +/** Double precision exponential + * + * Compute exponential value. + * + * @param arg Exponential argument. + * + * @return Exponential value. + * + */ +float64_t float64_exp(float64_t arg) +{ + float64_t f; + float64_u i; + float64_u m; + float64_u r; + + /* + * e^a = (2 ^ log2(e))^a = 2 ^ (log2(e) * a) + * log2(e) * a = i + f | f in [0, 1] + * e ^ a = 2 ^ (i + f) = 2^f * 2^i = (e ^ log(2))^f * 2^i = + * e^(log(2)*f) * 2^i + */ + + m.val = arg * M_LOG2E; + i.data = trunc_float64(m.data); + f = arg * M_LOG2E - i.val; + + r.val = taylor_exp_64(M_LN2 * f); + r.data.parts.exp += i.val; + return r.val; +} + +/** @} + */ Index: uspace/lib/math/generic/floor.c =================================================================== --- uspace/lib/math/generic/floor.c (revision bae1e1f6864fcff51386ca7dda92f703816ffeef) +++ uspace/lib/math/generic/floor.c (revision 992ffa61925d507ddd9ca906e90fb41d1ddbaacb) @@ -49,8 +49,8 @@ float64_u v; float64_u r; - + v.data = val; t.data = trunc_float64(val); - + if (val.parts.sign == 0 || v.val == t.val) { r = t; @@ -58,5 +58,5 @@ r.val = t.val - 1.0; } - + return r.data; } Index: uspace/lib/math/generic/log.c =================================================================== --- uspace/lib/math/generic/log.c (revision 992ffa61925d507ddd9ca906e90fb41d1ddbaacb) +++ uspace/lib/math/generic/log.c (revision 992ffa61925d507ddd9ca906e90fb41d1ddbaacb) @@ -0,0 +1,160 @@ +/* + * Copyright (c) 2015 Jiri Svoboda + * Copyright (c) 2014 Martin Decky + * All rights reserved. + * + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions + * are met: + * + * - Redistributions of source code must retain the above copyright + * notice, this list of conditions and the following disclaimer. + * - Redistributions in binary form must reproduce the above copyright + * notice, this list of conditions and the following disclaimer in the + * documentation and/or other materials provided with the distribution. + * - The name of the author may not be used to endorse or promote products + * derived from this software without specific prior written permission. + * + * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR + * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES + * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. + * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT, + * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT + * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, + * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY + * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT + * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF + * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + */ + +/** @addtogroup libmath + * @{ + */ +/** @file + */ + +#include +#include + +#define TAYLOR_DEGREE_32 31 +#define TAYLOR_DEGREE_64 63 + +/** Single precision log(1 - arg) approximation by Taylor series + * + * Compute the approximation of log(1 - arg) by a Taylor + * series (using the first TAYLOR_DEGREE terms). + * arg must be within [-1, 1]. + * + * @param arg Argument. + * + * @return log(1 - arg) + * + */ +static float32_t taylor_log_32(float32_t arg) +{ + float32_t ret = 0; + float32_t num = 1; + + for (unsigned int i = 1; i <= TAYLOR_DEGREE_32; i++) { + num *= arg; + + if ((i % 2) == 0) + ret += num / i; + else + ret -= num / i; + } + + return ret; +} + +/** Double precision log(1 - arg) approximation by Taylor series + * + * Compute the approximation of log(1 - arg) by a Taylor + * series (using the first TAYLOR_DEGREE terms). + * arg must be within [-1, 1]. + * + * @param arg Argument. + * + * @return log(1 - arg) + * + */ +static float64_t taylor_log_64(float64_t arg) +{ + float64_t ret = 0; + float64_t num = 1; + + for (unsigned int i = 1; i <= TAYLOR_DEGREE_64; i++) { + num *= arg; + + if ((i % 2) == 0) + ret += num / i; + else + ret -= num / i; + } + + return ret; +} + +/** Single precision logarithm + * + * @param arg Argument. + * + * @return Logarithm. + * + */ +float32_t float32_log(float32_t arg) +{ + float32_u m; + int e; + + m.val = arg; + /* + * Factor arg into m * 2^e where m has exponent -1, + * which means it is in [1.0000..e-1, 1.1111..e-1] = [0.5, 1.0] + * so the argument to taylor_log_32 will be in [0, 0.5] + * ensuring that we get at least one extra bit of precision + * in each iteration. + */ + e = m.data.parts.exp - (FLOAT32_BIAS - 1); + m.data.parts.exp = FLOAT32_BIAS - 1; + + /* + * arg = m * 2^e ; log(arg) = log(m) + log(2^e) = + * log(m) + log2(2^e) / log2(e) = log(m) + e / log2(e) + */ + return - taylor_log_32(m.val - 1.0) + e / M_LOG2E; +} + +/** Double precision logarithm + * + * @param arg Argument. + * + * @return Logarithm. + * + */ +float64_t float64_log(float64_t arg) +{ + float64_u m; + int e; + + m.val = arg; + + /* + * Factor arg into m * 2^e where m has exponent -1, + * which means it is in [1.0000..e-1, 1.1111..e-1] = [0.5, 1.0] + * so the argument to taylor_log_32 will be in [0, 0.5] + * ensuring that we get at least one extra bit of precision + * in each iteration. + */ + e = m.data.parts.exp - (FLOAT64_BIAS - 1); + m.data.parts.exp = FLOAT64_BIAS - 1; + + /* + * arg = m * 2^e ; log(arg) = log(m) + log(2^e) = + * log(m) + log2(2^e) / log2(e) = log(m) + e / log2(e) + */ + return - taylor_log_64(m.val - 1.0) + e / M_LOG2E; +} + +/** @} + */ Index: uspace/lib/math/generic/pow.c =================================================================== --- uspace/lib/math/generic/pow.c (revision 992ffa61925d507ddd9ca906e90fb41d1ddbaacb) +++ uspace/lib/math/generic/pow.c (revision 992ffa61925d507ddd9ca906e90fb41d1ddbaacb) @@ -0,0 +1,74 @@ +/* + * Copyright (c) 2015 Jiri Svoboda + * Copyright (c) 2014 Martin Decky + * All rights reserved. + * + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions + * are met: + * + * - Redistributions of source code must retain the above copyright + * notice, this list of conditions and the following disclaimer. + * - Redistributions in binary form must reproduce the above copyright + * notice, this list of conditions and the following disclaimer in the + * documentation and/or other materials provided with the distribution. + * - The name of the author may not be used to endorse or promote products + * derived from this software without specific prior written permission. + * + * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR + * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES + * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. + * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT, + * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT + * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, + * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY + * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT + * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF + * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + */ + +/** @addtogroup libmath + * @{ + */ +/** @file + */ + +#include +#include +#include +#include + +/** Single precision power + * + * Compute power value. + * + * @param x Base + * @param y Exponent + * + * @return Cosine value. + * + */ +float32_t float32_pow(float32_t x, float32_t y) +{ + /* x^y = (e ^ log(x))^y = e ^ (log(x) * y) */ + return float32_exp(float32_log(x) * y); +} + +/** Double precision power + * + * Compute power value. + * + * @param x Base + * @param y Exponent + * + * @return Cosine value. + * + */ +float64_t float64_pow(float64_t x, float64_t y) +{ + /* x^y = (e ^ log(x))^y = e ^ (log(x) * y) */ + return float64_exp(float64_log(x) * y); +} + +/** @} + */ Index: uspace/lib/math/generic/trig.c =================================================================== --- uspace/lib/math/generic/trig.c (revision bae1e1f6864fcff51386ca7dda92f703816ffeef) +++ uspace/lib/math/generic/trig.c (revision 992ffa61925d507ddd9ca906e90fb41d1ddbaacb) @@ -36,10 +36,13 @@ #include -#define TAYLOR_DEGREE 13 +#define TAYLOR_DEGREE_32 13 +#define TAYLOR_DEGREE_64 21 /** Precomputed values for factorial (starting from 1!) */ -static float64_t factorials[TAYLOR_DEGREE] = { +static float64_t factorials[TAYLOR_DEGREE_64] = { 1, 2, 6, 24, 120, 720, 5040, 40320, 362880, 3628800, 39916800, - 479001600, 6227020800 + 479001600, 6227020800.0L, 87178291200.0L, 1307674368000.0L, + 20922789888000.0L, 355687428096000.0L, 6402373705728000.0L, + 121645100408832000.0L, 2432902008176640000.0L, 51090942171709440000.0L }; @@ -61,5 +64,5 @@ float64_t nom = 1; - for (unsigned int i = 0; i < TAYLOR_DEGREE; i++) { + for (unsigned int i = 0; i < TAYLOR_DEGREE_64; i++) { nom *= arg; @@ -90,5 +93,5 @@ float64_t nom = 1; - for (unsigned int i = 0; i < TAYLOR_DEGREE; i++) { + for (unsigned int i = 0; i < TAYLOR_DEGREE_64; i++) { nom *= arg; Index: uspace/lib/math/generic/trunc.c =================================================================== --- uspace/lib/math/generic/trunc.c (revision bae1e1f6864fcff51386ca7dda92f703816ffeef) +++ uspace/lib/math/generic/trunc.c (revision 992ffa61925d507ddd9ca906e90fb41d1ddbaacb) @@ -35,4 +35,42 @@ #include #include + +/** Truncate fractional part (round towards zero) + * + * Truncate the fractional part of IEEE 754 single + * precision floating point number by zeroing fraction + * bits, effectively rounding the number towards zero + * to the nearest whole number. + * + * If the argument is infinity or NaN, an exception + * should be indicated. This is not implemented yet. + * + * @param val Floating point number. + * + * @return Number rounded towards zero. + * + */ +float32 trunc_float32(float32 val) +{ + int32_t exp = val.parts.exp - FLOAT32_BIAS; + + if (exp < 0) { + /* -1 < val < 1 => result is +0 or -0 */ + val.parts.exp = 0; + val.parts.fraction = 0; + } else if (exp >= FLOAT32_FRACTION_SIZE) { + if (exp == 1024) { + /* val is +inf, -inf or NaN => trigger an exception */ + // FIXME TODO + } + + /* All bits in val are relevant for the result */ + } else { + /* Truncate irrelevant fraction bits */ + val.parts.fraction &= ~(UINT32_C(0x007fffff) >> exp); + } + + return val; +} /** Truncate fractional part (round towards zero) Index: uspace/lib/math/include/exp.h =================================================================== --- uspace/lib/math/include/exp.h (revision 992ffa61925d507ddd9ca906e90fb41d1ddbaacb) +++ uspace/lib/math/include/exp.h (revision 992ffa61925d507ddd9ca906e90fb41d1ddbaacb) @@ -0,0 +1,46 @@ +/* + * Copyright (c) 2015 Jiri Svoboda + * All rights reserved. + * + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions + * are met: + * + * - Redistributions of source code must retain the above copyright + * notice, this list of conditions and the following disclaimer. + * - Redistributions in binary form must reproduce the above copyright + * notice, this list of conditions and the following disclaimer in the + * documentation and/or other materials provided with the distribution. + * - The name of the author may not be used to endorse or promote products + * derived from this software without specific prior written permission. + * + * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR + * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES + * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. + * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT, + * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT + * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, + * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY + * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT + * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF + * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + */ + +/** @addtogroup libmath + * @{ + */ +/** @file + */ + +#ifndef LIBMATH_EXP_H_ +#define LIBMATH_EXP_H_ + +#include + +extern float32_t float32_exp(float32_t); +extern float64_t float64_exp(float64_t); + +#endif + +/** @} + */ Index: uspace/lib/math/include/log.h =================================================================== --- uspace/lib/math/include/log.h (revision 992ffa61925d507ddd9ca906e90fb41d1ddbaacb) +++ uspace/lib/math/include/log.h (revision 992ffa61925d507ddd9ca906e90fb41d1ddbaacb) @@ -0,0 +1,46 @@ +/* + * Copyright (c) 2015 Jiri Svoboda + * All rights reserved. + * + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions + * are met: + * + * - Redistributions of source code must retain the above copyright + * notice, this list of conditions and the following disclaimer. + * - Redistributions in binary form must reproduce the above copyright + * notice, this list of conditions and the following disclaimer in the + * documentation and/or other materials provided with the distribution. + * - The name of the author may not be used to endorse or promote products + * derived from this software without specific prior written permission. + * + * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR + * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES + * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. + * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT, + * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT + * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, + * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY + * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT + * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF + * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + */ + +/** @addtogroup libmath + * @{ + */ +/** @file + */ + +#ifndef LIBMATH_LOG_H_ +#define LIBMATH_LOG_H_ + +#include + +extern float32_t float32_log(float32_t); +extern float64_t float64_log(float64_t); + +#endif + +/** @} + */ Index: uspace/lib/math/include/math.h =================================================================== --- uspace/lib/math/include/math.h (revision bae1e1f6864fcff51386ca7dda92f703816ffeef) +++ uspace/lib/math/include/math.h (revision 992ffa61925d507ddd9ca906e90fb41d1ddbaacb) @@ -38,5 +38,7 @@ #include -#define M_PI 3.14159265358979323846 +#define M_LN2 0.69314718055994530942 +#define M_LOG2E 1.4426950408889634074 +#define M_PI 3.14159265358979323846 #endif Index: uspace/lib/math/include/pow.h =================================================================== --- uspace/lib/math/include/pow.h (revision 992ffa61925d507ddd9ca906e90fb41d1ddbaacb) +++ uspace/lib/math/include/pow.h (revision 992ffa61925d507ddd9ca906e90fb41d1ddbaacb) @@ -0,0 +1,46 @@ +/* + * Copyright (c) 2015 Jiri Svoboda + * All rights reserved. + * + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions + * are met: + * + * - Redistributions of source code must retain the above copyright + * notice, this list of conditions and the following disclaimer. + * - Redistributions in binary form must reproduce the above copyright + * notice, this list of conditions and the following disclaimer in the + * documentation and/or other materials provided with the distribution. + * - The name of the author may not be used to endorse or promote products + * derived from this software without specific prior written permission. + * + * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR + * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES + * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. + * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT, + * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT + * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, + * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY + * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT + * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF + * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + */ + +/** @addtogroup libmath + * @{ + */ +/** @file + */ + +#ifndef LIBMATH_POW_H_ +#define LIBMATH_POW_H_ + +#include + +extern float32_t float32_pow(float32_t, float32_t); +extern float64_t float64_pow(float64_t, float64_t); + +#endif + +/** @} + */ Index: uspace/lib/math/include/trunc.h =================================================================== --- uspace/lib/math/include/trunc.h (revision bae1e1f6864fcff51386ca7dda92f703816ffeef) +++ uspace/lib/math/include/trunc.h (revision 992ffa61925d507ddd9ca906e90fb41d1ddbaacb) @@ -38,4 +38,5 @@ #include +extern float32 trunc_float32(float32); extern float64 trunc_float64(float64);