[flint-devel] factorise

[Richard Howell-P. <ric...@gm...>]

I include here 3 things, the improved square-free algorithm, some
documentation for it and some test code which it seems to pass ok. I have
noticed however that some polynomials cause it to put in way too many powers
in the exponents part, I think it might be to do with this part of the code:

*            if (res->num_factors) res->exponents[res->num_factors-1] *= i;
            i++;

*which I am not even sure what it is doing or why. I think it is a problem
when you input a polynomial with factors of high multiplicity, which maybe
some of the exponents get multiplied when they shouldn't but I'm not really
sure. I've also written a zmod_poly_faactorise function which takes care of
factoring and seems to work whenever square-free does.




/**
 *    This function takes an arbitary polynomial and factorises it. It first
performs a square-free factorisation, then
 * factorises all of the square free polynomails and returns the leading
coefficient, all the factors will be monic.
 */
unsigned long zmod_poly_factorise(zmod_poly_factor_t result, zmod_poly_t
input)
{
    //This will store all of the square-free factors
    zmod_poly_factor_t sqfree_factors;
    zmod_poly_factor_init(sqfree_factors);
    zmod_poly_t monic_input;
    zmod_poly_init(monic_input, zmod_poly_modulus(input));
    //Now we must make sure the input polynomial is monic. Get the highest
coeff and store it then call make monic
    unsigned long leading_coeff = zmod_poly_get_coeff_ui(input,
zmod_poly_degree(input));
    zmod_poly_make_monic(monic_input,input);

    //Now we do the square-free factorisation of the input polynomial
    zmod_poly_factor_square_free(sqfree_factors, monic_input);

    //space to store factors
    zmod_poly_factor_t* factors;
    factors = (zmod_poly_factor_t*) malloc (sizeof (*factors) *
sqfree_factors->num_factors );


    //Run berlekamp on each of the square-free factors
    for(unsigned long i = 0; i < sqfree_factors->num_factors; ++i)
    {
        zmod_poly_factor_init(factors[i]);
        zmod_poly_get_coeff_ui(sqfree_factors->factors[i],
zmod_poly_degree(sqfree_factors->factors[i]));


zmod_poly_make_monic(sqfree_factors->factors[i],sqfree_factors->factors[i]);
        zmod_poly_factor_berlekamp(factors[i], sqfree_factors->factors[i]);
        zmod_poly_factor_pow(factors[i], sqfree_factors->exponents[i]);

        //Now add all the factors to the final array
        zmod_poly_factor_concat(result, factors[i]);
        //clean up the memory being used
        zmod_poly_factor_clear(factors[i]);
    }

//clean up memory
   zmod_poly_factor_clear(sqfree_factors);
   free(factors);
    return leading_coeff;

}





/**
 * Computes the square free factorisation of the given polynomial. Input
must be a monic polynomial. All the factors will be monic.
 * @param f         The polynomial to factorise.
 * @param fac        The factorisation into monic square free factors of
<code>f</code>.
 */



/**
 * Square-Free Algorithm, takes an arbitary polynomial in F_p[X] and returns
an array of square free factors
 * LOW MULTIPLICITIES
 */
void zmod_poly_factor_square_free(zmod_poly_factor_t res, zmod_poly_t f)
{
    if (f->length <= 1)
    {
        res->num_factors = 0;
        return;
    }

    if (f->length == 2)
    {
        zmod_poly_factor_add(res, f);
        return;
    }
    unsigned long p = zmod_poly_modulus(f);        //order of the field
    unsigned long deg = zmod_poly_degree(f);    //degree of the polynomial

    //Step 1, look at f', if it is zero then we are done since f = h(x)^p
    //for some particular h(x), clearly f(x) = sum a_k x^kp, k <= deg(f)
    zmod_poly_t f_d, g, g_1;

    zmod_poly_init(g_1, p);
    zmod_poly_init(f_d, p);
    zmod_poly_init(g, p);

    zmod_poly_derivative(f_d, f);

    //CASE 1:
    if(zmod_poly_is_zero(f_d))
    {
        zmod_poly_t h;
               zmod_poly_init(h, p);
            for(unsigned long i = 0; i <= deg/p; ++i)    //this will be an
integer since f'=0
            {
                zmod_poly_set_coeff_ui(h, i, zmod_poly_get_coeff_ui(f,
i*p));
            }
        //now run square-free on h, and return it to the pth power
            zmod_poly_factor_t new_res;
            zmod_poly_factor_init(new_res);
        //h will be monic because the input is monic
            zmod_poly_factor_square_free(new_res, h);
            //now raise it to the power of p
            zmod_poly_factor_pow(new_res, p);

            zmod_poly_factor_concat(res, new_res);    //note, concatenating
is equivalent to multiplying
    }
    else
    {
            zmod_poly_gcd(g, f, f_d);
            zmod_poly_div(g_1, f, g);
            unsigned long i = 1;

            //CASE 2:
            while (!zmod_poly_is_one(g_1))
        {
            zmod_poly_t h, z;
                    zmod_poly_init(h, p);
                    zmod_poly_init(z, p);

            zmod_poly_gcd(h, g_1, g);
            zmod_poly_div(z, g_1, h);
            //Make sure the factor is monic
            zmod_poly_make_monic(z,z);
                    //now add it to the factor array
            zmod_poly_factor_add(res, z);
            if (res->num_factors) res->exponents[res->num_factors-1] *= i;
            i++;
            zmod_poly_set(g_1, h);
            zmod_poly_div(g, g, h);
            }

            if(!zmod_poly_is_one(g))
            {
            //so now we multiply res with square-free(g^1/p) ^ p
            zmod_poly_t g_p; //g^(1/p)
            zmod_poly_init(g_p, p);

            for(unsigned long i = 0; i <= zmod_poly_degree(g)/p; i++)
                zmod_poly_set_coeff_ui(g_p, i, zmod_poly_get_coeff_ui(g,
i*p));

            zmod_poly_factor_t new_res_2;
            zmod_poly_factor_init(new_res_2);
            //make sure g_p is monic
                    zmod_poly_make_monic(g_p,g_p);
            //square-free(g^(1/p))
            zmod_poly_factor_square_free(new_res_2, g_p);
            //now raise it to the power of p
            zmod_poly_factor_pow(new_res_2, p);
            zmod_poly_factor_concat(res, new_res_2);
            }
    }
}





int test_zmod_poly_factor_square_free()
{
   int result = 1;
   zmod_poly_t pol1,prod;
   zmod_poly_factor_t res;
   unsigned long bits;
   unsigned long modulus;

   for (unsigned long count1 = 0; (count1 < 10000) && (result == 1);
count1++)
   {
      bits = randint(FLINT_BITS-2)+2;

      do {modulus = randprime(bits);} while (modulus < 2);

 #if DEBUG
      printf("bits = %ld, modulus = %ld\n", bits, modulus);
#endif

    zmod_poly_init(pol1, modulus);
    zmod_poly_init(prod, modulus);

    zmod_poly_factor_init(res);

    zmod_poly_set_coeff_ui(pol1, randint(20), 1L);

    unsigned long k = zmod_poly_get_coeff_ui(pol1, zmod_poly_degree(pol1));
    zmod_poly_factor_square_free(res, pol1);
    //check the result
    zmod_poly_factor_product(prod, res);
    zmod_poly_scalar_mul(prod, prod, k);
    result &= zmod_poly_equal(prod, pol1);

    zmod_poly_clear(pol1);
    zmod_poly_clear(prod);
    zmod_poly_factor_clear(res);
   }
   return result;
}



-- 
Richard