LibBF – a small library to handle arbitrary precision floating point numbers(bellard.org)
bellard.org
LibBF – a small library to handle arbitrary precision floating point numbers
https://bellard.org/libbf/
62 comments
Yes, a comparison with GMP is the first link: https://bellard.org/libbf/benchmark.html
Is there a comparison for anything other than multiplication? Or is multiplication fairly representative of other operations on AP numbers?
God damn, as if the others were not enough... Can someone confirm that this dude is a human?
Can confirm that he is an alien--we're just not sure if extraterrestrial or extradimensional.
He also built a home brew 3G and then LTE base station, iirc.
Yes, his company Amarisoft sells a SDR LTE based station with EPC (the core network part) support too:
https://www.amarisoft.com/
And he generated MPEG2 DVB streams (actual RF signal) from clever GPU coding.
Fabrice Bellard is so outrageously awesome, it makes me wonder why I bother to exist.
I like the stuff he did for software radio decoding too.
I like the stuff he did for software radio decoding too.
I've just noticed that one of the jslinux configurations now boots to X Windows(!) for either RISC-V(!) 32 or 64 bit.
GET OUT OF TOWN.
GET OUT OF TOWN.
Mostly off topic but maybe this post attracts people who know the answer to a question that I ran into when I implemented a big integer based rational number library some time ago.
The decimal representation of a rational number p / q has the general form ±<integer part>.<non-repeating fractional part>(<repeating fractional part>), for example 41,111,111 / 333,000 equals 123.456(789). Is there an efficient algorithm to determine the lengths of the fractional parts or one that does at least provides reasonably tight bounds on them?
I did quite a bit of research at the time but only with limited success. I am reasonably sure that finding the exact lengths is comparable in hardness to factoring p and q in some cases and possibly harder in the general case. But I did not come across much with regard to bounds on the lengths.
Maybe there was somewhere an implicit answer to my question which I did not recognize or fully appreciate due to my limited knowledge of number theory but I never came across an explicit statement that said that printing the decimal representation of a rational number or estimating its length is a hard or unsolved problem.
I ended up just performing the expansion digit by digit until I detected a cycle or reached a specified maximum length, but it would by nice to know in advance how long it will be or at least to learn that this is actually an hard or unsolved problem.
The decimal representation of a rational number p / q has the general form ±<integer part>.<non-repeating fractional part>(<repeating fractional part>), for example 41,111,111 / 333,000 equals 123.456(789). Is there an efficient algorithm to determine the lengths of the fractional parts or one that does at least provides reasonably tight bounds on them?
I did quite a bit of research at the time but only with limited success. I am reasonably sure that finding the exact lengths is comparable in hardness to factoring p and q in some cases and possibly harder in the general case. But I did not come across much with regard to bounds on the lengths.
Maybe there was somewhere an implicit answer to my question which I did not recognize or fully appreciate due to my limited knowledge of number theory but I never came across an explicit statement that said that printing the decimal representation of a rational number or estimating its length is a hard or unsolved problem.
I ended up just performing the expansion digit by digit until I detected a cycle or reached a specified maximum length, but it would by nice to know in advance how long it will be or at least to learn that this is actually an hard or unsolved problem.
The length of the repeated part is the multiplicative order of 10 mod q [1].
Look at the Wikipedia page on repeated decimals and you can see many of the repeated sequences have length (q-1) ((q-1) is phi(q) when q is prime) [2].
There's a SO answer that talks about this [3]. I think I can derive it here, though take it with a grain of salt as I might have screwed something up.
Say you have a rational number p / q, and w.l.o.g., p < q. Consider an integer, R, that corresponds to the repeated portion of the rational number (in base 10), with length s (in decimal). This implies:
So to make the above equation true, let's be coarse and choose s = phi(q) in the below:
There's also the case of when 10 and q aren't relatively prime that needs to be worked out.
You can also see how this argument could be modified to work with other bases.
[1] https://softwareengineering.stackexchange.com/a/192077
[2] https://en.wikipedia.org/wiki/Repeating_decimal#Decimal_expa...
[3] https://softwareengineering.stackexchange.com/a/192081
Look at the Wikipedia page on repeated decimals and you can see many of the repeated sequences have length (q-1) ((q-1) is phi(q) when q is prime) [2].
There's a SO answer that talks about this [3]. I think I can derive it here, though take it with a grain of salt as I might have screwed something up.
Say you have a rational number p / q, and w.l.o.g., p < q. Consider an integer, R, that corresponds to the repeated portion of the rational number (in base 10), with length s (in decimal). This implies:
p/q = R sum_{k=1}^{\infty} 10^{-s k} = R / ( 10^s - 1 )
p (10^s - 1) - q R = 0
So to make the above equation true, let's be coarse and choose s = phi(q) in the below:
.. p (10^{ phi(q) } - 1 ) - q R = 0
-> p ( C q ) - q R = 0
-> q ( p C - R ) = 0
There's also the case of when 10 and q aren't relatively prime that needs to be worked out.
You can also see how this argument could be modified to work with other bases.
[1] https://softwareengineering.stackexchange.com/a/192077
[2] https://en.wikipedia.org/wiki/Repeating_decimal#Decimal_expa...
[3] https://softwareengineering.stackexchange.com/a/192081
This made me look at multiplicative order again and I finally found [1] what I could not find at the time, someone explicitly stating that this is a hard problem. Now it seems so obvious that I should have looked at what cryptography has to say about it that I can not imagine why I did not. Anyway, case finally closed, I got my peace of mind now.
[1] https://math.stackexchange.com/questions/837489/algorithms-f...
[1] https://math.stackexchange.com/questions/837489/algorithms-f...
While I’m sure you know this, when computing the fractional part you will have at most q-1 digits in your repeating cycle.
When you divide by q, you have at most q unique remainders (0 to q-1). Also when computing the fractional part, if you end up with the same remainder as you’ve previously seen then you’ve hit the cycle.
That is exactly what I did, I implemented some cycle finding algorithm, a tortoise and hare variant maybe but I can not tell from the top of my head, on the reminder sequence giving up after a user-specified number of steps. But it irks me that this may abort just a hand full of steps before finally finding the cycle.
And unfortunately q is just way to loose as a bound if it even deserves the name bound in this case. I initially started looking into this because I wanted to print something like »0.577 215 664 <901,532 more digits>« but in that case q is essentially useless because you either have to really search the cycle which quickly becomes impractical once q reaches millions or billions or you have to live with »0.577 215 664 <up to 999,999,991 more digits but maybe also just one>«.
And unfortunately q is just way to loose as a bound if it even deserves the name bound in this case. I initially started looking into this because I wanted to print something like »0.577 215 664 <901,532 more digits>« but in that case q is essentially useless because you either have to really search the cycle which quickly becomes impractical once q reaches millions or billions or you have to live with »0.577 215 664 <up to 999,999,991 more digits but maybe also just one>«.
> I implemented some cycle finding algorithm, a tortoise and hare variant maybe but I can not tell from the top of my head, on the reminder sequence giving up after a user-specified number of steps.
What's the cycle finding for? Isn't it enough to record the remainder sequence with indices? As soon as any remainder reoccurs, that's the cycle.
What's the cycle finding for? Isn't it enough to record the remainder sequence with indices? As soon as any remainder reoccurs, that's the cycle.
You could of course just maintain the set of remainders seen so far and then you would detect the cycle by reaching a remainder that is already in your set of remainders seen so far but that requires space proportional to the number of digits expanded and if you are expanding a lot of digits and have a large denominator, this will add substantially to the amount of memory consumed. The digits themselves require of course a linear amount of space, too, but only about one byte per digit compared to the potentially big remainder.
There are some closely related difficult problems. Consider the fraction a/b, and let r be the length of the repeating fractional part. Then r <= b - 1. Moreover, if r == b - 1, then b is prime. In general, r divides phi(b), where phi is the Euler phi function.
As you said, this is equivalent to factorisation in terms of complexity. The Wikipedia article[0] on the subject does a good job of summarising the details of this in the section "Other properties of repeatend lengths".
For a simple answer, the lengths are bounded above by Euler's totient function [1] (the sum of factors) of the denominator. Any tighter bounds depend on the specific primes in the factorisation.
Note that Euler's totient function is "always 'nearly n'" in magnitude (Hardy and Littlewood), so that representation will be very expensive for most fractions of large denominator.
[0]: https://en.wikipedia.org/wiki/Repeating_decimal [1]: https://en.wikipedia.org/wiki/Euler%27s_totient_function
For a simple answer, the lengths are bounded above by Euler's totient function [1] (the sum of factors) of the denominator. Any tighter bounds depend on the specific primes in the factorisation.
Note that Euler's totient function is "always 'nearly n'" in magnitude (Hardy and Littlewood), so that representation will be very expensive for most fractions of large denominator.
[0]: https://en.wikipedia.org/wiki/Repeating_decimal [1]: https://en.wikipedia.org/wiki/Euler%27s_totient_function
> As you said, this is equivalent to factorisation in terms of complexity.
This can't be right. Factorisation is, legendarily, a hard problem, unsolved (in realistic time frames) for large input. But more than that, it's harder than division is.
Finding the length of the nonrepeating and repeating parts of a rational number's decimal expansion can be done just by doing the division. That would appear to be evidence that it is not harder than division?
Edit: from the wikipedia section you link:
> If k = 2^a·5^b·n where n > 1 and n is not divisible by 2 or 5, then the length of the transient of 1/k is max(a,b) and the period [of the repetend, presumably] is r, where r is the smallest integer such that 10^r ≡ 1 (mod n)
All integers (except those whose only prime factors are 2 and 5) satisfy this condition, so this would appear to be exactly the method desired? Repeatedly dividing by 5 is much, much easier than factoring is. Taking the discrete log of 1... could be equivalent to factoring.
This can't be right. Factorisation is, legendarily, a hard problem, unsolved (in realistic time frames) for large input. But more than that, it's harder than division is.
Finding the length of the nonrepeating and repeating parts of a rational number's decimal expansion can be done just by doing the division. That would appear to be evidence that it is not harder than division?
Edit: from the wikipedia section you link:
> If k = 2^a·5^b·n where n > 1 and n is not divisible by 2 or 5, then the length of the transient of 1/k is max(a,b) and the period [of the repetend, presumably] is r, where r is the smallest integer such that 10^r ≡ 1 (mod n)
All integers (except those whose only prime factors are 2 and 5) satisfy this condition, so this would appear to be exactly the method desired? Repeatedly dividing by 5 is much, much easier than factoring is. Taking the discrete log of 1... could be equivalent to factoring.
That would appear to be evidence that it is not harder than division?
From your edit I guess that you already realized this, but the hard part here is that for 1 / q you would have to perform up to q divisions of integers of size n = log₂(q) which of course takes time O(2ⁿ) even assuming the divisions are O(1).
From your edit I guess that you already realized this, but the hard part here is that for 1 / q you would have to perform up to q divisions of integers of size n = log₂(q) which of course takes time O(2ⁿ) even assuming the divisions are O(1).
That Wikipedia article on repeating decimals is probably where I started looking into it, but the article is somewhat weird. For many problems Wikipedia articles say we know this and that but those other things are still unknown, this article has a lot of information but does not talk about unknown things at all. The length of the repeating fractional part of 1 / n divides φ(n)...and then nothing. Is this all we know or is there more but just not in this article?
I still remember that quite well because that was probably the first and only time that I was unable to find someone stating that we do not know any better, and not only the Wikipedia article but everywhere I looked there was no such statement to be found. I eventually gave up after a couple of evenings but it obviously kept stuck somewhere in my mind and that is why I wrote the initial comment.
I still remember that quite well because that was probably the first and only time that I was unable to find someone stating that we do not know any better, and not only the Wikipedia article but everywhere I looked there was no such statement to be found. I eventually gave up after a couple of evenings but it obviously kept stuck somewhere in my mind and that is why I wrote the initial comment.
The repeating digits in the fractional part are essentially the output of a linear congruence random number generator. Finding the period is almost the same as determining the order of an element of the state space, but by chance you might get a shorter repetition that fits into the longer one.
Finding the order of an element in a group is a hard problem on par with factoring http://mathworld.wolfram.com/GroupOrder.html
I'm not sure about cases where the state space doesn't form a multiplicative group, but I doubt that it gets much easier.
Finding the order of an element in a group is a hard problem on par with factoring http://mathworld.wolfram.com/GroupOrder.html
I'm not sure about cases where the state space doesn't form a multiplicative group, but I doubt that it gets much easier.
MathWorld is often a decent resource for finding this kind of information: <http://mathworld.wolfram.com/DecimalPeriod.html>, but it's often difficult to arrive at those pages if you don't have any idea where to start.
I certainly also looked through MathWorld but I failed to make some - in hindsight - rather obvious connections or did not pay enough attention. I certainly knew what the discrete logarithm problem is and that it is hard but I did not realize until a few minutes ago that finding the multiplicative order is just that only with a different name. And now being aware of that the MathWorld article spells exactly this out in a way that it seems impossible to miss.
There's an extra > at the end of your link.
Perhaps a side digression to this conversation, but... in spite of all the cool languages we talk about that have many on-paper advantages of either safety, convenience, or support of paradigms, is there really any option besides C for this kind of thing?
It still seems that, if you want to write Core Infrastructual Code that can be run anywhere, on anything, that links with any other software without any caveats, disclaimers, or compromise, C is still king.
Could this have practically been written in rust? haskell? ...and still have the same degree of embeddability and reach?
It still seems that, if you want to write Core Infrastructual Code that can be run anywhere, on anything, that links with any other software without any caveats, disclaimers, or compromise, C is still king.
Could this have practically been written in rust? haskell? ...and still have the same degree of embeddability and reach?
Fortran, assembler
In my limited experience, no. All of the languages you mention put developer experience above portability, so usually require a basic runtime. C on the other hand is not interested in saving you from the machine, which allows for extremely minimal programs. Even C++ isn't great here compared to C. Then again, debugging segfaults isn't for everyone, so most are fine giving up C.
> All of the languages you mention put developer experience above portability, so usually require a basic runtime.
Rust's runtime is comparable to C's. It is less portable because it uses LLVM as a backend, not because of the runtime.
> Then again, debugging segfaults isn't for everyone, so most are fine giving up C.
Depending on the application, having vulnerabilities due to undefined behavior is a much bigger problem than having to debug segfaults.
Rust's runtime is comparable to C's. It is less portable because it uses LLVM as a backend, not because of the runtime.
> Then again, debugging segfaults isn't for everyone, so most are fine giving up C.
Depending on the application, having vulnerabilities due to undefined behavior is a much bigger problem than having to debug segfaults.
+1 for mentioning Rust. Which in contrast to C++ has all the goodies without the runtime dependencies.
I wouldn't say LLVM is such a great restriction, except when talking about embedded (too few architectures supported here). The major desktop and server architectures are supported (AFAIK), even WASM.
I wouldn't say LLVM is such a great restriction, except when talking about embedded (too few architectures supported here). The major desktop and server architectures are supported (AFAIK), even WASM.
then port it to gcc?
We’d love to see a gcc front end for Rust. Someone attempted one in the pre-1.0 days, which was... not simple. Someone has to step up and do it though.
You could use D in -betterC mode. It allows compiling code that doesn't rely on runtime (has some limitations though, most notably lack of garbage collector, forcing you to use manual allocation).
The cool languages themselves tend to have this kind of feature written in C "under the hood".
GHC Haskell Runtime:
https://github.com/ghc/ghc/tree/master/rts
GHC Haskell Runtime:
https://github.com/ghc/ghc/tree/master/rts
For Mac users trying this at home, I made the following changes to build and run the tinypi tests on Sierra:
1. Comment out the malloc.h include in tinypi.c.
2. Add -Wno-unused-function to CFLAGS.
3. Use shasum, rather than sha1sum, in the test target.
1. Comment out the malloc.h include in tinypi.c.
2. Add -Wno-unused-function to CFLAGS.
3. Use shasum, rather than sha1sum, in the test target.
Having malloc.h as an include is bizarre, anyway. The standard location for malloc and related functions is stdlib.h, and removing that include should let it still compile on any reasonable platform.
> The basic arithmetic operations (addition, subtraction, multiplication, division, square root) have a near linear running time
I thought the best known multiplication algorithms (the ones based on Fourier transforms) for large numbers are a little slower than O(N log(N)). Has something better been found?
I thought the best known multiplication algorithms (the ones based on Fourier transforms) for large numbers are a little slower than O(N log(N)). Has something better been found?
"near linear" is sometimes used to refer to O(n log(n)^O(1)), that is, things that are linear other than a polynomial of log(n). This includes Schönhage–Strassen multiplication, with its O(n log(n) log(log(n))) complexity.
There is also BSDNT[0] (which is obviously BSD licensed) in case anyone is still interested... I think this takes the cake though (as it is MIT licensed and superior).
[0]: https://github.com/wbhart/bsdnt
[0]: https://github.com/wbhart/bsdnt
I also just realized that this is only for big numbers, not arbitrary precision.
MIT license makes this great for WebAssembly!
Looking forward to benchmarking this in a wasm context.
Looking forward to benchmarking this in a wasm context.
You can already try LibBF online at http://numcalc.com
Looking at the code this should be quite simple since libc is the only dep. You can Emscripten at first to be easy, but to wrap it up as a JS lib, you can probably implement the libc calls yourself easily and remove the Emscripten dep. The author even broke off bf_realloc to put your own alloc impl (which can be written in WASM, but I suggest writing it in C before the compilation).
cheers, sadly I won't be exploring JS lib wrapper approach.
we're already deep down the EMSCRIPTEN_BINDINGS() rabbit-hole :)
originally, the fixed size heap (/w optimisations) meant unbound use from a JS lib was too hard to reason about; hence avoided entirely in favour of a emscripten blackbox.
we're already deep down the EMSCRIPTEN_BINDINGS() rabbit-hole :)
originally, the fixed size heap (/w optimisations) meant unbound use from a JS lib was too hard to reason about; hence avoided entirely in favour of a emscripten blackbox.
Do you have a link for this ?
This needs some major work to be packagable into a library, esp. for non amd64 hosts. There's also no src repo, as with all his previous hacks.
Who will do the autotooling and repo hosting? I've added some sugar in my fork on github https://github.com/rurban/libbf/tree/my but I'm not sure yet if I will use it.
But the low constant overhead for small numbers, the small size, no assembler tricks and the MIT license makes it very attractive.
Who will do the autotooling and repo hosting? I've added some sugar in my fork on github https://github.com/rurban/libbf/tree/my but I'm not sure yet if I will use it.
But the low constant overhead for small numbers, the small size, no assembler tricks and the MIT license makes it very attractive.
I think the name is also already taken, for brainfuck and bloom filter.
This is perfect for jq[0].
[0] https://github.com/stedolan/jq
[0] https://github.com/stedolan/jq
...How is it perfect for jq?
It has a suitable license (MIT), is rather complete, and has semantics that will not surprise jq users.
Why does jq need arbitrary precision floating-point numbers?
Because JSON's spec defines a "number" as an arbitrary-precision floating point number, no limit on how many digits can be before or after the decimal or on what positive or negative integers you can put after the "e". Despite the name "JavaScript Object Notation," it does not inherit JavaScript's traditional interpretation that a "number" is an IEEE 754 double.
(JSON also does not permit infinities or NaN, which JavaScript does permit.)
http://json.org/
(JSON also does not permit infinities or NaN, which JavaScript does permit.)
http://json.org/
Wouldn't jq be better served by a library supporting arbitrary precision decimal values? LibBF is for base-2 floating point.
Probably yes, because you always want to convert to and from base 10 strings. mpdecimal [1] (which was developed for Python) might be a better choice, because it uses base 10 internally.
[1] http://www.bytereef.org/mpdecimal/
[1] http://www.bytereef.org/mpdecimal/
[deleted]
not that it's needed, but JSONschema has this kind of crazy part where it lets you do multipleOf for floating point numbers. Of course 2.7 is not a multiple of 0.3 in IEEE-64-bit FP.
The most FAQ / most filed issue for jq may well be that IEEE754 sucks. Users are frequently surprised by IEEE754 issues.
[deleted]
Does anyone know how LibBF stacks up against GMP? [5]
[1] http://ffmpeg.org/
[2] https://bellard.org/tcc/
[3] https://bellard.org/jslinux/
[4] https://bellard.org/
[5] https://en.wikipedia.org/wiki/GNU_Multiple_Precision_Arithme...