Implement sqrt function

[Pzixel]

Related to rust-num/num-bigint#55

Should it be implemented on top of BigInts or some better approach is possible?

Ratio will be trickier to deal with precision, because you probably don't want to just round √(a/b) by operating on the numerator and denominator separately. That is, ⌊√a⌋/⌊√b⌋ would be quite lossy.

[cuviper]

I'm not really sure what's the right answer here. Approximating irrational numbers as a rational requires some level of control over how much accuracy is desired -- we can't be perfect, but we can get arbitrarily close, especially with BigInt that can grow as much as you like.

I have some code I was playing with for #32, to generate the most accurate approximations for a few interesting constants, given any particular Ratio<T> limit. For √2 it looks like:

# √2 ≈ 1.4142135624, f64 error 9.667e-17
================================
    type |      error | rational
---------|------------|---------
      i8 |   7.215e-5 | 99/70
      u8 |  -1.238e-5 | 239/169
     i16 |  -1.840e-9 | 27720/19601
     u16 | -3.158e-10 | 47321/33461
     i32 | -2.055e-19 | 1855077841/1311738121
     u32 | -2.055e-19 | 1855077841/1311738121
     i64 |  1.493e-38 | 6882627592338442563/4866752642924153522
     u64 | -2.561e-39 | 16616132878186749607/11749380235262596085
    i128 | -7.878e-77 | 133984184101103275326877813426364627544/94741125149636933417873079920900017937
    u128 | -1.352e-77 | 228725309250740208744750893347264645481/161733217200188571081311986634082331709

These values and their errors were computed based on 100 digits of the constant from Wolfram Alpha. It's possible to derive that basis mathematically, but I didn't bother because I was doing this for a variety of constants.

Note that even the 32-bit ratios are more accurate than f64 can represent! Which makes some sense, since a pair of 32-bit values has 64 significant bits together, or 62 w/ signs, vs. a 53-bit mantissa in f64.

And yet, it's probably not so useful to have ratios that are so near their type limits, since a lot of arithmetic requires multiplication that would overflow. That's not a problem for BigRational, of course.

[Pzixel]

I actually tried to implement superpi using this crate. However, it wasn't a big success: https://github.com/Pzixel/superpi-rs

It gives 3.1410... for 10 iterations (should give 10^2 correct digits). It's not directly related but it provides some experience.

P.S. It also shows how much clones you need to implement some algorithms.

https://github.com/Pzixel/superpi-rs/blob/2b8874c0bff2f480cdf299bb8e17d58218f2dc18/src/main.rs#L24-L29

It also seems having almost O(exp(n)) complexity:

[cuviper]

let new_b = BigRat::new(new_b_square.numer().sqrt(),new_b_square.denom().sqrt());

That's exactly the ⌊√a⌋/⌊√b⌋ which I said would be lossy.

I think even your initial ratio b for 1/√2 might be too inaccurate. I don't have experience with this algorithm, but it doesn't appear to be good at self-correcting for errors. If I set the initial b using the u32 ratio I gave above, your code then converges on 3.1415591..., or with my u64 ratio it gives 3.141592653787... (vs the true 3.141592653589...)

Maybe it would fare better with a Ratio::sqrt that doesn't round as much, but I'm not sure.

[cuviper]

P.S. It also shows how much clones you need to implement some algorithms.

You can do most operations by reference, e.g. let new_t = t - &p * &a_diff * a_diff;

[Pzixel]

That's exactly the ⌊√a⌋/⌊√b⌋ which I said would be lossy.

Yep, I know, but I don't know any better approach. Error is too high to make this algorithm work fine, because it should work with any precision, while we have only 11 right digits on u128. It's even worse than f64 which often has 16 significant digits.

You can do most operations by reference, e.g. let new_t = t - &p * &a_diff * a_diff;

Thank you for pointing that out, gonna try that.

[cuviper]

We could just use that sqrt-by-parts as an initial guess for the Babylonian method, but the question remains of how many iterations to compute. Maybe even just one, because the ratio part will quickly get large!

[Pzixel]

@cuviper maybe sqrt just should take some kind of precision as parameter?..

[ExpHP]

I don't get it, what is the use case of having a square root function on a rational that does approximations?

[cuviper]

Any sqrt that's not represented symbolically will be an approximation, including those for floating point. A rational sqrt has a chance to improve accuracy. (Whether that's needed is a separate question.)

[ExpHP]

But as already noted in the discussion, there are many problems with this:

• The problem is ill specified; rationals are not closed under the square root operation
  • Nor is there the existence of a "nearest" rational, under any possible rounding scheme (at least integers have that!)
• If you get a big denominator on a fixed-size integer type, then you can't perform virtually any math on it at all without getting an integer overflow. So what on earth are you going to do with it?

The fact that Rational32 can represent sqrt(2) more accurately than f64 doesn't mean much if any further operation you perform on it causes the world to explode.

This seems like a very niche problem that would be better serviced by a niche solution. Something that e.g. gives approximate rational roots to any rational polynomial.

(or at the very least, please don't call it sqrt!)

[ExpHP]

N.B. why I care so much:

Just earlier today, working on a diophantine equation solver, I almost wrote the following:

trait SqrtExt {
    /// Compute the (exact) square root, or panic if there isn't one
    fn sqrt(self) -> Self;
}

impl SqrtExt for Rational64 {
    ...
}

but then stopped and decided to make it a free function instead when I considered that it could get overridden by an inherent method.

[Pzixel]

I have been talking with guys on Math.StackOverflow and they convinced me it doesn't make sense to perform any kind of sqrt for rationals.

That's highly inefficient, because you'll end up with very large integers a,b, that may be far larger than needed for the desired precision level; and it'll be imprecise, because of the ⌊√a⌋/⌊√b⌋ problem mentioned.

So we probably need to implement arbitrary-sized fixed-point arithmetic in new repo like num-fp and implement all such operations there. I with to help if I can.

[cuviper]

For fixpoint, consider existing crates like https://crates.io/crates/bigdecimal

[liborty]

Or you could implement continued fractions, whereby all square roots of rationals are known to have repeating continued fraction terms. So you only need to evaluate the continued fraction of the square root up to where you detect repetition and the rest of arbitrary precision you get for free.

[1011X]

Maybe the sqrt method could return an Option<Ratio> instead, where if the numerator and denominator are perfect squares it returns their square root and None otherwise?

[ExpHP]

• Continued fractions (including the function to obtain one from the square root of a rational) are probably best implemented in another crate.
• I am sympathetic to the &Self -> Option<Self> use case and could maybe see it being called sqrt_checked.
• Both of these ideas are nonetheless far away from the spirit of the original idea in this thread, and it seems odd to be discussing them all in one space.

[1011X]

Not sure what you mean by "far away from the spirit of the original idea" when the reason other ideas are being brought up is because the original post asked

Should it be implemented on top of BigInts or some better approach is possible?

[ckaran]

Quick hack, as yet untested. YMMV.

EDIT. @noahdahlman pointed out that there wasn't a license on the original code, now added.

//! Copyright 2019, Cem Karan
//! 
//! Licensed under the Apache License, Version 2.0 (the "License");
//! you may not use this file except in compliance with the License.
//! You may obtain a copy of the License at
//! 
//!     http://www.apache.org/licenses/LICENSE-2.0
//! 
//! Unless required by applicable law or agreed to in writing, software
//! distributed under the License is distributed on an "AS IS" BASIS,
//! WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
//! See the License for the specific language governing permissions and
//! limitations under the License.

/// # Calculates the approximate square root of the value
///
/// Calculates the approximate square root of `value`.  If the returned value is
/// `Ok(_)`, then it is guaranteed to be within `epsilon` of the actual 
/// answer.  If `epsilon <= 0.0`, then `Err` is returned (the reason for the
/// bound of `0.0` is because the approximation algorithm is unable to return an
/// exact answer).  If `value < 0.0`, then `Err` is returned (`BigRational` is
/// a real valued object; it cannot represent complex values).  In both `Err`
/// cases, the value will be a `String` explaining what the error actually is.
///
/// # Parameters
///
/// - `value` - The value whose approximate square root you wish to obtain.  If
///     this is less than `0.0`, then `Err` will be returned.
/// - `epsilon` - The maximum acceptable difference between the returned value
///     and the actual value.  The returned value is in the range 
///     `[actual - epsilon, actual + epsilon]`.
///
/// # Returns
///
/// If everything went as expected, then `Ok(_)` will be returned, containing
/// a value that is within `± epsilon` of the actual value.  If anything went 
/// wrong, then `Err(_)` will be returned, containing a `String` outlining what
/// the problem was.
pub fn approx_square_root(value: BigRational, epsilon: BigRational) -> Result<BigRational, String> {
    if value < BigRational::zero() {
        return Err(format!("approx_square_root() cannot calculate the square \
                            root of negative values.  value = {}", value).to_owned());
    } else if epsilon <= BigRational::zero() {
        return Err(format!("approx_square_root() cannot calculate the square \
                            root with a non-positive epsilon.  \
                            epsilon = {}", epsilon).to_owned());
    }

    // I'm going to use the Babylonian method to find the square root.  This is
    // described at 
    // https://en.wikipedia.org/wiki/Methods_of_computing_square_roots#Babylonian_method
    // To do so, I need to have an initial seed value that is the approximate
    // square root.  This estimate will be refined until it is within 
    // epsilon of the real value.

    // Calculates seed values for all values >= 1.0.  This is used below when
    // calculating the seed value.
    #[inline]
    fn calc_seed(value: &BigRational) -> BigRational {
        let bits = value.ceil().to_integer().bits();
        let half_bits = bits / 2;
        let approximate = BigInt::one() << half_bits;
        BigRational::from_integer(approximate)
    };

    let mut x = if value >= BigRational::one() {
        calc_seed(&value)
    } else {
        // Because the value is less than one, I can't use the trick above 
        // directly.  Instead, I'm going to find the reciprocal, and then do the
        // trick above, and then use the reciprocal of that as the seed.
        calc_seed(&(value.recip())).recip()
    };

    // We now have an initial seed.  Time to refine it until it is within
    // epsilon of the real value.  Inlined functions could probably be really
    // inlined, but this is slightly easier for me to read.

    #[inline]
    fn calc_next_x(value: BigRational, x: BigRational) -> BigRational {
        let two = BigRational::one() + BigRational::one();
        (x + (value / x)) / two
    };

    #[inline]
    fn calc_approx_error(value: BigRational, x: BigRational) -> BigRational {
        let two = BigRational::one() + BigRational::one();
        ((value - (x * x)) / (x * two)).abs()
    }

    while calc_approx_error(value, x) > epsilon {
        x = calc_next_x(value, x);
    }

    Ok(x)
}

[nbraud]

Maybe the sqrt method could return an Option<Ratio> instead, where if the numerator and denominator are perfect squares it returns their square root and None otherwise?

FWIW, I have a highly-optimized implementation of perfect square detection (for integers) that I wrote while working on uutils/coreutils' implementation of integer factorization: factor/src/numeric/sqrt.rs
Given an (exact) integer square root function, getting exact square roots of some x: Ratio<_> would be simply let y = x.reduced(); Ratio::new_raw(exact_sqrt(y.numer())?, exact_sqrt(y.denom())?).

It's not yet upstreamed (mostly because, for that integer factoring application, there were some I/O-related optimizations I'd need to get in first) but I'd be entirely happy to see that code used for your usecase; if there's interest, I could submit a PR to include a (generalized) version in num-integer.

[cmpute]

FWIW as well, I'm implementing a crate for irrational numbers. The QuadraticSurd struct can exactly meet the demand for this, as it can represent any irrational numbers in form (a+b√r)/c.
I already implemented a very simple to_rational method for it (basically just rounded down the square root part), and the square root of a rational number can be achieved by QuadraticSurd<T>::from_sqrt(t: Ratio<T>). The missing part might be a method to produce more accurate rational estimation of the surd.

What I have in mind now is to make all irrational types have a method that takes an additional parameter to determine the precision of rational estimation. It might looks like

trait Irrational {
    fn to_accurate_rational(&self, iterations: u32) -> BigRational;
}

With these methods implemented, one can get the rational approximation of the square root by QuadraticSurd::from_sqrt(num).to_accurate_rational(10)

@nbraud num-integer already has a sqrt implementation that has been optimized.

[noahdahlman]

Quick hack, as yet untested. YMMV.

/// # Calculates the approximate square root of the value
///
/// Calculates the approximate square root of `value`.  If the returned value is
/// `Ok(_)`, then it is guaranteed to be within `epsilon` of the actual 
/// answer.  If `epsilon <= 0.0`, then `Err` is returned (the reason for the
/// bound of `0.0` is because the approximation algorithm is unable to return an
/// exact answer).  If `value < 0.0`, then `Err` is returned (`BigRational` is
/// a real valued object; it cannot represent complex values).  In both `Err`
/// cases, the value will be a `String` explaining what the error actually is.
///
/// # Parameters
///
/// - `value` - The value whose approximate square root you wish to obtain.  If
///     this is less than `0.0`, then `Err` will be returned.
/// - `epsilon` - The maximum acceptable difference between the returned value
///     and the actual value.  The returned value is in the range 
///     `[actual - epsilon, actual + epsilon]`.
///
/// # Returns
///
/// If everything went as expected, then `Ok(_)` will be returned, containing
/// a value that is within `± epsilon` of the actual value.  If anything went 
/// wrong, then `Err(_)` will be returned, containing a `String` outlining what
/// the problem was.
pub fn approx_square_root(value: BigRational, epsilon: BigRational) -> Result<BigRational, String> {
    if value < BigRational::zero() {
        return Err(format!("approx_square_root() cannot calculate the square \
                            root of negative values.  value = {}", value).to_owned());
    } else if epsilon <= BigRational::zero() {
        return Err(format!("approx_square_root() cannot calculate the square \
                            root with a non-positive epsilon.  \
                            epsilon = {}", epsilon).to_owned());
    }

    // I'm going to use the Babylonian method to find the square root.  This is
    // described at 
    // https://en.wikipedia.org/wiki/Methods_of_computing_square_roots#Babylonian_method
    // To do so, I need to have an initial seed value that is the approximate
    // square root.  This will estimate will be refined until it is within 
    // epsilon of the real value.

    // Calculates seed values for all values >= 1.0.  This is used below when
    // calculating the seed value.
    #[inline]
    fn calc_seed(value: &BigRational) -> BigRational {
        let bits = value.ceil().to_integer().bits();
        let half_bits = bits / 2;
        let approximate = BigInt::one() << half_bits;
        BigRational::from_integer(approximate)
    };

    let mut x = if value >= BigRational::one() {
        calc_seed(&value)
    } else {
        // Because the value is less than one, I can't use the trick above 
        // directly.  Instead, I'm going to find the reciprocal, and then do the
        // trick above, and then use the reciprocal of that as the seed.
        calc_seed(&(value.recip())).recip()
    };

    // We now have an initial seed.  Time to refine it until it is within
    // epsilon of the real value.  Inlined functions could probably be really
    // inlined, but this is slightly easier for me to read.

    #[inline]
    fn calc_next_x(value: BigRational, x: BigRational) -> BigRational {
        let two = BigRational::one() + BigRational::one();
        (x + (value / x)) / two
    };

    #[inline]
    fn calc_approx_error(value: BigRational, x: BigRational) -> BigRational {
        let two = BigRational::one() + BigRational::one();
        ((value - (x * x)) / (x * two)).abs()
    }

    while calc_approx_error(value, x) > epsilon {
        x = calc_next_x(value, x);
    }

    Ok(x)
}

I wanted to use this code since there doesn't appear to be a good solution available, but I noticed you didn't include a license. Is it safe to assume this is MIT? I don't want to steal your work, could you specify a license? Thanks!

[ckaran]

@noahdahlman My apologies! It is Apache 2.0. Be warned! That was really a quick hack without any significant testing whatsoever! Before relying on it, please spend some time verifying that it actually works as expected!

[liborty]

I am not sure if this was mentioned but it is possible to compute and store square roots of unlimited precision, using continued fractions, which will have (guaranteed) recurring group of terms. We only need to store this recurring set. Such continued fraction can then be evaluated to any required precision, on demand.

[ckaran]

@liborty, you actually mentioned it earlier in this very same conversation! 😉

IMHO, @cmpute's crate is the solution that solves @cuviper and @ExpHP's concerns with regards to irrationality. Perhaps @cmpute and @cuviper could work together to combine forces in some way? Maybe bringing num-irrational into the num family of crates?

[cmpute]

@liborty, you actually mentioned it earlier in this very same conversation! 😉

IMHO, @cmpute's crate is the solution that solves @cuviper and @ExpHP's concerns with regards to irrationality. Perhaps @cmpute and @cuviper could work together to combine forces in some way? Maybe bringing num-irrational into the num family of crates?

I'm glad that you are interested in this. However, I'm actually planning to integrate the num-irrational trait with my dashu crates, because modeling the quadratic surd in a generic way is pretty troublesome.

[ckaran]

@cmpute Ah, I see. I actually thought it was going to be a part of the num family because of it's name; num-irrational really sounds like it belongs there.

[cuviper]

There are no formal namespaces, so it's free to be named num-* without being part of this @rust-num org.

[cmpute]

@cmpute Ah, I see. I actually thought it was going to be a part of the num family because of it's name; num-irrational really sounds like it belongs there.

Yes, it's indeed intended to interop with the num crates (hence the name), and it does support conversions pretty well now. I'm open to integrate it into the num ecosystem, but I don't plan to improve the num-irrational crate now.

[robs-zeynet]

Hi all - thanks for the great crate. What's the status of this discussion? FYI: my use-case is CAD software where I'm hoping to store coordinates in an arbitrary-sized format and I'm considering using this crate. But, like in many CAD systems, I'll often have to compute the distance between two points which requires a sqrt() function. Looking over the thread, @ckaran 's solution would work for me as it's easy for me to pick an epsilon that's smaller than physical machines can handle. Any thoughts about merging that solution into the code? If this whole thread has gone cold, I can put it up as a pull request...

Please let me know and thanks in advance.

[ckaran]

@robs-zeynet I, for one, completely forgot about this until you brought it back up!

I suspect that it won't get rolled into the num-rational codebase because it really violates the spirit of the num-rational crate (the operations are closed under the rationals). You might be better off just copying my code into your own code. The advantage there is that while the code I wrote does seem to work, anybody that looks at it for more than 5 seconds can immediately see that there are a lot of opportunities for performance improvements! In short, take it as a starting point and tune it until it works for you.