Is there a difference between Divide[a,b] and a/b?

Question

In this comment it was asserted that Divide[a,b] and a/b are different, though the documentation indicates that they are the same. In particular, it was asserted that a/b is evaluated as a * 1/b, whereas Divide[a,b] performs the division directly. It was further asserted that this could result in an observable difference in behavior for machine-precision numbers.

The question is: is there a difference? If so, what is it exactly? How does it manifest itself? Is there an example where Divide gives a different result than /?

I tried a few edge cases with $MaxMachineNumber and $MinMachineNumber, but didn't find any differences in behavior between the two.

Answer 1

Oleksandr is correct about the way evaluation works. a/b seems to be interpreted (parsed) directly as Times[a, Power[b,-1]], or more readably: $a\times b^{-1}$. Divide[a,b] is interpreted as is. Evaluation then proceeds from these forms, and the arithmetic is carried out differently for the two cases: either $a\times (1/b)$ or $a/b$.

Here are some examples that illustrate the evaluation sequence:

In[95]:= 
On[]
Divide[4,8]
Off[]

During evaluation of In[95]:= On::trace: On[] --> Null. >>
During evaluation of In[95]:= Divide::trace: 4/8 --> 1/2. >>

Out[96]= 1/2

In[98]:= 
On[]
4/8
Off[]

During evaluation of In[98]:= On::trace: On[] --> Null. >>
During evaluation of In[98]:= Power::trace: 1/8 --> 1/8. >>
During evaluation of In[98]:= Times::trace: 4/8 --> 4/8. >>
During evaluation of In[98]:= Times::trace: 4/8 --> 1/2. >>

Out[99]= 1/2

This can indeed theoretically lead to different machine precision results. Let's find out if it really does! We are going to compare the complete binary representation of the results, and we won't use == or === (which both have some tolerance).

Table[{k, RealDigits[k/137., 2] === RealDigits[Divide[k, 137.], 2]}, {k, 1, 20}]

(* ==> {{1, True}, {2, True}, {3, True}, {4, True}, {5, True}, {6, True}, {7, True}, 
        {8, True}, {9, True}, {10, True}, {11, True}, {12, True}, {13, True}, 
        {14, True}, {15, False}, {16, True}, {17, True}, {18, True}, {19, True}, {20, True}} 
 *)

So Divide[15, 137.] and 15/137. really do lead to different results. Conclusion: yes, there is an observable difference.

Again, keep in mind that even === has some tolerance (Internal`$SameQTolerance) when comparing machine precision numbers (though less than ==, Internal`$EqualTolerance) and 15/137. === Divide[15, 137.] returns True. The difference is there though as evidenced by the full 53-bit binary representation.

So will you ever see the effects of this in practice? Theoretically the error may accumulate, and there are some functions which do not honour these tolerances and perform strict comparisons (try e.g. Union[{15/137., Divide[15, 137.]}], which returns a list of length 2).

Answer 2

I have just a little to add to @Szabolcs' answer, and this seemed an appropriate place rather than a separate Q&A. Since there is a theorem involved (see below), it should be pointed out that this theorem tends to fail to hold in Mathematica. Well, wait, it's a theorem: it must true! Then let us say that it tends to be difficult to apply this theorem to Mathematica code.

As pointed out already, a / b is short for Times[a, Power[b, -1]], which involves two roundings. It also happens that sometimes Divide[a, b] evaluates to Times[a, Power[b, -1]], for example, symbolically and in some cases presented below. In fact, I sometimes find it difficult to get Divide to Divide[] and not to multiply by the reciprocal.

There's a remarkable theorem for base-10 addicts due to W. Kahan (see David Goldberg, "What every computer scientist should know about floating-point arithmetic", ACM Computing Surveys, Vol 23, No 1, March 1991, 5–48, Theorem 7):

Suppose $x$ is an integer with $|x|<2^{p-1}$ in a binary floating-point format with precision $p$ bits and floating-point operations that are exactly rounded. Then if $n$ is an integer of the form $n = 2^i + 2^j$ (such as $n=10$), then $(x \,/\, n) \times n = x$ (exactly).

Unsurprisingly, this theorem does not apply if $x \,/\, n$ is computed as Times[x, Power[n, -1]].

According to the theorem, all the totals below should be zero.

(*** N.B.: Vectorized FAILS ***)
(Divide[#, 10.]*10. - #) &[
   N@Range@32] // Unitize // Total
(*** WORKS on short packed arrays?  ***)
(Divide[#, 10.]*10. - #) &[
   N@Range@31] // Unitize // Total
(*** No, wait! Vectorized sometimes WORKS!!  ***)
(Divide[#, ConstantArray[10., Length@#]]*10. - #) &[
   N@Range@32] // Unitize // Total
(*** WORKS on unpacked arrays? ***)
(Divide[#, 10.]*10. - #) &[
   N@Range@32 // Developer`FromPackedArray] // Unitize // Total
(*** Nope, FAILS on large arrays! Maybe gets repacked/compiled? ***)
(Divide[#, 10.]*10. - #) &[
   N@Range@251 // Developer`FromPackedArray] // Unitize // Total
(*
  9   -- 32 numbers:    FAILURE
  0   -- 31 numbers:    SUCCESS
  0   -- Double array:  SUCCESS
  0   -- Unpacked 32:   SUCCESS
  65  -- Unpacked >250: FAILURE
*)

Perhaps we're seeing some patterns, but the Trace the command changes the computational environment and results:

ReleaseHold@Last@Trace[Divide[#, 10.]] - Divide[#, 10.] &[
   N@Range@32] // Unitize // Total
(*  11  *)

ReleaseHold@Last@Trace[Divide[#, 10.]*10. - #] &@[
   N@Range@32] // Unitize // Total
(*  0 *)

These are not consistent with the first result, although both are computed from the same input N@Range@32. Apparently Trace prevents Divide from converting to multiplying by the reciprocal.

On my system (Intel i7/Macbook Pro), we can address the two failures with system options:

SetSystemOptions[
 "ParallelOptions" -> (* originally {128, 32, 32} *)
  "VectorVendorLengthThresholds" -> {128, 32, Infinity}]

SetSystemOptions[
 "PackedArrayOptions" -> (* originally 250 *)
  "ListableAutoPackLength" -> Infinity]

It's been noted elsewhere that Divide[a, b] is faster than a/b:

First@AbsoluteTiming[Divide[##]] & @@ RandomReal[1, {2, 10^7}]
(*  0.026914  *)

First@AbsoluteTiming[#1/#2] & @@ RandomReal[1, {2, 10^7}]
(*  0.05279  *)

It seems completely consistent with Divide[a, b] having half as many FLOPS as a/b. It's too bad there isn't a reliable way to implement it.

---

Addendum: A criticism of x * (1/y).

The recently well-publicized flaw in floating-point division on the Pentium chip has prompted the question whether $\tt x/y$ may simply be replaced by $\tt x*(1/y)$. This note focuses on the example $y=x$ illustrating that this proposed fix would no longer conform to the IEEE standard. --- Edelman, "When is $x*(1/x)\ne1$" (1994)