Index: uspace/lib/math/generic/ceil.c =================================================================== --- uspace/lib/math/generic/ceil.c (revision bae1e1f6864fcff51386ca7dda92f703816ffeef) +++ uspace/lib/math/generic/ceil.c (revision 01cdd5ae7da116060201cbd32623d045eaf00d62) @@ -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 01cdd5ae7da116060201cbd32623d045eaf00d62) +++ uspace/lib/math/generic/exp.c (revision 01cdd5ae7da116060201cbd32623d045eaf00d62) @@ -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 01cdd5ae7da116060201cbd32623d045eaf00d62) @@ -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 01cdd5ae7da116060201cbd32623d045eaf00d62) +++ uspace/lib/math/generic/log.c (revision 01cdd5ae7da116060201cbd32623d045eaf00d62) @@ -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 01cdd5ae7da116060201cbd32623d045eaf00d62) +++ uspace/lib/math/generic/pow.c (revision 01cdd5ae7da116060201cbd32623d045eaf00d62) @@ -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 01cdd5ae7da116060201cbd32623d045eaf00d62) @@ -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 01cdd5ae7da116060201cbd32623d045eaf00d62) @@ -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)