1 /* $OpenBSD: factor.c,v 1.11 2001/10/24 14:32:56 deraadt Exp $ */ 2 /* $NetBSD: factor.c,v 1.5 1995/03/23 08:28:07 cgd Exp $ */ 3 4 /* 5 * Copyright (c) 1989, 1993 6 * The Regents of the University of California. All rights reserved. 7 * 8 * This code is derived from software contributed to Berkeley by 9 * Landon Curt Noll. 10 * 11 * Redistribution and use in source and binary forms, with or without 12 * modification, are permitted provided that the following conditions 13 * are met: 14 * 1. Redistributions of source code must retain the above copyright 15 * notice, this list of conditions and the following disclaimer. 16 * 2. Redistributions in binary form must reproduce the above copyright 17 * notice, this list of conditions and the following disclaimer in the 18 * documentation and/or other materials provided with the distribution. 19 * 3. All advertising materials mentioning features or use of this software 20 * must display the following acknowledgement: 21 * This product includes software developed by the University of 22 * California, Berkeley and its contributors. 23 * 4. Neither the name of the University nor the names of its contributors 24 * may be used to endorse or promote products derived from this software 25 * without specific prior written permission. 26 * 27 * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND 28 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 29 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE 30 * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE 31 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL 32 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS 33 * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) 34 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT 35 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY 36 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF 37 * SUCH DAMAGE. 38 */ 39 40 #ifndef lint 41 static char copyright[] = 42 "@(#) Copyright (c) 1989, 1993\n\ 43 The Regents of the University of California. All rights reserved.\n"; 44 #endif /* not lint */ 45 46 #ifndef lint 47 #if 0 48 static char sccsid[] = "@(#)factor.c 8.4 (Berkeley) 5/4/95"; 49 #else 50 static char rcsid[] = "$OpenBSD: factor.c,v 1.11 2001/10/24 14:32:56 deraadt Exp $"; 51 #endif 52 #endif /* not lint */ 53 54 /* 55 * factor - factor a number into primes 56 * 57 * By: Landon Curt Noll chongo@toad.com, ...!{sun,tolsoft}!hoptoad!chongo 58 * 59 * chongo <for a good prime call: 391581 * 2^216193 - 1> /\oo/\ 60 * 61 * usage: 62 * factor [number] ... 63 * 64 * The form of the output is: 65 * 66 * number: factor1 factor1 factor2 factor3 factor3 factor3 ... 67 * 68 * where factor1 < factor2 < factor3 < ... 69 * 70 * If no args are given, the list of numbers are read from stdin. 71 */ 72 73 #include <sys/types.h> 74 #include <err.h> 75 #include <ctype.h> 76 #include <errno.h> 77 #include <limits.h> 78 #include <math.h> 79 #include <stdio.h> 80 #include <stdlib.h> 81 #include <string.h> 82 #include <unistd.h> 83 84 #include "primes.h" 85 86 /* 87 * prime[i] is the (i+1)th prime. 88 * 89 * We are able to sieve 2^32-1 because this byte table yields all primes 90 * up to 65537 and 65537^2 > 2^32-1. 91 */ 92 extern const ubig prime[]; 93 extern const ubig *pr_limit; /* largest prime in the prime array */ 94 extern const char pattern[]; 95 extern const int pattern_size; 96 97 void pr_fact __P((u_int64_t)); /* print factors of a value */ 98 void pr_bigfact __P((u_int64_t)); 99 void usage __P((void)); 100 101 int 102 main(argc, argv) 103 int argc; 104 char *argv[]; 105 { 106 u_int64_t val; 107 int ch; 108 char *p, buf[100]; /* > max number of digits. */ 109 110 /* revoke privs */ 111 setegid(getgid()); 112 setgid(getgid()); 113 114 while ((ch = getopt(argc, argv, "")) != -1) 115 switch (ch) { 116 case '?': 117 default: 118 usage(); 119 } 120 argc -= optind; 121 argv += optind; 122 123 /* No args supplied, read numbers from stdin. */ 124 if (argc == 0) 125 for (;;) { 126 if (fgets(buf, sizeof(buf), stdin) == NULL) { 127 if (ferror(stdin)) 128 err(1, "stdin"); 129 exit (0); 130 } 131 if (*(p = buf + strlen(buf) - 1) == '\n') 132 *p = '\0'; 133 for (p = buf; isblank(*p); ++p); 134 if (*p == '\0') 135 continue; 136 if (*p == '-') 137 errx(1, "negative numbers aren't permitted."); 138 errno = 0; 139 val = strtouq(buf, &p, 10); 140 if (errno) 141 err(1, "%s", buf); 142 for (; isblank(*p); ++p); 143 if (*p != '\0') 144 errx(1, "%s: illegal numeric format.", buf); 145 pr_fact(val); 146 } 147 /* Factor the arguments. */ 148 else 149 for (; *argv != NULL; ++argv) { 150 if (argv[0][0] == '-') 151 errx(1, "negative numbers aren't permitted."); 152 errno = 0; 153 val = strtouq(argv[0], &p, 10); 154 if (errno) 155 err(1, "%s", argv[0]); 156 if (*p != '\0') 157 errx(1, "%s: illegal numeric format.", argv[0]); 158 pr_fact(val); 159 } 160 exit(0); 161 } 162 163 /* 164 * pr_fact - print the factors of a number 165 * 166 * If the number is 0 or 1, then print the number and return. 167 * If the number is < 0, print -1, negate the number and continue 168 * processing. 169 * 170 * Print the factors of the number, from the lowest to the highest. 171 * A factor will be printed multiple times if it divides the value 172 * multiple times. 173 * 174 * Factors are printed with leading tabs. 175 */ 176 void 177 pr_fact(val) 178 u_int64_t val; /* Factor this value. */ 179 { 180 const ubig *fact; /* The factor found. */ 181 182 /* Firewall - catch 0 and 1. */ 183 if (val == 0) /* Historical practice; 0 just exits. */ 184 exit(0); 185 if (val == 1) { 186 (void)printf("1: 1\n"); 187 return; 188 } 189 190 /* Factor value. */ 191 (void)printf("%llu:", val); 192 fflush(stdout); 193 for (fact = &prime[0]; val > 1; ++fact) { 194 /* Look for the smallest factor. */ 195 do { 196 if (val % (long)*fact == 0) 197 break; 198 } while (++fact <= pr_limit); 199 200 /* Watch for primes larger than the table. */ 201 if (fact > pr_limit) { 202 if (val > BIG) 203 pr_bigfact(val); 204 else 205 (void)printf(" %llu", val); 206 break; 207 } 208 209 /* Divide factor out until none are left. */ 210 do { 211 (void)printf(" %lu", (unsigned long) *fact); 212 val /= (long)*fact; 213 } while ((val % (long)*fact) == 0); 214 215 /* Let the user know we're doing something. */ 216 (void)fflush(stdout); 217 } 218 (void)putchar('\n'); 219 } 220 221 222 /* At this point, our number may have factors greater than those in primes[]; 223 * however, we can generate primes up to 32 bits (see primes(6)), which is 224 * sufficient to factor a 64-bit quad. 225 */ 226 void 227 pr_bigfact(val) 228 u_int64_t val; /* Factor this value. */ 229 { 230 ubig start, stop, factor; 231 char *q; 232 const ubig *p; 233 ubig fact_lim, mod; 234 char *tab_lim; 235 char table[TABSIZE]; /* Eratosthenes sieve of odd numbers */ 236 237 start = *pr_limit + 2; 238 stop = (ubig)sqrt((double)val); 239 if ((stop & 0x1) == 0) 240 stop++; 241 /* 242 * Following code barely modified from that in primes(6) 243 * 244 * we shall sieve a bytemap window, note primes and move the window 245 * upward until we pass the stop point 246 */ 247 while (start < stop) { 248 /* 249 * factor out 3, 5, 7, 11 and 13 250 */ 251 /* initial pattern copy */ 252 factor = (start%(2*3*5*7*11*13))/2; /* starting copy spot */ 253 memcpy(table, &pattern[factor], pattern_size-factor); 254 /* main block pattern copies */ 255 for (fact_lim = pattern_size - factor; 256 fact_lim + pattern_size <= TABSIZE; fact_lim += pattern_size) { 257 memcpy(&table[fact_lim], pattern, pattern_size); 258 } 259 /* final block pattern copy */ 260 memcpy(&table[fact_lim], pattern, TABSIZE-fact_lim); 261 262 if (stop-start > TABSIZE+TABSIZE) { 263 tab_lim = &table[TABSIZE]; /* sieve it all */ 264 fact_lim = (int)sqrt( 265 (double)(start)+TABSIZE+TABSIZE+1.0); 266 } else { 267 tab_lim = &table[(stop - start)/2]; /* partial sieve */ 268 fact_lim = (int)sqrt((double)(stop) + 1.0); 269 } 270 /* sieve for factors >= 17 */ 271 factor = 17; /* 17 is first prime to use */ 272 p = &prime[7]; /* 19 is next prime, pi(19)=7 */ 273 do { 274 /* determine the factor's initial sieve point */ 275 mod = start % factor; 276 if (mod & 0x1) 277 q = &table[(factor-mod)/2]; 278 else 279 q = &table[mod ? factor-(mod/2) : 0]; 280 /* sieve for our current factor */ 281 for ( ; q < tab_lim; q += factor) { 282 *q = '\0'; /* sieve out a spot */ 283 } 284 } while ((factor=(ubig)(*(p++))) <= fact_lim); 285 286 /* 287 * use generated primes 288 */ 289 for (q = table; q < tab_lim; ++q, start+=2) { 290 if (*q) { 291 if (val % start == 0) { 292 do { 293 (void)printf(" %lu", (unsigned long) start); 294 val /= start; 295 } while ((val % start) == 0); 296 (void)fflush(stdout); 297 stop = (ubig)sqrt((double)val); 298 if ((stop & 0x1) == 0) 299 stop++; 300 } 301 } 302 } 303 } 304 if (val > 1) 305 printf(" %llu", val); 306 } 307 308 309 void 310 usage() 311 { 312 (void)fprintf(stderr, "usage: factor [value ...]\n"); 313 exit (0); 314 } 315