Golomb Triangle Coordinates

[Tom Sirgedas]

I wrote some code to search for Golomb triangles, given 3 points. A Golomb triangle is a unit-distance equilateral triangle, where each of its 3 vertices is a unit distance away from 3 given points. (In total there are 6 points, and 6 unit distances involved).

In the picture below, the 3 given points are shown as small circles, and two Golomb triangles are drawn.

Sometimes there are a lot of solutions. It seems that up to 12 are possible, and they come in pairs.

[Tom Sirgedas]

My approach:

Start at point `A`. We'll search for an `angle`. The `angle` determines where `a` is. Now there are at most two possibilities for `b`. Next we find `c` based on `a` and `b` (2 locations are possible).

Every distance is now correct, except between `c` and `C`.

`f( angle, which ) = distance(c,C) - 1` is a function that quantifies the error. The zeroes of this function correspond to solutions. `which` is a parameter with just 4 values to indicate which "possibilities" should be chosen in the previous step.

My simple technique for finding the zeroes is to step through the angles .01 radians at a time and find sign-changes in `f`. When a sign change is detected, do a binary search to find the `angle` that makes `f` zero.

This can miss solutions where `f` only briefly changes signs, or solutions near where `f` is undefined. It works fairly well in practice, though. A simple improvement might be to scan over angles coarsely at first, and then make a finer scan over chaotic regions.

[Tom Sirgedas]

here's the relevant code:

const double PI = acos( 0. )*2;


// given `f(lo)` < 0 and `f(hi)` > 0, find `x` such that `f(x)=0`
double binarySearchRoot_( double lo, double hi, std::function<double( double )> f )
{
   double mid = (lo+hi)/2;
   if ( lo == mid || hi == mid )
      return mid;
   double fmid = f( mid );
   return fmid < 0 ? binarySearchRoot_( mid, hi, f ) : binarySearchRoot_( lo, mid, f );
}


bool binarySearchRoot( double lo, double hi, double& root, std::function<double( double )> f )
{
   double flo = f( lo );
   double fhi = f( hi );
   bool sameSign = (flo < 0) == (fhi < 0);
   if ( sameSign )
      return false;
   if ( flo > fhi )
      std::swap( lo, hi );
   root = binarySearchRoot_( lo, hi, f );
   return abs( f( root ) ) < 1e-6;
}


// given 2-d points A, B, C
// find 2-d points a, b, c such that the pairs (A,a),(B,b),(C,c),(a,b),(b,c),(c,a) are all 1 unit apart
class GolombFinder
{
public:
   struct Solution
   {
      XYf a, b, c;
   };


public:
   GolombFinder( const XYf& A, const XYf& B, const XYf& C ) : _A( A ), _B( B ), _C( C ) {}


   // return a point a unit away from both `p` and `q`
   // the last bit of `which` controls which of the 2 solutions to return
   XYf findPointDist1From( const XYf& p, const XYf& q, int which, bool& success ) const
   {
      XYf mid = (p+q)/2;
      XYf perpendicular = XYf( p.y - q.y, q.x - p.x );
      double dist = (p-q).len();
      success = dist <= 2;
      int sign = which&1 ? 1 : -1;
      return mid + perpendicular * (sign * sqrt( 1 - dist*dist/4 ) / dist);
   }




   // `angle` is the direction from `_A` to `a`
   // `which` is one of {0,1,2,3} -- four separate scenarios
   bool solveGivenAngle( double angle, int which, Solution& sol ) const
   {
      sol.a = _A + XYf( cos( angle ), sin( angle ) );


      // `b` is 1 unit away from `a` and `_B`
      bool success;
      sol.b = findPointDist1From( sol.a, _B, which, success );
      if ( !success )
         return false;


      // `c` is 1 unit away from `a` and `b`
      sol.c = findPointDist1From( sol.a, sol.b, which/2, success );
      return success;
   }


   bool eval( double angle, int which, double& result ) const
   {
      Solution sol;
      if ( !solveGivenAngle( angle, which, sol ) )
         return false;


      result = (_C-sol.c).len() - 1; // the distance between `_C` and `c` should be 1 -- return the error amount
      return true;
   }


   std::vector<double> findAngles( int which, double angleStepSize )
   {
      std::vector<double> ret;
      for ( double loAngle = 0; loAngle < 2*PI; loAngle += angleStepSize )
      {
         double hiAngle = loAngle + angleStepSize;
         double root;
         bool foundRoot = binarySearchRoot( loAngle, std::min( hiAngle, 2*PI ), root, [&]( double angle ) {
            double f;
            bool success = eval( angle, which, f );
            return f;
         } );
         if ( foundRoot )
            ret.push_back( root );
      }
      return ret;
   }


   bool isUnitDistance( double dist ) { return fabs( dist-1 ) < 1e-7; }
   bool isValidSolution( const Solution& sol )
   {
      return isUnitDistance( (sol.a - _A).len() )
          && isUnitDistance( (sol.b - _B).len() )
          && isUnitDistance( (sol.c - _C).len() )
          && isUnitDistance( (sol.a - sol.b).len() )
          && isUnitDistance( (sol.b - sol.c).len() )
          && isUnitDistance( (sol.c - sol.a).len() );
   }


   std::vector<Solution> findSolutions()
   {
      const double angleStepSize = .005; // controls how granularly to search for solutions


      std::vector<Solution> ret;
      for ( int which = 0; which < 4; which++ )
      {
         std::vector<double> angles = findAngles( which, angleStepSize );
         for ( double angle : angles )
         {
            Solution sol;
            if ( !solveGivenAngle( angle, which, sol ) )
               continue;


            if ( !isValidSolution( sol ) )
               throw std::runtime_error( "invalid solution" );
            ret.push_back( sol );
         }
      }
      return ret;
   }


private:
   XYf _A, _B, _C;
};




std::vector<std::vector<XYf>> golomb( const XYf& A, const XYf& B, const XYf& C )
{
   GolombFinder finder( A, B, C );
   std::vector<GolombFinder::Solution> sols = finder.findSolutions();


   std::vector<std::vector<XYf>> ret;
   for ( auto sol : sols )
      ret.push_back( { sol.a, sol.b, sol.c } );  


   return ret;
}

[Tom Sirgedas]

Here's a video that shows the Golomb Triangles as I drag the 3 given points around: https://www.youtube.com/watch?v=aCRrzjDD5Bw

[Tom Sirgedas]

BTW, this was motivated by the discussion here: https://dustingmixon.wordpress.com/2019/08/05/polymath16-fourteenth-thread-automated-graph-minimization/#comment-24410