/* factor.c -- factorization of very smooth numbers.
 *
 * Copyright (C) 2012, 2013, 2021 INRIA
 *
 * This file is part of CMH.
 *
 * CMH 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 3 of the License, or (at your
 * option) any later version.
 *
 * CMH 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 this program. If not, see http://www.gnu.org/licenses/ .
 */

#include <stdlib.h>
#include <stdio.h>
#include <limits.h>
#include <string.h>
#include <pari/pari.h>
#include "parsepari.h"
#include "factor.h"

static factor_matrix factor_matrix_init (int size)
   /* Allocate a matrix f to hold size factors, where f [i, 0] is the i-th
      prime factor and f [i, 1] its multiplicity; the last entry is
      indicated by a prime factor 0 */
{
   factor_matrix f;

   f = (factor_matrix) malloc ((size + 1) * 2 * sizeof (unsigned long));

   return f;
}


factor_matrix factor_matrix_dup(factor_matrix f)
{
    int size;
    for (size = 0; f [size][0] ; size++);
    factor_matrix g = factor_matrix_init (size);
    memcpy(g, f, (size + 1) * 2 * sizeof (unsigned long));
    return g;
}


void factor_matrix_free (factor_matrix f)
{
    free (f);
}


factor_matrix factor_trialdiv (mpz_srcptr n, unsigned long int bound)
   /* factors n by trial division with numbers <= bound, or indefinitely
         if bound is 0.
      The return value is a newly allocated factor list if the number is
      smooth, NULL otherwise.  */
{
   pari_sp av;
   factor_matrix f;
   GEN n_gen, fact;
   mpz_t tmp;
   unsigned long int last;
   int i, m, ok;

   if (mpz_cmp_ui (n, 1) == 0) {
      f = factor_matrix_init (0);
      f [0][0] = 0;
      return f;
   }

   if (bound == 0)
      bound = ULONG_MAX;

   pari_init (500000, 0);
   av = avma;

   n_gen = mpz_get_GEN (n);
   fact = boundfact (n_gen, bound);
   m = lg (gel (fact, 1)) - 1;
   last = itou_or_0 (gel (gel (fact, 1), m));
   if (last > bound || last == 0)
      ok = 0;
   else {
      ok = 1;
      mpz_init (tmp);
      f = factor_matrix_init (m);
      for (i = 0; i < m; i++) {
         mpz_set_GEN (tmp, gel (gel (fact, 1), i+1));
         f [i][0] = mpz_get_ui (tmp);
         mpz_set_GEN (tmp, gel (gel (fact, 2), i+1));
         f [i][1] = mpz_get_ui (tmp);
      }
      f [m][0] = 0;
      mpz_clear (tmp);
   }

   avma = av;
   pari_close ();

   return (ok ? f : NULL);
}


void factor_matrix_asprintf (char ** pres, factor_matrix f)
{
   size_t alloc = 10;
   size_t size = 0;
   char * res = malloc (alloc);
   int d, i;
   
#define PUSH(fmt__, arg__) do {                                         \
    for( ; (d = snprintf(res + size, alloc-size, fmt__, arg__)) >= (int) (alloc-size) ; ) {     \
        alloc = size+d+1 + alloc / 4;                                   \
        res = realloc(res, alloc);                                      \
    }                                                                   \
    size+=d;                                                            \
} while (0)

   if (f[0][0] == 0) {
       snprintf(res, alloc, "1");
       *pres = res;
       return;
   }
   for (i = 0; f [i][0] != 0; i++) {
      if (i > 0)
         PUSH ("%c", '*');
      PUSH ("%lu", f [i][0]);
      if (f [i][1] > 1)
         PUSH ("^%lu", f [i][1]);
   }
#undef PUSH

   *pres = res;
}


int factor_matrix_fprintf (FILE *file, factor_matrix f)
{
   int ok, i, size;

   for (size = 0; f [size][0] != 0; size++);
   ok = fprintf (file, "%i", size);
   for (i = 0; i < size; i++)
      ok += fprintf (file, " %lu %lu", f [i][0], f [i][1]);

   return (ok > 0);
}


static void factor_matrix_fold_inner(mpz_ptr z, factor_matrix f, int i0, int i1)
{
    if (i1 - i0 == 1) {
        mpz_ui_pow_ui(z, f[i0][0], f[i0][1]);
    } else {
        int j = (i0 + i1) / 2;
        mpz_t l,r;
        mpz_init(l);
        mpz_init(r);
        factor_matrix_fold_inner(l, f, i0, j);
        factor_matrix_fold_inner(r, f, j, i1);
        mpz_mul(z, l, r);
        mpz_clear(l);
        mpz_clear(r);
    }
}

void factor_matrix_fold(mpz_ptr z, factor_matrix f)
{
   int size;
   for (size = 0; f [size][0] != 0; size++);
   if (size == 0) {
       mpz_set_ui(z, 1);
   }
   else
      factor_matrix_fold_inner(z, f, 0, size);
}


void mpz_remove_squares(mpz_ptr D)
{
   factor_matrix f;
   int i;

   f = factor_trialdiv (D, 0);
   for (i = 0; f [i][0] != 0; i++)
      f [i][1] %= 2;

   factor_matrix_fold (D, f);

   factor_matrix_free (f);
}