Re: complex - direct access to members?

Re: complex - direct access to members?

[Gabriel Dos Reis]

[-- Warning: decoded text below may be mangled, UTF-8 assumed --]
[-- Attachment #1: Type: text/plain, Size: 784 bytes --]

>>>>> «Mike», Mike Stump <mrs@wrs.com> wrote:

> You can lobby ANSI to `clarify' this point, if you want...  See
> comp.std.c++...  Best bet.  If you want, you can also to this in user
> code, with apparently no loss in performace.


[snip proxy class idiom]


> Simple, easy, reasonably direct, all brought to you by the power of
> C++.  Yes, it is a bit of a pain, but I think you can do it.  The
> above obviously isn't complete, and you might find other ways to do
> it, the above is only a starting point, but enough to prove the idea.

This scheme works fine in an application programmer code; but as
a libstdc++ contributor it is helpless in order to make the library
efficient *and* conforming. I would have expected an optimization
support from the compiler. 

-- Gaby

Re: complex - direct access to members?

[Mike Stump]

You can lobby ANSI to `clarify' this point, if you want...  See
comp.std.c++...  Best bet.  If you want, you can also to this in user
code, with apparently no loss in performace.

// You want:

//	c.imag() = 12;

// You can only do:

//	c = C(c.real(), 12)

// First, we define helper classes to adapt the syntax:

template <class T>
class Hi {
  T& val;
public:
  Hi(T& t) : val(t) { }
  template <class Y>
  Hi& operator = (const Y& y) {
    val = T(val.real (), y);
    return *this;
  }
  template <class Y>
  operator Y() {
    return val.imag();
  }
};

template <class T>
class Hr {
  T& val;
public:
  Hr(T& t) : val(t) { }
  template <class Y>
  Hr& operator = (const Y& y) {
    val = T(y, val.imag());
    return *this;
  }
  template <class Y>
  operator Y() {
    return val.real();
  }
};

#include <complex.h>

template <class T>
class MyComplex {
  complex<T> i;
public:
  const Hi<complex<T> > imag() { return Hi<complex<T> >(i); }
};

void foo1();
void foo2();
void foo3();

int main() {
  complex<float> i;
  MyComplex<float> j;

  foo1();

  i = complex<float>(i.real(),12);

  foo2();

  j.imag() = 12;

  foo3();
}


Simple, easy, reasonably direct, all brought to you by the power of
C++.  Yes, it is a bit of a pain, but I think you can do it.  The
above obviously isn't complete, and you might find other ways to do
it, the above is only a starting point, but enough to prove the idea.

RE: complex - direct access to members?

[Joe Buck]

> What about an implementation that stores polar values R,Theta?  The
> following code sets real(a), and generates a lot assembly code (egcs g++
> 1.1).  In particular it appears that the inlined a.imag=imag(a)
> assignment does not get optimized away.

This optimization failure is a more general problem in egcs, in that the
compiler decides too soon that it has a store to memory.  However, as you
suggest, in principle if this flaw were corrected, we could generate as
good code from

template<class T>
inline void setReal(complex<T>& z, T newReal)
{
	z = complex<T>(newReal, imag(z));
}

as we could if imag(z) were changed to return a reference.  Certainly
Neal could use this setReal, though at present he will get worse code.

If the speed penalty is too much, then this setReal can be made a
friend of complex<T> in his local copy and implemented as

	z.re = newReal;

retaining the other definition for portability.

Re: complex - direct access to members?

[Michael Hayes]

Craig Burley writes:

 > In case you didn't notice somebody else's comment, I'll reiterate it:
 > there apparently is usefulness to an alternate representation of
 > complex numbers in some kinds of computation, specifically, polar
 > notation.  In this case, there'd still be two double values, but for
 > angle and length ("rho" is the angle IIRC, can't recall the geek name
 > for length).

r or rho is usually used for length; theta for angle.  (geek greek!)

Michael.

Re: complex - direct access to members?

[nbecker]

[-- Attachment #1: Type: text/plain, Size: 492 bytes --]

>>>>> "Gabriel" == Gabriel Dos Reis <Gabriel.Dos-Reis@dptmaths.ens-cachan.fr> writes:

    Gabriel> If you have any performance degradation due to the complex type as
    Gabriel> currently implemented it would be helpfull if you could precise your
    Gabriel> architecture, provide testcase and the optimization you use. Then, we
    Gabriel> could see what we can do.

OK, here's an example.

I'm writing algorithms that interface via STL iterators.  For example, 
here is an FIR filter:


[-- Attachment #2: ComplexAdapt.H --]
[-- Type: text/x-c++, Size: 2392 bytes --]

#ifndef _ComplexAdapt_H
#define _ComplexAdapt_H

#include <complex>
#include <cstddef>

template<class T, class ComplexIterator>
class ComplexToRealAdapt {
private:
  ComplexIterator I;
public:
  typedef T value_type;
  typedef ComplexToRealAdapt self;
  typedef ptrdiff_t difference_type;
  typedef size_t size_type;

  ComplexToRealAdapt (const ComplexIterator& C) :
    I (C)
  {}

  value_type operator*() const {
    return real (*I);
  }

  self& operator++ () {
    I++;
    return *this;
  }

  self operator++ (int) {
    self tmp = *this;
    ++*this;
    return tmp;
  }

  self& operator-- () {
    I--;
    return *this;
  }

  self operator-- (int) {
    self tmp = *this;
    --*this;
    return tmp;
  }

  self& operator+=(difference_type n) {
    I += n;
    return *this;
  }

  self& operator-=(difference_type n) { return *this += -n; }

  self operator+(difference_type n) const {
    self tmp = *this;
    return tmp += n;
  }

  value_type operator[](difference_type n) const { return *(*this + n); }

  bool operator==(const self& x) const { return I == x.I; }

  bool operator!=(const self& x) const { return !(*this == x); }

  difference_type operator- (const self& x) const {
      return I - x.I;
  }

};

template<class T, class ComplexIterator>
class ComplexToImagAdapt {
private:
  ComplexIterator I;
public:
  typedef T value_type;
  typedef ComplexToImagAdapt self;
  typedef ptrdiff_t difference_type;
  typedef size_t size_type;

  ComplexToImagAdapt (const ComplexIterator& C) :
    I (C)
  {}

  value_type operator*() const {
    return imag (*I);
  }

  self& operator++ () {
    I++;
    return *this;
  }

  self operator++ (int) {
    self tmp = *this;
    ++*this;
    return tmp;
  }

  self& operator-- () {
    I--;
    return *this;
  }

  self operator-- (int) {
    self tmp = *this;
    --*this;
    return tmp;
  }

  self& operator+=(difference_type n) {
    I += n;
    return *this;
  }

  self& operator-=(difference_type n) { return *this += -n; }

  self operator+(difference_type n) const {
    self tmp = *this;
    return tmp += n;
  }

  value_type operator[](difference_type n) const { return *(*this + n); }

  bool operator==(const self& x) const { return I == x.I; }

  bool operator!=(const self& x) const { return !(*this == x); }

  difference_type operator- (const self& x) const {
      return I - x.I;
  }

};

#endif


[-- Attachment #3: QPSKDemod.H --]
[-- Type: text/x-c++, Size: 1279 bytes --]

#ifndef QPSKDemod_H
#define QPSKDemod_H

#include <complex>
#include "FIRRcvPulse.H"
#include "PulseShape.H"
#include "ComplexAdapt.H"

template<class T>
class QPSKDemod {

protected:

  FIRRcvPulse<T> IrcvFP;
  FIRRcvPulse<T> QrcvFP;
  const int SPS;
  const int decim;

public:

  QPSKDemod (int _SPS, int _SI, const PulseShape<T>& _P, int _decim = 0);

  template<class InputIterator, class OutputIterator>
  void Compute (const InputIterator& instart, const InputIterator& inend, const OutputIterator& Iout, const OutputIterator& Qout);

};

template<class T>
QPSKDemod<T>::QPSKDemod (int _SPS, int _SI, const PulseShape<T>& _P, int _decim) :
  IrcvFP (_P, _SPS, _SI, (_decim == 0) ? _SPS : _decim, 0), 
  QrcvFP (_P, _SPS, _SI, (_decim == 0) ? _SPS : _decim, 0),
  SPS (_SPS),
  decim (_decim)
{}


template<class T>
template<class InputIterator, class OutputIterator>
void QPSKDemod<T>::Compute (const InputIterator& instart, const InputIterator& inend, const OutputIterator& Iout, const OutputIterator& Qout) {

  for (InputIterator in = instart; in != inend; in += decim) {
    ComplexToRealAdapt<T, InputIterator> Rin (in);
    ComplexToImagAdapt<T, InputIterator> Iin (in);

    IrcvFP.Go (Rin+0, Rin+decim, Iout);
    QrcvFP.Go (Iin+0, Iin+decim, Qout);
  }
}

#endif

[-- Attachment #4: SimpleFIR.H --]
[-- Type: text/x-c++, Size: 9870 bytes --]

#ifndef _SimpleFIR_H
#define _SimpleFIR_H

#include <vector>
#include "RingBuffer.H"

// FIR assumes the buffering of the input is done externally.

template<class InputType, class OutputType = InputType, class CoefType = InputType> class FIRbase {

protected:

  vector<CoefType> coef;

public:

  template<class InputIterator>
  FIRbase (InputIterator firsttap, InputIterator lasttap) :
    coef (firsttap, lasttap) {}

  FIRbase (const vector<CoefType>& _coef) :
    coef (_coef) {}

  const CoefType* Coef() const { return coef; }

  int Len() const { return coef.size(); }

  void Print (ostream& os) const;

};

template<class InputType, class OutputType = InputType, class CoefType = InputType> class FIR2base :
  public FIRbase<InputType, OutputType, CoefType> {

protected:

  RingBuffer<InputType> input;

  template<class InputIterator>
  void ShiftIn (const InputIterator& in);

  template<class InputIterator>
  void ShiftIn (const InputIterator& instart, const InputIterator& inend);

public:

  template<class InputIterator>
  FIR2base (InputIterator firsttap, InputIterator lasttap, int BufSize = lasttap - firsttap + 1) :
    FIRbase<InputType, OutputType, CoefType> (firsttap, lasttap),
    input (BufSize) {}

  FIR2base (const vector<CoefType>& _coef, int BufSize = 0) : 
    FIRbase<InputType, OutputType, CoefType> (_coef), 
    input ((BufSize == 0) ? _coef.size() : BufSize)
  {}

};


template<class InputType, class OutputType = InputType, class CoefType = InputType> class SimpleFIR :
  public FIRbase<InputType, OutputType, CoefType> {

public:

  template<class InputIterator>
  SimpleFIR (InputIterator firsttap, InputIterator lasttap) :
    FIRbase (firsttap, lasttap) {}

  SimpleFIR (istream& is);

  template<class InputIterator, class OutputIterator>
  void Go (const InputIterator& instart, const InputIterator& inend, const OutputIterator& outbuff);

};

// FIR2 provides a ring buffer for the input internally

template<class InputType, class OutputType = InputType, class CoefType = InputType> class SimpleFIR2 : 
  public FIR2base<InputType, OutputType, CoefType> {

public:

  template<class InputIterator>
  SimpleFIR2 (InputIterator firsttap, InputIterator lasttap) :
    FIR2base<InputType, OutputType, CoefType> (firsttap, lasttap) {}

  SimpleFIR2 (istream& is);

  template<class InputIterator, class OutputIterator>
  void Go (const InputIterator& instart, const InputIterator& inend, const OutputIterator& outbuff);
  
  OutputType Go (InputType in);

};

template<class InputType, class OutputType = InputType, class CoefType = InputType> class InterpFIR : public FIRbase<InputType, OutputType, CoefType> {

protected:
    
  int interp;  
    
public:

  template<class InputIterator>
  InterpFIR (InputIterator firsttap, InputIterator lasttap, int _interp) :
    FIRbase (firsttap, lasttap),
    interp (_interp)
    {}

  InterpFIR (vector<CoefType> V, int _interp) :
    FIRbase (V),
    interp (_interp)
    {}

  InterpFIR (istream& is, int _interp) :
    FIRbase (V),
    interp (_interp)
    {}

  template<class InputIterator, class OutputIterator>
  void Go (const InputIterator& instart, const InputIterator& inend, const OutputIterator& outbuff);

};

template<class InputType, class OutputType = InputType, class CoefType = InputType> class InterpFIR2 :
  public FIR2base<InputType, OutputType, CoefType> {

protected:

  int interp;

public:

  template<class InputIterator>
  InterpFIR2 (InputIterator firsttap, InputIterator lasttap, int _interp) :
    FIR2base<InputType, OutputType, CoefType> (firsttap, lasttap),
    interp (_interp) {}

  InterpFIR2 (istream& is, int _interp);

  template<class InputIterator, class OutputIterator>
  void Go (const InputIterator& instart, const InputIterator& inend, const OutputIterator& outbuff);

};


template<class InputType, class OutputType = InputType, class CoefType = InputType> class DecimFIR :
public FIRbase<InputType, OutputType, CoefType> {

protected:

  int decim;
  int phase;

public:

  template<class InputIterator>
  DecimFIR (InputIterator firsttap, InputIterator lasttap, int _decim, int _phase) :
    FIRbase (firsttap, lasttap),
    decim (_decim),
    phase (_phase)
    {}

  DecimFIR (istream& is, int _decim, int _phase) :
    FIRbase (V),
    decim (_decim),
    phase (_phase)
    {}

  template<class InputIterator, class OutputIterator>
  void Go (const InputIterator& instart, const InputIterator& inend, const OutputIterator& outbuff);


};

template<class InputType, class OutputType = InputType, class CoefType = InputType>
class DecimFIR2 : 
  public FIR2base<InputType, OutputType, CoefType> {

private:

  int decim;
  int phase;

public:

  template<class InputIterator>
  DecimFIR2 (InputIterator firsttap, InputIterator lasttap, int _decim, int _phase, int BufSize = lasttap - firsttap + 1 + _decim) :
    FIR2base<InputType, OutputType, CoefType> (firsttap, lasttap, BufSize),
    decim (_decim),
    phase (_phase) {}

  DecimFIR2 (const vector<CoefType>& _coef, int _decim, int _phase, int BufSize = 0) :
    FIR2base<InputType, OutputType, CoefType> (_coef, (BufSize == 0) ? _coef.size() + _decim : BufSize),
    decim (_decim),
    phase (_phase)
  {}

  //  DecimFIR2 (istream& is, int _decim, int _phase);

  int Decim() const { return decim; }

  template<class InputIterator, class OutputIterator>
  void Go (const InputIterator& instart, const InputIterator& inend, const OutputIterator& outbuff);

};

template<class InputType, class OutputType, class CoefType> 
SimpleFIR<InputType, OutputType, CoefType>::SimpleFIR (istream& is)
{
  while (is) {
    CoefType x;
    is >> x;
    if (is)
      coef->push_back (x);
  }
}

template<class InputType, class OutputType, class CoefType> 
template<class InputIterator, class OutputIterator> 
void SimpleFIR<InputType, OutputType, CoefType>::Go (const InputIterator& instart, const InputIterator& inend, const OutputIterator& outstart) {
  OutputIterator O = outstart;

  for (InputIterator in = instart; in != inend; in++) {

    OutputType sum = 0;
    for (int s = 0; s < Len(); s++)
      sum += coef[s] * in[-s];

    *O++ = sum;
  }
}

#include <iostream.h>

template <class InputType, class OutputType, class CoefType> 
void FIRbase<InputType, OutputType, CoefType>::Print( ostream& os ) const {
  for( int i = 0; i < coef.size(); i++ )
    os << i << '\t' << coef[ i ] << endl;

}

template<class InputType, class OutputType, class CoefType> 
template<class InputIterator>
void FIR2base<InputType, OutputType, CoefType>::ShiftIn (const InputIterator& in) {
  input.Shift ();
  input[0] = in[0];
}

template<class InputType, class OutputType, class CoefType> 
template<class InputIterator>
void FIR2base<InputType, OutputType, CoefType>::ShiftIn (const InputIterator& instart, const InputIterator& inend) {
  for (InputIterator in = instart; in != inend; in++) {
    input.Shift ();
    input[0] = *in;
  }
}

template<class InputType, class OutputType, class CoefType> 
template<class InputIterator, class OutputIterator> 
void SimpleFIR2<InputType, OutputType, CoefType>::Go (const InputIterator& instart, const InputIterator& inend, const OutputIterator& outbuff) {
  OutputIterator O = outbuff;

  for (InputIterator in = instart; in != inend; in++) {
    ShiftIn (in);

    OutputType sum = 0;
    for (int s = 0; s < Len(); s++)
      sum += coef[s] * input[-s];

    *O++ = sum;
  }
}

template<class InputType, class OutputType, class CoefType> 
OutputType SimpleFIR2<InputType, OutputType, CoefType>::Go (InputType in) {
  ShiftIn (&in);

  OutputType sum = 0;
  for (int s = 0; s < Len(); s++)
    sum += coef[s] * input[-s];

  return sum;
}

template<class InputType, class OutputType, class CoefType> 
template<class InputIterator, class OutputIterator> 
void InterpFIR<InputType, OutputType, CoefType>::Go (const InputIterator& instart, const InputIterator& inend, const OutputIterator& outbuff) {
  OutputIterator O = outbuff;

  for (InputIterator in = instart; in != inend; in++) {
    for (int phase = 0; phase < interp; phase++) {
      OutputType sum = 0;
      //      for (int i = 0; i < Len(); i++)
      int i = 0;
      while ((i * interp) + phase < Len()) {
	sum += coef[(i * interp) + phase] * in[-i];
	i++;
      }
      *O++ = sum;
    }
  }
}

template<class InputType, class OutputType, class CoefType> 
template<class InputIterator, class OutputIterator> 
void InterpFIR2<InputType, OutputType, CoefType>::Go (const InputIterator& instart, const InputIterator& inend, const OutputIterator& outbuff) {

  for (InputIterator in = instart; in != inend; in++) {
    ShiftIn (in);

    OutputIterator O = outbuff;
    for (int phase = 0; phase < interp; phase++) {
      OutputType sum = 0;
      int i = 0;
      while ((i * interp) + phase < Len()) {
	sum += coef[(i * interp) + phase] * input[-i];
	i++;
      }
      *O++ = sum;
    }
  }
}

template<class InputType, class OutputType, class CoefType> 
template<class InputIterator, class OutputIterator> 
void DecimFIR<InputType, OutputType, CoefType>::Go (const InputIterator& instart, const InputIterator& inend, const OutputIterator& outbuff) {
  OutputIterator O = outbuff;

  for (InputIterator in = instart; in != inend; in += decim) {
    OutputType sum = 0;
    for (int i = 0; i < Len(); i++)
      sum += coef[i] * in[-(i+phase)];
    *O++ = sum;
  }
}

template<class InputType, class OutputType, class CoefType> 
template<class InputIterator, class OutputIterator> 
void DecimFIR2<InputType, OutputType, CoefType>::Go (const InputIterator& instart, const InputIterator& inend, const OutputIterator& outbuff) {
  OutputIterator O = outbuff;

  for (InputIterator in = instart; in != inend; in += decim) {
    ShiftIn (in, in+decim);
    OutputType sum = 0;
    for (int i = 0; i < Len(); i++)
      sum += coef[i] * input[-(i+phase)];
    *O++ = sum;
  }
}

#endif

RE: complex - direct access to members?

[Jan Reimers]

What about an implementation that stores polar values R,Theta?  The
following code sets real(a), and generates a lot assembly code (egcs g++
1.1).  In particular it appears that the inlined a.imag=imag(a)
assignment does not get optimized away.

#include <complex>

complex<double> a(1.0,1.0);

void foo(double myreal)
{
  a=complex<double>(myreal,imag(a)); //Reset real(a).
}

Assembly after -O2 compilation:

_foo__Fd:
	pushl %ebp
	movl %esp,%ebp
	subl $16,%esp
	fldl _a+8
	fldl 8(%ebp)
	fstpl -16(%ebp)
	movl -16(%ebp),%eax
	fstpl -8(%ebp)
	movl %eax,_a
	movl $_a,%edx
	movl -12(%ebp),%eax
	movl %eax,4(%edx)
	movl -8(%ebp),%eax
	movl %eax,8(%edx)
	movl -4(%ebp),%eax
	movl %eax,12(%edx)
	movl %ebp,%esp
	popl %ebp
	ret



JR

> ----------
> From: 	nbecker@fred.net[SMTP:nbecker@fred.net]
> Sent: 	Tuesday, January 12, 1999 5:52 AM
> To: 	egcs@cygnus.com
> Subject: 	complex - direct access to members?
> 
> As a heavy user of c++ for numeric computation, I think it would be a
> great improvement if we made a simple modification to complex, so that
> 
> the member functions real() and imag() returned a reference.  This
> would allow modification of the values directly.  I really can't see
> any argument against it.  I don't think this violates any
> encapsulation.  I can't imagine any implementation of complex<double>
> in which there is not somewhere a double value for the real and for
> the imaginary parts.
> 

Re: complex - direct access to members?

[Gabriel Dos Reis]

>>>>> Jan Reimers <janr@molienergy.bc.ca> wrote:

> What about an implementation that stores polar values R,Theta?  The
> following code sets real(a), and generates a lot assembly code (egcs g++
> 1.1).  In particular it appears that the inlined a.imag=imag(a)
> assignment does not get optimized away.

[...]

Thanks. Your snippet reveals flaw in the EGCS optimization
machinery. But it seems to me that EGCS does not bet on a large number
crunching community so (basic) memory reference optimizations don't
appear to be on the hot TODO list :-(. Maybe some day they will.

-- Gaby

Re: complex - direct access to members?

[Joe Buck]

> As a heavy user of c++ for numeric computation, I think it would be a
> great improvement if we made a simple modification to complex, so that 
> the member functions real() and imag() returned a reference.  This
> would allow modification of the values directly.  I really can't see
> any argument against it.

If the standard were not final, this would be a good argument.
The argument against it now is "this is not what the standard says".

> I don't think this violates any
> encapsulation.  I can't imagine any implementation of complex<double>
> in which there is not somewhere a double value for the real and for
> the imaginary parts.

If we are going to use the class name specified in ISO C++, we should
implement the class specified by ISO C++.

If you want a class that permits direct access, you should name it
something other than complex<T>.

Re: complex - direct access to members?

[Craig Burley]

>I can't imagine any implementation of complex<double>
>in which there is not somewhere a double value for the real and for
>the imaginary parts.

In case you didn't notice somebody else's comment, I'll reiterate it:
there apparently is usefulness to an alternate representation of
complex numbers in some kinds of computation, specifically, polar
notation.  In this case, there'd still be two double values, but for
angle and length ("rho" is the angle IIRC, can't recall the geek name
for length).

Dunno whether the *design* of that class mandates real/imaginary,
though.  In some ways, Fortran effectively does this in its COMPLEX-
type support, by allowing aliasing with floating-point types and
mandating that the aliasing accesses the real and imaginary parts,
in that order (in memory).  But I've lost track of whether that
particular mandate is continued beyond just COMPLEX in Fortran 90
(i.e. whether COMPLEX of an arbitrary, vendor-supplied KIND must
meet that mandate).  My *guess* is that it does not, so a vendor
could provide COMPLEX types implemented in polar notation, corresponding
to each of the real/imaginary types it supports.

        tq vm, (burley)

complex - direct access to members?

[nbecker]

As a heavy user of c++ for numeric computation, I think it would be a
great improvement if we made a simple modification to complex, so that 
the member functions real() and imag() returned a reference.  This
would allow modification of the values directly.  I really can't see
any argument against it.  I don't think this violates any
encapsulation.  I can't imagine any implementation of complex<double>
in which there is not somewhere a double value for the real and for
the imaginary parts.

Re: complex - direct access to members?

[Gabriel Dos Reis]

>>>>>   <nbecker@fred.net> wrote:

 > As a heavy user of c++ for numeric computation, I think it would be a
 > great improvement if we made a simple modification to complex, so that 
 > the member functions real() and imag() returned a reference.  This
 > would allow modification of the values directly.  I really can't see
 > any argument against it.  I don't think this violates any
 > encapsulation. 

This violates the ISO C++ definition of std::complex<>. Furthermore
you can change any the real and imaginary part of any complex<> by using
operator= ; also it not obvious at all whether your proposal is
optimal (for all types and on the platforms currently supported).

 > ... I can't imagine any implementation of complex<double>
 > in which there is not somewhere a double value for the real and for
 > the imaginary parts.

If you have any performance degradation due to the complex type as
currently implemented it would be helpfull if you could precise your
architecture, provide testcase and the optimization you use. Then, we
could see what we can do.

-- Gaby