DGtal 1.3.0
Loading...
Searching...
No Matches
testMPolynomial.cpp
Go to the documentation of this file.
1
31#include <iostream>
32#include <iomanip>
33#include <sstream>
34#include "DGtal/base/Common.h"
35#include "DGtal/math/MPolynomial.h"
36#include "DGtal/io/readers/MPolynomialReader.h"
38
39using namespace std;
40using namespace DGtal;
41
43// Functions for testing class Signal.
45
49template <typename Ring>
52{
54 = mmonomial<Ring>( 3, 1, 0 )
55 + mmonomial<Ring>( 1, 0, 3 )
56 + mmonomial<Ring>( 0, 3, 1 )
57 + mmonomial<Ring>( 0, 0, 3 )
58 + ((Ring)5) * mmonomial<Ring>( 0, 0, 1 );
59 return p;
60}
61
62double
63durchblickC( const double & x, const double & y, const double z )
64{
65 return x*x*x*y + x*z*z*z + y*y*y*z + z*z*z + 5.0 * z;
66}
67
72bool testMPolynomialSpeed( double step = 0.01 )
73{
74 unsigned int nbok = 0;
75 unsigned int nb = 0;
76
77 trace.beginBlock ( "Testing block ... Evaluation speed of mpolynomials (naive)" );
78 trace.info() << setprecision( 15 ) << "step is " << step << std::endl;
79 trace.info() << "approximately " << 8.0/(step*step*step) << " computations." << std::endl;
80 MPolynomial<3, double> P = durchblick<double>();
81 double total = 0.0;
82 for ( double x = -1.0; x < 1.0; x += step )
83 {
84 for ( double y = -1.0; y < 1.0; y += step )
85 {
86 for ( double z = -1.0; z < 1.0; z += step )
87 total += P(x)(y)(z);
88 }
89 }
90 trace.info() << "Total = " << total << std::endl;
92
93 trace.beginBlock ( "Testing block ... Evaluation speed of mpolynomials" );
94 //MPolynomial<3, double> P = durchblick<double>();
95 double total1 = 0.0;
96 for ( double x = -1.0; x < 1.0; x += step )
97 {
98 MPolynomial<2, double> PX = P( x );
99 for ( double y = -1.0; y < 1.0; y += step )
100 {
101 MPolynomial<1, double> PXY = PX( y );
102 for ( double z = -1.0; z < 1.0; z += step )
103 total1 += PXY( z );
104 }
105 }
106 trace.info() << "Total1 = " << total1 << std::endl;
107 trace.endBlock();
108
109 trace.beginBlock ( "Testing block ... Same computation in C." );
110 double total2 = 0.0;
111 for ( double x = -1.0; x < 1.0; x += step )
112 {
113 for ( double y = -1.0; y < 1.0; y += step )
114 {
115 for ( double z = -1.0; z < 1.0; z += step )
116 total2 += durchblickC( x, y, z );
117 }
118 }
119 trace.info() << "Total2 = " << total2 << std::endl;
120 trace.endBlock();
121 nbok += fabs( total1 - total ) < 1e-8 ? 1 : 0;
122 nb++;
123 trace.info() << "(" << nbok << "/" << nb << ") "
124 << "fabs( total1 - total ) < 1e-8" << std::endl;
125 nbok += fabs( total2 - total ) < 1e-8 ? 1 : 0;
126 nb++;
127 trace.info() << "(" << nbok << "/" << nb << ") "
128 << "fabs( total2 - total ) < 1e-8" << std::endl;
129
130 trace.info() << "For information, ImaGene::Polynomial3 takes 164ms for step=0.01 and 1604ms for step = 0.005." << std::endl;
131 return nbok == nb;
132}
133
139{
140 unsigned int nbok = 0;
141 unsigned int nb = 0;
142
143 trace.beginBlock ( "Testing block ..." );
144
145 MPolynomial<2, int> f = mmonomial<int>(1, 2) + 3 * mmonomial<int>(4, 5);
146 trace.info() << f << std::endl;
147 int e = f(4)(2);
148 trace.info() << e << std::endl;
149 nbok += e == 24592 ? 1 : 0;
150 nb++;
151 trace.info() << derivative<1>(f) << std::endl;
152 MPolynomial<1, int> g = f(2), h = f(3);
153 trace.info() << g << " and " << h << std::endl;
154 trace.info() << gcd<double>(g, h) << std::endl; // cast polynomials to MPolynomial<1, double> first; otherwise,
155 // result will be incorrect since int is not a field
156 trace.info() << gcd(g, h) << std::endl; // to prove our point, check for yourself that the result is
157 // X_0^5 instead of X_0^2
158 nbok += gcd<double>(g,h) == mmonomial<double>(2) ? 1 : 0;
159 nb++;
160 MPolynomial<1, double> l = f(mmonomial<double>(1) - 1)(mmonomial<double>(1) - 2);
161 trace.info() << l << std::endl;
162 trace.info() << derivative<0>(l) << std::endl;
163 trace.info() << gcd(l, derivative<0>(l)) << " (gcd of two previous polys is 1)" << std::endl;
164 nbok += gcd(l, derivative<0>(l)) == 1.0 * mmonomial<double>(0) ? 1 : 0;
165 nb++;
166 trace.info() << "Durchblick (x3y+xz3+y3z+z3+5z)= " << durchblick<double>() << std::endl;
167 trace.endBlock();
168 MPolynomial<3, double> Q = mmonomial<double>( 0, 0, 0 )
169 + mmonomial<double>( 1, 2, 0 ) + mmonomial<double>( 4, 1, 1 );
170 std::cout << "Q(x,y,z)=1+xy^2+x^4yz = " << Q << std::endl;
171 std::cout << " degree = " << Q.degree() << std::endl;
172 std::cout << " leading = " << Q.leading() << std::endl;
173 std::cout << " Q[0] = " << Q[ 0 ] << std::endl;
174 std::cout << " Q[1] = " << Q[ 1 ] << std::endl;
175 std::cout << " Q[2] = " << Q[ 2 ] << std::endl;
176 std::cout << " Q[3] = " << Q[ 3 ] << std::endl;
177 std::cout << " Q[4] = " << Q[ 4 ] << std::endl;
178 std::cout << " dQ/dx = " << derivative<0>(Q) << std::endl;
179 std::cout << " dQ/dy = " << derivative<1>(Q) << std::endl;
180 std::cout << " dQ/dz = " << derivative<2>(Q) << std::endl;
182 P = Xe_k<3,double>( 2, 7 ) + Xe_k<3,double>( 1, 3 );
183 trace.info() << "P=" << P << std::endl;
184 return nbok == nb;
185}
186
187
189// Standard services - public :
190
191int main( int /*argc*/, char** /*argv*/ )
192{
193 trace.beginBlock ( "Testing class MPolynomial" );
194
195 bool res = testMPolynomial()
196 && testMPolynomialSpeed( 0.05 );
197 trace.emphase() << ( res ? "Passed." : "Error." ) << endl;
198 trace.endBlock();
199 return res ? 0 : 1;
200}
201// //
Aim: Represents a multivariate polynomial, i.e. an element of , where K is some ring or field.
Definition: MPolynomial.h:955
int degree() const
Definition: MPolynomial.h:1109
const MPolyNM1 & leading() const
Definition: MPolynomial.h:1119
void beginBlock(const std::string &keyword="")
std::ostream & emphase()
std::ostream & info()
double endBlock()
DGtal is the top-level namespace which contains all DGtal functions and types.
Trace trace
Definition: Common.h:154
MPolynomial< 1, Ring, Alloc > gcd(const MPolynomial< 1, Ring, Alloc > &f, const MPolynomial< 1, Ring, Alloc > &g)
Definition: MPolynomial.h:2047
STL namespace.
int main()
Definition: testBits.cpp:56
bool testMPolynomial()
MPolynomial< 3, Ring, std::allocator< Ring > > durchblick()
bool testMPolynomialSpeed(double step=0.01)
double durchblickC(const double &x, const double &y, const double z)