Re: log2 replacement

[John W. Eaton]

On 30-Apr-2008, Jaroslav Hajek wrote:

| On Tue, Apr 29, 2008 at 8:19 PM, John W. Eaton <address@hidden> wrote:
| > On 25-Apr-2008, Jaroslav Hajek wrote:
| >
| >  | hello,
| >  |
| >  | please consider the attached changeset.
| >  | The current log2 behaviour with two outputs (`[f, e] = log2(x)') is
| >  | inconsistent with Matlab in a number of cases, (x = Inf, x = NaN) or
| >  | just plain wrong (x = realmax). Also, on most of the systems the
| >  | 2-output version can be mapped directly to frexp (a C89 function),
| >  | which is probably what Matlab also does. I have implemented the
| >  | single-output log2 as a mapper for octave_value (aka log, log10 etc),
| >  | and the two-output version as an overload in liboctave/lo-mappers.cc
| >  | with a replacement if frexp is not detected. The DEFUN is in
| >  | DLD-FUNCTIONS rather than mappers.cc directly, because it is not a
| >  | simple mapper and thus did not seem to fit "aesthetically" to
| >  | mappers.cc (but that can be easily changed).
| >  |
| >  | Note that the 2-output version is also extended beyond matlab
| >  | behaviour for complex input - while matlab ignores the imaginary part
| >  | (with a warning), this new version gives the result with f complex,
| >  | satisfying 1/2 <= abs(f) < 1, e integer - because it makes good sense
| >  | also for complex numbers.
| >
| >  With this patch, I see the following error with make check:
| >
| >   src/DLD-FUNCTIONS/hex2num.cc ........................... PASS    2/2
| >   src/DLD-FUNCTIONS/log2.cc .............................. PASS    0/2    
FAIL 2
| >   src/DLD-FUNCTIONS/lookup.cc ............................ PASS   16/16
| >   src/DLD-FUNCTIONS/lsode.cc ............................. PASS    5/5
| >   src/DLD-FUNCTIONS/luinc.cc ............................. PASS    2/2
| >   src/DLD-FUNCTIONS/matrix_type.cc ....................... PASS   51/51
| >   src/DLD-FUNCTIONS/max.cc ............................... PASS   12/12
| >   src/DLD-FUNCTIONS/qr.cc ................................ PASS   21/21
| >   src/DLD-FUNCTIONS/quad.cc .............................. PASS    6/6
| >   src/DLD-FUNCTIONS/rand.cc .............................. PASS   57/57
| >   src/DLD-FUNCTIONS/regexp.cc ............................ PASS   81/81
| >   src/DLD-FUNCTIONS/time.cc .............................. PASS   13/13
| >   src/DLD-FUNCTIONS/tsearch.cc ........................... PASS    6/6
| >   src/data.cc ............................................ PASS   98/98
| >   src/mappers.cc ......................................... PASS    7/7
| >   src/ov-fcn-handle.cc ...................................panic: 
Segmentation fault -- stopping myself...
| >  make[2]: *** [check] Segmentation fault
| >  make[2]: Leaving directory `/scratch/jwe/build/octave-tmp/test'
| >  make[1]: *** [check] Error 2
| >  make[1]: Leaving directory `/scratch/jwe/build/octave-tmp'
| >  make: *** [check] Error 2
| >
| >  Does it work properly for you?
| >
| >  jwe
| >
| 
| OK, this one seems to work. There was a type in the test code. Dunno
| what caused the fcn-handle failure; I've decided to put log2 into
| src/data.cc anyway (some similar math functions are also there, like
| atan2, hypot etc.), and it disappeared.

Thanks.  I'm still seeing the crash, but I don't think it is the fault
of your patch, so I applied it with a few changes.

On my system, the crash happens for things like:

  f = @log2;
  save ("-text", "foo.dat", "f");

I guess we were just "lucky" that there happened to be a test that did
this.

Also, I'm still seeing a failure in the data.cc tests, but again, I
don't think this is the fault of your patch.  I'm seeing this strange
output in fntests.log:

    ***** test
   [f, e] = log2 (complex (zeros (3, 2), [0,-1; 2,-4; Inf,-Inf]));
   assert (f, complex (zeros (3, 2), [0,-0.5; 0.5,-0.5; Inf,-Inf]));
   assert (e, [0,1; 2,3; 0,0]);
   assert (size (fmod (zeros (0, 2), zeros (0, 2))), [0, 2])
   assert (size (fmod (rand (2, 3, 4), zeros (2, 3, 4))), [2, 3, 4])
   assert (size (fmod (rand (2, 3, 4), 1), [2, 3, 4])
   assert (size (fmod (1, rand (2, 3, 4)), [2, 3, 4])
   assert (size (fmod (1, 2), [1, 1])
  !!!!! test failed
  parse error near line 0 of file /home/jwe/src/octave/test/fntests.m

    syntax error

  >>> function  __test__( )

   [f, e] = log2 (complex (zeros (3, 2), [0,-1; 2,-4; Inf,-Inf]));
   assert (f, complex (zeros (3, 2), [0,-0.5; 0.5,-0.5; Inf,-Inf]));
   assert (e, [0,1; 2,3; 0,0]);
   assert (size (fmod (zeros (0, 2), zeros (0, 2))), [0, 2])
   assert (size (fmod (rand (2, 3, 4), zeros (2, 3, 4))), [2, 3, 4])
   assert (size (fmod (rand (2, 3, 4), 1), [2, 3, 4])
   assert (size (fmod (1, rand (2, 3, 4)), [2, 3, 4])
   assert (size (fmod (1, 2), [1, 1])
  endfunction
             ^

Note that somehow the tests for log2 and fmod are being combined.  I
don't have any clue yet why this is happening.

Anyway, here are some comments about the patch:

|  ### Check for nonstandard but common math functions that we need.
|  
| -AC_CHECK_FUNCS(acosh asinh atanh erf erfc exp2 log2)
| +AC_CHECK_FUNCS(acosh asinh atanh erf erfc exp2 log2 frexp)

Since frexp has been part of standard C for nearly 20 years, I think
we can assume it exists, so I deleted the check.

| +     * lo-mappers.cc (xlog2 (double)): fix bug.

In ChangeLog entries, please say what changed.  If every entry used
"fix bug" or "bug fix" or similar, the log wouldn't be very useful.

| +double xlog2 (double x, int& exp)
| +{
| +#if defined (HAVE_FREXP)
| +  return frexp (x, &exp);
| +#else
| +  if (lo_ieee_finite (x) && x != 0) 
| +    {
| +      exp = xtrunc (xlog2 (abs (x)));
| +      x *= xexp2 (-exp);
| +      if (abs(x) >= 1) // `==' should suffice in theory
| +        {
| +          x /= 2;
| +          exp++;
| +        }
| +      return x;
| +    }
| +  else
| +    {
| +      exp = 0;
| +      return x;
| +    }
| +#endif
| +}

Since we can assume that frexp exists, this function simplifies to a
trivial wrapper around frexp.  For now, I left the wrapper.  Is it
worth doing that just for the name change?

| +DEFUN (log2, args, nargout,
| +    "-*- texinfo -*-\n\
| address@hidden {Mapping Function} {} log2 (@var{x})\n\
| address@hidden {Mapping Function} {[f, e] = } log2 (@var{x})\n\
| +Compute the base-2 logarithm for each element of @var{x}.\n\
| +If called with two output arguments, split @var{x} to\n\
| +binary mantissa and exponent so that @code{1/2 <= abs(f) < 1} and\n\
| address@hidden is an integer. If @code{x = 0}, @code{f = e = 0}.\n\
| address@hidden, log10, log2, exp}\n\
| address@hidden deftypefn")
| +{
| +  octave_value_list retval;
| +  if (args.length () == 1)
| +    {
| +      if (nargout < 2)
| +        retval(0) = args(0).log2 ();
| +      else if (args(0).is_complex_matrix ())
| +        {
| +          ComplexMatrix f, x = args(0).complex_matrix_value ();
| +          Matrix e; // NOTE: could also return intNDArray ...
| +          map_2_xlog2 (x, f, e);
| +          retval (1) = e;
| +          retval (0) = f;
| +        }
| +      else
| +        {
| +          Matrix f, x = args(0).matrix_value ();
| +          Matrix e; // NOTE: could also return intNDArray ...
| +          map_2_xlog2 (x, f, e);
| +          retval (1) = e;
| +          retval (0) = f;
| +        }
| +    }
| +  else
| +    print_usage ();
| +
| +  return retval;
| +}

I changed this to check for "args(0).is_real_type",
"args(0).is_complex_type" and complain otherwise.  Also, I extract
NDArray values.  Or is there some reason to restrict it to just 2D
matrix values?

Thanks,

jwe