# Minimally promoting division

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

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

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.

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

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.