Re: log2 replacement

[Jaroslav Hajek]

On Wed, Apr 30, 2008 at 9:30 AM, John W. Eaton <address@hidden> wrote:
>
> 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.
>

OK, no problem. I think the same is true for log2 and perhpas other
functions, isn't it? Perhaps a "cleanup"
can be done.


>  | +     * 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.
>

Sorry. I just had no more inspiration, I think this was just a
redundant letter (log2 instead of log).

>  | +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?

I'm not sure. I think it's better if the complex and real version have
the same name. Can you do that with frexp, which is external "C"?
Also, I think some functions in <math.h> can actually be macros, in
which case the overload would not be possible.

>
>  | +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?
>

No, there isn't. I forgot that "Matrix" means just 2D array.


>  Thanks,
>
>  jwe
>



-- 
RNDr. Jaroslav Hajek
computing expert
Aeronautical Research and Test Institute (VZLU)
Prague, Czech Republic
url: www.highegg.matfyz.cz