about summary refs log tree commit diff
path: root/src/math
Commit message (Collapse)AuthorAgeFilesLines
* move lgamma-related internal declarations to libm.hRich Felker2018-09-124-12/+3
|
* add support for m68k 80-bit long double variantRich Felker2018-06-141-2/+10
| | | | | | | | | since x86 and m68k are the only archs with 80-bit long double and each has mandatory endianness, select the variant via endianness. differences are minor: apparently just byte order and representation of infinities. the m68k format is not well-documented anywhere I could find, so if other differences are found they may require additional changes later.
* fix fmaf wrong resultSzabolcs Nagy2018-04-021-1/+1
| | | | | | | | | | | | | if double precision r=x*y+z is not a half way case between two single precision floats or it is an exact result then fmaf returns (float)r. however the exactness check was wrong when |x*y| < |z| and could cause incorrectly rounded result in nearest rounding mode when r is a half way case. fmaf(-0x1.26524ep-54, -0x1.cb7868p+11, 0x1.d10f5ep-29) was incorrectly rounded up to 0x1.d117ap-29 instead of 0x1.d1179ep-29. (exact result is 0x1.d1179efffffffecp-29, r is 0x1.d1179fp-29)
* math: rewrite fma with mostly int arithmeticsSzabolcs Nagy2017-10-131-431/+154
| | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | the freebsd fma code failed to raise underflow exception in some cases in nearest rounding mode (affects fmal too) e.g. fma(-0x1p-1000, 0x1.000001p-74, 0x1p-1022) and the inexact exception may be raised spuriously since the fenv is not saved/restored around the exact multiplication algorithm (affects x86 fma too). another issue is that the underflow behaviour when the rounded result is the minimal normal number is target dependent, ieee754 allows two ways to raise underflow for inexact results: raise if the result before rounding is in the subnormal range (e.g. aarch64, arm, powerpc) or if the result after rounding with infinite exponent range is in the subnormal range (e.g. x86, mips, sh). to avoid all these issues the algorithm was rewritten with mostly int arithmetics and float arithmetics is only used to get correct rounding and raise exceptions according to the behaviour of the target without any fenv.h dependency. it also unifies x86 and non-x86 fma. fmaf is not affected, fmal need to be fixed too. this algorithm depends on a_clz_64 and it required a few spurious instructions to make sure underflow exception is raised in a particular corner case. (normally FORCE_EVAL(tiny*tiny) would be used for this, but on i386 gcc is broken if the expression is constant https://gcc.gnu.org/bugzilla/show_bug.cgi?id=57245 and there is no easy portable fix for the macro.)
* powerpc64: add single-instruction math functionsRich Felker2017-06-2322-0/+290
| | | | | | | | | | while the official elfv2 abi for "powerpc64le" sets power8 as the baseline isa, we use it for both little and big endian powerpc64 targets and need to maintain compatibility with pre-power8 models. the instructions for sqrt, fabs, and fma are in the baseline isa; support for the rest is conditional via predefined isa-level macros. patch by David Edelsohn.
* s390x: add single-instruction math functionsRich Felker2017-06-2324-0/+360
| | | | | | | | these were introduced in z196 and not available in the baseline (z900) ISA level. use __HTM__ as an alternate indicator for ISA level, since gcc did not define __ARCH__ until 7.x. patch by David Edelsohn.
* fix scalbn when result is in the subnormal rangeSzabolcs Nagy2017-04-213-12/+14
| | | | | | | | | | | | | | | | in nearest rounding mode scalbn could introduce double rounding error when an intermediate value and the final result were both in the subnormal range e.g. scalbn(0x1.7ffffffffffffp-1, -1073) returned 0x1p-1073 instead of 0x1p-1074, because the intermediate computation got rounded to 0x1.8p-1023. with the fix an intermediate value can only be in the subnormal range if the final result is 0 which is correct even after double rounding. (there still can be two roundings so signals may be raised twice, but that's only observable with trapping exceptions which is not supported.)
* aarch64: add single instruction math functionsSzabolcs Nagy2017-03-2134-24/+226
| | | | | | | | | | | | | | | | | | | | | | | | this should increase performance and reduce code size on aarch64. the compiled code was checked against using __builtin_* instead of inline asm with gcc-6.2.0. lrint is two instructions. c with inline asm is used because it is safer than a pure asm implementation, this prevents ll{rint,round} to be an alias of l{rint,round} (because the types don't match) and depends on gcc style inline asm support. ceil, floor, round, trunc can either raise inexact on finite non-integer inputs or not raise any exceptions. the new implementation does not raise exceptions while the generic c code does. on aarch64, the underflow exception is signaled before rounding (ieee 754 allows both before and after rounding, but it must be consistent), the generic fma c code signals it after rounding so using single instruction fixes a slight conformance issue too.
* fix threshold constants in j0f, y0f, j1f, y1fSzabolcs Nagy2017-03-152-13/+12
| | | | | | | | | | | | | partly following freebsd rev 279491 https://svnweb.freebsd.org/base?view=revision&revision=279491 (musl had some of the fixes before freebsd). the change should not matter much for j0f, y0f, but it improves j1f and y1f in [2.5,~3.75] (that is [0x40200000,~0x40700000]). near roots (e.g. around 3.8317 for j1f) there are still large ulp errors. dropped code that tried to raise inexact.
* math: fix pow signed shift ubSzabolcs Nagy2016-10-201-2/+2
| | | | | | | | | | | | j is int32_t and thus j<<31 is undefined if j==1, so j is changed to uint32_t locally as a quick fix, the generated code is not affected. (this is a strict conformance fix, future c standard may allow 1<<31, see DR 463. the bug was inherited from freebsd fdlibm, the proper fix is to use uint32_t for all bit hacks, but that requires more intrusive changes.) reported by Daniel Sabogal
* math: fix 128bit long double inverse trigonometric functionsSzabolcs Nagy2016-08-301-1/+1
| | | | | | | | | there was a copy paste error that could cause large ulp errors in atan2l, atanl, asinl and acosl on aarch64, mips64 and mipsn32. (the implementation is from freebsd fdlibm, but the tail end of the polynomial was wrong. 128 bit long double functions are not yet tested so this went undetected.)
* math: fix expf(-NAN) and exp2f(-NAN) to return -NAN instead of 0Szabolcs Nagy2016-03-042-0/+4
| | | | | | expf(-NAN) was treated as expf(-large) which unconditionally returns +0, so special case +-NAN. reported by Petr Hosek.
* work around regression building for armhf with clang (compiler bug)Rich Felker2016-02-192-2/+2
| | | | | | | | | commit e4355bd6bec89688e8c739cd7b4c76e675643dca moved the math asm from external source files to inline asm, but unfortunately, all current releases of clang use the wrong inline asm constraint codes for float and double ("w" and "P" instead of "t" and "w", respectively). this patch adds detection for the bug in configure, and, for now, just disables the affected asm on broken clang versions.
* improve macro logic for enabling arm math asmRich Felker2016-02-182-2/+2
| | | | | | | | | in order to take advantage of the fpu in -mfloat-abi=softfp mode, the __VFP_FP__ (presence of vfp fpu) was checked instead of checking for __ARM_PCS_VFP (hardfloat EABI variant). however, the latter macro is the one that's actually specified by the ABI documents rather than being compiler-specific, and should also be checked in case __VFP_FP__ is not defined on some compilers or some configurations.
* replace armhf math asm source files with inline asmRich Felker2016-01-2016-40/+60
| | | | | | | | | | | | this makes it possible to inline them with LTO, and is the simplest approach to eliminating the use of .sub files. this also makes VFP sqrt available for use with the standard EABI (plain arm rather than armhf subarch) when libc is built with -mfloat-abi=softfp. the same could have been done for fabs, but when the argument and return value are in integer registers, moving to VFP registers and back is almost certainly more costly than a simple integer operation.
* math: explicitly promote expressions to excess-precision typesRich Felker2015-11-213-4/+4
| | | | | | | | | | | | | | | | | | | | | | | | | | | | | a conforming compiler for an arch with excess precision floating point (FLT_EVAL_METHOD!=0; presently i386 is the only such arch supported) computes all intermediate results in the types float_t and double_t rather than the nominal type of the expression. some incorrect compilers, however, only keep excess precision in registers, and convert down to the nominal type when spilling intermediate results to memory, yielding unpredictable results that depend on the compiler's choices of what/when to spill. in particular, this happens on old gcc versions with -ffloat-store, which we need in order to work around bugs where the compiler wrongly keeps explicitly-dropped excess precision. by explicitly converting to double_t where expressions are expected be be evaluated in double_t precision, we can avoid depending on the compiler to get types correct when spilling; the nominal and intermediate precision now match. this commit should not change the code generated by correct compilers, or by old ones on non-i386 archs where double_t is defined as double. this fixes a serious bug in argument reduction observed on i386 with gcc 4.2: for values of x outside the unit circle, sin(x) was producing results outside the interval [-1,1]. changes made in commit 0ce946cf808274c2d6e5419b139e130c8ad4bd30 were likely responsible for breaking compatibility with this and other old gcc versions. patch by Szabolcs Nagy.
* explicitly assemble all arm asm sources as UALRich Felker2015-11-104-0/+4
| | | | | | | | these files are all accepted as legacy arm syntax when producing arm code, but legacy syntax cannot be used for producing thumb2 with access to the full ISA. even after switching to UAL, some asm source files contain instructions which are not valid in thumb mode, so these will need to be addressed separately.
* declare fpu usage to the assembler in arm hard-float asm filesSzabolcs Nagy2015-10-194-0/+4
| | | | | | | Some armhf gcc toolchains (built with --with-float=hard but without --with-fpu=vfp*) do not pass -mfpu=vfp to the assembler and then binutils rejects the UAL mnemonics for VFP unless there is an .fpu vfp directive in the asm source.
* fix regression in x86_64 math asm with old binutilsRich Felker2015-04-232-6/+6
| | | | | | | | | the implicit-operand form of fucomip is rejected by binutils 2.19 and perhaps other versions still in use. writing both operands explicitly fixes the issue. there is no change to the resulting output. commit a732e80d33b4fd6f510f7cec4f5573ef5d89bc4e was the source of this regression.
* remove potentially PIC-incompatible relocations from x86_64 and x32 asmRich Felker2015-04-182-2/+2
| | | | analogous to commit 8ed66ecbcba1dd0f899f22b534aac92a282f42d5 for i386.
* remove the last of possible-textrels from i386 asmRich Felker2015-04-182-1/+5
| | | | | | | | | | | | none of these are actual textrels because of ld-time binding performed by -Bsymbolic-functions, but I'm changing them with the goal of making ld-time binding purely an optimization rather than relying on it for semantic purposes. in the case of memmove's call to memcpy, making it explicit that the memmove asm is assuming the forward-copying behavior of the memcpy asm is desirable anyway; in case memcpy is ever changed, the semantic mismatch would be apparent while editing memmcpy.s.
* math: fix pow(+-0,-inf) not to raise divbyzero flagSzabolcs Nagy2015-04-183-3/+3
| | | | | | this reverts the commit f29fea00b5bc72d4b8abccba2bb1e312684d1fce which was based on a bug in C99 and POSIX and did not match IEEE-754 http://www.open-std.org/jtc1/sc22/wg14/www/docs/n1515.pdf
* add aarch64 portSzabolcs Nagy2015-03-114-0/+24
| | | | | | | | | | This adds complete aarch64 target support including bigendian subarch. Some of the long double math functions are known to be broken otherwise interfaces should be fully functional, but at this point consider this port experimental. Initial work on this port was done by Sireesh Tripurari and Kevin Bortis.
* math: add dummy implementations of 128 bit long double functionsSzabolcs Nagy2015-03-1116-4/+97
| | | | | | | | This is in preparation for the aarch64 port only to have the long double math symbols available on ld128 platforms. The implementations should be fixed up later once we have proper tests for these functions. Added bigendian handling for ld128 bit manipulations too.
* math: add ld128 exp2l based on the freebsd implementationSzabolcs Nagy2015-03-111-1/+366
| | | | | Changed the special case handling and bit manipulation to better match the double version.
* math: fix fmodl for IEEE binary128Szabolcs Nagy2015-02-091-1/+1
| | | | | This trivial copy-paste bug went unnoticed due to lack of testing. No currently supported target archs are affected.
* math: fix __fpclassifyl(-0.0) for IEEE binary128Szabolcs Nagy2015-02-081-3/+2
| | | | | The sign bit was not cleared before checking for 0 so -0.0 was misclassified as FP_SUBNORMAL instead of FP_ZERO.
* add parenthesis in fma.c to clarify intent and silence warningsSzabolcs Nagy2015-02-081-1/+1
|
* math: use fnstsw consistently instead of fstsw in x87 asmSzabolcs Nagy2014-11-0511-11/+11
| | | | | | fnstsw does not wait for pending unmasked x87 floating-point exceptions and it is the same as fstsw when all exceptions are masked which is the only environment libc supports.
* math: fix x86_64 and x32 asm not to use sahf instructionSzabolcs Nagy2014-11-056-28/+14
| | | | | | | | | | | Some early x86_64 cpus (released before 2006) did not support sahf/lahf instructions so they should be avoided (intel manual says they are only supported if CPUID.80000001H:ECX.LAHF-SAHF[bit 0] = 1). The workaround simplifies exp2l and expm1l because fucomip can be used instead of the fucomp;fnstsw;sahf sequence copied from i386. In fmodl and remainderl sahf is replaced by a simple bit test.
* math: use the rounding idiom consistentlySzabolcs Nagy2014-10-3113-58/+89
| | | | | | | | | | | | | | | | | | | | the idiomatic rounding of x is n = x + toint - toint; where toint is either 1/EPSILON (x is non-negative) or 1.5/EPSILON (x may be negative and nearest rounding mode is assumed) and EPSILON is according to the evaluation precision (the type of toint is not very important, because single precision float can represent the 1/EPSILON of ieee binary128). in case of FLT_EVAL_METHOD!=0 this avoids a useless store to double or float precision, and the long double code became cleaner with 1/LDBL_EPSILON instead of ifdefs for toint. __rem_pio2f and __rem_pio2 functions slightly changed semantics: on i386 a double-rounding is avoided so close to half-way cases may get evaluated differently eg. as sin(pi/4-eps) instead of cos(pi/4+eps)
* fix rint.c and rintf.c when FLT_EVAL_METHOD!=0Szabolcs Nagy2014-10-312-4/+22
| | | | | | | | | | | | | | | | | | | | The old code used the rounding idiom incorrectly: y = (double)(x + 0x1p52) - 0x1p52; the cast is useless if FLT_EVAL_METHOD==0 and causes a second rounding if FLT_EVAL_METHOD==2 which can give incorrect result in nearest rounding mode, so the correct idiom is to add/sub a power-of-2 according to the characteristics of double_t. This did not cause actual bug because only i386 is affected where rint is implemented in asm. Other rounding functions use a similar idiom, but they give correct results because they only rely on getting a neighboring integer result and the rounding direction is fixed up separately independently of the current rounding mode. However they should be fixed to use the idiom correctly too.
* always provide __fpclassifyl and __signbitl definitionsRich Felker2014-10-082-1/+9
| | | | | | | | | previously the external definitions of these functions were omitted on archs where long double is the same as double, since the code paths in the math.h macros which would call them are unreachable. however, even if they are unreachable, the definitions are still mandatory. omitting them is invalid C, and in the case of a non-optimizing compiler, will result in a link error.
* math: fix exp10 not to raise invalid exception on NaNSzabolcs Nagy2014-09-183-4/+13
| | | | | This was not caught earlier because gcc incorrectly generates quiet relational operators that never raise exceptions.
* fix exp10l.c to include float.hSzabolcs Nagy2014-09-081-0/+1
| | | | | | the previous commit was a no op in exp10l because LDBL_* macros were implicitly 0 (the preprocessor does not warn about undefined symbols).
* prune math code on archs with binary64 long doubleSzabolcs Nagy2014-09-082-0/+10
| | | | | | __polevll, __p1evll and exp10l were provided on archs when long double is the same as double. The first two were completely unused and exp10l can be a wrapper around exp10.
* math: fix aliasing violation in long double wrappersSzabolcs Nagy2014-04-112-2/+10
| | | | | | | | | modfl and sincosl were passing long double* instead of double* to the wrapped double precision functions (on archs where long double and double have the same size). This is fixed now by using temporaries (this is not optimized to a single branch so the generated code is a bit bigger). Found by Morten Welinder.
* x32 port (diff against vanilla x86_64)rofl0r2014-02-2318-69/+69
|
* import vanilla x86_64 code as x32rofl0r2014-02-2330-0/+396
|
* math: add drem and dremf weak aliases to i386 remainder asmSzabolcs Nagy2014-01-082-0/+6
| | | | | weak_alias was only in the c code, so drem was missing on platforms where remainder is implemented in asm.
* math: define _GNU_SOURCE when implementing non-standard math functionsSzabolcs Nagy2013-12-126-0/+6
| | | | | this makes the prototypes in math.h are visible so they are checked agaist the function definitions
* math: clean up __rem_pio2Szabolcs Nagy2013-11-243-71/+53
| | | | | | | | | | | - remove the HAVE_EFFICIENT_IRINT case: fn is an exact integer, so it can be converted to int32_t a bit more efficiently than with a cast (the rounding mode change can be avoided), but musl does not support this case on any arch. - __rem_pio2: use double_t where possible - __rem_pio2f: use less assignments to avoid stores on i386 - use unsigned int bit manipulation (and union instead of macros) - use hexfloat literals instead of named constants
* math: add (obsolete) bsd drem and finite functionsSzabolcs Nagy2013-11-214-0/+20
|
* math: lgamma cleanup (simpler sin(pi*x) for the negative case)Szabolcs Nagy2013-11-214-202/+110
| | | | | | | | | * simplify sin_pi(x) (don't care about inexact here, the result is inexact anyway, and x is not so small to underflow) * in lgammal add the previously removed special case for x==1 and x==2 (to fix the sign of zero in downward rounding mode) * only define lgammal on supported long double platforms * change tgamma so the generated code is a bit smaller
* math: extensive log*.c cleanupSzabolcs Nagy2013-10-2814-583/+369
| | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | The log, log2 and log10 functions share a lot of code and to a lesser extent log1p too. A small part of the code was kept separately in __log1p.h, but since it did not capture much of the common code and it was inlined anyway, it did not solve the issue properly. Now the log functions have significant code duplication, which may be resolved later, until then they need to be modified together. logl, log10l, log2l, log1pl: * Fix the sign when the return value should be -inf. * Remove the volatile hack from log10l (seems unnecessary) log1p, log1pf: * Change the handling of small inputs: only |x|<2^-53 is special (then it is enough to return x with the usual subnormal handling) this fixes the sign of log1p(0) in downward rounding. * Do not handle the k==0 case specially (other than skipping the elaborate argument reduction) * Do not handle 1+x close to power-of-two specially (this code was used rarely, did not give much speed up and the precision wasn't better than the general) * Fix the correction term formula (c=1-(u-x) was used incorrectly when x<1 but (double)(x+1)==2, this was not a critical issue) * Use the exact same method for calculating log(1+f) as in log (except in log1p the c correction term is added to the result). log, logf, log10, log10f, log2, log2f: * Use double_t and float_t consistently. * Now the first part of log10 and log2 is identical to log (until the return statement, hopefully this makes maintainence easier). * Most special case formulas were removed (close to power-of-two and k==0 cases), they increase the code size without providing precision or performance benefits (and obfuscate the code). Only x==1 is handled specially so in downward rounding mode the sign of zero is correct (the general formula happens to give -0). * For x==0 instead of -1/0.0 or -two54/0.0, return -1/(x*x) to force raising the exception at runtime. * Arg reduction code is changed (slightly simplified) * The thresholds for arg reduction to [sqrt(2)/2,sqrt(2)] are now consistently the [0x3fe6a09e00000000,0x3ff6a09dffffffff] and the [0x3f3504f3,0x3fb504f2] intervals for double and float reductions respectively (the exact threshold values are not critical) * Remove the obsolete comment for the FLT_EVAL_METHOD!=0 case in log2f (The same code is used for all eval methods now, on i386 slightly simpler code could be used, but we have asm there anyway) all: * Fix signed int arithmetics (using unsigned for bitmanipulation) * Fix various comments
* math: fix rare underflow issue in fmaSzabolcs Nagy2013-10-073-13/+55
| | | | | the issue is described in commits 1e5eb73545ca6cfe8b918798835aaf6e07af5beb and ffd8ac2dd50f99c3c83d7d9d845df9874ec3e7d5
* math: use sqrtl if FLT_EVAL_METHOD==2 in acosh and acoshfSzabolcs Nagy2013-10-072-0/+13
| | | | this makes acosh slightly more precise around 1.0 on i386
* math: remove an unused variable from modflSzabolcs Nagy2013-10-061-1/+0
|
* math: remove code duplication in erfl found by clang analyzerSzabolcs Nagy2013-10-041-13/+2
| | | | | erfl had some superflous code left around after the last erf cleanup. the issue was reported by Alexander Monakov
* math: remove a useless assignment in lgammal found by clang analyzerSzabolcs Nagy2013-10-041-2/+2
| | | | the issue was reported by Alexander Monakov
generated by cgit-pink 1.4.1 (git 2.36.1) at 2024-11-10 01:13:05 +0000