<commit_before>/*=============================================================================

    This file is part of FLINT.

    FLINT is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 2 of the License, or
    (at your option) any later version.

    FLINT is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with FLINT; if not, write to the Free Software
    Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301 USA

=============================================================================*/
/******************************************************************************

   Copyright (C) 2006, 2011 <NAME>

******************************************************************************/

#define ulong ulongxx /* interferes with system includes */
#include <stdlib.h>
#include <stdio.h>
#undef ulong
#define ulong mp_limb_t

#include <gmp.h>
#include "flint.h"
#include "ulong_extras.h"
#include "qsieve.h"
#include "fmpz.h"

#include <time.h>

/* taken from FLINT2 */

void qsieve_do_sieving(qs_t qs_inf, unsigned char * sieve)
{
   slong num_primes = qs_inf->num_primes;
   mp_limb_t * soln1 = qs_inf->soln1;
   mp_limb_t * soln2 = qs_inf->soln2;
   prime_t * factor_base = qs_inf->factor_base;
   mp_limb_t p;
   char * end = sieve + qs_inf->sieve_size;
   register char * pos1;
   register char * pos2;
   register char * bound;
   slong size;
   slong diff;
   slong pind;

   memset(sieve, 0, qs_inf->sieve_size + sizeof(ulong));
   *end = (char) 255;

   for (pind = qs_inf->small_primes; pind < num_primes; pind++)
   {
      if (soln2[pind] == 0) continue; /* don't sieve with A factors */

      p = factor_base[pind].p;
      size = factor_base[pind].size;
      pos1 = sieve + soln1[pind];
      pos2 = sieve + soln2[pind];
      diff = pos2 - pos1;
      bound = end - 2*p;

      while (bound - pos1 > 0)
      {
         (*pos1) += size, (*(pos1 + diff)) += size, pos1 += p;
         (*pos1) += size, (*(pos1 + diff)) += size, pos1 += p;
      }

      while ((end - pos1 > 0) && (end - pos1 - diff > 0))
      {
         (*pos1) += size, (*(pos1 + diff)) += size, pos1 += p;
      }

      pos2 = pos1 + diff;

      if (end - pos2 > 0)
      {
         (*pos2) += size;
      }

      if (end - pos1 > 0)
      {
         (*pos1) += size;
      }
   }
}

void qsieve_do_sieving2(qs_t qs_inf, unsigned char * sieve)
{
    slong b, d1, d2, Bp, B, i, j;
    slong pind, size, r0;
    slong num_primes = qs_inf->num_primes;
    mp_limb_t x, p;
    mp_limb_t * soln1 = qs_inf->soln1;
    mp_limb_t * soln2 = qs_inf->soln2;
    mp_limb_t * xr1 = qs_inf->xr1;
    mp_limb_t * xr2 = qs_inf->xr2;
    prime_t * factor_base = qs_inf->factor_base;

    memset(sieve, 0, qs_inf->sieve_size + sizeof(ulong));

    for (i = 0; i < num_primes; i++)
    {
        xr1[i] = soln1[i];
        xr2[i] = soln2[i] - xr1[i];
    }

    for (b = 1; b <= qs_inf->sieve_size / BLOCK_SIZE; b++)
    {
        B = b * BLOCK_SIZE;

        for (pind = qs_inf->small_primes; pind < qs_inf->num_primes; pind++)
        {
            if (soln2[pind] == 0)
                continue;

            p = factor_base[pind].p;

            if (p > BLOCK_SIZE)
                break;

            size = factor_base[pind].size;
            d1 = xr2[pind];
            d2 = p - d1;
            Bp = B - d1;

            for (x = xr1[pind]; x <= Bp; )
            {
                sieve[x] += size, x += d1;
                sieve[x] += size, x += d2;
            }

            if (x <= B)
            {
                sieve[x] += size, x += d1;
                xr2[pind] = d2;
            }
            else {xr2[pind] = d1;}

            xr1[pind] = x;
        }
    }

    r0 = pind;

    for (b = 1; b <= qs_inf->sieve_size / BLOCK_SIZE; b++)
    {
       B = b * BLOCK_SIZE;

       for (pind = r0; pind < num_primes; pind++)
        {
            p = factor_base[pind].p;

            if (soln2[pind] == 0)
                continue;

            size = factor_base[pind].size;

            if ((x = xr1[pind]) <= B)
            {
                sieve[x] += size;
                x += xr2[pind];

                if (x <= B)
                {
                    sieve[x] += size;
                    x += p - xr2[pind];
                }
                else {xr2[pind] = p - xr2[pind];}

                xr1[pind] = x;
            }
            else {xr1[pind] = x;}
        }
    }
}


/* check position 'i' in sieve array for smoothness */

slong qsieve_evaluate_candidate(qs_t qs_inf, slong i, unsigned char * sieve)
{
   slong bits, exp, extra_bits;
   mp_limb_t modp, prime;
   slong num_primes = qs_inf->num_primes;
   prime_t * factor_base = qs_inf->factor_base;
   fac_t * factor = qs_inf->factor;
   mp_limb_t * soln1 = qs_inf->soln1;
   mp_limb_t * soln2 = qs_inf->soln2;
   slong * small = qs_inf->small;
   mp_limb_t pinv;<commit_msg>optimizing block sieve<commit_after>/*=============================================================================

    This file is part of FLINT.

    FLINT is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 2 of the License, or
    (at your option) any later version.

    FLINT is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with FLINT; if not, write to the Free Software
    Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301 USA

=============================================================================*/
/******************************************************************************

   Copyright (C) 2006, 2011 <NAME>

******************************************************************************/

#define ulong ulongxx /* interferes with system includes */
#include <stdlib.h>
#include <stdio.h>
#undef ulong
#define ulong mp_limb_t

#include <gmp.h>
#include "flint.h"
#include "ulong_extras.h"
#include "qsieve.h"
#include "fmpz.h"

#include <time.h>

/* taken from FLINT2 */

void qsieve_do_sieving(qs_t qs_inf, unsigned char * sieve)
{
   slong num_primes = qs_inf->num_primes;
   mp_limb_t * soln1 = qs_inf->soln1;
   mp_limb_t * soln2 = qs_inf->soln2;
   prime_t * factor_base = qs_inf->factor_base;
   mp_limb_t p;
   char * end = sieve + qs_inf->sieve_size;
   register char * pos1;
   register char * pos2;
   register char * bound;
   slong size;
   slong diff;
   slong pind;

   memset(sieve, 0, qs_inf->sieve_size + sizeof(ulong));
   *end = (char) 255;

   for (pind = qs_inf->small_primes; pind < num_primes; pind++)
   {
      if (soln2[pind] == 0) continue; /* don't sieve with A factors */

      p = factor_base[pind].p;
      size = factor_base[pind].size;
      pos1 = sieve + soln1[pind];
      pos2 = sieve + soln2[pind];
      diff = pos2 - pos1;
      bound = end - 2*p;

      while (bound - pos1 > 0)
      {
         (*pos1) += size, (*(pos1 + diff)) += size, pos1 += p;
         (*pos1) += size, (*(pos1 + diff)) += size, pos1 += p;
      }

      while ((end - pos1 > 0) && (end - pos1 - diff > 0))
      {
         (*pos1) += size, (*(pos1 + diff)) += size, pos1 += p;
      }

      pos2 = pos1 + diff;

      if (end - pos2 > 0)
      {
         (*pos2) += size;
      }

      if (end - pos1 > 0)
      {
         (*pos1) += size;
      }
   }
}

void qsieve_do_sieving2(qs_t qs_inf, unsigned char * sieve)
{
    slong b, d1, d2, i, j;
    slong pind, size, r0;
    slong num_primes = qs_inf->num_primes;
    mp_limb_t p;
    mp_limb_t * soln1 = qs_inf->soln1;
    mp_limb_t * soln2 = qs_inf->soln2;
    mp_limb_t * xr1 = qs_inf->xr1;
    mp_limb_t * xr2 = qs_inf->xr2;
    prime_t * factor_base = qs_inf->factor_base;

    unsigned char * Bp;
    unsigned char * B;
    register unsigned char * x;

    memset(sieve, 0, qs_inf->sieve_size + sizeof(ulong));

    for (i = 0; i < num_primes; i++)
    {
        xr1[i] = soln1[i];
        xr2[i] = soln2[i] - xr1[i];
    }

    for (i = 0; i < num_primes; i++)
    {
        if (factor_base[i].p > BLOCK_SIZE)
            break;
    }

    r0 = i;

    for (b = 1; b <= qs_inf->sieve_size / BLOCK_SIZE; b++)
    {
        B = sieve + b * BLOCK_SIZE;

        for (pind = qs_inf->small_primes; pind < r0; pind++)
        {
            if (soln2[pind] == 0)
                continue;

            p = factor_base[pind].p;
            size = factor_base[pind].size;
            d1 = xr2[pind];
            d2 = p - d1;
            Bp = (B - xr2[pind]);

            for (x = sieve + xr1[pind]; x <= Bp; )
            {
                (*x) += size, x += d1;
                (*x) += size, x += d2;
            }

            if (x <= B)
            {
                (*x) += size, x += d1;
                xr2[pind] = d2;
            }
            else {xr2[pind] = d1;}

            xr1[pind] = (x - sieve);
        }
    }

    for (b = 1; b <= qs_inf->sieve_size / BLOCK_SIZE; b++)
    {
       B = sieve + b * BLOCK_SIZE;

       for (pind = r0; pind < num_primes; pind++)
        {
            p = factor_base[pind].p;

            if (soln2[pind] == 0)
                continue;

            size = factor_base[pind].size;

            if ((x = sieve + xr1[pind]) <= B)
            {
                (*x) += size;
                 x += xr2[pind];

                if (x <= B)
                {
                    (*x) += size;
                    x += p - xr2[pind];
                }
                else {xr2[pind] = p - xr2[pind];}

                xr1[pind] = (x - sieve);
            }
            else {xr1[pind] = (x - sieve);}
        }
    }
}

/* check position 'i' in sieve array for smoothness */

slong qsieve_evaluate_candidate(qs_t qs_inf, slong i, unsigned char * sieve)
{
   slong bits, exp, extra_bits;
   mp_limb_t modp, prime;
   slong num_primes = qs_inf->num_primes;
   prime_t * factor_base = qs_inf->factor_base;
   fac_t * factor = qs_inf->factor;
   mp_limb_t * soln1 = qs_inf->soln1;
   mp_limb_t * soln2 = qs_inf->soln2;
   slong * small = qs_inf->small;
   mp_limb_t pinv;
{'fraction_non_alphanumeric': 0.18072740907386578, 'fraction_numerical': 0.021747281589801274, 'mean_word_length': 1.971866295264624, 'pattern_counts': {'":': 0, '<': 34, '<?xml version=': 0, '>': 72, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '6240759', 'n_tokens_mistral': 4036, 'n_tokens_neox': 3826, 'n_words': 999}