[Octave-bug-tracker] [bug #60786] mpower: negative scalar to power of matrix containing non-integers returns nan

[A.R. Burgers]

Follow-up Comment #7, bug #60786 (project octave):

Here is test script:


a_vals = [2, -2, 0];                                                          
               
a_vals = [0];                                                                 
               
a_vals = [-2];
for a_s = a_vals
                       a_fs = single(a_s);
  b_m  = 0.5*eye(2);   b_fm = single(b_m);
  b_cm = complex(b_m); b_fcm = complex(b_fm);
  b_sm = sparse(b_m);  b_scm = sparse(b_cm);

  % powers of float/double matrices
  pow_s_m = a_s ^ b_m
  pow_s_fm = a_s ^ b_fm
  pow_fs_m = a_fs ^ b_m
  pow_fs_fm = a_fs ^ b_fm

  % powers of float/double complex matrices
  pow_s_cm = a_s ^ b_cm
  pow_s_fcm = a_s ^ b_fcm
  pow_fs_cm = a_fs ^ b_cm
  pow_fs_cfm = a_fs ^ b_fcm

  % powers of sparse  matrices
  pow_s_sm = a_s ^ b_sm
  pow_s_scm = a_s ^ b_scm
  if 0
    % not implemented
    pow_fs_sm = a_fs ^ b_sm
  end
end


which produces this after your patch, showing the sparse and real matrices
have been fixed, but not the complex ones.


pow_s_m =                                                                     
                          
   0.0000 + 1.4142i        0 +      0i
        0 +      0i   0.0000 + 1.4142i

pow_s_fm =
  -0.0000 + 1.4142i        0 +      0i
        0 +      0i  -0.0000 + 1.4142i

pow_fs_m =
  -0.0000 + 1.4142i        0 +      0i
        0 +      0i  -0.0000 + 1.4142i

pow_fs_fm =
  -0.0000 + 1.4142i        0 +      0i
        0 +      0i  -0.0000 + 1.4142i

pow_s_cm =
   NaN - NaNi   NaN - NaNi
   NaN - NaNi   NaN - NaNi

pow_s_fcm =
   NaN - NaNi   NaN - NaNi
   NaN - NaNi   NaN - NaNi

pow_fs_cm =
   NaN - NaNi   NaN - NaNi
   NaN - NaNi   NaN - NaNi

pow_fs_cfm =
   NaN - NaNi   NaN - NaNi
   NaN - NaNi   NaN - NaNi

pow_s_sm =
   0.0000 + 1.4142i        0 +      0i
        0 +      0i   0.0000 + 1.4142i

pow_s_scm =
   0.0000 + 1.4142i        0 +      0i
        0 +      0i   0.0000 + 1.4142i


also there was already an inconsistency in xpow whether EIG is called in- or
outside the try block. In xpow (double a, const Matrix& b) it is called
inside, otherwise outside. On first glance outside seems OK, since what fails
and needs to be caught is the Q.inverse.

I wonder whether the separate code path for a>0 is necessary for the real
matrix case, as later on for a>0 the real part of the result is taken anyway
before returning.

Additionally, would it make sense to have early exits when a=0.



    _______________________________________________________

Reply to this item at:

  <https://savannah.gnu.org/bugs/?60786>

_______________________________________________
  Message sent via Savannah
  https://savannah.gnu.org/