100 lines
3.9 KiB
C
100 lines
3.9 KiB
C
// C99 / C11 complex.h — complex-number types and core helpers.
|
|
//
|
|
// clang implements `_Complex` as a built-in type that lowers to a
|
|
// struct-of-two-reals on the W65816 (`_Complex double` = 16 bytes,
|
|
// `_Complex float` = 8 bytes). Plain arithmetic (`a + b`, `a * b`,
|
|
// etc.) is handled by the compiler via softFloat / softDouble.
|
|
//
|
|
// **Supported surface:** the core component / conjugate / magnitude /
|
|
// argument helpers — `creal`, `cimag`, `conj`, `cabs`, `carg`,
|
|
// `cproj` — plus the `CMPLX` constructor macros.
|
|
//
|
|
// **NOT supported:** the transcendental complex routines (`csin`,
|
|
// `ccos`, `cexp`, `clog`, `cpow`, `csqrt`, etc.) — they would each
|
|
// require a real polynomial-expansion implementation; not worth the
|
|
// runtime cost for our IIgs target. Code that references those
|
|
// symbols will link-fail; if you need them, implement them in your
|
|
// project and link them in.
|
|
|
|
#ifndef _COMPLEX_H
|
|
#define _COMPLEX_H
|
|
|
|
#include <math.h>
|
|
|
|
// Per C11: `complex` and `_Complex_I` are macros provided by
|
|
// <complex.h>. Real-world code mostly uses `complex` rather than
|
|
// the underscore form.
|
|
#define complex _Complex
|
|
#define _Complex_I ((float _Complex){0.0f, 1.0f})
|
|
#define I _Complex_I
|
|
|
|
// CMPLX(real, imag) — C11 constructor. Avoid the type-pun trick;
|
|
// clang implements this as a compound literal.
|
|
#define CMPLX(r, i) ((double _Complex){ (r), (i) })
|
|
#define CMPLXF(r, i) ((float _Complex){ (r), (i) })
|
|
#define CMPLXL(r, i) ((double _Complex){ (r), (i) }) // long double = double here
|
|
|
|
// ---- Component access -----------------------------------------------
|
|
// clang provides `__real__` and `__imag__` lvalue extensions that map
|
|
// directly to the underlying real / imag slot of the complex struct.
|
|
// Wrapping them as inline functions avoids leaking the gcc-extension
|
|
// keyword into user code.
|
|
|
|
static inline double creal (double _Complex z) { return __real__ z; }
|
|
static inline double cimag (double _Complex z) { return __imag__ z; }
|
|
static inline float crealf(float _Complex z) { return __real__ z; }
|
|
static inline float cimagf(float _Complex z) { return __imag__ z; }
|
|
static inline double creall(double _Complex z) { return __real__ z; }
|
|
static inline double cimagl(double _Complex z) { return __imag__ z; }
|
|
|
|
// ---- Conjugate -------------------------------------------------------
|
|
// conj(a + b*I) = a - b*I. Implemented via CMPLX so the compiler can
|
|
// optimise away the temporary.
|
|
|
|
static inline double _Complex conj (double _Complex z) {
|
|
return CMPLX(__real__ z, -__imag__ z);
|
|
}
|
|
static inline float _Complex conjf(float _Complex z) {
|
|
return CMPLXF(__real__ z, -__imag__ z);
|
|
}
|
|
static inline double _Complex conjl(double _Complex z) {
|
|
return CMPLX(__real__ z, -__imag__ z);
|
|
}
|
|
|
|
// ---- Magnitude / argument / projection ------------------------------
|
|
// cabs uses hypot to avoid intermediate over/underflow. carg uses
|
|
// atan2. cproj returns z unchanged unless either part is infinite,
|
|
// in which case it returns (INFINITY, +-0).
|
|
|
|
static inline double cabs (double _Complex z) {
|
|
return hypot(__real__ z, __imag__ z);
|
|
}
|
|
static inline float cabsf(float _Complex z) {
|
|
return hypotf(__real__ z, __imag__ z);
|
|
}
|
|
static inline double cabsl(double _Complex z) {
|
|
return hypot(__real__ z, __imag__ z);
|
|
}
|
|
|
|
static inline double carg (double _Complex z) {
|
|
return atan2(__imag__ z, __real__ z);
|
|
}
|
|
static inline float cargf(float _Complex z) {
|
|
return atan2f(__imag__ z, __real__ z);
|
|
}
|
|
static inline double cargl(double _Complex z) {
|
|
return atan2(__imag__ z, __real__ z);
|
|
}
|
|
|
|
static inline double _Complex cproj(double _Complex z) {
|
|
if (__isinf_d(__real__ z) || __isinf_d(__imag__ z)) {
|
|
return CMPLX(HUGE_VAL, __imag__ z < 0.0 ? -0.0 : 0.0);
|
|
}
|
|
return z;
|
|
}
|
|
static inline float _Complex cprojf(float _Complex z) {
|
|
return (float _Complex)cproj((double _Complex)z);
|
|
}
|
|
static inline double _Complex cprojl(double _Complex z) { return cproj(z); }
|
|
|
|
#endif
|