LCOV - code coverage report
Current view: top level - ecm - eval.c (source / functions) Hit Total Coverage
Test: unnamed Lines: 344 500 68.8 %
Date: 2022-03-21 11:19:20 Functions: 14 15 93.3 %

          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 : }

Generated by: LCOV version 1.14