1 #ifndef _MPS_PRIVATE 2 #define _MPS_PRIVATE 3 #endif 4 5 #include <octave/oct.h> 6 #include <gmp.h> 7 #include <octave/error.h> 8 #include <mps/mps.h> 9 10 DEFUN_DLD(mps_regenerate, args, nargout, 11 "-*- texinfo -*- \n\ 12 @deftypefn {Loadable Function} {@var{a} =} mps_det (@var{b}, @var{P0}, ..., @var{PN})\n\ 13 Compute a secular equation equivalent to the polynomial \n\ 14 \n\ 15 P(x) = det(P0 + x*P1 + ... + x^n*Pn) \n\ 16 \n\ 17 This function will take that b_i chosen as nodes and return the respective a_i such that \n\ 18 \n\ 19 sum (a ./ (x - b)) - 1 = 0 if and only if P(x) = 0 \n\ 20 \n\ 21 @end deftypefn") 22 { 23 int nargin = args.length(); 24 mps_context * ctx = mps_context_new (); 25 26 // Input sanitizing goes here 27 28 ComplexColumnVector b = args(0).complex_column_vector_value(); 29 int degree = nargin - 2; 30 ComplexColumnVector a(b.length()); 31 32 cplx_t * fa = (cplx_t*) a.fortran_vec(); 33 cplx_t * fb = (cplx_t*) b.fortran_vec(); 34 35 cdpe_t ctmp; 36 37 // Handle some special cases here. Here we are assuming that the second argument 38 // is diagonal then it's the identity. 39 if (degree == 1 && args(2).is_diag_matrix()) 40 { 41 ComplexMatrix A = args(1).complex_matrix_value().transpose(); 42 cplx_t * A_data = (cplx_t*) A.fortran_vec(); 43 44 for (int i = 0; i < a.length(); i++) 45 { 46 long int exponent = 0; 47 cplx_t prod; 48 49 mps_fhessenberg_shifted_determinant (ctx, A_data, fb[i], A.rows(), fa[i], &exponent); 50 cplx_set (prod, cplx_one); 51 52 for (int j = 0; j < b.length(); j++) 53 { 54 if (i == j) 55 continue; 56 57 cplx_t diff; 58 cplx_sub (diff, fb[i],fb[j]); 59 cplx_mul_eq (prod, diff); 60 61 if (j % 50 == 0) 62 { 63 int pe; 64 double mod = cplx_mod (prod); 65 double new_mod = frexp (mod, &pe); 66 67 exponent -= pe; 68 cplx_div_eq_d (prod, mod / new_mod); 69 } 70 } 71 72 cplx_div_eq (fa[i], prod); 73 cdpe_set_2dl (ctmp, 1.0, exponent, 0.0, 1); 74 cdpe_get_x (prod, ctmp); 75 cplx_mul_eq (fa[i], prod); 76 } 77 } 78 79 mps_context_free (ctx); 80 return octave_value (- a); 81 } 82