How can I set the level of precision for Raku's sqrt?

Question

With Perl, one could use bignum to set the level of precision for all operators. As in:

use bignum ( p => -50 );

print sqrt(20);    # 4.47213595499957939281834733746255247088123671922305

With Raku I have no problems with rationals since I can use Rat / FatRat, but I don't know how to use a longer level of precision for sqrt

say 20.sqrt # 4.47213595499958

Here is a non-answer: 1) You can change the _printed_ precision with `sprintf`. But this does not increase the actual precision – `say sprintf('%.50f', 20.sqrt)` prints `4.47213595499958000000000000000000000000000000000000`. I am not aware of a way to change the _actual_ precision of the `Num` type in Raku, though I'd certainly be interested in learning otherwise. Thanks for the interesting question. — codesections, Aug 19 '20 at 10:19
How does Perl's `bignum` do that? If it's just altering the *printed* precision, see @codesections' comment above. If it's altering the *computed* precision, what numerical techniques/library is it using? — raiph, Aug 19 '20 at 12:20
@raiph, I don't know about the inners workings in perl, but I believe that is not just a simple *printed* precission, as I can set a precissión of `10,000` and the extra digits are not 'zeroed'. Not to mention that it takes a lot to compute (7 seconds for a sqrt(20) with 10,000 digits of precision) — Julio, Aug 19 '20 at 12:55
OK, so I checked perl's internals. `Math::BigInt::Calc;` provides a `_sqrt` method for bigints that reads `square-root of $x in place. Compute a guess of the result (by rule of thumb), then improve it via Newton's method.`. That `_sqrt` method is called on the package `Math::BigFloat` by the `bsqrt` method. That `bsqrt` method does (amon other things): `sqrt(2) = 1.4 because sqrt(2*100) = 1.4*10; so we can increase the accuracy of the result by multiplying the input by 100 and then divide the integer result of sqrt(input) by 10. Rounding afterwards returns the real result.` — Julio, Aug 19 '20 at 16:09
Links to sources (see my previous comment): [Math::BigInt::Calc](https://metacpan.org/release/Math-BigInt/source/lib/Math/BigInt/Calc.pm) and [Math::BigFloat](https://metacpan.org/release/Math-BigInt/source/lib/Math/BigFloat.pm) — Julio, Aug 19 '20 at 16:12

Julio · Accepted Answer · 2020-09-09T16:28:37.070

As stated on Elizabeth's answer, sqrt returns a Num type, thus it has limited precision. See Elizabeth's answer for more detail.

For that reason I created a raku class: BigRoot, which uses newton's method and FatRat types to calculate the roots. You may use it like this:

use BigRoot;

# Can change precision level (Default precision is 30)
BigRoot.precision = 50;

my $root2 = BigRoot.newton's-sqrt: 2;
# 1.41421356237309504880168872420969807856967187537695

say $root2.WHAT;
# (FatRat)

# Can use other root numbers
say BigRoot.newton's-root: root => 3, number => 30;
# 3.10723250595385886687766242752238636285490682906742

# Numbers can be Int, Rational and Num:
say BigRoot.newton's-sqrt: 2.123;
# 1.45705181788431944566113502812562734420538186940001

# Can use other rational roots
say BigRoot.newton's-root: root => FatRat.new(2, 3), number => 30;
# 164.31676725154983403709093484024064018582340849939498

# Results are rounded:

BigRoot.precision = 8;
say BigRoot.newton's-sqrt: 2;
# 1.41421356

BigRoot.precision = 7;
say BigRoot.newton's-sqrt: 2;
# 1.4142136

In general it seems to be pretty fast (at least compared to Perl's bigfloat)

Benchmarks:

|---------------------------------------|-------------|------------|
|  sqrt with 10_000 precision digits    |    Raku     |    Perl    |
|---------------------------------------|-------------|------------|
|  20000000000                          |    0.714    |   3.713    |
|---------------------------------------|-------------|------------|
|  200000.1234                          |    1.078    |   4.269    |
|---------------------------------------|-------------|------------|
|  π                                    |    0.879    |   3.677    |
|---------------------------------------|-------------|------------|
|  123.9/12.29                          |    0.871    |   9.667    |
|---------------------------------------|-------------|------------|
|  999999999999999999999999999999999    |    1.208    |   3.937    |
|---------------------------------------|-------------|------------|
|  302187301.3727 / 123.30219380928137  |    1.528    |   7.587    |
|---------------------------------------|-------------|------------|
|  2 + 999999999999 ** 10               |    2.193    |   3.616    |
|---------------------------------------|-------------|------------|
|  91200937373737999999997301.3727 / π  |    1.076    |   7.419    |
|---------------------------------------|-------------|------------|

If want to implement your own sqrt using newton's method, this the basic idea behind it:

sub newtons-sqrt(:$number, :$precision) returns FatRat {
    my FatRat $error = FatRat.new: 1, 10 ** ($precision + 1);
    my FatRat $guess = (sqrt $number).FatRat;
    my FatRat $input = $number.FatRat;
    my FatRat $diff = $input;

    while $diff > $error {
        my FatRat $new-guess = $guess - (($guess ** 2 - $input) / (2 * $guess));
        $diff = abs($new-guess - $guess);
        $guess = $new-guess;
    }

    return $guess.round: FatRat.new: 1, 10 ** $precision;
}

I'm pretty sure this will only work for small integers. In my tests it works great for 10,000 digits of the square root of `2`, is OK for `20`, but is awfully slow for `200`. — raiph, Aug 23 '20 at 17:18
But it's still nice to see this. :) I've created tinyurl.com/raku-newtons-method as a permanent URL linking to a variant of your code. I didn't change Newton's method (!) but just altered the last couple lines of code to compare the number calculated with a set of digits slurped from Input. I chose 2 as the number to square root because there's a convenient public source for the correct first 10,000 digits. — raiph, Aug 23 '20 at 17:21
Hehe, I did same test on local with first 100_000 sqrt(2) digits. I guess newthon's method works!! :-) You are right it is not really fast with big numbers. But at least it is comparable (if not faster) to Perl's version. — Julio, Aug 23 '20 at 17:57
Hi @raiph, I made some improvements. It seems that the slow bigger ints issue was because of my poor choice for the first guess number (which was `number/2`) I knew it was not very accurate, but I expected the algorithm to cath the real number 'soon'. But that was not the case for bigger numbers, since the bigger the number, the bigger is the difference from `number/2` to the final sqrt. So now I'm using raku's original sqrt as the first guess and now there is no difference from small numbers to big ones. I also changed the code to being able to accept a Rational number aswell. — Julio, Aug 23 '20 at 22:39
\o/ Smokin'! How large is "large numbers"? What about a 20,000 digit integer or a `FatRat` with a random 1,000 digit numerator and denominator? I guess the main thing is that the iterative aspect does no division, and instead just addition and subtraction of integers, plus one integer square. — raiph, Aug 23 '20 at 23:07
@raiph, I added some benchmarks and changed the code to allow for input `Num`, `Rational` and `Int` — Julio, Aug 24 '20 at 00:41
Would consider putting this logic into a module? Or allow someone else to put it in a module? — Elizabeth Mattijsen, Aug 24 '20 at 07:54
@ElizabethMattijsen Sure! I was planning to add a module for this. Perhaps make it more generic to allow for any type of root, not just square roots. Perhaps even higher precision constants such `pi`, `e`, ... would be welcome. Do you think we have any chance a module like this would make into core? If so, what changes what we need to do? — Julio, Aug 27 '20 at 11:43
@Julio: so far, only one module made it from the ecosystem to the core: `NativeCall`. But having it in the ecosystem would already be pretty good! And inclusion in a package such as Rakudo Star would be much more likely. — Elizabeth Mattijsen, Aug 27 '20 at 18:09
I finally created the module, thank you for your suggestions and hints — Julio, Sep 09 '20 at 16:30

score 6 · Answer 2 · answered Aug 19 '20 at 10:51

6

In Rakudo, sqrt is implemented using the sqrt_n NQP opcode. Which indicates it only supports native nums (because of the _n suffix). Which implies limited precision.

Internally, I'm pretty sure this just maps to the sqrt functionality of one of the underlying math libraries that MoarVM uses.

I guess what we need is an ecosystem module that would export a sqrt function based on Rational arithmetic. That would give you the option to use higher precision sqrt implementations at the expense of performance. Which then in turn, might turn out to be interesting enough to integrate in core.

answered Aug 19 '20 at 10:51

Elizabeth Mattijsen

25,654
3
75
105

2

I'm pretty sure a rational `sqrt` would almost never help, and would usually make things a lot worse. ❶ Per math, a rational `sqrt` *can't* help with *any* irrational numbers (which would of course be stored as `Num`s) and *can't* return an accurate result for *almost all* rational numbers, even using `FatRat`s. ❷ Most `Num`s where float precision matters are already inaccurate, so converting them to `Rational` equivalents wouldn't help. ❸ For the sorts of numbers where precision of floats (`Num`s) matters, using a rational `sqrt` would almost always be *pathologically slow*. – raiph Aug 19 '20 at 12:29
@raiph I tried a rational sqrt (using newthon's method) and It is pretty fast (at least compared to Perl's bigfloat) – Julio Aug 23 '20 at 13:17
Heh. I will leave my original comment up as yet another reminder to myself and others of just how wrong I can be. (I was right about the math but thoroughly wrong about the computer science / reality.) – raiph Aug 24 '20 at 00:12

How can I set the level of precision for Raku's sqrt?

2 Answers2