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
19eval 'exec /usr/bin/perl -S $0 ${1+"$@"}'
20  if 0;
21
22use strict;
23(my $ME = $0) =~ s|.*/||;
24
25# A global destructor to close standard output with error checking.
26sub 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
36sub 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.  */
87EOF
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}
115