LCOV - code coverage report
Current view: top level - ecm - auxlib.c (source / functions) Hit Total Coverage
Test: unnamed Lines: 76 89 85.4 %
Date: 2022-03-21 11:19:20 Functions: 12 12 100.0 %

          Line data    Source code
       1             : /* Auxiliary routines for the ecm library.
       2             : 
       3             : Copyright 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2011 Paul Zimmermann,
       4             : Alexander Kruppa.
       5             : 
       6             : This file is part of the ECM Library.
       7             : 
       8             : The ECM Library is free software; you can redistribute it and/or modify
       9             : it under the terms of the GNU Lesser General Public License as published by
      10             : the Free Software Foundation; either version 3 of the License, or (at your
      11             : option) any later version.
      12             : 
      13             : The ECM Library is distributed in the hope that it will be useful, but
      14             : WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
      15             : or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
      16             : License for more details.
      17             : 
      18             : You should have received a copy of the GNU Lesser General Public License
      19             : along with the ECM Library; see the file COPYING.LIB.  If not, see
      20             : http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
      21             : 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */
      22             : 
      23             : /* need stdio.h and stdarg.h for gmp.h to declare gmp_vfprintf */
      24             : #include <stdio.h>
      25             : #include <stdarg.h>
      26             : #include <gmp.h>
      27             : #include "ecm-impl.h"
      28             : 
      29             : #if TIME_WITH_SYS_TIME
      30             : # include <sys/time.h>
      31             : # include <time.h>
      32             : #else
      33             : # if HAVE_SYS_TIME_H
      34             : #  include <sys/time.h>
      35             : # else
      36             : #  include <time.h>
      37             : # endif
      38             : #endif
      39             : 
      40             : #ifdef HAVE_LIMITS_H
      41             : # include <limits.h>
      42             : #else
      43             : # ifndef ULONG_MAX
      44             : #  define LONG_MAX (__GMP_ULONG_MAX / 2)
      45             : # endif
      46             : #endif
      47             : 
      48             : #ifdef HAVE_STDINT
      49             : #include <stdint.h>
      50             : #else
      51             : /* size_t is an unsigned integer so this ought to work */
      52             : #ifndef SIZE_MAX
      53             : #define SIZE_MAX (~((size_t) 0))
      54             : #endif
      55             : #endif
      56             : 
      57             : #define VERBOSE __ECM(verbose)
      58             : static int VERBOSE = OUTPUT_NORMAL;
      59             : 
      60             : void 
      61        2034 : mpz_add_si (mpz_t r, mpz_t s, long i)
      62             : {
      63        2034 :   if (i >= 0)
      64        1336 :     mpz_add_ui (r, s, (unsigned long) i);
      65             :   else
      66         698 :     mpz_sub_ui (r, s, (unsigned long) (-i));
      67        2034 : }
      68             : 
      69             : void 
      70       30752 : mpz_sub_si (mpz_t r, mpz_t s, long i)
      71             : {
      72       30752 :   if (i >= 0)
      73         400 :     mpz_sub_ui (r, s, (unsigned long) i);
      74             :   else
      75       30352 :     mpz_add_ui (r, s, (unsigned long) (-i));
      76       30752 : }
      77             : 
      78             : /* Divide RS by 3 */
      79             : void
      80    24318463 : mpz_divby3_1op (mpz_t RS)
      81             : {
      82    24318463 :   mp_size_t abssize = mpz_size (RS);
      83             :   
      84    24318463 :   if (abssize == 0)
      85           0 :     return;
      86             :   
      87    24318463 :   mpn_divexact_by3 (RS->_mp_d, RS->_mp_d, abssize);
      88             : 
      89    24318463 :   if (RS->_mp_d[abssize - 1] == 0)
      90      502202 :     RS->_mp_size -= mpz_sgn (RS);
      91             : }
      92             : 
      93             : /* Convert a double d to a size_t.
      94             :    Assumes d >= 0. If d > SIZE_MAX, returns SIZE_MAX. */
      95             : size_t
      96          58 : double_to_size (double d)
      97             : {
      98             :   ASSERT(d >= 0.0);
      99          58 :   return (d > (double) SIZE_MAX) ? SIZE_MAX : (size_t) d;
     100             : }
     101             : 
     102             : /* cputime () gives the elapsed time in milliseconds */
     103             : 
     104             : #if defined (_WIN32)
     105             : /* First case - GetProcessTimes () is the only known way of getting process
     106             :  * time (as opposed to calendar time) under mingw32 */
     107             : 
     108             : #include <windows.h>
     109             : 
     110             : long
     111             : cputime ()
     112             : {
     113             :   FILETIME lpCreationTime, lpExitTime, lpKernelTime, lpUserTime;
     114             :   ULARGE_INTEGER n;
     115             : 
     116             :   HANDLE hProcess = GetCurrentProcess();
     117             :   
     118             :   GetProcessTimes (hProcess, &lpCreationTime, &lpExitTime, &lpKernelTime,
     119             :       &lpUserTime);
     120             : 
     121             :   /* copy FILETIME to a ULARGE_INTEGER as recommended by MSDN docs */
     122             :   n.u.LowPart = lpUserTime.dwLowDateTime;
     123             :   n.u.HighPart = lpUserTime.dwHighDateTime;
     124             : 
     125             :   /* lpUserTime is in units of 100 ns. Return time in milliseconds */
     126             :   return (long) (n.QuadPart / 10000);
     127             : }
     128             : 
     129             : #elif defined (HAVE_GETRUSAGE)
     130             : /* Next case: getrusage () has higher resolution than clock () and so is
     131             :    preferred. */
     132             : 
     133             : #ifdef HAVE_SYS_TYPES_H
     134             : # include <sys/types.h>
     135             : #endif
     136             : #ifdef HAVE_SYS_RESOURCE_H
     137             : # include <sys/resource.h>
     138             : #endif
     139             : 
     140             : long
     141      463639 : cputime ()
     142             : {
     143             :   struct rusage rus;
     144             : 
     145      463639 :   getrusage (RUSAGE_SELF, &rus);
     146             :   /* This overflows a 32 bit signed int after 2147483s = 24.85 days */
     147      463639 :   return rus.ru_utime.tv_sec * 1000L + rus.ru_utime.tv_usec / 1000L;
     148             : }
     149             : 
     150             : #else
     151             : /* Resort to clock (), which on some systems may return calendar time. */
     152             : 
     153             : long
     154             : cputime ()
     155             : {
     156             :   /* Return time in milliseconds */
     157             :   return (long) (clock () * (1000. / (double) CLOCKS_PER_SEC));
     158             : }
     159             : 
     160             : #endif /* defining cputime () */
     161             : 
     162             : /* ellapsed time (in milliseconds) between st0 and st1 (values of cputime) */
     163             : long
     164      422338 : elltime (long st0, long st1)
     165             : {
     166      422338 :   return st1 - st0; /* assumes no wrap around */
     167             : }
     168             : 
     169             : /* Get real (wall-clock) time in milliseconds */
     170             : long
     171        4406 : realtime ()
     172             : {
     173             : #ifdef HAVE_GETTIMEOFDAY
     174             :   struct timeval tv;
     175        4406 :   int ret = gettimeofday (&tv, NULL);
     176        4406 :   ASSERT_ALWAYS(ret == 0); /* if HAVE_GETTIMEOFDAY, it should be functional */
     177        4406 :   return (long) tv.tv_sec * 1000L + (long) tv.tv_usec / 1000L;
     178             : #else
     179             :   return 0L;
     180             : #endif
     181             : }
     182             : 
     183             : /* Tests if loglevel gets printed with the current verbose setting */
     184             : 
     185             : int 
     186     6730584 : test_verbose (int loglevel)
     187             : {
     188     6730584 :   return (loglevel <= VERBOSE);
     189             : }
     190             : 
     191             : void 
     192        2129 : set_verbose (int v)
     193             : {
     194        2129 :   VERBOSE = v;
     195        2129 : }
     196             : 
     197             : int
     198    29953752 : outputf (int loglevel, const char *format, ...)
     199             : {
     200             :   va_list ap;
     201    29953752 :   int n = 0;
     202             :   
     203    29953752 :   va_start (ap, format);
     204             : 
     205    29953752 :   if (loglevel != OUTPUT_ERROR && loglevel <= VERBOSE)
     206             :     {
     207      499553 :       n = gmp_vfprintf (ECM_STDOUT, format, ap);
     208      499553 :       fflush (ECM_STDOUT);
     209             :     }
     210    29454199 :   else if (loglevel == OUTPUT_ERROR)
     211          72 :     n = gmp_vfprintf (ECM_STDERR, format, ap);
     212             :   
     213    29953752 :   va_end (ap);
     214             :   
     215    29953752 :   return n;
     216             : }
     217             : 
     218             : /* for P-1 and P+1 we have A = y = z = NULL */
     219             : void
     220          20 : writechkfile (char *chkfilename, int method, double p, mpmod_t modulus, 
     221             :               mpres_t A, mpres_t x, mpres_t y, mpres_t z)
     222             : {
     223             :   FILE *chkfile;
     224             :   char *methodname;
     225             :   mpz_t t;
     226             : 
     227          20 :   outputf (OUTPUT_VERBOSE, "Writing checkpoint to %s at p = %.0f\n",
     228             :            chkfilename, p);
     229             : 
     230          20 :   switch (method)
     231             :     {
     232          10 :     case ECM_ECM : methodname = "ECM"; break;
     233           5 :     case ECM_PM1 : methodname = "P-1"; break;
     234           5 :     case ECM_PP1 : methodname = "P+1"; break;
     235           0 :     default: 
     236           0 :       outputf (OUTPUT_ERROR, "writechkfile: Invalid method\n");
     237           0 :       return;
     238             :     }
     239             : 
     240          20 :   chkfile = fopen (chkfilename, "w");
     241          20 :   ASSERT_ALWAYS(chkfile != NULL);
     242             : 
     243          20 :   mpz_init (t);
     244             : 
     245          20 :   gmp_fprintf (chkfile, "METHOD=%s; B1=%.0f; N=%Zd;", 
     246          20 :                methodname, p, modulus->orig_modulus);
     247          20 :   mpres_get_z (t, x, modulus);
     248          20 :   gmp_fprintf (chkfile, " X=0x%Zx;", t);
     249          20 :   if (method == ECM_ECM)
     250             :     {
     251          10 :         if (y != NULL) /* this should mean Weierstrass form */
     252             :           {
     253             :             /* actually, we want to print (x:y:1) */
     254           0 :             mpres_get_z (t, y, modulus);
     255           0 :             gmp_fprintf (chkfile, " Y=0x%Zx;", t);
     256           0 :             fprintf (chkfile, " Z=0x1;");
     257             :           }
     258             :         else /* one day, we could have some homogeneous form to deal with */
     259             :           {
     260          10 :             mpres_get_z (t, z, modulus);
     261          10 :             gmp_fprintf (chkfile, " Z=0x%Zx;", t);
     262             :           }
     263          10 :       mpres_get_z (t, A, modulus);
     264          10 :       gmp_fprintf (chkfile, " A=0x%Zx;", t);
     265             :     }
     266          20 :   fprintf (chkfile, "\n");
     267          20 :   mpz_clear (t);
     268          20 :   fflush (chkfile);
     269          20 :   fclose (chkfile);
     270             : }
     271             : 
     272             : #if 0 /* currently unused (only used in listz_handle.c, currently inactive) */
     273             : int 
     274             : aux_fseek64(FILE *f, const int64_t offset, const int whence)
     275             : {
     276             : #ifdef HAVE__FSEEKI64
     277             :   return _fseeki64(f, offset, whence);
     278             : #endif
     279             : #if LONG_MAX == INT64_MAX
     280             :   return fseek (f, (long) offset, whence);
     281             : #endif
     282             :   ASSERT_ALWAYS (offset <= LONG_MAX);
     283             :   return fseek (f, (long) offset, whence);
     284             : }
     285             : #endif
     286             : 
     287             : int
     288    21151796 : ecm_tstbit (mpz_srcptr u, ecm_uint bit_index)
     289             : {
     290    21151796 :   mp_srcptr     u_ptr      = PTR(u);
     291    21151796 :   ecm_int       size       = SIZ(u);
     292    21151796 :   ecm_uint      abs_size   = ABS(size);
     293    21151796 :   ecm_uint      limb_index = bit_index / GMP_NUMB_BITS;
     294    21151796 :   mp_srcptr     p          = u_ptr + limb_index;
     295             :   mp_limb_t     limb;
     296             : 
     297    21151796 :   if (limb_index >= abs_size)
     298      216328 :     return (size < 0);
     299             : 
     300    20935468 :   limb = *p;
     301    20935468 :   if (size < 0)
     302             :     {
     303           0 :       limb = -limb;     /* twos complement */
     304             : 
     305           0 :       while (p != u_ptr)
     306             :         {
     307           0 :           p--;
     308           0 :           if (*p != 0)
     309             :             {
     310           0 :               limb--;   /* make it a ones complement instead */
     311           0 :               break;
     312             :             }
     313             :         }
     314             :     }
     315             : 
     316    20935468 :   return (limb >> (bit_index % GMP_NUMB_BITS)) & 1;
     317             : }

Generated by: LCOV version 1.14