At least we understand what behavior needs to be implemented.  At a minimum, Octave needs to check whether each operand contains a NaN and it needs to specifically guarantee that those locations result in NaN at the output.  The subsequent decision, about whether to convert to full if it would save space, should probably be left up to Octave, or to the sparse_auto_mutate code.

Pseudo-code which might work for Z = X op Y:

  list1 = isnan (X);
  list2 = isnan (Y);
  tmp = X op Y;
  tmp(list1) = NaN;
  tmp(list2) = NaN;
  Z = tmp;

Of course, all of that would need to be in C++.

Another choice would be to look carefully at the code which performs the operation 'op'.  For sparse matrices, I believe the code is equivalent to

nonzero1 = nonzeros (X);
nonzero2 = nonzeros (Y);
elements_to_operate_on = intersection (nonzero1, nonzero2);
Z = X(elements_to_operate_on) op Y(elements_to_operate_on);

The "intersection" function for sparse matrices in C++ could be changed to include elements where only one of the matrices had a NaN.

@Carlo: The memory problem already exists.  After one of the operations involving NaN, every single element will become NaN and require storage.  In fact, the situation will be worse than a full matrix because there will be the extra overhead of the sparse matrix implementation.

In general, even without sparse_auto_mutate set to true, Octave will convert to full for operations that are likely to result in a full matrix.  From the documentation,

22.1.4.2 Return Types of Operators and Functions

The two basic reasons to use sparse matrices are to reduce the memory usage and to not have
to do calculations on zero elements. The two are closely related in that the computation
time on a sparse matrix operator or function is roughly linear with the number of nonzero
elements.

Therefore, there is a certain density of nonzero elements of a matrix where it no longer
makes sense to store it as a sparse matrix, but rather as a full matrix. For this reason
operators and functions that have a high probability of returning a full matrix will always
return one. For example adding a scalar constant to a sparse matrix will almost always
make it a full matrix, and so the example,

speye (3) + 0
⇒
 1 0 0
0 1 0
0 0 1
returns a full matrix as can be seen.

For addition and subtraction, the sparse matrices are already converted to full and the NaN is added/subtracted and the operation follows IEEE guidelines.  For Sparse Matrix OPERATION Scalar NaN, the code would be easy to modify to check for NaN and then return

NaN (size (Xsparse))

For Matrix OPERATION Matrix operations it becomes a bit annoying because one would have to check for isnan() on every element.  This may be what Matlab does.  Can you test the following?

xs = sparse ([1 0; 0 2]);
ys = sparse ([1 NaN; 0 0]);
xs .* ys

Any obvious solution would be for NaN to force all matrices into a full implementation.

What does Matlab do for

Msparse = sparse ([1 0; 0 3]);
NaN .* Msparse

Confirmed on a recent (5/31/12) development tip.

The issue is more general than just divide so I have changed the bug Summary.

At present, I can verify that multiply and divide are affected while addition and subtraction behaves correctly.

The issue appears to be a side effect of the way that sparse matrices improve calculation performance.  In general, math operations like multiply and divide only need to be performed on the non-zero elements of a sparse matrix.  However, a NaN is different and should affect all elements of the matrix.

Below is some simple code to illustrate the problem.

Mfull = [1 0; 0 3];
Msparse = sparse (Mfull);

NaN .* Mfull
ans =

   NaN   NaN
   NaN   NaN

NaN .* Msparse
ans =

Compressed Column Sparse (rows = 2, cols = 2, nnz = 2 [50%])

  (1, 1) -> NaN
  (2, 2) -> NaN

NaN ./ Mfull
ans =

   NaN   NaN
   NaN   NaN

NaN ./ Msparse

ans =

   NaN  -Inf
  -Inf   NaN