src/wheel-gen.pl

   1 #!/usr/bin/perl -w
   2 # Generate the spokes of a wheel, for wheel factorization.
   3
   4 # Copyright (C) 2001, 2005, 2009-2010 Free Software Foundation, Inc.
   5
   6 # This program is free software: you can redistribute it and/or modify
   7 # it under the terms of the GNU General Public License as published by
   8 # the Free Software Foundation, either version 3 of the License, or
   9 # (at your option) any later version.
  10
  11 # This program is distributed in the hope that it will be useful,
  12 # but WITHOUT ANY WARRANTY; without even the implied warranty of
  13 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  14 # GNU General Public License for more details.
  15
  16 # You should have received a copy of the GNU General Public License
  17 # along with this program.  If not, see <http://www.gnu.org/licenses/>.
  18
  19 eval 'exec /usr/bin/perl -S $0 ${1+"$@"}'
  20   if 0;
  21
  22 use strict;
  23 (my $ME = $0) =~ s|.*/||;
  24
  25 # A global destructor to close standard output with error checking.
  26 sub END
  27 {
  28   defined fileno STDOUT
  29     or return;
  30   close STDOUT
  31     and return;
  32   warn "$ME: closing standard output: $!\n";
  33   $? ||= 1;
  34 }
  35
  36 sub is_prime ($)
  37 {
  38   my ($n) = @_;
  39   use integer;
  40
  41   $n == 2
  42     and return 1;
  43
  44   my $d = 2;
  45   my $w = 1;
  46   while (1)
  47     {
  48       my $q = $n / $d;
  49       $n == $q * $d
  50         and return 0;
  51       $d += $w;
  52       $q < $d
  53         and last;
  54       $w = 2;
  55     }
  56   return 1;
  57 }
  58
  59 {
  60   @ARGV == 1
  61     or die "$ME: missing argument\n";
  62
  63   my $wheel_size = $ARGV[0];
  64
  65   my @primes = (2);
  66   my $product = $primes[0];
  67   my $n_primes = 1;
  68   for (my $i = 3; ; $i += 2)
  69     {
  70       if (is_prime $i)
  71         {
  72           push @primes, $i;
  73           $product *= $i;
  74           ++$n_primes == $wheel_size
  75             and last;
  76         }
  77     }
  78
  79   my $ws_m1 = $wheel_size - 1;
  80   print <<EOF;
  81 /* The first $ws_m1 elements correspond to the incremental offsets of the
  82    first $wheel_size primes (@primes).  The $wheel_size(th) element is the
  83    difference between that last prime and the next largest integer
  84    that is not a multiple of those primes.  The remaining numbers
  85    define the wheel.  For more information, see
  86    http://www.utm.edu/research/primes/glossary/WheelFactorization.html.  */
  87 EOF
  88
  89   my @increments;
  90   my $prev = 2;
  91   for (my $i = 3; ; $i += 2)
  92     {
  93       my $rel_prime = 1;
  94       foreach my $divisor (@primes)
  95         {
  96           $i != $divisor && $i % $divisor == 0
  97             and $rel_prime = 0;
  98         }
  99
 100       if ($rel_prime)
 101         {
 102           #warn $i, ' ', $i - $prev, "\n";
 103           push @increments, $i - $prev;
 104           $prev = $i;
 105
 106           $product + 1 < $i
 107             and last;
 108         }
 109     }
 110
 111   print join (",\n", @increments), "\n";
 112
 113   exit 0;
 114 }