1 /*@ Implementation of prime.h.
2 *
3 * Copyright (c) 2001 - 2020 Steffen (Daode) Nurpmeso <steffen@sdaoden.eu>.
4 * SPDX-License-Identifier: ISC
5 *
6 * Permission to use, copy, modify, and/or distribute this software for any
7 * purpose with or without fee is hereby granted, provided that the above
8 * copyright notice and this permission notice appear in all copies.
9 *
10 * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
11 * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
12 * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
13 * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
14 * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
15 * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
16 * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
17 */
18 #undef su_FILE
19 #define su_FILE su_prime
20 #define su_SOURCE
21 #define su_SOURCE_PRIME
22
23 #include "su/code.h"
24
25 #include "su/prime.h"
26 #include "su/code-in.h"
27
28 /* Collected & merged from 'GLIB 1's 'gprimes.c' and 'GNU STL's 'hashtable'
29 * (around Y2K+1 or so) */
30 static u32 const a_prime_lookup[] = {
31 0x00000002, 0x00000005, 0x0000000B,
32 0x00000017, 0x0000002F, 0x00000061, 0x0000009D,
33 0x0000011B, 0x0000022D, 0x00000337, 0x000004D5, 0x00000741, 0x00000AD9,
34 0x00001051, 0x00001867, 0x0000249B, 0x000036E9, 0x00005261, 0x00007B8B,
35 0x0000B947,
36 0x000115E7, 0x0001A0E1, 0x00027149, 0x0003A9E5, 0x00057EE3, 0x00083E39,
37 0x000C5D67,
38 0x00128C09, 0x001BD1FF, 0x0029BB13, 0x003E988B, 0x005DE4C1, 0x008CD721,
39 0x00D342AB,
40 0x01800013, 0x03000005, 0x06000017, 0x0C000013,
41 0x18000005
42 };
43 #define a_PRIME_FIRST &a_prime_lookup[0]
44 #define a_PRIME_MIDDLE &a_prime_lookup[NELEM(a_prime_lookup) / 2]
45 #define a_PRIME_LAST &a_prime_lookup[NELEM(a_prime_lookup) - 1]
46
47 static u64 const a_prime_max = su_U64_C(0xFFFFFFFFFFFFFFFD);
48
49 /* */
50 static boole a_prime_is_pseudo(u64 no);
51 static boole a_prime_is_real(u64 no);
52
53 /* */
54 static u64 a_prime_calc(u64 no, sz add, boole pseudo_ok);
55
56 static boole
a_prime_is_pseudo(u64 no)57 a_prime_is_pseudo(u64 no){
58 boole rv;
59 NYD_IN;
60
61 switch(no){
62 case 2:
63 case 3:
64 case 5:
65 case 7:
66 rv = TRU1; break;
67 case 0:
68 case 1:
69 case 4:
70 case 6:
71 rv = FAL0; break;
72 default:
73 rv = ((no & 1) && (no % 3) && (no % 5) && (no % 7));
74 break;
75 }
76 NYD_OU;
77 return rv;
78 }
79
80 static boole
a_prime_is_real(u64 no)81 a_prime_is_real(u64 no){ /* TODO brute force yet (at least Miller-Rabin?) */
82 /* no is pseudo! */
83 union {uz x; u64 x64; boole rv;} u;
84 NYD_IN;
85
86 if(no < 2)
87 goto jfal;
88 if(no != 2)
89 for(u.x64 = 3; u.x64 < no; u.x64 += 2)
90 if(no % u.x64 == 0)
91 goto jfal;
92
93 u.rv = TRU1;
94 jleave:
95 NYD_OU;
96 return u.rv;
97 jfal:
98 u.rv = FAL0;
99 goto jleave;
100 }
101
102 static u64
a_prime_calc(u64 no,sz add,boole pseudo_ok)103 a_prime_calc(u64 no, sz add, boole pseudo_ok){
104 NYD_IN;
105
106 /* Primes are all odd (except 2 but that is NEVER evaluated in here) */
107 if(!(no & 1)){
108 no += add;
109 add <<= 1;
110 goto jdiver;
111 }
112
113 add <<= 1;
114 for(;;){
115 no += add;
116 jdiver:
117 if(!a_prime_is_pseudo(no))
118 continue;
119 if(pseudo_ok || a_prime_is_real(no))
120 break;
121 }
122 NYD_OU;
123 return no;
124 }
125
126 boole
su_prime_is_prime(u64 no,boole allowpseudo)127 su_prime_is_prime(u64 no, boole allowpseudo){
128 boole rv;
129 NYD_IN;
130
131 rv = a_prime_is_pseudo(no);
132 if(rv && !allowpseudo)
133 rv = a_prime_is_real(no);
134 NYD_OU;
135 return rv;
136 }
137
138 u64
su_prime_get_former(u64 no,boole allowpseudo)139 su_prime_get_former(u64 no, boole allowpseudo){
140 NYD_IN;
141
142 if(no <= 2 + 1)
143 no = 2;
144 else if(no > a_prime_max)
145 no = a_prime_max;
146 else
147 no = a_prime_calc(no, -1, allowpseudo);
148 NYD_OU;
149 return no;
150 }
151
152 u64
su_prime_get_next(u64 no,boole allowpseudo)153 su_prime_get_next(u64 no, boole allowpseudo){
154 NYD_IN;
155
156 if(no < 2)
157 no = 2;
158 else if(no >= a_prime_max - 2)
159 no = a_prime_max;
160 else
161 no = a_prime_calc(no, +1, allowpseudo);
162 NYD_OU;
163 return no;
164 }
165
166 u32
su_prime_lookup_former(u32 no)167 su_prime_lookup_former(u32 no){
168 u32 const *cursor;
169 NYD_IN;
170
171 cursor = ((no < *a_PRIME_MIDDLE) ? a_PRIME_MIDDLE - 1 : a_PRIME_LAST);
172 while(*cursor >= no && --cursor > a_PRIME_FIRST)
173 ;
174 NYD_OU;
175 return *cursor;
176 }
177
178 u32
su_prime_lookup_next(u32 no)179 su_prime_lookup_next(u32 no){
180 u32 const *cursor;
181 NYD_IN;
182
183 cursor = ((no > *a_PRIME_MIDDLE) ? a_PRIME_MIDDLE + 1 : a_PRIME_FIRST);
184 while(*cursor <= no && ++cursor < a_PRIME_LAST)
185 ;
186 NYD_OU;
187 return *cursor;
188 }
189
190 #include "su/code-ou.h"
191 /* s-it-mode */
192