Perl Advent Calendar 2008-12-03

Sleighing big numbers

by Bill Ricker

How Santa's IT Dept Saves Xmas (installment 2)

ATTN Comics Editors -- embargo Dec 3

Santa's Route planning department and his payload planning department have to work together. He does not have time to search through multiple sacks at each chimney for the right items, the sleigh needs to be loaded with the precision of a FedEx or UPS truck - the next present must be next to hand. So load and route sequence are inextricably linked.

We suspect the sleigh and reindeer must have elven teleport abilities (whether wormhole or other transdimensional passage or some anti-relativistic quirk), since the Classical acceleration required to make all the required stops in a mere 30 to 50 hours (from International Date Line (IDL) date start to IDL date stop) is horrendous, and the sonic booms would have been reported. Besides, his time on station at each stop improves dramatically if both up-the-chimney and roof to roof are instantaneous.

Let's check the calculations on Scientific Santa's explanation. Absolute precision is not required for optimization, they need plan and replan … fast; and "the best" is the enemy of good enough. So the Elven planners and their Neptunian outsourced IT staff use Math::BigApprox as an electronic sliderule for back-of-the envelope estimates when they fear blowout of floating point in their scripts (and Frink for interactive use).

Following Scientific Santa, we start with the UNICEF world child census and deduce the physics equation by equation.

2.106000e+09 children per unicef
8.424000e+08 stops (children/2.5)
3.141593e+00 pi
3.986000e+03 Re Earth radius, miles. 
1.996570e+08 area sq. miles. 
5.790052e+07 populated area (dry 29%) sq. miles. 
6.873281e-02 area each, sq.miles / houses
2.621694e-01 mean distance, miles (sqrt area each)
2.208515e+08 mileage= houses * mean distance; miles
1.728000e+05 seconds total flight
2.051282e-04 Time between each, s
1.278076e+03 speed avg, miles/second
4.601074e+06 speed avg, miles/hr
4.601074e+06 speed avg, m/s (metric)
6.134765e+03 Mach - web page off x1000 !
3.000000e+08 C, speed of light, metric, m/s
1.458531e+02 ratio C:spd :: _:1
4.010885e+10 acceleration m/s/s
4.092740e+09 accel G's
1.914545e+12 Cargo grams
4.212000e+09 Cargo Lbs (2 Lbs/Child)
7.679021e+19 F=ma: Force, Newtons
1.7e+7153006998 traveling salesman, houses!

We find the Scientific Santa slipped a decimal on Mach number, but otherwise holds up with assumptions given.

(Of course, out-sourcing delivery to various cultural gift bringers - Befanna, Three Kings/The Magi, Grandfather Winter, Pere Noel, St Basil, Christ Child, Ste Lucia, Father Christmas, Old Man Frost, Tio Nadal (Uncle Xmas), etc, greatly aids Santa's actual workload. And their team being able to do some households on St.Stephens Day or Three Kings Day or Solstice Eve or Lucia's Day also helps them all.)

Alas, even calculating the size of their un-structured Traveling Salesman problem as Fact($houses) takes 13 hours (47,212 seconds).


fact 1e1 = 3628800 [0s]
fact 1e2 = 9.33262154e+157 [0s]
fact 1e3 = 4.0238726 e+2567 [0s]
fact 1e4 = 2.84626   e+35659 [0s]
fact 1e5 = 2.82423   e+456573 [1s]
fact 1e6 = 8.2639    e+5565708 [6s]
fact 1e7 = 1.202     e+65657059 [70s]
fact 1e8 = 1.62      e+756570556 [646s]
2.106000e+09     children
8.424000e+08     houses
traveling salesman
fact 842400000 = 1.6e+7153006998 [47212s]
=end pre

This is hyperexponential, I would have expected 2 hours.

So they have to replace Math::BigApprox's Fact() with Stirling's Approximation when N is large enough that M:BA is barely keeping precision beyond order of magnitude anyway.


fact 1e1 = 3598695.6187 [0s]
fact 1e2 = 9.32484763e+157 [0s]
fact 1e3 = 4.0235373 e+2567 [0s]
fact 1e4 = 2.846236  e+35659 [0s]
fact 1e5 = 2.82423   e+456573 [0s]
fact 1e6 = 8.2639    e+5565708 [0s]
fact 1e7 = 1.202     e+65657059 [0s]
fact 1e8 = 1.62      e+756570556 [0s]
fact 1e9  = 1        e+8565705523 [0s]
fact 1e10 = 2        e+95657055186 [0s]
fact 1e11 = 1        e+1056570551820 [0s]
fact 1e12 = 1359     e+11565705518100 [0s]
fact 1e13 = 1        e+125657055181000 [0s]
fact 1e14 = 1        e +1.35657055181 e+15 [0s]
2.106000e+09     children
8.424000e+08     houses
traveling salesman
fact 842400000 = 1.7e+7153006998 [0s]
=end pre

This performs in constant time, which is pretty good.

1.7e+7_153_006_998. Ten to the 7 Billion. That's large. That many routes to optimize between is obviously not going to be precisely solved! Can the Eleven Planners and Neptunian IT find a fast enough good enough approximation in time to save Christmas?

(To Be Continued)

mod3.pl

   1 #! perl -lw
   2 
   3 use warnings;
   4 use strict;
   5 
   6 # require Math::BigApprox;
   7 use Math::BigApprox qw( c Prod );
   8 
   9 # Pi -
  10 # if it was BigFloat, I'd use
  11 # 3.14159_26535_89793_23846_26433_83279_50288_41971_69399;
  12 # but with Math::BigApprox
  13 my $pi = c( 355 / 113 ); # is good enough, better than 22/7, better than 3.14159
  14 
  15 # oddly,
  16 #  $pi = c(sqrt sqrt (9**2 + 19**2/22)); # is also as good.  http://xkcd.com/217/
  17 # http://use.perl.org/~n1vux/journal/32292
  18 # http://use.perl.org/~n1vux/journal/32686
  19 
  20 # Factorial -
  21 # Override Factorial if Stirling approximation good enough
  22 
  23 sub Stirling {
  24 
  25     # Stirlings approximate formula for Factorial
  26     # http://en.wikipedia.org/wiki/Stirling%27s_formula
  27     my $x = shift;
  28     my $n = ref $x eq "Math::BigApprox" ? $x : c($x);
  29     return sqrt( 2 * $pi * $n ) * ( $n / exp(1) )**$n;
  30 }
  31 
  32 sub Fact { return $_[0] > 1e4 ? Stirling shift : Math::BigApprox::Fact shift; }
  33 
  34 # -----------------------
  35 
  36 my $children = c( 2.106 * 10**9 );    # billion children aged under 18  UNICEF
  37 
  38 #when? update?
  39 printf "%e %s\n", $children, q{children per unicef};
  40 
  41 my $houses = $children / 2.5;         #Assume there are 2.5 children a house
  42 
  43 printf "%e stops (children/2.5)\n",
  44   $houses;                            # 842 million stops on Christmas Eve
  45 
  46 printf "%e %s\n", $pi, q{pi};
  47 
  48 my $radius = c(3_986);                # earth, miles
  49 printf "%e %s\n", $radius, q{Re Earth radius, miles. };
  50 
  51 my $area = ( $radius**2 ) * $pi * 4;
  52 printf "%e %s\n", $area, q{area sq. miles. };
  53 
  54 $area = $area * 0.29;                 # % dry
  55 
  56 printf "%e %s\n", $area, q{populated area (dry 29%) sq. miles. };
  57 $area = $area / $houses;
  58 printf "%e %s\n", $area, q{area each, sq.miles / houses};
  59 
  60 my $distPer = sqrt $area;
  61 
  62 printf "%e %s\n", $distPer, q{mean distance, miles (sqrt area each)};
  63 
  64 my $dist = $distPer * $houses;
  65 
  66 printf "%e %s\n", $dist, q{mileage= houses * mean distance; miles};
  67 
  68 my $time = c(48)    # hours IDL to IDL
  69   * 60 * 60;        # sec
  70 printf "%e %s\n", $time, q{seconds total flight};    #
  71 
  72 my $timePer = $time / $houses;                       #s
  73 printf "%e %s\n", $timePer, q{Time between each, s}; #
  74 
  75 my $speed = $dist / $time;                           #mi p s
  76 
  77 printf "%e %s\n", $speed, q{speed avg, miles/second};    #
  78 
  79 my $milePerMetre = 1 / (
  80     5_280                                                # ft/mi
  81       * 12                                               # in/ft
  82       * 2.54                                             # cm/in
  83       * 1 / 100                                          # m/cm
  84   )                                                      # 1/(m/mi) = mi/m
  85   ;
  86 
  87 my $speedMetric = $speed / $milePerMetre;                # mps
  88 
  89 $speed = $speed * 3600;                                  #MPH
  90 
  91 printf "%e %s\n", $speed, q{speed avg, miles/hr};
  92 printf "%e %s\n", $speed, q{speed avg, m/s (metric)};
  93 
  94 printf "%e %s\n", $speed / 750, q{Mach - web page off x1000 !};
  95 
  96 # web page slipped 3 decimals here !!
  97 
  98 my $C = 3E8;    # speed of light, metric, mps
  99 printf "%e %s\n", $C, q{C, speed of light, metric, m/s};
 100 printf "%e %s\n", 1 / ( $speedMetric / $C ),
 101   q{ratio C:spd :: _:1};    # speed/C, in mps
 102 
 103 #########
 104 # acceleration
 105 # Article is off by at least a factor of two here, must accell and then decell.
 106 # to average S, must accelerate to twice that and brake just as hard.
 107 
 108 # Real math is
 109 # R = 0.5 a*t**2. in this case,
 110 # R = 2 * (0.5 a* (dt/2)**2)
 111 # $dist = $a (timePer/2)**2
 112 
 113 my $accell = $distPer /     # miles
 114   (
 115     $milePerMetre * ( ( $timePer / 2 )**2 )    # s.s
 116   );
 117 printf "%e %s\n", $accell, q{acceleration m/s/s};
 118 
 119 my $g = c(9.8);                                # g=9.8m/s/s at surface
 120 printf "%e %s\n", $accell / $g, q{accel G's};
 121 
 122 ##########
 123 #  Cargo
 124 
 125 my $Cargo = ( 2 / 2.2 )                        # kg
 126   * $children;
 127 printf "%e %s\n", $Cargo * 1000, q{Cargo grams};
 128 printf "%e %s\n", $Cargo * 2.2,  q{Cargo Lbs (2 Lbs/Child)};
 129 
 130 printf "%e %s\n", $Cargo * $accell, q{F=ma: Force, Newtons};
 131 
 132 # work
 133 # food callories / 12 reindeer
 134 
 135 print Fact($houses), q{ traveling salesman, houses!};
 136 
 137 # printf "%e %s\n", c(1) << 31, q{2^31 for comparison};
 138 # 2.147484e+09 2^31 for comparison
 139 # use POSIX qw(FLT_MAX DBL_MAX);
 140 # printf "%e %s\n", DBL_MAX, q{DBL_MAX for comparison};
 141 # 1.797693e+308 DBL_MAX for comparison
 142 # printf "%e %s\n", FLT_MAX, q{FLT_MAX for comparison};
 143 # 3.402823e+38 FLT_MAX for comparison
 144 

Math::BigApprox certainly has its uses. While for most purposes 1.797693e+308 DBL_MAX for 32 bit perl is quite adequate, there are times when the extended range of Math::BigWhatever is needed but not the costly precision. With the Stirling approximation, it did let me compute factorial 10^30 which is sometimes necessary.

POD