internal.h

浏览该文件的文档.

#ifndef LMMC_INTERNAL_H
#define LMMC_INTERNAL_H
 
#include "lmmc/config.h"
#include <stddef.h>
#include <stdint.h>
#include "lmmc/dense.h"
#include <math.h>

static inline double lmmc_two_sum(double a, double b, double* residual)
{
    const double sum = a + b;
    const double virtual_b = sum - a;
    *residual = (a - (sum - virtual_b)) + (b - virtual_b);
    return sum;
}

static inline void lmmc_compensated_add(double value, double* sum, double* correction)
{
    double residual;
    *sum = lmmc_two_sum(*sum, value, &residual);
    *correction += residual;
}
 
typedef struct lmmc_rng_t lmmc_rng_t;

lmmc_rng_t* lmmc_rng_default_get(void);

void lmmc_rng_default_reset(void);
 
static inline void lmmc_fill_identity_unchecked(lmmc_mat_t* matrix)
{
    lmmc_real_t zero = 0.0;
    lmmc_real_t one = 1.0;
    for (size_t row = 0; row < matrix->rows; ++row) {
        for (size_t column = 0; column < matrix->cols; ++column) {
            LMMC_REAL_SET(&matrix->data[row * matrix->stride + column],
                          row == column ? &one : &zero);
        }
    }
}
 
static inline int lmmc_is_power_of_four(
    size_t value, unsigned* out_digits)
{
    size_t remaining = value;
    unsigned digits = 0;
    if (value == 0) {
        return 0;
    }
    while ((remaining & (size_t)3) == 0) {
        remaining >>= 2;
        ++digits;
    }
    if (remaining != (size_t)1) {
        return 0;
    }
    if (out_digits != NULL) {
        *out_digits = digits;
    }
    return 1;
}
static inline int lmmc_safe_mul_size(size_t a, size_t b, size_t *result)
{
    if (a == 0 || b == 0) {
        *result = 0;
        return 1;
    }
    if (a > SIZE_MAX / b) {
        return 0;
    }
    *result = a * b;
    return 1;
}

static inline int lmmc_safe_add_size(size_t a, size_t b, size_t *result)
{
    if (a > SIZE_MAX - b) {
        return 0;
    }
    *result = a + b;
    return 1;
}
 
 
typedef struct lmmc_storage_envelope {
    uintptr_t begin;
    uintptr_t end;
} lmmc_storage_envelope_t;

static inline int lmmc_storage_envelope_checked(
    const void *data,
    size_t rows,
    size_t cols,
    size_t stride,
    size_t element_size,
    lmmc_storage_envelope_t *out)
{
    size_t last_row;
    size_t elements;
    size_t bytes;
    uintptr_t begin;
 
    if (data == NULL || out == NULL) {
        return 0;
    }
    begin = (uintptr_t)data;
    out->begin = begin;
    out->end = begin;
    if (rows == 0 || cols == 0) {
        return 1;
    }
    if (!lmmc_safe_mul_size(rows - 1, stride, &last_row)) {
        return 0;
    }
    if (!lmmc_safe_add_size(last_row, cols, &elements)) {
        return 0;
    }
    if (!lmmc_safe_mul_size(elements, element_size, &bytes)) {
        return 0;
    }
    if ((uintptr_t)bytes > UINTPTR_MAX - begin) {
        return 0;
    }
    out->end = begin + (uintptr_t)bytes;
    return 1;
}

static inline int lmmc_mat_descriptor_is_valid(const lmmc_mat_t *matrix)
{
    lmmc_storage_envelope_t envelope;
    return matrix != NULL && matrix->data != NULL &&
           matrix->rows != 0 && matrix->cols != 0 &&
           matrix->stride >= matrix->cols &&
           lmmc_storage_envelope_checked(
               matrix->data, matrix->rows, matrix->cols, matrix->stride,
               sizeof(lmmc_real_t), &envelope);
}

static inline int lmmc_vec_descriptor_is_valid(const lmmc_vec_t *vector)
{
    lmmc_storage_envelope_t envelope;
    return vector != NULL && vector->data != NULL && vector->size != 0 &&
           lmmc_storage_envelope_checked(
               vector->data, 1, vector->size, vector->size,
               sizeof(lmmc_real_t), &envelope);
}

static inline int lmmc_storage_envelopes_overlap(
    const lmmc_storage_envelope_t *left,
    const lmmc_storage_envelope_t *right)
{
    return left->begin < right->end && right->begin < left->end;
}

static inline int lmmc_is_finite(const lmmc_real_t *x)
{
    return LMMC_REAL_IS_FINITE(x) ? 1 : 0;
}
 
static inline int lmmc_real_range_is_finite(
    const lmmc_real_t* values, size_t count)
{
    if (values == NULL) {
        return 0;
    }
    for (size_t index = 0; index < count; ++index) {
        if (!lmmc_is_finite(&values[index])) {
            return 0;
        }
    }
    return 1;
}
 
static inline void lmmc_interval_center_half_width(
    lmmc_real_t left, lmmc_real_t right,
    lmmc_real_t* center, lmmc_real_t* half_width)
{
    if ((left < 0.0 && right > 0.0) ||
        (left > 0.0 && right < 0.0)) {
        *center = left * 0.5 + right * 0.5;
        *half_width = right * 0.5 - left * 0.5;
        return;
    }
    *half_width = (right - left) * 0.5;
    *center = left + *half_width;
}
 
static inline lmmc_real_t lmmc_interval_midpoint(
    lmmc_real_t left, lmmc_real_t right)
{
    lmmc_real_t center;
    lmmc_real_t half_width;
    lmmc_interval_center_half_width(
        left, right, &center, &half_width);
    return center;
}

typedef struct lmmc_scaled_sumsq {
    lmmc_real_t scale;
    lmmc_real_t sumsq;
} lmmc_scaled_sumsq_t;
 
static inline void lmmc_scaled_sumsq_init(lmmc_scaled_sumsq_t *acc)
{
    acc->scale = 0.0;
    acc->sumsq = 1.0;
}
 
static inline void lmmc_scaled_sumsq_add(
    lmmc_scaled_sumsq_t *acc, lmmc_real_t value)
{
    lmmc_real_t magnitude = fabs(value);
    if (isnan(acc->scale)) {
        return;
    }
    if (isnan(magnitude)) {
        acc->scale = magnitude;
        acc->sumsq = 1.0;
        return;
    }
    if (magnitude == 0.0 || isinf(acc->scale)) {
        return;
    }
    if (isinf(magnitude)) {
        acc->scale = magnitude;
        acc->sumsq = 1.0;
        return;
    }
    if (acc->scale < magnitude) {
        lmmc_real_t ratio = acc->scale / magnitude;
        acc->sumsq = 1.0 + acc->sumsq * ratio * ratio;
        acc->scale = magnitude;
    } else {
        lmmc_real_t ratio = magnitude / acc->scale;
        acc->sumsq += ratio * ratio;
    }
}
 
static inline lmmc_real_t lmmc_scaled_sumsq_norm(
    const lmmc_scaled_sumsq_t *acc)
{
    return acc->scale == 0.0 ? 0.0 : acc->scale * sqrt(acc->sumsq);
}

static inline lmmc_real_t lmmc_abs(lmmc_real_t x)
{
    return fabs(x);
}

static inline lmmc_real_t lmmc_max(lmmc_real_t a, lmmc_real_t b)
{
    return (a >= b) ? a : b;
}

static inline lmmc_real_t lmmc_clamp(lmmc_real_t x, lmmc_real_t lo, lmmc_real_t hi)
{
    if (x < lo) {
        return lo;
    }
    if (x > hi) {
        return hi;
    }
    return x;
}

static inline void lmmc_swap(lmmc_real_t *a, lmmc_real_t *b)
{
    lmmc_real_t tmp = *a;
    *a = *b;
    *b = tmp;
}
 
#endif