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 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 2, or (at your option)
   9 # 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, write to the Free Software Foundation,
  18 # Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
  19
  20 eval 'exec /usr/bin/perl -S $0 ${1+"$@"}'
  21   if 0;
  22
  23 use strict;
  24 (my $ME = $0) =~ s|.*/||;
  25
  26 # A global destructor to close standard output with error checking.
  27 sub END
  28 {
  29   defined fileno STDOUT
  30     or return;
  31   close STDOUT
  32     and return;
  33   warn "$ME: closing standard output: $!\n";
  34   $? ||= 1;
  35 }
  36
  37 sub is_prime ($)
  38 {
  39   my ($n) = @_;
  40   use integer;
  41
  42   $n == 2
  43     and return 1;
  44
  45   my $d = 2;
  46   my $w = 1;
  47   while (1)
  48     {
  49       my $q = $n / $d;
  50       $n == $q * $d
  51         and return 0;
  52       $d += $w;
  53       $q < $d
  54         and last;
  55       $w = 2;
  56     }
  57   return 1;
  58 }
  59
  60 {
  61   @ARGV == 1
  62     or die "$ME: missing argument\n";
  63
  64   my $wheel_size = $ARGV[0];
  65
  66   my @primes = (2);
  67   my $product = $primes[0];
  68   my $n_primes = 1;
  69   for (my $i = 3; ; $i += 2)
  70     {
  71       if (is_prime $i)
  72         {
  73           push @primes, $i;
  74           $product *= $i;
  75           ++$n_primes == $wheel_size
  76             and last;
  77         }
  78     }
  79
  80   my $ws_m1 = $wheel_size - 1;
  81   print <<EOF;
  82 /* The first $ws_m1 elements correspond to the incremental offsets of the
  83    first $wheel_size primes (@primes).  The $wheel_size(th) element is the
  84    difference between that last prime and the next largest integer
  85    that is not a multiple of those primes.  The remaining numbers
  86    define the wheel.  For more information, see
  87    http://www.utm.edu/research/primes/glossary/WheelFactorization.html.  */
  88 EOF
  89
  90   my @increments;
  91   my $prev = 2;
  92   for (my $i = 3; ; $i += 2)
  93     {
  94       my $rel_prime = 1;
  95       foreach my $divisor (@primes)
  96         {
  97           $i != $divisor && $i % $divisor == 0
  98             and $rel_prime = 0;
  99         }
 100
 101       if ($rel_prime)
 102         {
 103           #warn $i, ' ', $i - $prev, "\n";
 104           push @increments, $i - $prev;
 105           $prev = $i;
 106
 107           $product + 1 < $i
 108             and last;
 109         }
 110     }
 111
 112   print join (",\n", @increments), "\n";
 113
 114   exit 0;
 115 }