From e2a77e4fec3e283a07614b26c46b7a0bf52fe54b Mon Sep 17 00:00:00 2001 From: feihoo87 Date: Fri, 28 Aug 2026 18:11:15 +0800 Subject: [PATCH 1/4] Fix symbolic frequency filtering in 3.3.1 --- tests/test_core.py | 72 ++++++++ waveforms/_cwaveform.c | 362 ++++++++++++++++++++++++++++++++--------- waveforms/version.py | 2 +- 3 files changed, 357 insertions(+), 79 deletions(-) diff --git a/tests/test_core.py b/tests/test_core.py index ccc30de..2a00d78 100644 --- a/tests/test_core.py +++ b/tests/test_core.py @@ -80,6 +80,78 @@ def test_c_symbolic_frequency_filter(): ) +def test_c_symbolic_filter_reduces_trigonometric_products_and_powers(): + x = np.linspace(-1.0, 1.0, 4097) + envelope = wf.gaussian(2.0) + first_frequency = 9.0 + second_frequency = 7.0 + first_phase = 0.31 + second_phase = -0.27 + cutoff = 5.0 + + first_cos = wf.cos(first_frequency, first_phase) + first_sin = wf.sin(first_frequency, first_phase) + second_cos = wf.cos(second_frequency, second_phase) + second_sin = wf.sin(second_frequency, second_phase) + difference = first_frequency - second_frequency + phase_difference = first_phase - second_phase + + assert np.allclose( + (2 * envelope * first_cos * second_cos).filter(high=cutoff)(x), + (envelope * wf.cos(difference, phase_difference))(x), + atol=2e-11, + ) + assert np.allclose( + (2 * envelope * first_cos * second_sin).filter(high=cutoff)(x), + (-envelope * wf.sin(difference, phase_difference))(x), + atol=2e-11, + ) + assert np.allclose( + (2 * envelope * first_sin * second_sin).filter(high=cutoff)(x), + (envelope * wf.cos(difference, phase_difference))(x), + atol=2e-11, + ) + + powered = wf.cos(first_frequency, first_phase) ** 2 + assert np.allclose(powered.filter(high=cutoff)(x), 0.5, atol=2e-11) + assert np.allclose( + powered.filter(2 * first_frequency, np.inf)(x), + (0.5 * wf.cos(2 * first_frequency, 2 * first_phase))(x), + atol=2e-11, + ) + + retained = (1.5e-15 * first_cos).filter(eps=1e-15) + discarded = (0.5e-15 * first_cos).filter(eps=1e-15) + assert np.max(np.abs(retained(x))) > 1e-15 + assert np.array_equal(discarded(x), np.zeros_like(x)) + + +def test_down_conversion_filter_matches_historic_symbolic_behavior(): + x = np.linspace(-100.0, 100.0, 10001) + envelope = wf.gaussian(100.0) + radio_frequency = 92.0451 + local_frequency = 92.0 + phase = 0.32 + rf, _ = wf.mixing( + envelope, freq=radio_frequency, phase=phase, DRAGScaling=0.0, + ) + + i = (2 * rf * wf.cos(-2 * np.pi * local_frequency)).filter( + high=2 * np.pi * local_frequency, + ) + q = (2 * rf * wf.sin(-2 * np.pi * local_frequency)).filter( + high=2 * np.pi * local_frequency, + ) + difference = 2 * np.pi * (radio_frequency - local_frequency) + + assert np.allclose( + i(x), (envelope * wf.cos(difference, -phase))(x), atol=2e-9, + ) + assert np.allclose( + q(x), (envelope * wf.sin(difference, -phase))(x), atol=2e-9, + ) + + def test_wave_block_roundtrip_pickle_and_numerics(): actual = _pulse() actual.start = -20e-9 diff --git a/waveforms/_cwaveform.c b/waveforms/_cwaveform.c index 51b2961..dd9d279 100644 --- a/waveforms/_cwaveform.c +++ b/waveforms/_cwaveform.c @@ -956,6 +956,19 @@ static uint32_t wf_clone_affine(const cwaveform_wave *source, wf_node *target, uint32_t offset, uint32_t parameter_offset, int64_t delay, double scale, uint32_t *next); +typedef struct wf_spectral_term { + cwaveform_wave *envelope; + double frequency; + double real; + double imag; +} wf_spectral_term; + +typedef struct wf_spectrum { + wf_spectral_term *terms; + size_t count; + size_t capacity; +} wf_spectrum; + static void wf_make_zero_node(wf_node *node) { memset(node, 0, sizeof(*node)); node->op = WF_OP_CONSTANT; @@ -963,105 +976,298 @@ static void wf_make_zero_node(wf_node *node) { node->upper = INT64_MIN; } -static int wf_carrier_count(const wf_node *nodes, uint32_t index, - int *counts, double *frequencies) { - const wf_node *node; - int count; - if (counts[index] >= 0) return counts[index]; - node = nodes + index; +static void wf_spectrum_clear(wf_spectrum *spectrum) { + size_t index; + if (spectrum == NULL) return; + for (index = 0; index < spectrum->count; ++index) + cwaveform_wave_release(spectrum->terms[index].envelope); + free(spectrum->terms); + memset(spectrum, 0, sizeof(*spectrum)); +} + +/* Append one signed complex-exponential component. Equal envelope/frequency + * pairs are combined as they are produced so powers of one carrier grow + * linearly rather than exponentially. The envelope reference is consumed. */ +static int wf_spectrum_append(wf_spectrum *spectrum, + cwaveform_wave *envelope, + double frequency, double real, double imag) { + size_t index; + wf_spectral_term *terms; + size_t capacity; + if (spectrum == NULL || envelope == NULL || !isfinite(frequency) + || !isfinite(real) || !isfinite(imag)) { + cwaveform_wave_release(envelope); + return -1; + } + if (frequency == 0.0) frequency = 0.0; + if (real == 0.0 && imag == 0.0) { + cwaveform_wave_release(envelope); + return 0; + } + for (index = 0; index < spectrum->count; ++index) { + wf_spectral_term *term = spectrum->terms + index; + if (term->frequency == frequency + && cwaveform_wave_equal(term->envelope, envelope)) { + term->real += real; + term->imag += imag; + cwaveform_wave_release(envelope); + if (!isfinite(term->real) || !isfinite(term->imag)) return -1; + if (term->real == 0.0 && term->imag == 0.0) { + cwaveform_wave_release(term->envelope); + spectrum->terms[index] = spectrum->terms[spectrum->count - 1]; + --spectrum->count; + } + return 0; + } + } + if (spectrum->count == spectrum->capacity) { + capacity = spectrum->capacity == 0 ? 8 : spectrum->capacity * 2; + if (capacity < spectrum->capacity + || capacity > SIZE_MAX / sizeof(*terms)) { + cwaveform_wave_release(envelope); + return -1; + } + terms = (wf_spectral_term *)realloc( + spectrum->terms, capacity * sizeof(*terms)); + if (terms == NULL) { + cwaveform_wave_release(envelope); + return -1; + } + spectrum->terms = terms; + spectrum->capacity = capacity; + } + spectrum->terms[spectrum->count].envelope = envelope; + spectrum->terms[spectrum->count].frequency = frequency; + spectrum->terms[spectrum->count].real = real; + spectrum->terms[spectrum->count].imag = imag; + ++spectrum->count; + return 0; +} + +static int wf_spectrum_move(wf_spectrum *target, wf_spectrum *source) { + size_t index; + for (index = 0; index < source->count; ++index) { + wf_spectral_term *term = source->terms + index; + cwaveform_wave *envelope = term->envelope; + term->envelope = NULL; + if (wf_spectrum_append(target, envelope, term->frequency, + term->real, term->imag) != 0) return -1; + } + return 0; +} + +static int wf_spectrum_product(const wf_spectrum *left, + const wf_spectrum *right, + wf_spectrum *result) { + size_t i; + size_t j; + for (i = 0; i < left->count; ++i) { + for (j = 0; j < right->count; ++j) { + const wf_spectral_term *a = left->terms + i; + const wf_spectral_term *b = right->terms + j; + cwaveform_wave *left_envelope = a->envelope; + cwaveform_wave *right_envelope = b->envelope; + cwaveform_wave *envelope; + double real = a->real * b->real - a->imag * b->imag; + double imag = a->real * b->imag + a->imag * b->real; + cwaveform_wave_retain(left_envelope); + cwaveform_wave_retain(right_envelope); + envelope = wf_mul_owned(left_envelope, right_envelope); + if (wf_spectrum_append(result, envelope, + a->frequency + b->frequency, + real, imag) != 0) return -1; + } + } + return 0; +} + +static int wf_spectrum_node(const cwaveform_wave *wave, uint32_t index, + wf_spectrum *result); + +static int wf_spectrum_power(const cwaveform_wave *wave, uint32_t index, + unsigned exponent, wf_spectrum *result) { + wf_spectrum base = {0}; + wf_spectrum power = {0}; + cwaveform_wave *one = cwaveform_wave_constant(1.0); + int status = -1; + if (one == NULL || wf_spectrum_append(&power, one, 0.0, 1.0, 0.0) != 0) + goto done; + if (wf_spectrum_node(wave, index, &base) != 0) goto done; + while (exponent != 0) { + if (exponent & 1u) { + wf_spectrum next = {0}; + if (wf_spectrum_product(&power, &base, &next) != 0) { + wf_spectrum_clear(&next); + goto done; + } + wf_spectrum_clear(&power); + power = next; + } + exponent >>= 1; + if (exponent != 0) { + wf_spectrum next = {0}; + if (wf_spectrum_product(&base, &base, &next) != 0) { + wf_spectrum_clear(&next); + goto done; + } + wf_spectrum_clear(&base); + base = next; + } + } + if (wf_spectrum_move(result, &power) != 0) goto done; + status = 0; +done: + wf_spectrum_clear(&base); + wf_spectrum_clear(&power); + return status; +} + +/* Convert the real expression tree into signed complex-exponential terms. + * Non-trigonometric nodes are retained as envelopes, matching the historic + * symbolic filter which treats their spectrum as baseband. */ +static int wf_spectrum_node(const cwaveform_wave *wave, uint32_t index, + wf_spectrum *result) { + const wf_node *node = wave->nodes + index; + wf_spectrum left = {0}; + wf_spectrum right = {0}; + wf_spectrum product = {0}; + cwaveform_wave *envelope; + size_t term_index; + int status = -1; switch (node->op) { + case WF_OP_CONSTANT: + envelope = cwaveform_wave_constant(1.0); + return wf_spectrum_append(result, envelope, 0.0, node->p0, 0.0); case WF_OP_COS: - case WF_OP_SIN: - counts[index] = 1; - frequencies[index] = fabs(node->p0); - return 1; - case WF_OP_ADD: - case WF_OP_MUL: { - int left = wf_carrier_count(nodes, node->left, - counts, frequencies); - int right = wf_carrier_count(nodes, node->right, - counts, frequencies); - count = left + right; - if (count > 2) count = 2; - counts[index] = count; - if (count == 1) - frequencies[index] = left == 1 - ? frequencies[node->left] : frequencies[node->right]; - return count; + case WF_OP_SIN: { + const double half = 0.5; + const double two_pi = 6.28318530717958647692; + double phase = remainder( + -node->p0 * (double)node->shift + / (double)wf_ticks_per_second, + two_pi); + double cosine = cos(phase); + double sine = sin(phase); + envelope = cwaveform_wave_constant(1.0); + if (node->op == WF_OP_COS) { + if (wf_spectrum_append(result, envelope, node->p0, + half * cosine, half * sine) != 0) + return -1; + envelope = cwaveform_wave_constant(1.0); + return wf_spectrum_append(result, envelope, -node->p0, + half * cosine, -half * sine); + } + if (wf_spectrum_append(result, envelope, node->p0, + half * sine, -half * cosine) != 0) + return -1; + envelope = cwaveform_wave_constant(1.0); + return wf_spectrum_append(result, envelope, -node->p0, + half * sine, half * cosine); } + case WF_OP_ADD: + if (wf_spectrum_node(wave, node->left, &left) != 0 + || wf_spectrum_node(wave, node->right, &right) != 0 + || wf_spectrum_move(result, &left) != 0 + || wf_spectrum_move(result, &right) != 0) goto done; + status = 0; + break; + case WF_OP_MUL: + if (wf_spectrum_node(wave, node->left, &left) != 0 + || wf_spectrum_node(wave, node->right, &right) != 0 + || wf_spectrum_product(&left, &right, &product) != 0 + || wf_spectrum_move(result, &product) != 0) goto done; + status = 0; + break; case WF_OP_SCALE: + if (wf_spectrum_node(wave, node->left, &left) != 0) goto done; + for (term_index = 0; term_index < left.count; ++term_index) { + left.terms[term_index].real *= node->p0; + left.terms[term_index].imag *= node->p0; + } + if (wf_spectrum_move(result, &left) != 0) goto done; + status = 0; + break; case WF_OP_POWER: - case WF_OP_WINDOW: - count = wf_carrier_count(nodes, node->left, - counts, frequencies); - counts[index] = count; - frequencies[index] = frequencies[node->left]; - return count; + if (node->p0 >= 0.0 && node->p0 <= (double)UINT_MAX) + return wf_spectrum_power( + wave, node->left, (unsigned)node->p0, result); + envelope = wf_subwave(wave, index); + return wf_spectrum_append(result, envelope, 0.0, 1.0, 0.0); + case WF_OP_WINDOW: { + int64_t upper; + if (wf_spectrum_node(wave, node->left, &left) != 0) goto done; + memcpy(&upper, &node->p0, sizeof(upper)); + for (term_index = 0; term_index < left.count; ++term_index) { + cwaveform_wave *windowed = cwaveform_wave_window( + left.terms[term_index].envelope, node->shift, upper); + if (windowed == NULL) goto done; + cwaveform_wave_release(left.terms[term_index].envelope); + left.terms[term_index].envelope = windowed; + } + if (wf_spectrum_move(result, &left) != 0) goto done; + status = 0; + break; + } default: - counts[index] = 0; - frequencies[index] = 0.0; - return 0; + envelope = wf_subwave(wave, index); + return wf_spectrum_append(result, envelope, 0.0, 1.0, 0.0); } -} - -static void wf_filter_context(wf_node *nodes, uint32_t index, - int *counts, double *frequencies, - double low, double high) { - wf_node *node = nodes + index; - if (node->op == WF_OP_ADD) { - wf_filter_context(nodes, node->left, counts, frequencies, low, high); - wf_filter_context(nodes, node->right, counts, frequencies, low, high); - return; - } - if (node->op == WF_OP_SCALE || node->op == WF_OP_WINDOW) { - wf_filter_context(nodes, node->left, counts, frequencies, low, high); - return; - } - { - int count = wf_carrier_count(nodes, index, counts, frequencies); - if ((count == 0 && low > 0.0) - || (count == 1 && !(low <= frequencies[index] - && frequencies[index] < high))) - wf_make_zero_node(node); +done: + wf_spectrum_clear(&left); + wf_spectrum_clear(&right); + wf_spectrum_clear(&product); + return status; +} + +static cwaveform_wave *wf_spectral_term_wave(const wf_spectral_term *term) { + cwaveform_wave *envelope = term->envelope; + cwaveform_wave *carrier; + double amplitude; + if (term->frequency == 0.0) { + if (term->real == 0.0) + return cwaveform_wave_constant(0.0); + cwaveform_wave_retain(envelope); + return wf_scaled(envelope, term->real); } + /* Only the positive half of a real signal's conjugate spectrum is + * reconstructed; double it here to recover the real carrier amplitude. */ + amplitude = 2.0 * hypot(term->real, term->imag); + if (amplitude == 0.0) return cwaveform_wave_constant(0.0); + carrier = cwaveform_wave_cos( + term->frequency, atan2(term->imag, term->real)); + carrier = wf_scaled(carrier, amplitude); + cwaveform_wave_retain(envelope); + return wf_mul_owned(envelope, carrier); } cwaveform_wave *cwaveform_wave_filter(const cwaveform_wave *wave, double low, double high, double epsilon) { - wf_node *nodes; - int *counts; - double *frequencies; - cwaveform_wave *result; - uint32_t index; - (void)epsilon; + wf_spectrum spectrum = {0}; + cwaveform_wave *result = NULL; + size_t index; if (wave == NULL || !isfinite(low) || isnan(high) || low < 0.0 - || high < low) return NULL; - nodes = (wf_node *)malloc((size_t)wave->node_count * sizeof(*nodes)); - counts = (int *)malloc((size_t)wave->node_count * sizeof(*counts)); - frequencies = (double *)calloc(wave->node_count, sizeof(*frequencies)); - if (nodes == NULL || counts == NULL || frequencies == NULL) { - free(nodes); free(counts); free(frequencies); - return NULL; - } - memcpy(nodes, wave->nodes, (size_t)wave->node_count * sizeof(*nodes)); - for (index = 0; index < wave->node_count; ++index) counts[index] = -1; - wf_filter_context(nodes, wave->root, counts, frequencies, low, high); - for (index = 0; index < wave->node_count; ++index) { - if (wf_prepare_decoded_node(nodes + index, nodes, index, - wave->parameter_count) != 0) { - free(nodes); free(counts); free(frequencies); - return NULL; + || high < low || !isfinite(epsilon) || epsilon < 0.0) return NULL; + if (wf_spectrum_node(wave, wave->root, &spectrum) != 0) goto done; + result = cwaveform_wave_constant(0.0); + if (result == NULL) goto done; + for (index = 0; index < spectrum.count; ++index) { + const wf_spectral_term *term = spectrum.terms + index; + double frequency = fabs(term->frequency); + if (term->frequency < 0.0) continue; + if (low <= frequency && frequency < high) { + cwaveform_wave *component = wf_spectral_term_wave(term); + result = wf_add_owned(result, component); + if (result == NULL) goto done; } } - result = wf_wave_from_parts(nodes, wave->node_count, wave->root, - wave->parameters, wave->parameter_count); - free(nodes); free(counts); free(frequencies); if (result != NULL) { cwaveform_wave *compact = cwaveform_wave_simplify(result, epsilon); cwaveform_wave_release(result); result = compact; } +done: + wf_spectrum_clear(&spectrum); return result; } diff --git a/waveforms/version.py b/waveforms/version.py index 20a9f63..6f292e0 100644 --- a/waveforms/version.py +++ b/waveforms/version.py @@ -1,2 +1,2 @@ """Define version number here and read it from setup.py automatically""" -__version__ = "3.3.0" +__version__ = "3.3.1" From f35904b0bfc8d02a4dce0cb96b30b7b71b45d539 Mon Sep 17 00:00:00 2001 From: feihoo87 Date: Sat, 29 Aug 2026 11:46:39 +0800 Subject: [PATCH 2/4] Optimize waveform sampling with cross-platform SIMD --- .github/workflows/workflow.yml | 4 + benchmarks/benchmark_simd.py | 80 ++++ benchmarks/benchmark_simd_c.c | 139 +++++++ setup.py | 4 + tests/test_waveform.py | 24 ++ waveforms/_cwaveform.c | 655 +++++++++++++++++++++++++++++++-- waveforms/_cwaveform.h | 3 + waveforms/_waveform.pyx | 73 ++-- 8 files changed, 893 insertions(+), 89 deletions(-) create mode 100644 benchmarks/benchmark_simd.py create mode 100644 benchmarks/benchmark_simd_c.c diff --git a/.github/workflows/workflow.yml b/.github/workflows/workflow.yml index 266ea43..6949e79 100644 --- a/.github/workflows/workflow.yml +++ b/.github/workflows/workflow.yml @@ -48,6 +48,10 @@ jobs: - name: Test with pytest run: | coverage run --source=waveforms -m pytest --verbose tests/ + - name: Benchmark SIMD hot paths + if: ${{ matrix.python-version == '3.12' }} + run: | + python benchmarks/benchmark_simd.py --count 1000000 --repeats 7 - name: Coveralls if: ${{ matrix.os == 'ubuntu-latest' && matrix.python-version == '3.12' }} continue-on-error: true diff --git a/benchmarks/benchmark_simd.py b/benchmarks/benchmark_simd.py new file mode 100644 index 0000000..8a83b3b --- /dev/null +++ b/benchmarks/benchmark_simd.py @@ -0,0 +1,80 @@ +"""Benchmark the loops targeted by native SIMD and batched-math kernels.""" + +from __future__ import annotations + +import argparse +import gc +import json +import time + +import numpy as np + +import waveforms as wf +from waveforms._waveform import quantize_samples + + +def best_time(function, repeats): + function() + timings = [] + for _ in range(repeats): + gc.collect() + start = time.perf_counter() + function() + timings.append(time.perf_counter() - start) + return min(timings) + + +def benchmark(count, rate, repeats): + duration = count / rate + start = -duration / 2 + stop = start + duration + width = duration * 1.8 + pulse = 0.8 * wf.gaussian(width) * wf.cos(2 * np.pi * 100e6) + powered = pulse**2 + positions = np.linspace(start, stop, count, endpoint=False) + values = 0.999 * np.sin(np.linspace(-100, 100, count)) + + pulse.start = start + pulse.stop = stop + pulse.sample_rate = rate + + cases = { + "evaluate_supported": lambda: pulse(positions), + "evaluate_power_fallback": lambda: powered(positions), + "sample_float64": lambda: pulse.sample(dtype=np.float64), + "sample_int16": lambda: pulse.sample(dtype=np.int16), + "sample_int32": lambda: pulse.sample(dtype=np.int32), + "quantize_int16": lambda: quantize_samples(values, 16), + "quantize_int32": lambda: quantize_samples(values, 32), + } + outputs = {name: function() for name, function in cases.items()} + if not np.array_equal( + outputs["sample_int16"], quantize_samples(outputs["sample_float64"], 16) + ): + raise AssertionError("int16 direct sampling differs from float quantization") + if not np.array_equal( + outputs["sample_int32"], quantize_samples(outputs["sample_float64"], 32) + ): + raise AssertionError("int32 direct sampling differs from float quantization") + + return { + "count": count, + "rate": rate, + "milliseconds": { + name: 1e3 * best_time(function, repeats) + for name, function in cases.items() + }, + } + + +def main(): + parser = argparse.ArgumentParser() + parser.add_argument("--count", type=int, default=1_000_000) + parser.add_argument("--rate", type=int, default=2_400_000_000) + parser.add_argument("--repeats", type=int, default=7) + args = parser.parse_args() + print(json.dumps(benchmark(args.count, args.rate, args.repeats), indent=2)) + + +if __name__ == "__main__": + main() diff --git a/benchmarks/benchmark_simd_c.c b/benchmarks/benchmark_simd_c.c new file mode 100644 index 0000000..e32b60a --- /dev/null +++ b/benchmarks/benchmark_simd_c.c @@ -0,0 +1,139 @@ +#define _POSIX_C_SOURCE 200809L + +#include "../waveforms/_cwaveform.h" + +#include +#include +#include +#include +#include +#include + +static double now_seconds(void) { + struct timespec value; + clock_gettime(CLOCK_MONOTONIC, &value); + return (double)value.tv_sec + 1e-9 * (double)value.tv_nsec; +} + +static double best_evaluate(const cwaveform_wave *wave, const double *positions, + size_t count, double *output, int repeats) { + double best = INFINITY; + int repeat; + for (repeat = -1; repeat < repeats; ++repeat) { + double start = now_seconds(); + double elapsed; + if (cwaveform_wave_evaluate( + wave, positions, count, 0, 1.0, -INFINITY, INFINITY, + output) != 0) abort(); + elapsed = now_seconds() - start; + if (repeat >= 0 && elapsed < best) best = elapsed; + } + return best; +} + +static double best_sample(const cwaveform_wave *wave, int64_t start_tick, + size_t count, int dtype, void *output, int repeats) { + double best = INFINITY; + int repeat; + for (repeat = -1; repeat < repeats; ++repeat) { + double start = now_seconds(); + double elapsed; + if (cwaveform_wave_sample( + wave, start_tick, count, 50, 1, 0, 1.0, + -INFINITY, INFINITY, dtype, 1.0, output) != 0) abort(); + elapsed = now_seconds() - start; + if (repeat >= 0 && elapsed < best) best = elapsed; + } + return best; +} + +static double best_quantize(const double *values, size_t count, int dtype, + void *output, int repeats) { + double best = INFINITY; + int repeat; + for (repeat = -1; repeat < repeats; ++repeat) { + double start = now_seconds(); + double elapsed; + if (cwaveform_quantize(values, count, dtype, 1.0, output) != 0) abort(); + elapsed = now_seconds() - start; + if (repeat >= 0 && elapsed < best) best = elapsed; + } + return best; +} + +int main(int argc, char **argv) { + const size_t count = argc > 1 ? (size_t)strtoull(argv[1], NULL, 10) + : 1000000u; + const int repeats = argc > 2 ? atoi(argv[2]) : 15; + const double rate = 2400000000.0; + const double duration = (double)count / rate; + const double start = -duration / 2.0; + const double pi = 3.14159265358979323846; + cwaveform_wave *gaussian = cwaveform_wave_gaussian(duration * 1.8); + cwaveform_wave *carrier = cwaveform_wave_cos(2.0 * pi * 100e6, 0.0); + cwaveform_wave *pulse = cwaveform_wave_mul_affine( + gaussian, 0, 0.8, carrier, 0, 1.0); + cwaveform_wave *powered = cwaveform_wave_power(pulse, 2); + double *positions = malloc(count * sizeof(*positions)); + double *values = malloc(count * sizeof(*values)); + double *samples = malloc(count * sizeof(*samples)); + int16_t *samples16 = malloc(count * sizeof(*samples16)); + int32_t *samples32 = malloc(count * sizeof(*samples32)); + int16_t *expected16 = malloc(count * sizeof(*expected16)); + int32_t *expected32 = malloc(count * sizeof(*expected32)); + size_t index; + double evaluate_supported; + double evaluate_power; + double sample_float64; + double sample_int16; + double sample_int32; + double quantize_int16; + double quantize_int32; + if (pulse == NULL || powered == NULL || positions == NULL || values == NULL + || samples == NULL || samples16 == NULL || samples32 == NULL + || expected16 == NULL || expected32 == NULL) abort(); + for (index = 0; index < count; ++index) { + positions[index] = start + (double)index / rate; + values[index] = 0.999 * sin(-100.0 + 200.0 * (double)index + / (double)count); + } + evaluate_supported = best_evaluate( + pulse, positions, count, samples, repeats); + evaluate_power = best_evaluate(powered, positions, count, samples, repeats); + sample_float64 = best_sample( + pulse, llround(start * 120000000000.0), count, + CWAVEFORM_FLOAT64, samples, repeats); + sample_int16 = best_sample( + pulse, llround(start * 120000000000.0), count, + CWAVEFORM_INT16, samples16, repeats); + sample_int32 = best_sample( + pulse, llround(start * 120000000000.0), count, + CWAVEFORM_INT32, samples32, repeats); + quantize_int16 = best_quantize( + values, count, CWAVEFORM_INT16, expected16, repeats); + quantize_int32 = best_quantize( + values, count, CWAVEFORM_INT32, expected32, repeats); + if (cwaveform_wave_sample( + pulse, llround(start * 120000000000.0), count, 50, 1, 0, 1.0, + -INFINITY, INFINITY, CWAVEFORM_FLOAT64, 1.0, samples) != 0 + || cwaveform_quantize( + samples, count, CWAVEFORM_INT16, 1.0, expected16) != 0 + || cwaveform_quantize( + samples, count, CWAVEFORM_INT32, 1.0, expected32) != 0 + || memcmp(samples16, expected16, count * sizeof(*samples16)) != 0 + || memcmp(samples32, expected32, count * sizeof(*samples32)) != 0) + abort(); + printf("{\n \"count\": %zu,\n \"milliseconds\": {\n", count); + printf(" \"evaluate_supported\": %.6f,\n", 1e3 * evaluate_supported); + printf(" \"evaluate_power\": %.6f,\n", 1e3 * evaluate_power); + printf(" \"sample_float64\": %.6f,\n", 1e3 * sample_float64); + printf(" \"sample_int16\": %.6f,\n", 1e3 * sample_int16); + printf(" \"sample_int32\": %.6f,\n", 1e3 * sample_int32); + printf(" \"quantize_int16\": %.6f,\n", 1e3 * quantize_int16); + printf(" \"quantize_int32\": %.6f\n }\n}\n", 1e3 * quantize_int32); + free(positions); free(values); free(samples); free(samples16); + free(samples32); free(expected16); free(expected32); + cwaveform_wave_release(powered); cwaveform_wave_release(pulse); + cwaveform_wave_release(carrier); cwaveform_wave_release(gaussian); + return 0; +} diff --git a/setup.py b/setup.py index 6c8bf82..a6d8526 100644 --- a/setup.py +++ b/setup.py @@ -39,15 +39,19 @@ def get_extensions(): sources = [os.path.join(dirpath, filename)] include_dirs = [] extra_link_args = [] + libraries = [] if filename == '_waveform.pyx': sources.append(os.path.join( 'waveforms', '_cwaveform.c')) include_dirs.append('waveforms') if sys.platform == 'darwin': extra_link_args.extend(['-framework', 'Accelerate']) + elif sys.platform.startswith('linux'): + libraries.append('m') extensions.append( Extension(module_name(dirpath, filename), sources, include_dirs=include_dirs, + libraries=libraries, extra_link_args=extra_link_args)) return extensions diff --git a/tests/test_waveform.py b/tests/test_waveform.py index 94c9243..3aca468 100644 --- a/tests/test_waveform.py +++ b/tests/test_waveform.py @@ -473,6 +473,15 @@ def test_fixed_width_quantization_supported_sampling_and_fallback(): 1073741824, 2147483647, 2147483647], ) + # Exercise the fused SIMD path, including its historic half-away-from-zero + # rounding rule and multidimensional ``out`` handling. + half_steps = np.array([0.5, -0.5, 1.5, -1.5] * 8) / 32768.0 + target = np.empty((8, 4), dtype=np.int16) + assert quantize_samples(half_steps.reshape(8, 4), 16, out=target) is target + assert np.array_equal( + target.reshape(-1), np.array([1, -1, 2, -2] * 8, dtype=np.int16) + ) + wav = 0.8 * wf.gaussian(20e-9) wav.start = -20e-9 wav.stop = 20e-9 @@ -493,6 +502,21 @@ def test_fixed_width_quantization_supported_sampling_and_fallback(): assert np.array_equal(wav.sample(odd_rate), wav(legacy_grid)) +def test_simd_node_evaluation_matches_scalar_boundaries(): + positions = np.linspace(-2.0, 2.0, 513) + waves = ( + wf.t(), + wf.exp(0.2), + wf.sinc(1.3), + wf.cosh(0.3), + wf.sinh(0.3), + (wf.gaussian(4.0, 1.5) * wf.cos(2.2)) ** 2, + ) + for wave in waves: + expected = np.array([wave(float(position)) for position in positions]) + assert np.allclose(wave(positions), expected, rtol=1e-12, atol=1e-13) + + def test_integer_sampling_fast_paths_are_bit_exact_and_pickle_safe(): rate = 2_400_000_000 pulse = 0.8 * wf.gaussian(20e-9) * wf.cos(2 * np.pi * 100e6) diff --git a/waveforms/_cwaveform.c b/waveforms/_cwaveform.c index dd9d279..1fcd0b2 100644 --- a/waveforms/_cwaveform.c +++ b/waveforms/_cwaveform.c @@ -7,6 +7,31 @@ #include #include +#if defined(__aarch64__) && defined(__ARM_NEON) +#include +#define WF_HAVE_ARM64_NEON 1 +#endif + +#if ((defined(__x86_64__) || defined(__i386__)) \ + && (defined(__GNUC__) || defined(__clang__))) \ + || (defined(_MSC_VER) && (defined(_M_X64) || defined(_M_IX86))) +#include +#if defined(_MSC_VER) +#include +#endif +#define WF_HAVE_X86_SIMD 1 +#endif + +#if defined(WF_HAVE_X86_SIMD) \ + && (defined(__GNUC__) || defined(__clang__)) +#define WF_TARGET_AVX2 __attribute__((target("avx2"))) +#define WF_TARGET_AVX512 \ + __attribute__((target("avx512f,avx512dq,avx512bw,avx512vl"))) +#else +#define WF_TARGET_AVX2 +#define WF_TARGET_AVX512 +#endif + #if defined(__APPLE__) #include #endif @@ -123,6 +148,22 @@ struct cwaveform_sample_plan { static double wf_node_parameter(const cwaveform_wave *wave, const wf_node *node, size_t index); +static double wf_hypot(double x, double y) { + double high = fabs(x); + double low = fabs(y); + double ratio; + if (isinf(high) || isinf(low)) return INFINITY; + if (isnan(high) || isnan(low)) return NAN; + if (low > high) { + double temporary = high; + high = low; + low = temporary; + } + if (high == 0.0) return 0.0; + ratio = low / high; + return high * sqrt(1.0 + ratio * ratio); +} + static void wf_put_u16(uint8_t *p, uint16_t value) { p[0] = (uint8_t)value; p[1] = (uint8_t)(value >> 8); @@ -1231,7 +1272,7 @@ static cwaveform_wave *wf_spectral_term_wave(const wf_spectral_term *term) { } /* Only the positive half of a real signal's conjugate spectrum is * reconstructed; double it here to recover the real carrier amplitude. */ - amplitude = 2.0 * hypot(term->real, term->imag); + amplitude = 2.0 * wf_hypot(term->real, term->imag); if (amplitude == 0.0) return cwaveform_wave_constant(0.0); carrier = cwaveform_wave_cos( term->frequency, atan2(term->imag, term->real)); @@ -2285,7 +2326,7 @@ static double wf_drag_sin_value(const cwaveform_wave *wave, peak_components[0] += transform[index * 4 + 0] * values[index]; peak_components[1] += transform[index * 4 + 2] * values[index]; } - normalization = hypot(peak_components[0], peak_components[1]); + normalization = wf_hypot(peak_components[0], peak_components[1]); /* Restore the actual derivatives after the peak calculation. */ { double midpoint = t0 + width / 2.0; @@ -2509,6 +2550,8 @@ static double wf_evaluate_one(const cwaveform_wave *wave, double position, return values[wave->root]; } +#if (defined(__APPLE__) || defined(WF_HAVE_X86_SIMD) || defined(_WIN32)) \ + && !defined(WF_DISABLE_BATCH_EVALUATOR) #if defined(__APPLE__) static void wf_vector_exp(double *output, const double *input, size_t count) { while (count != 0) { @@ -2530,9 +2573,156 @@ static void wf_vector_cos(double *output, const double *input, size_t count) { } } -/* Evaluate one node at a time so vForce can process transcendental functions - * in wide batches. The portable path below deliberately stays scalar. */ -static int wf_evaluate_many_apple( +static void wf_vector_sin(double *output, const double *input, size_t count) { + while (count != 0) { + int batch = count > (size_t)INT_MAX ? INT_MAX : (int)count; + vvsin(output, input, &batch); + output += batch; + input += batch; + count -= (size_t)batch; + } +} + +static void wf_vector_cosh(double *output, const double *input, size_t count) { + while (count != 0) { + int batch = count > (size_t)INT_MAX ? INT_MAX : (int)count; + vvcosh(output, input, &batch); + output += batch; + input += batch; + count -= (size_t)batch; + } +} + +static void wf_vector_sinh(double *output, const double *input, size_t count) { + while (count != 0) { + int batch = count > (size_t)INT_MAX ? INT_MAX : (int)count; + vvsinh(output, input, &batch); + output += batch; + input += batch; + count -= (size_t)batch; + } +} + +static void wf_vector_power_scalar(double *output, const double *input, + double exponent, size_t count) { + if (exponent == 2.0) { + vDSP_vsqD(input, 1, output, 1, (vDSP_Length)count); + return; + } + while (count != 0) { + int batch = count > (size_t)INT_MAX ? INT_MAX : (int)count; + vvpows(output, &exponent, input, &batch); + output += batch; + input += batch; + count -= (size_t)batch; + } +} +#else +#if defined(__GLIBC__) && defined(__x86_64__) \ + && defined(WF_HAVE_X86_SIMD) \ + && !defined(WF_DISABLE_X86_SIMD) +typedef double wf_v4df __attribute__((vector_size(32))); +typedef double wf_v8df __attribute__((vector_size(64))); +#if defined(WF_DISABLE_AVX512) +#define WF_AVX512_ENABLED 0 +#else +#define WF_AVX512_ENABLED 1 +#endif + +#define WF_DEFINE_VECTOR_MATH(name) \ + extern wf_v4df _ZGVdN4v_##name(wf_v4df); \ + extern wf_v8df _ZGVeN8v_##name(wf_v8df); \ + __attribute__((target("avx2"))) \ + static void wf_vector_##name##_avx2( \ + double *output, const double *input, size_t count) { \ + size_t index = 0; \ + for (; index + 4 <= count; index += 4) { \ + wf_v4df value = (wf_v4df)_mm256_loadu_pd(input + index); \ + _mm256_storeu_pd(output + index, \ + (__m256d)_ZGVdN4v_##name(value)); \ + } \ + for (; index < count; ++index) output[index] = name(input[index]); \ + } \ + __attribute__((target("avx512f"))) \ + static void wf_vector_##name##_avx512( \ + double *output, const double *input, size_t count) { \ + size_t index = 0; \ + for (; index + 8 <= count; index += 8) { \ + wf_v8df value = (wf_v8df)_mm512_loadu_pd(input + index); \ + _mm512_storeu_pd(output + index, \ + (__m512d)_ZGVeN8v_##name(value)); \ + } \ + for (; index < count; ++index) output[index] = name(input[index]); \ + } \ + static void wf_vector_##name( \ + double *output, const double *input, size_t count) { \ + if (count >= 8) { \ + if (WF_AVX512_ENABLED \ + && __builtin_cpu_supports("avx512f")) { \ + wf_vector_##name##_avx512(output, input, count); \ + return; \ + } \ + if (__builtin_cpu_supports("avx2")) { \ + wf_vector_##name##_avx2(output, input, count); \ + return; \ + } \ + } \ + { \ + size_t index; \ + for (index = 0; index < count; ++index) \ + output[index] = name(input[index]); \ + } \ + } + +WF_DEFINE_VECTOR_MATH(exp) +WF_DEFINE_VECTOR_MATH(cos) +WF_DEFINE_VECTOR_MATH(sin) +#undef WF_DEFINE_VECTOR_MATH +#undef WF_AVX512_ENABLED + +/* glibc only added vector cosh/sinh in 2.35, while Linux wheels target + * manylinux_2_34. Keep these uncommon nodes scalar to preserve that ABI. */ +static void wf_vector_cosh(double *output, const double *input, size_t count) { + size_t index; + for (index = 0; index < count; ++index) output[index] = cosh(input[index]); +} + +static void wf_vector_sinh(double *output, const double *input, size_t count) { + size_t index; + for (index = 0; index < count; ++index) output[index] = sinh(input[index]); +} +#else +#define WF_DEFINE_SCALAR_MATH(name) \ + static void wf_vector_##name( \ + double *output, const double *input, size_t count) { \ + size_t index; \ + for (index = 0; index < count; ++index) \ + output[index] = name(input[index]); \ + } +WF_DEFINE_SCALAR_MATH(exp) +WF_DEFINE_SCALAR_MATH(cos) +WF_DEFINE_SCALAR_MATH(sin) +WF_DEFINE_SCALAR_MATH(cosh) +WF_DEFINE_SCALAR_MATH(sinh) +#undef WF_DEFINE_SCALAR_MATH +#endif + +static void wf_vector_power_scalar(double *output, const double *input, + double exponent, size_t count) { + size_t index; + if (exponent == 2.0) { + for (index = 0; index < count; ++index) + output[index] = input[index] * input[index]; + } else { + for (index = 0; index < count; ++index) + output[index] = pow(input[index], exponent); + } +} +#endif + +/* Evaluate one node at a time so platform vector-math libraries can process + * transcendental functions in wide batches. */ +static int wf_evaluate_chunk_vector( const cwaveform_wave *wave, const double *positions, size_t count, double delay, double scale, double lower_clip, double upper_clip, double *output) { @@ -2540,7 +2730,9 @@ static int wf_evaluate_many_apple( double *scratch; uint32_t node_index; size_t index; - if (count != 0 && (size_t)wave->node_count > SIZE_MAX / count) return -2; + if (count > SIZE_MAX / sizeof(double) + || (count != 0 && (size_t)wave->node_count + > (SIZE_MAX / sizeof(double)) / count)) return -2; matrix = (double *)malloc((size_t)wave->node_count * count * sizeof(double)); scratch = (double *)malloc(count * sizeof(double)); if (matrix == NULL || scratch == NULL) { @@ -2558,16 +2750,22 @@ static int wf_evaluate_many_apple( double shift = (double)node->shift / (double)wf_ticks_per_second; switch (node->op) { case WF_OP_CONSTANT: - for (index = 0; index < count; ++index) { - double position = positions[index] - delay; - row[index] = position >= lower && position < upper - ? node->p0 : 0.0; + if (lower == -DBL_MAX && upper == DBL_MAX) { + for (index = 0; index < count; ++index) + row[index] = node->p0; + } else { + for (index = 0; index < count; ++index) { + double position = positions[index] - delay; + row[index] = position >= lower && position < upper + ? node->p0 : 0.0; + } } break; case WF_OP_GAUSSIAN: for (index = 0; index < count; ++index) { double position = positions[index] - delay; - double local = (position - shift) / node->p1; + double std = node->flags ? node->p0 : node->p1; + double local = (position - shift) / std; scratch[index] = position >= lower && position < upper ? -(local * local) : -INFINITY; } @@ -2582,14 +2780,7 @@ static int wf_evaluate_many_apple( if (node->op == WF_OP_COS) { wf_vector_cos(row, scratch, count); } else { - int batch = count > (size_t)INT_MAX ? INT_MAX : (int)count; - size_t cursor = 0; - while (cursor < count) { - batch = count - cursor > (size_t)INT_MAX - ? INT_MAX : (int)(count - cursor); - vvsin(row + cursor, scratch + cursor, &batch); - cursor += (size_t)batch; - } + wf_vector_sin(row, scratch, count); } if (lower != -DBL_MAX || upper != DBL_MAX) { for (index = 0; index < count; ++index) { @@ -2598,6 +2789,44 @@ static int wf_evaluate_many_apple( } } break; + case WF_OP_LINEAR: + for (index = 0; index < count; ++index) { + double position = positions[index] - delay; + row[index] = position >= lower && position < upper + ? position - shift : 0.0; + } + break; + case WF_OP_ERF: + for (index = 0; index < count; ++index) + scratch[index] = (positions[index] - delay - shift) / node->p0; + for (index = 0; index < count; ++index) + row[index] = erf(scratch[index]); + break; + case WF_OP_SINC: + for (index = 0; index < count; ++index) + scratch[index] = 3.14159265358979323846 * node->p0 + * (positions[index] - delay - shift); + wf_vector_sin(row, scratch, count); + for (index = 0; index < count; ++index) + row[index] = scratch[index] == 0.0 + ? 1.0 : row[index] / scratch[index]; + break; + case WF_OP_EXP: + for (index = 0; index < count; ++index) + scratch[index] = node->p0 + * (positions[index] - delay - shift); + wf_vector_exp(row, scratch, count); + break; + case WF_OP_COSH: + case WF_OP_SINH: + for (index = 0; index < count; ++index) + scratch[index] = node->p0 + * (positions[index] - delay - shift); + if (node->op == WF_OP_COSH) + wf_vector_cosh(row, scratch, count); + else + wf_vector_sinh(row, scratch, count); + break; case WF_OP_SQUARE: for (index = 0; index < count; ++index) { double position = positions[index] - delay; @@ -2625,6 +2854,26 @@ static int wf_evaluate_many_apple( row[index] = node->p0 * source[index]; break; } + case WF_OP_POWER: { + const double *source = matrix + (size_t)node->left * count; + wf_vector_power_scalar(row, source, node->p0, count); + if (lower != -DBL_MAX || upper != DBL_MAX) { + for (index = 0; index < count; ++index) { + double position = positions[index] - delay; + if (position < lower || position >= upper) row[index] = 0.0; + } + } + break; + } + case WF_OP_WINDOW: { + const double *source = matrix + (size_t)node->left * count; + for (index = 0; index < count; ++index) { + double position = positions[index] - delay; + row[index] = position >= lower && position < upper + ? source[index] : 0.0; + } + break; + } default: free(matrix); free(scratch); @@ -2633,6 +2882,8 @@ static int wf_evaluate_many_apple( } { const double *root = matrix + (size_t)wave->root * count; + /* Keep scaling and clipping fused in one compiler-vectorized pass. + * Two separate vDSP calls cost an extra full memory traversal. */ for (index = 0; index < count; ++index) { double value = scale * root[index]; if (value < lower_clip) value = lower_clip; @@ -2644,6 +2895,36 @@ static int wf_evaluate_many_apple( free(scratch); return 0; } + +static size_t wf_vector_chunk_count(const cwaveform_wave *wave, size_t count) { +#if defined(__APPLE__) + (void)wave; + return count; +#else + size_t chunk = (64u * 1024u) + / ((size_t)wave->node_count * sizeof(double)); + if (chunk < 256) chunk = 256; + if (chunk > 8192) chunk = 8192; + return chunk < count ? chunk : count; +#endif +} + +static int wf_evaluate_many_vector( + const cwaveform_wave *wave, const double *positions, size_t count, + double delay, double scale, double lower_clip, double upper_clip, + double *output) { + size_t cursor = 0; + size_t chunk = wf_vector_chunk_count(wave, count); + while (cursor < count) { + size_t batch = count - cursor < chunk ? count - cursor : chunk; + int status = wf_evaluate_chunk_vector( + wave, positions + cursor, batch, delay, scale, + lower_clip, upper_clip, output + cursor); + if (status != 0) return status; + cursor += batch; + } + return 0; +} #endif int cwaveform_wave_evaluate( @@ -2656,9 +2937,10 @@ int cwaveform_wave_evaluate( if (wave == NULL || positions == NULL || output == NULL || !isfinite(scale)) { return -1; } -#if defined(__APPLE__) +#if (defined(__APPLE__) || defined(WF_HAVE_X86_SIMD) || defined(_WIN32)) \ + && !defined(WF_DISABLE_BATCH_EVALUATOR) if (count >= 256) { - int status = wf_evaluate_many_apple( + int status = wf_evaluate_many_vector( wave, positions, count, (double)delay_tick / (double)wf_ticks_per_second, scale, lower_clip, upper_clip, output); @@ -2685,7 +2967,7 @@ static int16_t wf_quantize16(double value, double full_scale) { double scaled; if (value <= -full_scale) return INT16_MIN; if (value >= full_scale) return INT16_MAX; - scaled = nearbyint(value * (32768.0 / full_scale)); + scaled = round(value * (32768.0 / full_scale)); if (scaled <= -32768.0) return INT16_MIN; if (scaled >= 32767.0) return INT16_MAX; return (int16_t)scaled; @@ -2695,12 +2977,264 @@ static int32_t wf_quantize32(double value, double full_scale) { double scaled; if (value <= -full_scale) return INT32_MIN; if (value >= full_scale) return INT32_MAX; - scaled = nearbyint(value * (2147483648.0 / full_scale)); + scaled = round(value * (2147483648.0 / full_scale)); if (scaled <= -2147483648.0) return INT32_MIN; if (scaled >= 2147483647.0) return INT32_MAX; return (int32_t)scaled; } +static int wf_values_are_finite(const double *values, size_t count) { + size_t index; + for (index = 0; index < count; ++index) + if (!isfinite(values[index])) return 0; + return 1; +} + +#if defined(WF_HAVE_ARM64_NEON) +static int wf_quantize_array_neon(const double *values, size_t count, int dtype, + double full_scale, void *output) { + const float64x2_t finite_limit = vdupq_n_f64(DBL_MAX); + const float64x2_t multiplier = vdupq_n_f64( + dtype == CWAVEFORM_INT16 + ? 32768.0 / full_scale : 2147483648.0 / full_scale); + const float64x2_t minimum = vdupq_n_f64( + dtype == CWAVEFORM_INT16 ? -32768.0 : -2147483648.0); + const float64x2_t maximum = vdupq_n_f64( + dtype == CWAVEFORM_INT16 ? 32767.0 : 2147483647.0); + uint64x2_t finite = vdupq_n_u64(UINT64_MAX); + size_t index = 0; + for (; index + 4 <= count; index += 4) { + float64x2_t first = vld1q_f64(values + index); + float64x2_t second = vld1q_f64(values + index + 2); + int64x2_t first_integer; + int64x2_t second_integer; + finite = vandq_u64(finite, + vcleq_f64(vabsq_f64(first), finite_limit)); + finite = vandq_u64(finite, + vcleq_f64(vabsq_f64(second), finite_limit)); + first = vmulq_f64(first, multiplier); + second = vmulq_f64(second, multiplier); + first = vminq_f64(vmaxq_f64(first, minimum), maximum); + second = vminq_f64(vmaxq_f64(second, minimum), maximum); + first_integer = vcvtaq_s64_f64(first); + second_integer = vcvtaq_s64_f64(second); + if (dtype == CWAVEFORM_INT16) { + int32x4_t packed32 = vcombine_s32( + vmovn_s64(first_integer), vmovn_s64(second_integer)); + vst1_s16((int16_t *)output + index, vmovn_s32(packed32)); + } else { + vst1q_s32((int32_t *)output + index, vcombine_s32( + vmovn_s64(first_integer), vmovn_s64(second_integer))); + } + } + for (; index < count; ++index) { + if (!isfinite(values[index])) return -3; + if (dtype == CWAVEFORM_INT16) + ((int16_t *)output)[index] = wf_quantize16(values[index], full_scale); + else + ((int32_t *)output)[index] = wf_quantize32(values[index], full_scale); + } + return (vgetq_lane_u64(finite, 0) == UINT64_MAX + && vgetq_lane_u64(finite, 1) == UINT64_MAX) ? 0 : -3; +} +#endif + +#if defined(WF_HAVE_X86_SIMD) && !defined(WF_DISABLE_X86_SIMD) +static int wf_cpu_supports_avx2(void) { +#if defined(_MSC_VER) + int registers[4]; + unsigned __int64 xcr0; + __cpuid(registers, 0); + if (registers[0] < 7) return 0; + __cpuidex(registers, 1, 0); + if ((registers[2] & ((1 << 27) | (1 << 28))) + != ((1 << 27) | (1 << 28))) return 0; + xcr0 = _xgetbv(0); + if ((xcr0 & 0x6) != 0x6) return 0; + __cpuidex(registers, 7, 0); + return (registers[1] & (1 << 5)) != 0; +#else + return __builtin_cpu_supports("avx2"); +#endif +} + +#if !defined(WF_DISABLE_AVX512) && !defined(_MSC_VER) +static int wf_cpu_supports_avx512(void) { + return __builtin_cpu_supports("avx512f") + && __builtin_cpu_supports("avx512dq") + && __builtin_cpu_supports("avx512bw") + && __builtin_cpu_supports("avx512vl"); +} +#endif + +WF_TARGET_AVX2 +static int wf_quantize_array_avx2(const double *values, size_t count, int dtype, + double full_scale, void *output) { + const __m256d sign_mask = _mm256_set1_pd(-0.0); + const __m256d half = _mm256_set1_pd(0.5); + const __m256d finite_limit = _mm256_set1_pd(DBL_MAX); + const __m256d multiplier = _mm256_set1_pd( + dtype == CWAVEFORM_INT16 + ? 32768.0 / full_scale : 2147483648.0 / full_scale); + const __m256d minimum = _mm256_set1_pd( + dtype == CWAVEFORM_INT16 ? -32768.0 : -2147483648.0); + const __m256d maximum = _mm256_set1_pd( + dtype == CWAVEFORM_INT16 ? 32767.0 : 2147483647.0); + int finite = 1; + size_t index = 0; + for (; index + 4 <= count; index += 4) { + __m256d value = _mm256_loadu_pd(values + index); + __m256d absolute = _mm256_andnot_pd(sign_mask, value); + __m256d adjustment = _mm256_or_pd( + half, _mm256_and_pd(sign_mask, value)); + __m128i integer; + finite &= _mm256_movemask_pd(_mm256_cmp_pd( + absolute, finite_limit, _CMP_LE_OQ)) == 0xf; + value = _mm256_mul_pd(value, multiplier); + value = _mm256_min_pd(_mm256_max_pd(value, minimum), maximum); + integer = _mm256_cvttpd_epi32(_mm256_add_pd(value, adjustment)); + if (dtype == CWAVEFORM_INT16) { + __m128i packed = _mm_packs_epi32(integer, _mm_setzero_si128()); + _mm_storel_epi64((__m128i *)((int16_t *)output + index), packed); + } else { + _mm_storeu_si128((__m128i *)((int32_t *)output + index), integer); + } + } + for (; index < count; ++index) { + if (!isfinite(values[index])) return -3; + if (dtype == CWAVEFORM_INT16) + ((int16_t *)output)[index] = wf_quantize16(values[index], full_scale); + else + ((int32_t *)output)[index] = wf_quantize32(values[index], full_scale); + } + _mm256_zeroupper(); + return finite ? 0 : -3; +} + +#if !defined(WF_DISABLE_AVX512) && !defined(_MSC_VER) +WF_TARGET_AVX512 +static int wf_quantize_array_avx512( + const double *values, size_t count, int dtype, + double full_scale, void *output) { + const __m512d sign_mask = _mm512_set1_pd(-0.0); + const __m512d half = _mm512_set1_pd(0.5); + const __m512d finite_limit = _mm512_set1_pd(DBL_MAX); + const __m512d multiplier = _mm512_set1_pd( + dtype == CWAVEFORM_INT16 + ? 32768.0 / full_scale : 2147483648.0 / full_scale); + const __m512d minimum = _mm512_set1_pd( + dtype == CWAVEFORM_INT16 ? -32768.0 : -2147483648.0); + const __m512d maximum = _mm512_set1_pd( + dtype == CWAVEFORM_INT16 ? 32767.0 : 2147483647.0); + __mmask8 finite = (__mmask8)0xff; + size_t index = 0; + for (; index + 8 <= count; index += 8) { + __m512d value = _mm512_loadu_pd(values + index); + __m512d absolute = _mm512_andnot_pd(sign_mask, value); + __m512d adjustment = _mm512_or_pd( + half, _mm512_and_pd(sign_mask, value)); + __m512i integer; + __m256i packed32; + finite &= _mm512_cmp_pd_mask(absolute, finite_limit, _CMP_LE_OQ); + value = _mm512_mul_pd(value, multiplier); + value = _mm512_min_pd(_mm512_max_pd(value, minimum), maximum); + integer = _mm512_cvttpd_epi64(_mm512_add_pd(value, adjustment)); + packed32 = _mm512_cvtepi64_epi32(integer); + if (dtype == CWAVEFORM_INT16) { + __m128i packed16 = _mm256_cvtepi32_epi16(packed32); + _mm_storeu_si128( + (__m128i *)((int16_t *)output + index), packed16); + } else { + _mm256_storeu_si256( + (__m256i *)((int32_t *)output + index), packed32); + } + } + for (; index < count; ++index) { + if (!isfinite(values[index])) return -3; + if (dtype == CWAVEFORM_INT16) + ((int16_t *)output)[index] = wf_quantize16(values[index], full_scale); + else + ((int32_t *)output)[index] = wf_quantize32(values[index], full_scale); + } + return finite == (__mmask8)0xff ? 0 : -3; +} +#endif +#endif + +static int wf_quantize_array(const double *values, size_t count, int dtype, + double full_scale, void *output) { + size_t index; +#if defined(WF_HAVE_ARM64_NEON) + if (count >= 16) + return wf_quantize_array_neon( + values, count, dtype, full_scale, output); +#endif +#if defined(WF_HAVE_X86_SIMD) && !defined(WF_DISABLE_X86_SIMD) + if (count >= 16) { +#if !defined(WF_DISABLE_AVX512) && !defined(_MSC_VER) + if (wf_cpu_supports_avx512()) + return wf_quantize_array_avx512( + values, count, dtype, full_scale, output); +#endif + if (wf_cpu_supports_avx2()) + return wf_quantize_array_avx2( + values, count, dtype, full_scale, output); + } +#endif + if (!wf_values_are_finite(values, count)) return -3; +#if defined(__APPLE__) + if (count >= 256) { + double *scaled = (double *)malloc(count * sizeof(*scaled)); + double scale; + double minimum; + double maximum; + if (scaled == NULL) return -2; + if (dtype == CWAVEFORM_INT16) { + scale = 32768.0 / full_scale; + minimum = -32768.0; + maximum = 32767.0; + } else if (dtype == CWAVEFORM_INT32) { + scale = 2147483648.0 / full_scale; + minimum = -2147483648.0; + maximum = 2147483647.0; + } else { + free(scaled); + return -1; + } + vDSP_vsmulD(values, 1, &scale, scaled, 1, (vDSP_Length)count); + vDSP_vclipD(scaled, 1, &minimum, &maximum, scaled, 1, + (vDSP_Length)count); + if (dtype == CWAVEFORM_INT16) + vDSP_vfixr16D(scaled, 1, (int16_t *)output, 1, + (vDSP_Length)count); + else + vDSP_vfixr32D(scaled, 1, (int32_t *)output, 1, + (vDSP_Length)count); + free(scaled); + return 0; + } +#endif + if (dtype == CWAVEFORM_INT16) { + for (index = 0; index < count; ++index) + ((int16_t *)output)[index] = wf_quantize16(values[index], full_scale); + } else if (dtype == CWAVEFORM_INT32) { + for (index = 0; index < count; ++index) + ((int32_t *)output)[index] = wf_quantize32(values[index], full_scale); + } else { + return -1; + } + return 0; +} + +int cwaveform_quantize(const double *values, size_t count, int dtype, + double full_scale, void *output) { + if ((count != 0 && (values == NULL || output == NULL)) + || (dtype != CWAVEFORM_INT16 && dtype != CWAVEFORM_INT32) + || !isfinite(full_scale) || full_scale <= 0.0) return -1; + if (count == 0) return 0; + return wf_quantize_array(values, count, dtype, full_scale, output); +} + /* Form the local coordinate before converting to binary64. Subtracting two * already-rounded seconds values loses enough precision to move an exact * pulse boundary outside its half-open support. */ @@ -2729,6 +3263,53 @@ int cwaveform_wave_sample( || !isfinite(full_scale) || full_scale <= 0.0) { return -1; } + if (dtype != CWAVEFORM_FLOAT64 && dtype != CWAVEFORM_INT16 + && dtype != CWAVEFORM_INT32) return -1; +#if (defined(__APPLE__) || defined(WF_HAVE_X86_SIMD) || defined(_WIN32)) \ + && !defined(WF_DISABLE_BATCH_EVALUATOR) + if (count >= 256) { + size_t cursor = 0; + size_t capacity = wf_vector_chunk_count(wave, count); + size_t item_size = dtype == CWAVEFORM_INT16 + ? sizeof(int16_t) : dtype == CWAVEFORM_INT32 + ? sizeof(int32_t) : sizeof(double); + double *positions = (double *)malloc(capacity * sizeof(*positions)); + double *samples = dtype == CWAVEFORM_FLOAT64 + ? NULL : (double *)malloc(capacity * sizeof(*samples)); + if (positions != NULL + && (dtype == CWAVEFORM_FLOAT64 || samples != NULL)) { + int status = 0; + while (cursor < count) { + size_t batch = count - cursor < capacity + ? count - cursor : capacity; + double *chunk_samples = dtype == CWAVEFORM_FLOAT64 + ? (double *)output + cursor : samples; + for (index = 0; index < batch; ++index) + positions[index] = wf_local_grid_position( + start_tick, cursor + index, + step_numerator, step_denominator, delay_tick); + status = wf_evaluate_many_vector( + wave, positions, batch, 0.0, scale, + lower_clip, upper_clip, chunk_samples); + if (status != 0) break; + if (dtype == CWAVEFORM_FLOAT64) + status = wf_values_are_finite(chunk_samples, batch) + ? 0 : -3; + else + status = wf_quantize_array( + chunk_samples, batch, dtype, full_scale, + (uint8_t *)output + cursor * item_size); + if (status != 0) break; + cursor += batch; + } + free(positions); + free(samples); + return status; + } + free(positions); + free(samples); + } +#endif values = (double *)malloc((size_t)wave->node_count * sizeof(*values)); if (values == NULL) return -2; for (index = 0; index < count; ++index) { @@ -2747,9 +3328,6 @@ int cwaveform_wave_sample( ((int16_t *)output)[index] = wf_quantize16(value, full_scale); } else if (dtype == CWAVEFORM_INT32) { ((int32_t *)output)[index] = wf_quantize32(value, full_scale); - } else { - free(values); - return -1; } } free(values); @@ -3723,16 +4301,13 @@ int cwaveform_sample_plan_sample( if (!plan->non_overlapping) { double *values = plan->count == 0 ? NULL : (double *)malloc(plan->count * sizeof(double)); + int status; if (plan->count != 0 && values == NULL) return -2; wf_plan_sample_float(plan, offset, values); - for (index = 0; index < plan->count; ++index) { - if (dtype == CWAVEFORM_INT16) - ((int16_t *)output)[index] = wf_quantize16(values[index], full_scale); - else - ((int32_t *)output)[index] = wf_quantize32(values[index], full_scale); - } + status = wf_quantize_array( + values, plan->count, dtype, full_scale, output); free(values); - return 0; + return status; } if (dtype == CWAVEFORM_INT16) { int16_t base = wf_quantize16(offset, full_scale); @@ -3833,12 +4408,14 @@ int cwaveform_stack_sample( * wf_evaluate_one(wave, position, values); } } - if (dtype == CWAVEFORM_INT16) { - for (index = 0; index < count; ++index) - ((int16_t *)output)[index] = wf_quantize16(float_output[index], full_scale); - } else if (dtype == CWAVEFORM_INT32) { - for (index = 0; index < count; ++index) - ((int32_t *)output)[index] = wf_quantize32(float_output[index], full_scale); + if (dtype == CWAVEFORM_INT16 || dtype == CWAVEFORM_INT32) { + int status = wf_quantize_array( + float_output, count, dtype, full_scale, output); + if (status != 0) { + free(float_output); + free(values); + return status; + } } else if (dtype != CWAVEFORM_FLOAT64) { if (float_output != output) free(float_output); free(values); diff --git a/waveforms/_cwaveform.h b/waveforms/_cwaveform.h index bd75058..9d9041f 100644 --- a/waveforms/_cwaveform.h +++ b/waveforms/_cwaveform.h @@ -118,6 +118,9 @@ CWAVEFORM_API int cwaveform_wave_sample( int64_t step_numerator, int64_t step_denominator, int64_t delay_tick, double scale, double lower_clip, double upper_clip, int dtype, double full_scale, void *output); +CWAVEFORM_API int cwaveform_quantize( + const double *values, size_t count, int dtype, + double full_scale, void *output); CWAVEFORM_API cwaveform_stack *cwaveform_stack_create( cwaveform_wave *const *templates, const uint32_t *template_ids, diff --git a/waveforms/_waveform.pyx b/waveforms/_waveform.pyx index 449750c..13ed5cd 100644 --- a/waveforms/_waveform.pyx +++ b/waveforms/_waveform.pyx @@ -134,26 +134,22 @@ def quantize_samples(values, bits, full_scale=1.0, out=None): """Saturating real-signal quantizer for signed 16- and 32-bit DAC data.""" cdef object source cdef object target - cdef double[::1] source_view - cdef int16_t[::1] target16 - cdef int32_t[::1] target32 - cdef Py_ssize_t index, size - cdef double value, scaled - cdef double scale - cdef double maximum - cdef double minimum + cdef Py_ssize_t size cdef double full_scale_value + cdef const double *source_pointer + cdef void *target_pointer + cdef int status + cdef int bit_count full_scale_value = float(full_scale) if not np.isfinite(full_scale_value) or full_scale_value <= 0: raise ValueError("full_scale must be a finite positive number") if bits not in (16, 32): raise ValueError("bits must be 16 or 32") + bit_count = bits if np.iscomplexobj(values): raise TypeError("integer quantization requires a real signal") source = np.ascontiguousarray(values, dtype=np.float64) - if not np.all(np.isfinite(source)): - raise ValueError("cannot quantize non-finite samples") dtype = np.dtype(np.int16 if bits == 16 else np.int32) if out is None: target = np.empty(source.shape, dtype=dtype) @@ -166,46 +162,21 @@ def quantize_samples(values, bits, full_scale=1.0, out=None): if not target.flags.c_contiguous or not target.flags.writeable: raise ValueError("out must be a writable C-contiguous array") - source_view = source.reshape(-1) - size = source_view.shape[0] - if bits == 16: - target16 = target.reshape(-1) - scale = 32768.0 / full_scale_value - minimum = -32768.0 - maximum = 32767.0 - with nogil: - for index in range(size): - value = source_view[index] - if value <= -full_scale_value: - target16[index] = -32768 - elif value >= full_scale_value: - target16[index] = 32767 - else: - scaled = c_round(value * scale) - if scaled < minimum: - scaled = minimum - elif scaled > maximum: - scaled = maximum - target16[index] = scaled - else: - target32 = target.reshape(-1) - scale = 2147483648.0 / full_scale_value - minimum = -2147483648.0 - maximum = 2147483647.0 - with nogil: - for index in range(size): - value = source_view[index] - if value <= -full_scale_value: - target32[index] = -2147483648 - elif value >= full_scale_value: - target32[index] = 2147483647 - else: - scaled = c_round(value * scale) - if scaled < minimum: - scaled = minimum - elif scaled > maximum: - scaled = maximum - target32[index] = scaled + size = source.size + if size == 0: + return target + source_pointer = source.ctypes.data + target_pointer = target.ctypes.data + with nogil: + status = cwaveform_quantize( + source_pointer, size, bit_count, full_scale_value, target_pointer, + ) + if status == -3: + raise ValueError("cannot quantize non-finite samples") + if status == -2: + raise MemoryError("C quantization allocation failed") + if status != 0: + raise RuntimeError(f"C quantization failed with status {status}") return target @@ -525,6 +496,8 @@ cdef extern from "_cwaveform.h": int cwaveform_wave_sample( const cwaveform_wave *, int64_t, size_t, int64_t, int64_t, int64_t, double, double, double, int, double, void *) noexcept nogil + int cwaveform_quantize( + const double *, size_t, int, double, void *) noexcept nogil cwaveform_stack *cwaveform_stack_create( cwaveform_wave *const *, const uint32_t *, const int64_t *, From 95b7a38bdbb69335ec418e0b925053e371975fb3 Mon Sep 17 00:00:00 2001 From: feihoo87 Date: Sat, 29 Aug 2026 11:53:42 +0800 Subject: [PATCH 3/4] Make native SIMD benchmark portable to Windows --- benchmarks/benchmark_simd_c.c | 15 +++++++++++++++ 1 file changed, 15 insertions(+) diff --git a/benchmarks/benchmark_simd_c.c b/benchmarks/benchmark_simd_c.c index e32b60a..e8dcb9c 100644 --- a/benchmarks/benchmark_simd_c.c +++ b/benchmarks/benchmark_simd_c.c @@ -1,4 +1,6 @@ +#if !defined(_WIN32) #define _POSIX_C_SOURCE 200809L +#endif #include "../waveforms/_cwaveform.h" @@ -7,12 +9,25 @@ #include #include #include +#if defined(_WIN32) +#define WIN32_LEAN_AND_MEAN +#include +#else #include +#endif static double now_seconds(void) { +#if defined(_WIN32) + LARGE_INTEGER counter; + LARGE_INTEGER frequency; + QueryPerformanceFrequency(&frequency); + QueryPerformanceCounter(&counter); + return (double)counter.QuadPart / (double)frequency.QuadPart; +#else struct timespec value; clock_gettime(CLOCK_MONOTONIC, &value); return (double)value.tv_sec + 1e-9 * (double)value.tv_nsec; +#endif } static double best_evaluate(const cwaveform_wave *wave, const double *positions, From 6e867937a6c0995921a4a7986894232dafabbf58 Mon Sep 17 00:00:00 2001 From: feihoo87 Date: Sat, 29 Aug 2026 12:04:08 +0800 Subject: [PATCH 4/4] Release 3.3.2 --- waveforms/version.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/waveforms/version.py b/waveforms/version.py index 6f292e0..c52555c 100644 --- a/waveforms/version.py +++ b/waveforms/version.py @@ -1,2 +1,2 @@ """Define version number here and read it from setup.py automatically""" -__version__ = "3.3.1" +__version__ = "3.3.2"