src/wheel-gen.pl

   1 #!/usr/bin/perl -w
   2 # -*- perl -*-
   3
   4 eval 'exec /usr/bin/perl -S $0 ${1+"$@"}'
   5     if 0;
   6
   7 use strict;
   8 (my $program_name = $0) =~ s|.*/||;
   9
  10 sub END
  11 {
  12   use POSIX qw (_exit);
  13   # This is required if the code might send any output to stdout
  14   # E.g., even --version or --help.  So it's best to do it unconditionally.
  15   close STDOUT
  16     or (warn "$program_name: closing standard output: $!\n"), _exit (1);
  17 }
  18
  19 sub is_prime ($)
  20 {
  21   my ($n) = @_;
  22   use integer;
  23
  24   $n == 2
  25     and return 1;
  26
  27   my $d = 2;
  28   my $w = 1;
  29   while (1)
  30     {
  31       my $q = $n / $d;
  32       $n == $q * $d
  33         and return 0;
  34       $d += $w;
  35       $q < $d
  36         and last;
  37       $w = 2;
  38     }
  39   return 1;
  40 }
  41
  42 {
  43   @ARGV == 1
  44     or die "$program_name: missing argument\n";
  45
  46   my $wheel_size = $ARGV[0];
  47
  48   my @primes = (2);
  49   my $product = $primes[0];
  50   my $n_primes = 1;
  51   for (my $i = 3; ; $i += 2)
  52     {
  53       if (is_prime $i)
  54         {
  55           push @primes, $i;
  56           $product *= $i;
  57           ++$n_primes == $wheel_size
  58             and last;
  59         }
  60     }
  61
  62   my $ws_m1 = $wheel_size - 1;
  63   print <<EOF;
  64 /* The first $ws_m1 elements correspond to the incremental offsets of the
  65    first $wheel_size primes (@primes).  The $wheel_size(th) element is the
  66    difference between that last prime and the next largest integer
  67    that is not a multiple of those primes.  The remaining numbers
  68    define the wheel.  For more information, see
  69    http://www.utm.edu/research/primes/glossary/WheelFactorization.html.  */
  70 EOF
  71
  72   my @increments;
  73   my $prev = 2;
  74   for (my $i = 3; ; $i += 2)
  75     {
  76       my $rel_prime = 1;
  77       foreach my $divisor (@primes)
  78         {
  79           $i != $divisor && $i % $divisor == 0
  80             and $rel_prime = 0;
  81         }
  82
  83       if ($rel_prime)
  84         {
  85           #warn $i, ' ', $i - $prev, "\n";
  86           push @increments, $i - $prev;
  87           $prev = $i;
  88
  89           $product + 1 < $i
  90             and last;
  91         }
  92     }
  93
  94   print join (",\n", @increments), "\n";
  95
  96   exit 0;
  97 }