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);
 

Community Search:
MacTech Search:

Software Updates via MacUpdate

How to manage your time in Bakery Blitz
It can be tricky, especially when you risk burning your kitchen to the ground if you forget a cake in the oven, so make sure to use these time management tricks to keep your bakery running smoothly. Don’t collect the money right away [Read more] | Read more »
Model 15 (Music)
Model 15 1.0 Device: iOS iPhone Category: Music Price: $29.99, Version: 1.0 (iTunes) Description: The Moog Model 15 App is the first Moog modular synthesizer and synthesis educational tool created exclusively for iPad, iPhone and... | Read more »
How to deal with wind in Angry Birds Act...
Angry Birds Action! is a physics-based puzzler in which you're tasked with dragging and launching birds around an obstacle-littered field to achieve a set objective. It's simple enough at first, but when wind gets introduced things can get pretty... | Read more »
How to get three stars in every level of...
Angry Birds Action! is, essentially, a pinball-style take on the pull-and-fling action of the original games. When you first boot it up, you'll likely be wondering exactly what it is you have to do to get a good score. Well, never fear as 148Apps... | Read more »
The beginner's guide to Warbits
Warbits is a turn-based strategy that's clearly inspired by Nintendo's Advance Wars series. Since turn-based strategy games can be kind of tricky to dive into, see below for a few tips to help you in the beginning. Positioning is crucial [Read... | Read more »
How to upgrade your character in Spellsp...
So you’ve mastered the basics of Spellspire. By which I mean you’ve realised it’s all about spelling things in a spire. What next? Well you’re going to need to figure out how to toughen up your character. It’s all well and good being able to spell... | Read more »
5 slither.io mash-ups we'd love to...
If there's one thing that slither.io has proved, it's that the addictive gameplay of Agar.io can be transplanted onto basically anything and it will still be good fun. It wouldn't be surprising if we saw other developers jumping on the bandwagon,... | Read more »
How to navigate the terrain in Sky Charm...
Sky Charms is a whimsical match-'em up adventure that uses creative level design to really ramp up the difficulty. [Read more] | Read more »
Victorious Knight (Games)
Victorious Knight 1.3 Device: iOS Universal Category: Games Price: $1.99, Version: 1.3 (iTunes) Description: New challenges awaits you! Experience fresh RPG experience with a unique combat mechanic, packed with high quality 3D... | Read more »
Agent Gumball - Roguelike Spy Game (Gam...
Agent Gumball - Roguelike Spy Game 1.0 Device: iOS Universal Category: Games Price: $2.99, Version: 1.0 (iTunes) Description: Someone’s been spying on Gumball. What the what?! Two can play at that game! GO UNDERCOVERSneak past enemy... | Read more »

Price Scanner via MacPrices.net

13-inch 2.5GHz MacBook Pro on sale for $999,...
B&H Photo has the 13″ 2.5GHz MacBook Pro on sale for $999 including free shipping plus NY sales tax only. Their price is $100 off MSRP. Read more
Apple refurbished 2015 iMacs available for up...
Apple now has a full line of Certified Refurbished 2015 21″ & 27″ iMacs available for up to $350 off MSRP. Apple’s one-year warranty is standard, and shipping is free. The following models are... Read more
Indian Smartphone Market Grows Annually by 12...
India’s smartphone market grew by 12 percent year-over-year, with 24.4 million units shipping in Q1 2016. The top five vendors stayed the same, with Samsung in the lead, followed by Micromax, Intex... Read more
Get Notifications When Your Friend’s Phone Ba...
Calgary, Canada based Stonelight Pictures has announced the release of Battery Share 1.0.1, its new utility for iOS 9 supported devices. The company notes that people are spending more time on their... Read more
11-inch 1.6GHz/128GB MacBook Air on sale for...
Amazon has the current-generation 11″ 1.6GHz/128GB MacBook Air (sku MJVM2LL/A) on sale for $749.99 for a limited time. Their price is $150 off MSRP, and it’s the lowest price available for this model... Read more
Price drops on clearance 2015 13-inch MacBook...
B&H Photo has dropped prices on clearance 2015 13″ MacBook Airs by up to $250. Shipping is free, and B&H charges NY sales tax only: - 13″ 1.6GHz/4GB/128GB MacBook Air (MJVE2LL/A): $799, $200... Read more
Mac minis on sale for up to $100 off MSRP
B&H Photo has Mac minis on sale for up to $100 off MSRP including free shipping plus NY sales tax only: - 1.4GHz Mac mini: $449 $50 off MSRP - 2.6GHz Mac mini: $649 $50 off MSRP - 2.8GHz Mac mini... Read more
13-inch Retina MacBook Pros on sale for up to...
B&H Photo has 13″ Retina MacBook Pros on sale for $130-$200 off MSRP. Shipping is free, and B&H charges NY tax only: - 13″ 2.7GHz/128GB Retina MacBook Pro: $1169 $130 off MSRP - 13″ 2.7GHz/... Read more
Apple price trackers, updated continuously
Scan our Apple Price Trackers for the latest information on sales, bundles, and availability on systems from Apple’s authorized internet/catalog resellers. We update the trackers continuously: - 15″... Read more
SanDisk Half-Terabyte SSD Optimized for Every...
SanDisk Corporation has announced the SanDisk Z410 SSD, a cost-competitive, half-terabyte solid state drive (SSD) that enables manufacturers to design for a broad range of desktop PCs and laptops.... Read more

Jobs Board

Restaurant Manager (Neighborhood Captain) - A...
…in every aspect of daily operation WHY YOU LL LIKE IT You ll be the Big Apple You ll solve problems You ll get to show your ability to handle the stress and Read more
*Apple* Retail - Multiple Positions (US) - A...
Job Description:SalesSpecialist - Retail Customer Service and SalesTransform Apple Store visitors into loyal Apple customers. When customers enter the store, Read more
Restaurant Manager (Neighborhood Captain) - A...
…in every aspect of daily operation. WHY YOU'LL LIKE IT: You'll be the Big Apple . You'll solve problems. You'll get to show your ability to handle the stress and Read more
*Apple* Subject Matter Expert - NTT Data, In...
…in Owings Mills, MD has a 6+ month contract position available for an Apple Subject Matter Expert. TITLE: Apple Subject Matter Expert LOCATION: Owings Mills, Read more
*Apple* Retail - Multiple Positions - Apple,...
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.