Решение проблем округления с плавающей запятой C ++

Я разрабатываю научное приложение (моделирование движения хромосом в ядре клетки). Хромосомы разделены на небольшие фрагменты, которые вращаются вокруг случайной оси с использованием матриц вращения 4x4.

Проблема в том, что моделирование выполняет сотни миллиардов поворотов, поэтому ошибки округления с плавающей запятой складываются и растут экспоненциально, поэтому фрагменты имеют тенденцию «уплывать» и отделяться от остальной хромосомы с течением времени.

Я использую двойную точность в C ++. На данный момент программа работает на ЦП, но будет портирована для CUDA, и моделирование может длиться максимум 1 месяц.

Я понятия не имею, как я мог бы каким-то образом перенормировать хромосому, потому что все фрагменты связаны вместе (вы можете видеть это как двусвязный список), но я думаю, что это будет лучшая идея, если возможно.

Есть ли у вас предложения? Я чувствую себя немного потерянным.

Большое спасибо,

Х.

РЕДАКТИРОВАТЬ: Добавлен упрощенный пример кода. Вы можете предположить, что вся матричная математика является классической реализацией.

// Rotate 1000000 times
for (int i = 0; i < 1000000; ++i)
{
    // Pick a random section start
    int istart = rand() % chromosome->length;

    // Pick the end 20 segments further (cyclic)
    int iend = (istart + 20) % chromosome->length;

    // Build rotation axis
    Vector4 axis = chromosome->segments[istart].position - chromosome->segments[iend].position;
    axis.normalize();

    // Build rotation matrix and translation vector
    Matrix4 rotm(axis, rand() / float(RAND_MAX));
    Vector4 oldpos = chromosome->segments[istart].position;

    // Rotate each segment between istart and iend using rotm
    for (int j = (istart + 1) % chromosome->length; j != iend; ++j, j %= chromosome->length)
    {
        chromosome->segments[j].position -= oldpos;
        chromosome->segments[j].position.transform(rotm);
        chromosome->segments[j].position += oldpos;
    }
}