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