Line data Source code
1 : /* Simple expression parser for GMP-ECM.
2 :
3 : Copyright 2003, 2004, 2005, 2006, 2007, 2008, 2012 Jim Fougeron, Paul Zimmermann and Alexander Kruppa.
4 :
5 : This program is free software; you can redistribute it and/or modify
6 : it under the terms of the GNU General Public License as published by
7 : the Free Software Foundation; either version 3 of the License, or (at your
8 : option) any later version.
9 :
10 : This program is distributed in the hope that it will be useful, but
11 : WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
12 : or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for
13 : more details.
14 :
15 : You should have received a copy of the GNU General Public License
16 : along with this program; see the file COPYING. If not, see
17 : http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
18 : 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */
19 :
20 : #include <stdio.h>
21 : #include <stdlib.h>
22 : #include <string.h>
23 : #include <time.h>
24 : #include "ecm-ecm.h"
25 : #include "getprime_r.h"
26 :
27 : #ifdef HAVE_STRINGS_H
28 : # include <strings.h> /* for strncasecmp */
29 : #endif
30 :
31 : #ifdef HAVE_CTYPE_H
32 : # include <ctype.h>
33 : #endif
34 :
35 :
36 : /*****************************************************************
37 : * Syntax for this VERY simple recursive expression parser: *
38 : * *
39 : * ( or [ or { along with ) or ] or } are valid for grouping *
40 : * Normal "simple" operators: + - * / (. can be used for *) *
41 : * Module: n%m 345%11 *
42 : * Unary minus is supported: -n -500 *
43 : * Exponentation: n^m 2^500 *
44 : * Simple factorial: n! 53! == 1*2*3*4...*52*53 *
45 : * Multi-factorial: n!m 15!3 == 15.12.9.6.3 *
46 : * Simple Primorial: n# 11# == 2*3*5*7*11 *
47 : * Reduced Primorial: n#m 17#5 == 5.7.11.13.17 *
48 : * *
49 : * Supported functions: (case insensitive) *
50 : * Phi(n,x) *
51 : * GCD(m,n) *
52 : * U(p,q,n) *
53 : * primU(p,q,n) *
54 : * TODO: PhiL(k,n), PhiM(k,n) *
55 : * only for bases 2,3,5,6,7,10,11 (times a square) *
56 : * Note for developers: *
57 : * First k-1 arguments are passed as an mpz_t array *
58 : * *
59 : * NOTE Lines ending in a \ character are "joined" *
60 : * NOTE C++ // single line comments (rest of line is a comment) *
61 : * *
62 : ****************************************************************/
63 :
64 : /* value only used by the expression parser */
65 : static mpz_t t, mpOne;
66 : static char *expr_str;
67 :
68 : static void eval_power (mpz_t prior_n, mpz_t n,char op);
69 : static void eval_product (mpz_t prior_n, mpz_t n,char op);
70 : static void eval_sum (mpz_t prior_n, mpz_t n,char op);
71 : static int eval_Phi (mpz_t *params, mpz_t n);
72 : static int eval_PhiL (mpz_t *params, mpz_t n);
73 : static int eval_PhiM (mpz_t *params, mpz_t n);
74 : // static int eval_gcd (mpz_t *params, mpz_t n);
75 : static int eval_U (mpz_t *params, mpz_t n);
76 : static int eval_primU (mpz_t *params, mpz_t n);
77 : static int eval_2 (int bInFuncParams);
78 : static int aurif (mpz_t output, mpz_t n, mpz_t base, int sign);
79 :
80 : #if 0 /* strncasecmp is a required function in configure.in */
81 : #if defined (_MSC_VER) || defined (__MINGW32__)
82 : #define strncasecmp strnicmp
83 : #endif
84 : #endif
85 :
86 : /***************************************/
87 : /* Main expression evaluation function */
88 : /* This is the function that the app */
89 : /* calls to read the expression line */
90 : /***************************************/
91 2873 : int eval (mpcandi_t *n, FILE *fd, int primetest)
92 : {
93 : int ret;
94 2873 : int nMaxSize = 2000, nCurSize = 0;
95 : int c;
96 2873 : char *expr = (char *) malloc (nMaxSize + 1);
97 :
98 2873 : ASSERT_ALWAYS (expr != NULL);
99 2883 : JoinLinesLoop:
100 2883 : c = fgetc (fd);
101 : if (0)
102 : {
103 10 : ChompLine:
104 : do
105 90 : c = fgetc (fd);
106 90 : while (c != EOF && !IS_NEWLINE(c));
107 10 : if (IS_NEWLINE(c))
108 10 : goto JoinLinesLoop;
109 : }
110 :
111 130000 : while (c != EOF && !IS_NEWLINE(c) && c != ';')
112 : {
113 127127 : if (c == '/')
114 : {
115 : /* This might be a C++ // comment or it might be a / division operator.
116 : Check it out, and if it is a comment, then "eat it" */
117 148 : int peek_c = fgetc (fd);
118 148 : if (peek_c == '/')
119 : /* Got a C++ single line comment, so Chomp the line */
120 10 : goto ChompLine;
121 :
122 : /* Put the char back on the file, then allow the code to add the '/' char to the buffer */
123 138 : ungetc (peek_c, fd);
124 : }
125 :
126 : /* strip space and tabs out here, and then we DON'T have to mess with them in the rest of the parser */
127 127117 : if (!isspace (c) && c != '"' && c != '\'')
128 127107 : expr[nCurSize++] = (char) c;
129 :
130 127117 : if (nCurSize == nMaxSize)
131 : {
132 : char *cp;
133 10 : nMaxSize += nMaxSize / 2;
134 10 : cp = (char *) realloc (expr, nMaxSize + 1);
135 10 : ASSERT_ALWAYS (cp != NULL);
136 10 : expr = cp;
137 : }
138 127117 : c = fgetc (fd);
139 : }
140 2873 : expr[nCurSize] = 0;
141 2873 : if (!nCurSize)
142 0 : ret = 0;
143 : else
144 : {
145 2873 : if (c == ';')
146 157 : ungetc (c, fd);
147 2873 : mpz_init (t);
148 2873 : expr_str = expr;
149 2873 : ret = eval_2 (0);
150 2603 : if (ret)
151 : {
152 : char *s;
153 2603 : char *cpTmpExpr = expr;
154 2603 : s = mpz_get_str (NULL, 10, t);
155 2603 : if (!strcmp(s, cpTmpExpr))
156 2009 : cpTmpExpr = NULL;
157 2603 : ret = mpcandi_t_add_candidate (n, t, cpTmpExpr, primetest);
158 2603 : free (s); /* size strlen (s) + 1 */
159 : }
160 2603 : mpz_clear(t);
161 : }
162 2603 : free(expr);
163 2603 : return ret;
164 : }
165 :
166 181 : int eval_str (mpcandi_t *n, char *cp, int primetest, char **EndChar)
167 : {
168 : int ret;
169 181 : int nMaxSize=2000, nCurSize=0;
170 : char *c;
171 181 : char *expr = (char *) malloc(nMaxSize+1);
172 :
173 181 : ASSERT_ALWAYS (expr != NULL);
174 181 : c = cp;
175 181 : JoinLinesLoop:
176 181 : if (*c == '#')
177 : {
178 : do
179 0 : ++c;
180 0 : while (*c && !IS_NEWLINE(*c));
181 0 : if (IS_NEWLINE(*c))
182 0 : goto JoinLinesLoop;
183 : }
184 :
185 9533 : while (*c && !IS_NEWLINE(*c) && *c != ';')
186 : {
187 : /* strip space and tabs out here, and then we DON'T have to mess with them in the rest of the parser */
188 9352 : if (!isspace((int) *c) && *c != '"' && *c != '\'')
189 9352 : expr[nCurSize++] = *c;
190 9352 : if (nCurSize == nMaxSize)
191 : {
192 : char *cp;
193 0 : nMaxSize += nMaxSize / 2;
194 0 : cp = (char *) realloc (expr, nMaxSize + 1);
195 0 : ASSERT_ALWAYS (cp != NULL);
196 0 : expr = cp;
197 : }
198 9352 : ++c;
199 : }
200 181 : expr[nCurSize] = 0;
201 181 : if (!nCurSize)
202 0 : ret = 0;
203 : else
204 : {
205 181 : if (*c != ';')
206 181 : ++c;
207 181 : mpz_init(t);
208 181 : expr_str = expr;
209 181 : ret = eval_2(0);
210 171 : if (ret)
211 : {
212 : char *s;
213 171 : char *cpTmpExpr = expr;
214 171 : s = mpz_get_str (NULL, 10, t);
215 171 : if (!strcmp(s, cpTmpExpr))
216 152 : cpTmpExpr = NULL;
217 171 : ret = mpcandi_t_add_candidate(n, t, cpTmpExpr, primetest);
218 171 : free (s); /* size strlen (s) + 1 */
219 : }
220 171 : mpz_clear(t);
221 : }
222 171 : free(expr);
223 171 : if (EndChar && *EndChar)
224 0 : *EndChar = c;
225 171 : return ret;
226 : }
227 :
228 4806 : void eval_power (mpz_t prior_n, mpz_t n,char op)
229 : {
230 : #if defined (DEBUG_EVALUATOR)
231 : if ('#'==op || '^'==op || '!'==op || '@'==op || '$'==op)
232 : {
233 : fprintf (stderr, "eval_power ");
234 : mpz_out_str(stderr, 10, prior_n);
235 : fprintf (stderr, "%c", op);
236 : mpz_out_str(stderr, 10, n);
237 : fprintf (stderr, "\n");
238 : }
239 : #endif
240 :
241 4806 : if ('^'==op)
242 516 : mpz_pow_ui(n,prior_n,mpz_get_ui(n));
243 4290 : else if ('!'==op) /* simple factorial (syntax n! example: 7! == 1*2*3*4*5*6*7) */
244 10 : mpz_fac_ui(n,mpz_get_ui(n));
245 4280 : else if ('@'==op) /* Multi factorial (syntax n!prior_n
246 : Example: 15!3 == 15*12*9*6*3
247 : Note: 15!3 is substituted into 15@3 by the parser */
248 : {
249 : long nCur;
250 : unsigned long nDecr;
251 10 : nCur = mpz_get_si(prior_n);
252 10 : nDecr = mpz_get_ui(n);
253 10 : mpz_set_ui(n,1);
254 350 : while (nCur > 1)
255 : {
256 : /* This could be done much more efficiently (bunching mults using smaller "built-ins"), but I am not going to bother for now */
257 340 : mpz_mul_ui(n,n,nCur);
258 340 : nCur -= nDecr;
259 : }
260 : }
261 4270 : else if ('#'==op) /* simple primorial (syntax n# example: 11# == 2*3*5*7*11 */
262 : {
263 : unsigned long nMax;
264 : unsigned long p;
265 : prime_info_t prime_info;
266 :
267 10 : prime_info_init (prime_info);
268 10 : nMax = mpz_get_ui (n);
269 10 : mpz_set_ui (n, 1);
270 270 : for (p = 2; p <= nMax; p = getprime_mt (prime_info))
271 : /* This could be done much more efficiently (bunching mults using smaller "built-ins"), but I am not going to bother for now */
272 260 : mpz_mul_ui (n, n, p);
273 10 : prime_info_clear (prime_info); /* free the prime table */
274 : }
275 4260 : else if ('$'==op) /* reduced primorial (syntax n#prior_n example: 13#5 == (5*7*11*13) */
276 : {
277 : unsigned long p;
278 : unsigned long nMax;
279 : unsigned long nStart;
280 : prime_info_t prime_info;
281 :
282 10 : prime_info_init (prime_info);
283 10 : nMax = mpz_get_ui (prior_n);
284 10 : nStart = mpz_get_ui (n);
285 10 : mpz_set_ui (n, 1);
286 : /*printf ("Reduced-primorial %ld#%ld\n", nMax, nStart);*/
287 270 : for (p = 2; p <= nMax; p = getprime_mt (prime_info))
288 : {
289 260 : if (p >= nStart)
290 : /* This could be done much more efficiently (bunching mults using smaller "built-ins"), but I am not going to bother for now */
291 250 : mpz_mul_ui (n, n, p);
292 : }
293 10 : prime_info_clear (prime_info); /* free the prime table */
294 : }
295 4806 : }
296 :
297 : void
298 4260 : eval_product (mpz_t prior_n, mpz_t n, char op)
299 : {
300 : #if defined (DEBUG_EVALUATOR)
301 : if ('*'==op || '.'==op || '/'==op || '%'==op)
302 : {
303 : fprintf (stderr, "eval_product ");
304 : mpz_out_str(stderr, 10, prior_n);
305 : fprintf (stderr, "%c", op);
306 : mpz_out_str(stderr, 10, n);
307 : fprintf (stderr, "\n");
308 : }
309 : #endif
310 4260 : if ('*' == op || '.' == op)
311 52 : mpz_mul (n, prior_n, n);
312 4208 : else if ('/' == op)
313 : {
314 : mpz_t r;
315 138 : mpz_init (r);
316 138 : mpz_tdiv_qr (n, r, prior_n, n);
317 138 : if (mpz_cmp_ui (r, 0) != 0)
318 : {
319 10 : fprintf (stderr, "Parsing Error: inexact division\n");
320 10 : exit (EXIT_FAILURE);
321 : }
322 128 : mpz_clear (r);
323 : }
324 4070 : else if ('%' == op)
325 10 : mpz_tdiv_r (n, prior_n, n);
326 4250 : }
327 :
328 4050 : void eval_sum (mpz_t prior_n, mpz_t n,char op)
329 : {
330 : #if defined (DEBUG_EVALUATOR)
331 : if ('+'==op || '-'==op)
332 : {
333 : fprintf (stderr, "eval_sum ");
334 : mpz_out_str(stderr, 10, prior_n);
335 : fprintf (stderr, "%c", op);
336 : mpz_out_str(stderr, 10, n);
337 : fprintf (stderr, "\n");
338 : }
339 : #endif
340 :
341 4050 : if ('+' == op)
342 130 : mpz_add(n,prior_n,n);
343 3920 : else if ('-' == op)
344 437 : mpz_sub(n,prior_n,n);
345 4050 : }
346 :
347 150 : int eval_Phi (mpz_t* params, mpz_t n)
348 : {
349 : /* params[0]=exp, n=base */
350 : int factors[200];
351 150 : unsigned dwFactors=0, dw;
352 : unsigned long B;
353 : unsigned long p;
354 : mpz_t D, T, org_n;
355 : prime_info_t prime_info;
356 :
357 : /* deal with trivial cases first */
358 150 : if (mpz_cmp_ui (params[0], 0) == 0)
359 : {
360 0 : mpz_set_ui (n, 1);
361 0 : return 1;
362 : }
363 150 : if (mpz_cmp_ui (params[0], 0) < 0)
364 10 : return 0;
365 140 : if (mpz_cmp_ui (params[0], 1) == 0)
366 : {
367 10 : mpz_sub_ui (n, n, 1);
368 10 : return 1;
369 : }
370 130 : if (mpz_cmp_ui (params[0], 2) == 0)
371 : {
372 0 : mpz_add_ui (n, n, 1);
373 0 : return 1;
374 : }
375 130 : if (mpz_cmp_ui (n, 0) < 0)
376 : /* Convert to positive base; this is always valid when exp>=3 */
377 : {
378 20 : mpz_neg (n, n);
379 20 : if (mpz_congruent_ui_p (params[0], 1, 2))
380 : {
381 10 : mpz_mul_ui(params[0], params[0], 2);
382 : }
383 10 : else if (mpz_congruent_ui_p (params[0], 2, 4))
384 : {
385 10 : mpz_divexact_ui(params[0], params[0], 2);
386 : }
387 : }
388 130 : if (mpz_cmp_ui (n, 1) == 0)
389 : {
390 : /* return value is p if params[0] is prime power p^k, or 1 otherwise */
391 30 : int maxpower=mpz_sizeinbase(params[0], 2)+1;
392 30 : mpz_init (T);
393 350 : for (int power=maxpower; power>=1; --power)
394 : {
395 350 : if ( mpz_root (T, params[0], power) ) break;
396 : }
397 30 : int isPrime = mpz_probab_prime_p (T, PROBAB_PRIME_TESTS);
398 30 : mpz_set (n, isPrime ? T : mpOne);
399 30 : mpz_clear(T);
400 30 : return 1;
401 : }
402 :
403 :
404 : /* Ok, do the real h_primative work, since we are not one of the trivial case */
405 :
406 100 : if (mpz_fits_ulong_p (params[0]) == 0)
407 10 : return 0;
408 :
409 90 : B = mpz_get_ui (params[0]);
410 :
411 : /* Obtain the factors of B */
412 90 : prime_info_init (prime_info);
413 1040 : for (p = 2; p <= B; p = getprime_mt (prime_info))
414 : {
415 950 : if (B % p == 0)
416 : {
417 : /* Add the factor one time */
418 190 : factors[dwFactors++] = p;
419 : /* but be sure to totally remove it */
420 280 : do { B /= p; } while (B % p == 0);
421 : }
422 : }
423 90 : prime_info_clear (prime_info); /* free the prime tables */
424 90 : B = mpz_get_si (params[0]);
425 :
426 90 : mpz_init_set (org_n, n);
427 90 : mpz_set_ui (n, 1);
428 90 : mpz_init_set_ui (D, 1);
429 90 : mpz_init (T);
430 :
431 550 : for(dw=0;(dw<(1U<<dwFactors)); dw++)
432 : {
433 : /* for all Mobius terms */
434 460 : int iPower=B;
435 460 : int iMobius=0;
436 460 : unsigned dwIndex=0;
437 460 : unsigned dwMask=1;
438 :
439 1640 : while(dwIndex < dwFactors)
440 : {
441 1180 : if(dw&dwMask)
442 : {
443 : /* printf ("iMobius = %d iPower = %d, dwIndex = %d ", iMobius, iPower, dwIndex); */
444 590 : iMobius++;
445 590 : iPower/=factors[dwIndex];
446 : /* printf ("Then iPower = %d\n", iPower); */
447 : }
448 1180 : dwMask<<=1;
449 1180 : ++dwIndex;
450 : }
451 : /*
452 : fprintf (stderr, "Taking ");
453 : mpz_out_str(stderr, 10, org_n);
454 : fprintf (stderr, "^%d-1\n", iPower);
455 : */
456 460 : mpz_pow_ui(T, org_n, iPower);
457 460 : mpz_sub_ui(T, T, 1);
458 :
459 460 : if(iMobius&1)
460 : {
461 : /*
462 : fprintf (stderr, "Muling D=D*T ");
463 : mpz_out_str(stderr, 10, D);
464 : fprintf (stderr, "*");
465 : mpz_out_str(stderr, 10, T);
466 : fprintf (stderr, "\n");
467 : */
468 230 : mpz_mul(D, D, T);
469 : }
470 : else
471 : {
472 : /*
473 : fprintf (stderr, "Muling n=n*T ");
474 : mpz_out_str(stderr, 10, n);
475 : fprintf (stderr, "*");
476 : mpz_out_str(stderr, 10, T);
477 : fprintf (stderr, "\n");
478 : */
479 230 : mpz_mul(n, n, T);
480 : }
481 : }
482 90 : mpz_divexact(n, n, D);
483 90 : mpz_clear(T);
484 90 : mpz_clear(org_n);
485 90 : mpz_clear(D);
486 :
487 90 : return 1;
488 : }
489 :
490 70 : int aurif (mpz_t output, mpz_t n, mpz_t base, int sign) // Evaluate Aurifeullian polynomials
491 : {
492 70 : int b,k=mpz_get_ui(n);
493 : mpz_t orig_base;
494 : mpz_t C,D,l,m;
495 : // Find a proper base
496 70 : mpz_init_set(orig_base,base);
497 70 : mpz_inits(C,D,l,m,NULL);
498 480 : for(b=2;b<=11;b++)
499 : {
500 450 : mpz_set(base,orig_base);
501 450 : mpz_mul_ui(base,base,b);
502 450 : if(mpz_perfect_square_p(base)) break;
503 : }
504 70 : if(b==12) // not found
505 : {
506 30 : gmp_fprintf (stderr, "Error: base %Zd not supported for Aurifeullian factorization yet\n", orig_base);
507 30 : return 0;
508 : }
509 40 : if(k%((b==5)?b:(2*b))!=0)
510 : {
511 20 : gmp_fprintf (stderr, "Error: exponent %Zd does not make sense for base %Zd\n", n, orig_base);
512 20 : return 0;
513 : }
514 20 : k/=((b==5)?b:(2*b));
515 20 : if(k%2==0)
516 : {
517 20 : gmp_fprintf (stderr, "Error: exponent %Zd does not make sense for base %Zd\n", n, orig_base);
518 20 : return 0;
519 : }
520 0 : mpz_set(base,orig_base);
521 0 : mpz_pow_ui(m, base, k);
522 0 : mpz_mul_ui(l, m, b);
523 0 : mpz_sqrt(l, l);
524 0 : switch(b)
525 : {
526 0 : case 2:
527 : case 3:
528 0 : mpz_add_ui(C, m, 1);
529 0 : mpz_set_ui(D, 1);
530 0 : break;
531 0 : case 5:
532 : case 6:
533 0 : mpz_add_ui(C, m, 3);
534 0 : mpz_mul(C, C, m);
535 0 : mpz_add_ui(C, C, 1);
536 0 : mpz_add_ui(D, m, 1);
537 0 : break;
538 0 : case 7:
539 0 : mpz_add_ui(C, m, 1);
540 0 : mpz_pow_ui(C, C, 3);
541 0 : mpz_add_ui(D, m, 1);
542 0 : mpz_mul(D, D, m);
543 0 : mpz_add_ui(D, D, 1);
544 0 : break;
545 0 : case 10:
546 0 : mpz_add_ui(C, m, 5);
547 0 : mpz_mul(C, C, m);
548 0 : mpz_add_ui(C, C, 7);
549 0 : mpz_mul(C, C, m);
550 0 : mpz_add_ui(C, C, 5);
551 0 : mpz_mul(C, C, m);
552 0 : mpz_add_ui(C, C, 1);
553 0 : mpz_add_ui(D, m, 2);
554 0 : mpz_mul(D, D, m);
555 0 : mpz_add_ui(D, D, 2);
556 0 : mpz_mul(D, D, m);
557 0 : mpz_add_ui(D, D, 1);
558 0 : break;
559 0 : case 11:
560 0 : mpz_add_ui(C, m, 5);
561 0 : mpz_mul(C, C, m);
562 0 : mpz_sub_ui(C, C, 1);
563 0 : mpz_mul(C, C, m);
564 0 : mpz_sub_ui(C, C, 1);
565 0 : mpz_mul(C, C, m);
566 0 : mpz_add_ui(C, C, 5);
567 0 : mpz_mul(C, C, m);
568 0 : mpz_add_ui(C, C, 1);
569 0 : mpz_add_ui(D, m, 1);
570 0 : mpz_mul(D, D, m);
571 0 : mpz_sub_ui(D, D, 1);
572 0 : mpz_mul(D, D, m);
573 0 : mpz_add_ui(D, D, 1);
574 0 : mpz_mul(D, D, m);
575 0 : mpz_add_ui(D, D, 1);
576 0 : break;
577 0 : default: // not supposed to arrive here
578 0 : break;
579 : }
580 0 : mpz_set(output, C);
581 0 : (sign>0 ? mpz_addmul : mpz_submul)(output, D, l);
582 : // gmp_fprintf(stderr, "Calculated base=%Zd, exp=%Zd, C=%Zd, D=%Zd, output=%Zd\n",base,n,C,D,output);
583 0 : mpz_clears(orig_base,C,D,l,m,NULL);
584 0 : return 1;
585 : }
586 50 : int eval_PhiL (mpz_t *params, mpz_t n)
587 : {
588 : mpz_t aur;
589 : int err1,err2;
590 50 : mpz_init(aur);
591 50 : err1=aurif(aur,params[0],n,-1);
592 50 : err2=eval_Phi(params,n); // n now holds Phi(params[0],n)
593 50 : mpz_gcd(n,n,aur);
594 50 : mpz_clear(aur);
595 50 : return err1*err2;
596 : }
597 20 : int eval_PhiM (mpz_t *params, mpz_t n)
598 : {
599 : mpz_t aur;
600 : int err1,err2;
601 20 : mpz_init(aur);
602 20 : err1=aurif(aur,params[0],n,1);
603 20 : err2=eval_Phi(params,n); // n now holds Phi(params[0],n)
604 20 : mpz_gcd(n,n,aur);
605 20 : mpz_clear(aur);
606 20 : return err1*err2;
607 : }
608 :
609 0 : int eval_gcd (mpz_t *params, mpz_t n)
610 : {
611 0 : mpz_gcd(n, n, params[0]);
612 0 : return 1;
613 : }
614 :
615 20 : int eval_U (mpz_t *params, mpz_t n)
616 : /* params[0]=P, params[1]=Q */
617 : {
618 : unsigned long N;
619 : mpz_t U1,U0,org_n,D,T; /* At each step U1 holds U(k), and U0 holds U(k-1) */
620 : long k,l;
621 :
622 20 : if (mpz_cmp_si (n, 0) < 0)
623 10 : return 0;
624 10 : if (mpz_cmp_ui (n, 1) == 0)
625 : {
626 0 : mpz_set_ui (n, 1);
627 0 : return 1;
628 : }
629 10 : if (mpz_cmp_ui (n, 0) == 0)
630 : {
631 0 : mpz_set_ui (n, 0);
632 0 : return 1;
633 : }
634 10 : if (mpz_fits_ulong_p (n) == 0)
635 10 : return 0;
636 :
637 0 : N = mpz_get_ui (n);
638 0 : if (mpz_cmp_ui (params[0], 0) == 0)
639 : {
640 0 : if( N%2==0 )
641 : {
642 0 : mpz_set_ui (n, 0);
643 : }
644 : else
645 : {
646 0 : mpz_neg (params[1], params[1]);
647 0 : mpz_pow_ui (n, params[1], (N-1)/2);
648 0 : mpz_neg (params[1], params[1]);
649 : }
650 0 : return 1;
651 : }
652 :
653 :
654 0 : mpz_init_set (org_n, n);
655 0 : mpz_init_set_ui (U1, 1);
656 0 : mpz_init_set_ui (U0, 0);
657 0 : mpz_init (D);
658 0 : mpz_init (T);
659 0 : mpz_mul (D, params[0], params[0]);
660 0 : mpz_submul_ui (D, params[1], 4);
661 0 : k=1;
662 :
663 0 : for(l=mpz_sizeinbase(org_n,2)-2;l>=0;l--)
664 : {
665 0 : mpz_mul (U0, U0, U0);
666 0 : mpz_mul (U1, U1, U1);
667 0 : mpz_mul (U0, U0, params[1]);
668 0 : mpz_sub (U0, U1, U0); // U(2k-1)=U(k)^2-QU(k-1)^2
669 0 : mpz_pow_ui (T, params[1], k);
670 0 : mpz_mul (U1, U1, D);
671 0 : mpz_addmul_ui (U1, T, 2);
672 0 : mpz_addmul (U1, params[1], U0); // U(2k+1)=DU(k)^2+2Q^k+QU(2k-1)
673 0 : if (mpz_tstbit (org_n, l) )
674 : {
675 0 : k=2*k+1;
676 0 : mpz_mul (U0,U0,params[1]); // U0 is 2k, U1 is 2k+1
677 0 : mpz_add (U0,U1,U0);
678 0 : mpz_divexact (U0,U0,params[0]);
679 : }
680 : else
681 : {
682 0 : k=2*k;
683 0 : mpz_addmul (U1,U0,params[1]); // U0 is 2k-1, U1 is 2k
684 0 : mpz_divexact (U1,U1,params[0]);
685 : }
686 : /* gmp_printf("%d %Zd %Zd\n",k,U0,U1); */
687 : }
688 0 : mpz_set(n, U1);
689 :
690 0 : mpz_clear(U0);
691 0 : mpz_clear(U1);
692 0 : mpz_clear(org_n);
693 0 : mpz_clear(D);
694 0 : mpz_clear(T);
695 :
696 0 : return 1;
697 : }
698 :
699 40 : int eval_primU (mpz_t* params, mpz_t n)
700 : {
701 : int factors[200];
702 40 : unsigned dwFactors=0, dw;
703 : unsigned long N;
704 : unsigned long p;
705 : mpz_t D, T;
706 :
707 40 : if (mpz_cmp_ui (n, 0) <= 0)
708 10 : return 0;
709 30 : if (mpz_cmp_ui (n, 1) == 0)
710 : {
711 0 : mpz_set_ui (n, 1);
712 0 : return 1;
713 : }
714 :
715 : /* Ignore the special cases where P^2=0,Q or 4Q*/
716 30 : if (mpz_cmp_ui (params[0], 0) == 0)
717 : {
718 10 : return 0;
719 : }
720 20 : mpz_init(D);
721 20 : mpz_mul(D, params[0], params[0]);
722 20 : if (mpz_cmp (D, params[1]) == 0)
723 : {
724 10 : return 0;
725 : }
726 10 : mpz_submul_ui(D, params[1], 4);
727 10 : if (mpz_cmp_ui (D, 0) == 0)
728 : {
729 10 : return 0;
730 : }
731 :
732 :
733 0 : if (mpz_fits_ulong_p (n) == 0)
734 0 : return 0;
735 :
736 0 : N = mpz_get_ui (n);
737 :
738 : /* Obtain the factors of N */
739 0 : for (p = 2; p <= N; p++)
740 : {
741 0 : if (N % p == 0)
742 : {
743 : /* Add the factor one time */
744 0 : factors[dwFactors++] = p;
745 : /* but be sure to totally remove it */
746 0 : do { N /= p; } while (N % p == 0);
747 : }
748 : }
749 :
750 :
751 0 : N = mpz_get_ui (n);
752 :
753 0 : mpz_set_ui (n, 1);
754 0 : mpz_set_ui (D, 1);
755 0 : mpz_init (T);
756 :
757 0 : for(dw=0;(dw<(1U<<dwFactors)); dw++)
758 : {
759 : /* for all Mobius terms */
760 0 : int iPower=N;
761 0 : int iMobius=0;
762 0 : unsigned dwIndex=0;
763 0 : unsigned dwMask=1;
764 :
765 0 : while(dwIndex < dwFactors)
766 : {
767 0 : if(dw&dwMask)
768 : {
769 : /* printf ("iMobius = %d iPower = %d, dwIndex = %d ", iMobius, iPower, dwIndex); */
770 0 : iMobius++;
771 0 : iPower/=factors[dwIndex];
772 : /* printf ("Then iPower = %d\n", iPower); */
773 : }
774 0 : dwMask<<=1;
775 0 : ++dwIndex;
776 : }
777 : /*
778 : gmp_fprintf (stderr, "Taking U(%Zd,%Zd,%d)\n", P,Q,iPower);
779 : */
780 0 : mpz_set_ui(T,iPower);
781 0 : if(eval_U(params, T)==0)
782 : {
783 0 : return 0;
784 : }
785 :
786 0 : if(iMobius&1)
787 : {
788 : /*
789 : fprintf (stderr, "Muling D=D*T ");
790 : mpz_out_str(stderr, 10, D);
791 : fprintf (stderr, "*");
792 : mpz_out_str(stderr, 10, T);
793 : fprintf (stderr, "\n");
794 : */
795 0 : mpz_mul(D, D, T);
796 : }
797 : else
798 : {
799 : /*
800 : fprintf (stderr, "Muling n=n*T ");
801 : mpz_out_str(stderr, 10, n);
802 : fprintf (stderr, "*");
803 : mpz_out_str(stderr, 10, T);
804 : fprintf (stderr, "\n");
805 : */
806 0 : mpz_mul(n, n, T);
807 : }
808 : }
809 0 : mpz_divexact(n, n, D);
810 0 : mpz_clear(T);
811 0 : mpz_clear(D);
812 :
813 0 : return 1;
814 : }
815 :
816 : /* A simple partial-recursive decent parser */
817 3793 : int eval_2 (int bInFuncParams)
818 : {
819 : mpz_t n_stack[5];
820 : mpz_t n;
821 : mpz_t param_stack[5];
822 : int i,j;
823 : int num_base;
824 : char op_stack[5];
825 : char op;
826 : char negate;
827 : typedef int (*fptr)(mpz_t *,mpz_t);
828 3793 : const int num_of_funcs=6;
829 3793 : const char *func_names[]={"Phi","PhiL","PhiM","U","primU","gcd"};
830 3793 : const int func_num_params[]={2,2,2,3,3,2};
831 3793 : const fptr func_ptrs[]={eval_Phi,eval_PhiL,eval_PhiM,eval_U,eval_primU,eval_gcd};
832 : char *paren_position;
833 : char tentative_func_name[20];
834 : int func_id;
835 22758 : for (i=0;i<5;i++)
836 : {
837 18965 : op_stack[i]=0;
838 18965 : mpz_init(n_stack[i]);
839 18965 : mpz_init(param_stack[i]);
840 : }
841 3793 : mpz_init(n);
842 3793 : op = 0;
843 3793 : negate = 0;
844 :
845 : for (;;)
846 : {
847 5116 : if ('-'==(*expr_str))
848 : {
849 100 : expr_str++;
850 100 : negate=1;
851 : }
852 : else
853 5016 : negate=0;
854 5116 : if ('('==(*expr_str) || '['==(*expr_str) || '{'==(*expr_str)) /* Number is subexpression */
855 : {
856 139 : expr_str++;
857 139 : eval_2 (bInFuncParams);
858 129 : mpz_set(n, t);
859 : }
860 : else /* Number is decimal value */
861 : {
862 133792 : for (i=0;isdigit((int) expr_str[i]);i++)
863 : ;
864 4977 : if (!i) /* No digits found */
865 : {
866 : /* check for a valid "function" */
867 320 : paren_position=strchr(&expr_str[i], '(');
868 320 : if (NULL==paren_position)
869 : {
870 : /* No parentheses found */
871 20 : fprintf (stderr, "\nError - invalid number [%c]\n", expr_str[i]);
872 20 : exit (EXIT_FAILURE);
873 : }
874 300 : strncpy(tentative_func_name,&expr_str[i],paren_position-&expr_str[i]);
875 300 : tentative_func_name[paren_position-&expr_str[i]]='\0';
876 870 : for (func_id=0;func_id<num_of_funcs;func_id++)
877 : {
878 860 : if (!strcasecmp (tentative_func_name, func_names[func_id]))
879 290 : break;
880 : }
881 300 : if(func_id==num_of_funcs) /* No matching function name found */
882 : {
883 10 : fprintf (stderr, "Error, Unknown function %s()\n", tentative_func_name);
884 10 : exit (EXIT_FAILURE);
885 : }
886 : /* Now we can actually process existing functions */
887 290 : expr_str=paren_position+1;
888 600 : for(j=0;j<func_num_params[func_id]-1;j++)
889 : {
890 : /* eval the first parameters. NOTE we pass a 1 since we ARE in parameter mode,
891 : and this causes the ',' character to act as the end of expression */
892 360 : if(eval_2 (1) != 2)
893 : {
894 50 : fprintf (stderr, "Error, Function %s() requires %d parameters\n", func_names[func_id], func_num_params[func_id]);
895 50 : exit (EXIT_FAILURE);
896 : }
897 310 : mpz_set(param_stack[j], t);
898 : }
899 : /* Now eval the last parameter NOTE we pass a 0 since we are NOT expecting a ','
900 : character to end the expression, but are expecting a ) character to end the function */
901 240 : if (eval_2 (0))
902 : {
903 210 : mpz_set(n, t);
904 210 : if( (func_ptrs[func_id])(param_stack, n) == 0 )
905 : {
906 150 : fprintf (stderr, "\nParsing Error - Invalid "
907 : "parameter passed to the %s function\n", func_names[func_id]);
908 150 : exit (EXIT_FAILURE);
909 : }
910 : }
911 60 : goto MONADIC_SUFFIX_LOOP;
912 : }
913 : /* Now check for a hex number. If so, handle it as such */
914 4657 : num_base=10; /* assume base 10 */
915 4657 : if (i == 1 && !strncasecmp(expr_str, "0x", 2))
916 : {
917 22 : num_base = 16; /* Kick up to hex */
918 22 : expr_str += 2; /* skip the 0x string of the number */
919 3853 : for (i=0;isxdigit((int) expr_str[i]);i++)
920 : ;
921 : }
922 4657 : op=expr_str[i];
923 4657 : expr_str[i]=0;
924 4657 : mpz_set_str(n,expr_str,num_base);
925 4657 : expr_str+=i;
926 4657 : *expr_str=op;
927 : }
928 4786 : if (negate)
929 100 : mpz_neg(n,n);
930 :
931 : /* This label is needed for "normal" primorials and factorials, since they are evaluated Monadic suffix
932 : expressions. Most of this parser assumes Dyadic operators with the only exceptino being the
933 : Monadic prefix operator of "unary minus" which is handled by simply ignoring it (but remembering),
934 : and then fixing the expression up when completed. */
935 : /* This is ALSO where functions should be sent. A function should "act" like a stand alone number.
936 : We should NOT start processing, and expecting a number, but we should expect an operator first */
937 4846 : MONADIC_SUFFIX_LOOP:
938 4866 : op=*expr_str++;
939 :
940 4866 : if (0==op || ')'==op || ']'==op || '}'==op || (','==op&&bInFuncParams))
941 : {
942 3483 : eval_power (n_stack[2],n,op_stack[2]);
943 3483 : eval_product (n_stack[1],n,op_stack[1]);
944 3473 : eval_sum (n_stack[0],n,op_stack[0]);
945 3473 : mpz_set(t, n);
946 3473 : mpz_clear(n);
947 20838 : for (i=0;i<5;i++)
948 : {
949 17365 : mpz_clear(n_stack[i]);
950 17365 : mpz_clear(param_stack[i]);
951 : }
952 : /* Hurray! a valid expression (or sub-expression) was parsed! */
953 3473 : return ','==op?2:1;
954 : }
955 : else
956 : {
957 1383 : if ('^' == op)
958 : {
959 526 : eval_power (n_stack[2],n,op_stack[2]);
960 526 : mpz_set(n_stack[2], n);
961 526 : op_stack[2]='^';
962 : }
963 857 : else if ('!' == op)
964 : {
965 20 : if (!isdigit((int) *expr_str))
966 : {
967 : /* If the next char is not a digit, then this is a simple factorial, and not a "multi" factorial */
968 10 : mpz_set(n_stack[2], n);
969 10 : op_stack[2]='!';
970 10 : goto MONADIC_SUFFIX_LOOP;
971 : }
972 10 : eval_power (n_stack[2],n,op_stack[2]);
973 10 : mpz_set(n_stack[2], n);
974 10 : op_stack[2]='@';
975 : }
976 837 : else if ('#' == op)
977 : {
978 20 : if (!isdigit((int) *expr_str))
979 : {
980 : /* If the next char is not a digit, then this is a simple primorial, and not a "reduced" primorial */
981 10 : mpz_set(n_stack[2], n);
982 10 : op_stack[2]='#';
983 10 : goto MONADIC_SUFFIX_LOOP;
984 : }
985 10 : eval_power (n_stack[2],n,op_stack[2]);
986 10 : mpz_set(n_stack[2], n);
987 10 : op_stack[2]='$';
988 : }
989 : else
990 : {
991 817 : if ('.'==op || '*'==op || '/'==op || '%'==op)
992 : {
993 200 : eval_power (n_stack[2],n,op_stack[2]);
994 200 : op_stack[2]=0;
995 200 : eval_product (n_stack[1],n,op_stack[1]);
996 200 : mpz_set(n_stack[1], n);
997 200 : op_stack[1]=op;
998 : }
999 : else
1000 : {
1001 617 : if ('+'==op || '-'==op)
1002 : {
1003 577 : eval_power (n_stack[2],n,op_stack[2]);
1004 577 : op_stack[2]=0;
1005 577 : eval_product (n_stack[1],n,op_stack[1]);
1006 577 : op_stack[1]=0;
1007 577 : eval_sum (n_stack[0],n,op_stack[0]);
1008 577 : mpz_set(n_stack[0], n);
1009 577 : op_stack[0]=op;
1010 : }
1011 : else /* Error - invalid operator */
1012 : {
1013 40 : fprintf (stderr, "\nError - unknown operator: '%c'\n", op);
1014 40 : exit (EXIT_FAILURE);
1015 : }
1016 : }
1017 : }
1018 : }
1019 : }
1020 : }
1021 :
1022 2963 : void init_expr(void)
1023 : {
1024 2963 : mpz_init_set_ui(mpOne, 1);
1025 2963 : }
1026 :
1027 2474 : void free_expr(void)
1028 : {
1029 2474 : mpz_clear(mpOne);
1030 2474 : }
|