Simbody  3.6
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 
64 using 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
73 namespace SimTK {
74 // complex<float> with int, double
75 inline complex<float> operator*(const complex<float>& c,int r) {return c*(float)r;}
76 inline complex<float> operator*(int r,const complex<float>& c) {return (float)r*c;}
77 inline complex<double> operator*(const complex<float>& c,const double& r) {return complex<double>(c)*r;}
78 inline complex<double> operator*(const double& r,const complex<float>& c) {return r*complex<double>(c);}
79 
80 inline complex<float> operator/(const complex<float>& c,int r) {return c/(float)r;}
81 inline complex<float> operator/(int r,const complex<float>& c) {return (float)r/c;}
82 inline complex<double> operator/(const complex<float>& c,const double& r) {return complex<double>(c)/r;}
83 inline complex<double> operator/(const double& r,const complex<float>& c) {return r/complex<double>(c);}
84 
85 inline complex<float> operator+(const complex<float>& c,int r) {return c+(float)r;}
86 inline complex<float> operator+(int r,const complex<float>& c) {return (float)r+c;}
87 inline complex<double> operator+(const complex<float>& c,const double& r) {return complex<double>(c)+r;}
88 inline complex<double> operator+(const double& r,const complex<float>& c) {return r+complex<double>(c);}
89 
90 inline complex<float> operator-(const complex<float>& c,int r) {return c-(float)r;}
91 inline complex<float> operator-(int r,const complex<float>& c) {return (float)r-c;}
92 inline complex<double> operator-(const complex<float>& c,const double& r) {return complex<double>(c)-r;}
93 inline complex<double> operator-(const double& r,const complex<float>& c) {return r-complex<double>(c);}
94 
95 // complex<double> with int, float
96 inline complex<double> operator*(const complex<double>& c,int r) {return c*(double)r;}
97 inline complex<double> operator*(int r,const complex<double>& c) {return (double)r*c;}
98 inline complex<double> operator*(const complex<double>& c,const float& r) {return c*(double)r;}
99 inline complex<double> operator*(const float& r,const complex<double>& c) {return (double)r*c;}
100 
101 inline complex<double> operator/(const complex<double>& c,int r) {return c/(double)r;}
102 inline complex<double> operator/(int r,const complex<double>& c) {return (double)r/c;}
103 inline complex<double> operator/(const complex<double>& c,const float& r) {return c/(double)r;}
104 inline complex<double> operator/(const float& r,const complex<double>& c) {return (double)r/c;}
105 
106 inline complex<double> operator+(const complex<double>& c,int r) {return c+(double)r;}
107 inline complex<double> operator+(int r,const complex<double>& c) {return (double)r+c;}
108 inline complex<double> operator+(const complex<double>& c,const float& r) {return c+(double)r;}
109 inline complex<double> operator+(const float& r,const complex<double>& c) {return (double)r+c;}
110 
111 inline complex<double> operator-(const complex<double>& c,int r) {return c-(double)r;}
112 inline complex<double> operator-(int r,const complex<double>& c) {return (double)r-c;}
113 inline complex<double> operator-(const complex<double>& c,const float& r) {return c-(double)r;}
114 inline complex<double> operator-(const float& r,const complex<double>& c) {return (double)r-c;}
115 } // namespace SimTK
116 #endif
117 
118 namespace SimTK {
119 
120 template <class R> class conjugate; // Only defined for float, double
121 template <> class conjugate<float>;
122 template <> 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.
132 template <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  */
139 template <class R1, class R2> struct Wider {/* Only defined for built-ins. */};
140 template <> struct Wider<float,float> {
141  typedef float WReal;
142  typedef complex<float> WCplx;
144 };
145 template <> struct Wider<float,double> {
146  typedef double WReal;
147  typedef complex<double> WCplx;
149 };
150 template <> struct Wider<double,float> {
151  typedef double WReal;
152  typedef complex<double> WCplx;
154 };
155 template <> struct Wider<double,double> {
156  typedef double WReal;
157  typedef complex<double> WCplx;
159 };
160 
161 
178 template <class R> class conjugate {/*Only defined for float, double*/};
179 
181 // Specialization for conjugate<float> //
183 
184 template <> class conjugate<float> {
185 public:
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 
300 private:
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 
312 template <> class conjugate<double> {
313 public:
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 
464 private:
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.
482 inline const float& real(const conjugate<float>& c) { return c.real(); }
483 inline const negator<float>& imag(const conjugate<float>& c) { return c.imag(); }
484 inline const complex<float>& conj(const conjugate<float>& c) { return c.conj(); }
485 inline float abs (const conjugate<float>& c) { return std::abs(c.conj()); }
486 inline float norm(const conjugate<float>& c) { return std::norm(c.conj()); }
487 
488 inline const double& real(const conjugate<double>& c) { return c.real(); }
489 inline const negator<double>& imag(const conjugate<double>& c) { return c.imag(); }
490 inline const complex<double>& conj(const conjugate<double>& c) { return c.conj(); }
491 inline double abs (const conjugate<double>& c) { return std::abs(c.conj()); }
492 inline 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 
504 template <class R, class CHAR, class TRAITS> inline std::basic_istream<CHAR,TRAITS>&
505 operator>>(std::basic_istream<CHAR,TRAITS>& is, conjugate<R>& c) {
506  complex<R> z; is >> z; c=z;
507  return is;
508 }
509 template <class R, class CHAR, class TRAITS> inline std::basic_ostream<CHAR,TRAITS>&
510 operator<<(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
526 template <class R> inline conjugate<R> operator+(const conjugate<R>& a, const float& b)
527  { return conjugate<R>(a) += b; }
528 template <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
532 template <class R> inline conjugate<R> operator+(const float& a, const conjugate<R>& b) {return b+a;}
533 template <class R> inline typename Wider<R,double>::WConj operator+(const double& a, const conjugate<R>& b) {return b+a;}
534 
535 // conjugate = conjugate * real
536 template <class R> inline conjugate<R> operator*(const conjugate<R>& a, const float& b)
537  { return conjugate<R>(a) *= b; }
538 template <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
542 template <class R> inline conjugate<R> operator*(const float& a, const conjugate<R>& b) {return b*a;}
543 template <class R> inline typename Wider<R,double>::WConj operator*(const double& a, const conjugate<R>& b) {return b*a;}
544 
545 // bool = conjugate==real
546 template <class R> inline bool operator==(const conjugate<R>& a, const float& b)
547  { return a.isReal() && a.real()==b; }
548 template <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
552 template <class R> inline bool operator==(const float& a, const conjugate<R>& b) {return b==a;}
553 template <class R> inline bool operator==(const double& a, const conjugate<R>& b) {return b==a;}
554 template <class R> inline bool operator!=(const conjugate<R>& a, const float& b) {return !(a==b);}
555 template <class R> inline bool operator!=(const conjugate<R>& a, const double& b) {return !(a==b);}
556 template <class R> inline bool operator!=(const float& a, const conjugate<R>& b) {return !(a==b);}
557 template <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
562 template <class R> inline conjugate<R> operator-(const conjugate<R>& a, const float& b)
563  { return conjugate<R>(a) -= b; }
564 template <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.
569 template <class R> inline complex<R> operator-(const float& a, const conjugate<R>& b)
570  { return complex<R>(a-b.real(), -b.imag()); }
571 template <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
575 template <class R> inline conjugate<R> operator/(const conjugate<R>& a, const float& b)
576  { return conjugate<R>(a) /= b; }
577 template <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.
583 template <class R> inline complex<R> operator/(const float& a, const conjugate<R>& b)
584  { return (R)a/complex<R>(b); }
585 template <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
596 template <class R, class S> inline typename Wider<R,S>::WConj
597 operator+(const conjugate<R>& a, const conjugate<S>& r) {
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
603 template <class R, class S> inline typename Wider<R,S>::WCplx
604 operator+(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 }
608 template <class R, class S> inline typename Wider<R,S>::WCplx
609 operator+(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
612 template <class R, class S> inline typename Wider<R,S>::WCplx
613 operator-(const conjugate<R>& a, const conjugate<S>& r) {
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]
619 template <class R, class S> inline negator<typename Wider<R,S>::WCplx>
620 operator-(const conjugate<R>& a, const complex<S>& r) {
621  return negator<typename Wider<R,S>::WCplx>::recast
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
627 template <class R, class S> inline typename Wider<R,S>::WCplx
628 operator-(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]
639 template <class R, class S> inline negator<typename Wider<R,S>::WCplx>
640 operator*(const conjugate<R>& a, const conjugate<S>& r) {
641  return negator<typename Wider<R,S>::WCplx>::recast
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
647 template <class R, class S> inline typename Wider<R,S>::WCplx
648 operator*(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 
653 template <class R, class S> inline typename Wider<R,S>::WCplx
654 operator*(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.
660 template <class R, class S> inline typename Wider<R,S>::WCplx
661 operator*(const negator< complex<R> >& a, const conjugate<S>& r)
662  { return (-a)*(-r); } // -a is free here
663 template <class R, class S> inline typename Wider<R,S>::WCplx
664 operator*(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 
671 template <class R, class S> inline typename Wider<R,S>::WCplx
672 operator/(const conjugate<R>& a, const conjugate<S>& r) {
673  return std::conj(a.conj()/r.conj());
674 }
675 
676 template <class R, class S> inline typename Wider<R,S>::WCplx
677 operator/(const conjugate<R>& a, const complex<S>& r) {
678  return std::conj(a.conj()/std::conj(r));
679 }
680 
681 template <class R, class S> inline typename Wider<R,S>::WCplx
682 operator/(const complex<R>& a, const conjugate<S>& r) {
683  return std::conj(std::conj(a)/r.conj());
684 }
685 
686 
687 template <class R, class S> inline bool
688 operator==(const conjugate<R>& a, const conjugate<S>& r) {
689  return a.real() == r.real() && a.negImag() == r.negImag();
690 }
691 
692 template <class R, class S> inline bool
693 operator==(const conjugate<R>& a, const complex<S>& r) {
694  return a.real() == r.real() && -a.negImag() == r.imag();
695 }
696 
697 template <class R, class S> inline bool
698 operator==(const complex<R>& a, const conjugate<S>& r) {return r==a;}
699 
700 template <class R, class S> inline bool
701 operator!=(const conjugate<R>& a, const conjugate<S>& r) {return !(a==r);}
702 
703 template <class R, class S> inline bool
704 operator!=(const conjugate<R>& a, const complex<S>& r) {return !(a==r);}
705 
706 template <class R, class S> inline bool
707 operator!=(const complex<R>& a, const conjugate<S>& r) {return !(a==r);}
708 
709 
710 } // namespace SimTK
711 
712 #endif //SimTK_SIMMATRIX_CONJUGATE_H_
bool isReal() const
Definition: conjugate.h:298
const float & negImag() const
Definition: conjugate.h:296
Matrix_< E > operator/(const MatrixBase< E > &l, const typename CNT< E >::StdNumber &r)
Definition: BigMatrix.h:613
conjugate & operator*=(const complex< float > &t)
Definition: conjugate.h:268
double & real()
Definition: conjugate.h:451
conjugate< double > WConj
Definition: conjugate.h:153
conjugate & operator=(const double &r)
Definition: conjugate.h:363
conjugate & operator*=(const double &r)
Definition: conjugate.h:369
conjugate(const float &real)
Implicit conversion from float to conjugate<float>.
Definition: conjugate.h:200
conjugate(const double &real, int i)
Definition: conjugate.h:323
conjugate & operator-=(const conjugate< double > &c)
Definition: conjugate.h:394
complex< float > & conj()
Definition: conjugate.h:293
negator< float > & imag()
Definition: conjugate.h:290
complex< double > operator-() const
Definition: conjugate.h:354
conjugate & operator-=(const float &r)
Definition: conjugate.h:378
negator< double > & imag()
Definition: conjugate.h:454
conjugate(const complex< float > &x)
Definition: conjugate.h:341
double WReal
Definition: conjugate.h:151
conjugate & operator-=(const double &r)
Definition: conjugate.h:367
const float & real() const
Definition: conjugate.h:286
complex< double > & conj()
Definition: conjugate.h:457
This is the top-level SimTK namespace into which all SimTK names are placed to avoid collision with o...
Definition: Assembler.h:37
SimTK::conjugate<R> should be instantiated only for float, double.
Definition: String.h:45
conjugate & operator*=(const float &r)
Definition: conjugate.h:380
conjugate & operator+=(const double &r)
Definition: conjugate.h:365
conjugate & operator+=(const conjugate< float > &c)
Definition: conjugate.h:397
Matrix_< typename CNT< E1 >::template Result< E2 >::Sub > operator-(const MatrixBase< E1 > &l, const MatrixBase< E2 > &r)
Definition: BigMatrix.h:584
const double & negImag() const
Definition: conjugate.h:460
float WReal
Definition: conjugate.h:141
Definition: conjugate.h:184
conjugate & operator*=(int i)
Definition: conjugate.h:389
Definition: conjugate.h:312
const negator< double > & imag() const
Definition: conjugate.h:453
conjugate & operator=(const float &r)
Definition: conjugate.h:235
conjugate< double > WConj
Definition: conjugate.h:148
std::basic_istream< CHAR, TRAITS > & operator>>(std::basic_istream< CHAR, TRAITS > &is, conjugate< R > &c)
Definition: conjugate.h:505
complex< double > WCplx
Definition: conjugate.h:152
conjugate(int r, int i)
Definition: conjugate.h:325
negator<N>, where N is a number type (real, complex, conjugate), is represented in memory identically...
Definition: String.h:44
conjugate & operator/=(const conjugate< double > &d)
Definition: conjugate.h:436
conjugate & operator-=(const complex< float > &c)
Definition: conjugate.h:255
conjugate(const float &rf)
Definition: conjugate.h:335
conjugate< float > WConj
Definition: conjugate.h:143
complex< float > WCplx
Definition: conjugate.h:142
float & real()
Definition: conjugate.h:287
bool operator==(const PhiMatrix &p1, const PhiMatrix &p2)
Definition: SpatialAlgebra.h:791
conjugate & operator+=(const complex< float > &c)
Definition: conjugate.h:253
double WReal
Definition: conjugate.h:156
conjugate(const complex< double > &x)
Definition: conjugate.h:343
conjugate & operator=(const float &r)
Definition: conjugate.h:374
conjugate & operator+=(const conjugate< float > &c)
Definition: conjugate.h:246
conjugate & operator*=(const conjugate< float > &c)
Definition: conjugate.h:431
conjugate & operator*=(const conjugate< double > &c)
Definition: conjugate.h:422
conjugate(int r, const double &imag)
Definition: conjugate.h:324
conjugate(const conjugate< float > &cf)
Definition: conjugate.h:333
const conjugate & operator+() const
Definition: conjugate.h:229
conjugate & operator/=(const complex< double > &d)
Definition: conjugate.h:441
conjugate(int r, const float &imag)
Definition: conjugate.h:196
double & negImag()
Definition: conjugate.h:461
conjugate(int r)
Definition: conjugate.h:329
Matrix_< E > operator*(const MatrixBase< E > &l, const typename CNT< E >::StdNumber &r)
Definition: BigMatrix.h:605
const negator< float > & imag() const
Definition: conjugate.h:289
complex< double > WCplx
Definition: conjugate.h:157
Definition: conjugate.h:139
conjugate & operator*=(const complex< double > &t)
Definition: conjugate.h:426
conjugate(const complex< float > &x)
Definition: conjugate.h:213
conjugate()
Definition: conjugate.h:186
complex< double > WCplx
Definition: conjugate.h:147
conjugate & operator+=(int i)
Definition: conjugate.h:387
conjugate & operator+=(const float &r)
Definition: conjugate.h:237
conjugate & operator-=(const complex< float > &c)
Definition: conjugate.h:413
conjugate(int r, int i)
Definition: conjugate.h:197
conjugate(int r)
Definition: conjugate.h:201
conjugate(const float &real, const float &imag)
Construction from reals. Note that the numeric result is (real-imag*i).
Definition: conjugate.h:194
float norm(const conjugate< float > &c)
Definition: conjugate.h:486
conjugate & operator-=(const conjugate< float > &c)
Definition: conjugate.h:248
conjugate & operator/=(const float &r)
Definition: conjugate.h:382
const float & real(const conjugate< float > &c)
Definition: conjugate.h:482
RowVectorBase< typename CNT< ELEM >::TAbs > abs(const RowVectorBase< ELEM > &v)
Definition: VectorMath.h:120
double norm(const conjugate< double > &c)
Definition: conjugate.h:492
double WReal
Definition: conjugate.h:146
conjugate & operator/=(const double &r)
Definition: conjugate.h:371
conjugate & operator*=(const complex< float > &c)
Definition: conjugate.h:432
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 complex< double > &c)
Definition: conjugate.h:402
const complex< double > & conj(const conjugate< double > &c)
Definition: conjugate.h:490
conjugate & operator/=(const complex< float > &c)
Definition: conjugate.h:448
conjugate & operator+=(const complex< float > &c)
Definition: conjugate.h:411
const complex< float > & conj() const
Definition: conjugate.h:292
conjugate & operator*=(const conjugate< float > &c)
Definition: conjugate.h:264
conjugate(const double &rd)
Definition: conjugate.h:207
conjugate & operator/=(int i)
Definition: conjugate.h:390
conjugate(const float &real, int i)
Definition: conjugate.h:195
conjugate()
Definition: conjugate.h:314
conjugate & operator=(const complex< float > &c)
Definition: conjugate.h:251
const complex< double > & conj() const
Definition: conjugate.h:456
Matrix_< typename CNT< E1 >::template Result< E2 >::Add > operator+(const MatrixBase< E1 > &l, const MatrixBase< E2 > &r)
Definition: BigMatrix.h:568
bool isReal() const
Definition: conjugate.h:462
const negator< float > & imag(const conjugate< float > &c)
Definition: conjugate.h:483
conjugate(const double &real)
Implicit conversion from double to conjugate<double>.
Definition: conjugate.h:328
float & negImag()
Definition: conjugate.h:297
bool operator!=(const L &left, const R &right)
Definition: SimTKcommon/include/SimTKcommon/internal/common.h:641
const double & real() const
Definition: conjugate.h:450
conjugate< double > WConj
Definition: conjugate.h:158
const conjugate & operator+() const
Definition: conjugate.h:357
conjugate & operator-=(int i)
Definition: conjugate.h:388
const complex< float > & conj(const conjugate< float > &c)
Definition: conjugate.h:484
conjugate & operator-=(const complex< double > &c)
Definition: conjugate.h:406
conjugate & operator-=(const conjugate< float > &c)
Definition: conjugate.h:399
conjugate & operator/=(const conjugate< float > &c)
Definition: conjugate.h:447
conjugate & operator/=(const float &r)
Definition: conjugate.h:243
conjugate & operator=(const complex< float > &c)
Definition: conjugate.h:409
conjugate & operator-=(const float &r)
Definition: conjugate.h:239
complex< float > operator-() const
Definition: conjugate.h:226
conjugate & operator+=(const complex< double > &c)
Definition: conjugate.h:404
conjugate & operator/=(const conjugate< float > &d)
Definition: conjugate.h:275
conjugate(const complex< double > &x)
Definition: conjugate.h:215
conjugate & operator+=(const float &r)
Definition: conjugate.h:376
conjugate & operator/=(const complex< float > &d)
Definition: conjugate.h:280
conjugate & operator+=(const conjugate< double > &c)
Definition: conjugate.h:392
conjugate & operator*=(const float &r)
Definition: conjugate.h:241
double abs(const conjugate< double > &c)
Definition: conjugate.h:491