| |
 | Posted by H. S. Teoh
in reply to Geroge Little |  Permalink Reply |
|
H. S. Teoh 
Posted in reply to Geroge Little
| On Mon, Apr 10, 2017 at 06:54:54PM +0000, Geroge Little via Digitalmars-d-learn wrote:
> Is there support for BigFloat in phobos or any other package? I was playing around with D and wrote some code that calculates a Fibonacci sequence (iterative) with overflow detection that upgrades ulong to BigInt. I also wanted to use Binet's formula which requires sqrt(5) but it only works up to n=96 or so due to floating point precision loss. The code is here:
>
> https://gist.github.com/ggl/38458b57b1eb3945ce447c8bf1d4e458
There is no BigFloat in phobos, you could try looking at code.dlang.org to see if there's anything that you could use.
One way to use sqrt(5) in your calculations might be to use a quadratic
rational representation, i.e., as a triple (x,y,z) representing (x +
y*√5)/z where x, y, z are integers (can be BigInt).
The main observation is that numbers of this form are closed under +, -, *, and / (excluding divison by 0 of course), and thus form a field. Even nicer, is that all calculations are exact (it can be reduced to purely integer manipulations). So you could ostensibly implement a QuadRational type based on BigInt that you can use to freely compute with √5.
The key to implementing division is to note that (x + y*√5)*(x - y*√5) =
x^2 - 5*y^2, so you can use this fact to eliminate √5 from the
denominator so the rest of the computation can be carried out purely
using +, -, and *. I.e.:
(x1 + y1*√5)
------------
(x2 + y2*√5)
(x1 + y1*√5) (x2 - y2*√5)
= ------------ * ------------
(x2 + y2*√5) (x2 - y2*√5)
(x1 + y1*√5) * (x2 - y2*√5)
= ---------------------------
x2^2 - 5*y2^2
so you just multiply out the top, which will be of the form p+q*√5, and with the denominator you again have the representation (x + y*√5)/z. You can then reduce the representation by dividing each of x, y, z with gcd(x,y,z).
Note that this scheme works with both BigInt and built-in types like long/ulong, but BigInt is recommended because the square terms in the denominator means that performing divisions tend to overflow long/ulong pretty quickly.
Some time ago I wanted to implement a QuadRational library that
basically does the above, plus a neat algorithm for exact comparison,
but didn't finish it because I got a bit too ambitious and wanted to
support numbers of the form x + y*√a + z*√b + ... as well. Turns out
that was too much because the required representation would be
exponential in length in the number of different radicals you wish to
support. Perhaps a less ambitious library would be one that supports
numbers of the form (x + y*√r)/z, where r is fixed. It would then
encompass the complex numbers (by setting r=-1, though it loses
orderability in that case), combinations of √5 like you have. This can
be pretty useful for implementing exact arithmetic in certain
geometrical computations (mainly r=√2 for things containing octagons and
r=√5 for things involving pentagons, r=√3 for certain triangular
constructions).
T
--
It is not the employer who pays the wages. Employers only handle the money. It is the customer who pays the wages. -- Henry Ford
|
