flint-devel — Development list for the FLINT project

Re: [flint-devel] improved differentiation function

[David H. <dmh...@ma...>]

Also I'm guessing it probably won't work as is, since the arguments  
to mulmod_precomp probably need to be less than p (which j is not  
necessarily).

david

On Jul 16, 2008, at 10:56 AM, Richard Howell-Peak wrote:

> Can we do better than this?
>
>
>
> void zmod_poly_diff(zmod_poly_t x_primed, zmod_poly_t x)
> {
>     //get the degree and base of the polynomial we are working with
>     unsigned long degree = zmod_poly_length (x);
>     unsigned long p = zmod_poly_modulus(x);
>     unsigned long coeff;
>
>     pre_inv_t p_inv = z_precompute_inverse(p);
>     for (unsigned long j = 1; j < degree ; ++j)
>     {
>         coeff = z_mulmod_precomp ( zmod_poly_get_coeff_ui ( x ,  
> j) ,  j , p , p_inv );
>         zmod_poly_set_coeff_ui ( x_primed, j-1, coeff );
>     }
>
> }

[flint-devel] improved differentiation function

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

oh yeah, I know how to do that

-- 
Richard

Re: [flint-devel] improved differentiation function

[David H. <dmh...@ma...>]

You should pre-allocate space for x_primed. I don't know how to do  
that; perhaps Bill knows.

david

On Jul 16, 2008, at 10:56 AM, Richard Howell-Peak wrote:

> Can we do better than this?
>
>
>
> void zmod_poly_diff(zmod_poly_t x_primed, zmod_poly_t x)
> {
>     //get the degree and base of the polynomial we are working with
>     unsigned long degree = zmod_poly_length (x);
>     unsigned long p = zmod_poly_modulus(x);
>     unsigned long coeff;
>
>     pre_inv_t p_inv = z_precompute_inverse(p);
>     for (unsigned long j = 1; j < degree ; ++j)
>     {
>         coeff = z_mulmod_precomp ( zmod_poly_get_coeff_ui ( x ,  
> j) ,  j , p , p_inv );
>         zmod_poly_set_coeff_ui ( x_primed, j-1, coeff );
>     }
>
> }

[flint-devel] improved differentiation function

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

Can we do better than this?



void zmod_poly_diff(zmod_poly_t x_primed, zmod_poly_t x)
{
    //get the degree and base of the polynomial we are working with
    unsigned long degree = zmod_poly_length (x);
    unsigned long p = zmod_poly_modulus(x);
    unsigned long coeff;

    pre_inv_t p_inv = z_precompute_inverse(p);
    for (unsigned long j = 1; j < degree ; ++j)
    {
        coeff = z_mulmod_precomp ( zmod_poly_get_coeff_ui ( x , j) ,  j , p
, p_inv );
        zmod_poly_set_coeff_ui ( x_primed, j-1, coeff );
    }

}

-- 
Richard

Re: [flint-devel] Gaussian Elimination

[Bill H. <goo...@go...>]

Oh, yes you are right!!

In the given example, you have nullity 3, which means there should be
3 things in the basis of the null space.

I'll presume that the polynomials have constant coefficients to the
right of the matrix.

It's probably easier if you go all the way to Gauss Jordan
elimination, i.e. zeroes above all the pivots (1's on the diagonals).

You can take your free variables to be the ones corresponding to the
last three columns (since they don't have pivots).

You can set these to anything you like, e.g. 0,0,1 or 0,1,0, or 1,0,0.
Then you work back through the rows for each of these three
combinations and solve for the other variables. This gives you three
vectors which form the basis you are after.

You might not get the same basis as this guy, but you'll get a basis,
and that'll do.

I think if you reduce to Gauss-Jordan canonical form it is easier. I
think if you negate the last three columns mod p and put

1 0 0
0 1 0
0 0 1

in the bottom right corner of the matrix, then the last three columns
of the matrix are the vectors you are after.

Bill.

2008/7/16 Richard Howell-Peak <ric...@gm...>:
> I thought and example might help, the one in the paper gives the matrix
>
> 0    7
> 8    7  0 6 0 0 0 0 0 1
> 10  7  6 6 4 3 1 2 6 2 3 1
> 8    7  4 6 6 0 5 4 3 4
> 10  7  3 5 3 4 2 5 5 2 2 5
> 10  7  0 6 1 1 6 6 0 2 6 6
> 10  7  4 6 3 1 0 4 0 5 2 5
> 10  7  4 2 0 1 4 3 6 1 2 4
> 10  7  2 2 4 1 6 6 2 4 1 5
> 10  7  6 3 5 5 1 6 6 4 6 5
>
> And after Gaussian Elimination I get this:
>
> 10  7  1 1 3 4 6 5 1 5 4 6
> 8    7  0 1 0 0 0 0 0 6
> 10  7  0 0 1 4 5 1 2 3 0 6
> 10  7  0 0 0 1 1 6 1 3 0 2
> 10  7  0 0 0 0 1 1 1 3 5 1
> 10  7  0 0 0 0 0 1 5 6 0 6
> 10  7  0 0 0 0 0 0 1 2 1 4
> 0    7
> 0    7
> 0    7
>
> He then gives the basis as this:
>
> b1 = 1
> b2 = 2x^3 + x^4 + 3x^6 + 6x^7 + x^8
> b3 = 6x + 2x^2 + 6x^3 + 6x^4 + x^5 + 3x^6 + 3x^7 + x^9
>
> I can see why there are 3 of them since there are 3 zero rows, i.e. the
> nullity is 3, but I can't see how to calculate the basis.
> --
> Richard
> -------------------------------------------------------------------------
> This SF.Net email is sponsored by the Moblin Your Move Developer's challenge
> Build the coolest Linux based applications with Moblin SDK & win great
> prizes
> Grand prize is a trip for two to an Open Source event anywhere in the world
> http://moblin-contest.org/redirect.php?banner_id=100&url=/
> _______________________________________________
> flint-devel mailing list
> fli...@li...
> https://lists.sourceforge.net/lists/listinfo/flint-devel
>
>

[flint-devel] Gaussian Elimination

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

I thought and example might help, the one in the paper gives the matrix

0    7
8    7  0 6 0 0 0 0 0 1
10  7  6 6 4 3 1 2 6 2 3 1
8    7  4 6 6 0 5 4 3 4
10  7  3 5 3 4 2 5 5 2 2 5
10  7  0 6 1 1 6 6 0 2 6 6
10  7  4 6 3 1 0 4 0 5 2 5
10  7  4 2 0 1 4 3 6 1 2 4
10  7  2 2 4 1 6 6 2 4 1 5
10  7  6 3 5 5 1 6 6 4 6 5

And after Gaussian Elimination I get this:

10  7  1 1 3 4 6 5 1 5 4 6
8    7  0 1 0 0 0 0 0 6
10  7  0 0 1 4 5 1 2 3 0 6
10  7  0 0 0 1 1 6 1 3 0 2
10  7  0 0 0 0 1 1 1 3 5 1
10  7  0 0 0 0 0 1 5 6 0 6
10  7  0 0 0 0 0 0 1 2 1 4
0    7
0    7
0    7

He then gives the basis as this:

b1 = 1
b2 = 2x^3 + x^4 + 3x^6 + 6x^7 + x^8
b3 = 6x + 2x^2 + 6x^3 + 6x^4 + x^5 + 3x^6 + 3x^7 + x^9

I can see why there are 3 of them since there are 3 zero rows, i.e. the
nullity is 3, but I can't see how to calculate the basis.
-- 
Richard

[flint-devel] Gaussian Elimination

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

Well, I thought the whole point was to compute the kernel of the map, which
would mean we are more interested in the null space and the nullity, opposed
to the rank and the image or am I wrong?

-- 
Richard

Re: [flint-devel] Gaussian Elimination

[Bill H. <goo...@go...>]

Hi Richard,

Can you do the Gaussian elimination itself? It's just like GE over Q
except you are working over F_p (there is a typo in that thesis which
says you are working over F_q[x]). The elementary row operations are
to add two rows together, for which you need z_addmod, to multiply a
row by a scalar, for which you need z_mulmod_precomp and finally to
divide a row through by a scalar, for which you use z_invert and
z_mulmod_precomp.

The number of non-zero rows is the rank. The non-zero rows form the
basis you are after. Each row represents the coefficients of a
polynomial in the Berlekamp basis that you are after.

If this doesn't make sense, come back with more questions.

Bill.

2008/7/16 Richard Howell-Peak <ric...@gm...>:
> I was wrong, I don't see how to do it at all. I know the rank of the matrix
> is to do with the non-zero rows, but I can't see how to get a basis from
> that.
>
> Richard
> -------------------------------------------------------------------------
> This SF.Net email is sponsored by the Moblin Your Move Developer's challenge
> Build the coolest Linux based applications with Moblin SDK & win great
> prizes
> Grand prize is a trip for two to an Open Source event anywhere in the world
> http://moblin-contest.org/redirect.php?banner_id=100&url=/
> _______________________________________________
> flint-devel mailing list
> fli...@li...
> https://lists.sourceforge.net/lists/listinfo/flint-devel
>
>

Re: [flint-devel] Progress

[Bill H. <goo...@go...>]

Excellent work. The figures you are giving look quite promising.

I'll commit the changes you've made to pocklington to the FLINT repository soon.

Currently I'm working on some code for a collaborator which computes
theta functions (essentially power series like Sum q^{ax^2+bx+c} for
various values of a,, b, c). As soon as I've got somewhere with it,
I'll put some more time into the primality testing code.

Bill.

2008/7/16 Peter Shrimpton <ps...@gm...>:
> Hello.
>
> I think I have pretty much finished the pocklington test. The partial
> factorization is done by adjusting the code for the z_factor function that
> was already there. This seams to be pretty quick. Also I have combined it
> with the SPRP test already in flint and the results seam quite promising. I
> have checked the numbers between 10^16 and 10^16+10^9 and the combined SPRP
> and Pocklington tese took on average about 2*10^-6 seconds per number to
> check with 100 iterations of the pocklington test. The number of iterations
> in the pocklington seams to affect the number of potential SPRP's we get
> quite a bit:
>
> number of iterations           number of primes proved        number of
> 2-SPRP's       number of 2-SPRP's that we couldn't prove are prime
> 200
> 1033994                                   1033994                       0
> 130
> 1033986                                   1033994                       8
> 120
> 1033982                                   1033994                       12
> 100
> 1033938                                   1033994                       56
> 90
> 1033897                                   1033994                       97
>
> As the number of composite SPRP's found was 0 in this sample the amount of
> time taken to execute the code was affected very little by the number of
> iterations of the Pocklington test. This is because quite often the first
> few numbers are enough to find a b such that gcd(b^((n-1)/p),n). However I
> am not sure if this trend will continue. The code for this is attached.
>
>
> I am now working on the Lucas psudoprime test. I have implemented a slow
> version that works up to about 100 but then the numbers in the lucas chain
> get too big. There is an algorithm in the book "Prime Numbers. A
> Computational perspective" (one of the authors is in fact Carl Pomerance)
> which keeps everything reduced mod n and should be quite quick but I am
> having difficulty in implementing it.
>
> Peter
>
> -------------------------------------------------------------------------
> This SF.Net email is sponsored by the Moblin Your Move Developer's challenge
> Build the coolest Linux based applications with Moblin SDK & win great
> prizes
> Grand prize is a trip for two to an Open Source event anywhere in the world
> http://moblin-contest.org/redirect.php?banner_id=100&url=/
> _______________________________________________
> flint-devel mailing list
> fli...@li...
> https://lists.sourceforge.net/lists/listinfo/flint-devel
>
>

[flint-devel] Progress

[Peter S. <ps...@gm...>]

Hello.

I think I have pretty much finished the pocklington test. The partial
factorization is done by adjusting the code for the z_factor function that
was already there. This seams to be pretty quick. Also I have combined it
with the SPRP test already in flint and the results seam quite promising. I
have checked the numbers between 10^16 and 10^16+10^9 and the combined SPRP
and Pocklington tese took on average about 2*10^-6 seconds per number to
check with 100 iterations of the pocklington test. The number of iterations
in the pocklington seams to affect the number of potential SPRP's we get
quite a bit:

number of iterations           number of primes proved        number of
2-SPRP's       number of 2-SPRP's that we couldn't prove are prime
200
1033994                                   1033994                       0
130
1033986                                   1033994                       8
120
1033982                                   1033994                       12
100
1033938                                   1033994                       56
90
1033897                                   1033994                       97

As the number of composite SPRP's found was 0 in this sample the amount of
time taken to execute the code was affected very little by the number of
iterations of the Pocklington test. This is because quite often the first
few numbers are enough to find a b such that gcd(b^((n-1)/p),n). However I
am not sure if this trend will continue. The code for this is attached.


I am now working on the Lucas psudoprime test. I have implemented a slow
version that works up to about 100 but then the numbers in the lucas chain
get too big. There is an algorithm in the book "Prime Numbers. A
Computational perspective" (one of the authors is in fact Carl Pomerance)
which keeps everything reduced mod n and should be quite quick but I am
having difficulty in implementing it.

Peter

[flint-devel] Gaussian Elimination

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

I was wrong, I don't see how to do it at all. I know the rank of the matrix
is to do with the non-zero rows, but I can't see how to get a basis from
that.

Richard

[flint-devel] Gaussian Elimination

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

I just remembered some linear algebra from the first year. I think I can do
it now.

-- 
Richard

[flint-devel] Gaussian Elimination

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

I have implemented a Gaussian Elimination algorithm from wikipedia which
seems to work ok, but somehow from it I need to obtain the basis for the
Berlekamp Sub-algebra. I understand this is something to do with the columns
of the matrix but I can't quite see what's going on here.

-- 
Richard

Re: [flint-devel] Square-Free Algorithm

[Bill H. <goo...@go...>]

Ah, I see how it all works now.

By the way, the function zmod_poly_set_ui doesn't reduce mod p. You
have to use the function in long_extras called z_mulmod_precomp in
your diff function. The other problem of course is that if p is nearly
64 bits then coeff*(i%p) might be well above 64 bits, so it won't fit
in a limb. z_mulmod_precomp takes care of that for you (or one of the
versions does anyway).

The other thing is that i%p takes a lot of cycles in your diff
function. Better to start j = 1, increment j every time and have, if
(j == p) j = 0; in your loop. That way, every time it gets to p it
will reset to 0.

There's probably no need to check if coeff == 0 or not. This probably
happens in z_mulmod_precomp and if it doesn't, it might not be faster.

Bill.

2008/7/15 Richard Howell-Peak <ric...@gm...>:
> I agree, the goto's are not very nice it's just how the algorithm was
> written in the paper. I've basically got some code for doing the other
> stuff, you can have a look at it here. Obviously none of it's optimised but
> I'd quite like to get a working version of Berlekamp as soon as possible
> then worry about optimisation later.
>
>
> /**
>  * Initialises an array of zmod_poly's
>  */
> void zmod_poly_arr_init(zmod_poly_arr_t arr, unsigned long p)
> {
>     arr->len = 1;
>     arr->used = 0;
>     arr->factors = (zmod_poly_t*) malloc(sizeof(zmod_poly_t));
>     zmod_poly_init(arr->factors[0],p);
> }
>
> /**
>  * Initialises one of specified length
>  */
> void zmod_poly_arr_init_len(zmod_poly_arr_t arr, unsigned long p, unsigned
> long length)
> {
>     arr->len = length;
>     arr->used = 0;
>     arr->factors = (zmod_poly_t*) malloc(sizeof(zmod_poly_t) * length);
>     for(unsigned long i=0; i< length; i++)
>         zmod_poly_init(arr->factors[i],p);
> }
>
> /**
>  * Frees up memory being used
>  */
> void zmod_poly_arr_free(zmod_poly_arr_t arr)
> {
>     free(arr->factors);
>     arr->len = 0;
>     arr->used = 0;
> }
>
> /**
>  * Adds an extra element to the array
>  */
> void zmod_poly_arr_add(zmod_poly_arr_t arr, zmod_poly_t add)
> {
>     //how much space left in the array?, if none make a new one twice as big
> (for efficiency) and copy contents across
>     if(arr->len == arr->used)
>     {
>         zmod_poly_arr_t temp;
>         zmod_poly_arr_init_len(temp, zmod_poly_modulus(add), arr->len * 2);
>         //copy it across
>         for(unsigned long i=0; i < arr->len; i++)
>         {
>             zmod_poly_set(temp->factors[i], arr->factors[i]);
>         }
>         //add in the extra one
>         temp->len = arr->len*2;
>         zmod_poly_set(temp->factors[arr->used], add);
>         temp->used = 1 + arr->used;
>         zmod_poly_arr_free(arr);
>         arr[0] = temp[0];
>     } //plenty of space left, just add it to the end
>     else
>     {
>         zmod_poly_set(arr->factors[arr->used], add);
>         ++arr->used;
>     }
> }
>
> /**
>  * Concatenates array res to equal array res + add
>  */
> void zmod_poly_arr_cat(zmod_poly_arr_t res, zmod_poly_arr_t add)
> {
>     for(unsigned long i=0; i < add->used; ++i)
>         zmod_poly_arr_add(res, add->factors[i]);
> }
>
> //Getters and Setters
> /**
>  * Get the number of elements currently in the array
>  */
> unsigned long zmod_poly_arr_get_size(zmod_poly_arr_t arr)
> {
>     return arr->used;
> }
>
> /**
>  * Gets the product of all the factors
>  */
> void zmod_poly_arr_product(zmod_poly_t result, zmod_poly_arr_t arr)
> {
>     zmod_poly_init(result, zmod_poly_modulus(arr->factors[0]) );
>     zmod_poly_set_coeff_ui(result, 0 , 1);
>     for(unsigned long i=0; i < arr->used; ++i)
>         zmod_poly_mul(result, result, arr->factors[i]);
> }
>
> /**
>  * Dumps the array to stdout
>  */
> void zmod_poly_arr_print(zmod_poly_arr_t output)
> {
>     printf("\n");
>     for(unsigned long i=0; i < output->used; ++i)
>     {
>         zmod_poly_print(output->factors[i]);
>         printf("\t");
>     }
>     printf("\n");
> }
>
> /**
>  * Prints the product of all elements in the array
>  */
> void zmod_poly_arr_print_res(zmod_poly_arr_t output)
> {
>     zmod_poly_t res;
>     zmod_poly_init(res, zmod_poly_modulus(output->factors[0]));
>     zmod_poly_set_coeff_ui(res, 0, 1);
>     for(unsigned long i=0; i < output->used; ++i)
>     {
>         zmod_poly_mul(res, res, output->factors[i]);
>
>     }
>     zmod_poly_print(res);
> }
>
> /**************************** Utility functions for zmod_poly's
> *****************/
>
>
> /**
>  * Returns the derivative of a polynomial mod p, where p is the prime
>  * that is already associated with the polynomial x.
>  * <p>
>  * This is done in a very naive way, simply going through the polynomial
>  * differentiating term by term reducing the product mod p each time
>  *
>  * @param  zmod_poly_t The polynomial to be derived
>  * @return      The derivative of the original polynomial
>  */
> void naiveDiff(zmod_poly_t x_primed, zmod_poly_t x)
> {
> //Now we compute the derivative in a very naive way :)
>
>     unsigned long degree = zmod_poly_length (x);
>     long length = zmod_poly_degree (x);
>
>     zmod_poly_init( x_primed , zmod_poly_modulus(x));
>     for (unsigned long i=1; i<degree; i++)
>     {
>         unsigned long coeff = zmod_poly_get_coeff_ui ( x , i);
>         unsigned long newCoeff = i*coeff;
>         zmod_poly_set_coeff_ui ( x_primed, i-1, newCoeff);
>     }
> }
>
>
> /**
>  * This is slightly more efficient but basically does the same thing
>  */
> void zmod_poly_diff(zmod_poly_t x_primed, zmod_poly_t x)
> {
>     //get the degree of the polynomial we are working with
>     unsigned long degree = zmod_poly_length (x);
>     //make sure new poly is in the same base
>     unsigned long p = zmod_poly_modulus(x);
>     zmod_poly_init( x_primed , p);
>     unsigned long coeff;
>     //Loop through and diff term by term
>     for (unsigned long i=1; i<degree; i++)
>     {
>         coeff = zmod_poly_get_coeff_ui ( x , i);
>         if(coeff!=0)
>             zmod_poly_set_coeff_ui ( x_primed, i-1, (i%p) *coeff);
>     }
> }
>
>
> /**
>  * Repeated square algorithm, efficiently computes a to the exp mod <f>
> using the repeated square method
>  */
> void zmod_poly_repeated_square(zmod_poly_t res, zmod_poly_t a, unsigned long
> exp, zmod_poly_t f)
> {
>     //unefficient way, quick fix:
>     zmod_poly_t a_exp,temp;
>     zmod_poly_init(a_exp, zmod_poly_modulus(f) );
>     zmod_poly_init(temp, zmod_poly_modulus(f) );
>     zmod_poly_init(res , zmod_poly_modulus(f) );
>
>     zmod_poly_set(a_exp,a);
>     for(unsigned long i=1; i<exp; i++)         //calculate a_exp = a^exp
>         zmod_poly_mul(a_exp ,a_exp ,a);
>     //for A = QB + R we get divrem(Q, R, A, B);
>     zmod_poly_divrem( temp , res, a_exp , f );
> }
>
> /**
>  * Raises a polynomial to the nth power
>  * Sets res to f^n
>  */
> void zmod_poly_power(zmod_poly_t res, zmod_poly_t f, unsigned long n)
> {
>     zmod_poly_set(res, f);
>     for(unsigned long i=1; i < n; ++i)
>         zmod_poly_mul(res, res, f);
> }
>
> --
> Richard
> -------------------------------------------------------------------------
> This SF.Net email is sponsored by the Moblin Your Move Developer's challenge
> Build the coolest Linux based applications with Moblin SDK & win great
> prizes
> Grand prize is a trip for two to an Open Source event anywhere in the world
> http://moblin-contest.org/redirect.php?banner_id=100&url=/
> _______________________________________________
> flint-devel mailing list
> fli...@li...
> https://lists.sourceforge.net/lists/listinfo/flint-devel
>
>

[flint-devel] Square-Free Algorithm

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

This is what I've done so far. It almost works perfectly, except for some
reason it is always out by some factor of x, i.e. 4x or something like that.
I have no idea why this is happening and I can't work it out at all. I'm
currently moving on to doing Berlekamp's algorithm which in order to do that
I'll need to do some Gaussian Elimination but it doesn't look too hard so
I'll get on with that for the time being.

*<https://sourceforge.net/mailarchive/forum.php?thread_name=d42dea700807150320r57bb77fj25c6706c971e94e%40mail.gmail.com&forum_name=flint-devel>
/**
 * Square-Free Algorithm, takes an arbitary polynomial in Fq[X] and returns
an array of square free factors
 * LOW MULTIPLICITIES
 */
void zmod_poly_square_free(zmod_poly_arr_t res, zmod_poly_t f)
{
    if(DBG) {NL printf("STARTING \n");}
    unsigned long q = zmod_poly_modulus(f);    //order of the field
    unsigned long deg = zmod_poly_degree(f); //degree of the polynomial
    //Let's create some constant polynomials that will be very useful later
on:
    zmod_poly_t zero, one;
    zmod_poly_init(zero,q);
    zmod_poly_init(one,q);
    zmod_poly_zero(zero);
    zmod_poly_set_coeff_ui(one, 0 , 1 );

    //set the result to the constant polynomial 1
    //zmod_poly_arr_add(res,one);

    //Step 1, look at f', if it is zero then we are done since f = h(x)^q
    //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,q);
    zmod_poly_init(f_d,q);
    zmod_poly_init(g,q);
    zmod_poly_diff(f_d, f);
    //CASE 1:
    if( zmod_poly_equal(f_d, zero) )
    {
        if(DBG) {printf("f'=0, f = "); zmod_poly_print(f); NL}
        unsigned long coeff[deg/q];
        zmod_poly_t h;
        zmod_poly_init(h,q);
        for(unsigned long i=0; i <= deg/q; ++i)    //this will be an integer
since f'=0
        {
            zmod_poly_set_coeff_ui(h, i, zmod_poly_get_coeff_ui(f, i*q));
        }
        if(DBG) {printf("h = ");zmod_poly_print(h); NL}
        //now run square-free on h, and return it to the pth power
        zmod_poly_arr_t new_res;
        zmod_poly_arr_init(new_res,q);
        zmod_poly_square_free(new_res, h);

        for(unsigned long j=0; j < q ; ++j)
        {
            zmod_poly_arr_cat(res, new_res);    //note, concatenating is
equivalent to multiplying
        }
    }
    else
    {
        zmod_poly_gcd(g, f, f_d);
        zmod_poly_div(g_1, f, g);

        //CASE 2:
    in:
        if(DBG) printf("Is g_1 = 1? \n");
        if( zmod_poly_equal(g_1, one))
        {
            goto out;
        }
        else
        {
            if(DBG){ printf("no, g_1 = "); zmod_poly_print(g_1); NL }
            zmod_poly_t h, z;
            zmod_poly_init(h,q);
            zmod_poly_init(z,q);
            //moar!!
            zmod_poly_gcd(h, g_1, g);
            zmod_poly_div(z, g_1, h);
            //out <- out.z
            if(DBG) { printf("z = "); zmod_poly_print(z); NL }
            zmod_poly_arr_add(res, z);
            zmod_poly_set(g_1, h);
            zmod_poly_div(g, g, h);
            if(DBG){ printf("also, g = "); zmod_poly_print(g); NL }
        }
        goto in;

    out:
        if(DBG) printf("g_1 = 1 \n");

        if( ! zmod_poly_equal(g,one) )
        {
            if(DBG) printf("g != 1 \n");
            //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, q);
            for(unsigned long i=0; i <= zmod_poly_degree(g) / q; ++i)
                zmod_poly_set_coeff_ui(g_p, i, zmod_poly_get_coeff_ui(g,
i*q));

            zmod_poly_arr_t new_res_2;
            zmod_poly_arr_init(new_res_2,q);
            //square-free(g^(1/p))
            zmod_poly_square_free(new_res_2, g_p);
            //now raise it to the power of p
            for(unsigned long i=0; i < q; ++i)
                zmod_poly_arr_cat(res,new_res_2);
        }
    }
    if(DBG){ N: printf("FINISH! \n"); NL}
}*

-- 
Richard

Re: [flint-devel] Square-Free Algorithm

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

I agree, the goto's are not very nice it's just how the algorithm was
written in the paper. I've basically got some code for doing the other
stuff, you can have a look at it here. Obviously none of it's optimised but
I'd quite like to get a working version of Berlekamp as soon as possible
then worry about optimisation later.


/**
 * Initialises an array of zmod_poly's
 */
void zmod_poly_arr_init(zmod_poly_arr_t arr, unsigned long p)
{
    arr->len = 1;
    arr->used = 0;
    arr->factors = (zmod_poly_t*) malloc(sizeof(zmod_poly_t));
    zmod_poly_init(arr->factors[0],p);
}

/**
 * Initialises one of specified length
 */
void zmod_poly_arr_init_len(zmod_poly_arr_t arr, unsigned long p, unsigned
long length)
{
    arr->len = length;
    arr->used = 0;
    arr->factors = (zmod_poly_t*) malloc(sizeof(zmod_poly_t) * length);
    for(unsigned long i=0; i< length; i++)
        zmod_poly_init(arr->factors[i],p);
}

/**
 * Frees up memory being used
 */
void zmod_poly_arr_free(zmod_poly_arr_t arr)
{
    free(arr->factors);
    arr->len = 0;
    arr->used = 0;
}

/**
 * Adds an extra element to the array
 */
void zmod_poly_arr_add(zmod_poly_arr_t arr, zmod_poly_t add)
{
    //how much space left in the array?, if none make a new one twice as big
(for efficiency) and copy contents across
    if(arr->len == arr->used)
    {
        zmod_poly_arr_t temp;
        zmod_poly_arr_init_len(temp, zmod_poly_modulus(add), arr->len * 2);
        //copy it across
        for(unsigned long i=0; i < arr->len; i++)
        {
            zmod_poly_set(temp->factors[i], arr->factors[i]);
        }
        //add in the extra one
        temp->len = arr->len*2;
        zmod_poly_set(temp->factors[arr->used], add);
        temp->used = 1 + arr->used;
        zmod_poly_arr_free(arr);
        arr[0] = temp[0];
    } //plenty of space left, just add it to the end
    else
    {
        zmod_poly_set(arr->factors[arr->used], add);
        ++arr->used;
    }
}

/**
 * Concatenates array res to equal array res + add
 */
void zmod_poly_arr_cat(zmod_poly_arr_t res, zmod_poly_arr_t add)
{
    for(unsigned long i=0; i < add->used; ++i)
        zmod_poly_arr_add(res, add->factors[i]);
}

//Getters and Setters
/**
 * Get the number of elements currently in the array
 */
unsigned long zmod_poly_arr_get_size(zmod_poly_arr_t arr)
{
    return arr->used;
}

/**
 * Gets the product of all the factors
 */
void zmod_poly_arr_product(zmod_poly_t result, zmod_poly_arr_t arr)
{
    zmod_poly_init(result, zmod_poly_modulus(arr->factors[0]) );
    zmod_poly_set_coeff_ui(result, 0 , 1);
    for(unsigned long i=0; i < arr->used; ++i)
        zmod_poly_mul(result, result, arr->factors[i]);
}

/**
 * Dumps the array to stdout
 */
void zmod_poly_arr_print(zmod_poly_arr_t output)
{
    printf("\n");
    for(unsigned long i=0; i < output->used; ++i)
    {
        zmod_poly_print(output->factors[i]);
        printf("\t");
    }
    printf("\n");
}

/**
 * Prints the product of all elements in the array
 */
void zmod_poly_arr_print_res(zmod_poly_arr_t output)
{
    zmod_poly_t res;
    zmod_poly_init(res, zmod_poly_modulus(output->factors[0]));
    zmod_poly_set_coeff_ui(res, 0, 1);
    for(unsigned long i=0; i < output->used; ++i)
    {
        zmod_poly_mul(res, res, output->factors[i]);

    }
    zmod_poly_print(res);
}

/**************************** Utility functions for zmod_poly's
*****************/


/**
 * Returns the derivative of a polynomial mod p, where p is the prime
 * that is already associated with the polynomial x.
 * <p>
 * This is done in a very naive way, simply going through the polynomial
 * differentiating term by term reducing the product mod p each time
 *
 * @param  zmod_poly_t The polynomial to be derived
 * @return      The derivative of the original polynomial
 */
void naiveDiff(zmod_poly_t x_primed, zmod_poly_t x)
{
//Now we compute the derivative in a very naive way :)

    unsigned long degree = zmod_poly_length (x);
    long length = zmod_poly_degree (x);

    zmod_poly_init( x_primed , zmod_poly_modulus(x));
    for (unsigned long i=1; i<degree; i++)
    {
        unsigned long coeff = zmod_poly_get_coeff_ui ( x , i);
        unsigned long newCoeff = i*coeff;
        zmod_poly_set_coeff_ui ( x_primed, i-1, newCoeff);
    }
}


/**
 * This is slightly more efficient but basically does the same thing
 */
void zmod_poly_diff(zmod_poly_t x_primed, zmod_poly_t x)
{
    //get the degree of the polynomial we are working with
    unsigned long degree = zmod_poly_length (x);
    //make sure new poly is in the same base
    unsigned long p = zmod_poly_modulus(x);
    zmod_poly_init( x_primed , p);
    unsigned long coeff;
    //Loop through and diff term by term
    for (unsigned long i=1; i<degree; i++)
    {
        coeff = zmod_poly_get_coeff_ui ( x , i);
        if(coeff!=0)
            zmod_poly_set_coeff_ui ( x_primed, i-1, (i%p) *coeff);
    }
}


/**
 * Repeated square algorithm, efficiently computes a to the exp mod <f>
using the repeated square method
 */
void zmod_poly_repeated_square(zmod_poly_t res, zmod_poly_t a, unsigned long
exp, zmod_poly_t f)
{
    //unefficient way, quick fix:
    zmod_poly_t a_exp,temp;
    zmod_poly_init(a_exp, zmod_poly_modulus(f) );
    zmod_poly_init(temp, zmod_poly_modulus(f) );
    zmod_poly_init(res , zmod_poly_modulus(f) );

    zmod_poly_set(a_exp,a);
    for(unsigned long i=1; i<exp; i++)         //calculate a_exp = a^exp
        zmod_poly_mul(a_exp ,a_exp ,a);
    //for A = QB + R we get divrem(Q, R, A, B);
    zmod_poly_divrem( temp , res, a_exp , f );
}

/**
 * Raises a polynomial to the nth power
 * Sets res to f^n
 */
void zmod_poly_power(zmod_poly_t res, zmod_poly_t f, unsigned long n)
{
    zmod_poly_set(res, f);
    for(unsigned long i=1; i < n; ++i)
        zmod_poly_mul(res, res, f);
}

-- 
Richard

Re: [flint-devel] Square-Free Algorithm

[Bill H. <goo...@go...>]

Oh wait, I see you zmod_poly_arr_add is just the same as zmod_poly_set
for now. You are actually right, we don't want zmod_poly_arr_cat
there.

Bill.

2008/7/15 Bill Hart <goo...@go...>:
> I think the second part of the algorithm is OK, except you probably
> want zmod_poly_arr_cat(res, z) not zmod_poly_arr_add(res, z).
>
> If I've understood correctly, the way you've implemented this it is
> self testing, i.e. the output should return the same thing as the
> input. That's quite clever.
>
> Later on we'll set up a struct for containing the factors and we'll
> add some functions for catenating factors to such a struct.
>
> One thing I'd like to get rid of is the goto statements. I think what
> you have is equivalent to somethink like:
>
> while (g1 == 1)
> {
>   blah;
> }
>
> blah blah;
>
> Bill.
>
> 2008/7/15 Bill Hart <goo...@go...>:
>> As you worked out, to take f(x)^p, simply raise each term to the power
>> of p. That means raising each coefficient to the power p modulo p (use
>> Fermat's little theorem for that) and multiplying the exponent of each
>> term by p.
>>
>> When trying to go the other way, i.e. compute f(x)^(1/p),  it's just
>> the reverse. I think you got the reverse of the second part OK, but
>> not the first.
>>
>> Bill.
>>
>> 2008/7/15 Richard Howell-Peak <ric...@gm...>:
>>> /**
>>>  * Square-Free Algorithm, takes an arbitary polynomial in Fq[X] and returns
>>> an array of square free factors
>>>  * LOW MULTIPLICITIES
>>>  */
>>> #define DBG 1
>>> void zmod_poly_square_free(zmod_poly_arr_t res, zmod_poly_t f)
>>> {
>>>     unsigned long q = zmod_poly_modulus(f);    //order of the field
>>>     unsigned long deg = zmod_poly_degree(f); //degree of the polynomial
>>>     //Let's create some constant polynomials that will be very useful later
>>> on:
>>>     zmod_poly_t zero, one;
>>>     zmod_poly_init(zero,q);
>>>     zmod_poly_init(one,q);
>>>     zmod_poly_zero(zero);
>>>     zmod_poly_set_coeff_ui(one, 0 , 1 );
>>>
>>>     //set the result to the constant polynomial 1
>>>     zmod_poly_arr_init(res,q);    //? Not sure if we should assume it is
>>> pre-initialised (NOTE - probably won't)
>>>     zmod_poly_arr_add(res,one);
>>>
>>>     //Step 1, look at f', if it is zero then we are done since f = h(x)^q
>>>     //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,q);
>>>     zmod_poly_init(f_d,q);
>>>     zmod_poly_init(g,q);
>>>     zmod_poly_diff(f_d, f);
>>>     //CASE 1:
>>>     if( zmod_poly_equal(f_d, zero) )
>>>     {
>>>         if(DBG) printf("f'=0, f = "); zmod_poly_print(f); NL
>>>         unsigned long coeff[deg/q];
>>>         zmod_poly_t h;
>>>         zmod_poly_init(h,q);
>>>         for(unsigned long i=0; i <= deg/q; ++i)    //this will be an integer
>>> since f'=0
>>>         {
>>>             if(DBG) printf("Changing new coefficients\n");
>>>             coeff[i] = zmod_poly_get_coeff_ui(f, i*q);
>>>             //change them slightly
>>>             coeff[i] = coeff[i] ^ (i/q);    //clearly this is incorrect,
>>> need to do math.pow but for finite field
>>>             //now we set up the h poly, where f = h^q
>>>             zmod_poly_set_coeff_ui(h, i, coeff[i]);
>>>         }
>>>         if(DBG) printf("h = ");zmod_poly_print(h); NL
>>>         //now run square-free on h, and return it to the pth power
>>>         zmod_poly_arr_t new_res;
>>>         zmod_poly_square_free(new_res, h);
>>>
>>>         for(unsigned long j=0; j < q ; ++j)
>>>         {
>>>             zmod_poly_arr_cat(res, new_res);    //note, concatenating is
>>> equivalent to multiplying
>>>         }
>>>
>>>     }
>>>     else
>>>     {
>>>         zmod_poly_gcd(g, f, f_d);
>>>         zmod_poly_div(g_1, f, g);
>>>
>>>         //CASE 2:
>>>     in:
>>>         if(DBG) printf("Is g_1 = 1? \n");
>>>         if( zmod_poly_equal(g_1, one))
>>>         {
>>>             goto out;
>>>         }
>>>         else
>>>         {
>>>             if(DBG){ printf("no, g_1 = "); zmod_poly_print(g_1); NL }
>>>             zmod_poly_t h, z;
>>>             zmod_poly_init(h,q);
>>>             zmod_poly_init(z,q);
>>>             //moar!!
>>>             zmod_poly_gcd(h, g_1, g);
>>>             zmod_poly_div(z, g_1, h);
>>>             //out <- out.z
>>>             zmod_poly_arr_add(res, z);
>>>             zmod_poly_set(g_1, h);
>>>             zmod_poly_div(g, g, h);
>>>             if(DBG){ printf("also, g = "); zmod_poly_print(g); NL }
>>>         }
>>>         goto in;
>>>
>>>     out:
>>>         if(DBG) printf("g_1 = 1 \n");
>>>
>>>         if( ! zmod_poly_equal(g,one) )
>>>         {
>>>             if(DBG) printf("g != 1 \n");
>>>             //so now we multiply res with square-free(g^1/p) ^ p
>>>             zmod_poly_arr_t new_res_2;
>>>             zmod_poly_square_free(new_res_2, g);
>>>             zmod_poly_arr_cat(res,new_res_2);
>>>         }
>>>         if(DBG) printf("FINISH! \n");
>>>     }
>>> }
>>>
>>>
>>> --
>>> Richard
>>> -------------------------------------------------------------------------
>>> This SF.Net email is sponsored by the Moblin Your Move Developer's challenge
>>> Build the coolest Linux based applications with Moblin SDK & win great
>>> prizes
>>> Grand prize is a trip for two to an Open Source event anywhere in the world
>>> http://moblin-contest.org/redirect.php?banner_id=100&url=/
>>> _______________________________________________
>>> flint-devel mailing list
>>> fli...@li...
>>> https://lists.sourceforge.net/lists/listinfo/flint-devel
>>>
>>>
>>
>

Re: [flint-devel] Square-Free Algorithm

[Bill H. <goo...@go...>]

I think the second part of the algorithm is OK, except you probably
want zmod_poly_arr_cat(res, z) not zmod_poly_arr_add(res, z).

If I've understood correctly, the way you've implemented this it is
self testing, i.e. the output should return the same thing as the
input. That's quite clever.

Later on we'll set up a struct for containing the factors and we'll
add some functions for catenating factors to such a struct.

One thing I'd like to get rid of is the goto statements. I think what
you have is equivalent to somethink like:

while (g1 == 1)
{
   blah;
}

blah blah;

Bill.

2008/7/15 Bill Hart <goo...@go...>:
> As you worked out, to take f(x)^p, simply raise each term to the power
> of p. That means raising each coefficient to the power p modulo p (use
> Fermat's little theorem for that) and multiplying the exponent of each
> term by p.
>
> When trying to go the other way, i.e. compute f(x)^(1/p),  it's just
> the reverse. I think you got the reverse of the second part OK, but
> not the first.
>
> Bill.
>
> 2008/7/15 Richard Howell-Peak <ric...@gm...>:
>> /**
>>  * Square-Free Algorithm, takes an arbitary polynomial in Fq[X] and returns
>> an array of square free factors
>>  * LOW MULTIPLICITIES
>>  */
>> #define DBG 1
>> void zmod_poly_square_free(zmod_poly_arr_t res, zmod_poly_t f)
>> {
>>     unsigned long q = zmod_poly_modulus(f);    //order of the field
>>     unsigned long deg = zmod_poly_degree(f); //degree of the polynomial
>>     //Let's create some constant polynomials that will be very useful later
>> on:
>>     zmod_poly_t zero, one;
>>     zmod_poly_init(zero,q);
>>     zmod_poly_init(one,q);
>>     zmod_poly_zero(zero);
>>     zmod_poly_set_coeff_ui(one, 0 , 1 );
>>
>>     //set the result to the constant polynomial 1
>>     zmod_poly_arr_init(res,q);    //? Not sure if we should assume it is
>> pre-initialised (NOTE - probably won't)
>>     zmod_poly_arr_add(res,one);
>>
>>     //Step 1, look at f', if it is zero then we are done since f = h(x)^q
>>     //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,q);
>>     zmod_poly_init(f_d,q);
>>     zmod_poly_init(g,q);
>>     zmod_poly_diff(f_d, f);
>>     //CASE 1:
>>     if( zmod_poly_equal(f_d, zero) )
>>     {
>>         if(DBG) printf("f'=0, f = "); zmod_poly_print(f); NL
>>         unsigned long coeff[deg/q];
>>         zmod_poly_t h;
>>         zmod_poly_init(h,q);
>>         for(unsigned long i=0; i <= deg/q; ++i)    //this will be an integer
>> since f'=0
>>         {
>>             if(DBG) printf("Changing new coefficients\n");
>>             coeff[i] = zmod_poly_get_coeff_ui(f, i*q);
>>             //change them slightly
>>             coeff[i] = coeff[i] ^ (i/q);    //clearly this is incorrect,
>> need to do math.pow but for finite field
>>             //now we set up the h poly, where f = h^q
>>             zmod_poly_set_coeff_ui(h, i, coeff[i]);
>>         }
>>         if(DBG) printf("h = ");zmod_poly_print(h); NL
>>         //now run square-free on h, and return it to the pth power
>>         zmod_poly_arr_t new_res;
>>         zmod_poly_square_free(new_res, h);
>>
>>         for(unsigned long j=0; j < q ; ++j)
>>         {
>>             zmod_poly_arr_cat(res, new_res);    //note, concatenating is
>> equivalent to multiplying
>>         }
>>
>>     }
>>     else
>>     {
>>         zmod_poly_gcd(g, f, f_d);
>>         zmod_poly_div(g_1, f, g);
>>
>>         //CASE 2:
>>     in:
>>         if(DBG) printf("Is g_1 = 1? \n");
>>         if( zmod_poly_equal(g_1, one))
>>         {
>>             goto out;
>>         }
>>         else
>>         {
>>             if(DBG){ printf("no, g_1 = "); zmod_poly_print(g_1); NL }
>>             zmod_poly_t h, z;
>>             zmod_poly_init(h,q);
>>             zmod_poly_init(z,q);
>>             //moar!!
>>             zmod_poly_gcd(h, g_1, g);
>>             zmod_poly_div(z, g_1, h);
>>             //out <- out.z
>>             zmod_poly_arr_add(res, z);
>>             zmod_poly_set(g_1, h);
>>             zmod_poly_div(g, g, h);
>>             if(DBG){ printf("also, g = "); zmod_poly_print(g); NL }
>>         }
>>         goto in;
>>
>>     out:
>>         if(DBG) printf("g_1 = 1 \n");
>>
>>         if( ! zmod_poly_equal(g,one) )
>>         {
>>             if(DBG) printf("g != 1 \n");
>>             //so now we multiply res with square-free(g^1/p) ^ p
>>             zmod_poly_arr_t new_res_2;
>>             zmod_poly_square_free(new_res_2, g);
>>             zmod_poly_arr_cat(res,new_res_2);
>>         }
>>         if(DBG) printf("FINISH! \n");
>>     }
>> }
>>
>>
>> --
>> Richard
>> -------------------------------------------------------------------------
>> This SF.Net email is sponsored by the Moblin Your Move Developer's challenge
>> Build the coolest Linux based applications with Moblin SDK & win great
>> prizes
>> Grand prize is a trip for two to an Open Source event anywhere in the world
>> http://moblin-contest.org/redirect.php?banner_id=100&url=/
>> _______________________________________________
>> flint-devel mailing list
>> fli...@li...
>> https://lists.sourceforge.net/lists/listinfo/flint-devel
>>
>>
>

Re: [flint-devel] Square-Free Algorithm

[Bill H. <goo...@go...>]

As you worked out, to take f(x)^p, simply raise each term to the power
of p. That means raising each coefficient to the power p modulo p (use
Fermat's little theorem for that) and multiplying the exponent of each
term by p.

When trying to go the other way, i.e. compute f(x)^(1/p),  it's just
the reverse. I think you got the reverse of the second part OK, but
not the first.

Bill.

2008/7/15 Richard Howell-Peak <ric...@gm...>:
> /**
>  * Square-Free Algorithm, takes an arbitary polynomial in Fq[X] and returns
> an array of square free factors
>  * LOW MULTIPLICITIES
>  */
> #define DBG 1
> void zmod_poly_square_free(zmod_poly_arr_t res, zmod_poly_t f)
> {
>     unsigned long q = zmod_poly_modulus(f);    //order of the field
>     unsigned long deg = zmod_poly_degree(f); //degree of the polynomial
>     //Let's create some constant polynomials that will be very useful later
> on:
>     zmod_poly_t zero, one;
>     zmod_poly_init(zero,q);
>     zmod_poly_init(one,q);
>     zmod_poly_zero(zero);
>     zmod_poly_set_coeff_ui(one, 0 , 1 );
>
>     //set the result to the constant polynomial 1
>     zmod_poly_arr_init(res,q);    //? Not sure if we should assume it is
> pre-initialised (NOTE - probably won't)
>     zmod_poly_arr_add(res,one);
>
>     //Step 1, look at f', if it is zero then we are done since f = h(x)^q
>     //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,q);
>     zmod_poly_init(f_d,q);
>     zmod_poly_init(g,q);
>     zmod_poly_diff(f_d, f);
>     //CASE 1:
>     if( zmod_poly_equal(f_d, zero) )
>     {
>         if(DBG) printf("f'=0, f = "); zmod_poly_print(f); NL
>         unsigned long coeff[deg/q];
>         zmod_poly_t h;
>         zmod_poly_init(h,q);
>         for(unsigned long i=0; i <= deg/q; ++i)    //this will be an integer
> since f'=0
>         {
>             if(DBG) printf("Changing new coefficients\n");
>             coeff[i] = zmod_poly_get_coeff_ui(f, i*q);
>             //change them slightly
>             coeff[i] = coeff[i] ^ (i/q);    //clearly this is incorrect,
> need to do math.pow but for finite field
>             //now we set up the h poly, where f = h^q
>             zmod_poly_set_coeff_ui(h, i, coeff[i]);
>         }
>         if(DBG) printf("h = ");zmod_poly_print(h); NL
>         //now run square-free on h, and return it to the pth power
>         zmod_poly_arr_t new_res;
>         zmod_poly_square_free(new_res, h);
>
>         for(unsigned long j=0; j < q ; ++j)
>         {
>             zmod_poly_arr_cat(res, new_res);    //note, concatenating is
> equivalent to multiplying
>         }
>
>     }
>     else
>     {
>         zmod_poly_gcd(g, f, f_d);
>         zmod_poly_div(g_1, f, g);
>
>         //CASE 2:
>     in:
>         if(DBG) printf("Is g_1 = 1? \n");
>         if( zmod_poly_equal(g_1, one))
>         {
>             goto out;
>         }
>         else
>         {
>             if(DBG){ printf("no, g_1 = "); zmod_poly_print(g_1); NL }
>             zmod_poly_t h, z;
>             zmod_poly_init(h,q);
>             zmod_poly_init(z,q);
>             //moar!!
>             zmod_poly_gcd(h, g_1, g);
>             zmod_poly_div(z, g_1, h);
>             //out <- out.z
>             zmod_poly_arr_add(res, z);
>             zmod_poly_set(g_1, h);
>             zmod_poly_div(g, g, h);
>             if(DBG){ printf("also, g = "); zmod_poly_print(g); NL }
>         }
>         goto in;
>
>     out:
>         if(DBG) printf("g_1 = 1 \n");
>
>         if( ! zmod_poly_equal(g,one) )
>         {
>             if(DBG) printf("g != 1 \n");
>             //so now we multiply res with square-free(g^1/p) ^ p
>             zmod_poly_arr_t new_res_2;
>             zmod_poly_square_free(new_res_2, g);
>             zmod_poly_arr_cat(res,new_res_2);
>         }
>         if(DBG) printf("FINISH! \n");
>     }
> }
>
>
> --
> Richard
> -------------------------------------------------------------------------
> This SF.Net email is sponsored by the Moblin Your Move Developer's challenge
> Build the coolest Linux based applications with Moblin SDK & win great
> prizes
> Grand prize is a trip for two to an Open Source event anywhere in the world
> http://moblin-contest.org/redirect.php?banner_id=100&url=/
> _______________________________________________
> flint-devel mailing list
> fli...@li...
> https://lists.sourceforge.net/lists/listinfo/flint-devel
>
>

Re: [flint-devel] Square-Free Algorithm

[Bill H. <goo...@go...>]

Hi Richard,

Do you also want to send me the code for computing the derivative. Did
you already write that?

Also did you also already write test code for it in zmod_poly-test? If
not, I would do that first. See some of the other examples in
zmod_poly-test. The test code should probably compute the derivative
in a really simple way and compare, e.g. set up an mpz_poly (see the
flint documentation for mpz_poly), multiply each coefficient by the
required multiple, then shift the whole polynomial right by 1, then
reduce all the coefficients mod p.

To convert the coefficients from unsigned longs (which is what they
are in a zmod_poly) to mpz's, just use mpz_set_ui. Once they are
unsigned longs again after you've reduced mod p, you can get their
values as unsigned longs using mpz_get_ui.

Note you can act on coefficients of an mpz_poly with any mpz function
from GMP (see the online GNU MP documentation).

I also noticed you had a zmod_poly_arr_add function. Were you looking
for zmod_add(res, res, one)?

Finally, there is no need to initialise the output poly. One assumes
that the user passes in an initialised poly.

Bill.

2008/7/15 Richard Howell-Peak <ric...@gm...>:
> /**
>  * Square-Free Algorithm, takes an arbitary polynomial in Fq[X] and returns
> an array of square free factors
>  * LOW MULTIPLICITIES
>  */
> #define DBG 1
> void zmod_poly_square_free(zmod_poly_arr_t res, zmod_poly_t f)
> {
>     unsigned long q = zmod_poly_modulus(f);    //order of the field
>     unsigned long deg = zmod_poly_degree(f); //degree of the polynomial
>     //Let's create some constant polynomials that will be very useful later
> on:
>     zmod_poly_t zero, one;
>     zmod_poly_init(zero,q);
>     zmod_poly_init(one,q);
>     zmod_poly_zero(zero);
>     zmod_poly_set_coeff_ui(one, 0 , 1 );
>
>     //set the result to the constant polynomial 1
>     zmod_poly_arr_init(res,q);    //? Not sure if we should assume it is
> pre-initialised (NOTE - probably won't)
>     zmod_poly_arr_add(res,one);
>
>     //Step 1, look at f', if it is zero then we are done since f = h(x)^q
>     //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,q);
>     zmod_poly_init(f_d,q);
>     zmod_poly_init(g,q);
>     zmod_poly_diff(f_d, f);
>     //CASE 1:
>     if( zmod_poly_equal(f_d, zero) )
>     {
>         if(DBG) printf("f'=0, f = "); zmod_poly_print(f); NL
>         unsigned long coeff[deg/q];
>         zmod_poly_t h;
>         zmod_poly_init(h,q);
>         for(unsigned long i=0; i <= deg/q; ++i)    //this will be an integer
> since f'=0
>         {
>             if(DBG) printf("Changing new coefficients\n");
>             coeff[i] = zmod_poly_get_coeff_ui(f, i*q);
>             //change them slightly
>             coeff[i] = coeff[i] ^ (i/q);    //clearly this is incorrect,
> need to do math.pow but for finite field
>             //now we set up the h poly, where f = h^q
>             zmod_poly_set_coeff_ui(h, i, coeff[i]);
>         }
>         if(DBG) printf("h = ");zmod_poly_print(h); NL
>         //now run square-free on h, and return it to the pth power
>         zmod_poly_arr_t new_res;
>         zmod_poly_square_free(new_res, h);
>
>         for(unsigned long j=0; j < q ; ++j)
>         {
>             zmod_poly_arr_cat(res, new_res);    //note, concatenating is
> equivalent to multiplying
>         }
>
>     }
>     else
>     {
>         zmod_poly_gcd(g, f, f_d);
>         zmod_poly_div(g_1, f, g);
>
>         //CASE 2:
>     in:
>         if(DBG) printf("Is g_1 = 1? \n");
>         if( zmod_poly_equal(g_1, one))
>         {
>             goto out;
>         }
>         else
>         {
>             if(DBG){ printf("no, g_1 = "); zmod_poly_print(g_1); NL }
>             zmod_poly_t h, z;
>             zmod_poly_init(h,q);
>             zmod_poly_init(z,q);
>             //moar!!
>             zmod_poly_gcd(h, g_1, g);
>             zmod_poly_div(z, g_1, h);
>             //out <- out.z
>             zmod_poly_arr_add(res, z);
>             zmod_poly_set(g_1, h);
>             zmod_poly_div(g, g, h);
>             if(DBG){ printf("also, g = "); zmod_poly_print(g); NL }
>         }
>         goto in;
>
>     out:
>         if(DBG) printf("g_1 = 1 \n");
>
>         if( ! zmod_poly_equal(g,one) )
>         {
>             if(DBG) printf("g != 1 \n");
>             //so now we multiply res with square-free(g^1/p) ^ p
>             zmod_poly_arr_t new_res_2;
>             zmod_poly_square_free(new_res_2, g);
>             zmod_poly_arr_cat(res,new_res_2);
>         }
>         if(DBG) printf("FINISH! \n");
>     }
> }
>
>
> --
> Richard
> -------------------------------------------------------------------------
> This SF.Net email is sponsored by the Moblin Your Move Developer's challenge
> Build the coolest Linux based applications with Moblin SDK & win great
> prizes
> Grand prize is a trip for two to an Open Source event anywhere in the world
> http://moblin-contest.org/redirect.php?banner_id=100&url=/
> _______________________________________________
> flint-devel mailing list
> fli...@li...
> https://lists.sourceforge.net/lists/listinfo/flint-devel
>
>

[flint-devel] z_intsqrt

[William H. <ha...@ya...>]

I finally got around to implementing a z_intsqrt
function which works for arbitrary sized unsigned
longs, and correctly.

It basically computes floor(sqrt(n)).

I got it to an average of 54 cycles, but was surprised
to find I couldn't get it lower than this without
lookup tables. 

I left the other code I tried, in long_extras.c, but
commented out. It uses a clever algorithm which takes
advantage of the IEEE format to compute an
approximation to 1/sqrt(n), then refines it.

The first approximation is accurate to 4 bits and
subsequent refinements more than double the number of
bits. Probably I only need 4 refinements to be able to
recognize the square root (which is an integer of at
most 32 bits), so long as I check that I am not one
out at the end. 

But it doesn't seem to go any faster than just using
the floating point square root, which is what I
eventually ended up using. 

I think lookup tables could speed up the approximation
method. If I looked up an 8 bit approximation, then
only two refinements would be needed to get the square
root.

I'll leave it for now, since we don't need every last
cycle at this stage, but I'll revisit it some time and
speed it up.

Bill.


      

Re: [flint-devel] Square-Free Algorithm

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

/**
 * Square-Free Algorithm, takes an arbitary polynomial in Fq[X] and returns
an array of square free factors
 * LOW MULTIPLICITIES
 */
#define DBG 1
void zmod_poly_square_free(zmod_poly_arr_t res, zmod_poly_t f)
{
    unsigned long q = zmod_poly_modulus(f);    //order of the field
    unsigned long deg = zmod_poly_degree(f); //degree of the polynomial
    //Let's create some constant polynomials that will be very useful later
on:
    zmod_poly_t zero, one;
    zmod_poly_init(zero,q);
    zmod_poly_init(one,q);
    zmod_poly_zero(zero);
    zmod_poly_set_coeff_ui(one, 0 , 1 );

    //set the result to the constant polynomial 1
    zmod_poly_arr_init(res,q);    //? Not sure if we should assume it is
pre-initialised (NOTE - probably won't)
    zmod_poly_arr_add(res,one);

    //Step 1, look at f', if it is zero then we are done since f = h(x)^q
    //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,q);
    zmod_poly_init(f_d,q);
    zmod_poly_init(g,q);
    zmod_poly_diff(f_d, f);
    //CASE 1:
    if( zmod_poly_equal(f_d, zero) )
    {
        if(DBG) printf("f'=0, f = "); zmod_poly_print(f); NL
        unsigned long coeff[deg/q];
        zmod_poly_t h;
        zmod_poly_init(h,q);
        for(unsigned long i=0; i <= deg/q; ++i)    //this will be an integer
since f'=0
        {
            if(DBG) printf("Changing new coefficients\n");
            coeff[i] = zmod_poly_get_coeff_ui(f, i*q);
            //change them slightly
            coeff[i] = coeff[i] ^ (i/q);    //clearly this is incorrect,
need to do math.pow but for finite field
            //now we set up the h poly, where f = h^q
            zmod_poly_set_coeff_ui(h, i, coeff[i]);
        }
        if(DBG) printf("h = ");zmod_poly_print(h); NL
        //now run square-free on h, and return it to the pth power
        zmod_poly_arr_t new_res;
        zmod_poly_square_free(new_res, h);

        for(unsigned long j=0; j < q ; ++j)
        {
            zmod_poly_arr_cat(res, new_res);    //note, concatenating is
equivalent to multiplying
        }

    }
    else
    {
        zmod_poly_gcd(g, f, f_d);
        zmod_poly_div(g_1, f, g);

        //CASE 2:
    in:
        if(DBG) printf("Is g_1 = 1? \n");
        if( zmod_poly_equal(g_1, one))
        {
            goto out;
        }
        else
        {
            if(DBG){ printf("no, g_1 = "); zmod_poly_print(g_1); NL }
            zmod_poly_t h, z;
            zmod_poly_init(h,q);
            zmod_poly_init(z,q);
            //moar!!
            zmod_poly_gcd(h, g_1, g);
            zmod_poly_div(z, g_1, h);
            //out <- out.z
            zmod_poly_arr_add(res, z);
            zmod_poly_set(g_1, h);
            zmod_poly_div(g, g, h);
            if(DBG){ printf("also, g = "); zmod_poly_print(g); NL }
        }
        goto in;

    out:
        if(DBG) printf("g_1 = 1 \n");

        if( ! zmod_poly_equal(g,one) )
        {
            if(DBG) printf("g != 1 \n");
            //so now we multiply res with square-free(g^1/p) ^ p
            zmod_poly_arr_t new_res_2;
            zmod_poly_square_free(new_res_2, g);
            zmod_poly_arr_cat(res,new_res_2);
        }
        if(DBG) printf("FINISH! \n");
    }
}


-- 
Richard

Re: [flint-devel] Square-Free Algorithm

[Bill H. <goo...@go...>]

Do you want to attach your code and I'll take a look.

Bill.

2008/7/15 Richard Howell-Peak <ric...@gm...>:
> In trying to implement the square-free algorithm I have run into a couple of
> troubles.
>
> In the first case you look to see if f'=0, which implies f = h^p where you
> compute the coefficients of h by raising them to powers i/p but I have no
> idea how to compute such powers in Z/pZ. I'm not sure if this is causing the
> other problem in which you end up calling recursively a constant polynomial
> forever until the stack overflows but this doesn't always happen.
>
> The other part of the algorithm does not work either but I haven't had
> sufficient time to go through it and see why.
>
> --
> Richard
> -------------------------------------------------------------------------
> This SF.Net email is sponsored by the Moblin Your Move Developer's challenge
> Build the coolest Linux based applications with Moblin SDK & win great
> prizes
> Grand prize is a trip for two to an Open Source event anywhere in the world
> http://moblin-contest.org/redirect.php?banner_id=100&url=/
> _______________________________________________
> flint-devel mailing list
> fli...@li...
> https://lists.sourceforge.net/lists/listinfo/flint-devel
>
>

Re: [flint-devel] OpenMP and latest build

[Bill H. <goo...@go...>]

Oh yeah, you are right, we need to install gcc 4.3.

You might be able to install gcc 4.3 in your home directory, but I
guess I should do it globally some time. I'm probably not going to get
to it today. Can you survive without it for a few days?

Bill.

On 15/07/2008, Peter Shrimpton <ps...@gm...> wrote:
> Hey.
>
> It appears that OpenMP is currently not supported on the main server. I have
> tried the following code:
>
> .
> .
> .
> #include <omp.h>
> int main (int argc, char const *argv[])
>  {
>     int i,array[20];
>     unsigned long j;
>
>     for (i=0; i<20; i++){
>         array[i]=i%7;
>     }
>
> #ifdef _OPENMP
>     printf("WE HAVE MP\n");
> #endif
>
> #pragma omp parallel private(i,j) num_threads(4)
>      {
> #pragma omp for  schedule(dynamic,1) nowait
>         for (i=0; i<20; i++){
>             if (array[i]==0){
>                 printf("%d = 0 mod 7\n",i);
>                 for (j=0; j<10000000000; j++){
>                  }
>             }
>             //printf("hello\n");
>         }
>     }
>
> }
>
>
> gcc says it can't find omp.h, when I coment out that line it dosn't print
> "WE HAVE MP" and the for loop doesn't appear to be executed in parallel. As
> openMP is a set of compiler pragma's I have no idea how to install it or if
> I would have permission. I think we may need gcc 4.2 but we currently only
> have gcc 4.1.2.
>
> Also after checking out the latest build I have suddenly come accross
> "undefined reference to `z_intsqrt'" errors. I have tried to make the
> library but it dosn't compile.
>
> Peter
>
>
>
> -------------------------------------------------------------------------
>  This SF.Net email is sponsored by the Moblin Your Move Developer's
> challenge
>  Build the coolest Linux based applications with Moblin SDK & win great
> prizes
>  Grand prize is a trip for two to an Open Source event anywhere in the world
> http://moblin-contest.org/redirect.php?banner_id=100&url=/
> _______________________________________________
>  flint-devel mailing list
>  fli...@li...
>  https://lists.sourceforge.net/lists/listinfo/flint-devel
>
>

[flint-devel] Square-Free Algorithm

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

In trying to implement the square-free algorithm I have run into a couple of
troubles.

In the first case you look to see if f'=0, which implies f = h^p where you
compute the coefficients of h by raising them to powers i/p but I have no
idea how to compute such powers in Z/pZ. I'm not sure if this is causing the
other problem in which you end up calling recursively a constant polynomial
forever until the stack overflows but this doesn't always happen.

The other part of the algorithm does not work either but I haven't had
sufficient time to go through it and see why.

-- 
Richard