Simbody 3.7
Loading...
Searching...
No Matches
conjugate.h
Go to the documentation of this file.
1#ifndef SimTK_SIMMATRIX_CONJUGATE_H_
2#define SimTK_SIMMATRIX_CONJUGATE_H_
3
4/* -------------------------------------------------------------------------- *
5 * Simbody(tm): SimTKcommon *
6 * -------------------------------------------------------------------------- *
7 * This is part of the SimTK biosimulation toolkit originating from *
8 * Simbios, the NIH National Center for Physics-Based Simulation of *
9 * Biological Structures at Stanford, funded under the NIH Roadmap for *
10 * Medical Research, grant U54 GM072970. See https://simtk.org/home/simbody. *
11 * *
12 * Portions copyright (c) 2005-12 Stanford University and the Authors. *
13 * Authors: Michael Sherman *
14 * Contributors: *
15 * *
16 * Licensed under the Apache License, Version 2.0 (the "License"); you may *
17 * not use this file except in compliance with the License. You may obtain a *
18 * copy of the License at http://www.apache.org/licenses/LICENSE-2.0. *
19 * *
20 * Unless required by applicable law or agreed to in writing, software *
21 * distributed under the License is distributed on an "AS IS" BASIS, *
22 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. *
23 * See the License for the specific language governing permissions and *
24 * limitations under the License. *
25 * -------------------------------------------------------------------------- */
26
60#include <complex>
61#include <iostream>
62#include <limits>
63
64using std::complex;
65
66
67// These mixed-precision real/complex operators should not be needed
68// and may have to be removed on some systems. However neither
69// gcc 3.4.4 nor VC++ 7 provide them.
70// Define the symbol below if these conflict with standard library operators.
71
72#ifndef SimTK_MIXED_PRECISION_REAL_COMPLEX_ALREADY_DEFINED
73namespace SimTK {
74// complex<float> with int, double
75inline complex<float> operator*(const complex<float>& c,int r) {return c*(float)r;}
76inline complex<float> operator*(int r,const complex<float>& c) {return (float)r*c;}
77inline complex<double> operator*(const complex<float>& c,const double& r) {return complex<double>(c)*r;}
78inline complex<double> operator*(const double& r,const complex<float>& c) {return r*complex<double>(c);}
79
80inline complex<float> operator/(const complex<float>& c,int r) {return c/(float)r;}
81inline complex<float> operator/(int r,const complex<float>& c) {return (float)r/c;}
82inline complex<double> operator/(const complex<float>& c,const double& r) {return complex<double>(c)/r;}
83inline complex<double> operator/(const double& r,const complex<float>& c) {return r/complex<double>(c);}
84
85inline complex<float> operator+(const complex<float>& c,int r) {return c+(float)r;}
86inline complex<float> operator+(int r,const complex<float>& c) {return (float)r+c;}
87inline complex<double> operator+(const complex<float>& c,const double& r) {return complex<double>(c)+r;}
88inline complex<double> operator+(const double& r,const complex<float>& c) {return r+complex<double>(c);}
89
90inline complex<float> operator-(const complex<float>& c,int r) {return c-(float)r;}
91inline complex<float> operator-(int r,const complex<float>& c) {return (float)r-c;}
92inline complex<double> operator-(const complex<float>& c,const double& r) {return complex<double>(c)-r;}
93inline complex<double> operator-(const double& r,const complex<float>& c) {return r-complex<double>(c);}
94
95// complex<double> with int, float
96inline complex<double> operator*(const complex<double>& c,int r) {return c*(double)r;}
97inline complex<double> operator*(int r,const complex<double>& c) {return (double)r*c;}
98inline complex<double> operator*(const complex<double>& c,const float& r) {return c*(double)r;}
99inline complex<double> operator*(const float& r,const complex<double>& c) {return (double)r*c;}
100
101inline complex<double> operator/(const complex<double>& c,int r) {return c/(double)r;}
102inline complex<double> operator/(int r,const complex<double>& c) {return (double)r/c;}
103inline complex<double> operator/(const complex<double>& c,const float& r) {return c/(double)r;}
104inline complex<double> operator/(const float& r,const complex<double>& c) {return (double)r/c;}
105
106inline complex<double> operator+(const complex<double>& c,int r) {return c+(double)r;}
107inline complex<double> operator+(int r,const complex<double>& c) {return (double)r+c;}
108inline complex<double> operator+(const complex<double>& c,const float& r) {return c+(double)r;}
109inline complex<double> operator+(const float& r,const complex<double>& c) {return (double)r+c;}
110
111inline complex<double> operator-(const complex<double>& c,int r) {return c-(double)r;}
112inline complex<double> operator-(int r,const complex<double>& c) {return (double)r-c;}
113inline complex<double> operator-(const complex<double>& c,const float& r) {return c-(double)r;}
114inline complex<double> operator-(const float& r,const complex<double>& c) {return (double)r-c;}
115} // namespace SimTK
116#endif
117
118namespace SimTK {
119
120template <class R> class conjugate; // Only defined for float, double
121template <> class conjugate<float>;
122template <> class conjugate<double>;
123
124// This is an adaptor for number types which negates the apparent values. A
125// negator<N> has exactly the same internal representation as a number
126// type N, but it is to be interpreted has having the negative of the value
127// it would have if interpreted as an N. This permits negation to be done
128// by reinterpretation rather than computation. A full set of arithmetic operators
129// are provided involving negator<N>'s and N's. Sometimes we can save an op or
130// two this way. For example negator<N>*negator<N> can be performed as an N*N
131// since the negations cancel, and we saved two floating point negations.
132template <class N> class negator; // Only defined for numbers
133
134/*
135 * This class is specialized for all 9 combinations of built-in real types
136 * and contains a typedef for the appropriate "widened" real, complex, or
137 * conjugate type for use when R1 & R2 appear in an operation together.
138 */
139template <class R1, class R2> struct Wider {/* Only defined for built-ins. */};
140template <> struct Wider<float,float> {
141 typedef float WReal;
142 typedef complex<float> WCplx;
144};
145template <> struct Wider<float,double> {
146 typedef double WReal;
147 typedef complex<double> WCplx;
149};
150template <> struct Wider<double,float> {
151 typedef double WReal;
152 typedef complex<double> WCplx;
154};
155template <> struct Wider<double,double> {
156 typedef double WReal;
157 typedef complex<double> WCplx;
159};
160
161
178template <class R> class conjugate {/*Only defined for float, double*/};
179
181// Specialization for conjugate<float> //
183
184template <> class conjugate<float> {
185public:
187 #ifndef NDEBUG
188 re = negIm = std::numeric_limits<float>::quiet_NaN();
189 #endif
190 }
191 // default copy constructor, copy assignment, destructor
192
194 conjugate(const float& real, const float& imag) { re = real; negIm = imag; }
195 conjugate(const float& real, int i) { re = real; negIm = float(i); }
196 conjugate(int r, const float& imag) { re = float(r); negIm = imag; }
197 conjugate(int r, int i) { re = float(r); negIm = float(i); }
198
200 conjugate(const float& real) { re = real; negIm = 0.f; }
201 conjugate(int r) { re = float(r); negIm = 0.f; }
202
203 // No implicit conversions from double because precision will be lost. Some
204 // definitions must be deferred until conjugate<double> is defined below.
205 inline explicit conjugate(const conjugate<double>& cd);
206
207 explicit conjugate(const double& rd)
208 { re = float(rd); negIm = 0.f; }
209
210 // Conversions from complex are always explicit. Note that the value
211 // represented by the conjugate must be identical to that represented by
212 // the complex, which means we must negate the imaginary part.
213 explicit conjugate(const complex<float>& x)
214 { re = x.real(); negIm = -x.imag(); }
215 explicit conjugate(const complex<double>& x)
216 { re = float(x.real()); negIm = float(-x.imag()); }
217
220 operator complex<float>() const
221 { return complex<float>(re,-negIm); }
222
223 // Can't defer here by casting to negator<conjugate> -- this must act
224 // like a built-in. But ... we can use this as a chance to convert
225 // to complex and save one negation.
226 complex<float> operator-() const { return complex<float>(-re,negIm); }
227
228 // Useless.
229 const conjugate& operator+() const { return *this; }
230
231 // Computed assignment operators. We don't depend on implicit conversions
232 // from reals to conjugates here because we can save a few flops by handling
233 // the reals explicitly. Note that we only provide operators for implicitly
234 // convertible precisions, though, which in this case means only floats.
235 conjugate& operator=(const float& r)
236 { re = r; negIm = 0.f; return *this; }
237 conjugate& operator+=(const float& r)
238 { re += r; return *this; }
239 conjugate& operator-=(const float& r)
240 { re -= r; return *this; }
241 conjugate& operator*=(const float& r)
242 { re *= r; negIm *= r; return *this; }
243 conjugate& operator/=(const float& r)
244 { re /= r; negIm /= r; return *this; }
245
247 { re += c.re; negIm += c.negIm; return *this; }
249 { re -= c.re; negIm -= c.negIm; return *this; }
250
251 conjugate& operator=(const complex<float>& c)
252 { re = c.real(); negIm = -c.imag(); return *this; }
253 conjugate& operator+=(const complex<float>& c)
254 { re += c.real(); negIm -= c.imag(); return *this; }
255 conjugate& operator-=(const complex<float>& c)
256 { re -= c.real(); negIm += c.imag(); return *this; }
257
258 // It is pleasant to note that we can self-multiply by either a complex or
259 // a conjugate (leaving a conjugate result) in six flops which is the same
260 // cost as an ordinary complex multiply:
261 // cplx=cplx*cplx: (a+bi)(r+si) = (ar-bs)+(as+br)i
262 // conj=conj*conj: (a-bi)(r-si) = (ar-bs)-(as+br)i
263 // conj=conj*cplx: (a-bi)(r+si) = (ar+bs)-(br-as)i
265 const float r=(re*c.re - negIm*c.negIm);
266 negIm=(re*c.negIm + negIm*c.re); re=r; return *this;
267 }
268 conjugate& operator*=(const complex<float>& t) {
269 const float r=(re*t.real() + negIm*t.imag());
270 negIm=(negIm*t.real() - re*t.imag()); re=r; return *this;
271 }
272
273 // Complex divide is messy and slow anyway so we'll convert to complex and back here,
274 // making use of the fact that for complex c and d, c/d=conj(conj(c)/conj(d)).
276 const complex<float> t = conj()/d.conj();
277 re = t.real(); negIm = t.imag(); // conjugating!
278 return *this;
279 }
280 conjugate& operator/=(const complex<float>& d) {
281 const complex<float> t = conj()/std::conj(d);
282 re = t.real(); negIm = t.imag(); // conjugating!
283 return *this;
284 }
285
286 const float& real() const { return re; }
287 float& real() { return re; }
288
289 const negator<float>& imag() const { return reinterpret_cast<const negator<float>&>(negIm); }
290 negator<float>& imag() { return reinterpret_cast<negator<float>&>(negIm); }
291
292 const complex<float>& conj() const { return reinterpret_cast<const complex<float>&>(*this); }
293 complex<float>& conj() { return reinterpret_cast<complex<float>&>(*this); }
294
295 // Special conjugate methods of use primarily in operator implementations.
296 const float& negImag() const { return negIm; }
297 float& negImag() { return negIm; }
298 bool isReal() const { return negIm==0.f; }
299
300private:
301 float re; // The value represented here is re - negIm*i.
302 float negIm;
303};
304
305
306
307
309// Specialization for conjugate<double> //
311
312template <> class conjugate<double> {
313public:
315 #ifndef NDEBUG
316 re = negIm = std::numeric_limits<double>::quiet_NaN();
317 #endif
318 }
319 // default copy constructor, copy assignment, destructor
320
322 conjugate(const double& real, const double& imag) { re = real; negIm = imag; }
323 conjugate(const double& real, int i) { re = real; negIm = double(i); }
324 conjugate(int r, const double& imag) { re = double(r); negIm = imag; }
325 conjugate(int r, int i) { re = double(r); negIm = double(i); }
326
328 conjugate(const double& real) { re = real; negIm = 0.; }
329 conjugate(int r) { re = double(r); negIm = 0.; }
330
331 // Implicit conversions from float are allowed since
332 // there is no loss in going to double.
334 { re = double(cf.real()); negIm = double(cf.negImag()); }
335 conjugate(const float& rf)
336 { re = double(rf); negIm = 0.; }
337
338 // Conversions from complex are always explicit. Note that the value
339 // represented by the conjugate must be identical to that represented by
340 // the complex, which means we must negate the imaginary part.
341 explicit conjugate(const complex<float>& x)
342 { re = double(x.real()); negIm = double(-x.imag()); }
343 explicit conjugate(const complex<double>& x)
344 { re = x.real(); negIm = -x.imag(); }
345
348 operator complex<double>() const
349 { return complex<double>(re,-negIm); }
350
351 // Can't defer here by casting to negator<conjugate> -- this must act
352 // like a built-in. But ... we can use this as a chance to convert
353 // to complex and save one negation.
354 complex<double> operator-() const { return complex<double>(-re,negIm); }
355
356 // Useless.
357 const conjugate& operator+() const { return *this; }
358
359 // Computed assignment operators. We don't depend on implicit conversions
360 // from reals to conjugates here because we can save a few flops by handling
361 // the reals explicitly. Note that we only provide operators for implicitly
362 // convertible precisions, though, which in this case means floats and doubles.
363 conjugate& operator=(const double& r)
364 { re = r; negIm = 0.; return *this; }
365 conjugate& operator+=(const double& r)
366 { re += r; return *this; }
367 conjugate& operator-=(const double& r)
368 { re -= r; return *this; }
369 conjugate& operator*=(const double& r)
370 { re *= r; negIm *= r; return *this; }
371 conjugate& operator/=(const double& r)
372 { re /= r; negIm /= r; return *this; }
373
374 conjugate& operator=(const float& r)
375 { re = r; negIm = 0.; return *this; }
376 conjugate& operator+=(const float& r)
377 { re += r; return *this; }
378 conjugate& operator-=(const float& r)
379 { re -= r; return *this; }
380 conjugate& operator*=(const float& r)
381 { re *= r; negIm *= r; return *this; }
382 conjugate& operator/=(const float& r)
383 { re /= r; negIm /= r; return *this; }
384
385 // Disambiguate int to be a double.
386 conjugate& operator =(int i) {*this =(double)i; return *this;}
387 conjugate& operator+=(int i) {*this+=(double)i; return *this;}
388 conjugate& operator-=(int i) {*this-=(double)i; return *this;}
389 conjugate& operator*=(int i) {*this*=(double)i; return *this;}
390 conjugate& operator/=(int i) {*this/=(double)i; return *this;}
391
393 { re += c.re; negIm += c.negIm; return *this; }
395 { re -= c.re; negIm -= c.negIm; return *this; }
396
398 { re += c.real(); negIm += c.negImag(); return *this; }
400 { re -= c.real(); negIm -= c.negImag(); return *this; }
401
402 conjugate& operator=(const complex<double>& c)
403 { re = c.real(); negIm = -c.imag(); return *this; }
404 conjugate& operator+=(const complex<double>& c)
405 { re += c.real(); negIm -= c.imag(); return *this; }
406 conjugate& operator-=(const complex<double>& c)
407 { re -= c.real(); negIm += c.imag(); return *this; }
408
409 conjugate& operator=(const complex<float>& c)
410 { re = c.real(); negIm = -c.imag(); return *this; }
411 conjugate& operator+=(const complex<float>& c)
412 { re += c.real(); negIm -= c.imag(); return *this; }
413 conjugate& operator-=(const complex<float>& c)
414 { re -= c.real(); negIm += c.imag(); return *this; }
415
416 // It is pleasant to note that we can self-multiply by either a complex or
417 // a conjugate (leaving a conjugate result) in six flops which is the same
418 // cost as an ordinary complex multiply:
419 // cplx=cplx*cplx: (a+bi)(r+si) = (ar-bs)+(as+br)i
420 // conj=conj*conj: (a-bi)(r-si) = (ar-bs)-(as+br)i
421 // conj=conj*cplx: (a-bi)(r+si) = (ar+bs)-(br-as)i
423 const double r=(re*c.re - negIm*c.negIm);
424 negIm=(re*c.negIm + negIm*c.re); re=r; return *this;
425 }
426 conjugate& operator*=(const complex<double>& t) {
427 const double r=(re*t.real() + negIm*t.imag());
428 negIm=(negIm*t.real() - re*t.imag()); re=r; return *this;
429 }
430
432 conjugate& operator*=(const complex<float>& c) { return operator*=(complex<double>(c)); }
433
434 // Complex divide is messy and slow anyway so we'll convert to complex and back here,
435 // making use of the fact that for complex c and d, c/d=conj(conj(c)/conj(d)).
437 const complex<double> t = conj()/d.conj();
438 re = t.real(); negIm = t.imag(); // conjugating!
439 return *this;
440 }
441 conjugate& operator/=(const complex<double>& d) {
442 const complex<double> t = conj()/std::conj(d);
443 re = t.real(); negIm = t.imag(); // conjugating!
444 return *this;
445 }
446
448 conjugate& operator/=(const complex<float>& c) { return operator/=(complex<double>(c)); }
449
450 const double& real() const { return re; }
451 double& real() { return re; }
452
453 const negator<double>& imag() const { return reinterpret_cast<const negator<double>&>(negIm); }
454 negator<double>& imag() { return reinterpret_cast<negator<double>&>(negIm); }
455
456 const complex<double>& conj() const { return reinterpret_cast<const complex<double>&>(*this); }
457 complex<double>& conj() { return reinterpret_cast<complex<double>&>(*this); }
458
459 // Special conjugate methods of use primarily in operator implementations.
460 const double& negImag() const { return negIm; }
461 double& negImag() { return negIm; }
462 bool isReal() const { return negIm==0.; }
463
464private:
465 double re; // The value represented here is re - negIm*i.
466 double negIm;
467};
468
469
470
471// These definitions had to be deferred until all the specializations have been declared.
473 re = float(cd.real()); negIm = float(cd.negImag());
474}
475
476// Global functions real(),imag(), conj(), abs(), and norm() are overloaded here
477// for efficiency (e.g., abs(c)==abs(conj(c))). The others still work through
478// the implicit conversion from conjugate<T> to complex<T>, which costs
479// one negation. The non-overloaded functions defined for complex are
480// arg(), which returns a real, and a set of functions which return complex:
481// polar,cos,cosh,exp,log,log10,pow,sin,sinh,sqrt,tan,tanh.
482inline const float& real(const conjugate<float>& c) { return c.real(); }
483inline const negator<float>& imag(const conjugate<float>& c) { return c.imag(); }
484inline const complex<float>& conj(const conjugate<float>& c) { return c.conj(); }
485inline float abs (const conjugate<float>& c) { return std::abs(c.conj()); }
486inline float norm(const conjugate<float>& c) { return std::norm(c.conj()); }
487
488inline const double& real(const conjugate<double>& c) { return c.real(); }
489inline const negator<double>& imag(const conjugate<double>& c) { return c.imag(); }
490inline const complex<double>& conj(const conjugate<double>& c) { return c.conj(); }
491inline double abs (const conjugate<double>& c) { return std::abs(c.conj()); }
492inline double norm(const conjugate<double>& c) { return std::norm(c.conj()); }
493
494
495
496
497
498
499// Binary operators with conjugate as one of the operands, and the other any
500// numerical type (real, complex, conjugate) but NOT a negator type. Each operator
501// will silently work with operands of mixed precision, widening the result as
502// necessary. We try to return complex rather than conjugate whenever possible.
503
504template <class R, class CHAR, class TRAITS> inline std::basic_istream<CHAR,TRAITS>&
505operator>>(std::basic_istream<CHAR,TRAITS>& is, conjugate<R>& c) {
506 complex<R> z; is >> z; c=z;
507 return is;
508}
509template <class R, class CHAR, class TRAITS> inline std::basic_ostream<CHAR,TRAITS>&
510operator<<(std::basic_ostream<CHAR,TRAITS>& os, const conjugate<R>& c) {
511 return os << complex<R>(c);
512}
513
514// Operators involving only conjugate and complex can be templatized reliably, as can
515// operators which do not mix precision. But we have to deal explicitly with mixes
516// of conjugate<R> and some other real type S, because the 'class S' template
517// argument can match anything and create ambiguities.
518
519// conjugate<R> with float, double. With 'float' we can be sure that R
520// is the right width for the return value. 'double' is trickier and we have
521// to use the Wider<R,...> helper class to give us the right return type.
522
523// Commutative ops need be done only once: +, *, ==, and != is defined in terms of ==.
524
525// conjugate = conjugate + real
526template <class R> inline conjugate<R> operator+(const conjugate<R>& a, const float& b)
527 { return conjugate<R>(a) += b; }
528template <class R> inline typename Wider<R,double>::WConj operator+(const conjugate<R>& a, const double& b)
529 { return typename Wider<R,double>::WConj(a) += b; }
530
531// conjugate = real + conjugate
532template <class R> inline conjugate<R> operator+(const float& a, const conjugate<R>& b) {return b+a;}
533template <class R> inline typename Wider<R,double>::WConj operator+(const double& a, const conjugate<R>& b) {return b+a;}
534
535// conjugate = conjugate * real
536template <class R> inline conjugate<R> operator*(const conjugate<R>& a, const float& b)
537 { return conjugate<R>(a) *= b; }
538template <class R> inline typename Wider<R,double>::WConj operator*(const conjugate<R>& a, const double& b)
539 { return typename Wider<R,double>::WConj(a) *= b; }
540
541// conjugate = real * conjugate
542template <class R> inline conjugate<R> operator*(const float& a, const conjugate<R>& b) {return b*a;}
543template <class R> inline typename Wider<R,double>::WConj operator*(const double& a, const conjugate<R>& b) {return b*a;}
544
545// bool = conjugate==real
546template <class R> inline bool operator==(const conjugate<R>& a, const float& b)
547 { return a.isReal() && a.real()==b; }
548template <class R> inline bool operator==(const conjugate<R>& a, const double& b)
549 { return a.isReal() && a.real()==b; }
550
551// bool = real==conjugate, bool = conjugate!=real, bool = real!=conjugate
552template <class R> inline bool operator==(const float& a, const conjugate<R>& b) {return b==a;}
553template <class R> inline bool operator==(const double& a, const conjugate<R>& b) {return b==a;}
554template <class R> inline bool operator!=(const conjugate<R>& a, const float& b) {return !(a==b);}
555template <class R> inline bool operator!=(const conjugate<R>& a, const double& b) {return !(a==b);}
556template <class R> inline bool operator!=(const float& a, const conjugate<R>& b) {return !(a==b);}
557template <class R> inline bool operator!=(const double& a, const conjugate<R>& b) {return !(a==b);}
558
559// Non-commutative ops are a little messier.
560
561// conjugate = conjugate - real
562template <class R> inline conjugate<R> operator-(const conjugate<R>& a, const float& b)
563 { return conjugate<R>(a) -= b; }
564template <class R> inline typename Wider<R,double>::WConj operator-(const conjugate<R>& a, const double& b)
565 { return typename Wider<R,double>::WConj(a) -= b; }
566
567// complex = real - conjugate
568// This is nice because -conjugate.imag() is free.
569template <class R> inline complex<R> operator-(const float& a, const conjugate<R>& b)
570 { return complex<R>(a-b.real(), -b.imag()); }
571template <class R> inline typename Wider<R,double>::WCplx operator-(const double& a, const conjugate<R>& b)
572 { return typename Wider<R,double>::WCplx(a-b.real(), -b.imag()); }
573
574// conjugate = conjugate / real
575template <class R> inline conjugate<R> operator/(const conjugate<R>& a, const float& b)
576 { return conjugate<R>(a) /= b; }
577template <class R> inline typename Wider<R,double>::WConj operator/(const conjugate<R>& a, const double& b)
578 { return typename Wider<R,double>::WConj(a) /= b; }
579
580// complex = real / conjugate
581// Division by complex is tricky and slow anyway so we'll just convert to complex
582// at the cost of one negation.
583template <class R> inline complex<R> operator/(const float& a, const conjugate<R>& b)
584 { return (R)a/complex<R>(b); }
585template <class R> inline typename Wider<R,double>::WCplx operator/(const double& a, const conjugate<R>& b)
586 { return (typename Wider<R,double>::WReal)a/(typename Wider<R,double>::WCplx(b)); }
587
588
589// That's it for (conjugate, real) combinations. Now we need to do all the (conjugate, conjugate) and
590// (conjugate, complex) combinations which are safer to templatize fully. There are many more opportunities
591// to return complex rather than conjugate here. Keep in mind here that for a conjugate number c=a-bi,
592// a=c.real() and b=-c.imag() are available for free. conjugate provides negImag() to return b directly.
593// Below we'll capitalize components that should be accessed with negImag().
594
595// conjugate = conjugate + conjugate: (a-Bi)+(r-Si) = (a+r)-(B+S)i
596template <class R, class S> inline typename Wider<R,S>::WConj
598 return typename Wider<R,S>::WConj(a.real()+r.real(),
599 a.negImag()+r.negImag());
600}
601
602// complex = conjugate + complex = complex + conjugate: (a-Bi)+(r+si) = (a+r)+(s-B)i
603template <class R, class S> inline typename Wider<R,S>::WCplx
604operator+(const conjugate<R>& a, const complex<S>& r) {
605 return typename Wider<R,S>::WCplx(a.real()+r.real(),
606 r.imag()-a.negImag());
607}
608template <class R, class S> inline typename Wider<R,S>::WCplx
609operator+(const complex<R>& a, const conjugate<S>& r) { return r+a; }
610
611// complex = conjugate - conjugate: (a-Bi)-(r-Si) = (a-r)+(S-B)i
612template <class R, class S> inline typename Wider<R,S>::WCplx
614 return typename Wider<R,S>::WCplx(a.real()-r.real(),
615 r.negImag()-a.negImag());
616}
617
618// Neg<complex> = conjugate - complex: (a-Bi)-(r+si) = (a-r)+(-B-s)i = -[(r-a)+(B+s)i]
619template <class R, class S> inline negator<typename Wider<R,S>::WCplx>
620operator-(const conjugate<R>& a, const complex<S>& r) {
622 (typename Wider<R,S>::WCplx(r.real()-a.real(),
623 a.negImag()+r.imag()));
624}
625
626// complex = complex - conjugate: (a+bi)-(r-Si) = (a-r)+(b+S)i
627template <class R, class S> inline typename Wider<R,S>::WCplx
628operator-(const complex<R>& a, const conjugate<S>& r) {
629 return typename Wider<R,S>::WCplx(a.real()-r.real(),
630 a.imag()+r.negImag());
631}
632
633// We can multiply by either a complex or a conjugate (leaving a complex
634// or negated complex result) in six flops which is the same cost as an
635// ordinary complex multiply.
636// (cplx=cplx*cplx: (a+bi)(r+si) = (ar-bs)+(as+br)i)
637
638// Neg<cplx>=conj*conj: (a-Bi)(r-Si) = -[(BS-ar)+(aS+Br)i]
639template <class R, class S> inline negator<typename Wider<R,S>::WCplx>
642 (typename Wider<R,S>::WCplx(a.negImag()*r.negImag() - a.real()*r.real(),
643 a.real()*r.negImag() + a.negImag()*r.real()));
644}
645
646// cplx=conj*cplx: (a-Bi)(r+si) = (ar+Bs)+(as-Br)i
647template <class R, class S> inline typename Wider<R,S>::WCplx
648operator*(const conjugate<R>& a, const complex<S>& r) {
649 return typename Wider<R,S>::WCplx(a.real()*r.real() + a.negImag()*r.imag(),
650 a.real()*r.imag() - a.negImag()*r.real());
651}
652
653template <class R, class S> inline typename Wider<R,S>::WCplx
654operator*(const complex<R>& a, const conjugate<S>& r)
655 { return r*a; }
656
657// If there's a negator on the complex number, move it to the conjugate
658// one which will change to complex with just one flop; then this
659// is an ordinary complex mutiply.
660template <class R, class S> inline typename Wider<R,S>::WCplx
661operator*(const negator< complex<R> >& a, const conjugate<S>& r)
662 { return (-a)*(-r); } // -a is free here
663template <class R, class S> inline typename Wider<R,S>::WCplx
664operator*(const conjugate<R>& a, const negator< complex<S> >& r)
665 { return (-a)*(-r); } // -r is free here
666
667// Division is tricky and there is little to gain by trying to exploit
668// the conjugate class, so we'll just convert to complex. (Remember that
669// conj() is free for Conjugates.)
670
671template <class R, class S> inline typename Wider<R,S>::WCplx
673 return std::conj(a.conj()/r.conj());
674}
675
676template <class R, class S> inline typename Wider<R,S>::WCplx
677operator/(const conjugate<R>& a, const complex<S>& r) {
678 return std::conj(a.conj()/std::conj(r));
679}
680
681template <class R, class S> inline typename Wider<R,S>::WCplx
682operator/(const complex<R>& a, const conjugate<S>& r) {
683 return std::conj(std::conj(a)/r.conj());
684}
685
686
687template <class R, class S> inline bool
689 return a.real() == r.real() && a.negImag() == r.negImag();
690}
691
692template <class R, class S> inline bool
693operator==(const conjugate<R>& a, const complex<S>& r) {
694 return a.real() == r.real() && -a.negImag() == r.imag();
695}
696
697template <class R, class S> inline bool
698operator==(const complex<R>& a, const conjugate<S>& r) {return r==a;}
699
700template <class R, class S> inline bool
701operator!=(const conjugate<R>& a, const conjugate<S>& r) {return !(a==r);}
702
703template <class R, class S> inline bool
704operator!=(const conjugate<R>& a, const complex<S>& r) {return !(a==r);}
705
706template <class R, class S> inline bool
707operator!=(const complex<R>& a, const conjugate<S>& r) {return !(a==r);}
708
709
710} // namespace SimTK
711
712#endif //SimTK_SIMMATRIX_CONJUGATE_H_
conjugate & operator+=(const conjugate< double > &c)
Definition conjugate.h:392
conjugate & operator*=(const conjugate< float > &c)
Definition conjugate.h:431
conjugate & operator/=(const complex< double > &d)
Definition conjugate.h:441
conjugate & operator+=(const double &r)
Definition conjugate.h:365
conjugate & operator=(const complex< float > &c)
Definition conjugate.h:409
conjugate & operator=(const float &r)
Definition conjugate.h:374
conjugate & operator*=(const conjugate< double > &c)
Definition conjugate.h:422
conjugate & operator/=(const conjugate< float > &c)
Definition conjugate.h:447
conjugate & operator/=(const complex< float > &c)
Definition conjugate.h:448
conjugate(int r, const double &imag)
Definition conjugate.h:324
conjugate & operator+=(const complex< float > &c)
Definition conjugate.h:411
const double & real() const
Definition conjugate.h:450
conjugate(int r, int i)
Definition conjugate.h:325
double & real()
Definition conjugate.h:451
conjugate(const complex< float > &x)
Definition conjugate.h:341
conjugate & operator=(const double &r)
Definition conjugate.h:363
double & negImag()
Definition conjugate.h:461
negator< double > & imag()
Definition conjugate.h:454
conjugate & operator+=(const float &r)
Definition conjugate.h:376
const negator< double > & imag() const
Definition conjugate.h:453
const complex< double > & conj() const
Definition conjugate.h:456
conjugate & operator*=(const double &r)
Definition conjugate.h:369
conjugate(const complex< double > &x)
Definition conjugate.h:343
conjugate & operator-=(const double &r)
Definition conjugate.h:367
complex< double > & conj()
Definition conjugate.h:457
conjugate(const double &real, const double &imag)
Construction from reals. Note that the numeric result is (real-imag*i).
Definition conjugate.h:322
conjugate & operator/=(const conjugate< double > &d)
Definition conjugate.h:436
conjugate & operator/=(const float &r)
Definition conjugate.h:382
bool isReal() const
Definition conjugate.h:462
conjugate & operator-=(int i)
Definition conjugate.h:388
conjugate & operator*=(const float &r)
Definition conjugate.h:380
conjugate & operator-=(const complex< double > &c)
Definition conjugate.h:406
conjugate & operator/=(const double &r)
Definition conjugate.h:371
conjugate & operator-=(const complex< float > &c)
Definition conjugate.h:413
conjugate(const double &real, int i)
Definition conjugate.h:323
conjugate & operator*=(const complex< double > &t)
Definition conjugate.h:426
conjugate & operator*=(int i)
Definition conjugate.h:389
conjugate & operator-=(const conjugate< double > &c)
Definition conjugate.h:394
conjugate & operator+=(int i)
Definition conjugate.h:387
conjugate(int r)
Definition conjugate.h:329
const conjugate & operator+() const
Definition conjugate.h:357
conjugate & operator*=(const complex< float > &c)
Definition conjugate.h:432
conjugate & operator+=(const conjugate< float > &c)
Definition conjugate.h:397
const double & negImag() const
Definition conjugate.h:460
complex< double > operator-() const
Definition conjugate.h:354
conjugate & operator+=(const complex< double > &c)
Definition conjugate.h:404
conjugate()
Definition conjugate.h:314
conjugate & operator=(const complex< double > &c)
Definition conjugate.h:402
conjugate(const double &real)
Implicit conversion from double to conjugate<double>.
Definition conjugate.h:328
conjugate & operator/=(int i)
Definition conjugate.h:390
conjugate(const conjugate< float > &cf)
Definition conjugate.h:333
conjugate & operator-=(const conjugate< float > &c)
Definition conjugate.h:399
conjugate(const float &rf)
Definition conjugate.h:335
conjugate & operator-=(const float &r)
Definition conjugate.h:378
conjugate(int r, int i)
Definition conjugate.h:197
conjugate(const float &real, int i)
Definition conjugate.h:195
conjugate & operator+=(const complex< float > &c)
Definition conjugate.h:253
conjugate & operator-=(const conjugate< float > &c)
Definition conjugate.h:248
conjugate(const float &real, const float &imag)
Construction from reals. Note that the numeric result is (real-imag*i).
Definition conjugate.h:194
float & negImag()
Definition conjugate.h:297
conjugate(int r, const float &imag)
Definition conjugate.h:196
conjugate & operator/=(const complex< float > &d)
Definition conjugate.h:280
const float & negImag() const
Definition conjugate.h:296
conjugate(int r)
Definition conjugate.h:201
const float & real() const
Definition conjugate.h:286
conjugate & operator*=(const float &r)
Definition conjugate.h:241
negator< float > & imag()
Definition conjugate.h:290
complex< float > & conj()
Definition conjugate.h:293
conjugate(const complex< double > &x)
Definition conjugate.h:215
conjugate & operator/=(const float &r)
Definition conjugate.h:243
conjugate & operator+=(const conjugate< float > &c)
Definition conjugate.h:246
conjugate & operator-=(const float &r)
Definition conjugate.h:239
conjugate & operator=(const float &r)
Definition conjugate.h:235
float & real()
Definition conjugate.h:287
bool isReal() const
Definition conjugate.h:298
conjugate(const double &rd)
Definition conjugate.h:207
conjugate(const float &real)
Implicit conversion from float to conjugate<float>.
Definition conjugate.h:200
conjugate & operator*=(const conjugate< float > &c)
Definition conjugate.h:264
const negator< float > & imag() const
Definition conjugate.h:289
conjugate & operator+=(const float &r)
Definition conjugate.h:237
conjugate & operator=(const complex< float > &c)
Definition conjugate.h:251
conjugate(const complex< float > &x)
Definition conjugate.h:213
conjugate & operator/=(const conjugate< float > &d)
Definition conjugate.h:275
conjugate & operator*=(const complex< float > &t)
Definition conjugate.h:268
const complex< float > & conj() const
Definition conjugate.h:292
conjugate()
Definition conjugate.h:186
const conjugate & operator+() const
Definition conjugate.h:229
conjugate & operator-=(const complex< float > &c)
Definition conjugate.h:255
complex< float > operator-() const
Definition conjugate.h:226
SimTK::conjugate<R> should be instantiated only for float, double.
Definition conjugate.h:178
negator<N>, where N is a number type (real, complex, conjugate), is represented in memory identically...
Definition negator.h:75
This is the top-level SimTK namespace into which all SimTK names are placed to avoid collision with o...
Definition Assembler.h:37
RowVectorBase< typename CNT< ELEM >::TAbs > abs(const RowVectorBase< ELEM > &v)
Definition VectorMath.h:120
Matrix_< E > operator*(const MatrixBase< E > &l, const typename CNT< E >::StdNumber &r)
Definition BigMatrix.h:605
Matrix_< E > operator/(const MatrixBase< E > &l, const typename CNT< E >::StdNumber &r)
Definition BigMatrix.h:613
std::basic_istream< CHAR, TRAITS > & operator>>(std::basic_istream< CHAR, TRAITS > &is, conjugate< R > &c)
Definition conjugate.h:505
Matrix_< typename CNT< E1 >::template Result< E2 >::Add > operator+(const MatrixBase< E1 > &l, const MatrixBase< E2 > &r)
Definition BigMatrix.h:568
std::ostream & operator<<(std::ostream &o, const ContactForce &f)
Definition CompliantContactSubsystem.h:387
Matrix_< typename CNT< E1 >::template Result< E2 >::Sub > operator-(const MatrixBase< E1 > &l, const MatrixBase< E2 > &r)
Definition BigMatrix.h:584
const float & real(const conjugate< float > &c)
Definition conjugate.h:482
float norm(const conjugate< float > &c)
Definition conjugate.h:486
const complex< float > & conj(const conjugate< float > &c)
Definition conjugate.h:484
bool operator==(const PhiMatrix &p1, const PhiMatrix &p2)
Definition SpatialAlgebra.h:791
bool operator!=(const L &left, const R &right)
Definition SimTKcommon/include/SimTKcommon/internal/common.h:641
const negator< float > & imag(const conjugate< float > &c)
Definition conjugate.h:483
double WReal
Definition conjugate.h:156
complex< double > WCplx
Definition conjugate.h:157
conjugate< double > WConj
Definition conjugate.h:158
conjugate< double > WConj
Definition conjugate.h:153
double WReal
Definition conjugate.h:151
complex< double > WCplx
Definition conjugate.h:152
double WReal
Definition conjugate.h:146
conjugate< double > WConj
Definition conjugate.h:148
complex< double > WCplx
Definition conjugate.h:147
complex< float > WCplx
Definition conjugate.h:142
float WReal
Definition conjugate.h:141
conjugate< float > WConj
Definition conjugate.h:143
Definition conjugate.h:139