/*
Copyright (C) 2015 FLINT authors
This file is part of FLINT.
FLINT is free software: you can redistribute it and/or modify it under
the terms of the GNU Lesser General Public License (LGPL) as published
by the Free Software Foundation; either version 2.1 of the License, or
(at your option) any later version. See .
*/
#include
#include "flint.h"
#include "fmpz.h"
#include "fmpz_vec.h"
#define fr_node_mref(x) (&(x)->m)
typedef struct fr_node_struct
{
fmpz m;
ulong e;
struct fr_node_struct * next;
} fr_node_struct;
typedef fr_node_struct * fr_node_ptr;
/* functions related to factor refinement nodes */
void fr_node_init(fr_node_ptr x);
void fr_node_init_fmpz_ui(fr_node_ptr x, const fmpz_t m, ulong e);
void fr_node_clear(fr_node_ptr x);
int fr_node_is_one(fr_node_ptr x);
void fr_node_set_fmpz_ui(fr_node_ptr x, const fmpz_t m, ulong e);
void fr_node_get_fmpz_ui(fmpz_t m, ulong * e, fr_node_ptr x);
int fr_node_base_cmp(const fr_node_ptr x, const fr_node_ptr y);
int fr_node_base_pcmp(const void * a, const void * b);
/* functions related to lists of factor refinement nodes */
slong fr_node_list_length(fr_node_ptr x);
void fr_node_list_pop_front(fr_node_ptr *phead, fr_node_ptr *ptail);
void fr_node_list_concat(fr_node_ptr *phead, fr_node_ptr *ptail,
fr_node_ptr rhead, fr_node_ptr rtail);
void fr_node_list_clear(fr_node_ptr head);
void fr_node_list_print(fr_node_ptr head);
/* fmpz_factor_t convenience function */
int _fmpz_factor_sgn(const fmpz_factor_t f);
/* functions related to the actual algorithms of interest */
void pair_refine_unreduced(fr_node_ptr *phead,
fmpz_t m1, ulong e1, fmpz_t m2, ulong e2);
void remove_ones(fr_node_ptr *phead, fr_node_ptr *ptail, fr_node_ptr ohead);
void pair_refine(fr_node_ptr *phead, fr_node_ptr *ptail,
fmpz_t m1, ulong e1, fmpz_t m2, ulong e2);
void augment_refinement(fr_node_ptr *phead, fr_node_ptr *ptail,
const fmpz_t m_jp1, ulong e_jp1,
fr_node_ptr L_j, fr_node_ptr L_j_tail);
void
fr_node_init(fr_node_ptr x)
{
fmpz_init(fr_node_mref(x));
x->e = 0;
x->next = NULL;
}
void
fr_node_init_fmpz_ui(fr_node_ptr x, const fmpz_t m, ulong e)
{
fmpz_init_set(fr_node_mref(x), m);
x->e = e;
x->next = NULL;
}
void
fr_node_clear(fr_node_ptr x)
{
fmpz_clear(fr_node_mref(x));
x->e = 0;
x->next = NULL;
}
int
fr_node_is_one(fr_node_ptr x)
{
/* follow the fmpz_pow_ui convention 0^0 = 1 */
return (x->e == WORD(0) || fmpz_is_one(fr_node_mref(x)));
}
void
fr_node_set_fmpz_ui(fr_node_ptr x, const fmpz_t m, ulong e)
{
fmpz_set(fr_node_mref(x), m);
x->e = e;
}
void
fr_node_get_fmpz_ui(fmpz_t m, ulong * e, fr_node_ptr x)
{
fmpz_set(m, fr_node_mref(x));
*e = x->e;
}
int
fr_node_base_cmp(const fr_node_ptr x, const fr_node_ptr y)
{
return fmpz_cmp(fr_node_mref(x), fr_node_mref(y));
}
int
fr_node_base_pcmp(const void * a, const void * b)
{
/* intended for qsorting an array of node pointers */
fr_node_ptr x = *( (fr_node_ptr *) a);
fr_node_ptr y = *( (fr_node_ptr *) b);
return fr_node_base_cmp(x, y);
}
slong
fr_node_list_length(fr_node_ptr x)
{
slong count;
count = 0;
while (x)
{
count++;
x = x->next;
}
return count;
}
void
fr_node_list_pop_front(fr_node_ptr *phead, fr_node_ptr *ptail)
{
fr_node_ptr tmp;
if (phead == ptail)
{
flint_printf("aliasing issue...\n");
abort();
}
if (*phead)
{
if (*phead == *ptail)
{
*ptail = NULL;
}
tmp = (*phead)->next;
fr_node_clear(*phead);
flint_free(*phead);
*phead = tmp;
}
}
void
fr_node_list_concat(fr_node_ptr *phead, fr_node_ptr *ptail,
fr_node_ptr rhead, fr_node_ptr rtail)
{
/* if the second list is empty then do nothing */
if (!rhead)
{
return;
}
/* if the first list is empty then use the second list */
if (!(*phead))
{
*phead = rhead;
*ptail = rtail;
return;
}
/* neither list is empty so connect the tail to the head */
(*ptail)->next = rhead;
*ptail = rtail;
}
void
fr_node_list_clear(fr_node_ptr head)
{
fr_node_ptr curr, next;
curr = head;
while (curr)
{
next = curr->next;
fr_node_clear(curr);
flint_free(curr);
curr = next;
}
}
void
fr_node_list_print(fr_node_ptr head)
{
fr_node_ptr curr;
curr = head;
while (curr)
{
fmpz_print(fr_node_mref(curr));
flint_printf("^%wu ", curr->e);
curr = curr->next;
}
flint_printf("\n");
}
void
pair_refine_unreduced(fr_node_ptr *phead,
fmpz_t m1, ulong e1,
fmpz_t m2, ulong e2)
{
fr_node_ptr head, tail, curr, next, neo;
fmpz_t d;
int boring;
if (fmpz_is_one(m1) && fmpz_is_one(m2))
{
*phead = NULL;
return;
}
fmpz_init(d);
head = flint_malloc(sizeof(fr_node_struct));
fr_node_init_fmpz_ui(head, m1, e1);
tail = flint_malloc(sizeof(fr_node_struct));
fr_node_init_fmpz_ui(tail, m2, e2);
head->next = tail;
boring = 0;
while (!boring)
{
curr = head;
next = curr->next;
boring = 1;
while (next)
{
if (!fr_node_is_one(curr) && !fr_node_is_one(next))
{
fmpz_gcd(d, fr_node_mref(curr), fr_node_mref(next));
fmpz_divexact(fr_node_mref(curr), fr_node_mref(curr), d);
fmpz_divexact(fr_node_mref(next), fr_node_mref(next), d);
neo = flint_malloc(sizeof(fr_node_struct));
fr_node_init(neo);
fmpz_set(fr_node_mref(neo), d);
neo->e = curr->e + next->e;
curr->next = neo;
neo->next = next;
next = neo;
boring = 0;
}
else
{
curr = next;
next = next->next;
}
}
}
fmpz_clear(d);
*phead = head;
}
void
remove_ones(fr_node_ptr *phead, fr_node_ptr *ptail, fr_node_ptr ohead)
{
fr_node_ptr head, curr, ocurr, onext;
if (!ohead)
{
*phead = NULL;
*ptail = NULL;
return;
}
head = NULL;
curr = NULL;
ocurr = ohead;
while (ocurr)
{
onext = ocurr->next;
if (fr_node_is_one(ocurr))
{
fr_node_clear(ocurr);
flint_free(ocurr);
}
else
{
if (!head)
{
head = ocurr;
curr = ocurr;
}
else
{
curr->next = ocurr;
curr = curr->next;
}
}
ocurr = onext;
}
curr->next = NULL;
*phead = head;
*ptail = curr;
}
void
pair_refine(fr_node_ptr *phead, fr_node_ptr *ptail,
fmpz_t m1, ulong e1,
fmpz_t m2, ulong e2)
{
pair_refine_unreduced(phead, m1, e1, m2, e2);
remove_ones(phead, ptail, *phead);
}
void
augment_refinement(fr_node_ptr *phead, fr_node_ptr *ptail,
const fmpz_t m_jp1, ulong e_jp1,
fr_node_ptr L_j, fr_node_ptr L_j_tail)
{
/*
* m_jp1 must have absolute value greater than 1,
* and its sign will be discarded.
* This function is destructive to the existing refinement L_j
*/
fr_node_ptr L_jp1, L_jp1_tail, L_prime, L_prime_tail, neo;
fmpz_t m;
ulong e;
/* initialize (m, e) <- (m_{j+1}, 1), L_{j+1} <- empty list */
fmpz_init(m);
fmpz_abs(m, m_jp1);
e = e_jp1;
L_jp1 = NULL;
L_jp1_tail = NULL;
L_prime = NULL;
L_prime_tail = NULL;
while (L_j && !fmpz_is_one(m))
{
if (!fr_node_is_one(L_j))
{
/* L' <- Pair-Refine((m, e), First(L_j)) */
pair_refine(&L_prime, &L_prime_tail, m, e,
fr_node_mref(L_j), L_j->e);
/* (m, e) <- First(L') */
fr_node_get_fmpz_ui(m, &e, L_prime);
/* L' <- Rest(L') */
fr_node_list_pop_front(&L_prime, &L_prime_tail);
/* L_{j+1} <- Concat(L_{j+1}, L') */
fr_node_list_concat(&L_jp1, &L_jp1_tail, L_prime, L_prime_tail);
}
/* L_j <- Rest(L_j) */
fr_node_list_pop_front(&L_j, &L_j_tail);
}
/* create a single-element list like [(m, e)] */
neo = flint_malloc(sizeof(fr_node_struct));
fr_node_init_fmpz_ui(neo, m, e);
/* L_{j+1} <- Concat(L_{j+1}, Rest(L_j), [(m, e)]) */
fr_node_list_pop_front(&L_j, &L_j_tail);
fr_node_list_concat(&L_jp1, &L_jp1_tail, L_j, L_j_tail);
fr_node_list_concat(&L_jp1, &L_jp1_tail, neo, neo);
/* output list of pairs (n_i, e_i) in L_{j+1} with n_i != 1 */
remove_ones(phead, ptail, L_jp1);
fmpz_clear(m);
}
int
_fmpz_factor_sgn(const fmpz_factor_t f)
{
int i, s;
ulong e, neg;
if (!f->sign)
{
return 0;
}
neg = f->sign < 0;
for (i = 0; i < f->num; i++)
{
/* follow the fmpz_pow_ui convention 0^0 = 1 */
e = f->exp[i];
if (e != WORD(0))
{
s = fmpz_sgn(f->p+i);
if (!s)
{
return 0;
}
else if (s < 0)
{
neg = (neg + e) % 2;
}
}
}
return neg ? -1 : 1;
}
void
fmpz_factor_refine(fmpz_factor_t res, const fmpz_factor_t f)
{
int s;
fr_node_ptr L, L_tail, curr;
slong i, len;
fmpz * b;
ulong e;
fr_node_ptr * qsort_arr;
/* check the sign of f without requiring canonical form */
s = _fmpz_factor_sgn(f);
if (!s)
{
_fmpz_factor_set_length(res, 0);
res->sign = 0;
return;
}
/* compute the refinement as a linked list */
for (L = L_tail = NULL, i = 0; i < f->num; i++)
{
b = f->p+i;
e = f->exp[i];
if (e != WORD(0) && !fmpz_is_pm1(b))
{
augment_refinement(&L, &L_tail, b, e, L, L_tail);
}
}
/* this is the length of the list of refined factors */
len = fr_node_list_length(L);
/* make a sorted array of pointers to the refined factor nodes */
qsort_arr = flint_malloc(sizeof(fr_node_ptr) * len);
for (i = 0, curr = L; i < len; i++, curr = curr->next)
{
qsort_arr[i] = curr;
}
qsort(qsort_arr, len, sizeof(fr_node_ptr), fr_node_base_pcmp);
/* set the length, sign, bases, and exponents of the output structure */
_fmpz_factor_fit_length(res, len);
_fmpz_factor_set_length(res, len);
res->sign = s;
for (i = 0; i < len; i++)
{
curr = qsort_arr[i];
fmpz_set(res->p+i, fr_node_mref(curr));
res->exp[i] = curr->e;
}
flint_free(qsort_arr);
fr_node_list_clear(L);
}