summaryrefslogtreecommitdiff
path: root/src/wheel-gen.pl
blob: 5ed9841b0c668a77ee4cbdaaa231a5cb4991cdb6 (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
#!/usr/bin/perl -w
# Generate the spokes of a wheel, for wheel factorization.

# Copyright (C) 2001, 2005, 2009-2010 Free Software Foundation, Inc.

# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.

# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.

# You should have received a copy of the GNU General Public License
# along with this program.  If not, see <http://www.gnu.org/licenses/>.

eval 'exec /usr/bin/perl -S $0 ${1+"$@"}'
  if 0;

use strict;
(my $ME = $0) =~ s|.*/||;

# A global destructor to close standard output with error checking.
sub END
{
  defined fileno STDOUT
    or return;
  close STDOUT
    and return;
  warn "$ME: closing standard output: $!\n";
  $? ||= 1;
}

sub is_prime ($)
{
  my ($n) = @_;
  use integer;

  $n == 2
    and return 1;

  my $d = 2;
  my $w = 1;
  while (1)
    {
      my $q = $n / $d;
      $n == $q * $d
        and return 0;
      $d += $w;
      $q < $d
        and last;
      $w = 2;
    }
  return 1;
}

{
  @ARGV == 1
    or die "$ME: missing argument\n";

  my $wheel_size = $ARGV[0];

  my @primes = (2);
  my $product = $primes[0];
  my $n_primes = 1;
  for (my $i = 3; ; $i += 2)
    {
      if (is_prime $i)
        {
          push @primes, $i;
          $product *= $i;
          ++$n_primes == $wheel_size
            and last;
        }
    }

  my $ws_m1 = $wheel_size - 1;
  print <<EOF;
/* The first $ws_m1 elements correspond to the incremental offsets of the
   first $wheel_size primes (@primes).  The $wheel_size(th) element is the
   difference between that last prime and the next largest integer
   that is not a multiple of those primes.  The remaining numbers
   define the wheel.  For more information, see
   http://www.utm.edu/research/primes/glossary/WheelFactorization.html.  */
EOF

  my @increments;
  my $prev = 2;
  for (my $i = 3; ; $i += 2)
    {
      my $rel_prime = 1;
      foreach my $divisor (@primes)
        {
          $i != $divisor && $i % $divisor == 0
            and $rel_prime = 0;
        }

      if ($rel_prime)
        {
          #warn $i, ' ', $i - $prev, "\n";
          push @increments, $i - $prev;
          $prev = $i;

          $product + 1 < $i
            and last;
        }
    }

  print join (",\n", @increments), "\n";

  exit 0;
}