00001
00002 #include <basix/double.hpp>
00003 #include <basix/int.hpp>
00004 #include <basix/vector.hpp>
00005 #include <basix/port.hpp>
00006 #include <basix/literal.hpp>
00007 #include <numerix/integer.hpp>
00008 #include <numerix/modular.hpp>
00009 #include <numerix/modular_integer.hpp>
00010 #include <numerix/rational.hpp>
00011 #include <numerix/complex.hpp>
00012 #include <numerix/complex_double.hpp>
00013 #include <algebramix/vector_unrolled.hpp>
00014 #include <algebramix/vector_simd.hpp>
00015 #include <algebramix/vector_modular.hpp>
00016 #include <analyziz/vector_double.hpp>
00017 #include <basix/compound.hpp>
00018 #include <basix/mmx_syntax.hpp>
00019 #include <basix/lisp_syntax.hpp>
00020 #include <basix/cpp_syntax.hpp>
00021 #include <basix/syntactic.hpp>
00022 #include <algebramix/polynomial.hpp>
00023 #include <algebramix/polynomial_polynomial.hpp>
00024 #include <algebramix/polynomial_integer.hpp>
00025 #include <algebramix/polynomial_rational.hpp>
00026 #include <algebramix/polynomial_modular.hpp>
00027 #include <algebramix/polynomial_modular_integer.hpp>
00028 #include <algebramix/polynomial_complex.hpp>
00029 #include <algebramix/polynomial_schonhage.hpp>
00030 #include <analyziz/polynomial_numeric.hpp>
00031 #include <analyziz/polynomial_double.hpp>
00032 #include <analyziz/solver.hpp>
00033 #include <analyziz/solver_aberth.hpp>
00034 #include <algebramix/permutation.hpp>
00035 #include <basix/row_tuple.hpp>
00036 #include <algebramix/matrix.hpp>
00037 #include <algebramix/matrix_ring_naive.hpp>
00038 #include <algebramix/matrix_integer.hpp>
00039 #include <algebramix/matrix_modular_integer.hpp>
00040 #include <algebramix/matrix_quotient.hpp>
00041 #include <algebramix/matrix_complex.hpp>
00042 #include <analyziz/matrix_jacobian.hpp>
00043 #include <analyziz/matrix_double.hpp>
00044 #include <analyziz/eigen.hpp>
00045 #include <analyziz/eigen_homotopy.hpp>
00046 #include <multimix/multivariate_coordinates.hpp>
00047 #include <multimix/multivariate_monomial.hpp>
00048 #include <multimix/multivariate_polynomial.hpp>
00049 #include <multimix/sparse_polynomial_integer.hpp>
00050 #include <multimix/sparse_polynomial_rational.hpp>
00051 #include <multimix/sparse_polynomial_modular.hpp>
00052 #include <multimix/sparse_polynomial_modular_integer.hpp>
00053 #include <numerix/floating.hpp>
00054 #include <numerix/ball.hpp>
00055 #include <numerix/ball_complex.hpp>
00056 #include <analyziz/taylor.hpp>
00057 #include <continewz/ode_taylor.hpp>
00058 #include <basix/routine.hpp>
00059 #include <basix/alias.hpp>
00060 #include <basix/glue.hpp>
00061
00062 #define double_literal(x) as_double (as_string (x))
00063 #define int_literal(x) as_int (as_string (x))
00064 #define is_generic_literal is<literal>
00065 #define gen_literal_apply(f,v) gen (as<generic> (f), v)
00066 #define gen_literal_access(f,v) access (as<generic> (f), v)
00067 #define is_generic_compound is<compound>
00068 #define compound_arguments(x) cdr (as_vector (x))
00069 #define gen_compound_apply(f,v) gen (as<generic> (f), v)
00070 namespace mmx {
00071 template<typename C> polynomial<C>
00072 polynomial_reverse (const vector<C>& v) {
00073 return polynomial<C> (reverse (v)); }
00074
00075 template<typename C> polynomial<modular<modulus<C>, modular_local> >
00076 as_polynomial_modular (const polynomial<C>& f, const modulus<C>& p) {
00077 modular<modulus<C>, modular_local>::set_modulus (p);
00078 return as<polynomial<modular<modulus<C>, modular_local> > > (f); }
00079
00080 template<typename C> vector<generic>
00081 wrap_subresultants (const polynomial<C>& f, const polynomial<C>& g) {
00082 return as<vector<generic> > (subresultants (f, g)); }
00083
00084 }
00085 namespace mmx { POLYNOMIAL_GENERIC_USES_SCHONHAGE }
00086 #define identity_matrix_integer identity_matrix<integer>
00087 #define hilbert_matrix_rational hilbert_matrix<rational>
00088
00089 namespace mmx {
00090 template<typename C> matrix<C>
00091 matrix_new (const vector<row_tuple<C> >& t) {
00092 if (N(t) == 0) return matrix<C> ();
00093 nat i, j, rows= N(t), cols= N(t[0]);
00094 C dummy= zero_cst<C> ();
00095 matrix<C> r (dummy, rows, cols);
00096 for (i=0; i<rows; i++) {
00097 ASSERT (N(t[i]) == cols, "unequal row lengths");
00098 for (j=0; j<cols; j++)
00099 r(i,j)= t[i][j];
00100 }
00101 return r;
00102 }
00103
00104 template<typename C> matrix<C>
00105 matrix_new (const vector<C>& t) {
00106 nat j, rows= 1, cols= N(t);
00107 C dummy= zero_cst<C> ();
00108 matrix<C> r (dummy, rows, cols);
00109 for (j=0; j<cols; j++)
00110 r(0,j)= t[j];
00111 return r;
00112 }
00113
00114 template<typename C> vector<generic>
00115 wrap_column_reduced_echelon_with_permutation (const matrix<C>& m) {
00116 permutation permut;
00117 generic tp=as<generic> (column_reduced_echelon (m, permut));
00118 return vec (tp, as<generic> (permut));
00119 }
00120
00121 template<typename C> vector<generic>
00122 wrap_column_reduced_echelon_with_transform (const matrix<C>& m) {
00123 matrix<C> k;
00124 generic tp=as<generic> (column_reduced_echelon (m, k));
00125 return vec (tp, as<generic> (k));
00126 }
00127
00128 template<typename C> vector<generic>
00129 wrap_row_reduced_echelon_with_transform (const matrix<C>& m) {
00130 matrix<C> k;
00131 generic tp=as<generic> (row_reduced_echelon (m, k));
00132 return vec (tp, as<generic> (k));
00133 }
00134 }
00135
00136 #define mmx_coordinate multivariate_coordinate<>
00137 #define mmx_coordinates multivariate_coordinates<>
00138 #define mv_monomial multivariate<monomial<> >
00139
00140 #define mv_polynomial(C) multivariate<sparse_polynomial<C> >
00141
00142
00143 namespace mmx {
00144 #define Polynomial \
00145 multivariate<sparse_polynomial<modular<modulus<C>, modular_local> > >
00146 template<typename C> Polynomial
00147 as_mv_polynomial_modular (const mv_polynomial(C)& f, const modulus<C>& p) {
00148 modular<modulus<C>, modular_local>::set_modulus (p);
00149 return as<Polynomial> (f); }
00150 #undef Polynomial
00151 }
00152
00153 namespace mmx {
00154 static vector<integer>
00155 GLUE_1 (const vector<int> &arg_1) {
00156 return as<vector<integer> > (arg_1);
00157 }
00158
00159 static vector<rational>
00160 GLUE_2 (const vector<int> &arg_1) {
00161 return as<vector<rational> > (arg_1);
00162 }
00163
00164 static vector<complex<rational> >
00165 GLUE_3 (const vector<int> &arg_1) {
00166 return as<vector<complex<rational> > > (arg_1);
00167 }
00168
00169 static vector<double>
00170 GLUE_4 (const vector<integer> &arg_1) {
00171 return as<vector<double> > (arg_1);
00172 }
00173
00174 static vector<int>
00175 GLUE_5 (const vector<integer> &arg_1) {
00176 return as<vector<int> > (arg_1);
00177 }
00178
00179 static vector<complex<double> >
00180 GLUE_6 (const vector<integer> &arg_1) {
00181 return as<vector<complex<double> > > (arg_1);
00182 }
00183
00184 static vector<generic>
00185 GLUE_7 (const vector<mv_polynomial(rational) > &arg_1, const vector<mmx_coordinate> &arg_2, const vector<double> &arg_3) {
00186 return solve_ode_taylor (arg_1, arg_2, arg_3);
00187 }
00188
00189 static matrix<generic>
00190 GLUE_8 (const vector<mv_polynomial(rational) > &arg_1, const vector<mmx_coordinate> &arg_2, const vector<double> &arg_3) {
00191 return first_variation_taylor (arg_1, arg_2, arg_3);
00192 }
00193
00194 static vector<generic>
00195 GLUE_9 (const vector<mv_polynomial(rational) > &arg_1, const vector<mmx_coordinate> &arg_2, const vector<double> &arg_3, const double &arg_4) {
00196 return integrate_ode_taylor (arg_1, arg_2, arg_3, arg_4);
00197 }
00198
00199 static vector<generic>
00200 GLUE_10 (const vector<mv_polynomial(rational) > &arg_1, const vector<mmx_coordinate> &arg_2, const vector<complex<double> > &arg_3) {
00201 return solve_ode_taylor (arg_1, arg_2, arg_3);
00202 }
00203
00204 static matrix<generic>
00205 GLUE_11 (const vector<mv_polynomial(rational) > &arg_1, const vector<mmx_coordinate> &arg_2, const vector<complex<double> > &arg_3) {
00206 return first_variation_taylor (arg_1, arg_2, arg_3);
00207 }
00208
00209 static vector<generic>
00210 GLUE_12 (const vector<mv_polynomial(rational) > &arg_1, const vector<mmx_coordinate> &arg_2, const vector<complex<double> > &arg_3, const complex<double> &arg_4) {
00211 return integrate_ode_taylor (arg_1, arg_2, arg_3, arg_4);
00212 }
00213
00214 static vector<double>
00215 GLUE_13 (const vector<rational> &arg_1) {
00216 return as<vector<double> > (arg_1);
00217 }
00218
00219 static vector<complex<double> >
00220 GLUE_14 (const vector<rational> &arg_1) {
00221 return as<vector<complex<double> > > (arg_1);
00222 }
00223
00224 static vector<complex<double> >
00225 GLUE_15 (const vector<complex<rational> > &arg_1) {
00226 return as<vector<complex<double> > > (arg_1);
00227 }
00228
00229 void
00230 glue_ode_taylor_double () {
00231 static bool done = false;
00232 if (done) return;
00233 done = true;
00234 call_glue (string ("glue_matrix_double"));
00235 call_glue (string ("glue_mvpolynomial_rational"));
00236 static alias<int> mmx_significant_digits_alias = global_alias (((int&) mmx_significant_digits));
00237 define_constant<alias<int> > ("significant_digits", mmx_significant_digits_alias);
00238 static alias<int> mmx_bit_precision_alias = global_alias (((int&) mmx_bit_precision));
00239 define_constant<alias<int> > ("bit_precision", mmx_bit_precision_alias);
00240 static alias<int> mmx_discrepancy_alias = global_alias (((int&) mmx_discrepancy));
00241 define_constant<alias<int> > ("discrepancy", mmx_discrepancy_alias);
00242 static alias<bool> mmx_pretty_exponents_alias = global_alias (((bool&) mmx_pretty_exponents));
00243 define_constant<alias<bool> > ("pretty_exponents", mmx_pretty_exponents_alias);
00244 static alias<int> mmx_ode_order_alias = global_alias (((int&) mmx_ode_order));
00245 define_constant<alias<int> > ("ode_order", mmx_ode_order_alias);
00246 static alias<int> mmx_ode_precision_alias = global_alias (((int&) mmx_ode_precision));
00247 define_constant<alias<int> > ("ode_precision", mmx_ode_precision_alias);
00248 define_converter (":>", GLUE_1, PENALTY_INCLUSION);
00249 define_converter (":>", GLUE_2, PENALTY_INCLUSION);
00250 define_converter (":>", GLUE_3, PENALTY_HOMOMORPHISM);
00251 define_converter (":>", GLUE_4, PENALTY_INCLUSION);
00252 define_converter (":>", GLUE_5, PENALTY_INCLUSION);
00253 define_converter (":>", GLUE_6, PENALTY_HOMOMORPHISM);
00254 define ("solve_ode_taylor", GLUE_7);
00255 define ("first_variation_taylor", GLUE_8);
00256 define ("integrate_ode_taylor", GLUE_9);
00257 define ("solve_ode_taylor", GLUE_10);
00258 define ("first_variation_taylor", GLUE_11);
00259 define ("integrate_ode_taylor", GLUE_12);
00260 define_converter (":>", GLUE_13, PENALTY_INCLUSION);
00261 define_converter (":>", GLUE_14, PENALTY_HOMOMORPHISM);
00262 define_converter (":>", GLUE_15, PENALTY_INCLUSION);
00263 }
00264 }
