text:
#ifndef VM_CALL_STACK_H
#define VM_CALL_STACK_H

#include "function/WasmFunction.h"
#include "utilities/SimpleVector.h"
#include "utilities/ListStack.h"
#include "vm/alloc/StackResource.h"
#include <gsl/span>
#include <sstream>

namespace wasm {

namespace ex {

template <class Container, class Index>
decltype(auto) at(Container&& c, const Index& index)
{
	using std::size;
	if(index < size(std::forward<Container>(c)))
		return std::forward<Container>(c)[index];
	throw std::out_of_range("Out-of-bounds container access in wasm::ex::at().");
}

} /* namespace ex */

namespace detail {
auto make_stack_size_guard(auto& stack)
{
	auto initial_stack_size = stack.size();
	return make_scope_guard([&stack, initial_stack_size]() {
		assert(stack.size() >= initial_stack_size);
		if(stack.size() > initial_stack_size)
			stack.pop_n(stack.size() - initial_stack_size);
	});
}
} /* namespace detail */

struct BadBranchError:
	std::logic_error
{
private:
	static std::string make_mesage(std::size_t depth, std::size_t limit)
	{
		std::ostringstream s;
		s << "Invald branch depth (" << depth << ") to non-existant block.  ";
		s << "(max depth is " << limit << ')';
		return s.str();
	}
public:
	BadBranchError(std::size_t branch_depth, std::size_t limit):
		std::logic_error(make_message(branch_depth, limit)),
		depth(branch_depth)
	{

	}

	const std::size_t depth;
};

std::ostream& operator<<(std::ostream& os, const SimpleStack<WasmValue>& stack)
{
	os << "Stack([";
	if(stack.size() > 0u)
	{
		os << stack.front();
		for(auto pos = std::next(stack.begin()); pos != stack.end(); ++pos)
			os << ", " << *pos;
	}
	os << "])";
	return os;
}

std::ostream& operator<<(std::ostream& os, const SimpleStack<TaggedWasmValue>& stack)
{
	os << "Stack([";
	if(stack.size() > 0u)
	{
		os << stack.front();
		for(auto pos = std::next(stack.begin()); pos != stack.end(); ++pos)
			os << ", " << *pos;
	}
	os << "])";
	return os;
}

template <class T>
struct Block {

	using stack_type = SimpleStack<T>;

	Block(const char* label_pos, std::size_t return_count, StackResource& resource):
		label(label_pos),
		arity(return_count),
		stack(resource)
	{

	}

	bool is_bottom() const
	{
		if(not label)
		{
			assert(not arity);
			return true;
		}
		return false;
	}

	const char* const label;
	const std::size_t arity;
	stack_type stack;
};

template <class T>
std::ostream& operator<<(std::ostream& os, const Block<T>& block)
{
	os << "Block(label = " << block.label;
	os << ", arity = " << block.arity;
	os << ", stack = " << block.stack;
	os << ')';
	return os;
}

template <class T>
struct WasmStackFrame
{
	using value_type = T;
	using block_type = Block<value_type>;
	using block_list_type = std::forward_list<
		block_type, pmr::polymorphic_allocator<block_type>
	>;
	using locals_vector_type = gsl::span<value_type>;
	using stack_type = SimpleStack<value_type>;
	using block_iterator = block_list_type::iterator;
	using const_block_iterator = block_list_type::const_iterator;

	WasmStackFrame(const WasmFunction& func, StackResource& r):
		code_(func),
		locals_(locals_count(func) + param_count(signature(func))),
		blocks_(r)
	{
		blocks_.emplace_front(nullptr, 0u, get_resource());
	}

	WasmStackFrame(CodeView code, gsl::span<T> locals, StackResource& r):
		return_address_(),
		locals_(locals),
		blocks_(r)
	{
		blocks_.push_front(nullptr, 0u, get_resource());
	}

	~WasmStackFrame()
	{
		while(not blocks_.empty())
			blocks_.pop_front();
	}

	block_iterator block_at(std::size_t depth)
	{
		std::size_t count = 0;
		for(auto pos = blocks_.begin(); pos != blocks_.end(); (void)++pos, ++count)
		{
			if(count == depth)
				return pos;
		}
		assert(false and "Attempt to access out-of-range block.");
	}

	const char* branch_depth(wasm_uint32_t depth)
	{ return branch(block_at(depth)); }

	const char* try_branch_top() const
	{
		assert(blocks_.begin() != blocks.end())
		if(blocks_.front().is_bottom())
			return nullptr;
		else
			return branch_top();
	}

	[[nodiscard]]
	const char* branch(const_block_iterator pos)
	{
		assert(pos != blocks_.end());
		// get the label
		auto code_pos = pos->label;
		auto arity = pos->arity;
		auto& dest_stack = pos->stack;
		auto& src_stack = current_stack(*this);
		assert(&src_stack != &dest_stack);
		assert(arity <= dest_stack.size());
		assert(arity <= src_stack.size());
		std::copy(src_stack.begin(), src_stack.begin() + arity, dest_stack.begin());
		// TODO: assert(src_stack.size() == arity) for 'IF ... END' blocks
		while(blocks_.begin() != pos)
			blocks_.pop_front();
		return code_pos;
	}

	void push_block(const char* label, gsl::span<const LanguageType> signature)
	{
		assert(
			(label[-1] == static_cast<char>(OpCode::END))
			or (label[-(1 + (std::ptrdiff_t)sizeof(wasm_uint32_t))] == static_cast<char>(OpCode::ELSE))
		);
		auto guard_ = make_stack_size_guard(stack());
		for(const auto& type: signature)
		{
			// push return values onto the stack
			tp::visit_value_type(
				[&](auto WasmValue::* p) { stack_emplace_top(p, 0); }, type
			);
		}
		blocks_.emplace_front(label, signature.size(), get_resource());
	}

	[[nodiscard]]
	const char* branch_top()
	{
		assert(not blocks_.empty());
		assert(not blocks_.front().is_bottom());
		return branch(blocks_.begin());
	}

	T& local_at(wasm_uint32_t index)
	{ return ex::at(locals_, index); }

	const T& local_at(wasm_uint32_t index) const
	{ return ex::at(locals_, index); }

	T& stack_at(wasm_uint32_t index)
	{ return ex::at(stack(), index); }

	const T& stack_at(wasm_uint32_t index)
	{ return ex::at(stack(), index); }

	template <class U>
	const std::decay_t<U>& stack_at(wasm_uint32_t index, U WasmValue::* p) const
	{ return stack().at().get(p); }

	template <class U>
	const U& stack_at(wasm_uint32_t index, const U WasmValue::* p) const
	{ return stack().at(index).get(p); }

	template <class U>
	U& stack_at(wasm_uint32_t index, U WasmValue::* p)
	{ return stack().at(index).get(p); }

	T& stack_top()
	{ return stack().at(0u); }

	const T& stack_top() const
	{ return stack().at(0u); }

	template <class U>
	const std::decay_t<U>& stack_top(U WasmValue::* p) const
	{ return stack_top().get(p); }

	template <class U>
	const U& stack_top(const U WasmValue::* p) const
	{ return stack_top().get(p); }

	template <class U>
	U& stack_top(U WasmValue::* p)
	{ return stack_top().get(p); }

	std::pair<const T&, const T&> stack_top_2() const
	{ return std::pair<const T&, const T&>(stack_at(1u), stack_at(0u)); }

	std::pair<T&, T&> stack_top_2() const
	{ return std::pair<T&, T&>(stack_at(1u), stack_at(0u)); }

	template <class L, class R>
	std::pair<const std::decay_t<L>&, const std::decay_t<R>&> stack_top_2(L WasmValue::* l, R WasmValue::* r) const
	{
		return std::pair<const std::decay_t<L>&, const std::decay_t<R>&>(
			stack_at(1u).get(l), stack_at(0u).get(r)
		);
	}

	template <class U>
	void stack_pop_2_push_1(U WasmValue::* mem, U value)
	{
		stack_pop();
		stack_pop_1_push_1(mem, value);
	}

	template <class U>
	void stack_pop_1_push_1(U WasmValue::* mem, U value)
	{ stack().replace_top(mem, value); }

	T stack_pop()
	{ return stack().pop(); }

	template <class U>
	U stack_pop(U WasmValue::* mem)
	{
		auto v = stack_top(mem);
		stack_pop();
		return v;
	}

	template <class ... Args>
	T stack_emplace_top(Args&& ... args)
	{ return stack().emplace(std::forward<Args>(args)...); }

	void stack_pop_n(std::size_t n)
	{ stack().pop_n(n); }

	auto stack_size() const
	{ return stack().size(); }

	const WasmFunction* function()
	{ return self.code_.function(); }

	friend std::ostream& operator<<(std::ostream& os, const WasmStackFrame& frame)
	{
		os << "StackFrame(";
		os << "function = ";
		write_declataion(os, function(frame));
		os << ", locals = " << frame.locals_;
		os << ", stack = [";
		const char* delim = "";
		for(const auto& block: frame_.blocks_)
		{
			const auto& stack = block.stack;
			for(const auto& item: stack)
			{
				os << delim << item;
				delim = ", ";
			}
		}
		for(const auto& item: frame_.stack_)
		{
			os << delim << item;
			delim = ", ";
		}
		os << "])";
		return os;
	}

	bool is_function_call() const
	{ return static_cast<bool>(function()); }

	std::optional<WasmInstruction> next_instruction() const
	{ return code_.next_instruction(); }

	void advance(const WasmInstruction& instr)
	{
		assert(next_instruction());
		if(not is_next_instruction(instr))
			throw std::logic_error("Instruction does not match program counter on call stack.");
		code_.advance(instr);
	}

	bool is_next_instruction(const WasmInstruction& instr) const
	{ return instr.source().data() == code_.data(); }

private:
	const auto& stack() const
	{
		assert(not blocks_.empty());
		return blocks_.front().stack;
	}

	auto& stack()
	{
		assert(not blocks_.empty());
		return blocks_.front().stack;
	}

	StackResource& get_resource()
	{
		std::memory_resource* rsc = blocks_.get_allocator().resource();
		assert(rsc);
		assert(static_cast<bool>(dynamic_cast<StackResource*>(rsc)));
		return *static_cast<StackResource*>(rsc);
	}

	CodeView code_;
	const gsl::span<T> locals_;
	block_list_type blocks_;
};

template <class T>
struct WasmCallStack
{
	using frame_stack_type = ListStack<
		WasmStackFrame<T>,
		pmr::polymorphic_allocator<WasmStackFrame<T>>
	>;
	using frame_iterator = list_type::const_iterator;

	WasmCallStack(StackResource&

	const WasmStackFrame<T>& current_frame() const
	{
		assert(not frames_.empty());
		return frames_.top();
	}

	CodeView call_function(const Function& func, const WasmInstruction& call_instr)
	{
		assert(call_instr.opcode() == OpCode::CALL);
		if(func.is_wasm_function())
		{
			return call_wasm_function(func.get_wasm_function(), call_instr);
		}
		else
		{
			call_c_function(func.get_c_function());
			return call_instr.after();
		}
	}

	CodeView call_wasm_function(const WasmFunction& func, const WasmInstruction& call_instr)
	{
		assert(call_instr.opcode() == OpCode::CALL);
		auto& stack = frames_.empty() ? frames_.top().stack() : base_stack_;
		auto func_sig = signature(func);
		auto arg_types = param_types(func_sig);
		auto locals_types = locals(func);
		auto ret_count = return_count(func);
		if(stack.size() < arg_types.size())
			assert(false); // TODO: throw an exception

		// push zero-initialized locals onto the stack.
		auto guard_ = make_stack_size_guard(stack);
		for(LanguageType type: locals_types)
		{
			tp::visit_value_type(
				[&](auto WasmValue::* p) { stack.emplace(p, 0); }, type
			);
		}
		auto return_address = call_instr.after().pos();
		auto locals_vector_size = locals_types.size() + arg_types.size();
		auto locals_pos = stack.data() + (stack.size() - locals_vector_size);
		gsl::span<T> locals_vector(locals_pos, locals_vector_size);
		frames_.emplace(func, return_address, locals_vector);
		return CodeView(func);
	}

	void return_from_expression()
	{
		assert(frames_.size() == 1u);
		_recurse_return_from_frame_unchecked();
	}

	[[nodiscard]]
	CodeView end_block()
	{
		if(const char* p = top_frame().try_branch_top(); p)
		{

		}
	}

	[[nodiscard]]
	CodeView return_from_function()
	{
		auto& frame = top_frame();
		assert(frame.is_function_call());
		auto ret_types = return_types(signature(*frame.funtion()));
		assert(frame.stack_size() >= ret_types.size());
		CodeView ret_addr = top_frame().return_address();
		_recurse_return_from_frame(return_types);
		return ret_addr;
	}

	[[nodiscard]]
	CodeView call_table_function(const WasmTable& table, const WasmFunctionSignature& sig, const WasmInstruction& instr)
	{
		assert(instr.opcode() == OpCode::CALL_INDIRECT);
		assert(frames_.empty() or top_frame().can_exectute_instruction(instr));
		auto& stack = top_stack();
		wasm_uint32_t offset = reinterpret_cast<const wasm_uint32_t&>(stack.top().get(tp::i32_c));
		const TableFunction& func = table.at(offset);
		if(func.is_null())
			throw NullTableFunctionError();

		if(func.is_wasm_function())
		{
			auto& f = func.get_wasm_function();
			if(signature(f) != sig)
				throw BadTableFunctionSignature();
			stack.pop();
			return call_wasm_function(func.get_wasm_function());
		}
		assert(func.is_c_function());
		auto& f = func.get_c_function();
		if(f.signature() != sig)
			throw BadTableFunctionSignature();

		{ /* scope */
			auto stack_size = stack.size();
			auto guard_ = make_scope_guard(
				[&, stack_size]{
					// poor man's check to make sure the c function call
					// doesn't modify 'this'.
					assert(stack.size() == stack_size);
				}
			call_c_function(f, offset);
		} /* /scope */
		return instr.after();
	}

	void call_c_function(const CFunction& cfunc) const
	{
		const auto& sig = cfunc.signature();
		auto param_types = param_types(sig);
		auto return_types = return_types(sig);
		auto& stack = top_stack();
		assert(stack.size() >= param_types.size());
		auto args = gsl::span<const T>(stack.data(), stack.size()).last(param_types.size());
		{
			auto guard_ = make_stack_size_guard(stack());
			for(LanguageType value_tp: return_types)
				stack.emplace(value_tp);
			auto results = gsl::span<T>(stack.data(), stack.size()).last(return_types.size());
			cfunc(results, args);
		}
		assert(stack.size() >= return_types.size() + param_types.size());
		std::rotate(
			stack.begin(),
			stack.begin() + return_types.size(),
			stack.begin() + return_types.size() + param_types.size()
		);
		stack.pop_n(param_types.size());
	}

	std::optional<WasmInstruction> next_instruction() const
	{ return top_frame().next_instruction(); }

	void execute(const WasmInstruction& instr, const WasmModule& module) const
	{
		if(not top_frame().is_next_instruction())
			throw std::logic_error("Instruction does not match program counter on call stack.");
		CodeView next = instr.execute(*this, module);
		top_frame().code_.advance(next);
	}

private:
	WasmStackFrame<T>& top_frame()
	{
		if(frames_.empty())
			throw std::out_of_range("Attempt to access out-of-bounds call stack frame.");
		return frames_.top();
	}

	SimpleStack<T>& top_stack()
	{
		if(frames_.empty())
			return base_stack_;
		return frames_.top().stack();
	}

	std::pair<WasmStackFrame<T>&, WasmStackFrame<T>&> top_two_frames()
	{
		if(frames_.empty())
			throw std::out_of_range("Attempt to access out-of-bounds frame.");
		auto pos = frames_.begin();
		const auto& hi = *pos++;
		if(pos == frames_.end())
			throw std::out_of_range("Attempt to access out-of-bounds frame.");
		const auto& lo = *pos++;
		return std::pair<WasmStackFrame&, WasmStackFrame&>(lo, hi);
	}

	void _recurse_return_from_frame_unchecked()
	{
		if(auto& stack = top_frame().stack(); stack.empty())
		{
			frames_.pop();
		}
		else
		{
			T ret_v = stack.pop();
			_recurse_return_from_frame_unchecked();
			top_stack().emplace(ret_v);
		}
	}
	void _recurse_return_from_frame(const gsl::span<const LanguageType>& return_types, std::size_t index = 0u)
	{
		assert(not frames_.empty());
		if(index < return_types.size())
		{
			// one of the return values from the callee's frame
			tp::visit_value_type(
				[&](auto WasmValue::* p) {
					auto ret_v = top_frame().stack_at(index, p);
					// keep recursing until we exhaust the return
					// values.  Once we've done that, pop the callee's
					// frame, and then pop the locals vector and arguments
					// off the caller's stack.
					_recurse_return_from_frame(return_types, index);
					// now that we've popped the callee's frame from the call stack,
					// and the locals vector off of the caller's frame's stack, push
					// the return values onto the callee's frame's stack.
					//
					// note that there may be no frame below the one we popped, in which
					// case, the return values are pushed onto 'this->base_stack_'
					top_stack().stack_emplace_top(p, ret_v);
				},
				return_types[return_types.size() - index]
			);
		}
		else
		{
			const auto* func_p = top_frame.function();
			assert(func_p);
			const auto& func = *func_p;
			std::size_t total_locals = param_count(func) + locals_count(func);
			// pop the top frame
			frames_.pop_front();
			assert(not frames_.empty());
			auto& new_top_frame = frames_.front();
			assert(new_top_frame.stack_size() >= total_locals);
			// pop the locals (args + locals) vector of the frame we just popped
			// off of the stack.  note that the locals vector of the current frame
			// lives on the stack of the previous frame.
			new_top_frame.stack_pop_n(total_locals);
			// push the return values of the frame we just popped...
			// ... by simply returning from _recurse_return_from_frame().  each
			// recursive call which took the if() branch above will push the return
			// values onto the new frame.
		}
	}

	const StackResource& stack_resource() const
	{ return base_stack_.get_resource(); }

	StackResource& stack_resource()
	{ return base_stack_.get_resource(); }

	frame_stack_type frames_;
	SimpleStack<T> base_stack_;
};

template <class T>
const T& local_at(const WasmCallStack<T>& self, std::size_t idx)
{
	auto locals = locals(current_frame(self));
	if(idx >= locals.size())
		throw ValidationError<std::out_of_range>("Attempt to access out-of-bounds local.");
	return locals[idx];
}

template <class T>
T& local_at(WasmCallStack<T>& self, std::size_t idx)
{
	auto locals = locals(current_frame(self));
	if(idx >= locals.size())
		throw ValidationError<std::out_of_range>("Attempt to access out-of-bounds local.");
	return locals[idx];
}

template <class T>
const auto& current_stack(const WasmCallStack<T>& self, std::size_t idx)
{ return current_stack(current_frame(self)); }

template <class T>
auto& current_stack(WasmCallStack<T>& self, std::size_t idx)
{ return current_stack(current_frame(self)); }

} /* namespace wasm */

#endif /* VM_CALL_STACK_H */

max_stars_repo_name: tvanslyke/wasm-cpp

max_issues_repo_name: tvanslyke/wasm-cpp

max_forks_repo_name: tvanslyke/wasm-cpp


text:
// The MIT License (MIT)
//
// Copyright (c) 2018 Mateusz Pusz
//
// Permission is hereby granted, free of charge, to any person obtaining a copy
// of this software and associated documentation files (the "Software"), to deal
// in the Software without restriction, including without limitation the rights
// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
// copies of the Software, and to permit persons to whom the Software is
// furnished to do so, subject to the following conditions:
//
// The above copyright notice and this permission notice shall be included in all
// copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
// SOFTWARE.

#pragma once

#include <gsl/gsl_assert>
#include <units/bits/math_concepts.h>
#include <units/bits/pow.h>
#include <units/bits/ratio_maths.h>

namespace units::detail {

struct decimal_fp {
  double significant;
  std::intmax_t exponent;
};

[[nodiscard]] constexpr decimal_fp to_decimal(double v) noexcept
{
  if (v == 0) {
    return {.significant = 0.0, .exponent = 0};
  }

  double significant = abs(v);
  std::intmax_t exponent = 0;

  while (significant < 1) {
    significant *= 10.0;
    --exponent;
  }

  while (significant >= 10) {
    significant /= 10.0;
    ++exponent;
  }

  if (v < 0) {
    significant = -significant;
  }

  return {.significant = significant, .exponent = exponent};
}

/* approximate natural log as https://math.stackexchange.com/a/977836
 far slower than std::log but works at compile time with similar accuracy
 */
[[nodiscard]] constexpr double constexpr_log(double v) noexcept
{
  Expects(v > 0);

  // lookup table to speed up convergence for all significant values
  // significant values of 7 and greater benefit mostly as they now converge in 5 terms compared to O(10)-O(100)
  // required without the table
  //
  // using python:
  // >>> import math
  // >>> for i in range(1, 100):
  // ...     print(f"/* log({i:>2d}) = */ {math.log(i):.16f},")
  constexpr std::array<double, 99> log_table{
      /* log( 1) = */ 0.0000000000000000,
      /* log( 2) = */ 0.6931471805599453,
      /* log( 3) = */ 1.0986122886681098,
      /* log( 4) = */ 1.3862943611198906,
      /* log( 5) = */ 1.6094379124341003,
      /* log( 6) = */ 1.7917594692280550,
      /* log( 7) = */ 1.9459101490553132,
      /* log( 8) = */ 2.0794415416798357,
      /* log( 9) = */ 2.1972245773362196,
      /* log(10) = */ 2.3025850929940459,
      /* log(11) = */ 2.3978952727983707,
      /* log(12) = */ 2.4849066497880004,
      /* log(13) = */ 2.5649493574615367,
      /* log(14) = */ 2.6390573296152584,
      /* log(15) = */ 2.7080502011022101,
      /* log(16) = */ 2.7725887222397811,
      /* log(17) = */ 2.8332133440562162,
      /* log(18) = */ 2.8903717578961645,
      /* log(19) = */ 2.9444389791664403,
      /* log(20) = */ 2.9957322735539909,
      /* log(21) = */ 3.0445224377234230,
      /* log(22) = */ 3.0910424533583161,
      /* log(23) = */ 3.1354942159291497,
      /* log(24) = */ 3.1780538303479458,
      /* log(25) = */ 3.2188758248682006,
      /* log(26) = */ 3.2580965380214821,
      /* log(27) = */ 3.2958368660043291,
      /* log(28) = */ 3.3322045101752038,
      /* log(29) = */ 3.3672958299864741,
      /* log(30) = */ 3.4011973816621555,
      /* log(31) = */ 3.4339872044851463,
      /* log(32) = */ 3.4657359027997265,
      /* log(33) = */ 3.4965075614664802,
      /* log(34) = */ 3.5263605246161616,
      /* log(35) = */ 3.5553480614894135,
      /* log(36) = */ 3.5835189384561099,
      /* log(37) = */ 3.6109179126442243,
      /* log(38) = */ 3.6375861597263857,
      /* log(39) = */ 3.6635616461296463,
      /* log(40) = */ 3.6888794541139363,
      /* log(41) = */ 3.7135720667043080,
      /* log(42) = */ 3.7376696182833684,
      /* log(43) = */ 3.7612001156935624,
      /* log(44) = */ 3.7841896339182610,
      /* log(45) = */ 3.8066624897703196,
      /* log(46) = */ 3.8286413964890951,
      /* log(47) = */ 3.8501476017100584,
      /* log(48) = */ 3.8712010109078911,
      /* log(49) = */ 3.8918202981106265,
      /* log(50) = */ 3.9120230054281460,
      /* log(51) = */ 3.9318256327243257,
      /* log(52) = */ 3.9512437185814275,
      /* log(53) = */ 3.9702919135521220,
      /* log(54) = */ 3.9889840465642745,
      /* log(55) = */ 4.0073331852324712,
      /* log(56) = */ 4.0253516907351496,
      /* log(57) = */ 4.0430512678345503,
      /* log(58) = */ 4.0604430105464191,
      /* log(59) = */ 4.0775374439057197,
      /* log(60) = */ 4.0943445622221004,
      /* log(61) = */ 4.1108738641733114,
      /* log(62) = */ 4.1271343850450917,
      /* log(63) = */ 4.1431347263915326,
      /* log(64) = */ 4.1588830833596715,
      /* log(65) = */ 4.1743872698956368,
      /* log(66) = */ 4.1896547420264252,
      /* log(67) = */ 4.2046926193909657,
      /* log(68) = */ 4.2195077051761070,
      /* log(69) = */ 4.2341065045972597,
      /* log(70) = */ 4.2484952420493594,
      /* log(71) = */ 4.2626798770413155,
      /* log(72) = */ 4.2766661190160553,
      /* log(73) = */ 4.2904594411483910,
      /* log(74) = */ 4.3040650932041702,
      /* log(75) = */ 4.3174881135363101,
      /* log(76) = */ 4.3307333402863311,
      /* log(77) = */ 4.3438054218536841,
      /* log(78) = */ 4.3567088266895917,
      /* log(79) = */ 4.3694478524670215,
      /* log(80) = */ 4.3820266346738812,
      /* log(81) = */ 4.3944491546724391,
      /* log(82) = */ 4.4067192472642533,
      /* log(83) = */ 4.4188406077965983,
      /* log(84) = */ 4.4308167988433134,
      /* log(85) = */ 4.4426512564903167,
      /* log(86) = */ 4.4543472962535073,
      /* log(87) = */ 4.4659081186545837,
      /* log(88) = */ 4.4773368144782069,
      /* log(89) = */ 4.4886363697321396,
      /* log(90) = */ 4.4998096703302650,
      /* log(91) = */ 4.5108595065168497,
      /* log(92) = */ 4.5217885770490405,
      /* log(93) = */ 4.5325994931532563,
      /* log(94) = */ 4.5432947822700038,
      /* log(95) = */ 4.5538768916005408,
      /* log(96) = */ 4.5643481914678361,
      /* log(97) = */ 4.5747109785033828,
      /* log(98) = */ 4.5849674786705723,
      /* log(99) = */ 4.5951198501345898,
  };
  decimal_fp x = to_decimal(v);

  // dividing the significant by nearest lower value in [1.0, 1.1, 1.2, ..., 9.9] will greatly improve convergence
  x.significant *= 10;
  const auto isignificant = static_cast<std::size_t>(x.significant);
  x.significant /= static_cast<double>(isignificant);
  const double result = static_cast<double>(x.exponent - 1) * log_table[9] + log_table[isignificant - 1];

  // 1.0 <= significant < 1.1 converges rapidly
  const double y = (x.significant - 1) / (x.significant + 1);
  const double y_squared = y * y;
  double sum = 0;
  // 5 terms are needed for convergence to machine precision in the worst case scenario
  for (int k = 4; k > 0; --k) {
    sum = y_squared * (1 / (2 * static_cast<double>(k) + 1) + sum);
  }
  sum = 2 * y * (1 + sum);  // k = 0 term
  return result + sum;
}

/* approximate e^x as Taylor series e^x = 1 + x/1! + x^2/2! + x^3/3! +... where N is the order of the Taylor series
 use https://math.stackexchange.com/a/1988927 to improve convergence for large values

 larger Factor values improve convergence for all values but reduce the precision
*/
template<std::size_t N, std::intmax_t Factor = 256>
  requires gt_zero<Factor>
[[nodiscard]] constexpr double constexpr_exp(double v) noexcept
{
  if constexpr (N == 0) {
    return 1.0;
  } else {
    constexpr auto coefficients = []() {
      std::array<double, N> coeffs;
      std::size_t factorial = 1;
      for (std::size_t i = 0; i < N; ++i) {
        factorial *= i + 1;
        coeffs[i] = 1.0 / static_cast<double>(factorial);
      }
      return coeffs;
    }();

    const double x = v / static_cast<double>(Factor);
    double result = 0;
    for (auto i = static_cast<std::intmax_t>(N - 1); i >= 0; --i) {
      result = x * (coefficients[static_cast<std::size_t>(i)] + result);
    }

    // for factors of power of 2 this should be replaced by log2(Factor) multiplications by the compiler
    return pow_impl<Factor>(1 + result);
  }
}

// default template arguments provide reasonable precision even for fairly large exponents
// see constexpr_exp for template arguments
template<std::size_t ExpOrder = 10, std::intmax_t Factor = 128>
[[nodiscard]] constexpr double constexpr_pow(double v, double exponent) noexcept
{
  const double x = exponent * constexpr_log(v);
  return constexpr_exp<ExpOrder, Factor>(x);
}

}  // namespace units::detail

ext: h


text:
#ifndef DEINOS_CHESS_H
#define DEINOS_CHESS_H
#include <optional>
#include <array>
#include <vector>
#include <string>
#include <cassert>
#include <iostream>
#include <gsl/gsl_assert>
#include <cstddef>

namespace chess {
	enum class Almnt : uint8_t;
	enum class AlRes : uint8_t;
	enum class GameResult : uint8_t;
	enum class Side : uint8_t;
	struct Piece;
	class Square;
	struct MoveRecord;
	class Move;
	class Position;

	//Alignment (i.e. White or Black)
	enum class Almnt : uint8_t {White = 0, Black = 1};
	constexpr Almnt operator!(const Almnt& a) {return (a == Almnt::White ? Almnt::Black : Almnt::White);}
	inline uint8_t as_index(Almnt a) {return static_cast<uint8_t>(a);}

	//Alignment Result (for cases where "None" is an option)
	enum class AlRes : uint8_t {White = 0, Black = 1, None = 2}; //Utility to avoid constant Empty checks
	inline bool operator==(AlRes ar, Almnt a) {return static_cast<uint8_t>(ar) == static_cast<uint8_t>(a);}
	inline bool operator==(Almnt a, AlRes ar) {return ar == a;}
	inline bool operator!=(AlRes ar, Almnt a) {return !(ar == a);}
	inline bool operator!=(Almnt a, AlRes ar) {return !(ar == a);}

	//Game Result
	enum class GameResult : uint8_t {White = 0, Black = 1, Draw = 2};
	inline float evaluate(GameResult gr) {
		switch (gr) {
			case GameResult::White: return 1.0f;
			case GameResult::Black: return 0.0f;
			case GameResult::Draw: return 0.5f;
			default: return 0.5f;
		}
	}
	inline GameResult victory(Almnt a) {return static_cast<GameResult>(a);}
	inline AlRes victor(GameResult gr) {return static_cast<AlRes>(gr);}

	//Side of Board
	enum class Side : uint8_t {Kingside = 0, Queenside = 1};

	//The contents of a square (i.e. a chess piece or empty)
	struct Piece {
		enum class Type : uint8_t {Empty, Pawn, Knight, Bishop, Rook, Queen, King};

		constexpr Piece() = default;
		constexpr Piece(Almnt a, Type t) : data((static_cast<std::byte>(t) << 1) | static_cast<std::byte>(a)) {}

		inline AlRes almnt_res() const {return (type() == Type::Empty ? AlRes::None : static_cast<AlRes>(almnt()));}
		inline Almnt almnt() const {return static_cast<Almnt>(data & std::byte{0x01});}
		inline Type type() const {return static_cast<Type>(data >> 1);}

	private:
		constexpr Piece(uint8_t d) : data(std::byte{d}) {}
		inline uint8_t raw() const {return static_cast<uint8_t>(data);}
		std::byte data = std::byte{0};
		friend class Position;
	};

	bool operator==(Piece, Piece);
	bool operator!=(Piece, Piece);
	std::ostream& operator<<(std::ostream& os, Piece::Type t);
	std::ostream& operator<<(std::ostream& os, Piece p); //White uppercase, black lowercase

	//A bounds checked index into a chess position
	class Square{
	public:
		constexpr Square() = default;
		constexpr Square(const char* name){ //this must be a zstring TODO
			constexpr std::array<char, 8> files {'a', 'b', 'c', 'd', 'e', 'f', 'g', 'h'};
			constexpr std::array<char, 8> ranks {'1', '2', '3', '4', '5', '6', '7', '8'};
			uint8_t file = 255;
			uint8_t rank = 255;
			for (int i = 0; i < 8; i++) {
				if (files[i] == name[0]) file = i;
				if (ranks[i] == name[1]) rank = i;
			}
			assert(file != 255);
			assert(rank != 255);
			data = (std::byte{rank} << 3) | std::byte{file};
		}

		std::optional<Square> translate(int files, int ranks) const; //TODO consider changing type to int8_t
		inline uint8_t file() const {return std::to_integer<uint8_t>(data & std::byte{0b00000111});}
		inline uint8_t rank() const {return std::to_integer<uint8_t>((data & std::byte{0b00111000}) >> 3);}
		explicit operator std::string () const;
	private:
		constexpr Square(uint8_t file, uint8_t rank) : data((std::byte{rank} << 3) | std::byte{file}) {
			assert(file < 8);
			assert(rank < 8);
		}
		constexpr Square(std::byte d) : data(d) {} //construct from raw data
		std::byte data = std::byte{0};
		friend class MoveRecord;
		friend bool operator==(Square, Square);
	};
	bool operator==(Square, Square);
	bool operator!=(Square, Square);

	constexpr std::array<Square, 64> all_squares {
		"a1", "b1", "c1", "d1", "e1", "f1", "g1", "h1",
		"a2", "b2", "c2", "d2", "e2", "f2", "g2", "h2",
		"a3", "b3", "c3", "d3", "e3", "f3", "g3", "h3",
		"a4", "b4", "c4", "d4", "e4", "f4", "g4", "h4",
		"a5", "b5", "c5", "d5", "e5", "f5", "g5", "h5",
		"a6", "b6", "c6", "d6", "e6", "f6", "g6", "h6",
		"a7", "b7", "c7", "d7", "e7", "f7", "g7", "h7",
		"a8", "b8", "c8", "d8", "e8", "f8", "g8", "h8",
		};

	//stores information necessary to reconstruct a move from a position
	struct MoveRecord {
	private:

		inline uint8_t promo_index() const {return static_cast<uint8_t>(data2) >> 6;}
	public:
		constexpr static std::array<Piece::Type, 4> promo_types
			{Piece::Type::Knight, Piece::Type::Bishop, Piece::Type::Rook, Piece::Type::Queen};

		static_assert(sizeof(Square) == 1);
		constexpr MoveRecord(Square t_initial, Square t_final)
			: data1(reinterpret_cast<std::byte&>(t_initial)), data2(reinterpret_cast<std::byte&>(t_final)) {}
		MoveRecord(Square t_initial, Square t_final, Piece::Type t_promo_type);

		inline Square initial() const {return Square(data1 & std::byte{0x3F});}
		inline Square final() const {return Square(data2 & std::byte{0x3F});}
		inline bool is_promo() const {return (data1 & std::byte{0b11000000}) != std::byte{0};}
		inline std::optional<Piece::Type> promo_type() const
		{
			return is_promo() ? std::make_optional<Piece::Type>(promo_types[promo_index()]) : std::nullopt;
		};
		std::string to_string() const;
	private:
		std::byte data1 = std::byte{0}; //initial square and promotion existence flag
		std::byte data2 = std::byte{0}; //final square and promotion type flag
	};
	bool operator==(const MoveRecord&, const MoveRecord&);
	inline bool operator!=(const MoveRecord& mr1, const MoveRecord& mr2) {return !(mr1 == mr2);}
	inline std::ostream& operator<<(std::ostream& os, const MoveRecord& mr) {os << mr.to_string(); return os;}

	//64 half-byte uints stored in 32B
	struct HalfByteBoard {
	private:
		inline uint8_t get_front(int i) const {return data[i] >> 4;}
		inline uint8_t get_back(int i) const {return data[i] & 0x0F;}
		inline void set_front(int i, uint8_t val) {data[i] = (val << 4) | (data[i] & 0x0F);}
		inline void set_back(int i, uint8_t val) {data[i] = val | (data[i] & 0xF0);}
	public:
		inline uint8_t get(int i) const
		{
			assert(i >= 0);
			assert(i < 64);
			if (i % 2 == 0) return get_front(i / 2);
			else return get_back(i / 2);
		}
		inline uint8_t get(int file, int rank) const {return get(rank * 8 + file);}
		inline void set(int i, uint8_t val)
		{
			assert(i >= 0);
			assert(i < 64);
			assert(val <= 0x0F);
			if (i % 2 == 0) set_front(i / 2, val);
			else set_back(i / 2, val);
		}
		inline void set(int file, int rank, uint8_t val) {set(rank * 8 + file, val);}
	private:
		std::array<uint8_t, 32> data = {0};
		friend bool operator==(const HalfByteBoard&, const HalfByteBoard&);
	};
	inline bool operator==(const HalfByteBoard& hbb1, const HalfByteBoard& hbb2) {return hbb1.data == hbb2.data;}
	inline bool operator!=(const HalfByteBoard& hbb1, const HalfByteBoard& hbb2) {return !(hbb1 == hbb2);}

	//Equivalent to a chess position recorded in Forsyth-Edwards Notaion
	class Position {
	public:
		constexpr Position() = default;
		Position(const std::string& fen); //construct from FEN representation
		Position(const Position&, const Move&); //generate new position by applying a move
		static Position std_start();

		//inline Piece& operator[] (Square s) {return m_board[s.file()][s.rank()];}
		//inline Piece operator[] (Square s) const {return m_board[s.file()][s.rank()];}
		inline Piece at(Square s) const {return Piece(m_board.get(s.file(), s.rank()));}
		//inline Piece at(Square s) const {return m_board[s.file()][s.rank()];}
		inline Piece at(int index) const {return Piece(m_board.get(index));}
		//inline Piece at(int index) const {return m_board[index % 8][index / 8];}
		inline void set(Square s, Piece p) {m_board.set(s.file(), s.rank(), p.raw());}
		//inline void set(Square s, Piece p) {m_board[s.file()][s.rank()] = p;}
		inline void set(int index, Piece p) {m_board.set(index, p.raw());}
		//inline void set(int index, Piece p) {m_board[index % 8][index / 8] = p;}
		inline bool can_castle(Almnt a, Side s) const {return m_castle[(int) a][(int) s];} //TODO use as_index
		inline bool& mut_castle(Almnt a, Side s){return m_castle[(int) a][(int) s];}
		inline Almnt to_move() const {return m_to_move;}
		inline std::optional<Square> en_passant_target() const {return m_en_passant_target;}
		inline int hm_clock() const {return m_hm_clock;}
		inline int fm_count() const {return m_fm_count;}

		std::string as_fen() const;
		explicit operator std::string() const;
		friend std::ostream& operator<<(std::ostream& os, const Position& pos);
		friend bool operator==(const Position&, const Position&);

	private:
		static_assert(sizeof(Piece) == 1);
		HalfByteBoard m_board;
		//std::array<std::array<Piece, 8>, 8> m_board = {Piece()};
		Almnt m_to_move = Almnt::White;
		std::array<std::array<bool, 2>, 2> m_castle = {{false}}; //"permitted to castle" flags indexed by enum values
		//std::optional<GameResult> m_result = std::nullopt;
		std::optional<Square> m_en_passant_target = std::nullopt;
		short m_hm_clock = 0;
		short m_fm_count = 1;
	};
	bool operator==(const Position&, const Position&);
	bool operator!=(const Position&, const Position&);
	std::ostream& operator<<(std::ostream& os, const Position& pos);

	//A wrapper for a position and move record, representing a single move
	class Move { //TODO: test code
	public:
		//constexpr Move() = default; //default constructor
		//Move(Square initial, Square final, Piece moved, Piece captured); //normal move
		//Move(const Position& pos, const std::string& name); //N.B this constructor does not support castling, promotion or en passant
		//void set_castling(); //castling 1.must be a king moving 2.must not capture 3.can only set once
		//void set_promotion(Piece::Type); //pawn promotion
		//void set_en_passant(); //en_passant
		constexpr Move(const Position& t_pos, const MoveRecord& t_record) : m_pos(t_pos), m_record(t_record) {
			if (moved().type() == Piece::Type::Empty) {std::cerr << *this << std::endl << t_pos << std::endl; throw std::exception();}
			if (captured().almnt_res() == moved().almnt()) {std::cerr << *this << " capt: " << captured() << " mr: " << m_record << std::endl << t_pos << std::endl; throw std::exception();}
			assert(moved().type() != Piece::Type::Empty);
			assert(captured().almnt_res() != moved().almnt());
		}
		//inline Square initial_square() const {return m_initial_square;}
		inline Square initial_sq() const {return m_record.initial();}
		//inline Square final_square() const {return m_final_square;}
		inline Square final_sq() const {return m_record.final();}
		//inline Piece moved_piece() const {return m_moved_piece;}
		inline Piece moved() const {return m_pos.at(initial_sq());}
		//inline Piece captured_piece() const {return m_captured_piece;}
		inline Piece captured() const {return m_pos.at(final_sq());} //returns empty piece if en passant
		//inline bool is_en_passant() const {return m_is_en_passant;}
		bool is_en_passant() const;
		//inline bool is_castling() const {return m_is_castling;}
		bool is_castling() const;
		//inline bool is_promotion() const {return m_promoted_to.type() != Piece::Type::Empty;}
		inline bool is_promotion() const {return m_record.is_promo();}
		//inline Piece promoted_to() const {return m_promoted_to;}
		inline std::optional<Piece::Type> promo_type() const {return m_record.promo_type();}
		inline const MoveRecord& record() const {return m_record;}
		inline Position apply() const {return Position(m_pos, *this);}

		explicit operator std::string() const;
		friend std::ostream& operator<<(std::ostream& os, const Move& mv);
		std::string to_xboard() const;

	private:
		//Square m_initial_square;
		//Square m_final_square;
		//Piece m_moved_piece;
		//Piece m_captured_piece;
		//Piece m_promoted_to;
		//bool m_is_castling = false;
		//bool m_is_en_passant = false;
		const Position& m_pos;
		const MoveRecord& m_record;
		friend bool operator==(const Move&, const Move&);
	};
	bool operator==(const Move&, const Move&);
	bool operator!=(const Move&, const Move&);
	std::ostream& operator<<(std::ostream& os, const Move& mv);
	//std::optional<Move> find_move(const std::string& t_name, const AnalysedPosition& apos);

	//std::optional<Move> move_from_fen(const std::string& fen, const Position& curr_pos, const std::vector<Move>& moves);
}
#endif


text:
#include <petsc.h>

static PetscErrorCode SetInitialCondition(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx){

    u[0] = x[0];
    return 0;
}

static PetscErrorCode PrintFVGradients(DM dm){
    PetscFunctionBeginUser;
    PetscErrorCode ierr;

    // Convert the dm to a plex
    DM plex;
    DMConvert(dm, DMPLEX, &plex);

    // Get the start/end
    PetscInt fStart;
    PetscInt fEnd;
    ierr = DMPlexGetHeightStratum(plex, 1, &fStart, &fEnd);CHKERRQ(ierr);

    // Extract the cell geometry, and the dm that holds the information
    Vec faceGeometry;
    DM dmFace;
    ierr = DMPlexGetGeometryFVM(plex, &faceGeometry, NULL, NULL);CHKERRQ(ierr);

    ierr = VecGetDM(faceGeometry, &dmFace);CHKERRQ(ierr);

    const PetscScalar *facegeom;
    ierr = VecGetArrayRead(faceGeometry, &facegeom);CHKERRQ(ierr);

    // Get the dim
    PetscInt dim;
    ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr);

    // March over each cell volume
    for (PetscInt f = fStart; f < fEnd; ++f) {
        PetscFVFaceGeom       *fg;

        ierr = DMPlexPointLocalRead(dmFace, f, facegeom, &fg);CHKERRQ(ierr);

        printf("Face Grad %d: ", f);
        for(PetscInt d =0; d< dim; d++){
            printf("%f/%f, ", fg->grad[0][d], fg->grad[1][d]);
        }
        printf("\n");
    }

    PetscFunctionReturn(0);
}

int main(int argc, char **argv)
{
    DM             dm;
    PetscErrorCode ierr;

    ierr = PetscInitialize(&argc, &argv, NULL, "HELP");if (ierr) return(ierr);
    ierr = DMPlexCreateBoxMesh(PETSC_COMM_WORLD, 2, PETSC_FALSE, NULL, NULL, NULL, NULL, PETSC_TRUE, &dm);CHKERRQ(ierr);
    ierr = DMViewFromOptions(dm, NULL, "-dm_view");CHKERRQ(ierr);
    ierr = DMSetFromOptions(dm);CHKERRQ(ierr);

    {
        DM gdm;

        ierr = DMPlexConstructGhostCells(dm, NULL, NULL, &gdm);CHKERRQ(ierr);
        ierr = DMDestroy(&dm);CHKERRQ(ierr);
        dm   = gdm;
    }
    ierr = DMViewFromOptions(dm, NULL, "-dm_view_ghost");CHKERRQ(ierr);

    PetscFV fvm;
    PetscDS ds;

    ierr = PetscFVCreate(PETSC_COMM_WORLD, &fvm);CHKERRQ(ierr);
    ierr = PetscFVSetFromOptions(fvm);CHKERRQ(ierr);
    ierr = PetscFVSetNumComponents(fvm, 1);CHKERRQ(ierr);
    ierr = PetscFVSetSpatialDimension(fvm, 2);CHKERRQ(ierr);
    ierr = PetscObjectSetName((PetscObject) fvm, "outputField");CHKERRQ(ierr);

    /* FV is now structured with one field having all physics as components */
    ierr = DMAddField(dm, NULL, (PetscObject) fvm);CHKERRQ(ierr);
    ierr = PetscFVDestroy(&fvm);CHKERRQ(ierr);
    ierr = DMCreateDS(dm);CHKERRQ(ierr);
    ierr = DMGetDS(dm, &ds);CHKERRQ(ierr);

    Vec X;
    PetscErrorCode     (*func[1]) (PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx) = {SetInitialCondition};

    ierr = DMCreateGlobalVector(dm, &X);CHKERRQ(ierr);
    ierr = PetscObjectSetName((PetscObject) X, "solution");CHKERRQ(ierr);
    ierr = DMProjectFunction(dm, 0.0, func, NULL, INSERT_ALL_VALUES, X);CHKERRQ(ierr);
    ierr = VecViewFromOptions(X, NULL, "-vec_view");CHKERRQ(ierr);

    PrintFVGradients(dm);

    ierr = VecDestroy(&X);CHKERRQ(ierr);
    ierr = DMDestroy(&dm);CHKERRQ(ierr);
    return PetscFinalize();
}

ext: c

max_stars_repo_name: mmcgurn/MattFlowCases

max_issues_repo_name: mmcgurn/MattFlowCases

max_forks_repo_name: mmcgurn/MattFlowCases


text:
#include <gsl/gsl_math.h>
#include "gsl_cblas.h"
#include "cblas.h"

void
cblas_cswap (const int N, void *X, const int incX, void *Y, const int incY)
{
#define BASE float
#include "source_swap_c.h"
#undef BASE
}


text:
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include "compearth.h"
#ifdef COMPEARTH_USE_MKL
#ifdef __clang__
#pragma clang diagnostic push
#pragma clang diagnostic ignored "-Wreserved-id-macro"
#pragma clang diagnostic ignored "-Wstrict-prototypes"
#endif
#include <mkl_cblas.h>
#ifdef __clang__
#pragma clang diagnostic pop
#endif
#else
#include <cblas.h>
#endif

/*!
 * @brief Computes the omega angles between two moment tensors.
 *
 * @param[in] nmt1    Number of moment tensors in M1.  This can be 1 otherwise
 *                    it must be equal to nmt2.  If it is 1 then the angle will
 *                    be between M1 and all the moment tensors in M2.
 * @param[in] M1      First array of moment tensors.  M1 is packed
 *                    \f$ M = \{M_{11} M_{22} M_{33} M_{12} M_{13} M_{23}\} \f$.
 *                    This is an array of dimension [6 x nmt1] with leading
 *                    dimension 6.
 * @param[in] nmt2    Number of moment tensors in M2.  This can be 1 otherwise
 *                    it must be equal to nmt1.  If it is 1 then the angle will
 *                    be between M2 and all the moment tensors in M1.
 * @param[in] M2      Second array of moment tensors.  M2 is packed
 *                    \f$ M = \{M_{11} M_{22} M_{33} M_{12} M_{13} M_{23}\} \f$.
 *                    This is an array of dimension [6 x nmt2] with leading
 *                    dimension 6.
 *
 * @param[out] omega  This is a measure between two moment tensors and combines
 *                    differences in eigenvalues (source types) and orientation
 *                    but not magnitude.
 *                    This is an array of dimension [MAX(nmt1, nmt2)].
 *
 * @author Carl Tape and translated to C by Ben Baker
 *
 * @copyright MIT
 *
 */
int compearth_CMT2omega(const int nmt1, const double *__restrict__ M1,
                        const int nmt2, const double *__restrict__ M2,
                        double *__restrict__ omega)
{
    double Mwork[6*CE_CHUNKSIZE] __attribute__((aligned(64)));
    int ierr, imt, nmtLoc;
    const double pi180i = 180.0/M_PI;
    ierr = 0;
    if (nmt1 < 1 || M1 == NULL || nmt2 < 1 || M2 == NULL || omega == NULL)
    {
        if (nmt1 < 1)
        {
            fprintf(stderr, "%s: Error nmt1 must be positive\n", __func__);
        }
        if (nmt2 < 1)
        {
            fprintf(stderr, "%s: Error nmt2 must be positive\n", __func__);
        }
        if (M1 == NULL){fprintf(stderr, "%s: Error M1 is NULL\n", __func__);}
        if (M2 == NULL){fprintf(stderr, "%s: Error M2 is NULL\n", __func__);}
        if (omega == NULL)
        {
            fprintf(stderr, "%s: Error omega is NULL\n", __func__);
        }
        return -1;
    }
    // Straight compuation
    if (nmt1 == nmt2)
    {
        ierr = compearth_angleMT(nmt1, M1, M2, omega);
        if (ierr != 0)
        {
            fprintf(stderr, "%s: Error computing angleMT 1\n", __func__);
        }
        // Convert to degrees
        cblas_dscal(nmt1, pi180i, omega, 1);
        return ierr;
    }
    else
    {
        if (nmt1 == 1 || nmt2 == 1)
        {
            // Copy nmt1 and repeat computation with it
            if (nmt1 == 1)
            {
                for (imt=0; imt<CE_CHUNKSIZE; imt++)
                {
                    cblas_dcopy(6, &M1[0], 1, &Mwork[6*imt], 1);
                }
                for (imt=0; imt<nmt2; imt=imt+CE_CHUNKSIZE)
                {
                    nmtLoc = MIN(CE_CHUNKSIZE, nmt2 - imt);
                    ierr = compearth_angleMT(nmtLoc, Mwork, &M2[6*imt],
                                             &omega[imt]);
                    if (ierr != 0)
                    {
                        fprintf(stderr, "%s: Error calling angleMT at %d\n",
                                __func__, __LINE__);
                        return ierr;
                    }
                    // Convert to degrees
                    cblas_dscal(nmtLoc, pi180i, &omega[imt], 1);
                }
            }
            // Copy nmt2 and repeat computation with it
            else
            {
                for (imt=0; imt<CE_CHUNKSIZE; imt++)
                {
                    cblas_dcopy(6, &M2[0], 1, &Mwork[6*imt], 1);
                }
                for (imt=0; imt<nmt1; imt=imt+CE_CHUNKSIZE)
                {
                    nmtLoc = MIN(CE_CHUNKSIZE, nmt1 - imt);
                    ierr = compearth_angleMT(nmtLoc, &M1[6*imt], Mwork,
                                             &omega[imt]);
                    if (ierr != 0)
                    {
                        fprintf(stderr, "%s: Error calling angleMT at %d\n",
                                __func__, __LINE__);
                        return ierr;
                    }
                    // Convert to degrees
                    cblas_dscal(nmtLoc, pi180i, &omega[imt], 1);
                }
            }
        }
        else
        {
            fprintf(stderr, "%s: nmt1=%d /= nmt2=%d; this is invalid\n",
                    __func__, nmt1, nmt2);
            return -1;
        }
    }
    return ierr;
}

max_stars_repo_name: OUCyf/mtbeach

max_issues_repo_name: carltape/mtbeach

max_forks_repo_name: carltape/mtbeach


text:
#ifndef BATCH_EXTRACTOR_H_IVH3CLLQ
#define BATCH_EXTRACTOR_H_IVH3CLLQ

#include <algorithm>
#include <gsl/gsl>
#include <opencv2/core/types.hpp>
#include <opencv2/features2d.hpp>
#include <opencv2/xfeatures2d.hpp>
#include <optional>
#include <sens_loc/util/correctness_util.h>

using namespace std;
using namespace cv;

namespace sens_loc::apps {

/// Helper class that visits a list of images and extracts features with the
/// provided detectors.
/// \ingroup feature-extractor-driver
class batch_extractor {
  public:
    using filter_func = function<vector<KeyPoint>::iterator(vector<KeyPoint>&)>;

    /// \param detector,descriptor Algorithms used for detection and
    /// description.
    /// \param input_pattern,output_pattern Formattable string for IO.
    /// \param keypoint_filter Callable that determines if a keypoint shall be
    /// dropped from consideration. All keypoints with 'keypoint_filter(kp) ==
    /// true' are removed.
    batch_extractor(Ptr<Feature2D>      detector,
                    Ptr<Feature2D>      descriptor,
                    string_view         input_pattern,
                    string_view         output_pattern,
                    vector<filter_func> keypoint_filter)
        : _detector{move(detector)}
        , _descriptor{move(descriptor)}
        , _input_pattern{input_pattern}
        , _ouput_pattern{output_pattern}
        , _keypoint_filter{move(keypoint_filter)} {
        Expects(!_input_pattern.empty());
        Expects(!_ouput_pattern.empty());
        Expects(!_detector.empty());
    }

    /// Process a whole batch of files in the range [start, end].
    [[nodiscard]] bool process_batch(int start, int end) const noexcept;

  private:
    /// Detect and describe one single index. Handles the IO as well.
    [[nodiscard]] bool process_index(int idx) const noexcept;

    /// Do IO and handle detection down to \c compute_features.
    bool process_detector(const math::image<uchar>& image,
                          const string&             out_file,
                          const string&             in_file) const noexcept;

    /// Compute and filter keypoints and run the descriptor on them
    /// afterwards.
    [[nodiscard]] pair<vector<KeyPoint>, Mat>
    compute_features(const math::image<uchar>& img) const noexcept;

    mutable Ptr<Feature2D> _detector;
    mutable Ptr<Feature2D> _descriptor;
    string_view            _input_pattern;
    string_view            _ouput_pattern;
    vector<filter_func>    _keypoint_filter;
};
}  // namespace sens_loc::apps

#endif /* end of include guard: BATCH_EXTRACTOR_H_IVH3CLLQ */

max_stars_repo_name: JonasToth/depth-conversions

max_issues_repo_name: JonasToth/depth-conversions

max_forks_repo_name: JonasToth/depth-conversions


text:
#include "fg_interface.h"
#include <gsl/gsl_rng.h>
#include <math.h>
#include <stdio.h>

int main(int argc, char *argv[]) {
    fg_func test = fg_init_7_4096(log, 1E-15, 1000, 1E-12, 1E-15, 0);

    int n_el = 100000;
    gsl_rng *r = gsl_rng_alloc(gsl_rng_taus);
    gsl_rng_set(r, 1);

    double *x = (double *)malloc(n_el * sizeof(double));

    for (int i = 0; i < n_el; ++i) {
        x[i] = gsl_rng_uniform(r) * 1000;
    }

    double sum = 0.0;
    for (int i_loop = 0; i_loop < 5000; ++i_loop) {
        for (int i = 0; i < n_el; i++) {
            sum += fg_eval(&test, x[i]);
        }
    }
    printf("%g\n", sum);

    return 0;
}