Line data    Source code

       1             : /* Auxiliary arithmetic routines on unsigned long ints for the ecm library.
       2             : 
       3             : Copyright 2001, 2002, 2003, 2004, 2005, 2007, 2008 Paul Zimmermann and 
       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             : #include "config.h"
      24             : #include "ecm-impl.h"
      25             : 
      26             : /* Returns the gcd of a and b */
      27             : unsigned long
      28     6135223 : gcd (unsigned long a, unsigned long b)
      29             : {
      30             :   unsigned long t;
      31             : 
      32    34735960 :   while (b != 0UL)
      33             :     {
      34    28600737 :       t = a % b;
      35    28600737 :       a = b;
      36    28600737 :       b = t;
      37             :     }
      38             : 
      39     6135223 :   return a;
      40             : }
      41             : 
      42             : 
      43             : /* returns Euler's totient phi function */
      44             : unsigned long
      45      365796 : eulerphi (unsigned long n)
      46             : {
      47      365796 :   unsigned long phi = 1UL, p;
      48             : 
      49     3630400 :   for (p = 2UL; p * p <= n; p += 2UL)
      50             :     {
      51     3264604 :       if (n % p == 0UL)
      52             :         {
      53     2009764 :           phi *= p - 1UL;
      54     2009764 :           n /= p;
      55     2372894 :           while (n % p == 0UL)
      56             :             {
      57      363130 :               phi *= p;
      58      363130 :               n /= p;
      59             :             }
      60             :         }
      61             : 
      62     3264604 :       if (p == 2UL)
      63      363882 :         p--;
      64             :     }
      65             : 
      66             :   /* now n is prime or 1 */
      67             : 
      68      365796 :   return (n == 1UL) ? phi : phi * (n - 1UL);
      69             : }
      70             : 
      71             : 
      72             : /* returns ceil(log(n)/log(2)) */
      73             : unsigned int
      74       14459 : ceil_log2 (unsigned long n)
      75             : {
      76       14459 :   unsigned int k = 0;
      77             : 
      78             :   ASSERT (n > 0UL);
      79             : 
      80       14459 :   n--;
      81       97241 :   while (n)
      82             :     {
      83       82782 :       k++;
      84       82782 :       n >>= 1;
      85             :     }
      86             : 
      87       14459 :   return k;
      88             : }
      89             : 
      90             : /* Returns the smallest prime factor of N. If N == 1, return 1. */
      91             : unsigned long
      92     3495507 : find_factor (const unsigned long N)
      93             : {
      94             :   unsigned long i;
      95             : 
      96     3495507 :   ASSERT_ALWAYS (N != 0UL);
      97             : 
      98     3495507 :   if (N == 1UL)
      99           0 :     return 1UL;
     100             : 
     101     3495507 :   if (N % 2UL == 0UL)
     102           0 :     return 2UL;
     103             : 
     104     7369516 :   for (i = 3UL; i*i <= N; i += 2UL)
     105     6435376 :     if (N % i == 0UL)
     106     2561367 :       return i;
     107             : 
     108      934140 :   return N;
     109             : }