Minimally promoting division

General Usage
1

I could use a division function or operator that did minimal type promotion. E.g. a divide(x,y) function and associated operator that would return a Rational for rational or integer x,y and a floating-point number if either of x,y is floating-point. Neither of Bases’s / or // does the job. / returns a float for integer x,y, and // is not defined for floating-point types.

Any advice on the best approach? Would defining // to be / for floating-point types be reasonable? Wouldn’t it be good to have this in Base?

This is triggered for me by teaching polynomials in a numerical methods class. I’d like my Newton Divided differences algorithm for polynomial interpolation to return rational polynomial coefficients if the x,y datapoints are given as integers or rationals, and floating-point coefficients if either of x,y are floating point.

"""
function newtondivdiff(x,y)

Compute interpolating polynomial passing through points (x_i, y_i) 
Returns (c,b): the polynomial coeffs {c} and base points {b} for 
  interpolating polynomial p(x) = c[1] + c[2] (x-b[1]) + c[3] (x-b[1])(x-b[2]) + ...
""" 
function newtondivdiff(x,y)
    N = length(x)    
    length(y) == N || throw(ArgumentError("length(x)=$N ≠ length(y)=$length(y)"))

    # determine type of y/x. Would like this to be Rational for rationals and integers, floating-point otherwise
    T = typeof(one(eltype(y))/one(eltype(x))) 
    F = zeros(T,N,N) # matrix to store divided differences (wastes half of F allocation)
    F[:,1] = y
    
    for j=2:N
        for i=1:N+1-j
            F[i,j] = (F[i+1,j-1]-F[i,j-1])/(x[i+j-1]-x[i])  # would division with minimal promotion
        end
    end
    (F[1,:], x[1:N-1])
end

# give integer inputs
xdata = [-2 0 1 3]
ydata = [10 -11 1 0]
c,b = newtondivdiff(xdata, ydata)

@show c
@show b;

c = [10.0, -10.5, 7.5, -2.3333333333333335]  # want Rational output
b = [-2, 0, 1]
1 Like
2

How about defining your own function divide that has the wanted behavior?

3

That’s my plan for now. But it seems to me like a feature missing from the language, and I’m curious if others think so and what the best fix would be.

4

Note that you can also define an infix operator by appending a subscript, e.g. /ᵣ.

This doesn’t seem crazy to me. It would ideally have to be done in Base to avoid type piracy.

You could alternatively define:

newtondivdiff(x::AbstractVector{<:Integer}, y::AbstractVector{<:Integer}) =
    newtondivdiff(Rational.(x), Rational.(y))

and then just use / (which is exact for rationals or mixtures of rationals and integers).

1 Like
5

That does it. Thanks! I implemented it as

newtondivdiff(x::AbstractVector{<:Integer}, y::AbstractVector{<:Integer}) =
    newtondivdiff(x, Rational.(y))

so that the returned base points x[1:N-1] keep their original type.

I’ll experiment with defining // for floats as / and maybe make a PR.