I don't know how I got myself into such a heated discussion.
The following is a copy of my original posting. I don't understand
how it has led to the responses it has generated.
****************** ORIGINAL POSTING BEGIN *******************
I have discovered an extension to the Sieve of Eratosthenes
that not only finds the primes, but also finds the prime
factorization of the composites in canonical form.
The algorithm requires no explicit multiplication or
division, only loops and increments, the same as the "Sieve".
Is such an algorithm already known? I have never
been able to find a reference to such an algorithm.
I realize that the algorithm probably does not have
any practical application. But I think such an algorithm
would be of interest to the study of the "Fundamental
Theorem of Arithmetic" since it is so simple, just as
the "Sieve" is.
I just want to know if it is something new, which I can't
imagine it is. But I cannot find a reference to it.
Any reply would be appreciated.
----------
I have been a member of GIMPS since 1999 and have Prime95
running on 3 or 4 computers most of the time, but no "joy" yet.
- rlindley
******************** ORIGINAL POSTING END ************************
As is clearly stated in the original posting above, I just want to
know if such an algorithm is common knowledge. Apparently it
is not, so the following is my FORTRAN implementation. As I have
stated in other postings, I realize the FORTRAN code is not the best
and I know that there are many optimizations possible. And I have made
no claims that the algorithm has any practical value. It certainly will
not factor 2^p -1 for very large values of p. As is stated in the
comments of the code, the program, as written, is bounded by
machine constraints (basically 32 bit integers). It is nothing more
than a simple algorithm to find the prime factorization of all integers
up to some bound. Just as the "Sieve" finds primes up to some
bound.
The following FORTRAN source code compiles and executes ok for me
under Windows XP using the open source g77 compiler. (No guarantees
for you.)
Following the source code is a sample output file to demonstrate that
the algorithm does work, at least up to the parameterized bounds specified.
Increase them to 2^31 -1 or whatever to test it further. Or better yet,
provide a mathematical proof that it works ad infinitum, which I am
not capable of providing. Or equally as good, provide a counter example.
******************** FORTRAN SOURCE CODE BEGIN ************************
program primefactors
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Author: Roger Lindley
!
! Date: October 19, 2004
!
! Program to find the prime factorization of all integers up to some
! bound within machine constraints.
!
! The factorization is in the form
! N = (2**e0)*(p1**e1)*(p2**e2)*(p3**e3)... for even integers
! or
! N = (p1**e1)*(p2**e2)*(p3**e3)... for odd integers
! where N is a positive integer, the p's are odd primes
! and the the e's are positive integers.
!
! The algorithm is an extension of the Sieve of Eratosthenes and was
! unknown to the author.
!
! It may be interesting to note that the algorithm requires no explicit
! multiplication or division, only loops and increments.
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! The single array "factors" is used to contain the required data. (A single
! array (instead of multiple arrays) is used in an effort to reduce paging
! by keeping closely related data close together in memory.)
!
! The structure of "factors" is as follows:
!
! Dimension: factors(max_ints,max_factors,2) where max_ints is the
! number of integers to
! consider and max_factors
! is the maximum number of
! prime factors allowed for
! each integer.
!
! Upon completion of the program the array "factors" will contain
! the following:
!
! factors(i,1,1)(i=1,max_ints) : the integers up to some bound
! factors(i,1,2)(i=1,max_ints) : the index of the storage location
! of the largest prime factor of this
! integer.
! When this value is 1, then the
! integer is prime.
! factors(i,j,1)(j=2,factors(i,1,2)) : the prime factors of this
! integer
! factors(i,j,2)(j=2,factors(i,1,2)) : the exponents of each prime
! factor
!
! Example:
!
! Consider the integer 45;
! 45 is in the 45th "row"
! So,
! factors(45,1,1) = 45
! factors(45,1,2) = 3 the largest prime factor of 45 is stored
! in the third element
! factors(45,2,1) = 3 3 is a prime factor of 45
! factors(45,2,2) = 2 the exponent of 3 is 2, i.e.., 45
! has a factor of 3**2 = 9
! factors(45,3,1) = 5 5 is a prime factor of 45
! factors(45,3,2) = 1 the exponent of 5 is 1
!
! Thus, all this data means that
! 45 = (3**2)*(5**1)
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!!!!!!!!!!!!!!!!! Data definitions and declarations. !!!!!!!!!!!!!!!!!
!
implicit none
integer*4 max_ints,max_factors
parameter (max_ints = 2**8) ! these parameters are
parameter (max_factors = 3) ! purposely set low for
! demonstration purposes
integer*4 factors(max_ints,max_factors,2) ! described above
integer*4 i,j,k,ii,int,maxindex,primefactor ! local variables
integer*2 out_lun ! logical unit of output
parameter(out_lun=11) ! file
character*20 out_file ! output file name
parameter (out_file = 'factors.dat')
!!!!!!!!!!!!!!!!! Begin executable code. !!!!!!!!!!!!!!!!!!!!!!!
do i = 1,max_ints
factors(i,1,1) = i ! load integers into array
factors(i,1,2) = 1 ! mark all as being
! prime
end do
open(unit=out_lun,name=out_file) ! create the output file
write(out_lun,*)1,0 ! 1 is a special case
do i = 2,max_ints ! must look at each integer
! can't get away with just
! looking at those up to
! the square root of the
! largest
int = factors(i,1,1) ! save the integer
maxindex = factors(i,1,2) ! and the index of it's
! largest prime in local
! variables to avoid
! paging
if(maxindex .eq. 1) then ! if no factors then this
! integer is prime
do ii = i + int,max_ints,int ! sieve of eratosthenes
! loop over multiples of
! this prime
factors(ii,1,2) = factors(ii,1,2) + 1 ! bump factor index
! this also marks it
! composite
if (factors(ii,1,2).le. max_factors) then ! make sure enough space
factors(ii,factors(ii,1,2),1) = int ! load prime factor
factors(ii,factors(ii,1,2),2) = 1 ! set exponent
else
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! No empty slots to save factor in. This does not affect other results,
! so the program can continue. These integers could be written to a data
! base for further examination. No need to print an error message now.
! That can be done later.
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
end if
end do
else ! this integer is composite
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! This section is the extension of the Sieve of Eratosthenes.
!
! The reasoning here is that if the composite being examined only has one
! prime factor, then the composite is a power of that prime factor. Thus
! every multiple of this composite will have another factor of that prime
! factor. For example, the first such occurrence is 4. Since 4 is
! composite and only has one prime factor, 4 must be a power of that
! prime factor, so every multiple of 4, including 4 itself, will have
! another factor of that prime factor.
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
if (maxindex .eq. 2) then ! if this composite only
! has one prime factor
! then bump the exponents
primefactor = factors(i,2,1) ! save prime factor in
! local variable
do ii = i,max_ints,int ! loop over the multiples
! of this composite
do j = 2,factors(ii,1,2) ! look for prime factor in
! list
if (factors(ii,j,1) .eq. primefactor) then
factors(ii,j,2) = factors(ii,j,2) + 1 ! increment the exponent
exit ! done, so exit loop
end if
end do
end do
end if
end if
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! At this point the record for this integer can be written to a data base
! as it will not be needed any more. Then this memory need not be paged
! back in again.
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
if(factors(i,1,2) .le. max_factors) then
write(out_lun,*)((factors(i,j,k),k=1,2),j=1,factors(i,1,2))
else
write(out_lun,*) "Incomplete *******",&
((factors(i,j,k),k=1,2),j=1,max_factors)
end if
end do
close(unit=out_lun)
stop
end
******************** FORTRAN SOURCE CODE END ************************
I imported the above source code from a file created by the Crimson
editor and the tabs did not import exactly right, but it is still readable.
The comments just don't line up quite right.
The sample output file follows to demonstrate that the program works
for at least a little ways. Refer to the source code heading comments
for an explanation of the output.
******************** SAMPLE OUTPUT BEGIN ************************
1 0
2 1
3 1
4 2 2 2
5 1
6 3 2 1 3 1
7 1
8 2 2 3
9 2 3 2
10 3 2 1 5 1
11 1
12 3 2 2 3 1
13 1
14 3 2 1 7 1
15 3 3 1 5 1
16 2 2 4
17 1
18 3 2 1 3 2
19 1
20 3 2 2 5 1
21 3 3 1 7 1
22 3 2 1 11 1
23 1
24 3 2 3 3 1
25 2 5 2
26 3 2 1 13 1
27 2 3 3
28 3 2 2 7 1
29 1
Incomplete ******* 30 4 2 1 3 1
31 1
32 2 2 5
33 3 3 1 11 1
34 3 2 1 17 1
35 3 5 1 7 1
36 3 2 2 3 2
37 1
38 3 2 1 19 1
39 3 3 1 13 1
40 3 2 3 5 1
41 1
Incomplete ******* 42 4 2 1 3 1
43 1
44 3 2 2 11 1
45 3 3 2 5 1
46 3 2 1 23 1
47 1
48 3 2 4 3 1
49 2 7 2
50 3 2 1 5 2
51 3 3 1 17 1
52 3 2 2 13 1
53 1
54 3 2 1 3 3
55 3 5 1 11 1
56 3 2 3 7 1
57 3 3 1 19 1
58 3 2 1 29 1
59 1
Incomplete ******* 60 4 2 2 3 1
61 1
62 3 2 1 31 1
63 3 3 2 7 1
64 2 2 6
65 3 5 1 13 1
Incomplete ******* 66 4 2 1 3 1
67 1
68 3 2 2 17 1
69 3 3 1 23 1
Incomplete ******* 70 4 2 1 5 1
71 1
72 3 2 3 3 2
73 1
74 3 2 1 37 1
75 3 3 1 5 2
76 3 2 2 19 1
77 3 7 1 11 1
Incomplete ******* 78 4 2 1 3 1
79 1
80 3 2 4 5 1
81 2 3 4
82 3 2 1 41 1
83 1
Incomplete ******* 84 4 2 2 3 1
85 3 5 1 17 1
86 3 2 1 43 1
87 3 3 1 29 1
88 3 2 3 11 1
89 1
Incomplete ******* 90 4 2 1 3 2
91 3 7 1 13 1
92 3 2 2 23 1
93 3 3 1 31 1
94 3 2 1 47 1
95 3 5 1 19 1
96 3 2 5 3 1
97 1
98 3 2 1 7 2
99 3 3 2 11 1
100 3 2 2 5 2
101 1
Incomplete ******* 102 4 2 1 3 1
103 1
104 3 2 3 13 1
Incomplete ******* 105 4 3 1 5 1
106 3 2 1 53 1
107 1
108 3 2 2 3 3
109 1
Incomplete ******* 110 4 2 1 5 1
111 3 3 1 37 1
112 3 2 4 7 1
113 1
Incomplete ******* 114 4 2 1 3 1
115 3 5 1 23 1
116 3 2 2 29 1
117 3 3 2 13 1
118 3 2 1 59 1
119 3 7 1 17 1
Incomplete ******* 120 4 2 3 3 1
121 2 11 2
122 3 2 1 61 1
123 3 3 1 41 1
124 3 2 2 31 1
125 2 5 3
Incomplete ******* 126 4 2 1 3 2
127 1
128 2 2 7
129 3 3 1 43 1
Incomplete ******* 130 4 2 1 5 1
131 1
Incomplete ******* 132 4 2 2 3 1
133 3 7 1 19 1
134 3 2 1 67 1
135 3 3 3 5 1
136 3 2 3 17 1
137 1
Incomplete ******* 138 4 2 1 3 1
139 1
Incomplete ******* 140 4 2 2 5 1
141 3 3 1 47 1
142 3 2 1 71 1
143 3 11 1 13 1
144 3 2 4 3 2
145 3 5 1 29 1
146 3 2 1 73 1
147 3 3 1 7 2
148 3 2 2 37 1
149 1
Incomplete ******* 150 4 2 1 3 1
151 1
152 3 2 3 19 1
153 3 3 2 17 1
Incomplete ******* 154 4 2 1 7 1
155 3 5 1 31 1
Incomplete ******* 156 4 2 2 3 1
157 1
158 3 2 1 79 1
159 3 3 1 53 1
160 3 2 5 5 1
161 3 7 1 23 1
162 3 2 1 3 4
163 1
164 3 2 2 41 1
Incomplete ******* 165 4 3 1 5 1
166 3 2 1 83 1
167 1
Incomplete ******* 168 4 2 3 3 1
169 2 13 2
Incomplete ******* 170 4 2 1 5 1
171 3 3 2 19 1
172 3 2 2 43 1
173 1
Incomplete ******* 174 4 2 1 3 1
175 3 5 2 7 1
176 3 2 4 11 1
177 3 3 1 59 1
178 3 2 1 89 1
179 1
Incomplete ******* 180 4 2 2 3 2
181 1
Incomplete ******* 182 4 2 1 7 1
183 3 3 1 61 1
184 3 2 3 23 1
185 3 5 1 37 1
Incomplete ******* 186 4 2 1 3 1
187 3 11 1 17 1
188 3 2 2 47 1
189 3 3 3 7 1
Incomplete ******* 190 4 2 1 5 1
191 1
192 3 2 6 3 1
193 1
194 3 2 1 97 1
Incomplete ******* 195 4 3 1 5 1
196 3 2 2 7 2
197 1
Incomplete ******* 198 4 2 1 3 2
199 1
200 3 2 3 5 2
201 3 3 1 67 1
202 3 2 1 101 1
203 3 7 1 29 1
Incomplete ******* 204 4 2 2 3 1
205 3 5 1 41 1
206 3 2 1 103 1
207 3 3 2 23 1
208 3 2 4 13 1
209 3 11 1 19 1
Incomplete ******* 210 5 2 1 3 1
211 1
212 3 2 2 53 1
213 3 3 1 71 1
214 3 2 1 107 1
215 3 5 1 43 1
216 3 2 3 3 3
217 3 7 1 31 1
218 3 2 1 109 1
219 3 3 1 73 1
Incomplete ******* 220 4 2 2 5 1
221 3 13 1 17 1
Incomplete ******* 222 4 2 1 3 1
223 1
224 3 2 5 7 1
225 3 3 2 5 2
226 3 2 1 113 1
227 1
Incomplete ******* 228 4 2 2 3 1
229 1
Incomplete ******* 230 4 2 1 5 1
Incomplete ******* 231 4 3 1 7 1
232 3 2 3 29 1
233 1
Incomplete ******* 234 4 2 1 3 2
235 3 5 1 47 1
236 3 2 2 59 1
237 3 3 1 79 1
Incomplete ******* 238 4 2 1 7 1
239 1
Incomplete ******* 240 4 2 4 3 1
241 1
242 3 2 1 11 2
243 2 3 5
244 3 2 2 61 1
245 3 5 1 7 2
Incomplete ******* 246 4 2 1 3 1
247 3 13 1 19 1
248 3 2 3 31 1
249 3 3 1 83 1
250 3 2 1 5 3
251 1
Incomplete ******* 252 4 2 2 3 2
253 3 11 1 23 1
254 3 2 1 127 1
Incomplete ******* 255 4 3 1 5 1
256 2 2 8
******************** SAMPLE OUTPUT END ************************
Please don't hassle me about all the "incompletes" in the sample
output file. Those are for demonstration. The parameter "max_factors"
is only set to 3. It is not hard to precompute "max_factors" from
"max_ints" to make sure you don't run out of space. But nowadays
with all the new languages that dynamically allocate space it might
be possible to stretch the constraints that this program
is limited to.
Ok, there you have it. If it is all elementary and obvious, that is ok
with me. I have laid out what I know. I just had never seen it before.
So it was new to me.
- rlindley
_______________________________________________
Prime mailing list
[email protected]
http://hogranch.com/mailman/listinfo/prime
