URL https://opencores.org/ocsvn/openrisc/openrisc/trunk

Subversion Repositories openrisc

[/] [openrisc/] [trunk/] [gnu-dev/] [or1k-gcc/] [libquadmath/] [math/] [complex.c] - Blame information for rev 801

Go to most recent revision | Details | Compare with Previous | View Log


#include "quadmath-imp.h"
 
 
#define REALPART(z) (__real__(z)) 
#define IMAGPART(z) (__imag__(z)) 
#define COMPLEX_ASSIGN(z_, r_, i_) {__real__(z_) = (r_); __imag__(z_) = (i_);} 
 
 
// Horrible... GCC doesn't know how to multiply or divide these
// __complex128 things. We have to do it on our own.
// Protect it around macros so, some day, we can switch it on
 
#if 0
 
# define C128_MULT(x,y) ((x)*(y))
# define C128_DIV(x,y) ((x)/(y))
 
#else
 
#define C128_MULT(x,y) mult_c128(x,y)
#define C128_DIV(x,y) div_c128(x,y)
 
static inline __complex128 mult_c128 (__complex128 x, __complex128 y)
{
  __float128 r1 = REALPART(x), i1 = IMAGPART(x);
  __float128 r2 = REALPART(y), i2 = IMAGPART(y);
  __complex128 res;
  COMPLEX_ASSIGN(res, r1*r2 - i1*i2, i2*r1 + i1*r2);
  return res;
}
 
 
// Careful: the algorithm for the division sucks. A lot.
static inline __complex128 div_c128 (__complex128 x, __complex128 y)
{
  __float128 n = hypotq (REALPART (y), IMAGPART (y));
  __float128 r1 = REALPART(x), i1 = IMAGPART(x);
  __float128 r2 = REALPART(y), i2 = IMAGPART(y);
  __complex128 res;
  COMPLEX_ASSIGN(res, r1*r2 + i1*i2, i1*r2 - i2*r1);
  return res / n;
}
 
#endif
 
 
 
__float128
cabsq (__complex128 z)
{
  return hypotq (REALPART (z), IMAGPART (z));
}
 
 
__complex128
cexpq (__complex128 z)
{
  __float128 a, b;
  __complex128 v;
 
  a = REALPART (z);
  b = IMAGPART (z);
  COMPLEX_ASSIGN (v, cosq (b), sinq (b));
  return expq (a) * v;
}
 
 
__complex128
cexpiq (__float128 x)
{
  __complex128 v;
  COMPLEX_ASSIGN (v, cosq (x), sinq (x));
  return v;
}
 
 
__float128
cargq (__complex128 z)
{
  return atan2q (IMAGPART (z), REALPART (z));
}
 
 
__complex128
clogq (__complex128 z)
{
  __complex128 v;
  COMPLEX_ASSIGN (v, logq (cabsq (z)), cargq (z));
  return v;
}
 
 
__complex128
clog10q (__complex128 z)
{
  __complex128 v;
  COMPLEX_ASSIGN (v, log10q (cabsq (z)), cargq (z));
  return v;
}
 
 
__complex128
cpowq (__complex128 base, __complex128 power)
{
  return cexpq (C128_MULT(power, clogq (base)));
}
 
 
__complex128
csinq (__complex128 a)
{
  __float128 r = REALPART (a), i = IMAGPART (a);
  __complex128 v;
  COMPLEX_ASSIGN (v, sinq (r) * coshq (i), cosq (r) * sinhq (i));
  return v;
}
 
 
__complex128
csinhq (__complex128 a)
{
  __float128 r = REALPART (a), i = IMAGPART (a);
  __complex128 v;
  COMPLEX_ASSIGN (v, sinhq (r) * cosq (i), coshq (r) * sinq (i));
  return v;
}
 
 
__complex128
ccosq (__complex128 a)
{
  __float128 r = REALPART (a), i = IMAGPART (a);
  __complex128 v;
  COMPLEX_ASSIGN (v, cosq (r) * coshq (i), - (sinq (r) * sinhq (i)));
  return v;
}
 
 
__complex128
ccoshq (__complex128 a)
{
  __float128 r = REALPART (a), i = IMAGPART (a);
  __complex128 v;
  COMPLEX_ASSIGN (v, coshq (r) * cosq (i),  sinhq (r) * sinq (i));
  return v;
}
 
 
__complex128
ctanq (__complex128 a)
{
  __float128 rt = tanq (REALPART (a)), it = tanhq (IMAGPART (a));
  __complex128 n, d;
  COMPLEX_ASSIGN (n, rt, it);
  COMPLEX_ASSIGN (d, 1, - (rt * it));
  return C128_DIV(n,d);
}
 
 
__complex128
ctanhq (__complex128 a)
{
  __float128 rt = tanhq (REALPART (a)), it = tanq (IMAGPART (a));
  __complex128 n, d;
  COMPLEX_ASSIGN (n, rt, it);
  COMPLEX_ASSIGN (d, 1, rt * it);
  return C128_DIV(n,d);
}
 
 
/* Square root algorithm from glibc.  */
__complex128
csqrtq (__complex128 z)
{
  __float128 re = REALPART(z), im = IMAGPART(z);
  __complex128 v;
 
  if (im == 0)
  {
    if (re < 0)
    {
      COMPLEX_ASSIGN (v, 0, copysignq (sqrtq (-re), im));
    }
    else
    {
      COMPLEX_ASSIGN (v, fabsq (sqrtq (re)), copysignq (0, im));
    }
  }
  else if (re == 0)
  {
    __float128 r = sqrtq (0.5 * fabsq (im));
    COMPLEX_ASSIGN (v, r, copysignq (r, im));
  }
  else
  {
    __float128 d = hypotq (re, im);
    __float128 r, s;
 
    /* Use the identity   2  Re res  Im res = Im x
        to avoid cancellation error in  d +/- Re x.  */
    if (re > 0)
      r = sqrtq (0.5 * d + 0.5 * re), s = (0.5 * im) / r;
    else
      s = sqrtq (0.5 * d - 0.5 * re), r = fabsq ((0.5 * im) / s);
 
    COMPLEX_ASSIGN (v, r, copysignq (s, im));
  }
  return v;
}
 

Browse

Tools

Subversion Repositories openrisc

[/] [openrisc/] [trunk/] [gnu-dev/] [or1k-gcc/] [libquadmath/] [math/] [complex.c] - Blame information for rev 801

Line No.	Rev	Author	Line
1	740	jeremybenn	`#include "quadmath-imp.h"`
2
3
4			`#define REALPART(z) (__real__(z))`
5			`#define IMAGPART(z) (__imag__(z))`
6			`#define COMPLEX_ASSIGN(z_, r_, i_) {__real__(z_) = (r_); __imag__(z_) = (i_);}`
7
8
9			`// Horrible... GCC doesn't know how to multiply or divide these`
10			`// __complex128 things. We have to do it on our own.`
11			`// Protect it around macros so, some day, we can switch it on`
12
13			`#if 0`
14
15			`# define C128_MULT(x,y) ((x)*(y))`
16			`# define C128_DIV(x,y) ((x)/(y))`
17
18			`#else`
19
20			`#define C128_MULT(x,y) mult_c128(x,y)`
21			`#define C128_DIV(x,y) div_c128(x,y)`
22
23			`static inline __complex128 mult_c128 (__complex128 x, __complex128 y)`
24			`{`
25			`__float128 r1 = REALPART(x), i1 = IMAGPART(x);`
26			`__float128 r2 = REALPART(y), i2 = IMAGPART(y);`
27			`__complex128 res;`
28			`COMPLEX_ASSIGN(res, r1r2 - i1i2, i2r1 + i1r2);`
29			`return res;`
30			`}`
31
32
33			`// Careful: the algorithm for the division sucks. A lot.`
34			`static inline __complex128 div_c128 (__complex128 x, __complex128 y)`
35			`{`
36			`__float128 n = hypotq (REALPART (y), IMAGPART (y));`
37			`__float128 r1 = REALPART(x), i1 = IMAGPART(x);`
38			`__float128 r2 = REALPART(y), i2 = IMAGPART(y);`
39			`__complex128 res;`
40			`COMPLEX_ASSIGN(res, r1r2 + i1i2, i1r2 - i2r1);`
41			`return res / n;`
42			`}`
43
44			`#endif`
45
46
47
48			`__float128`
49			`cabsq (__complex128 z)`
50			`{`
51			`return hypotq (REALPART (z), IMAGPART (z));`
52			`}`
53
54
55			`__complex128`
56			`cexpq (__complex128 z)`
57			`{`
58			`__float128 a, b;`
59			`__complex128 v;`
60
61			`a = REALPART (z);`
62			`b = IMAGPART (z);`
63			`COMPLEX_ASSIGN (v, cosq (b), sinq (b));`
64			`return expq (a) * v;`
65			`}`
66
67
68			`__complex128`
69			`cexpiq (__float128 x)`
70			`{`
71			`__complex128 v;`
72			`COMPLEX_ASSIGN (v, cosq (x), sinq (x));`
73			`return v;`
74			`}`
75
76
77			`__float128`
78			`cargq (__complex128 z)`
79			`{`
80			`return atan2q (IMAGPART (z), REALPART (z));`
81			`}`
82
83
84			`__complex128`
85			`clogq (__complex128 z)`
86			`{`
87			`__complex128 v;`
88			`COMPLEX_ASSIGN (v, logq (cabsq (z)), cargq (z));`
89			`return v;`
90			`}`
91
92
93			`__complex128`
94			`clog10q (__complex128 z)`
95			`{`
96			`__complex128 v;`
97			`COMPLEX_ASSIGN (v, log10q (cabsq (z)), cargq (z));`
98			`return v;`
99			`}`
100
101
102			`__complex128`
103			`cpowq (__complex128 base, __complex128 power)`
104			`{`
105			`return cexpq (C128_MULT(power, clogq (base)));`
106			`}`
107
108
109			`__complex128`
110			`csinq (__complex128 a)`
111			`{`
112			`__float128 r = REALPART (a), i = IMAGPART (a);`
113			`__complex128 v;`
114			`COMPLEX_ASSIGN (v, sinq (r) * coshq (i), cosq (r) * sinhq (i));`
115			`return v;`
116			`}`
117
118
119			`__complex128`
120			`csinhq (__complex128 a)`
121			`{`
122			`__float128 r = REALPART (a), i = IMAGPART (a);`
123			`__complex128 v;`
124			`COMPLEX_ASSIGN (v, sinhq (r) * cosq (i), coshq (r) * sinq (i));`
125			`return v;`
126			`}`
127
128
129			`__complex128`
130			`ccosq (__complex128 a)`
131			`{`
132			`__float128 r = REALPART (a), i = IMAGPART (a);`
133			`__complex128 v;`
134			`COMPLEX_ASSIGN (v, cosq (r) * coshq (i), - (sinq (r) * sinhq (i)));`
135			`return v;`
136			`}`
137
138
139			`__complex128`
140			`ccoshq (__complex128 a)`
141			`{`
142			`__float128 r = REALPART (a), i = IMAGPART (a);`
143			`__complex128 v;`
144			`COMPLEX_ASSIGN (v, coshq (r) * cosq (i), sinhq (r) * sinq (i));`
145			`return v;`
146			`}`
147
148
149			`__complex128`
150			`ctanq (__complex128 a)`
151			`{`
152			`__float128 rt = tanq (REALPART (a)), it = tanhq (IMAGPART (a));`
153			`__complex128 n, d;`
154			`COMPLEX_ASSIGN (n, rt, it);`
155			`COMPLEX_ASSIGN (d, 1, - (rt * it));`
156			`return C128_DIV(n,d);`
157			`}`
158
159
160			`__complex128`
161			`ctanhq (__complex128 a)`
162			`{`
163			`__float128 rt = tanhq (REALPART (a)), it = tanq (IMAGPART (a));`
164			`__complex128 n, d;`
165			`COMPLEX_ASSIGN (n, rt, it);`
166			`COMPLEX_ASSIGN (d, 1, rt * it);`
167			`return C128_DIV(n,d);`
168			`}`
169
170
171			`/* Square root algorithm from glibc. */`
172			`__complex128`
173			`csqrtq (__complex128 z)`
174			`{`
175			`__float128 re = REALPART(z), im = IMAGPART(z);`
176			`__complex128 v;`
177
178			`if (im == 0)`
179			`{`
180			`if (re < 0)`
181			`{`
182			`COMPLEX_ASSIGN (v, 0, copysignq (sqrtq (-re), im));`
183			`}`
184			`else`
185			`{`
186			`COMPLEX_ASSIGN (v, fabsq (sqrtq (re)), copysignq (0, im));`
187			`}`
188			`}`
189			`else if (re == 0)`
190			`{`
191			`__float128 r = sqrtq (0.5 * fabsq (im));`
192			`COMPLEX_ASSIGN (v, r, copysignq (r, im));`
193			`}`
194			`else`
195			`{`
196			`__float128 d = hypotq (re, im);`
197			`__float128 r, s;`
198
199			`/* Use the identity 2 Re res Im res = Im x`
200			`to avoid cancellation error in d +/- Re x. */`
201			`if (re > 0)`
202			`r = sqrtq (0.5 * d + 0.5 * re), s = (0.5 * im) / r;`
203			`else`
204			`s = sqrtq (0.5 * d - 0.5 * re), r = fabsq ((0.5 * im) / s);`
205
206			`COMPLEX_ASSIGN (v, r, copysignq (s, im));`
207			`}`
208			`return v;`
209			`}`
210