Кирпичи Эйлера. Вторая итерация

Иван в комментариях к предыдущему посту придумал очень быстрый и эффективный способ вычисления пифагоровых троек к нечетному числу. Мне стало интересно, я тоже немного погрузился в тему статьи Матсона.

Правда, той эффективности, что у Ивана, мне достичь не удалось, я пока в процессе. И в этом процессе я обнаружил любопытную вещь.
Мне скучно просто искать совершенный кубоид, поэтому  я решил в процессе еще и сохранять все найденные обычные кирпичи Эйлера. Это чуть более сложная задача, потому должна работать чуть медленнее, чем поиск совершенного кубоида. И для её работы нужно уметь эффективно считать целочисленный корень для больших целых чисел, в моем случае 128-битных.
Я решил проверить, сможет ли FPU (или по-другому x87) выполнить эту задачу. Шансы у него были, так как во внутреннем представлении FPU хранит вещественные в расширенном формате, длиной не стандартные 64, а целых 80 бит.
В процессе обнаружил, что точности как бы не хватает. Начал копать и через непродолжительное время выяснил, что по какой-то странной причине точность FPU оказалась занижена до 64 бит. Да, в FPU с помощью флагов можно менять точность внутреннего представления чисел от 32 до 80 бит, что кроме точности влияет и на производительность, понятное дело.
В документации на сайте AMD утверждается, что по умолчанию FPU установлен на максимальную точность. Но возможно, это было давно и в современных версиях процессора уже не действует? Кроме того, точность может менять Windows, а так же компилятор Delphi. Кто из них повинен в такой ситуации, мне не ведомо.
Интереснее другое: зачем это вообще сделано? Лично я вижу лишь одну причину: что бы пользователи не заметили разницу в решениях своих задач на FPU и SSE, иначе кто-нибудь как-нибудь да непременно выскажет претензии. Но это лишь предположение.

Что касается целочисленного корня. При максимальной точности FPU может находить точное целое значения корня из 128-битового беззнакового целого. Но вот точности, что бы определить, целым ли получился результат, уже не хватает, а было бы очень кстати.

Теперь про мою реализацию. Во-первых, я отказался от проверки деления на 19. Конечно, такое деление снижает количество рассматриваемых вариантов, но, если я правильно понял, не так уж значительно по сравнению с делением на 16. В то же время, проверка делимости на 16 происходит гораздо быстрее.
Но возможно здесь я ошибаюсь, так как проверка каждого варианта на то, что он является кирпичом Эйлера стоит гораздо дороже, чем разница в скорости проверок делимости.

Во-вторых, я пока не сделал сохранение найденных делителей для последующего использования, что, конечно, существенно снижает производительность.

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

Все это, конечно, весьма далеко от  гораздо более эффективной реализации Ивана. Но, тем не менее, тоже показывает относительно неплохую производительность: от 100001 до 100000001 по нечетной стороне все варианты кирпичей Эйлера были найдены за час с небольшим, всего 3182002 варианта.
Я решил, что вместо того, что бы проверять на первичность решения в процессе поиска (мне показалось это достаточно затратной операцией), сохранять все найденные варианты. А уже потом быстренько отсеять все составные, ну и плюс тут же проверить и на совершенность.