Skip to content

Add SLEEF acceleration for trigonometric and power functions #595

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
bhaller merged 6 commits into from
Dec 18, 2025
Merged

Add SLEEF acceleration for trigonometric and power functions #595

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
64643b0
Add SLEEF acceleration for trigonometric functions
andrewkern Dec 17, 2025
4019376
Add SLEEF acceleration for pow() function
andrewkern Dec 17, 2025
dc3c4a3
update benchmark script
andrewkern Dec 17, 2025
17d0906
Fix whitespace: add tabs to blank lines inside functions
andrewkern Dec 18, 2025
dad94b7
forgot to alter tests to use allClose
andrewkern Dec 18, 2025
9d16245
Update VERSIONS for new SIMD functions
andrewkern Dec 18, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions VERSIONS
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -28,6 +28,7 @@ development head (in the master branch):
SIMD optimizations for spatial interaction strength calculations, thanks to Andy Kern, https://github.com/MesserLab/SLiM/pull/590
fix #593, Student's T distribution had a sign error in tdist(), thanks to Andy Kern, https://github.com/MesserLab/SLiM/issues/591
add comprehensive tests for spatial interaction kernel calculations; includes SIMD vs scalar consistency tests, C++ level tests, and script-level tests, thanks to Andy Kern, https://github.com/MesserLab/SLiM/pull/592
add SIMD acceleration for trig functions (sin(), cos(), tan(), asin(), acos(), atan(), atan2()), power functions (^ opperator, pow, SpatialMap.pow()), thanks to Andy Kern, https://github.com/MesserLab/SLiM/pull/595


version 5.1 (Eidos version 4.1):
Expand Down
18 changes: 9 additions & 9 deletions core/slim_test_other.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1668,9 +1668,9 @@ void _RunSpatialMapTests(void)
SLiMAssertScriptStop(prefix_1D + "m1.divide(mv2); if (identical(mv1 / mv2, m1.gridValues())) stop(); } ");
SLiMAssertScriptStop(prefix_1D + "m1.divide(m2); if (identical(mv1 / mv2, m1.gridValues())) stop(); } ");

SLiMAssertScriptStop(prefix_1D + "m1.power(0.25); if (identical(mv1 ^ 0.25, m1.gridValues())) stop(); } ");
SLiMAssertScriptStop(prefix_1D + "m1.power(mv2); if (identical(mv1 ^ mv2, m1.gridValues())) stop(); } ");
SLiMAssertScriptStop(prefix_1D + "m1.power(m2); if (identical(mv1 ^ mv2, m1.gridValues())) stop(); } ");
SLiMAssertScriptStop(prefix_1D + "m1.power(0.25); if (allClose(mv1 ^ 0.25, m1.gridValues())) stop(); } ");
SLiMAssertScriptStop(prefix_1D + "m1.power(mv2); if (allClose(mv1 ^ mv2, m1.gridValues())) stop(); } ");
SLiMAssertScriptStop(prefix_1D + "m1.power(m2); if (allClose(mv1 ^ mv2, m1.gridValues())) stop(); } ");

SLiMAssertScriptStop(prefix_1D + "m1.exp(); if (allClose(exp(mv1), m1.gridValues())) stop(); } ");

Expand Down Expand Up @@ -1783,9 +1783,9 @@ void _RunSpatialMapTests(void)
SLiMAssertScriptStop(prefix_2D + "m1.divide(mv2); if (identical(mv1 / mv2, m1.gridValues())) stop(); } ");
SLiMAssertScriptStop(prefix_2D + "m1.divide(m2); if (identical(mv1 / mv2, m1.gridValues())) stop(); } ");

SLiMAssertScriptStop(prefix_2D + "m1.power(0.25); if (identical(mv1 ^ 0.25, m1.gridValues())) stop(); } ");
SLiMAssertScriptStop(prefix_2D + "m1.power(mv2); if (identical(mv1 ^ mv2, m1.gridValues())) stop(); } ");
SLiMAssertScriptStop(prefix_2D + "m1.power(m2); if (identical(mv1 ^ mv2, m1.gridValues())) stop(); } ");
SLiMAssertScriptStop(prefix_2D + "m1.power(0.25); if (allClose(mv1 ^ 0.25, m1.gridValues())) stop(); } ");
SLiMAssertScriptStop(prefix_2D + "m1.power(mv2); if (allClose(mv1 ^ mv2, m1.gridValues())) stop(); } ");
SLiMAssertScriptStop(prefix_2D + "m1.power(m2); if (allClose(mv1 ^ mv2, m1.gridValues())) stop(); } ");

SLiMAssertScriptStop(prefix_2D + "m1.exp(); if (allClose(exp(mv1), m1.gridValues())) stop(); } ");

Expand Down Expand Up @@ -1899,9 +1899,9 @@ void _RunSpatialMapTests(void)
SLiMAssertScriptStop(prefix_3D + "m1.divide(mv2); if (identical(mv1 / mv2, m1.gridValues())) stop(); } ");
SLiMAssertScriptStop(prefix_3D + "m1.divide(m2); if (identical(mv1 / mv2, m1.gridValues())) stop(); } ");

SLiMAssertScriptStop(prefix_3D + "m1.power(0.25); if (identical(mv1 ^ 0.25, m1.gridValues())) stop(); } ");
SLiMAssertScriptStop(prefix_3D + "m1.power(mv2); if (identical(mv1 ^ mv2, m1.gridValues())) stop(); } ");
SLiMAssertScriptStop(prefix_3D + "m1.power(m2); if (identical(mv1 ^ mv2, m1.gridValues())) stop(); } ");
SLiMAssertScriptStop(prefix_3D + "m1.power(0.25); if (allClose(mv1 ^ 0.25, m1.gridValues())) stop(); } ");
SLiMAssertScriptStop(prefix_3D + "m1.power(mv2); if (allClose(mv1 ^ mv2, m1.gridValues())) stop(); } ");
SLiMAssertScriptStop(prefix_3D + "m1.power(m2); if (allClose(mv1 ^ mv2, m1.gridValues())) stop(); } ");

SLiMAssertScriptStop(prefix_3D + "m1.exp(); if (allClose(exp(mv1), m1.gridValues())) stop(); } ");

Expand Down
6 changes: 2 additions & 4 deletions core/spatial_map.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1654,8 +1654,7 @@ EidosValue_SP SpatialMap::ExecuteMethod_power(EidosGlobalStringID p_method_id, c
double power_scalar = x_value->NumericAtIndex_NOCAST(0, nullptr);

// FIXME: TO BE PARALLELIZED
for (int64_t i = 0; i < values_size_; ++i)
values_[i] = pow(values_[i], power_scalar);
Eidos_SIMD::pow_float64_scalar_exp(values_, power_scalar, values_, values_size_);
}
else
{
Expand All @@ -1666,8 +1665,7 @@ EidosValue_SP SpatialMap::ExecuteMethod_power(EidosGlobalStringID p_method_id, c
EIDOS_TERMINATION << "ERROR (SpatialMap::ExecuteMethod_power): power() requires the target SpatialMap to be compatible with the SpatialMap supplied in x (using the same spatiality and bounds, and having the same grid resolution)." << EidosTerminate();

// FIXME: TO BE PARALLELIZED
for (int64_t i = 0; i < values_size_; ++i)
values_[i] = pow(values_[i], power_map_values[i]);
Eidos_SIMD::pow_float64(values_, power_map_values, values_, values_size_);
}

_ValuesChanged();
Expand Down
150 changes: 124 additions & 26 deletions eidos/eidos_functions_math.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -111,12 +111,26 @@ EidosValue_SP Eidos_ExecuteFunction_acos(const std::vector<EidosValue_SP> &p_arg
EidosValue_SP result_SP(nullptr);

EidosValue *x_value = p_arguments[0].get();
EidosValueType x_type = x_value->Type();
int x_count = x_value->Count();
EidosValue_Float *float_result = (new (gEidosValuePool->AllocateChunk()) EidosValue_Float())->resize_no_initialize(x_count);
result_SP = EidosValue_SP(float_result);

for (int value_index = 0; value_index < x_count; ++value_index)
float_result->set_float_no_check(acos(x_value->NumericAtIndex_NOCAST(value_index, nullptr)), value_index);
if (x_type == EidosValueType::kValueInt)
{
EidosValue_Float *float_result = (new (gEidosValuePool->AllocateChunk()) EidosValue_Float())->resize_no_initialize(x_count);
result_SP = EidosValue_SP(float_result);

for (int value_index = 0; value_index < x_count; ++value_index)
float_result->set_float_no_check(acos(x_value->NumericAtIndex_NOCAST(value_index, nullptr)), value_index);
}
else if (x_type == EidosValueType::kValueFloat)
{
const double *float_data = x_value->FloatData();
EidosValue_Float *float_result = (new (gEidosValuePool->AllocateChunk()) EidosValue_Float())->resize_no_initialize(x_count);
double *float_result_data = float_result->data_mutable();
result_SP = EidosValue_SP(float_result);

Eidos_SIMD::acos_float64(float_data, float_result_data, x_count);
}

result_SP->CopyDimensionsFromValue(x_value);

Expand All @@ -129,12 +143,26 @@ EidosValue_SP Eidos_ExecuteFunction_asin(const std::vector<EidosValue_SP> &p_arg
EidosValue_SP result_SP(nullptr);

EidosValue *x_value = p_arguments[0].get();
EidosValueType x_type = x_value->Type();
int x_count = x_value->Count();
EidosValue_Float *float_result = (new (gEidosValuePool->AllocateChunk()) EidosValue_Float())->resize_no_initialize(x_count);
result_SP = EidosValue_SP(float_result);

for (int value_index = 0; value_index < x_count; ++value_index)
float_result->set_float_no_check(asin(x_value->NumericAtIndex_NOCAST(value_index, nullptr)), value_index);
if (x_type == EidosValueType::kValueInt)
{
EidosValue_Float *float_result = (new (gEidosValuePool->AllocateChunk()) EidosValue_Float())->resize_no_initialize(x_count);
result_SP = EidosValue_SP(float_result);

for (int value_index = 0; value_index < x_count; ++value_index)
float_result->set_float_no_check(asin(x_value->NumericAtIndex_NOCAST(value_index, nullptr)), value_index);
}
else if (x_type == EidosValueType::kValueFloat)
{
const double *float_data = x_value->FloatData();
EidosValue_Float *float_result = (new (gEidosValuePool->AllocateChunk()) EidosValue_Float())->resize_no_initialize(x_count);
double *float_result_data = float_result->data_mutable();
result_SP = EidosValue_SP(float_result);

Eidos_SIMD::asin_float64(float_data, float_result_data, x_count);
}

result_SP->CopyDimensionsFromValue(x_value);

Expand All @@ -147,12 +175,26 @@ EidosValue_SP Eidos_ExecuteFunction_atan(const std::vector<EidosValue_SP> &p_arg
EidosValue_SP result_SP(nullptr);

EidosValue *x_value = p_arguments[0].get();
EidosValueType x_type = x_value->Type();
int x_count = x_value->Count();
EidosValue_Float *float_result = (new (gEidosValuePool->AllocateChunk()) EidosValue_Float())->resize_no_initialize(x_count);
result_SP = EidosValue_SP(float_result);

for (int value_index = 0; value_index < x_count; ++value_index)
float_result->set_float_no_check(atan(x_value->NumericAtIndex_NOCAST(value_index, nullptr)), value_index);
if (x_type == EidosValueType::kValueInt)
{
EidosValue_Float *float_result = (new (gEidosValuePool->AllocateChunk()) EidosValue_Float())->resize_no_initialize(x_count);
result_SP = EidosValue_SP(float_result);

for (int value_index = 0; value_index < x_count; ++value_index)
float_result->set_float_no_check(atan(x_value->NumericAtIndex_NOCAST(value_index, nullptr)), value_index);
}
else if (x_type == EidosValueType::kValueFloat)
{
const double *float_data = x_value->FloatData();
EidosValue_Float *float_result = (new (gEidosValuePool->AllocateChunk()) EidosValue_Float())->resize_no_initialize(x_count);
double *float_result_data = float_result->data_mutable();
result_SP = EidosValue_SP(float_result);

Eidos_SIMD::atan_float64(float_data, float_result_data, x_count);
}

result_SP->CopyDimensionsFromValue(x_value);

Expand All @@ -166,8 +208,10 @@ EidosValue_SP Eidos_ExecuteFunction_atan2(const std::vector<EidosValue_SP> &p_ar

EidosValue *x_value = p_arguments[0].get();
int x_count = x_value->Count();
EidosValueType x_type = x_value->Type();
EidosValue *y_value = p_arguments[1].get();
int y_count = y_value->Count();
EidosValueType y_type = y_value->Type();

if (x_count != y_count)
EIDOS_TERMINATION << "ERROR (Eidos_ExecuteFunction_atan2): function atan2() requires arguments of equal length." << EidosTerminate(nullptr);
Expand All @@ -183,8 +227,20 @@ EidosValue_SP Eidos_ExecuteFunction_atan2(const std::vector<EidosValue_SP> &p_ar
EidosValue_Float *float_result = (new (gEidosValuePool->AllocateChunk()) EidosValue_Float())->resize_no_initialize(x_count);
result_SP = EidosValue_SP(float_result);

for (int value_index = 0; value_index < x_count; ++value_index)
float_result->set_float_no_check(atan2(x_value->NumericAtIndex_NOCAST(value_index, nullptr), y_value->NumericAtIndex_NOCAST(value_index, nullptr)), value_index);
// Use SIMD when both operands are float type
if ((x_type == EidosValueType::kValueFloat) && (y_type == EidosValueType::kValueFloat))
{
const double *x_data = x_value->FloatData();
const double *y_data = y_value->FloatData();
double *float_result_data = float_result->data_mutable();

Eidos_SIMD::atan2_float64(x_data, y_data, float_result_data, x_count);
}
else
{
for (int value_index = 0; value_index < x_count; ++value_index)
float_result->set_float_no_check(atan2(x_value->NumericAtIndex_NOCAST(value_index, nullptr), y_value->NumericAtIndex_NOCAST(value_index, nullptr)), value_index);
}

// Copy dimensions from whichever operand we chose at the beginning
result_SP->CopyDimensionsFromValue(result_dim_source.get());
Expand Down Expand Up @@ -225,12 +281,26 @@ EidosValue_SP Eidos_ExecuteFunction_cos(const std::vector<EidosValue_SP> &p_argu
EidosValue_SP result_SP(nullptr);

EidosValue *x_value = p_arguments[0].get();
EidosValueType x_type = x_value->Type();
int x_count = x_value->Count();
EidosValue_Float *float_result = (new (gEidosValuePool->AllocateChunk()) EidosValue_Float())->resize_no_initialize(x_count);
result_SP = EidosValue_SP(float_result);

for (int value_index = 0; value_index < x_count; ++value_index)
float_result->set_float_no_check(cos(x_value->NumericAtIndex_NOCAST(value_index, nullptr)), value_index);
if (x_type == EidosValueType::kValueInt)
{
EidosValue_Float *float_result = (new (gEidosValuePool->AllocateChunk()) EidosValue_Float())->resize_no_initialize(x_count);
result_SP = EidosValue_SP(float_result);

for (int value_index = 0; value_index < x_count; ++value_index)
float_result->set_float_no_check(cos(x_value->NumericAtIndex_NOCAST(value_index, nullptr)), value_index);
}
else if (x_type == EidosValueType::kValueFloat)
{
const double *float_data = x_value->FloatData();
EidosValue_Float *float_result = (new (gEidosValuePool->AllocateChunk()) EidosValue_Float())->resize_no_initialize(x_count);
double *float_result_data = float_result->data_mutable();
result_SP = EidosValue_SP(float_result);

Eidos_SIMD::cos_float64(float_data, float_result_data, x_count);
}

result_SP->CopyDimensionsFromValue(x_value);

Expand Down Expand Up @@ -2429,12 +2499,26 @@ EidosValue_SP Eidos_ExecuteFunction_sin(const std::vector<EidosValue_SP> &p_argu
EidosValue_SP result_SP(nullptr);

EidosValue *x_value = p_arguments[0].get();
EidosValueType x_type = x_value->Type();
int x_count = x_value->Count();
EidosValue_Float *float_result = (new (gEidosValuePool->AllocateChunk()) EidosValue_Float())->resize_no_initialize(x_count);
result_SP = EidosValue_SP(float_result);

for (int value_index = 0; value_index < x_count; ++value_index)
float_result->set_float_no_check(sin(x_value->NumericAtIndex_NOCAST(value_index, nullptr)), value_index);
if (x_type == EidosValueType::kValueInt)
{
EidosValue_Float *float_result = (new (gEidosValuePool->AllocateChunk()) EidosValue_Float())->resize_no_initialize(x_count);
result_SP = EidosValue_SP(float_result);

for (int value_index = 0; value_index < x_count; ++value_index)
float_result->set_float_no_check(sin(x_value->NumericAtIndex_NOCAST(value_index, nullptr)), value_index);
}
else if (x_type == EidosValueType::kValueFloat)
{
const double *float_data = x_value->FloatData();
EidosValue_Float *float_result = (new (gEidosValuePool->AllocateChunk()) EidosValue_Float())->resize_no_initialize(x_count);
double *float_result_data = float_result->data_mutable();
result_SP = EidosValue_SP(float_result);

Eidos_SIMD::sin_float64(float_data, float_result_data, x_count);
}

result_SP->CopyDimensionsFromValue(x_value);

Expand Down Expand Up @@ -2619,12 +2703,26 @@ EidosValue_SP Eidos_ExecuteFunction_tan(const std::vector<EidosValue_SP> &p_argu
EidosValue_SP result_SP(nullptr);

EidosValue *x_value = p_arguments[0].get();
EidosValueType x_type = x_value->Type();
int x_count = x_value->Count();
EidosValue_Float *float_result = (new (gEidosValuePool->AllocateChunk()) EidosValue_Float())->resize_no_initialize(x_count);
result_SP = EidosValue_SP(float_result);

for (int value_index = 0; value_index < x_count; ++value_index)
float_result->set_float_no_check(tan(x_value->NumericAtIndex_NOCAST(value_index, nullptr)), value_index);
if (x_type == EidosValueType::kValueInt)
{
EidosValue_Float *float_result = (new (gEidosValuePool->AllocateChunk()) EidosValue_Float())->resize_no_initialize(x_count);
result_SP = EidosValue_SP(float_result);

for (int value_index = 0; value_index < x_count; ++value_index)
float_result->set_float_no_check(tan(x_value->NumericAtIndex_NOCAST(value_index, nullptr)), value_index);
}
else if (x_type == EidosValueType::kValueFloat)
{
const double *float_data = x_value->FloatData();
EidosValue_Float *float_result = (new (gEidosValuePool->AllocateChunk()) EidosValue_Float())->resize_no_initialize(x_count);
double *float_result_data = float_result->data_mutable();
result_SP = EidosValue_SP(float_result);

Eidos_SIMD::tan_float64(float_data, float_result_data, x_count);
}

result_SP->CopyDimensionsFromValue(x_value);

Expand Down
13 changes: 7 additions & 6 deletions eidos/eidos_interpreter.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -23,6 +23,7 @@
#include "eidos_ast_node.h"
#include "eidos_rng.h"
#include "eidos_call_signature.h"
#include "eidos_simd.h"

#include <sstream>
#include <stdexcept>
Expand Down Expand Up @@ -3211,9 +3212,9 @@ EidosValue_SP EidosInterpreter::Evaluate_Exp(const EidosASTNode *p_node)
{
const double *first_child_data = first_child_value->FloatData();
const double *second_child_data = second_child_value->FloatData();
double *float_result_data = float_result->data_mutable();

for (int value_index = 0; value_index < first_child_count; ++value_index)
float_result->set_float_no_check(pow(first_child_data[value_index], second_child_data[value_index]), value_index);
Eidos_SIMD::pow_float64(first_child_data, second_child_data, float_result_data, first_child_count);
}
else if ((first_child_type == EidosValueType::kValueFloat) && (second_child_type == EidosValueType::kValueInt))
{
Expand Down Expand Up @@ -3258,9 +3259,9 @@ EidosValue_SP EidosInterpreter::Evaluate_Exp(const EidosASTNode *p_node)
else // (second_child_type == EidosValueType::kValueFloat)
{
const double *second_child_data = second_child_value->FloatData();
double *float_result_data = float_result->data_mutable();

for (int value_index = 0; value_index < second_child_count; ++value_index)
float_result->set_float_no_check(pow(singleton_float, second_child_data[value_index]), value_index);
Eidos_SIMD::pow_float64_scalar_base(singleton_float, second_child_data, float_result_data, second_child_count);
}

result_SP = float_result_SP;
Expand All @@ -3281,9 +3282,9 @@ EidosValue_SP EidosInterpreter::Evaluate_Exp(const EidosASTNode *p_node)
else // (first_child_type == EidosValueType::kValueFloat)
{
const double *first_child_data = first_child_value->FloatData();
double *float_result_data = float_result->data_mutable();

for (int value_index = 0; value_index < first_child_count; ++value_index)
float_result->set_float_no_check(pow(first_child_data[value_index], singleton_float), value_index);
Eidos_SIMD::pow_float64_scalar_exp(first_child_data, singleton_float, float_result_data, first_child_count);
}

result_SP = float_result_SP;
Expand Down
Loading
Loading