flint-devel — Development list for the FLINT project

Re: [flint-devel] Function Interface Updates

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

2008/8/5  <D.C...@wa...>:
> Hi,
>
>> Wow, that's detailed. Whoever goes to implement this will find it quite
>>easy.
>>
>> One thing that worries me is we decided to represent a p-adic as a
>>multiprecision integer, but if I am reading what you've written
>> correctly, you are thinking of a padic as a power series in number 6
>>below. Or am I just misunderstanding what you mean?
>
> I think you may be misunderstanding me, all the arithmetic I've planned is
> by representing a p-adic as an integer. My reasoning for (6) was in case a
> user wished to input a p-adic as you would by hand, say 3^-1 + 1 + 2*3 +
> 3^2 + O(3^3) (or in the notation below; 1 . 1 2 1) instead of inputting it
> more directly as an integer. This function would allow us to convert the
> input into a format that we would do arithmetic with. All the arithmetic
> I've planned is done with integers that represent the p-adic, theres no
> arithmetic planned using the power series representation.

OK, I was just being stupid. Yes I see what you mean now by number
(6). That's definitely a good idea to have such a function.

>
> I was thinking as to what would be the best way to deal with say 16/5
> which lies in Z_3. I was proposing that we store two integers to represent
> such a p-adic as (16,5) instead of storing (16/5)%(3^n); However will
> storing (16/5)%3^n be faster than storing two integers? The actual
> manipulation of the integers will take longer if we store two integers. We
> have to check if a power of p divides the denominator straight away
> otherwise the inverse modulo p^n wont exist. For example if we wanted
> 16/15 +O(3^5), we would have to take out 3^-1 straight away to get 16/5
> and then find 16/5%3^6. We need to check the numerator for powers of p
> also since if we perform a division we could end up with a power of p on
> the denominator, in which case the reductions modulo p^k will be
> impossible. This was one of the things I was trying to take care of in my
> last email outlining division.

Fractions are interesting because many algorithms work over Q. My
initial though is that there are going to be a whole load of extra
functions to do arithmetic with p-adic fractions if we store both a
numerator and denominator. That sounds like a lot of extra work.

But then again, as you discovered with those papers on hensel codes,
it may actually be possible to write an entire module for Q where
elements of Q are actually stored as p-adics. Of course that doesn't
mean we can't just use the original p-adic format that we came up
with.

I'm inclined to say it makes more sense to just use the original
format, i.e. don't store a separate denominator, just store the
"shift" (e.g. the 3^-1 as per the example you gave).

I guess if you store a fraction as a quotient of two p-adic integers,
you have to check the numerator and denominator for divisibility by p
anyway. So it isn't really going to cost you extra to have to do this
check if we want to use our {Zp, shift} model.

>
> In my previous email I was working on the basis that we store two integers
> for each p-adic, hence the a,b,c,d, but things will be simpler if we store
> a single integer. I'll have a look at this, re-write it and send it again
> later.
>

Yeah I don't think I see any reason to store two integers for each
p-adic fraction. If you do come across a reason, we can revisit it. I
remember when we were thinking about polynomials the same issue came
up, do you store a fraction for each coefficient, or do you just use
integers and a single common denominator. The big cost with fractions
is always having to do gcd's to figure out common denominators. I
remember there was a big discussion at one of the sage conferences
about this, and again online.

In the end, the single common denominator model won out, mainly due to
John Cremona saying it should work out.

>
>> By the way, I think you asked a few days ago about how to represent
>>fractions in Qp. I think it will be enough to represent them as an
>>element of Zp multiplied by the appropriate power of p (where all we
> store is >the exponent of the power), e.g:
>
> Yes, all the information about the structure of Q_p we need is stored in
> Z_p so we can reduce arithmetic of Q_p to Z_p in this way which is how
> I've been thinking of doing things. So for example if we had 17/6 +
> O(3^5), we would write this as (3^-1)*(17/2) + O(3^6).

Yep.

>
>
>
>> 2*3^-2+3^-1+1+2+O(3^2) would be represented as:
>>
>> -2, 2+3+3^2+2*3^3+O(3^4) which will really get stored as the quadruple
> of numbers:
>>
>> -2, 68, 3, 4 i.e. /3^-2)*(68 mod 3^4)
>>
>> Does that look like it works?
>>
>> Bill.
>>
>
>
> If x = 2*3^-2 +3^-1 + 1 +2*3 + O(3^2)
> then I would write x*3^2 = 2 + 3^1 + 3^2 + 2*3^3 + O(3^4)
>
> So we would need to store -2, 68, 3 and 4 as you say.
>

Great.

Bill.

Re: [flint-devel] Function Interface Updates

[<D.C...@wa...>]

Hi,

> Wow, that's detailed. Whoever goes to implement this will find it quite
>easy.
>
> One thing that worries me is we decided to represent a p-adic as a
>multiprecision integer, but if I am reading what you've written
> correctly, you are thinking of a padic as a power series in number 6
>below. Or am I just misunderstanding what you mean?

I think you may be misunderstanding me, all the arithmetic I've planned is
by representing a p-adic as an integer. My reasoning for (6) was in case a
user wished to input a p-adic as you would by hand, say 3^-1 + 1 + 2*3 +
3^2 + O(3^3) (or in the notation below; 1 . 1 2 1) instead of inputting it
more directly as an integer. This function would allow us to convert the
input into a format that we would do arithmetic with. All the arithmetic
I've planned is done with integers that represent the p-adic, theres no
arithmetic planned using the power series representation.

I was thinking as to what would be the best way to deal with say 16/5
which lies in Z_3. I was proposing that we store two integers to represent
such a p-adic as (16,5) instead of storing (16/5)%(3^n); However will
storing (16/5)%3^n be faster than storing two integers? The actual
manipulation of the integers will take longer if we store two integers. We
have to check if a power of p divides the denominator straight away
otherwise the inverse modulo p^n wont exist. For example if we wanted
16/15 +O(3^5), we would have to take out 3^-1 straight away to get 16/5
and then find 16/5%3^6. We need to check the numerator for powers of p
also since if we perform a division we could end up with a power of p on
the denominator, in which case the reductions modulo p^k will be
impossible. This was one of the things I was trying to take care of in my
last email outlining division.

In my previous email I was working on the basis that we store two integers
for each p-adic, hence the a,b,c,d, but things will be simpler if we store
a single integer. I'll have a look at this, re-write it and send it again
later.


> By the way, I think you asked a few days ago about how to represent
>fractions in Qp. I think it will be enough to represent them as an
>element of Zp multiplied by the appropriate power of p (where all we
store is >the exponent of the power), e.g:

Yes, all the information about the structure of Q_p we need is stored in
Z_p so we can reduce arithmetic of Q_p to Z_p in this way which is how
I've been thinking of doing things. So for example if we had 17/6 +
O(3^5), we would write this as (3^-1)*(17/2) + O(3^6).



> 2*3^-2+3^-1+1+2+O(3^2) would be represented as:
>
> -2, 2+3+3^2+2*3^3+O(3^4) which will really get stored as the quadruple
of numbers:
>
> -2, 68, 3, 4 i.e. /3^-2)*(68 mod 3^4)
>
> Does that look like it works?
>
> Bill.
>


If x = 2*3^-2 +3^-1 + 1 +2*3 + O(3^2)
then I would write x*3^2 = 2 + 3^1 + 3^2 + 2*3^3 + O(3^4)

So we would need to store -2, 68, 3 and 4 as you say.


Daniel Ellam

[flint-devel] Polynomial factorisation

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

I've been working this week on the code Richard has
been writing for polynomial factorisation.

By the way, has anyone seen Richard? I haven't seen
him for about a week!!

So far we've implemented the square-free factorisation
for factorising polynomials over Z/pZ into polynomials
without square factors, and Berlekamp's algorithm for
factoring polynomials over Z/pZ into irreducible
polynomials, which takes as input... you guessed it,
polynomials with no square factors.

Amazingly for polynomials of length up to about 40
with word sized primes (up to 64 bits), we are already
2-3 times faster than Magma!! 

Surely I've made a mistake in the timing somewhere
because we haven't even really optimised all that well
yet. But I've checked for some time now to find
something wrong in the timings, and I think they are
correct.

It wasn't the easiest thing to do as many of the
references on Berlekamp's algorithm are very poor
indeed. The main reference Richard and I had been
working from was completely unclear and contained
errors, though in the end, the algorithm contained
therein actually seems very good.

Overall an excellent result for 5 weeks work.

Bill.


      

Re: [flint-devel] Simultaneous Sieving and Trial factoring or Pocklington

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

Holy moly!! That was quick work. I thought it would take you at least
a week to write that code. That's amazing.

I only just figured out the best hash function to use last night, and
it is precisely the one you chose. Well done!

The idea I had doesn't work I think. It was basically this:

Suppose we want 2^K mod n1, n2, ..., nr

where n1, n2, ..., nr are all relatively close together (and so in our
case coprime).

My thought was to first compute A = 2^K mod n1*n2*....*nr using
multiprecision arithmetic.

Now we get 2^K mod n1*n2*...n{r/2} and 2^K mod n{r/2+1}*...*nr

simply by reducing A modulo both of the smaller products.

We keep splitting the products into two (recusively), until we are
finally getting 2^K mod n1, n2, ...., nr.

This definitely won't work with ordinary GMP division however, since
it uses an algorithm which is asymptotically too slow (O(M(n)log(n))
where M(n) is the time taken to do a multiplication. However, FLINT
has some code for Montgomery multiplication in fmpz.c/h which gets rid
of that extra log factor in precisely this situation. I can't recall
if we already implemented reduction modulo a whole pile of n's using
montgomery multiplication yet or not. If it is there, it is exactly
what you need. If not, I have to implement it anyway.

But, I'm still not convinced my idea works anyway. We could try it I
suppose. But I think I convinced myself yesterday it will actually be
slower than computing the individual 2^K mod ni's

Bill.

2008/8/5 Peter Shrimpton <ps...@gm...>:
> Hello Bill
>
> I have written some code to sieve and factor simultaneously
>
> sievefact.c sieves a chunk of numbers and uses the trial devision to factor
> n-1 for all n that get through the sieve. It uses a hash table to reduce the
> storage space required with the hash function being the 2nd to last 4 bits.
> (The last bit is alway the same). This seams to work very well indeed.
>
> sievfact2.c simultaneously sieves and checks to see if it can find a b such
> that gcd(b^((n-1)/p)-1,n)=1. It uses the same code as sievefact.c but it
> doesn't store the factors, only the cofactor. This could be useful if the
> Pocklington test is quicker then the pseudoprime tests the only trouble
> being to then decide which numbers to do the pseudoprime tests on to find a
> counter example to the Baillie-PSW test.
>
> Would it be possible for you to explain your idea about calculating powers
> of 2 mod N for lots of different numbers at once again? I can't quite
> remember exactly how it was to be done or how we could go about splitting
> the powers up once we had calculated the power of the product. An
> explanation in an email should suffice.
>
> I have also written a program that uses the list I emailed you before to
> look for counter examples to the Baillie-PSW test below 2^63. Unfortunately
> there weren't that many products in the list that were less then2^63 and I
> didn't find any counter examples. Maybe it could be worth implementing a
> Baillie-PSW test in fmpz's and having another go?
>
>
> 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
>
>

Re: [flint-devel] Function Interface Updates

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

Wow, that's detailed. Whoever goes to implement this will find it quite easy.

One thing that worries me is we decided to represent a p-adic as a
multiprecision integer, but if I am reading what you've written
correctly, you are thinking of a padic as a power series in number 6
below. Or am I just misunderstanding what you mean?

By the way, I think you asked a few days ago about how to represent
fractions in Qp. I think it will be enough to represent them as an
element of Zp multiplied by the appropriate power of p (where all we
store is the exponent of the power), e.g:

2*3^-2+3^-1+1+2+O(3^2) would be represented as:

-2, 2+3+3^2+2*3^3+O(3^4) which will really get stored as the quadruple
of numbers:

-2, 68, 3, 4 i.e. /3^-2)*(68 mod 3^4)

Does that look like it works?

Bill.

2008/8/5  <D.C...@wa...>:
> I give an outline of p-adic division. I've used some functions from long
> extras.
>
>
>
> --------------------------------------------------------------------------------------
> Functions
>
> (1) unsigned long z_powmod(unsigned long b, long exp, unsigned long n).
>
> (2) unsigned long z_gcd(long a, long b)
>
> (3) unsigned long z_pow(unsigned long a, unsigned long exp)
>
> (4) unsigned long mulmod2_precomp(unsigned long a, unsigned long b,
> unsigned long n, pre_inv2_t ninv)
>
> (5) A function for converting Q -> Q_p. Something like num_to_padic_conv.
> To compute the coefficients, repetedly use (1) to solve congruences modulo
> p^k for a and (1/b) and then a*(1/b). Store the elements in an array as
> they are computed and output something like 1 . 1 0 1 to represent the
> p-adic P^-1 + 1 + p^2.
>
> (6) A function to convert Q_p -> Q. Again something like
> padic_to_num_conv. We want to take an array of coefficients like 1 0 2 . 0
> 1 1 representing p^-3 + 2*p^-1 + p + p^2 and convert it to (a,b). I think
> in this case we should set a = 1 + 2*p^2 + p^4 + p^5 and b = p^3. More
> generally, consider the non-negative coefficient with the largest negative
> exponent of p, say -m and set b = pow(p,m) by using (3) and once we've
> essentially taken p^-m out as a factor we can sum the terms to find a.
>
> (7) unsigned long z_div2_precomp(unsigned long a, unsigned long n,
> pre_inv2_t ninv)
>
> (8) A function to find the number of times p divides b. Can we use
> unsigned long
> z_powmod(unsigned long b, long exp, unsigned long n) for this?
>
> (9) Function for reducing precision. This will be needed if we have two
> p-adics of different precision. Something like padic_prec_redu. Simply
> taking the p-adic modp^k for some suitable k should do.
>
> (10) A function to "shift" the p-adic, so for example adding or taking out
> a factor of p from the p-adic.
>
> --------------------------------------------------------------------------------------
>
>
> * Use (6) if necessary to convert a p-adic into an element of Q. So I can
> assume the p-adic is inputted in the form a/b together with a prime p and
> a precision n, say something like padic_num_t (a,b). Since we are likely
> going to need modulo arithmetic, we probably want to be working with the
> long extras module as opposed to fmpz when we deal with the two intgers a
> and b representing numerator and denominator. This also allows us to use
> unsigned long z_invert in a function for converting our result a/b back to
> a p-adic. I do realise that this limits us to unsigned longs.
>
>
> p-adic division.  We need a precision, prime and the input which will be
> of the form (a,b) representing a/b.
>
> * Assume we are given p-adics (a,b) and (c,d) with precision n and prime
> p. a may be signed, but we can assume b isnt signed.
>
> * The long extras module can deal with gcd, so given ints a and b, first
> find gcd(a,b) = D. Use (2) for this.
>        -> Set a = a/D, b = b/D (6).
>        -> Use (8) to find the p-adic valuation of a and b. If v_p(b) = N
>>1 where v_p denotes the p-adic valuation of b, then we know that
> a/b lies in Q_p, so letting b' = b/pow(p,N)  by using (3) to find
> pow(p, N) and (6) for the division. Similarly we want to take
> powers of p out of a (if any exist after the previous step) to
> avoid impossible divisions modulo p^k later. Let v_p(a) = M. By
> the last step, only one of M or N can be greater than 0. Let the
> resulting number be a'. We want to store a', b', M-N, p and n,
> then a'/b' lies in Z_p and the arithmetic in Q_p can be reduced to
> arithmetic in Z_p. Note that we are basically rewriting the p-adic
> a/b as p^(M-N)a'/b' and we shall perform the arithmetic on a'/b'.
>        -> Then we want something like padic_num_init (a', b', M-N, p, n).
>
> * Repeat the same process for (c, d).
>
> * The p-adic arithmetic is performed on the integers a', b', c', d' (where
> d' = d/pow(p,v_p(d) and similarly for c'.)
>
>        -> We want to find (a'/b')/(c'/d') = (a'*d')/(b'*c'). We don't
> want to reduce a'*d' and b'*c' modp^n after each step in say a while loop,
> it will be quicker to do one reduction mod p^n at the end rather than a
> reduction at each step. The end reduction (or if no loop is being used)
> can use (4) where we take the multiplications a'*d' and b'*c' modulo p^n.
> This will make converting a'*d'/b'*c' to its p-adic expansion quicker. We
> also cancel the powers of p we took out of a/b and c/d at this point. This
> requires an addition.
>
> The final step is converting (a'*d',b'*c') back into its p-adic expansion.
> We want to put the powers of p we took out of a, b, c and d (after
> cancellation has taken place; we have a/b = p^(M-N)a'/b' and something
> similar for c/d and when we perform the division we want to cancel the
> powers of p appropriately) back in at the last step, i.e. "shift" the
> resulting p-adic depending on what powers of p are left, (10). Use (5)
> (although we will want modp^(n+(N-M)) since we took factors of powers of p
> from a and b and so shifted the p-adic) and then perform a final shift as
> just described.
>
>
>
> Daniel Ellam
>
>
>
>
>>
>
>
> -------------------------------------------------------------------------
> 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] Simultaneous Sieving and Trial factoring or Pocklington

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

Hello Bill

I have written some code to sieve and factor simultaneously

sievefact.c sieves a chunk of numbers and uses the trial devision to factor
n-1 for all n that get through the sieve. It uses a hash table to reduce the
storage space required with the hash function being the 2nd to last 4 bits.
(The last bit is alway the same). This seams to work very well indeed.

sievfact2.c simultaneously sieves and checks to see if it can find a b such
that gcd(b^((n-1)/p)-1,n)=1. It uses the same code as sievefact.c but it
doesn't store the factors, only the cofactor. This could be useful if the
Pocklington test is quicker then the pseudoprime tests the only trouble
being to then decide which numbers to do the pseudoprime tests on to find a
counter example to the Baillie-PSW test.

Would it be possible for you to explain your idea about calculating powers
of 2 mod N for lots of different numbers at once again? I can't quite
remember exactly how it was to be done or how we could go about splitting
the powers up once we had calculated the power of the product. An
explanation in an email should suffice.

I have also written a program that uses the list I emailed you before to
look for counter examples to the Baillie-PSW test below 2^63. Unfortunately
there weren't that many products in the list that were less then2^63 and I
didn't find any counter examples. Maybe it could be worth implementing a
Baillie-PSW test in fmpz's and having another go?


Peter

Re: [flint-devel] Function Interface Updates

[<D.C...@wa...>]

I give an outline of p-adic division. I've used some functions from long
extras.



--------------------------------------------------------------------------------------
Functions

(1) unsigned long z_powmod(unsigned long b, long exp, unsigned long n).

(2) unsigned long z_gcd(long a, long b)

(3) unsigned long z_pow(unsigned long a, unsigned long exp)

(4) unsigned long mulmod2_precomp(unsigned long a, unsigned long b,
unsigned long n, pre_inv2_t ninv)

(5) A function for converting Q -> Q_p. Something like num_to_padic_conv.
To compute the coefficients, repetedly use (1) to solve congruences modulo
p^k for a and (1/b) and then a*(1/b). Store the elements in an array as
they are computed and output something like 1 . 1 0 1 to represent the
p-adic P^-1 + 1 + p^2.

(6) A function to convert Q_p -> Q. Again something like
padic_to_num_conv. We want to take an array of coefficients like 1 0 2 . 0
1 1 representing p^-3 + 2*p^-1 + p + p^2 and convert it to (a,b). I think
in this case we should set a = 1 + 2*p^2 + p^4 + p^5 and b = p^3. More
generally, consider the non-negative coefficient with the largest negative
exponent of p, say -m and set b = pow(p,m) by using (3) and once we've
essentially taken p^-m out as a factor we can sum the terms to find a.

(7) unsigned long z_div2_precomp(unsigned long a, unsigned long n,
pre_inv2_t ninv)

(8) A function to find the number of times p divides b. Can we use
unsigned long
z_powmod(unsigned long b, long exp, unsigned long n) for this?

(9) Function for reducing precision. This will be needed if we have two
p-adics of different precision. Something like padic_prec_redu. Simply
taking the p-adic modp^k for some suitable k should do.

(10) A function to "shift" the p-adic, so for example adding or taking out
a factor of p from the p-adic.

--------------------------------------------------------------------------------------


* Use (6) if necessary to convert a p-adic into an element of Q. So I can
assume the p-adic is inputted in the form a/b together with a prime p and
a precision n, say something like padic_num_t (a,b). Since we are likely
going to need modulo arithmetic, we probably want to be working with the
long extras module as opposed to fmpz when we deal with the two intgers a
and b representing numerator and denominator. This also allows us to use 
unsigned long z_invert in a function for converting our result a/b back to
a p-adic. I do realise that this limits us to unsigned longs.


p-adic division.  We need a precision, prime and the input which will be
of the form (a,b) representing a/b.

* Assume we are given p-adics (a,b) and (c,d) with precision n and prime
p. a may be signed, but we can assume b isnt signed.

* The long extras module can deal with gcd, so given ints a and b, first
find gcd(a,b) = D. Use (2) for this.
        -> Set a = a/D, b = b/D (6).
        -> Use (8) to find the p-adic valuation of a and b. If v_p(b) = N
>1 where v_p denotes the p-adic valuation of b, then we know that
a/b lies in Q_p, so letting b' = b/pow(p,N)  by using (3) to find
pow(p, N) and (6) for the division. Similarly we want to take
powers of p out of a (if any exist after the previous step) to
avoid impossible divisions modulo p^k later. Let v_p(a) = M. By
the last step, only one of M or N can be greater than 0. Let the
resulting number be a'. We want to store a', b', M-N, p and n,
then a'/b' lies in Z_p and the arithmetic in Q_p can be reduced to
arithmetic in Z_p. Note that we are basically rewriting the p-adic
a/b as p^(M-N)a'/b' and we shall perform the arithmetic on a'/b'.
        -> Then we want something like padic_num_init (a', b', M-N, p, n).

* Repeat the same process for (c, d).

* The p-adic arithmetic is performed on the integers a', b', c', d' (where
d' = d/pow(p,v_p(d) and similarly for c'.)

        -> We want to find (a'/b')/(c'/d') = (a'*d')/(b'*c'). We don't
want to reduce a'*d' and b'*c' modp^n after each step in say a while loop,
it will be quicker to do one reduction mod p^n at the end rather than a
reduction at each step. The end reduction (or if no loop is being used)
can use (4) where we take the multiplications a'*d' and b'*c' modulo p^n.
This will make converting a'*d'/b'*c' to its p-adic expansion quicker. We
also cancel the powers of p we took out of a/b and c/d at this point. This
requires an addition.

The final step is converting (a'*d',b'*c') back into its p-adic expansion.
We want to put the powers of p we took out of a, b, c and d (after
cancellation has taken place; we have a/b = p^(M-N)a'/b' and something
similar for c/d and when we perform the division we want to cancel the
powers of p appropriately) back in at the last step, i.e. "shift" the
resulting p-adic depending on what powers of p are left, (10). Use (5)
(although we will want modp^(n+(N-M)) since we took factors of powers of p
from a and b and so shifted the p-adic) and then perform a final shift as
just described.



Daniel Ellam




>

Re: [flint-devel] Fixed base exponention

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

That's interesting. Now the question is, can we adapt this to compute
powers mod n?

Bill.

2008/8/4 Peter Shrimpton <ps...@gm...>:
> Hello
>
> I was looking at powering algorithms and found a couple that only work for
> fixed bases but could be quicker as they use precomputed numbers.
>
> The first is in "Crandall and Pomerance: Prime Numbers a Computational
> Perspective". The psudo-code is:
>
> Computes x^y
> y is written in base B as y_0 y_1 .... y_(D-1) with y_(D-1)>0 the high
> digit.
> We assume that (B-1)*(D-1) numbers have been precalculated:
> x^(i*(B^j))   (x to the power of (i times (B to the power of j))    for i=1
> to B-1, j=1 to D-1
>
> Initialize:
>    z=1;
> Loop:
>    for (0<=j<D){
>       z=z*x^(y_j*(B^j))
>    }
>    return z;
>
> The other algorithm can be found near the bottom of page 15 on:
>
> www.daimi.au.dk/~ivan/FastExpproject.pdf
>
> According to the author the fixed base algorithm is 3 times faster then the
> general purpose one. Not much of an improvement but in the conclusion the
> author says he/she was using a programing language that was really slow at
> accessing arrays so we may be able to get more out of 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] Fixed base exponention

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

Hello

I was looking at powering algorithms and found a couple that only work for
fixed bases but could be quicker as they use precomputed numbers.

The first is in "Crandall and Pomerance: Prime Numbers a Computational
Perspective". The psudo-code is:

Computes x^y
y is written in base B as y_0 y_1 .... y_(D-1) with y_(D-1)>0 the high
digit.
We assume that (B-1)*(D-1) numbers have been precalculated:
x^(i*(B^j))   (x to the power of (i times (B to the power of j))    for i=1
to B-1, j=1 to D-1

Initialize:
   z=1;
Loop:
   for (0<=j<D){
      z=z*x^(y_j*(B^j))
   }
   return z;

The other algorithm can be found near the bottom of page 15 on:

www.daimi.au.dk/~ivan/FastExpproject.pdf

According to the author the fixed base algorithm is 3 times faster then the
general purpose one. Not much of an improvement but in the conclusion the
author says he/she was using a programing language that was really slow at
accessing arrays so we may be able to get more out of it

Peter

Re: [flint-devel] Function Interface Updates

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

Mate this is waaaay too vague and way too brief. Are these elements of Zp or Qp?

How will division be handled.

What about square roots?

What about functions to return the various fields of the padic
structure, to reduce the precision?

What about raising to a power, nth roots, etc etc etc.

You haven't even scratched the surface and I see no imdication of what
algoritms can be used.

I hope there is more to come, as this is the most important part of
the project. The outcome of all that reading you did.

Bill.
On 04/08/2008, D.C...@wa... <D.C...@wa...> wrote:
> Hi,
>
> This is my list so far;
>
> We need a p-adic data type, something like padic_num_t which I suppose
> will store things like the prime p, the precision and the integer. Then, I
> imagine a list of the basic functions will look something like:
>
> void  padic_num_add(padic_num_t   output, padic_num_t   padic1,
> padic_num_t   padic2)
> void  padic_num_sub(padic_num_t   output, padic_num_t   padic1,
> padic_num_t   padic2)
> void  padic_num_mul(padic_num_t   output, padic_num_t   padic1,
> padic_num_t   padic2)
> void  padic_num_div(padic_num_t   output, padic_num_t   padic1,
> padic_num_t   padic2)
>
> Then p-adic data type for polynomials over the p-adics would be something
> like padic_poly_t.
>
> I make this the 6th week of our projects and so I'm scheduled to leave
> campus at the end of this week. When I'm home I'll write the report and
> can tie up any loose ends by email. If I need to stay on campus longer it
> can easily be arranged.
>
> Thanks,
>
> Daniel Ellam
>
>
>
>>
>>
>> I've been considering the basic operations +,- ,* ,/  for the function
>> interface and converting between p-adic and integer/rational
>> representations.
>>
>
>
> -------------------------------------------------------------------------
> 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] Function Interface Updates

[<D.C...@wa...>]

Hi,

This is my list so far;

We need a p-adic data type, something like padic_num_t which I suppose
will store things like the prime p, the precision and the integer. Then, I
imagine a list of the basic functions will look something like:

void  padic_num_add(padic_num_t   output, padic_num_t   padic1,
padic_num_t   padic2)
void  padic_num_sub(padic_num_t   output, padic_num_t   padic1,
padic_num_t   padic2)
void  padic_num_mul(padic_num_t   output, padic_num_t   padic1,
padic_num_t   padic2)
void  padic_num_div(padic_num_t   output, padic_num_t   padic1,
padic_num_t   padic2)

Then p-adic data type for polynomials over the p-adics would be something
like padic_poly_t.

I make this the 6th week of our projects and so I'm scheduled to leave
campus at the end of this week. When I'm home I'll write the report and
can tie up any loose ends by email. If I need to stay on campus longer it
can easily be arranged.

Thanks,

Daniel Ellam



>
>
> I've been considering the basic operations +,- ,* ,/  for the function
> interface and converting between p-adic and integer/rational
> representations.
>

[flint-devel] host-57-44 back up

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

I've now reinstalled host-57-44 with KUbuntu. The only problem I had
so far is the KDE package manager crashes fairly often. But it's just
as easy to use apt-get anyway.

I've reinstalled mpir, mpfr, par/gp and am in the process of
reinstalling gcc-4.3.1 and magma.

It's definitely 64 bits now.

If I gave you a password a few days ago, it should be the same one
now. If not, I am about to email one to you.

Peter Shrimpton noted that to use gcc-4.3.1 you additionally need to
put mpfr in the LD_LIBRARY_PATH. It is currently installed in
/usr/lib, i.e. you need to do

export PATH=/usr/local/gcc-4.3.1/bin:$PATH
export LD_LIBRARY_PATH=/usr/local/mpir-0.9.0/:/usr/lib/:$LD_LIBRARY_PATH

in order to use the new gcc.

Bill.

Re: [flint-devel] host-57-44 reinstalled

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

Hi all,

I installed Kubuntu x86_64 on my laptop last night, and it is really
excellent!! So this morning I am going to have a go at upgrading the
server host-57-44 and see if it takes it. This means the server will
be inaccessible from about 9pm onwards for a few hours.

I'll back up anything that has been put into home directories in the
past few days to another machine or to disk so that we don't lose
anything this time!

Bill.

2008/7/29 William Hart <ha...@ya...>:
> Hi all,
>
> the SUSE linux I used installed a 32 bit kernel on
> host-57-44, which is a 64 bit Opteron machine.
>
> Thus it needs to be reinstalled again.
>
> -1 (again) for the SUSE developers for not telling me
> this was a 32 bit kernel only.
>
> Bill.
>
>
>
>
> -------------------------------------------------------------------------
> 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] Parallel Stuff

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

I should be able to upgrade to a 64 bit kernel tomorrow. That will
also speed things up slightly.

It's not clear to me that we actually need to use threads. We could
just run the same program 4 times starting at different points each
time. So I don't know if we really need OpenMP or not. If we don't use
it, we don't have to worry about threads having to wait for each
other.

It's up to you of course. Lot's of fun to play with OpenMP anyway.

The suggestions you've given do seem like they should speed things up
quite a bit.

Bill.

2008/8/1 Peter Shrimpton <ps...@gm...>:
> Hello
>
> I have written some code to do the 2-SPRP test Lucas test and Pocklington
> code in parallel (See attached). I am fairly certain it is working but
> unfortunately it is pretty slow (and currently only dealing with numbers up
> to 10^10 because it is only 32bits). I think the following will speed it up:
>
> 1)Do multiple Pocklington tests at once
> 2)Redo the parallel sieving part: When a thread reaches the end of the
> current sieved chunk it starts sieving the next chunk. However all of the
> other threads have to wait until it has finished sieving before they can
> continue.
> 3)Store part of the Lucas chain. Currently every time it does a Lucas test
> it calculates the Lucas chain from the beginning which is clearly not
> necessary.
>
> 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] Parallel Stuff

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

Hello

I have written some code to do the 2-SPRP test Lucas test and Pocklington
code in parallel (See attached). I am fairly certain it is working but
unfortunately it is pretty slow (and currently only dealing with numbers up
to 10^10 because it is only 32bits). I think the following will speed it up:

1)Do multiple Pocklington tests at once
2)Redo the parallel sieving part: When a thread reaches the end of the
current sieved chunk it starts sieving the next chunk. However all of the
other threads have to wait until it has finished sieving before they can
continue.
3)Store part of the Lucas chain. Currently every time it does a Lucas test
it calculates the Lucas chain from the beginning which is clearly not
necessary.

Peter

Re: [flint-devel] Square-Free Algorithm

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

Hi Richard,

I merged your squarefree factoring algorithm in. I made lot's of
cosmetic changes.

I think I have it working OK, and I spent a bit of time figuring out
ways to speed it up. I added some more thoughts about this at the end
of todo.txt in the FLINT source tree.

There were two major things wrong with it the way you had implemented it:

1) When you had to add f(x)^p to the factor array you were adding f(x)
over and over, p times. If p is a 16 digit number, this would never
finish and would take more memory that we'd ever have. So I added an
exponents array to the factor structure, and now I just store the p in
the exponent for that factor. It's instant and takes no memory.

2) There was an error in the version of the algorithm you were using.
The thesis has an error. Basically you have to raise z(x) to the power
of i, where i is the number of iterations you have done, before adding
it to the factor array. It seems to give correct results now, though
I've only tried it on a handful of polynomials. You can fiddle around
with it if you like by changing the polynomials in the test code,
which is not really written properly yet. We need is_irreducible
before we can do that.

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,

I noticed here that you have a function for raising a polynomial to a
power modulo another poly. I just added such a function to zmod_poly
which should be more efficient. It uses the repeated squaring trick.
It is called zmod_poly_powmod. I also added a function for multiplying
two polynomials modulo another, called zmod_poly_mulmod.

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
>
>

Re: [flint-devel] Function Interface Updates

[<D.C...@wa...>]

Hi,

I've been thinking a little more about how our p-adics could work in
FLINT. Am I right in thinking that fractions aren't yet implemented in
FLINT? If so, it probably wouldn't be difficult to get around the problem
of inputting a fraction, but its probably not best to bother at this stage
as everything would need redoing at a later stage. Although having said
that, I do remember you saying that when dealing with polynomials over
Q[X] you store an element of Z[X] and an integer representing the
denominator so maybe we should proceed in a similar way to this? Either
way, I've included my thoughts for both cases.


Assuming for now that we just input integers;

1) Input an integer say a and a precision n.

2) Performing arithmetic; Input integers a and b.

        * Store a and b and the precision n.
        * To perform multiplication,  find a*b
        * Find the p-adic expansion of a*b.
        * Print to screen.


I need to look at ways of converting integer -> p-adic expansion quickly.
I assume this conversion would want its own .c file in FLINT?

Alternatively, if a user inputs a p-adic expansion directly then we would
want to 'sum it up' to give us an integer to store. Again, would this need
its own .c file?


If we were to allow fractions as inputs, then this is how I imagine it
working;

1) Storing a p-adic; First input a rational, say a/b and a precision n.
        *Check hcf(a,b) = 1. If not then cancel their common factors.
        *Find v_p(b) where v_p denotes the p-adic valuation.
                --> If v_p(b) = N > 0, then write a/b = p^(-N)*a/b', so
that a/b' lies in Z_p.
                           Store a, b' and N.
                --> Otherwise store the integers a and b.



2) Performing arithmetic; Input rationals a/b and c/d, where we can assume
a/b and
c/d lie in Z_p.
        * Store a, b, c, d and the precision n.
        * To perform division, i.e. (a/b)/(c/d), find a*d and b*c, and
cancel common factors in ad/bc.
        * Find the p-adic expansion of ad/bc;
                --> Reduce a*d mod p^n == X and  b*c mod p^n == Y first to
make the following
                      divisions smaller
                --> Find the p-adic expansion of X/Y to precision n. (I
see that FLINT can do congruences of the form (1/Y)mod
p^k).
        * Print to screen.


I need to look through FLINT to see how everything would link together.

I've been looking into ramified extensions a little over the last few days
and will be spending some more time on this over the next few days. It may
turn out that in order to deal with finite extensions of p-adic fields we
may want to deal with p-adics a little differently, which is another
reason to be looking into it further.


--------------------------------------------------

This isn't relevant now, but I'll include it anyway;

I've also been thinking about how to deal with the p-adic expansion of an
irrational and whether or not this will be useful. We know that the field
Q_p is strictly bigger than Q, and it is easy to see when a number sqrt(x)
for x a natural number lies in Q_p; I believe that a necessary and
sufficient condition on x such that sqrt(x) lies in Q_p for some p is that
m^2 == x mod p has a solution for some m, 0=<m<p (then m is the first term
in the p-adic expansion of m). It is
not clear to me then how to proceed from here other than by the 'by-hand'
method
which wouldn't be easy at all to implement. Obviously we can't approximate
x by a
rational, i.e. sqrt(2) ~= 141/100, etc in the p-adics. So unless theres an
easier
way to compute the p-adic expansion of an irrational I think its going to
be hard to
do and with little benefit. I notice that PARI won't give the p-adic
expansion of
say sqrt(2) in Q_7, despite the face that sqrt(2) lies in Q_7.


--------------------------------------------------

Thanks,

Daniel Ellam





> Hi,
>
> I've been considering the basic operations +,- ,* ,/  for the function
> interface and converting between p-adic and integer/rational
> representations.

Re: [flint-devel] host-57-44 reinstalled

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

Hi all,

the SUSE linux I used installed a 32 bit kernel on
host-57-44, which is a 64 bit Opteron machine.

Thus it needs to be reinstalled again.

-1 (again) for the SUSE developers for not telling me
this was a 32 bit kernel only.

Bill.


      

Re: [flint-devel] Making FLINT

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

Looks like you do. Try

export LD_LIBRARY_PATH=/usr/local/mpfr-2.3.1:$LD_LIBRARY_PATH

Bill.

2008/7/29 Peter Shrimpton <ps...@gm...>:
> Hey
>
> I followed the instructions in the last email but when I try t make library
> using gcc 4.3 it comes up with the error:
>
> /usr/local/gcc-4.3.1/libexec/gcc/i686-pc-linux-gnu/4.3.1/cc1: error while
> loading shared libraries: libmpfr.so.1: cannot open shared object file: No
> such file or directory
>
> Do I need to tell it where mpfr is? If so how?
>
> 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] Function Interface Updates

[<D.C...@wa...>]

Hi,

It sounds like you're busy fixing the computer so theres no rush to get
back to me on any of this, Ive just outlined what I've been thinking about
with the p-adics.


I've been considering the basic operations +,- ,* ,/  for the function
interface and converting between p-adic and integer/rational
representations.

When the user creates a p-adic, they are either going to want to input it
as a rational/integer, say x = 16 + O(3^5) or directly as x = 1 + 2*3 +
3^2 +O(3^5) . Working with elements of Z_p = {a/b | p doesnt divide p} for
now, it is possible that the user may want the p-adic expansion of say
16/5. So when we convert an integer to its p-adic expansion, one way is
obviously just to solve congruences mod p^k (I seem to remember you saying
that there was a quicker way to do this. Theres always repeated squaring
but I don't think it would offer any advantage here?). This method for
integer -> p-adic extends for dealing with elements of the form a/b when p
doesnt divide b since we are able to solve the congruences 1/b mod p^k for
all k. So unless theres a better way or I'm missing something, to find the
p-adic expansion of a/b we would need to find a modp^k, (1/b)modp^k and
then reduce their product a*(1/b)mod p^k for eack k as necessary. Can
FLINT solve congruences like 1/b mod p^k already (as you'd expect, my C++
compiler complains)?

Suppose we were working with precision n. If you inputted a p-adic of the
form a/b +O(p^n), I'm not sure yet what would be better; either 1: Storing
the two integers a and b straight away, or 2: Just store the integer which
approximates the p-adic expansion of the quotient a/b to precision n. For
example,

Performing additions in Case 1;

Input two rationals, a/b and c/d. Store 4 integers a,b,c,d. To find a/b +
c/d we find X = a*d + b*c and Y = b*d. Then X and Y will be large, so
reduce them mod p^n to make the following divisions smaller; find X*(1/Y)
mod p^k by finding X mod p^k and 1/Y mod p^k for each k and then find
their product mod p^k to get its p-adic expansion.


Performing additions in Case 2;

Input two rationals, a/b and c/d and a precision n. Find the expansion of
a/b and c/d by finding a*(1/b)mod p^k and c*(1/d)mod p^k (again it looks
like we would have to find (1/b)modp^k, etc) to approximate a/b as P and
c/d as Q. We store two integers P and Q. Then to perform the addition a/b
+ c/d we just need P + Q and then convert this to its p-adic expansion to
precision n.


Unless I'm missing something, I think case 1, storing 4 integers straight
away will be quicker. However, if you were going to be performing the same
operation many times with say a while loop, case 2 would become quicker
after so many iterations since you need fewer integer operations. Which
would you think would be quicker specifically in FLINT?


One problem with case 1 is that the functions would have to return two
values, numberator and denominator when functions will only return one. I
know it is possible to get around that problem using references or
pointers, would we have to use this here? If we did need to do this, what
would be better for FLINT, pointers or references?


Thanks,

Daniel Ellam

[flint-devel] Making FLINT

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

Hey

I followed the instructions in the last email but when I try t make library
using gcc 4.3 it comes up with the error:

/usr/local/gcc-4.3.1/libexec/gcc/i686-pc-linux-gnu/4.3.1/cc1: error while
loading shared libraries: libmpfr.so.1: cannot open shared object file: No
such file or directory

Do I need to tell it where mpfr is? If so how?

Peter

[flint-devel] host-57-44 reinstalled

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

Hi all,

as most of you know, when I installed gcc-4.3.1 on
host-57-44 I did "make install" into the default
location which ended up mixing two gcc's with
disasterour consequences. -1 for the gcc developers.

Unfortunately gcc "doesn't support make uninstall". -1
for the gcc developers. 

So I tried to uninstall it myself. This is impossible,
so the machine died.

I had to reinstall SUSE linux. However, SUSE linux in
its infinite wisdom saw fit to *completely format* the
home directories of everyone as default whilst making
it completely opaque that this is what it was doing.
-1 for SUSE developers.

It's all back up and running now, minus anything you
had in your home directories. Very sorry about that.
We all lost everything. This is just not supposed to
happen, as the home directories are purposely
installed in a different partition to / But SUSE has
other ideas if it finds unformated unallocated space
on the drive. 

Anyhow, I have now:

1) Fixed the CDROM drive on the machine (it was not
connected - physically).

2) Reinstalled Open SUSE 10.2. I put pdflatex on this
time too. /home and /usr/local are separate partitions
from / now.

3) Installed a new fast version of GMP that I and some
colleagues have been working on, called MPIR. You'll
find it in /usr/local/mpir-0.9.0 It is currently
functionally equivalent to GMP. Most GNU software that
requires GMP allows you to specify where GMP is with
./configure --with-gmp=/usr/local/mpir-0.9.0/

To link FLINT against MPIR you must do:

export FLINT_GMP_LIB_DIR=/usr/local/mpir-0.9.0/lib
export
FLINT_GMP_INCLUDE_DIR=/usr/local/mpir-0.9.0/include
source flint_env
make library

to link your program against MPIR and FLINT you can do
-L/usr/local/mpir-0.9.0/lib -L/path_to_flint -lgmp
-lflint

where obviously path_to_flint is wherever you
installed FLINT.

If you use gcc-4.3.1 you *must* use MPIR, since the
standard system GMP 4.2.1 doesn't support gcc-4.3.1.

4) Installed mpfr in /usr/local/mpfr-2.3.1

5) Installed gcc-4.3.1 in /usr/local/gcc-4.3.1 To use
it instead of the system gcc (4.1.2) you can either
export CC=/usr/local/gcc-4.3.1/bin (in the case of
FLINT that should be FLINT_CC) or even easier still,
just do export PATH=/usr/local/gcc-4.3.1/bin:$PATH

gcc-4.3.1 now works fine (it builds FLINT - what more
could you want!!)

6) Installed pari/gp in /usr/local. If you just type
gp at the prompt it will run from anywhere. I built it
against MPIR and so it will be faster than if it were
built against GMP.
 
I have not installed Magma yet, for lack of a pass
file. But it is coming.

John and Samir I have not put home directories back on
for you, but the root password remains unchanged. Or I
can do it for you if you ask. 

Bill.


      

Re: [flint-devel] Computer

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

I think most of the functions that I had written are in the email with the
subject "Additions to long_extras". There was some other code but it was
fairly simple code testing combinations of the SPRP and Lucas tests that I
can re wright quite easily.

Peter

2008/7/28 Bill Hart <goo...@go...>

> Guy's, I am about 90% sure we've had a computing disaster!
>
> I think in my attempt to restore this machine it has completely wiped
> the home directories.
>
> I have most of Richards's code that he sent to me via email and some
> version of Peter's Pocklington code, hopefully the most recent.
>
> What else have we lost? Did you guys have anything backed up elsewhere
> by any chance.
>
> Really sorry about this. I'm doing what I can to get the machine back
> online now. It will probably be some time tomorrow before it is ready
> to log in again and I'll have to set up home directories for you
> again.
>
> Bill.
>
> -------------------------------------------------------------------------
> 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] Computer

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

Guy's, I am about 90% sure we've had a computing disaster!

I think in my attempt to restore this machine it has completely wiped
the home directories.

I have most of Richards's code that he sent to me via email and some
version of Peter's Pocklington code, hopefully the most recent.

What else have we lost? Did you guys have anything backed up elsewhere
by any chance.

Really sorry about this. I'm doing what I can to get the machine back
online now. It will probably be some time tomorrow before it is ready
to log in again and I'll have to set up home directories for you
again.

Bill.