Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension


Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
296 changes: 280 additions & 16 deletions src/dmd/backend/bcomplex.d
Original file line number Diff line number Diff line change
Expand Up @@ -4,14 +4,15 @@
*
* Copyright: public domain
* License: public domain
* Source: $(LINK2 https://github.com/dlang/dmd/blob/master/src/dmd/backend/bcomplex.d, backend/_bcomplex.d)
* Source: $(DMDSRC backend/_bcomplex.d)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

FYI underscore escaping is no longer needed. We enabled the opt-out about two months ago.
See: dlang/dlang.org#2307

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Also last year you made a point of using absolute links here as they are clickable (see #7030).
Though I was a fan of DMDSRC, so I don't mind ;-)

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

yes, it should be the full version.

*/

module dmd.backend.bcomplex;

// Online documentation: https://dlang.org/phobos/dmd_backend_bcomplex.html

import dmd.backend.cdef;
public import dmd.root.longdouble : targ_ldouble = longdouble;
import core.stdc.math : fabs, fabsl, sqrt;
version(CRuntime_Microsoft)
private import dmd.root.longdouble : fabsl, sqrt; // needed if longdouble is longdouble_soft

extern (C++):
@nogc:
Expand All @@ -24,28 +25,291 @@ struct Complex_f
{
float re, im;

static Complex_f div(ref Complex_f x, ref Complex_f y);
static Complex_f mul(ref Complex_f x, ref Complex_f y);
static targ_ldouble abs(ref Complex_f z);
static Complex_f sqrtc(ref Complex_f z);
static Complex_f div(ref Complex_f x, ref Complex_f y)
{
Complex_f q;
targ_ldouble r;
targ_ldouble den;

if (fabs(y.re) < fabs(y.im))
{
r = y.re / y.im;
den = y.im + r * y.re;
q.re = cast(float)((x.re * r + x.im) / den);
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Translation looks good, but I don't understand why the casts to float and double in this file are needed.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

they make it clearer what is happening

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

and the casts are necessary if targ_ldouble resolves to longdouble_soft, as there is no implicit conversion to float/double then.

q.im = cast(float)((x.im * r - x.re) / den);
}
else
{
r = y.im / y.re;
den = y.re + r * y.im;
q.re = cast(float)((x.re + r * x.im) / den);
q.im = cast(float)((x.im - r * x.re) / den);
}
return q;
}

static Complex_f mul(ref Complex_f x, ref Complex_f y)
{
Complex_f p;

p.re = x.re * y.re - x.im * y.im;
p.im = x.im * y.re + x.re * y.im;
return p;
}

static targ_ldouble abs(ref Complex_f z)
{
targ_ldouble x,y,ans,temp;

x = fabs(z.re);
y = fabs(z.im);
if (x == 0)
ans = y;
else if (y == 0)
ans = x;
else if (x > y)
{
temp = y / x;
ans = x * sqrt(1 + temp * temp);
}
else
{
temp = x / y;
ans = y * sqrt(1 + temp * temp);
}
return ans;
}

static Complex_f sqrtc(ref Complex_f z)
{
Complex_f c;
targ_ldouble x,y,w,r;

if (z.re == 0 && z.im == 0)
{
c.re = 0;
c.im = 0;
}
else
{
x = fabs(z.re);
y = fabs(z.im);
if (x >= y)
{
r = y / x;
w = sqrt(x) * sqrt(0.5 * (1 + sqrt(1 + r * r)));
}
else
{
r = x / y;
w = sqrt(y) * sqrt(0.5 * (r + sqrt(1 + r * r)));
}
if (z.re >= 0)
{
c.re = cast(float)w;
c.im = cast(float)(targ_ldouble(z.im) / (w + w));
}
else
{
c.im = cast(float)((z.im >= 0) ? w : -w);
c.re = z.im / (c.im + c.im);
}
}
return c;
}
}

struct Complex_d
{
double re, im;

static Complex_d div(ref Complex_d x, ref Complex_d y);
static Complex_d mul(ref Complex_d x, ref Complex_d y);
static targ_ldouble abs(ref Complex_d z);
static Complex_d sqrtc(ref Complex_d z);
static Complex_d div(ref Complex_d x, ref Complex_d y)
{
Complex_d q;
targ_ldouble r;
targ_ldouble den;

if (fabs(y.re) < fabs(y.im))
{
r = y.re / y.im;
den = y.im + r * y.re;
q.re = cast(double)((x.re * r + x.im) / den);
q.im = cast(double)((x.im * r - x.re) / den);
}
else
{
r = y.im / y.re;
den = y.re + r * y.im;
q.re = cast(double)((x.re + r * x.im) / den);
q.im = cast(double)((x.im - r * x.re) / den);
}
return q;
}

static Complex_d mul(ref Complex_d x, ref Complex_d y)
{
Complex_d p;
p.re = x.re * y.re - x.im * y.im;
p.im = x.im * y.re + x.re * y.im;
return p;
}

static targ_ldouble abs(ref Complex_d z)
{
targ_ldouble x,y,ans,temp;
x = fabs(z.re);
y = fabs(z.im);
if (x == 0)
ans = y;
else if (y == 0)
ans = x;
else if (x > y)
{
temp = y / x;
ans = x * sqrt(1 + temp * temp);
}
else
{
temp = x / y;
ans = y * sqrt(1 + temp * temp);
}
return ans;
}

static Complex_d sqrtc(ref Complex_d z)
{
Complex_d c;
targ_ldouble x,y,w,r;

if (z.re == 0 && z.im == 0)
{
c.re = 0;
c.im = 0;
}
else
{
x = fabs(z.re);
y = fabs(z.im);
if (x >= y)
{
r = y / x;
w = sqrt(x) * sqrt(0.5 * (1 + sqrt(1 + r * r)));
}
else
{
r = x / y;
w = sqrt(y) * sqrt(0.5 * (r + sqrt(1 + r * r)));
}
if (z.re >= 0)
{
c.re = cast(double)w;
c.im = cast(double)(targ_ldouble(z.im) / (w + w));
}
else
{
c.im = cast(double)((z.im >= 0) ? w : -w);
c.re = z.im / (2 * c.im);
}
}
return c;
}
}


struct Complex_ld
{
targ_ldouble re, im;

static Complex_ld div(ref Complex_ld x, Complex_ld y);
static Complex_ld mul(ref Complex_ld x, ref Complex_ld y);
static targ_ldouble abs(ref Complex_ld z);
static Complex_ld sqrtc(ref Complex_ld z);
static Complex_ld div(ref Complex_ld x, ref Complex_ld y)
{
Complex_ld q = void;
targ_ldouble r;
targ_ldouble den;

if (fabsl(y.re) < fabsl(y.im))
{
r = y.re / y.im;
den = y.im + r * y.re;
q.re = (x.re * r + x.im) / den;
q.im = (x.im * r - x.re) / den;
}
else
{
r = y.im / y.re;
den = y.re + r * y.im;
q.re = (x.re + r * x.im) / den;
q.im = (x.im - r * x.re) / den;
}
return q;
}

static Complex_ld mul(ref Complex_ld x, ref Complex_ld y)
{
Complex_ld p = void;

p.re = x.re * y.re - x.im * y.im;
p.im = x.im * y.re + x.re * y.im;
return p;
}

static targ_ldouble abs(ref Complex_ld z)
{
targ_ldouble x,y,ans,temp;

x = fabsl(z.re);
y = fabsl(z.im);
if (x == 0)
ans = y;
else if (y == 0)
ans = x;
else if (x > y)
{
temp = y / x;
ans = x * sqrt(1 + temp * temp);
}
else
{
temp = x / y;
ans = y * sqrt(1 + temp * temp);
}
return ans;
}

static Complex_ld sqrtc(ref Complex_ld z)
{
Complex_ld c = void;
targ_ldouble x,y,w,r;

if (z.re == 0 && z.im == 0)
{
c.re = 0;
c.im = 0;
}
else
{
x = fabsl(z.re);
y = fabsl(z.im);
if (x >= y)
{
r = y / x;
w = sqrt(x) * sqrt(0.5 * (1 + sqrt(1 + r * r)));
}
else
{
r = x / y;
w = sqrt(y) * sqrt(0.5 * (r + sqrt(1 + r * r)));
}
if (z.re >= 0)
{
c.re = w;
c.im = z.im / (w + w);
}
else
{
c.im = (z.im >= 0) ? w : -w;
c.re = z.im / (c.im + c.im);
}
}
return c;
}
}
Loading