1 #ifndef MPR_COMPLEX_H
2 #define MPR_COMPLEX_H
3 /****************************************
4 * Computer Algebra System SINGULAR *
5 ****************************************/
6
7 /*
8 * ABSTRACT - multipolynomial resultants - real floating-point numbers using gmp
9 * and complex numbers based on pairs of real floating-point numbers
10 *
11 */
12
13 //-> include & define stuff
14 // must have gmp version >= 2
15 #include "coeffs/si_gmp.h"
16 #include "coeffs/mpr_global.h"
17
18 #define ZTOF 1
19 #define QTOF 2
20 #define RTOF 3
21 #define CTOF 4
22
23 void setGMPFloatDigits( size_t digits, size_t rest );
24
25 //-> class gmp_float
26 /**
27 * @short wrapper class for GNU Multi Precision Floats
28 */
29 class gmp_float;
30 char *floatToStr( const gmp_float & r, const unsigned int oprec );
31 class gmp_float
32 {
33 public:
34 gmp_float( const int v = 0 )
35 {
36 mpf_init_set_si( t, (long)v );
37 }
gmp_float(const long v)38 gmp_float( const long v )
39 {
40 mpf_init_set_si( t, v );
41 }
gmp_float(const mprfloat v)42 gmp_float( const mprfloat v ) // double
43 {
44 mpf_init_set_d( t, v );
45 }
gmp_float(const mpf_t v)46 gmp_float( const mpf_t v )
47 {
48 mpf_init_set( t, v );
49 }
gmp_float(const mpz_t v)50 gmp_float( const mpz_t v ) // gnu mp Z
51 {
52 mpf_init( t );
53 mpf_set_z( t, v );
54 }
gmp_float(const gmp_float & v)55 gmp_float( const gmp_float & v ) // copy constructor
56 {
57 mpf_init_set( t, v.t );
58 }
59
~gmp_float()60 ~gmp_float()
61 {
62 mpf_clear( t );
63 }
64
65 inline gmp_float & operator = ( const gmp_float & a )
66 {
67 mpf_set( t, a.t );
68 return *this;
69 };
70 inline gmp_float & operator = ( const mpz_t & a )
71 {
72 mpf_set_z( t, a );
73 return *this;
74 };
75 inline gmp_float & operator = ( const mprfloat a )
76 {
77 mpf_set_d( t, (double) a );
78 return *this;
79 };
80 inline gmp_float & operator = ( const long a )
81 {
82 mpf_set_d( t, (double) a );
83 return *this;
84 };
85
86 gmp_float & operator += ( const gmp_float & a );
87 gmp_float & operator -= ( const gmp_float & a );
88 inline gmp_float & operator *= ( const gmp_float & a )
89 {
90 mpf_mul( t, t, a.t );
91 return *this;
92 };
93
94 inline gmp_float & operator /= ( const gmp_float & a )
95 {
96 mpf_div( t, t, a.t );
97 return *this;
98 };
99
neg()100 inline gmp_float & neg ( ) { mpf_neg(t,t); return *this; };
101
102 friend gmp_float operator + ( const gmp_float & a, const gmp_float & b );
103 friend gmp_float operator - ( const gmp_float & a, const gmp_float & b );
104 friend gmp_float operator * ( const gmp_float & a, const gmp_float & b );
105 friend gmp_float operator / ( const gmp_float & a, const gmp_float & b );
106
107 inline gmp_float operator ^ ( const int exp ) const
108 {
109 mpf_t b;
110 mpf_init(b);
111 mpf_pow_ui( b, this->t, (unsigned long)exp );
112 return gmp_float(b);
113 };
114
115 friend bool operator == ( const gmp_float & a, const gmp_float & b );
116 friend bool operator > ( const gmp_float & a, const gmp_float & b );
117 friend bool operator < ( const gmp_float & a, const gmp_float & b );
118 friend bool operator >= ( const gmp_float & a, const gmp_float & b );
119 friend bool operator <= ( const gmp_float & a, const gmp_float & b );
120
121 friend gmp_float operator - ( const gmp_float & a );
122
sign()123 inline int sign() // t>0:+1, t==0:0, t<0:-1
124 { return mpf_sgn( t ); };
125
126 bool isZero() const; // t == 0 ?
127 bool isOne() const; // t == 1 ?
128 bool isMOne() const; // t == -1 ?
129
130 void setFromStr(const char * in );
131
132 // access
mpfp()133 inline const mpf_t *mpfp() const { return &t; };
_mpfp()134 inline mpf_t *_mpfp() { return &t; };
135
136 inline operator double() { return mpf_get_d( t ); };
137 inline operator double() const { return mpf_get_d( t ); };
138
139 #if 0
140 inline operator int() { return (int)mpf_get_d( t ); };
141 inline operator int() const { return (int)mpf_get_d( t ); };
142 //#else
143 inline operator int() const
144 { if (mpf_fits_sint_p(t))
145 { return (int)mpf_get_si( t ); }
146 return 0;
147 };
148 #endif
149
150 private:
151 mpf_t t;
152 };
153
154
155 // built-in functions of GMP
156 gmp_float abs( const gmp_float & );
157 gmp_float sqrt( const gmp_float & );
158 gmp_float hypot( const gmp_float &, const gmp_float & );
159 //gmp_float pow( const gmp_float &, int & );
160
161 // simulated functions using double functions
162 gmp_float sin( const gmp_float & );
163 gmp_float cos( const gmp_float & );
164 gmp_float log( const gmp_float & );
165 gmp_float exp( const gmp_float & );
166
167 gmp_float max( const gmp_float &, const gmp_float & );
168
169 gmp_float numberToFloat( number num, const coeffs src );
170 gmp_float numberFieldToFloat( number num, int src );
171 //char *floatToStr( const gmp_float & r, const unsigned int oprec );
172 //<-
173
174 //-> class gmp_complex
175 /**
176 * @short gmp_complex numbers based on
177 */
178 class gmp_complex
179 {
180 private:
181 gmp_float r, i;
182
183 public:
184 gmp_complex( const gmp_float re= 0.0, const gmp_float im= 0.0 )
185 {
186 r= re;
187 i= im;
188 }
189 gmp_complex( const mprfloat re, const mprfloat im = 0.0 )
190 {
191 r= re;
192 i= im;
193 }
gmp_complex(const long re,const long im)194 gmp_complex( const long re, const long im )
195 {
196 r= re;
197 i= im;
198 }
gmp_complex(const gmp_complex & v)199 gmp_complex( const gmp_complex & v )
200 {
201 r= v.r;
202 i= v.i;
203 }
~gmp_complex()204 ~gmp_complex() {}
205
206 gmp_complex & neg ( );
207
208 friend gmp_complex operator + ( const gmp_complex & a, const gmp_complex & b );
209 friend gmp_complex operator - ( const gmp_complex & a, const gmp_complex & b );
210 friend gmp_complex operator * ( const gmp_complex & a, const gmp_complex & b );
211 friend gmp_complex operator / ( const gmp_complex & a, const gmp_complex & b );
212
213 // gmp_complex <operator> real
214 inline friend gmp_complex operator + ( const gmp_complex & a, const gmp_float b_d );
215 inline friend gmp_complex operator - ( const gmp_complex & a, const gmp_float b_d );
216 inline friend gmp_complex operator * ( const gmp_complex & a, const gmp_float b_d );
217 inline friend gmp_complex operator / ( const gmp_complex & a, const gmp_float b_d );
218
219 gmp_complex & operator += ( const gmp_complex & a );
220 gmp_complex & operator -= ( const gmp_complex & a );
221 gmp_complex & operator *= ( const gmp_complex & a );
222 gmp_complex & operator /= ( const gmp_complex & a );
223
224 inline friend bool operator == ( const gmp_complex & a, const gmp_complex & b );
225 inline friend bool operator > ( const gmp_complex & a, const gmp_complex & b );
226 inline friend bool operator < ( const gmp_complex & a, const gmp_complex & b );
227 inline friend bool operator >= ( const gmp_complex & a, const gmp_complex & b );
228 inline friend bool operator <= ( const gmp_complex & a, const gmp_complex & b );
229
230 inline gmp_complex & operator = ( const gmp_complex & a );
231 inline gmp_complex & operator = ( const gmp_float & f );
232
233 // access to real and imaginary part
real()234 inline gmp_float real() const { return r; }
imag()235 inline gmp_float imag() const { return i; }
236
real(gmp_float val)237 inline void real( gmp_float val ) { r = val; }
imag(gmp_float val)238 inline void imag( gmp_float val ) { i = val; }
239
240
isZero()241 inline bool isZero() { return (r.isZero() && i.isZero()); }
242 void SmallToZero();
243 };
244
245 // <gmp_complex> = <gmp_complex> operator <gmp_float>
246 //
247 inline gmp_complex operator + ( const gmp_complex & a, const gmp_float b_d )
248 {
249 return gmp_complex( a.r + b_d, a.i );
250 }
251 inline gmp_complex operator - ( const gmp_complex & a, const gmp_float b_d )
252 {
253 return gmp_complex( a.r - b_d, a.i );
254 }
255 inline gmp_complex operator * ( const gmp_complex & a, const gmp_float b_d )
256 {
257 return gmp_complex( a.r * b_d, a.i * b_d );
258 }
259 inline gmp_complex operator / ( const gmp_complex & a, const gmp_float b_d )
260 {
261 return gmp_complex( a.r / b_d, a.i / b_d );
262 }
263
264 // <gmp_complex> == <gmp_complex> ?
265 inline bool operator == ( const gmp_complex & a, const gmp_complex & b )
266 {
267 return ( b.real() == a.real() ) && ( b.imag() == a.imag() );
268 }
269 inline bool operator > ( const gmp_complex & a, const gmp_complex & b )
270 {
271 return ( a.real() > b.real() );
272 }
273 inline bool operator < ( const gmp_complex & a, const gmp_complex & b )
274 {
275 return ( a.real() < b.real() );
276 }
277 inline bool operator >= ( const gmp_complex & a, const gmp_complex & b )
278 {
279 return ( a.real() >= b.real() );
280 }
281 inline bool operator <= ( const gmp_complex & a, const gmp_complex & b )
282 {
283 return ( a.real() <= b.real() );
284 }
285
286
287 // <gmp_complex> = <gmp_complex>
288 inline gmp_complex & gmp_complex::operator = ( const gmp_complex & a )
289 {
290 r= a.r;
291 i= a.i;
292 return *this;
293 }
294
295 // <gmp_complex> = <gmp_complex>
296 inline gmp_complex & gmp_complex::operator = ( const gmp_float & f )
297 {
298 r= f;
299 i= (long int)0;
300 return *this;
301 }
302
303 // Returns absolute value of a gmp_complex number
304 //
abs(const gmp_complex & c)305 inline gmp_float abs( const gmp_complex & c )
306 {
307 return hypot(c.real(),c.imag());
308 }
309
310 gmp_complex sqrt( const gmp_complex & x );
311
numberToComplex(number num,const coeffs r)312 inline gmp_complex numberToComplex( number num, const coeffs r )
313 {
314 if (nCoeff_is_long_C(r))
315 {
316 return *(gmp_complex*)num;
317 }
318 else
319 {
320 return gmp_complex( numberToFloat(num, r) );
321 }
322 }
323
324 char *complexToStr( gmp_complex & c, const unsigned int oprec, const coeffs src );
325 //<-
326
327 bool complexNearZero( gmp_complex * c, int digits );
328
329 #endif /* MPR_COMPLEX_H */
330
331 // local Variables: ***
332 // folded-file: t ***
333 // compile-command-1: "make installg" ***
334 // compile-command-2: "make install" ***
335 // End: ***
336