Support rational-coefficient series for square root?

[Vectornaut]

When t is a monomial with rational coefficients, sqrt(1+t) ends up with floating point coefficients—even though they could be rational. As far as I can tell, the only obstacle to returning rational coefficients is the sqrt that gets applied to the leading term in the definition of sqrt!. The sqrt of a rational number is always a floating-point number, even if the square root happens to be rational.

When sqrt is applied to a series with rational coefficients, would it make sense to return a series with rational coefficients when the leading term is square? I think this could be done simply and efficiently using isqrt. (If desired, there could be an option that forces a floating-point result.)

I'd be happy to implement this and make a pull request. If you'd like me to do that, please let me know how to test the new code, and point out any conventions I should follow or pitfalls I should watch out for.

[dpsanders]

I agree that it would be nice to be able to do this; see https://github.com/JuliaDiff/TaylorSeries.jl/issues/51 for a similar request / idea and https://github.com/JuliaDiff/TaylorSeries.jl/pull/59 for a similar implementation.

However, there is the issue of type stability.

[lbenet]

This indeed has been proposed before, as David says.

Whenever I have encountered those cases, I proceed as follows:

julia> using  TaylorSeries

julia> t = Taylor1(Rational{Int}, 8)
 1//1 t + 𝒪(t⁹)

julia> p = sqrt(1+t)
 1.0 + 0.5 t - 0.125 t² + 0.0625 t³ - 0.0390625 t⁴ + 0.02734375 t⁵ - 0.0205078125 t⁶ + 0.01611328125 t⁷ - 0.013092041015625 t⁸ + 𝒪(t⁹)

julia> convert(Taylor1{Rational{Int}}, p)
 1//1 + 1//2 t - 1//8 t² + 1//16 t³ - 5//128 t⁴ + 7//256 t⁵ - 21//1024 t⁶ + 33//2048 t⁷ - 429//32768 t⁸ + 𝒪(t⁹)

julia> convert(Taylor1{Rational{Int}}, exp(t))
 1//1 + 1//1 t + 1//2 t² + 1//6 t³ + 1//24 t⁴ + 1//120 t⁵ + 1//720 t⁶ + 1//5040 t⁷ + 1//40320 t⁸ + 𝒪(t⁹)

This is maybe not the most elegant solution, but good enough.

In my opinion, to build an API we should distinguish among specific (user) cases, such as in sqrt(1//1+t) or exp(t), where the user knows the answer (or that the answer can be written nicely in some specific numeric types), from the more general use cases. As @Vectornaut points out, the conversion occurs is due to the fact that for sqrt(1) Julia returns a Float64, and that we use some recurrence relation whose initial condition are sqrt(1) (or exp(0)). As David also points out, this is related to type stability, which would then become an issue here.

Let me point out that there are some subtleties related to the conversion to Ratioinals:

julia> convert(Rational{Int}, 1/factorial(8))
28594283348384//1152921504606842881

julia> rationalize(1/factorial(8))
1//40320

We internally use rationalize (inside convert); maybe it would be better to overload rationalize, so there may be other issues aside from type stability.

[Vectornaut]

@dpsanders, thank you for pointing out the related issues! I can see better now how the type stability issues that inform the behavior of functions like sqrt and exp for numbers also apply to Taylor series.

@lbenet's comment about transcendental functions was enlightening for me, and it suggests to me that we should think about sqrt and exp differently. The graph of exp has only one rational point, so the rational-coefficient expansion there can be described as a constant. The graph of sqrt has rational points throughout it, so you really would need a function to access them all.

If coherence with functions of numbers is an important goal, would it make sense to define qsqrt for rational numbers by applying isqrt to the numerator and denominator, and then extend to rational series by analogy with sqrt? With this definition, I think qsqrt(f) would end up being sqrt(f) times the number qsqrt(f[1]) / sqrt(f[1]). Of course, this package isn't the right place to define qsqrt.

If the AlgebraicNumbers package gets published, there will be an elegant solution on the user end: just use series with AlgebraicNumber coefficients. It might not be as efficient as using rational coefficients, but it should work fine for casual calculations.

[lbenet]

Sorry for the delay to get back to this discussion...

Have you tested if TaylorSeries works fine with AlgebraicNumbers?

[lbenet]

TLDR: Sorry again for the long delay... I didn't know how to understand your question of qsqrt or isqrt, and finally it occurred to me how one could deal with it, instead of using convert. Performance is more or less the same as using convert, and there may be some issues with the rationals that i have not checked. Though the trick is nice, it doesn't always work.

Based on the implementation of sqrt here, you can define the following qsqrt:

#### I impose `a::Taylor1{Rational{T}}`, but this is not important...
function qsqrt(a::Taylor1{Rational{T}}) where {T<:Integer}
    ##### This first part deals with the possibility of factoring square terms; not essential but nice (but ugly code)
    # First non-zero coefficient
    l0nz = findfirst(a)
    aux = zero(rationalize(sqrt( constant_term(a) )))  ##### note the `rationalize` here (type stability)
    if l0nz < 0
        return Taylor1(aux, a.order)
    elseif l0nz%2 == 1 # l0nz must be pair
        throw(DomainError(a,
        """First non-vanishing Taylor1 coefficient must correspond
        to an **even power** in order to expand `sqrt` around 0."""))
    end

    # The last l0nz coefficients are set to zero.
    lnull = l0nz >> 1 # integer division by 2
    c_order = l0nz == 0 ? a.order : a.order >> 1
    c = Taylor1( aux, c_order )            #### This is the eventual output Taylor1{Rational{T}} series
    @inbounds c[lnull] = rationalize(sqrt( a[l0nz] ))   #### This "rationalizes" the first term of `c`
    #### This is standard for `sqrt`
    for k = lnull+1:c_order
        sqrt!(c, a, k, lnull)
    end

    return c
end

Two neat examples:

julia> using TaylorSeries

julia> tr = Taylor1(Rational{Int}, 15)
 1//1 t + 𝒪(t¹⁷)

julia> qsqrt(1//1 + tr)
 1//1 + 1//2 t - 1//8 t² + 1//16 t³ - 5//128 t⁴ + 7//256 t⁵ - 21//1024 t⁶ + 33//2048 t⁷ - 429//32768 t⁸ + 715//65536 t⁹ - 2431//262144 t¹⁰ + 4199//524288 t¹¹ - 29393//4194304 t¹² + 52003//8388608 t¹³ - 185725//33554432 t¹⁴ + 334305//67108864 t¹⁵ + 𝒪(t¹⁶)

julia> sqrt(1//1 + tr)   # This is the usual result
 1.0 + 0.5 t - 0.125 t² + 0.0625 t³ - 0.0390625 t⁴ + 0.02734375 t⁵ - 0.0205078125 t⁶ + 0.01611328125 t⁷ - 0.013092041015625 t⁸ + 0.0109100341796875 t⁹ - 0.009273529052734375 t¹⁰ + 0.008008956909179688 t¹¹ - 0.0070078372955322266 t¹² + 0.006199240684509277 t¹³ - 0.005535036325454712 t¹⁴ + 0.004981532692909241 t¹⁵ + 𝒪(t¹⁶)

julia> convert(Taylor1{Rational{Int}},sqrt(1//1+tr)) == TaylorSeries.qsqrt(1//1+tr)
true

### Another nice example
julia> qsqrt(1//4 + tr)
 1//2 + 1//1 t - 1//1 t² + 2//1 t³ - 5//1 t⁴ + 14//1 t⁵ - 42//1 t⁶ + 132//1 t⁷ - 429//1 t⁸ + 1430//1 t⁹ - 4862//1 t¹⁰ + 16796//1 t¹¹ - 58786//1 t¹² + 208012//1 t¹³ - 742900//1 t¹⁴ + 2674440//1 t¹⁵ + 𝒪(t¹⁶)

The whole idea is to make sure that rational arithmetic is used always, i.e., including the 0-th order term of c; then, the existent machinery works, though not always. For instance:

julia> qsqrt(1//2 + tr)
ERROR: OverflowError: -2982076586042449 * 54608393 overflowed for type Int64
...

As I mentioned above, there may also be issues related with the rational arithmetic, which I haven't investigated. For instance, the check with the converted series above doesn't work if i consider larger orders; not yet sure why.

Clearly, the same kind of idea can be used for other functions, e.g., the exp. Once said all this, I don't think the package should support these q-functions. When they work it's nice, but in general I think they will not.

[lbenet]

Closing. Feel free to reopen this issue.