github.com/aergoio/aergo@v1.3.1/libtool/src/gmp-6.1.2/nextprime.c (about)

     1  /* gmp_nextprime -- generate small primes reasonably efficiently for internal
     2     GMP needs.
     3  
     4     Contributed to the GNU project by Torbjorn Granlund.  Miscellaneous
     5     improvements by Martin Boij.
     6  
     7     THE FUNCTIONS IN THIS FILE ARE INTERNAL WITH MUTABLE INTERFACES.  IT IS ONLY
     8     SAFE TO REACH THEM THROUGH DOCUMENTED INTERFACES.  IN FACT, IT IS ALMOST
     9     GUARANTEED THAT THEY WILL CHANGE OR DISAPPEAR IN A FUTURE GNU MP RELEASE.
    10  
    11  Copyright 2009 Free Software Foundation, Inc.
    12  
    13  This file is part of the GNU MP Library.
    14  
    15  The GNU MP Library is free software; you can redistribute it and/or modify
    16  it under the terms of either:
    17  
    18    * the GNU Lesser General Public License as published by the Free
    19      Software Foundation; either version 3 of the License, or (at your
    20      option) any later version.
    21  
    22  or
    23  
    24    * the GNU General Public License as published by the Free Software
    25      Foundation; either version 2 of the License, or (at your option) any
    26      later version.
    27  
    28  or both in parallel, as here.
    29  
    30  The GNU MP Library is distributed in the hope that it will be useful, but
    31  WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
    32  or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
    33  for more details.
    34  
    35  You should have received copies of the GNU General Public License and the
    36  GNU Lesser General Public License along with the GNU MP Library.  If not,
    37  see https://www.gnu.org/licenses/.  */
    38  
    39  /*
    40    Optimisation ideas:
    41  
    42    1. Unroll the sieving loops.  Should reach 1 write/cycle.  That would be a 2x
    43       improvement.
    44  
    45    2. Separate sieving with primes p < SIEVESIZE and p >= SIEVESIZE.  The latter
    46       will need at most one write, and thus not need any inner loop.
    47  
    48    3. For primes p >= SIEVESIZE, i.e., typically the majority of primes, we
    49       perform more than one division per sieving write.  That might dominate the
    50       entire run time for the nextprime function.  A incrementally initialised
    51       remainder table of Pi(65536) = 6542 16-bit entries could replace that
    52       division.
    53  */
    54  
    55  #include "gmp.h"
    56  #include "gmp-impl.h"
    57  #include <string.h>		/* for memset */
    58  
    59  
    60  unsigned long int
    61  gmp_nextprime (gmp_primesieve_t *ps)
    62  {
    63    unsigned long p, d, pi;
    64    unsigned char *sp;
    65    static unsigned char addtab[] =
    66      { 2,4,2,4,6,2,6,4,2,4,6,6,2,6,4,2,6,4,6,8,4,2,4,2,4,8,6,4,6,2,4,6,2,6,6,4,
    67        2,4,6,2,6,4,2,4,2,10,2,10 };
    68    unsigned char *addp = addtab;
    69    unsigned long ai;
    70  
    71    /* Look for already sieved primes.  A sentinel at the end of the sieving
    72       area allows us to use a very simple loop here.  */
    73    d = ps->d;
    74    sp = ps->s + d;
    75    while (*sp != 0)
    76      sp++;
    77    if (sp != ps->s + SIEVESIZE)
    78      {
    79        d = sp - ps->s;
    80        ps->d = d + 1;
    81        return ps->s0 + 2 * d;
    82      }
    83  
    84    /* Handle the number 2 separately.  */
    85    if (ps->s0 < 3)
    86      {
    87        ps->s0 = 3 - 2 * SIEVESIZE; /* Tricky */
    88        return 2;
    89      }
    90  
    91    /* Exhausted computed primes.  Resieve, then call ourselves recursively.  */
    92  
    93  #if 0
    94    for (sp = ps->s; sp < ps->s + SIEVESIZE; sp++)
    95      *sp = 0;
    96  #else
    97    memset (ps->s, 0, SIEVESIZE);
    98  #endif
    99  
   100    ps->s0 += 2 * SIEVESIZE;
   101  
   102    /* Update sqrt_s0 as needed.  */
   103    while ((ps->sqrt_s0 + 1) * (ps->sqrt_s0 + 1) <= ps->s0 + 2 * SIEVESIZE - 1)
   104      ps->sqrt_s0++;
   105  
   106    pi = ((ps->s0 + 3) / 2) % 3;
   107    if (pi > 0)
   108      pi = 3 - pi;
   109    if (ps->s0 + 2 * pi <= 3)
   110      pi += 3;
   111    sp = ps->s + pi;
   112    while (sp < ps->s + SIEVESIZE)
   113      {
   114        *sp = 1, sp += 3;
   115      }
   116  
   117    pi = ((ps->s0 + 5) / 2) % 5;
   118    if (pi > 0)
   119      pi = 5 - pi;
   120    if (ps->s0 + 2 * pi <= 5)
   121      pi += 5;
   122    sp = ps->s + pi;
   123    while (sp < ps->s + SIEVESIZE)
   124      {
   125        *sp = 1, sp += 5;
   126      }
   127  
   128    pi = ((ps->s0 + 7) / 2) % 7;
   129    if (pi > 0)
   130      pi = 7 - pi;
   131    if (ps->s0 + 2 * pi <= 7)
   132      pi += 7;
   133    sp = ps->s + pi;
   134    while (sp < ps->s + SIEVESIZE)
   135      {
   136        *sp = 1, sp += 7;
   137      }
   138  
   139    p = 11;
   140    ai = 0;
   141    while (p <= ps->sqrt_s0)
   142      {
   143        pi = ((ps->s0 + p) / 2) % p;
   144        if (pi > 0)
   145  	pi = p - pi;
   146        if (ps->s0 + 2 * pi <= p)
   147  	  pi += p;
   148        sp = ps->s + pi;
   149        while (sp < ps->s + SIEVESIZE)
   150  	{
   151  	  *sp = 1, sp += p;
   152  	}
   153        p += addp[ai];
   154        ai = (ai + 1) % 48;
   155      }
   156    ps->d = 0;
   157    return gmp_nextprime (ps);
   158  }
   159  
   160  void
   161  gmp_init_primesieve (gmp_primesieve_t *ps)
   162  {
   163    ps->s0 = 0;
   164    ps->sqrt_s0 = 0;
   165    ps->d = SIEVESIZE;
   166    ps->s[SIEVESIZE] = 0;		/* sentinel */
   167  }