# Division by zero is slow

Suppose I have two arrays of data, which I want to divide:

n = 1000;
a = RandomInteger[{1, 9}, {n, n}];
b = RandomInteger[{1, 9}, {n, n}];
First @ AbsoluteTiming @ (a/b)


0.637064

But suppose that two thirds of the values are actually zero. I'd want the ones dividing by zero to just remain zero.

mask = RandomChoice[{0, 0, 1}, {n, n}];
First @ AbsoluteTiming @ Quiet[c/d /. Indeterminate -> 0]


8.568857

Now it has become very slow. How can we speed this up?

Edit: Thanks for all the answers guys. Below are some timings:

-

SparseArray can help, given the size and nature of the mask. It's slightly faster to convert c and d to sparse arrays than to convert a and b.

mask = SparseArray@RandomChoice[{0, 0, 1}, {n, n}];

First@AbsoluteTiming[
Quiet[foo2 = Block[{Indeterminate = 0}, SparseArray[c] / SparseArray[d]]]
]
(*
0.151565
*)


Compare:

First@AbsoluteTiming[
foo1 = Quiet[c/d /. Indeterminate -> 0]
]
(*
6.492534
*)

foo2 == foo1
(*
True
*)

-
@kguler I'm not sure why you deleted your Block trick, which I've used because it's faster than any other trick I've tried. But let me know if there's a problem. My use of SparseArray is a significant difference, but Replace does not work on it. Block and Map were the only tools I had, and Map was way too slow. –  Michael E2 Oct 28 '14 at 4:37
The trick will fail when c and d are using different masks. (Block won't help because the FullForm of the output of ComplexInfinity generated by 1/0 etc. is DirectedInfinity[].) Anyway, it's the fastest method for OP's specific example so far. –  xzczd Oct 28 '14 at 5:38
@xzczd Actually it seems very odd to have different masks. In the OP's Q, it's clear that Indeterminate is to go to 0. 1/0 does not become Indeterminate and there is no indication that ComplexInfinity should become 0. (I decided it probably should not.) –  Michael E2 Oct 28 '14 at 12:15
I like this answer the best. It's a bit wordy, but clear what it is doing. It's fast and flexible (if you want to do more than a simple divide). It's not quite so fast if you want to do a Normal afterwards - about the same speed as the case with no zeroes. –  wxffles Oct 28 '14 at 21:34
Michael, just saw your comment (somehow i think @name mechanism is broken; i don't get pinged when mentioned in comments). The only reason i deleted my answer was I forgot to refresh the kernel in a long/crowded mma session and got confused with conflicting results:) (+1) –  kguler Oct 29 '14 at 3:05
SetAttributes[f, Listable];
f@__ := 0
f[a_, b_] /; b != 0 := a/b
First@AbsoluteTiming[f[c, d];]


Edit

On the same spirit, this is 50% faster:

SetAttributes[g, Listable];
g[_, 0] := 0
g[a_, b_] := Divide[a, b]
First@AbsoluteTiming[g[c, d];]

-
Unfortunately this is not particularly fast. –  Mr.Wizard Oct 28 '14 at 11:29
@Mr.Wizard Unfortunately, it looks a lot like your code too... what have you done to belisarius! :P –  rm -rf Oct 29 '14 at 5:10
@rm-rf I am taking over the world one syntax at a time. :p –  Mr.Wizard Oct 29 '14 at 5:36

Simply turn all the 0/0s into 0/1s:

ans2 = c/(1 - Sign@d + d); // AbsoluteTiming

{0.7630000, Null}


The above answer only works for specific example in the question i.e. it only handles positive lists and identical mask. For more general cases one can use:

n = 1000;
(* a and b contain non-positive elements now *)
origin := RandomInteger[{-9, 9}, {n, n}]
a = origin; b = origin;

(* c and d now use different masks *)
mask := RandomChoice[{0, 0, 1}, {n, n}];

ans2 = With[{s = Unitize@d}, Divide[s c, 1 - s + d]]; // AbsoluteTiming

{0.606000, Null}


Notice that Divide is faster than /.

-
The best general answer so far IMHO. –  Mr.Wizard Oct 28 '14 at 8:55
Note that you can use Unitize@d in place of Abs@Sign@d. –  Mr.Wizard Oct 28 '14 at 8:59
@Mr.Wizard It makes sense that one might know beforehand that 0/0 ought to be 0 - in any case, Indeterminate -> 0 is specified in the question. The question leaves open what should happen to a result like 1/0 or ComplexInfinity. One might want infinity, zero or a different result, such as a maximum value. The general answer here effectively takes the union of the masks, which is a plausible (even common) use-case. But I think one ought to be careful in stating which generalized problem it is a solution to. - +1 long before now. :) –  Michael E2 Oct 28 '14 at 12:40
Clever and fast, but I'm going to wonder what the heck it is doing when I come back to read the code a few months later. Comments are for the weak. –  wxffles Oct 28 '14 at 21:39
@wxffles I think the most probable obstacle for one to understand the code is unfamiliarity with Unitize and vectorization :) …OK, I've added a bit explanation, just a bit though. –  xzczd Nov 1 '14 at 3:29
c/(d /. (0 | 0. -> Infinity));

-

For the given example where no values in b are zero before the mask one can use:

Divide[mask*a, b]


Note the use of explicit Divide for optimum performance. Timings:

(* your data *)

(r1 = Quiet[c/d /. Indeterminate -> 0]); // AbsoluteTiming // First
(r2 = Divide[mask*a, b]);                // AbsoluteTiming // First

r1 == r2

4.578262

0.212012

True


For the general problem I cannot beat xzczd's method.

-
+! I keep forgetting about Divide.... –  Michael E2 Oct 28 '14 at 12:44
How about Divide[maskb maska a, b]? It's only a bit faster. If all factors are converted to SparseArray`, it's faster still, even including the conversion time. –  Michael E2 Oct 28 '14 at 13:15