0
0
Fork 0
mirror of https://github.com/matrix-construct/construct synced 2024-12-30 17:34:04 +01:00
construct/include/ircd/math/norm.h
2021-03-17 18:55:56 -07:00

102 lines
2 KiB
C++

// Matrix Construct
//
// Copyright (C) Matrix Construct Developers, Authors & Contributors
// Copyright (C) 2016-2021 Jason Volk <jason@zemos.net>
//
// Permission to use, copy, modify, and/or distribute this software for any
// purpose with or without fee is hereby granted, provided that the above
// copyright notice and this permission notice is present in all copies. The
// full license for this software is available in the LICENSE file.
#pragma once
#define HAVE_IRCD_MATH_NORM_H
namespace ircd::math
{
template<class T,
class D>
void
norm(const vector_view<T> &out,
const vector_view<const T> &in,
const simd::lane_type<T> &eps,
const vector_view<D> &tmp);
}
/// Renormalization.
///
/// TODO: Figure out how to use less double
template<class T,
class D>
inline void
ircd::math::norm(const vector_view<T> &__restrict__ out,
const vector_view<const T> &__restrict__ in,
const simd::lane_type<T> &epsilon,
const vector_view<D> &tmp)
{
using singlep = simd::lane_type<T>;
using doublep = simd::lane_type<D>;
static_assert(simd::lanes<T>() == simd::lanes<D>());
assert(out.size() >= in.size());
assert(tmp.size() >= in.size());
const size_t num
{
in.size()
};
const size_t elems
{
num * simd::lanes<T>()
};
const doublep mean
(
//TODO: XXX
math::mean<singlep, doublep>(vector_view<const singlep>
{
reinterpret_cast<const singlep *>(in.data()), elems
})
);
for(uint i(0); i < num; ++i)
{
const auto wider
{
lane_cast<D>(in[i])
};
tmp[i] = wider - mean;
const auto res
{
pow(tmp[i], 2)
};
out[i] = lane_cast<T>(res);
}
const doublep s
(
//TODO: XXX
math::mean<singlep, doublep>(vector_view<const singlep>
{
reinterpret_cast<const singlep *>(out.data()), elems
})
);
const doublep divisor
(
sqrt(s + epsilon)
);
for(uint i(0); i < num; ++i)
{
const auto res
{
tmp[i] / divisor
};
out[i] = lane_cast<T>(res);
}
}