From 331c34ecbdd5af55b2d00a4510d0a7cc0157b103 Mon Sep 17 00:00:00 2001 From: "Gavin D. Howard" Date: Sat, 3 Jun 2023 22:53:35 -0600 Subject: Tweak comments Signed-off-by: Gavin D. Howard --- src/num.c | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/num.c b/src/num.c index 0a597072..e45aa62d 100644 --- a/src/num.c +++ b/src/num.c @@ -2929,14 +2929,14 @@ exit: #endif // BC_ENABLE_EXTRA_MATH /** - * Converts a number from limbs with base BC_BASE_POW to base @a pow, where - * @a pow is obase^N. + * Takes a number with limbs with base BC_BASE_POW and converts the limb at the + * given index to base @a pow, where @a pow is obase^N. * @param n The number to convert. * @param rem BC_BASE_POW - @a pow. * @param pow The power of obase we will convert the number to. * @param idx The index of the number to start converting at. Doing the * conversion is O(n^2); we have to sweep through starting at the - * least significant limb + * least significant limb. */ static void bc_num_printFixup(BcNum* restrict n, BcBigDig rem, BcBigDig pow, size_t idx) @@ -2998,8 +2998,8 @@ bc_num_printFixup(BcNum* restrict n, BcBigDig rem, BcBigDig pow, size_t idx) } /** - * Prepares a number for printing in a base that is not a divisor of - * BC_BASE_POW. This basically converts the number from having limbs of base + * Prepares a number for printing in a base that does not have BC_BASE_POW as a + * power. This basically converts the number from having limbs of base * BC_BASE_POW to limbs of pow, where pow is obase^N. * @param n The number to prepare for printing. * @param rem The remainder of BC_BASE_POW when divided by a power of the base. -- cgit v1.2.3 From 454be80933ed8a1e9bf90ddb316a55d6078fe841 Mon Sep 17 00:00:00 2001 From: "Gavin D. Howard" Date: Wed, 14 Jun 2023 23:57:16 -0600 Subject: Tweak fuzz prep Signed-off-by: Gavin D. Howard --- scripts/fuzz_prep.sh | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/scripts/fuzz_prep.sh b/scripts/fuzz_prep.sh index 2c57d94a..1198b7f4 100755 --- a/scripts/fuzz_prep.sh +++ b/scripts/fuzz_prep.sh @@ -83,6 +83,11 @@ if [ "$asan" -ne 0 ]; then CFLAGS="$CFLAGS -fsanitize=address" fi +# These are to get better instrumentation. +export AFL_LLVM_LAF_SPLIT_SWITCHES=1 +export AFL_LLVM_LAF_TRANSFORM_COMPARES=1 +export AFL_LLVM_LAF_SPLIT_COMPARES=1 + # We want a debug build because asserts are counted as crashes too. CC="$CC" CFLAGS="$CFLAGS" ./configure.sh -gO3 -z -- cgit v1.2.3 From 75cf2e3358b5a76780db0d448c02cf6f61065921 Mon Sep 17 00:00:00 2001 From: "Gavin D. Howard" Date: Sat, 29 Jul 2023 15:30:34 -0600 Subject: Fix a missed coverage file Signed-off-by: Gavin D. Howard --- .gitignore | 1 + Makefile.in | 2 +- 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/.gitignore b/.gitignore index ea0bf73e..d2b57124 100644 --- a/.gitignore +++ b/.gitignore @@ -70,6 +70,7 @@ perf.data.old *.gcno *.gcov *.html +*.css *.profraw core.* diff --git a/Makefile.in b/Makefile.in index 55e2e4a6..e1309cd6 100644 --- a/Makefile.in +++ b/Makefile.in @@ -554,7 +554,7 @@ clean_config: clean clean_benchmarks clean_coverage: @printf 'Cleaning coverage files...\n' @$(RM) -f *.gcov - @$(RM) -f *.html + @$(RM) -f *.html *.css @$(RM) -f *.gcda *.gcno @$(RM) -f *.profraw @$(RM) -f $(GCDA) $(GCNO) -- cgit v1.2.3 From 81744cfd21ca386be78ef685b0f0317dc5c44e3d Mon Sep 17 00:00:00 2001 From: "Gavin D. Howard" Date: Fri, 18 Aug 2023 00:37:46 -0600 Subject: Add a better p(x, y) from TediusTimmy in GitHub issue #69 This new algorithm gives much better results in many cases. I'll document in a later commit. Signed-off-by: Gavin D. Howard --- gen/lib2.bc | 28 ++++++++++++++++++++++++++-- tests/bc/lib2.txt | 1 + tests/bc/lib2_results.txt | 4 +++- 3 files changed, 30 insertions(+), 3 deletions(-) diff --git a/gen/lib2.bc b/gen/lib2.bc index ba3f76b1..2d508714 100644 --- a/gen/lib2.bc +++ b/gen/lib2.bc @@ -34,10 +34,34 @@ */ define p(x,y){ - auto a + auto a,i,s,z + if(y==0)return 1@scale + if(x==0){ + if(y>0)return 0 + return 1/0 + } a=y$ if(y==a)return(x^a)@scale - return e(y*l(x)) + z=0 + if(x<1){ + y=-y + a=-a + z=x + x=1/x + } + if(y<0){ + return e(y*l(x)) + } + i=x^a + s=scale + scale+=length(i) + if(z){ + x=1/z + i=x^a + } + i*=e((y-a)*l(x)) + scale=s + return i@scale } define r(x,p){ auto t,n diff --git a/tests/bc/lib2.txt b/tests/bc/lib2.txt index 0032da19..74e1256d 100644 --- a/tests/bc/lib2.txt +++ b/tests/bc/lib2.txt @@ -1,6 +1,7 @@ p(2, 8.0000) p(2, 8.0001) p(2, -8.0001) +p(1024,32.1) r(0, 0) r(0, 1) r(0, 100) diff --git a/tests/bc/lib2_results.txt b/tests/bc/lib2_results.txt index f0753aff..e5ddb516 100644 --- a/tests/bc/lib2_results.txt +++ b/tests/bc/lib2_results.txt @@ -1,6 +1,8 @@ 256.00000000000000000000 -256.01774518281640169821 +256.01774518281640171326 .00390597924876622489 +42719740718418201647900434123391042292054090447133055398940832156444\ +39451561281100045924173873151.99999999999999999999 0 0 0 -- cgit v1.2.3 From 84ab7982fa36d5cdb5dcb7e30f35608522310db7 Mon Sep 17 00:00:00 2001 From: "Gavin D. Howard" Date: Fri, 18 Aug 2023 01:06:25 -0600 Subject: Document the better pow algorithm by TediusTimmy Signed-off-by: Gavin D. Howard --- manuals/algorithms.md | 68 +++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 68 insertions(+) diff --git a/manuals/algorithms.md b/manuals/algorithms.md index 4d7a0edc..fa4c3511 100644 --- a/manuals/algorithms.md +++ b/manuals/algorithms.md @@ -193,6 +193,74 @@ The algorithm used is to use the formula `e(y*l(x))`. It has a complexity of `O(n^3)` because both `e()` and `l()` do. +However, there are details to this algorithm, described by the author, +TediusTimmy, in GitHub issue #69. + +First, check if the exponent is 0. If it is, return 1 at the appropriate +`scale`. + +Next, check if the number is 0. If so, check if the exponent is greater than +zero; if it is, return 0. If the exponent is less than 0, error (with a divide +by 0) because that is undefined. + +Next, check if the exponent is actually an integer, and if it is, use the +exponentiation operator. + +At the `z=0` line is the start of the meat of the new code. + +`z` is set to zero as a flag and as a value. What I mean by that will be clear +later. + +Then we check if the number is less than 0. If it is, we negate the exponent +(and the integer version of the exponent, which we calculated earlier to check +if it was an integer). We also save the number in `z`; being non-zero is a flag +for later and a value to be used. Then we store the reciprocal of the number in +itself. + +All of the above paragraph will not make sense unless you remember the +relationship `l(x) == -l(1/x)`; we negated the exponent, which is equivalent to +the negative sign in that relationship, and we took the reciprocal of the +number, which is equivalent to the reciprocal in the relationship. + +But what if the number is negative? We ignore that for now because we eventually +call `l(x)`, which will raise an error if `x` is negative. + +Now, we can keep going. + +If at this point, the exponent is negative, we need to use the original formula +(`e(y * l(x))`) and return that result because the result will go to zero +anyway. + +But if we did *not* return, we know the exponent is *not* negative, so we can +get clever. + +We then compute the integral portion of the power by computing the number to +power of the integral portion of the exponent. + +Then we have the most clever trick: we add the length of that integer power to +the `scale`. Why? Because this will ensure that the next part is calculated to +at least as many digits as should be in the integer *plus* any extra `scale` +that was wanted. + +Then we check `z`, which, if it is not zero, is the original value of the +number. If it is not zero, we need to take the take the reciprocal *again* +because now we have the correct `scale`. And we *also* have to calculate the +integer portion of the power again. + +Then we need to calculate the fractional portion of the number. We do this by +using the original formula, but we instead of calculating `e(y * l(x))`, we +calculate `e((y - a) * l(x))`, where `a` is the integer portion of `y`. It's +easy to see that `y - a` will be just the fractional portion of `y` (the +exponent), so this makes sense. + +But then we *multiply* it into the integer portion of the power. Why? Because +remember: we're dealing with an exponent and a power; the relationship is +`x^(y+z) == (x^y)*(x^z)`. + +So we multiply it into the integer portion of the power. + +Finally, we set the result to the `scale`. + ### Rounding (`bc` Math Library 2 Only) This is implemented in the function `r(x,p)`. -- cgit v1.2.3 From 30e671105414d48ba1007291db97572628e1e569 Mon Sep 17 00:00:00 2001 From: "Gavin D. Howard" Date: Tue, 22 Aug 2023 21:24:22 -0600 Subject: Make a tiny tweak to TediusTimmy's algorithm I just increased the scale a bit. Signed-off-by: Gavin D. Howard --- gen/lib2.bc | 2 +- manuals/algorithms.md | 8 ++++---- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/gen/lib2.bc b/gen/lib2.bc index 2d508714..1d7d48c6 100644 --- a/gen/lib2.bc +++ b/gen/lib2.bc @@ -54,7 +54,7 @@ define p(x,y){ } i=x^a s=scale - scale+=length(i) + scale+=length(i)+5 if(z){ x=1/z i=x^a diff --git a/manuals/algorithms.md b/manuals/algorithms.md index fa4c3511..e0a5e8a4 100644 --- a/manuals/algorithms.md +++ b/manuals/algorithms.md @@ -237,10 +237,10 @@ get clever. We then compute the integral portion of the power by computing the number to power of the integral portion of the exponent. -Then we have the most clever trick: we add the length of that integer power to -the `scale`. Why? Because this will ensure that the next part is calculated to -at least as many digits as should be in the integer *plus* any extra `scale` -that was wanted. +Then we have the most clever trick: we add the length of that integer power (and +a little extra) to the `scale`. Why? Because this will ensure that the next part +is calculated to at least as many digits as should be in the integer *plus* any +extra `scale` that was wanted. Then we check `z`, which, if it is not zero, is the original value of the number. If it is not zero, we need to take the take the reciprocal *again* -- cgit v1.2.3 From d465612195b194e620e1c12cac8db5da59e6ac7e Mon Sep 17 00:00:00 2001 From: "Gavin D. Howard" Date: Thu, 31 Aug 2023 22:38:52 -0600 Subject: Remove a done TODO Signed-off-by: Gavin D. Howard --- include/bcl.h | 3 --- 1 file changed, 3 deletions(-) diff --git a/include/bcl.h b/include/bcl.h index 0908e215..d3a9f42c 100644 --- a/include/bcl.h +++ b/include/bcl.h @@ -36,9 +36,6 @@ #ifndef BC_BCL_H #define BC_BCL_H -// TODO: Add a generation index when building with Valgrind to check for -// use-after-free's or double frees. - #include #include #include -- cgit v1.2.3 From 86b7343c079caf96139a70e90e9a22b3747f14cb Mon Sep 17 00:00:00 2001 From: "Gavin D. Howard" Date: Thu, 31 Aug 2023 22:40:08 -0600 Subject: Increment the version and update the NEWS Signed-off-by: Gavin D. Howard --- NEWS.md | 7 +++++++ include/version.h | 2 +- 2 files changed, 8 insertions(+), 1 deletion(-) diff --git a/NEWS.md b/NEWS.md index de3b3502..28639ea5 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,5 +1,12 @@ # News +## 6.6.1 + +This is a production release with an improved `p()` function in the [extended +math library][16]. + +Users who don't care do not need to upgrade. + ## 6.6.0 This is a production release with two bug fixes and one change. diff --git a/include/version.h b/include/version.h index a4df383e..c7d2449c 100644 --- a/include/version.h +++ b/include/version.h @@ -37,6 +37,6 @@ #define BC_VERSION_H /// The current version. -#define VERSION 6.6.0 +#define VERSION 6.6.1 #endif // BC_VERSION_H -- cgit v1.2.3