Floating point inaccuracies breaks idempotence in weird ways

Viewed 126

I was wondering how bad floating point inaccuracies can get when normalizing vectors, specifically a 3D one. Normalizing once should give the same value as normalizing twice, but, as expected, floating point inaccuracies don't guarantee us that. So I decided to test some other things: if you normalize indefinitely, will it always "converge" to a value that normalizes to itself? The answer is no. For example, a few seconds testing and I found a value that loops between 2: -0.60445204839058564 -0.54547899327834748 -0.58059485795902943 and -0.60445204839058575 -0.54547899327834759 -0.58059485795902954. Some "converge" after 3, 5 steps. Here's my boring code:

#include <iostream>
#include <numeric>
#include <random>

using namespace std;

struct Point
{
    double x, y, z;
    bool operator==(const Point& o) const
    {
        return
            x == o.x &&
            y == o.y &&
            z == o.z;
    }
};

auto square(auto x)
{
    return x * x;
}

Point normalize(const Point& p)
{
    auto l = sqrt(square(p.x) + square(p.y) + square(p.z));
    return { p.x / l, p.y / l, p.z / l };
}

std::mt19937 mt(24);
std::uniform_real_distribution d(-1., 1.);

int main()
{
    while (true)
    {
        auto
            x = d(mt),
            y = d(mt),
            z = d(mt);

        Point p{ x, y, z };

        auto
            n1 = normalize(p),
            n2 = normalize(normalize(p));
        int cnt = 1;
        while (n1 != n2)
        {
            n1 = n2;
            cnt++;
            n2 = normalize(n2);
            auto len = square(n2.x) + square(n2.y) + square(n2.z);
            if (len != 1.)
            {
                __debugbreak();
            }
        }
        if (cnt != 2)
        {
            cout << x << ' ' << y << ' ' << z
                << " took " << cnt << '\n';
        }
    }
}

So I have a few curiosities:

  1. What inputs create the longest loop/cycle and never converge?
  2. What inputs take longest to converge?
  3. After normalizing once, what input gives the maximum error from a length of 1?
0 Answers
Related