TweetFollow Us on Twitter

Apr 99 Challenge

Volume Number: 15 (1999)
Issue Number: 4
Column Tag: Programmer's Challenge

Apr 99 Challenge

by Bob Boonstra, Westford, MA

Shortest Network

This month's problem was suggested by Michael Kennedy, who wins two Challenge points for making the suggestion. The problem is to find the shortest network of line segments interconnecting a specified set of points. Shortest network algorithms have obvious practical application in constructing transportation and communications networks. In a January 1989, Scientific American article, Marshall Bern and Ronald Graham discussed the shortest network "Steiner" problem as one of a class of NP-hard problems. While no polynomial-time algorithm is known, the article (which, unfortunately, I have not been able to find online) discusses practical algorithms that produce networks slightly longer than the optimal one. Your Challenge for this month is to produce a near-optimal network in minimum time. Fortunately, we have been granted unlimited power of eminent domain, so there are no restrictions on where intermediate nodes may be placed or where connections may be routed.

The prototype for the code you should write is:

#if defined(__cplusplus)
extern "C" {
#endif

typedef struct Node {   /* node coordinates */
   double x;
   double y;
} Node;

typedef struct Connection {
         /* connection between Node[index1] and Node[index2] */
   long index1;
   long index2;
} Connection;

long /* numConnections */ ShortestNetwork(
   long numInitialNodes,         /* number of nodes to connect */
   long *numIntermediateNodes,   /* number of nodes added by ShortestNetwork */
   Node nodes[], 
      /* Nodes 0..numInitialNodes-1 are initialized on entry. */
      /* Nodes numInitialNodes..numInitialNodes+*numIntermediateNodes 
                  are added by ShortestNetwork */
   Connection connections[],   /* connections between nodes */
   long maxNodes,        /* number of entries allocated for nodes */
   long maxConnections   /* number of entries allocated for connections */
);

#if defined(__cplusplus)
}
#endif

Your ShortestNetwork routine will be given a list of numInitialNodes nodes to connect. You may add intermediate nodes to help you form a shorter network, and must produce as output a list of connections between pairs of nodes. The connections must provide a path between any pair of the initial nodes.

Your solution must return the number of intermediate nodes added to the network in *numIntermediateNodes, while storing the location of those nodes in nodes[numInitialNodes+k], k=0..*numIntermediateNodes-1. A connection is specified by storing the indices of the two nodes being connected into the connection array. Your ShortestNetwork routine should return the number of connections created.

The maxNodes and maxConnections parameters indicate how much storage has been allocated for nodes and connections. It is my intention to allocate enough storage for all the nodes and connections your solution might create, but if it turns out that there is not enough storage, your solution should return a value of -1 to indicate that storage was exhausted.

The winner will be the solution that generates the shortest network in the minimum amount of time. Specifically, your solution will be assigned a cost equal to the sum of the distances between nodes in your list of connections, plus a penalty of 10% for each second of execution time. Solutions that do not connect all of the initial nodes will be penalized with a very large cost. The solution with the lowest total cost over a series of networking problems will be the winner.

This will be a native PowerPC Challenge, using the latest CodeWarrior environment. Solutions may be coded in C, C++, or Pascal. Thanks to Michael for suggesting this Challenge.

Three Months Ago Winner

Congratulations to Tom Saxton for submitting the winning solution to the January Sphere Packing Challenge. You may recall that this Challenge was to pack a set of spheres of varying size into a box with minimum volume, and to do so in the shortest amount of time possible. Tom submitted one of only two solutions received for this Challenge, and his was the only one that performed correctly.

Tom's approach is to decide on a footprint for the box to contain the spheres, "drop" the spheres individually into the box until they hit another sphere or the bottom of the box, while attempting to move the dropped sphere around the obstacle without going outside the box footprint. The solution then iterates with random movements to try to converge to a better solution. Tom observed in his submission that the time penalty for this problem (1% per millisecond of execution time) was very severe, making it unproductive to let his algorithm iterate very long. Every tenth of a second of execution time requires a factor of 2 reduction in volume to be productive, a rate of improvement smaller than what Tom was able to achieve.

I evaluated the solutions using six test cases with between 200 and 2000 spheres per test case. As one might expect, execution time grew exponentially with the number of spheres. A test case with 1000 spheres took about 20 times as long to solve as a 200-sphere case, and a 2000-sphere case took about 4 times longer than the 1000-sphere case. Tom's solution generated solutions that, in aggregate, occupied between 1.3 and 3.9 times the volume of individual cubes containing the individual spheres, which suggests that better solutions could be achieved with a more relaxed time penalty.

The table below lists, for each of the solutions submitted, the total volume of the boxes containing the spheres, the total execution time, and the total score including the time penalty, as well as the code and data sizes for each entry. As usual, the number in parentheses after the entrant's name is the total number of Challenge points earned in all Challenges prior to this one.

Name Volume (x1.0E12) Time (secs) Score (x1.0e12) Code Size Data Size
Tom Saxton (79)65.3142.310107.25796372
A. D.***820104

Top Contestants

Listed here are the Top Contestants for the Programmer's Challenge, including everyone who has accumulated 20 or more points during the past two years. The numbers below include points awarded over the 24 most recent contests, including points earned by this month's entrants.

  1. Munter, Ernst 200
  2. Saxton, Tom 99
  3. Boring, Randy 56
  4. Mallett, Jeff 50
  5. Rieken, Willeke 47
  6. Maurer, Sebastian 40
  7. Heithcock, JG 37
  8. Cooper, Greg 34
  9. Murphy, ACC 34
  10. Lewis, Peter 31
  11. Nicolle, Ludovic 27
  12. Brown, Pat 20
  13. Day, Mark 20
  14. Higgins, Charles 20
  15. Hostetter, Mat 20

There are three ways to earn points: (1) scoring in the top 5 of any Challenge, (2) being the first person to find a bug in a published winning solution or, (3) being the first person to suggest a Challenge that I use. The points you can win are:

1st place20 points
2nd place10 points
3rd place7 points
4th place4 points
5th place2 points
finding bug2 points
suggesting Challenge2 points

Here is Tom's winning Sphere Packing solution:

Spheres.cpp
Copyright © 1999 Tom Saxton

#include "Spheres.h"
#include "VecUtil.h"

#include <math.h>
#include <stdlib.h>

enum {
   fFalse = 0,
   fTrue = 1
};

typedef unsigned long ulong;

// disable asserts
#define Assert(f)

// hard iteration limit
#define cIterLim   10000

// scoring an accepting solutions
#define _FAccept(volNew, volBest) ((volNew) < (volBest))
#define _Score(vol, dtick)      ((vol) * (1.0 + (dtick)*10.0/60.0))

// define this to ignore the time penalty
// #define KEEP_GOING

// time checking parameters
#define dtickSec         60
#define dtickCheckScore      (dtickSec/30)
#define dtickFirstCheck      (dtickSec/30)
#define dtickLastCheck      (10*dtickSec)

static const Position s_normalX = { 1.0, 0.0, 0.0 };
static const Position s_normalY = { 0.0, 1.0, 0.0 };
static const Position s_normalZ = { 0.0, 0.0, 1.0 };
static const Position s_normalXNeg = { -1.0, 0.0, 0.0 };
static const Position s_normalYNeg = { 0.0, -1.0, 0.0 };
static const Position s_normalZNeg = { 0.0, 0.0, -1.0 };

static void _InitStartingPos(
   const long csphere,
   long aisphere[],
   const double aradius[],
   double baseMin,
   double baseBest,
   double baseMax,
   double *pbase,
   Position aposStart[]);
static void _TweakStartingPos(
   const long csphere,
   long aisphere[],
   const double aradius[],
   double baseMin,
   double baseBest,
   double baseMax,
   double *pbase,
   Position aposStart[]);
static void _DropSpheres(
   long csphere,
   const long *paisphere,
   const double aradius[],
   const Position *paposStart,
   Position apos[],
   double base,
   double *pvolume);
static void _DropOneSphere(
   const Position &posStart,
   double radius,
   int csphere,
   const long *paisphere,
   const double aradius[],
   const Position apos[],
   Position * pposResult,
   long * pisphereHit);
static int _FFindObstruction(
   const Position normalMove,
   int fNear,
   const Position &posStart,
   double radius,
   int csphere,
   const long *paisphere,
   const double aradius[],
   const Position apos[],
   Position * pposResult,
   long * pisphereHit);

PackSpheres
void PackSpheres(
  long csphere,        /* input: number of spheres to pack */
  double aradius[],    /* input: radius of each of numSpheres spheres */
  Position aposBest[]  /* output: location of center of each sphere */
)
{
   int isphere;
   double volGuess, vol, volBest;
   double base, baseMin, baseMax, baseBest;
   double radiusLarge, radiusSum;
   ulong tickStart, tickCur;
   
   tickStart = LMGetTicks();
   radiusLarge = radiusSum = 0.0;
   for (isphere = 0, volGuess = 0.0; isphere < csphere; ++isphere)
   {
      double radius = aradius[isphere];
      volGuess += 8.0 * radius * radius * radius;
      
      if (radius > radiusLarge)
         radiusLarge = radius;
      radiusSum += radius;
   }
   
   baseMin = 2.0 * radiusLarge;
   baseMax = 2.0 * radiusSum;
   Assert(baseMin <= baseMax);
   
   baseBest = baseMin;
   _DropSpheres(csphere, NULL, aradius, NULL, aposBest, 
            baseBest, &volBest);
   
   base = baseMax;
   _DropSpheres(csphere, NULL, aradius, NULL, aposBest, 
            base, &vol);
   if (vol < volBest)
   {
      volBest = vol;
      baseBest = base;
   }
   
   base = sqrt(baseMin * baseMax);
   _DropSpheres(csphere, NULL, aradius, NULL, aposBest, 
            base, &vol);
   if (vol < volBest)
   {
      volBest = vol;
      baseBest = base;
   }
   
   char * pbBlock = NewPtr(csphere * sizeof(Position) + 
                  csphere * sizeof(Position) + csphere * sizeof(long));
   
   if (pbBlock != NULL)
   {
      long iIter;

      Position * aposStart = (Position *)pbBlock;
      Position * aposEnd = &aposStart[csphere];
      long * aisphere = (long *)&aposEnd[csphere];
      long tickNext = tickStart + dtickCheckScore;
   double scorePrev = _Score(volBest, LMGetTicks() - tickStart);
#ifdef KEEP_GOING
      double scoreBest = scorePrev;
      int iIterBest = 0;
#endif   
      
      for (iIter = 0; iIter < cIterLim; ++iIter)
      {
         tickCur = LMGetTicks();
         if (tickCur >= tickNext)
         {
            ulong dtickCur = tickCur - tickStart;
            if (dtickCur >= dtickFirstCheck)
            {
               if (dtickCur >= dtickLastCheck)
                  break;
                  
               double score = _Score(volBest, dtickCur);
#ifdef KEEP_GOING
               if (score < scoreBest)
               {
                  scoreBest = score;
                  iIterBest = iIter;
               }
#else
               if (score > scorePrev)
                  break;
#endif
               scorePrev = score;
            }
            while (tickNext < tickCur)
               tickNext += dtickCheckScore;
         }
         
         // pick a new scenario
         if (iIter == 0)
            _InitStartingPos(csphere, aisphere, aradius, 
                  baseMin, baseBest, baseMax, &base, aposStart);
         else
            _TweakStartingPos(csphere, aisphere, aradius, 
               baseMin, baseBest, baseMax, &base, aposStart);
         
         // try the new scenario
         _DropSpheres(csphere, aisphere, aradius, aposStart, 
               aposEnd, base, &vol);
         if (_FAccept(vol, volBest))
         {
            volBest = vol;
            baseBest = base;
         BlockMove(aposEnd, aposBest, csphere * sizeof(Position));
         }
         
         // if the largest sphere determined the height, then reduce baseMax
      if (vol <= 2.0 * (radiusLarge + epsilon) * base * base)
         {
            Assert(base <= baseMax);
            baseMax = base;
         }
      }
   }
   
   if (pbBlock != NULL)
      DisposePtr((Ptr) pbBlock);
}

_InitStartingPos
static void _InitStartingPos(
   const long csphere,
   long aisphere[],
   const double aradius[],
   double baseMin,
   double baseBest,
   double baseMax,
   double *pbase,
   Position aposStart[])
{
   long isphereCur;
   
   *pbase = baseBest;
   for (isphereCur = 0; isphereCur < csphere; ++isphereCur)
   {
      Position *ppos = &aposStart[isphereCur];
      double radiusCur = aradius[isphereCur];
      
      aisphere[isphereCur] = isphereCur;
      ppos->coordinate[0] = 
            GRandInRange(radiusCur, *pbase - radiusCur);
      ppos->coordinate[1] = 
            GRandInRange(radiusCur, *pbase - radiusCur);
      ppos->coordinate[2] = csphere * *pbase;
   }
}

_TweakStartingPos
static void _TweakStartingPos(
   const long csphere,
   long aisphere[],
   const double aradius[],
   double baseMin,
   double baseBest,
   double baseMax,
   double *pbase,
   Position aposStart[])
{
   long isphereCur;
   double dbase;
   
   // change the base size?
   if (GRandInRange(0.0, 1.0) < 0.1)
   {
      dbase = GRandInRange(-1.0, 1.0);
      dbase *= fabs(dbase);
      dbase *= 0.25 * (baseMax - baseMin);
      *pbase = baseBest + dbase;
      *pbase = fmax(baseMin, *pbase);
      *pbase = fmin(baseMax, *pbase);
   }
   
   // rearrange the drop order?
   if (GRandInRange(0.0, 1.0) < 4.0)
   {
      for (long index = csphere; - index > 0; )
      {
         long indexSwap;
         long isphereSav;
         
         indexSwap = ((unsigned long)LRand()) % index;
         Assert(0 <= indexSwap && indexSwap < index);
         isphereSav = aisphere[index];
         aisphere[index] = aisphere[indexSwap];
         aisphere[indexSwap] = isphereSav;
      }
   }
   
   // change the starting positions
   for (isphereCur = 0; isphereCur < csphere; ++isphereCur)
   {
      Position *ppos = &aposStart[isphereCur];
      double radiusCur = aradius[isphereCur];
      
      ppos->coordinate[0] = 
            GRandInRange(radiusCur, *pbase - radiusCur);
      ppos->coordinate[1] = 
            GRandInRange(radiusCur, *pbase - radiusCur);
      ppos->coordinate[2] = csphere * *pbase;
   }
}

_DropSpheres
static void _DropSpheres(
   const long csphere,
   const long *paisphere,
   const double aradius[],
   const Position *paposStart,
   Position aposEnd[],
   double base,
   double *pvol)
{
   long csphereDone;
   
   for (csphereDone = 0; csphereDone < csphere; ++csphereDone)
   {
      Position posStart, posLand;
      double radiusCur;
      long isphereHit;
      long isphereCur;

      isphereCur = paisphere == NULL ? csphereDone : 
            paisphere[csphereDone];

      radiusCur = aradius[isphereCur];
      
      // pick a starting point for the current sphere
      Assert(base >= radiusCur*2.0);
      if (paposStart == NULL)
      {
         posStart.coordinate[0] = 
            GRandInRange(radiusCur, base - radiusCur);
         posStart.coordinate[1] = 
            GRandInRange(radiusCur, base - radiusCur);
      }
      else
      {
         posStart.coordinate[0] = 
               paposStart[isphereCur].coordinate[0];
         posStart.coordinate[1] = 
               paposStart[isphereCur].coordinate[1];
      }
      
      // drop it
   _DropOneSphere(posStart, radiusCur, csphereDone, paisphere, 
         aradius, aposEnd, &aposEnd[isphereCur], &isphereHit);
      
      // try to move it around the sphere it hit
   for (int cIter = 0; isphereHit != -1 && cIter < isphereCur; 
            ++cIter)
      {
         Position vecMove, vecMoveXY, normalMove;
         Position posHit;
         double distH, distMove;
         int icoord;
         
         posHit = aposEnd[isphereHit];
         SubVec(aposEnd[isphereCur], posHit, &vecMove);
         vecMoveXY = vecMove;
         vecMoveXY.coordinate[2] = 0;
         distH = LengthVec(vecMoveXY);
         
         if (distH < epsilon)
            break;
            
         ScaleVec(1.0/distH, vecMoveXY, &normalMove);
         distMove = radiusCur + aradius[isphereHit];
         Assert(distMove > distH - epsilon);
         
         // don't move out of the box
         for (icoord = 0; icoord <= 1; ++icoord)
         {
            if (normalMove.coordinate[icoord] < -epsilon)
            {
               if (posHit.coordinate[icoord] + 
         distMove * normalMove.coordinate[icoord] < radiusCur)
            distMove = (radiusCur - posHit.coordinate[icoord]) / 
                     normalMove.coordinate[icoord];
            }
            else if (normalMove.coordinate[icoord] > epsilon)
            {
               if (posHit.coordinate[icoord] + distMove * 
               normalMove.coordinate[icoord] > base - radiusCur)
                  distMove = (base - radiusCur - 
                                             posHit.coordinate[icoord]) / 
                                          normalMove.coordinate[icoord];
            }
         }
         
         Assert(distMove >= distH - epsilon);
         if (distMove < distH + epsilon)
            break;
            
         AddScaleVec(posHit, distMove, normalMove, &posStart);
         
   _DropOneSphere(posStart, radiusCur, csphereDone, paisphere, 
                  aradius, aposEnd, &posLand, &isphereHit);
         
         if (posLand.coordinate[2] > 
                  aposEnd[isphereCur].coordinate[2] - epsilon)
            break;
         
         aposEnd[isphereCur] = posLand;
      
      }

      // try move it toward the edges
      int fImproved, cIter;
      for (fImproved = fTrue, cIter = 1; fImproved; ++cIter)
      {
         Assert(cIter < 15);
         fImproved = fFalse;
         for (int dir = 0; dir < 4; ++dir)
         {
            Position normalMove;
            int fHit;
            double sEdge;
            Position aposStart[2];
            int cposStart;
            
            switch (dir)
            {
            case 0:
               normalMove = s_normalX;
               sEdge = base - radiusCur;
               break;
            case 1:
               normalMove = s_normalY;
               sEdge = base - radiusCur;
               break;
            case 2:
               normalMove = s_normalXNeg;
               sEdge = -radiusCur;
               break;
            case 3:
               normalMove = s_normalYNeg;
               sEdge = -radiusCur;
               break;
            }
            
            fHit = _FFindObstruction(
                     normalMove,
                     fTrue/*fNear*/,
                     aposEnd[isphereCur],
                     radiusCur,
                     csphereDone,
                     paisphere,
                     aradius,
                     aposEnd,
                     &posLand,
                     &isphereHit);
            
            cposStart = 0;
            if (!fHit || DotVec(posLand, normalMove) > sEdge)
            {
               posLand = aposEnd[isphereCur];
               AddScaleVec(posLand, sEdge - 
                        DotVec(posLand, normalMove), normalMove, 
                                    &aposStart[cposStart++]);
               cposStart = 1;
            }
            else
            {
         LinearComboVec(0.5, posLand, 0.5, aposEnd[isphereCur], 
                                 &aposStart[cposStart++]);
               aposStart[cposStart++] = posLand;
            }
            
for (int iposStart = 0; iposStart < cposStart; ++iposStart)
            {
               _DropOneSphere(aposStart[iposStart], radiusCur, 
            csphereDone, paisphere, aradius, aposEnd, &posLand, 
                     &isphereHit);
               
               if (posLand.coordinate[2] < 
                     aposEnd[isphereCur].coordinate[2] + epsilon)
               {
                  if (aposEnd[isphereCur].coordinate[2] - 
                           posLand.coordinate[2] > radiusCur * 0.05)
                     fImproved = fTrue;
                  aposEnd[isphereCur] = posLand;
            
               }
            }
         }
      }
   }
   
   ComputeVol(csphere, NULL, aradius, aposEnd, base, pvol);
}

_DropOneSphere
static void _DropOneSphere(
   const Position &posStart,
   double radius,
   int csphere,
   const long *paisphere,
   const double aradius[],
   const Position apos[],
   Position * pposResult,
   long * pisphereHit)
{
   Position posBase;
   int fHit;
   
   posBase = posStart;
   posBase.coordinate[2] = 0.0;
   
   *pposResult = posBase;

   fHit = _FFindObstruction(
            s_normalZ,
            fFalse, /* fNear */
            posBase,
            radius,
            csphere,
            paisphere,
            aradius,
            apos,
            pposResult,
            pisphereHit
            );
   
   if (!fHit || pposResult->coordinate[2] < radius)
   {
      *pisphereHit = -1;
      pposResult->coordinate[2] = radius;
   }

   // add some fudge
   pposResult->coordinate[2] += epsilon;
   
#ifdef DEBUG
   for (long csphereChecked = 0; csphereChecked < csphere; 
               ++csphereChecked)
   {
      Position vecT;
      double dist, distGap;
      int isphere;
      
      isphere = paisphere == NULL ? csphereChecked : 
                                                paisphere[csphereChecked];
      
      SubVec(apos[isphere], *pposResult, &vecT);
      dist = LengthVec(vecT);
      distGap = dist - (radius + aradius[isphere]);
      Assert(distGap >= 0.0);
   }
#endif
}

_FFindObstruction
// moving a sphere with specifed radius from posStart in the direction normalMove,
// find the nearest or farthest obstruction
// If there is an obstruction, return the index to the obstructing sphere
// and the position to which the object can move.
static int _FFindObstruction(
   const Position normalMove,
   int fNear,
   const Position &posStart,
   double radius,
   int csphere,
   const long *paisphere,
   const double aradius[],
   const Position apos[],
   Position * pposResult,
   long * pisphereHit)
{
   double zBest;
   
   *pisphereHit = -1;

   for (int csphereChecked = 0; csphereChecked < csphere; 
               ++csphereChecked)
   {
      Position vecToOther, vecPerp, vecParallel;
      double distPerpSq, distSep, distSepSq;
      double z, dz;
      int isphere;
      
      isphere = paisphere == NULL ? csphereChecked : 
                                                paisphere[csphereChecked];
      SubVec(apos[isphere], posStart, &vecToOther);
      ProjectVec(vecToOther, normalMove, &vecParallel);
      SubVec(vecToOther, vecParallel, &vecPerp);
      
      distPerpSq = DotVec(vecPerp, vecPerp);
      distSep = radius + aradius[isphere];
      distSepSq = distSep * distSep;
      
      if (distPerpSq > distSepSq)
         continue;
      
      dz = sqrt(distSepSq - distPerpSq);
      if (fNear)
         dz = -dz;
      z = DotVec(vecParallel, normalMove) + dz;
      
      if (z >= 0.0 && (*pisphereHit == -1 || 
                                    (fNear ? z < zBest : z > zBest)))
      {
         zBest = z;
         *pisphereHit = isphere;
      }
   }
   
   if (*pisphereHit == -1)
      return fFalse;
      
   *pposResult = posStart;
   AddScaleVec(posStart, zBest, normalMove, pposResult);
   
   return fTrue;
}

VecUtil.cpp
#include "Spheres.h"
#include "VecUtil.h"

#include <math.h>
#include <stdlib.h>

enum {
   fFalse = 0,
   fTrue = 1
};

// disable asserts
#define Assert(f)

// math utilities

double GRandInRange(double gLow, double gHigh)
{
   double g;
   
   g = gLow + rand() * (gHigh - gLow) / RAND_MAX;
   Assert(gLow <= g && g <= gHigh);
   return g;
}

// return a long's worth of randomness
long LRand()
{
   long lw;
   
   Assert(RAND_MAX > 256);
   
   lw = 0;
   for (int ib = 0; ib < sizeof(long); ++ib)
      lw = (lw << 8) + (rand() & 0xFF);
   return lw;
}

// vector utilities

void SubVec(const Position &pos1, const Position &pos2, 
      Position *pposResult)
{
   for (int i = 0; i < 3; ++i)
      pposResult->coordinate[i] = pos1.coordinate[i] - 
                                                         pos2.coordinate[i];
}

double DotVec(const Position &pos1, const Position &pos2)
{
   double g = 0.0;
   for (int i = 0; i < 3; ++i)
      g += pos1.coordinate[i] * pos2.coordinate[i];
   return g;
}

double LengthVec(const Position &pos)
{
   return sqrt(DotVec(pos, pos));
}

void ScaleVec(double g, const Position &pos, 
   Position *pposResult)
{
   for (int i = 0; i < 3; ++i)
      pposResult->coordinate[i] = g * pos.coordinate[i];
}

void AddScaleVec(const Position &posBase, double g, 
   const Position &posAdd, Position *pposResult)
{
   for (int i = 0; i < 3; ++i)
      pposResult->coordinate[i] = posBase.coordinate[i] + 
                                                g * posAdd.coordinate[i];
}

void LinearComboVec(double g1, const Position &pos1, double g2, 
      const Position &pos2, Position *pposResult)
{
   for (int i = 0; i < 3; ++i)
      pposResult->coordinate[i] = g1 * pos1.coordinate[i] + 
                                                g2 * pos2.coordinate[i];
}

// project "vec" onto a "normal" vector
void ProjectVec(const Position &vec, const Position &normal, 
      Position *pvecResult)
{
   ScaleVec(DotVec(vec, normal), normal, pvecResult);
}

// sphere stuff

void ComputeVol(
   const long csphere,
   const long *paisphere,
   const double aradius[],
   const Position apos[],
   double base,
   double *pvol)
{
   Position posMin, posMax;
   long index;
   int icoord;
   double radius;
   const Position * ppos;

   posMin = posMax = apos[0];

   for (index = 0; index < csphere; ++index)
   {
      long isphere;
      
      isphere = paisphere == NULL ? index : paisphere[index];
      ppos = &apos[isphere];
      radius = aradius[isphere];
      for (icoord = 0; icoord < 3; ++icoord)
      {
         if (ppos->coordinate[icoord] - radius < 
                  posMin.coordinate[icoord])
            posMin.coordinate[icoord] = ppos->coordinate[icoord] - 
                                                            radius;

         if (ppos->coordinate[icoord] + radius > 
                  posMax.coordinate[icoord])
            posMax.coordinate[icoord] = ppos->coordinate[icoord] + 
                                                            radius;
      }
   }

   *pvol = 1.0;
   
   for (icoord = 0; icoord < 3; ++icoord)
   {
      Assert(posMin.coordinate[icoord] >= -epsilon);
   Assert(base == 0 || posMax.coordinate[icoord] <= base+epsilon 
                                       || icoord == 2);
      *pvol *= posMax.coordinate[icoord] - 
                        posMin.coordinate[icoord];
   }
}

Spheres.h

#if defined(__cplusplus)
extern "C" {
#endif

typedef struct Position {
  double coordinate[3];  /* coordinate[0]==X position, [1]==Y, [2]==Z */
} Position;

void PackSpheres(
  long numSpheres,       /* input: number of spheres to pack */
  double radius[],       /* input: radius of each of numSpheres spheres */
  Position location[]    /* output: location of center of each sphere */
);

#if defined (__cplusplus)
}
#endif

VecUtil.h
// error tolerance

const double epsilon (1.0e-10);

// math utilities

double GRandInRange(double gLow, double gHigh);
long LRand();

// vector utilities

void SubVec(const Position &pos1, const Position &pos2, 
   Position *pposResult);
double DotVec(const Position &pos1, const Position &pos2);
double LengthVec(const Position &pos);
void ScaleVec(double g, const Position &pos, 
   Position *pposResult);
void AddScaleVec(const Position &posBase, double g, 
   const Position &posAdd, Position *pposResult);
void ProjectVec(const Position &vec, const Position &normal, 
   Position *pvecResult);
void LinearComboVec(double g1, const Position &pos1, 
   double g2, const Position &pos2, Position *pposResult);

// sphere stuff

void ComputeVol(
   const long csphere,
   const long *paisphere,
   const double aradius[],
   const Position apos[],
   double base,
   double *pvol);
 
AAPL
$111.78
Apple Inc.
-0.87
MSFT
$47.66
Microsoft Corpora
+0.14
GOOG
$516.35
Google Inc.
+5.25

MacTech Search:
Community Search:

Software Updates via MacUpdate

LibreOffice 4.3.5.2 - Free Open Source o...
LibreOffice is an office suite (word processor, spreadsheet, presentations, drawing tool) compatible with other major office suites. The Document Foundation is coordinating development and... Read more
CleanApp 5.0.0 Beta 5 - Application dein...
CleanApp is an application deinstaller and archiver.... Your hard drive gets fuller day by day, but do you know why? CleanApp 5 provides you with insights how to reclaim disk space. There are... Read more
Monolingual 1.6.2 - Remove unwanted OS X...
Monolingual is a program for removing unnecesary language resources from OS X, in order to reclaim several hundred megabytes of disk space. It requires a 64-bit capable Intel-based Mac and at least... Read more
NetShade 6.1 - Browse privately using an...
NetShade is an Internet security tool that conceals your IP address on the web. NetShade routes your Web connection through either a public anonymous proxy server, or one of NetShade's own dedicated... Read more
calibre 2.13 - Complete e-library manage...
Calibre is a complete e-book library manager. Organize your collection, convert your books to multiple formats, and sync with all of your devices. Let Calibre be your multi-tasking digital librarian... Read more
Mellel 3.3.7 - Powerful word processor w...
Mellel is the leading word processor for OS X and has been widely considered the industry standard since its inception. Mellel focuses on writers and scholars for technical writing and multilingual... Read more
ScreenFlow 5.0.1 - Create screen recordi...
Save 10% with the exclusive MacUpdate coupon code: AFMacUpdate10 Buy now! ScreenFlow is powerful, easy-to-use screencasting software for the Mac. With ScreenFlow you can record the contents of your... Read more
Simon 4.0 - Monitor changes and crashes...
Simon monitors websites and alerts you of crashes and changes. Select pages to monitor, choose your alert options, and customize your settings. Simon does the rest. Keep a watchful eye on your... Read more
BBEdit 11.0.2 - Powerful text and HTML e...
BBEdit is the leading professional HTML and text editor for the Mac. Specifically crafted in response to the needs of Web authors and software developers, this award-winning product provides a... Read more
ExpanDrive 4.2.1 - Access cloud storage...
ExpanDrive builds cloud storage in every application, acts just like a USB drive plugged into your Mac. With ExpanDrive, you can securely access any remote file server directly from the Finder or... Read more

Latest Forum Discussions

See All

Make your own Tribez Figures (and More)...
Make your own Tribez Figures (and More) with Toyze Posted by Jessica Fisher on December 19th, 2014 [ permalink ] Universal App - Designed for iPhone and iPad | Read more »
So Many Holiday iOS Sales Oh My Goodness...
The holiday season is in full-swing, which means a whole lot of iOS apps and games are going on sale. A bunch already have, in fact. Naturally this means we’re putting together a hand-picked list of the best discounts and sales we can find in order... | Read more »
It’s Bird vs. Bird in the New PvP Mode f...
It’s Bird vs. Bird in the New PvP Mode for Angry Birds Epic Posted by Jessica Fisher on December 19th, 2014 [ permalink ] Universal App - Designed for iPhone and iPad | Read more »
Telltale Games and Mojang Announce Minec...
Telltale Games and Mojang Announce Minecraft: Story Mode – A Telltale Games Series Posted by Jessica Fisher on December 19th, 2014 [ permalink ] | Read more »
WarChest and Splash Damage Annouce Their...
WarChest and Splash Damage Annouce Their New Game: Tempo Posted by Jessica Fisher on December 19th, 2014 [ permalink ] WarChest Ltd and Splash Damage Ltd are teaming up again to work | Read more »
BulkyPix Celebrates its 6th Anniversary...
BulkyPix Celebrates its 6th Anniversary with a Bunch of Free Games Posted by Jessica Fisher on December 19th, 2014 [ permalink ] BulkyPix has | Read more »
Indulge in Japanese cuisine in Cooking F...
Indulge in Japanese cuisine in Cooking Fever’s new sushi-themed update Posted by Simon Reed on December 19th, 2014 [ permalink ] Lithuanian developer Nordcurrent has yet again updated its restaurant simulat | Read more »
Badland Daydream Level Pack Arrives to C...
Badland Daydream Level Pack Arrives to Celebrate 20 Million Downloads Posted by Ellis Spice on December 19th, 2014 [ permalink ] | Read more »
Far Cry 4, Assassin’s Creed Unity, Desti...
Far Cry 4, Assassin’s Creed Unity, Destiny, and Beyond – AppSpy Takes a Look at AAA Companion Apps Posted by Rob Rich on December 19th, 2014 [ permalink ] These day | Read more »
A Bunch of Halfbrick Games Are Going Fre...
A Bunch of Halfbrick Games Are Going Free for the Holidays Posted by Ellis Spice on December 19th, 2014 [ permalink ] Universal App - Designed for iPhone and iPad | Read more »

Price Scanner via MacPrices.net

Kodak Returns to CES With New Consumer Produ...
Former photography colossus Kodak is returning to CES for the first time in three years where the Kodak booth (#21818 South Hall 1) will showcase a wide range of innovative, imaging-related products... Read more
Invaluable Launches New Eponymously -Named A...
Invaluable, the world’s largest online live auction marketplace, hhas announced the official launch of the Invaluable app for iPad, now available for download in the iTunes App Store. Invaluable... Read more
IDC Reveals Worldwide Mobile Enterprise Appli...
International Data Corporation (IDC) last week hosted the IDC FutureScape: Worldwide Mobile Enterprise Applications and Solutions 2015 Predictions Web conference. The session provided organizations... Read more
Hello Vino Wine App Launches “Safe Ride Home”...
Hello Vino has announced addition of a new “Get a Safe Ride Home” feature in its Food & Drink app with a direct connection to Uber, the technology platform that connects users with rides. The... Read more
DEVON-technologies Releases DEVONthink To Go...
Coeur d’Alene, Idaho based DEVON-technologies, LLC has updated DEVONthink To Go, its mobile companion to DEVONthink, to version 1.5. The update includes an iOS 8 extension, compatibility with the... Read more
The Apple Store offering free next-day shippi...
The Apple Store is now offering free next-day shipping on all in stock items if ordered before 12/23/14 at 10:00am PT. Local store pickup is also available within an hour of ordering for any in stock... Read more
It’s 1992 Again At Sony Pictures, Except For...
Techcrunch’s John Biggs interviewed a Sony Pictures Entertainment (SPE) employee, who quite understandably wished to remain anonymous, regarding post-hack conditions in SPE’s L.A office, explaining “... Read more
Holiday sales this weekend: MacBook Pros for...
 B&H Photo has new MacBook Pros on sale for up to $300 off MSRP as part of their Holiday pricing. Shipping is free, and B&H charges NY sales tax only: - 15″ 2.2GHz Retina MacBook Pro: $1699... Read more
Holiday sales this weekend: MacBook Airs for...
B&H Photo has 2014 MacBook Airs on sale for up to $120 off MSRP, for a limited time, for the Thanksgiving/Christmas Holiday shopping season. Shipping is free, and B&H charges NY sales tax... Read more
Holiday sales this weekend: iMacs for up to $...
B&H Photo has 21″ and 27″ iMacs on sale for up to $200 off MSRP including free shipping plus NY sales tax only. B&H will also include a free copy of Parallels Desktop software: - 21″ 1.4GHz... Read more

Jobs Board

*Apple* Store Leader Program (US) - Apple, I...
…Summary Learn and grow as you explore the art of leadership at the Apple Store. You'll master our retail business inside and out through training, hands-on experience, Read more
Project Manager, *Apple* Financial Services...
**Job Summary** Apple Financial Services (AFS) offers consumers, businesses and educational institutions ways to finance Apple purchases. We work with national and Read more
*Apple* Retail - Multiple Positions (US) - A...
Sales Specialist - Retail Customer Service and Sales Transform Apple Store visitors into loyal Apple customers. When customers enter the store, you're also the Read more
*Apple* Retail - Multiple Positions (US) - A...
Sales Specialist - Retail Customer Service and Sales Transform Apple Store visitors into loyal Apple customers. When customers enter the store, you're also the Read more
*Apple* Retail - Multiple Positions (US) - A...
Job Description: Sales Specialist - Retail Customer Service and Sales Transform Apple Store visitors into loyal Apple customers. When customers enter the store, Read more
All contents are Copyright 1984-2011 by Xplain Corporation. All rights reserved. Theme designed by Icreon.