Complex question

[jacob navia]

hi

I am trying to use the complex data type.

Consider this code, taken from the cclib library:

struct complex csqrt(Cpx z)
{ double r;
r=sqrt(z.re*z.re+z.im*z.im);
r=sqrt(ldexp(r+fabs(z.re),-1));
if(r==0.) z.re=z.im=0.;
else{
if(z.re>=0.){ z.re=r; z.im=ldexp(z.im/r,-1);}
else{ z.re=ldexp(fabs(z.im)/r,-1);
if(z.im>=0.) z.im=r; else z.im= -r;
}
}
return z;
}

Now, this works in lcc-win easily because lcc-win
allows access to the real and imaginary parts of the complex
data type easily, using
z.re (real part)
and
ze.im

for the imaginary part.

I do not see how can be this written in standard C.
Has anyone an idea?

Thanks
--
jacob navia
jacob at jacob point remcomp point fr
logiciels/informatique
http://www.cs.virginia.edu/~lcc-win32

[user923005]

On Jan 24, 2:03*pm, jacob navia <[e-mail]:

hi

I am trying to use the complex data type.

Consider this code, taken from the cclib library:

struct complex csqrt(Cpx z)
{ double r;
* *r=sqrt(z.re*z.re+z.im*z.im);
* *r=sqrt(ldexp(r+fabs(z.re),-1));
* *if(r==0.) z.re=z.im=0.;
* *else{
* * *if(z.re>=0.){ z.re=r; z.im=ldexp(z.im/r,-1);}
* * *else{ z.re=ldexp(fabs(z.im)/r,-1);
* * * *if(z.im>=0.) z.im=r; else z.im= -r;
* * * }
* * }
* *return z;

}

Now, this works in lcc-win easily because lcc-win
allows access to the real and imaginary parts of the complex
data type easily, using
z.re (real part)
and
ze.im

for the imaginary part.

I do not see how can be this written in standard C.
Has anyone an idea?

To get the imaginary part:

#include <complex.h>
double cimag(double complex z);
float cimagf(float complex z);
long double cimagl(long double complex z);

But then again, to get the square root, you would just call the
function:

#include <complex.h>
double complex csqrt(double complex z);
float complex csqrtf(float complex z);
long double complex csqrtl(long double complex z);

So I guess you do not have the header complex.h?
In such a case, maybe this is helpful:
http://www.moshier.net/#Complex_variables

[user923005]

On Jan 24, 2:03*pm, jacob navia <[e-mail]:

hi

I am trying to use the complex data type.

Consider this code, taken from the cclib library:

struct complex csqrt(Cpx z)
{ double r;
* *r=sqrt(z.re*z.re+z.im*z.im);
* *r=sqrt(ldexp(r+fabs(z.re),-1));
* *if(r==0.) z.re=z.im=0.;
* *else{
* * *if(z.re>=0.){ z.re=r; z.im=ldexp(z.im/r,-1);}
* * *else{ z.re=ldexp(fabs(z.im)/r,-1);
* * * *if(z.im>=0.) z.im=r; else z.im= -r;
* * * }
* * }
* *return z;

}

Now, this works in lcc-win easily because lcc-win
allows access to the real and imaginary parts of the complex
data type easily, using
z.re (real part)
and
ze.im

for the imaginary part.

I do not see how can be this written in standard C.
Has anyone an idea?

Just guessing, since I do not have a C99 compiler handy right now, but
I would say:

#include <complex.h>
double complex my_csqrt(double complex z)
{
double r;
r = sqrt(creal(z) * creal(z) + (cimag(z) * cimag(z)) * I;
r = sqrt(ldexp(r + fabs(creal(z)), -1));
if (r == 0.)
z= 0.0 + 0.0 * I;
else {
if (creal(z) >= 0.) {
z = r + ldexp(cimag(z) / r, -1) * I;
} else {
z = ldexp(fabs(cimag(z)) / r, -1) + cimag(z) >= 0 ? r *
I : -r * I;
}
}
return z;
}

[Martin Ambuhl]

jacob navia wrote:

hi

I am trying to use the complex data type.

Consider this code, taken from the cclib library:

struct complex csqrt(Cpx z)
{ double r;
r=sqrt(z.re*z.re+z.im*z.im);
r=sqrt(ldexp(r+fabs(z.re),-1));
if(r==0.) z.re=z.im=0.;
else{
if(z.re>=0.){ z.re=r; z.im=ldexp(z.im/r,-1);}
else{ z.re=ldexp(fabs(z.im)/r,-1);
if(z.im>=0.) z.im=r; else z.im= -r;
}
}
return z;
}

Now, this works in lcc-win easily because lcc-win
allows access to the real and imaginary parts of the complex
data type easily, using
z.re (real part)
and
ze.im

for the imaginary part.

I do not see how can be this written in standard C.
Has anyone an idea?

Your function works just fine in standard C. Check the following:
#include <math.h>
#include <stdio.h>

/* if you don't have <complex.h>, remove or comment out the line below */
#define USE_COMPLEX_H

typedef struct cmplx
{
double re, im;
} Cpx;
/* 'cmplx' chosen for the name of the struct since 'complex' should be
<complex.h*/

Cpx Csqrt
/* renamed because csqrt should be in <complex.h*/
(Cpx z)
{
double r;
r = sqrt(z.re * z.re + z.im * z.im);
r = sqrt(ldexp(r + fabs(z.re), -1));
if (r == 0.)
z.re = z.im = 0.;
else {
if (z.re >= 0.) {
z.re = r;
z.im = ldexp(z.im / r, -1);
}
else {
z.re = ldexp(fabs(z.im) / r, -1);
if (z.im >= 0.)
z.im = r;
else
z.im = -r;
}
}
return z;
}

Cpx Csquare(Cpx z)
{
Cpx r;
r.re = z.re * z.re - z.im * z.im;
r.im = 2 * z.re * z.im;
return r;
}

#if defined(USE_COMPLEX_H)
#include <complex.h>

double complex Csqrt2(double complex z)
{
double r;
r = sqrt(creal(z) * creal(z) + cimag(z) * cimag(z));
r = sqrt(ldexp(r + fabs(creal(z)), -1));
if (r == 0.)
z = 0.;
else {
if (creal(z) >= 0.) {
z = r + ldexp(cimag(z) / r, -1) * I;
}
else {
z = ldexp(fabs(cimag(z)) / r,
-1) + (cimag(z) >= 0. ? r : -r) * I;
}
}
return z;
}
#endif

int main(void)
{
Cpx x, y, z;
#if defined(USE_COMPLEX_H)
double complex cx, cy, cz;
#endif
int i, j;
printf("Using the function Jacob Navia provided,\n"
"with a struct for complex types.\n\n");
for (i = -2; i < 3; i++)
for (j = -2; j < 3; j++) {
x.re = i;
x.im = j;
y = Csqrt(x);
z = Csquare(y);
printf("y = Csqrt(%g+%gi)) = %g+%gi, y*y=%g+%gi\n",
x.re, x.im, y.re, y.im, z.re, z.im);
}

#if defined(USE_COMPLEX_H)
printf("\nUsing the function Jacob Navia provided,\n"
"rewritten to use complex types.\n\n");
for (i = -2; i < 3; i++)
for (j = -2; j < 3; j++) {
cx = i + j * I;
cy = Csqrt2(cx);
cz = cy * cy;
printf("y = Csqrt2(%g+%gi)) = %g+%gi, y*y=%g+%gi\n",
creal(cx), cimag(cx),
creal(cy), cimag(cy),
creal(cz), cimag(cz));
}

#endif

return 0;
}

[Output]

Using the function Jacob Navia provided,
with a struct for complex types.

y = Csqrt(-2+-2i)) = 0.643594+-1.55377i, y*y=-2+-2i
y = Csqrt(-2+-1i)) = 0.343561+-1.45535i, y*y=-2+-1i
y = Csqrt(-2+0i)) = 0+1.41421i, y*y=-2+0i
y = Csqrt(-2+1i)) = 0.343561+1.45535i, y*y=-2+1i
y = Csqrt(-2+2i)) = 0.643594+1.55377i, y*y=-2+2i
y = Csqrt(-1+-2i)) = 0.786151+-1.27202i, y*y=-1+-2i
y = Csqrt(-1+-1i)) = 0.45509+-1.09868i, y*y=-1+-1i
y = Csqrt(-1+0i)) = 0+1i, y*y=-1+0i
y = Csqrt(-1+1i)) = 0.45509+1.09868i, y*y=-1+1i
y = Csqrt(-1+2i)) = 0.786151+1.27202i, y*y=-1+2i
y = Csqrt(0+-2i)) = 1+-1i, y*y=0+-2i
y = Csqrt(0+-1i)) = 0.707107+-0.707107i, y*y=1.5702e-16+-1i
y = Csqrt(0+0i)) = 0+0i, y*y=0+0i
y = Csqrt(0+1i)) = 0.707107+0.707107i, y*y=1.5702e-16+1i
y = Csqrt(0+2i)) = 1+1i, y*y=0+2i
y = Csqrt(1+-2i)) = 1.27202+-0.786151i, y*y=1+-2i
y = Csqrt(1+-1i)) = 1.09868+-0.45509i, y*y=1+-1i
y = Csqrt(1+0i)) = 1+0i, y*y=1+0i
y = Csqrt(1+1i)) = 1.09868+0.45509i, y*y=1+1i
y = Csqrt(1+2i)) = 1.27202+0.786151i, y*y=1+2i
y = Csqrt(2+-2i)) = 1.55377+-0.643594i, y*y=2+-2i
y = Csqrt(2+-1i)) = 1.45535+-0.343561i, y*y=2+-1i
y = Csqrt(2+0i)) = 1.41421+0i, y*y=2+0i
y = Csqrt(2+1i)) = 1.45535+0.343561i, y*y=2+1i
y = Csqrt(2+2i)) = 1.55377+0.643594i, y*y=2+2i

Using the function Jacob Navia provided,
rewritten to use __complex__ types.

y = Csqrt2(-2+-2i)) = 0.643594+-1.55377i, y*y=-2+-2i
y = Csqrt2(-2+-1i)) = 0.343561+-1.45535i, y*y=-2+-1i
y = Csqrt2(-2+0i)) = 0+1.41421i, y*y=-2+0i
y = Csqrt2(-2+1i)) = 0.343561+1.45535i, y*y=-2+1i
y = Csqrt2(-2+2i)) = 0.643594+1.55377i, y*y=-2+2i
y = Csqrt2(-1+-2i)) = 0.786151+-1.27202i, y*y=-1+-2i
y = Csqrt2(-1+-1i)) = 0.45509+-1.09868i, y*y=-1+-1i
y = Csqrt2(-1+0i)) = 0+1i, y*y=-1+0i
y = Csqrt2(-1+1i)) = 0.45509+1.09868i, y*y=-1+1i
y = Csqrt2(-1+2i)) = 0.786151+1.27202i, y*y=-1+2i
y = Csqrt2(0+-2i)) = 1+-1i, y*y=0+-2i
y = Csqrt2(0+-1i)) = 0.707107+-0.707107i, y*y=1.5702e-16+-1i
y = Csqrt2(0+0i)) = 0+0i, y*y=0+0i
y = Csqrt2(0+1i)) = 0.707107+0.707107i, y*y=1.5702e-16+1i
y = Csqrt2(0+2i)) = 1+1i, y*y=0+2i
y = Csqrt2(1+-2i)) = 1.27202+-0.786151i, y*y=1+-2i
y = Csqrt2(1+-1i)) = 1.09868+-0.45509i, y*y=1+-1i
y = Csqrt2(1+0i)) = 1+0i, y*y=1+0i
y = Csqrt2(1+1i)) = 1.09868+0.45509i, y*y=1+1i
y = Csqrt2(1+2i)) = 1.27202+0.786151i, y*y=1+2i
y = Csqrt2(2+-2i)) = 1.55377+-0.643594i, y*y=2+-2i
y = Csqrt2(2+-1i)) = 1.45535+-0.343561i, y*y=2+-1i
y = Csqrt2(2+0i)) = 1.41421+0i, y*y=2+0i
y = Csqrt2(2+1i)) = 1.45535+0.343561i, y*y=2+1i
y = Csqrt2(2+2i)) = 1.55377+0.643594i, y*y=2+2i

[Pierre Asselin]

jacob navia <ja***@nospam.comwrote:

In lcc-win it works since lcc-win defines
#define creal(z) z.re

Nitpick: you need
#define creal(z) ((z).re)

to avoid bad interactions with operators when z is an expression
and creal(z) is part of a larger expression.

(this apart from whether creal() can be implemented as a macro,
and whether the standard should provide direct accessors to the
members of a complex.)
--
pa at panix dot com

[Army1987]

jacob navia wrote:

Army1987 wrote:

>z = creal(z) + I * new_imag_part;
z = new_real_part + I * cimag(z);

Well in the chapter "Awful syntax anyone?" that would
take the first prize...

Or maybe that's too harsh.

fn = (int *(*)(int))foo;

is maybe better

Afwulness is in the eye of the beholder.

I think I will keep my "extensions"

z.im = new_imag_part;

Easy isn't it?

I think that you're allowed to do that, as that couldn't break any strictly
conforming program.

Why do we have creal()z) (a "getter" function) and not setreal(z) (the
"setter" part)?

Yes, maybe that'be useful. (Of course that'd not be setreal(z)...)
--
Army1987 (Replace "NOSPAM" with "email")

[Philip Potter]

Army1987 wrote:

jacob navia wrote:

>I think I will keep my "extensions"

z.im = new_imag_part;

Easy isn't it?

I think that you're allowed to do that, as that couldn't break any strictly
conforming program.

Not breaking strictly conforming programs is only one requirement of a
conforming implementation. Jacob's extension violates a constraint,
according to n1256 6.5.2.3p1. He's only allowed to do it provided he
issues a diagnostic in conforming mode.

>Why do we have creal()z) (a "getter" function) and not setreal(z) (the
"setter" part)?

Yes, maybe that'be useful. (Of course that'd not be setreal(z)...)

I can see the argument for it, though it's purely syntactic since

setreal(&z,real); /* or z = setreal(z,real); */

would be equivalent to

z = real + cimag(z)*I;

In fact I prefer the current version of doing things, because C's syntax
does not lend itself to setter functions (ie z.setreal(real) or
setreal(z,real) would not be possible.)

[jacob navia]

Philip Potter wrote:

Army1987 wrote:

>jacob navia wrote:

>>I think I will keep my "extensions"

z.im = new_imag_part;

Easy isn't it?

I think that you're allowed to do that, as that couldn't break any
strictly
conforming program.

Not breaking strictly conforming programs is only one requirement of a
conforming implementation. Jacob's extension violates a constraint,
according to n1256 6.5.2.3p1. He's only allowed to do it provided he
issues a diagnostic in conforming mode.

>>Why do we have creal()z) (a "getter" function) and not setreal(z) (the
"setter" part)?

Yes, maybe that'be useful. (Of course that'd not be setreal(z)...)

I can see the argument for it, though it's purely syntactic since

setreal(&z,real); /* or z = setreal(z,real); */

would be equivalent to

z = real + cimag(z)*I;

In fact I prefer the current version of doing things, because C's syntax
does not lend itself to setter functions (ie z.setreal(real) or
setreal(z,real) would not be possible.)

But why z.real and z.imaginary would be forbidden?

It is clear that complex are built from float/double/long double PAIRS.
Why couldn't we agree in the name of those elemnts of the pairs
and be done with it?

Writing

z.real = 5.887;

is Much clearer than

z = creal(z) + 5.887*I;

Besides, it requires a lot LESS implementation effort
than to recognize the first thing as its equivalent.

For a compiler is harder to recognize ALL the possible forms of the
first statement as equivalent

z = 5.887*I+creal(z);
z = (6.887-1.0)*I+creal(z);

etc etc etc!

--
jacob navia
jacob at jacob point remcomp point fr
logiciels/informatique
http://www.cs.virginia.edu/~lcc-win32

[Pierre Asselin]

Army1987 <ar******@nospam.itwrote:

Any function can be implemented as a macro, but it ought to be
#define creal(z) (((double _Complex)(z)).re)
as it should have the right type regardless of the type of z.

But then you lose the lvalue-preserving property. It looks like
there is no way to do what Jacob wanted without breaking the
conformance of creal(). Too bad.

Come to think of it, how is it done in Fortran ? The same way
as in C99, I think.
complex z
z= cmplx(new_real_val, imag(z))
z= cmplx(real(z), new_imag_val)

No direct accessors there either. To change one of the components
you have to synthesize a new complex number.

--
pa at panix dot com

[Martin Ambuhl]

user923005 wrote:

It could be implemented (for instance) as:

complex {
double parts[2];
};

or
complex {
double re;
double im;
};

or whatever underneath the covers, since that bit is not specified
clearly.

Of course it might be implemented as the standard specifies:
"Each complex type has the same representation and alignment
requirements as an array type containing exactly two elements of the
corresponding real type; the first element is equal to the real part,
and the second element to the imaginary part, of the complex
number."
When you think about that, then I'm sure you will realize that there are
a number of ways to accomplish what you want.

[user923005]

On Jan 25, 9:01*pm, Martin Ambuhl <[e-mail]:

user923005 wrote:

It could be implemented (for instance) as:

complex {
double parts[2];
};

or
complex {
double re;
double im;
};

or whatever underneath the covers, since that bit is not specified
clearly.

Of course it might be implemented as the standard specifies:
"Each complex type has the same representation and alignment
requirements as an array type containing exactly two elements of the
corresponding real type; the first element is equal to the real part,
and the second element to the imaginary part, of the complex
number."
When you think about that, then I'm sure you will realize that there are
a number of ways to accomplish what you want.

It's kind of sad if you are forced to do it in an underhanded way by
taking the address, though.

[jacob navia]

user923005 wrote:

On Jan 25, 9:01 pm, Martin Ambuhl <[e-mail]:

>user923005 wrote:

>>It could be implemented (for instance) as:
complex {
double parts[2];
};
or
complex {
double re;
double im;
};
or whatever underneath the covers, since that bit is not specified
clearly.

Of course it might be implemented as the standard specifies:
"Each complex type has the same representation and alignment
requirements as an array type containing exactly two elements of the
corresponding real type; the first element is equal to the real part,
and the second element to the imaginary part, of the complex
number."
When you think about that, then I'm sure you will realize that there are
a number of ways to accomplish what you want.

It's kind of sad if you are forced to do it in an underhanded way by
taking the address, though.

Specially if it could have been so easy to specify a couple of "setter"
functions:

// For double _Complex
double sreal(complex,double); // Returns the old value
double simag(complex,double); // Returns the old value

that would do this!

--
jacob navia
jacob at jacob point remcomp point fr
logiciels/informatique
http://www.cs.virginia.edu/~lcc-win32