|
From: Bill H. <goo...@go...> - 2008-11-07 03:08:08
|
I implemented the heuristic gcd for polynomials whose bitsize is up to 63 bits. The result looks like this: http://sage.math.washington.edu/home/wbhart/flint-trunk/graphing/gcd9.png Interestingly although it makes the bottom of the graph nice and blue, the big red patch on the left is still there. I think perhaps this is just a highly efficient euclidean algorithm, or perhaps a Lehmer GCD algorithm. The red on the right is nor due to FLINT. That is due to a lack of integer half-gcd in GMP. We're still working on fixing that. Bill. |
|
From: William H. <ha...@ya...> - 2008-11-07 17:03:54
|
I decided to remove the test for correctness in the Heuristic GCD (which it actually can't do without) and see if that sped things up. It did: http://sage.math.washington.edu/home/wbhart/flint-trunk/graphing/gcd10.png The red on the bottom left turned blue. So it looks like it is the function fmpz_poly_divides which checks if the purported GCD actually divides the original polynomials is taking most of the time (like 80%). One way of doing a quick check is to pack the coefficients of the polynomials into large integers and divide. If the integer divides exactly one gets the quotient, which one can unpack into a polynomial and then if that quotient times the divisor is the dividend, one can prove that the polynomials divide exactly. Clearly the brown at the bottom of the graph on the left is because Magma is packing (both in the GCD and the divides function) into 32 bits instead of 64. So I think it is clear now what has to be implemented. I'm quite surprised that the fmpz_poly_divides_modular function was not sufficiently fast. It seems one can do a lot better if exact division is expected. Bill. --- On Fri, 11/7/08, Bill Hart <goo...@go...> wrote: > From: Bill Hart <goo...@go...> > Subject: [flint-devel] Heuristic gcd > To: "Development list for FLINT" <fli...@li...> > Date: Friday, November 7, 2008, 2:07 PM > I implemented the heuristic gcd for polynomials whose > bitsize is up to > 63 bits. The result looks like this: > > http://sage.math.washington.edu/home/wbhart/flint-trunk/graphing/gcd9.png > > Interestingly although it makes the bottom of the graph > nice and blue, > the big red patch on the left is still there. > > I think perhaps this is just a highly efficient euclidean > algorithm, > or perhaps a Lehmer GCD algorithm. > > The red on the right is nor due to FLINT. That is due to a > lack of > integer half-gcd in GMP. We're still working on fixing > that. > > 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 |
|
From: Bill H. <goo...@go...> - 2008-11-08 15:58:48
|
I finally realised that if the number of bits one is packing into is bigger than the Mignotte bound, then one does not need to do a divides test at all. This is all that Magma is doing. Here is why: To do a divides, if one is under the Mignotte bound, one would pack A and B and pack the purported GCD and see if the packed A and B divide the packed GCD. This would be done by simply doing a multiprecision division and seeing if there is a remainder. If not, they divide precisely and the quotient is calculated exactly. But in computing the heuristic gcd, we already have the packed versions of A and B, as the gcd of these divides both of them by definition. But the packed version of G is precisely the gcd of the packed version of A and B. I've now inserted code to compute the Mignotte bound and omit the divides test in this case. The function passes the test code and I am now doing a profile run against Magma again: http://sage.math.washington.edu/home/wbhart/flint-trunk/graphing/gcd10.png The light blue layer at the bottom left is due to the fact that FLINT packs into 64 bits not 30 or 32 as Magma surely does. This means that we more easily exceed the Mignotte bound so that the divides test is not required. The browny red region at the bottom left is surely due to the fact that Magma is doing a heuristic gcd packing into 32 bits (or perhaps 30) instead of 64. Bill. 2008/11/7 William Hart <ha...@ya...>: > I decided to remove the test for correctness in the Heuristic GCD (which it actually can't do without) and see if that sped things up. It did: > > http://sage.math.washington.edu/home/wbhart/flint-trunk/graphing/gcd10.png > > The red on the bottom left turned blue. So it looks like it is the function fmpz_poly_divides which checks if the purported GCD actually divides the original polynomials is taking most of the time (like 80%). > > One way of doing a quick check is to pack the coefficients of the polynomials into large integers and divide. If the integer divides exactly one gets the quotient, which one can unpack into a polynomial and then if that quotient times the divisor is the dividend, one can prove that the polynomials divide exactly. > > Clearly the brown at the bottom of the graph on the left is because Magma is packing (both in the GCD and the divides function) into 32 bits instead of 64. > > So I think it is clear now what has to be implemented. > > I'm quite surprised that the fmpz_poly_divides_modular function was not sufficiently fast. It seems one can do a lot better if exact division is expected. > > Bill. > > --- On Fri, 11/7/08, Bill Hart <goo...@go...> wrote: > >> From: Bill Hart <goo...@go...> >> Subject: [flint-devel] Heuristic gcd >> To: "Development list for FLINT" <fli...@li...> >> Date: Friday, November 7, 2008, 2:07 PM >> I implemented the heuristic gcd for polynomials whose >> bitsize is up to >> 63 bits. The result looks like this: >> >> http://sage.math.washington.edu/home/wbhart/flint-trunk/graphing/gcd9.png >> >> Interestingly although it makes the bottom of the graph >> nice and blue, >> the big red patch on the left is still there. >> >> I think perhaps this is just a highly efficient euclidean >> algorithm, >> or perhaps a Lehmer GCD algorithm. >> >> The red on the right is nor due to FLINT. That is due to a >> lack of >> integer half-gcd in GMP. We're still working on fixing >> that. >> >> 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 > > > > > ------------------------------------------------------------------------- > 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 > |
|
From: Bill H. <goo...@go...> - 2008-11-09 00:19:46
|
I coded up the packing code to pack into 32 bits instead of 64. There were a couple of problems with this. When the Mignotte bound is bigger than 32 bits then a divides gets carried out, which for small gcd's dominates the time. Thus sometimes it is better to pack into 64 bits instead. Once that is sorted out one soon realises that just computing the Mignotte bound takes too long. So I implemented an approximation which is sometimes slightly too large but isn't too bad. Once these issues are fixed, here is the result: http://sage.math.washington.edu/home/wbhart/flint-trunk/graphing/gcd11.png Now the brown in the bottom left must be where packing is done to the bit. Here the Minotte bound is usually less than the number of bits that one is packing to. Bill. 2008/11/8 Bill Hart <goo...@go...>: > I finally realised that if the number of bits one is packing into is > bigger than the Mignotte bound, then one does not need to do a divides > test at all. This is all that Magma is doing. > > Here is why: > > To do a divides, if one is under the Mignotte bound, one would pack A > and B and pack the purported GCD and see if the packed A and B divide > the packed GCD. This would be done by simply doing a multiprecision > division and seeing if there is a remainder. If not, they divide > precisely and the quotient is calculated exactly. > > But in computing the heuristic gcd, we already have the packed > versions of A and B, as the gcd of these divides both of them by > definition. But the packed version of G is precisely the gcd of the > packed version of A and B. > > I've now inserted code to compute the Mignotte bound and omit the > divides test in this case. The function passes the test code and I am > now doing a profile run against Magma again: > > http://sage.math.washington.edu/home/wbhart/flint-trunk/graphing/gcd10.png > > The light blue layer at the bottom left is due to the fact that FLINT > packs into 64 bits not 30 or 32 as Magma surely does. This means that > we more easily exceed the Mignotte bound so that the divides test is > not required. > > The browny red region at the bottom left is surely due to the fact > that Magma is doing a heuristic gcd packing into 32 bits (or perhaps > 30) instead of 64. > > Bill. > > 2008/11/7 William Hart <ha...@ya...>: >> I decided to remove the test for correctness in the Heuristic GCD (which it actually can't do without) and see if that sped things up. It did: >> >> http://sage.math.washington.edu/home/wbhart/flint-trunk/graphing/gcd10.png >> >> The red on the bottom left turned blue. So it looks like it is the function fmpz_poly_divides which checks if the purported GCD actually divides the original polynomials is taking most of the time (like 80%). >> >> One way of doing a quick check is to pack the coefficients of the polynomials into large integers and divide. If the integer divides exactly one gets the quotient, which one can unpack into a polynomial and then if that quotient times the divisor is the dividend, one can prove that the polynomials divide exactly. >> >> Clearly the brown at the bottom of the graph on the left is because Magma is packing (both in the GCD and the divides function) into 32 bits instead of 64. >> >> So I think it is clear now what has to be implemented. >> >> I'm quite surprised that the fmpz_poly_divides_modular function was not sufficiently fast. It seems one can do a lot better if exact division is expected. >> >> Bill. >> >> --- On Fri, 11/7/08, Bill Hart <goo...@go...> wrote: >> >>> From: Bill Hart <goo...@go...> >>> Subject: [flint-devel] Heuristic gcd >>> To: "Development list for FLINT" <fli...@li...> >>> Date: Friday, November 7, 2008, 2:07 PM >>> I implemented the heuristic gcd for polynomials whose >>> bitsize is up to >>> 63 bits. The result looks like this: >>> >>> http://sage.math.washington.edu/home/wbhart/flint-trunk/graphing/gcd9.png >>> >>> Interestingly although it makes the bottom of the graph >>> nice and blue, >>> the big red patch on the left is still there. >>> >>> I think perhaps this is just a highly efficient euclidean >>> algorithm, >>> or perhaps a Lehmer GCD algorithm. >>> >>> The red on the right is nor due to FLINT. That is due to a >>> lack of >>> integer half-gcd in GMP. We're still working on fixing >>> that. >>> >>> 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 >> >> >> >> >> ------------------------------------------------------------------------- >> 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 >> > |
|
From: Bill H. <goo...@go...> - 2008-11-09 20:54:34
|
Hmm, I think I have gotten ahead of myself with my argument about the Mignotte bound. If f(x) = a(x)*b(x) and g(x) = a(x)*d(x) with b(x) and d(x) coprime. when I substitute x = 2^64 then certainly a(2^64) will be a common factor of f(2^64) and g(2^64). But it is also possible that even though b(x) and d(x) are coprime that b(2^64) and d(2^64) are not. But this means that we can't simply rely on the Mignotte bound for the size of the gcd of f(x) and g(x). Certainly if 2^64 is bigger than the Mignotte bound times the gcd of b(2^64) and d(2^64) then we are ok. The polynomial gcd that we are after will just end up being multiplied by some integer constant which we can remove by computing the content. But there is no way that I can see to bound gcd of b(2^64) and d(2^64). Thus it seems essential to me that one actually check that the purported poly gcd actually divides f(x) and g(x). But the time it takes to do this "divides" test is too great. I don't have a solution to this problem at present. Bill. 2008/11/9 Bill Hart <goo...@go...>: > I coded up the packing code to pack into 32 bits instead of 64. There > were a couple of problems with this. When the Mignotte bound is bigger > than 32 bits then a divides gets carried out, which for small gcd's > dominates the time. Thus sometimes it is better to pack into 64 bits > instead. > > Once that is sorted out one soon realises that just computing the > Mignotte bound takes too long. So I implemented an approximation which > is sometimes slightly too large but isn't too bad. > > Once these issues are fixed, here is the result: > > http://sage.math.washington.edu/home/wbhart/flint-trunk/graphing/gcd11.png > > Now the brown in the bottom left must be where packing is done to the > bit. Here the Minotte bound is usually less than the number of bits > that one is packing to. > > Bill. > > 2008/11/8 Bill Hart <goo...@go...>: >> I finally realised that if the number of bits one is packing into is >> bigger than the Mignotte bound, then one does not need to do a divides >> test at all. This is all that Magma is doing. >> >> Here is why: >> >> To do a divides, if one is under the Mignotte bound, one would pack A >> and B and pack the purported GCD and see if the packed A and B divide >> the packed GCD. This would be done by simply doing a multiprecision >> division and seeing if there is a remainder. If not, they divide >> precisely and the quotient is calculated exactly. >> >> But in computing the heuristic gcd, we already have the packed >> versions of A and B, as the gcd of these divides both of them by >> definition. But the packed version of G is precisely the gcd of the >> packed version of A and B. >> >> I've now inserted code to compute the Mignotte bound and omit the >> divides test in this case. The function passes the test code and I am >> now doing a profile run against Magma again: >> >> http://sage.math.washington.edu/home/wbhart/flint-trunk/graphing/gcd10.png >> >> The light blue layer at the bottom left is due to the fact that FLINT >> packs into 64 bits not 30 or 32 as Magma surely does. This means that >> we more easily exceed the Mignotte bound so that the divides test is >> not required. >> >> The browny red region at the bottom left is surely due to the fact >> that Magma is doing a heuristic gcd packing into 32 bits (or perhaps >> 30) instead of 64. >> >> Bill. >> >> 2008/11/7 William Hart <ha...@ya...>: >>> I decided to remove the test for correctness in the Heuristic GCD (which it actually can't do without) and see if that sped things up. It did: >>> >>> http://sage.math.washington.edu/home/wbhart/flint-trunk/graphing/gcd10.png >>> >>> The red on the bottom left turned blue. So it looks like it is the function fmpz_poly_divides which checks if the purported GCD actually divides the original polynomials is taking most of the time (like 80%). >>> >>> One way of doing a quick check is to pack the coefficients of the polynomials into large integers and divide. If the integer divides exactly one gets the quotient, which one can unpack into a polynomial and then if that quotient times the divisor is the dividend, one can prove that the polynomials divide exactly. >>> >>> Clearly the brown at the bottom of the graph on the left is because Magma is packing (both in the GCD and the divides function) into 32 bits instead of 64. >>> >>> So I think it is clear now what has to be implemented. >>> >>> I'm quite surprised that the fmpz_poly_divides_modular function was not sufficiently fast. It seems one can do a lot better if exact division is expected. >>> >>> Bill. >>> >>> --- On Fri, 11/7/08, Bill Hart <goo...@go...> wrote: >>> >>>> From: Bill Hart <goo...@go...> >>>> Subject: [flint-devel] Heuristic gcd >>>> To: "Development list for FLINT" <fli...@li...> >>>> Date: Friday, November 7, 2008, 2:07 PM >>>> I implemented the heuristic gcd for polynomials whose >>>> bitsize is up to >>>> 63 bits. The result looks like this: >>>> >>>> http://sage.math.washington.edu/home/wbhart/flint-trunk/graphing/gcd9.png >>>> >>>> Interestingly although it makes the bottom of the graph >>>> nice and blue, >>>> the big red patch on the left is still there. >>>> >>>> I think perhaps this is just a highly efficient euclidean >>>> algorithm, >>>> or perhaps a Lehmer GCD algorithm. >>>> >>>> The red on the right is nor due to FLINT. That is due to a >>>> lack of >>>> integer half-gcd in GMP. We're still working on fixing >>>> that. >>>> >>>> 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 >>> >>> >>> >>> >>> ------------------------------------------------------------------------- >>> 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 >>> >> > |
|
From: Bill H. <goo...@go...> - 2008-11-10 22:31:34
|
It turns out that Schoenhage already computed a bound that works for the heuristic gcd unconditionally. If A is the max norm of the polynomials f and g and n is the maximum degree then one can evaluate at an integer u with u > 4*(n+1)^n*A^(2n) and not have to run divides checks. If Magma were doing this in the red region of http://sage.math.washington.edu/home/wbhart/flint-trunk/graphing/gcd.png then we would have n = 32 and A = 2^16 say. The u would be 4*(33^32)*2^(16*64) > (2^64)^18, i.e. u would have to be bigger than 18 limbs!! So this proven bound is totally useless. If Magma is not doing divisions in this region, I would question whether the results would be correct. No doubt they do something else clever. Bill. 2008/11/9 Bill Hart <goo...@go...>: > Hmm, I think I have gotten ahead of myself with my argument about the > Mignotte bound. > > If f(x) = a(x)*b(x) and g(x) = a(x)*d(x) with b(x) and d(x) coprime. > > when I substitute x = 2^64 then certainly a(2^64) will be a common > factor of f(2^64) and g(2^64). > > But it is also possible that even though b(x) and d(x) are coprime > that b(2^64) and d(2^64) are not. > > But this means that we can't simply rely on the Mignotte bound for the > size of the gcd of f(x) and g(x). Certainly if 2^64 is bigger than the > Mignotte bound times the gcd of b(2^64) and d(2^64) then we are ok. > The polynomial gcd that we are after will just end up being multiplied > by some integer constant which we can remove by computing the content. > > But there is no way that I can see to bound gcd of b(2^64) and d(2^64). > > Thus it seems essential to me that one actually check that the > purported poly gcd actually divides f(x) and g(x). > > But the time it takes to do this "divides" test is too great. I don't > have a solution to this problem at present. > > Bill. > > 2008/11/9 Bill Hart <goo...@go...>: >> I coded up the packing code to pack into 32 bits instead of 64. There >> were a couple of problems with this. When the Mignotte bound is bigger >> than 32 bits then a divides gets carried out, which for small gcd's >> dominates the time. Thus sometimes it is better to pack into 64 bits >> instead. >> >> Once that is sorted out one soon realises that just computing the >> Mignotte bound takes too long. So I implemented an approximation which >> is sometimes slightly too large but isn't too bad. >> >> Once these issues are fixed, here is the result: >> >> http://sage.math.washington.edu/home/wbhart/flint-trunk/graphing/gcd11.png >> >> Now the brown in the bottom left must be where packing is done to the >> bit. Here the Minotte bound is usually less than the number of bits >> that one is packing to. >> >> Bill. >> >> 2008/11/8 Bill Hart <goo...@go...>: >>> I finally realised that if the number of bits one is packing into is >>> bigger than the Mignotte bound, then one does not need to do a divides >>> test at all. This is all that Magma is doing. >>> >>> Here is why: >>> >>> To do a divides, if one is under the Mignotte bound, one would pack A >>> and B and pack the purported GCD and see if the packed A and B divide >>> the packed GCD. This would be done by simply doing a multiprecision >>> division and seeing if there is a remainder. If not, they divide >>> precisely and the quotient is calculated exactly. >>> >>> But in computing the heuristic gcd, we already have the packed >>> versions of A and B, as the gcd of these divides both of them by >>> definition. But the packed version of G is precisely the gcd of the >>> packed version of A and B. >>> >>> I've now inserted code to compute the Mignotte bound and omit the >>> divides test in this case. The function passes the test code and I am >>> now doing a profile run against Magma again: >>> >>> http://sage.math.washington.edu/home/wbhart/flint-trunk/graphing/gcd10.png >>> >>> The light blue layer at the bottom left is due to the fact that FLINT >>> packs into 64 bits not 30 or 32 as Magma surely does. This means that >>> we more easily exceed the Mignotte bound so that the divides test is >>> not required. >>> >>> The browny red region at the bottom left is surely due to the fact >>> that Magma is doing a heuristic gcd packing into 32 bits (or perhaps >>> 30) instead of 64. >>> >>> Bill. >>> >>> 2008/11/7 William Hart <ha...@ya...>: >>>> I decided to remove the test for correctness in the Heuristic GCD (which it actually can't do without) and see if that sped things up. It did: >>>> >>>> http://sage.math.washington.edu/home/wbhart/flint-trunk/graphing/gcd10.png >>>> >>>> The red on the bottom left turned blue. So it looks like it is the function fmpz_poly_divides which checks if the purported GCD actually divides the original polynomials is taking most of the time (like 80%). >>>> >>>> One way of doing a quick check is to pack the coefficients of the polynomials into large integers and divide. If the integer divides exactly one gets the quotient, which one can unpack into a polynomial and then if that quotient times the divisor is the dividend, one can prove that the polynomials divide exactly. >>>> >>>> Clearly the brown at the bottom of the graph on the left is because Magma is packing (both in the GCD and the divides function) into 32 bits instead of 64. >>>> >>>> So I think it is clear now what has to be implemented. >>>> >>>> I'm quite surprised that the fmpz_poly_divides_modular function was not sufficiently fast. It seems one can do a lot better if exact division is expected. >>>> >>>> Bill. >>>> >>>> --- On Fri, 11/7/08, Bill Hart <goo...@go...> wrote: >>>> >>>>> From: Bill Hart <goo...@go...> >>>>> Subject: [flint-devel] Heuristic gcd >>>>> To: "Development list for FLINT" <fli...@li...> >>>>> Date: Friday, November 7, 2008, 2:07 PM >>>>> I implemented the heuristic gcd for polynomials whose >>>>> bitsize is up to >>>>> 63 bits. The result looks like this: >>>>> >>>>> http://sage.math.washington.edu/home/wbhart/flint-trunk/graphing/gcd9.png >>>>> >>>>> Interestingly although it makes the bottom of the graph >>>>> nice and blue, >>>>> the big red patch on the left is still there. >>>>> >>>>> I think perhaps this is just a highly efficient euclidean >>>>> algorithm, >>>>> or perhaps a Lehmer GCD algorithm. >>>>> >>>>> The red on the right is nor due to FLINT. That is due to a >>>>> lack of >>>>> integer half-gcd in GMP. We're still working on fixing >>>>> that. >>>>> >>>>> 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 >>>> >>>> >>>> >>>> >>>> ------------------------------------------------------------------------- >>>> 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 >>>> >>> >> > |
|
From: Bill H. <goo...@go...> - 2008-11-11 00:08:49
|
After heating up one of the cores on sage.math for a while I have come to the conclusion that Magma is not cheating. It gets the right answers for these small GCD problems. But I am mystified as to how? I recall when timing the division code in FLINT that in that region Magma was much more efficient. Something seems to reek of inefficiency when two tests for divisibility dominate the time to do a GCD. Perhaps I just haven't thought of an efficient enough way to test divisibility. I would have thought the multi-modular technique should be fast enough, but it doesn't seem to be. Bill. 2008/11/10 Bill Hart <goo...@go...>: > It turns out that Schoenhage already computed a bound that works for > the heuristic gcd unconditionally. > > If A is the max norm of the polynomials f and g and n is the maximum > degree then one can evaluate at an integer u with > > u > 4*(n+1)^n*A^(2n) > > and not have to run divides checks. > > If Magma were doing this in the red region of > http://sage.math.washington.edu/home/wbhart/flint-trunk/graphing/gcd.png > then we would have n = 32 and A = 2^16 say. > > The u would be 4*(33^32)*2^(16*64) > (2^64)^18, i.e. u would have to > be bigger than 18 limbs!! So this proven bound is totally useless. > > If Magma is not doing divisions in this region, I would question > whether the results would be correct. No doubt they do something else > clever. > > Bill. > > 2008/11/9 Bill Hart <goo...@go...>: >> Hmm, I think I have gotten ahead of myself with my argument about the >> Mignotte bound. >> >> If f(x) = a(x)*b(x) and g(x) = a(x)*d(x) with b(x) and d(x) coprime. >> >> when I substitute x = 2^64 then certainly a(2^64) will be a common >> factor of f(2^64) and g(2^64). >> >> But it is also possible that even though b(x) and d(x) are coprime >> that b(2^64) and d(2^64) are not. >> >> But this means that we can't simply rely on the Mignotte bound for the >> size of the gcd of f(x) and g(x). Certainly if 2^64 is bigger than the >> Mignotte bound times the gcd of b(2^64) and d(2^64) then we are ok. >> The polynomial gcd that we are after will just end up being multiplied >> by some integer constant which we can remove by computing the content. >> >> But there is no way that I can see to bound gcd of b(2^64) and d(2^64). >> >> Thus it seems essential to me that one actually check that the >> purported poly gcd actually divides f(x) and g(x). >> >> But the time it takes to do this "divides" test is too great. I don't >> have a solution to this problem at present. >> >> Bill. >> >> 2008/11/9 Bill Hart <goo...@go...>: >>> I coded up the packing code to pack into 32 bits instead of 64. There >>> were a couple of problems with this. When the Mignotte bound is bigger >>> than 32 bits then a divides gets carried out, which for small gcd's >>> dominates the time. Thus sometimes it is better to pack into 64 bits >>> instead. >>> >>> Once that is sorted out one soon realises that just computing the >>> Mignotte bound takes too long. So I implemented an approximation which >>> is sometimes slightly too large but isn't too bad. >>> >>> Once these issues are fixed, here is the result: >>> >>> http://sage.math.washington.edu/home/wbhart/flint-trunk/graphing/gcd11.png >>> >>> Now the brown in the bottom left must be where packing is done to the >>> bit. Here the Minotte bound is usually less than the number of bits >>> that one is packing to. >>> >>> Bill. >>> >>> 2008/11/8 Bill Hart <goo...@go...>: >>>> I finally realised that if the number of bits one is packing into is >>>> bigger than the Mignotte bound, then one does not need to do a divides >>>> test at all. This is all that Magma is doing. >>>> >>>> Here is why: >>>> >>>> To do a divides, if one is under the Mignotte bound, one would pack A >>>> and B and pack the purported GCD and see if the packed A and B divide >>>> the packed GCD. This would be done by simply doing a multiprecision >>>> division and seeing if there is a remainder. If not, they divide >>>> precisely and the quotient is calculated exactly. >>>> >>>> But in computing the heuristic gcd, we already have the packed >>>> versions of A and B, as the gcd of these divides both of them by >>>> definition. But the packed version of G is precisely the gcd of the >>>> packed version of A and B. >>>> >>>> I've now inserted code to compute the Mignotte bound and omit the >>>> divides test in this case. The function passes the test code and I am >>>> now doing a profile run against Magma again: >>>> >>>> http://sage.math.washington.edu/home/wbhart/flint-trunk/graphing/gcd10.png >>>> >>>> The light blue layer at the bottom left is due to the fact that FLINT >>>> packs into 64 bits not 30 or 32 as Magma surely does. This means that >>>> we more easily exceed the Mignotte bound so that the divides test is >>>> not required. >>>> >>>> The browny red region at the bottom left is surely due to the fact >>>> that Magma is doing a heuristic gcd packing into 32 bits (or perhaps >>>> 30) instead of 64. >>>> >>>> Bill. >>>> >>>> 2008/11/7 William Hart <ha...@ya...>: >>>>> I decided to remove the test for correctness in the Heuristic GCD (which it actually can't do without) and see if that sped things up. It did: >>>>> >>>>> http://sage.math.washington.edu/home/wbhart/flint-trunk/graphing/gcd10.png >>>>> >>>>> The red on the bottom left turned blue. So it looks like it is the function fmpz_poly_divides which checks if the purported GCD actually divides the original polynomials is taking most of the time (like 80%). >>>>> >>>>> One way of doing a quick check is to pack the coefficients of the polynomials into large integers and divide. If the integer divides exactly one gets the quotient, which one can unpack into a polynomial and then if that quotient times the divisor is the dividend, one can prove that the polynomials divide exactly. >>>>> >>>>> Clearly the brown at the bottom of the graph on the left is because Magma is packing (both in the GCD and the divides function) into 32 bits instead of 64. >>>>> >>>>> So I think it is clear now what has to be implemented. >>>>> >>>>> I'm quite surprised that the fmpz_poly_divides_modular function was not sufficiently fast. It seems one can do a lot better if exact division is expected. >>>>> >>>>> Bill. >>>>> >>>>> --- On Fri, 11/7/08, Bill Hart <goo...@go...> wrote: >>>>> >>>>>> From: Bill Hart <goo...@go...> >>>>>> Subject: [flint-devel] Heuristic gcd >>>>>> To: "Development list for FLINT" <fli...@li...> >>>>>> Date: Friday, November 7, 2008, 2:07 PM >>>>>> I implemented the heuristic gcd for polynomials whose >>>>>> bitsize is up to >>>>>> 63 bits. The result looks like this: >>>>>> >>>>>> http://sage.math.washington.edu/home/wbhart/flint-trunk/graphing/gcd9.png >>>>>> >>>>>> Interestingly although it makes the bottom of the graph >>>>>> nice and blue, >>>>>> the big red patch on the left is still there. >>>>>> >>>>>> I think perhaps this is just a highly efficient euclidean >>>>>> algorithm, >>>>>> or perhaps a Lehmer GCD algorithm. >>>>>> >>>>>> The red on the right is nor due to FLINT. That is due to a >>>>>> lack of >>>>>> integer half-gcd in GMP. We're still working on fixing >>>>>> that. >>>>>> >>>>>> 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 >>>>> >>>>> >>>>> >>>>> >>>>> ------------------------------------------------------------------------- >>>>> 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 >>>>> >>>> >>> >> > |
|
From: Bill H. <goo...@go...> - 2008-11-13 14:58:20
|
For those following the polynomial GCD thread, I've moved it across to: http://groups.google.co.uk/group/sage-nt/topics?hl=en&start= in particular to this thread: http://groups.google.co.uk/group/sage-nt/browse_thread/thread/2da0e434f09d1fe1?hl=en in case anyone there has any ideas. For now I've completely run out of ideas. Magma has got me completely stumped. Here is the latest graph: http://sage.math.washington.edu/home/wbhart/flint-trunk/graphing/gcd13.png A description of the new improvements is on the sage-nt site linked above. Whilst the right hand side now looks pretty good, there is one problem. If the polynomials are coprime the Heuristic GCD is not a good algorithm to use for long polynomials. The reason is that the modular algorithm will only need to use O(1) primes to determine that the polys are coprime, whereas the Heuristic algorithm will go ahead and pack right up to the bit size of the input polynomials, which is much more expensive. So I will have to determine some crossover points experimentally. The brown we had on the right of the graph will come back and I will only be able to get rid of it by packing more than one coefficient into a single limb in the zmod_poly's, as that problem is a caching one. Bill. 2008/11/11 Bill Hart <goo...@go...>: > After heating up one of the cores on sage.math for a while I have come > to the conclusion that Magma is not cheating. It gets the right > answers for these small GCD problems. > > But I am mystified as to how? I recall when timing the division code > in FLINT that in that region Magma was much more efficient. > > Something seems to reek of inefficiency when two tests for > divisibility dominate the time to do a GCD. Perhaps I just haven't > thought of an efficient enough way to test divisibility. I would have > thought the multi-modular technique should be fast enough, but it > doesn't seem to be. > > Bill. > > 2008/11/10 Bill Hart <goo...@go...>: >> It turns out that Schoenhage already computed a bound that works for >> the heuristic gcd unconditionally. >> >> If A is the max norm of the polynomials f and g and n is the maximum >> degree then one can evaluate at an integer u with >> >> u > 4*(n+1)^n*A^(2n) >> >> and not have to run divides checks. >> >> If Magma were doing this in the red region of >> http://sage.math.washington.edu/home/wbhart/flint-trunk/graphing/gcd.png >> then we would have n = 32 and A = 2^16 say. >> >> The u would be 4*(33^32)*2^(16*64) > (2^64)^18, i.e. u would have to >> be bigger than 18 limbs!! So this proven bound is totally useless. >> >> If Magma is not doing divisions in this region, I would question >> whether the results would be correct. No doubt they do something else >> clever. >> >> Bill. >> >> 2008/11/9 Bill Hart <goo...@go...>: >>> Hmm, I think I have gotten ahead of myself with my argument about the >>> Mignotte bound. >>> >>> If f(x) = a(x)*b(x) and g(x) = a(x)*d(x) with b(x) and d(x) coprime. >>> >>> when I substitute x = 2^64 then certainly a(2^64) will be a common >>> factor of f(2^64) and g(2^64). >>> >>> But it is also possible that even though b(x) and d(x) are coprime >>> that b(2^64) and d(2^64) are not. >>> >>> But this means that we can't simply rely on the Mignotte bound for the >>> size of the gcd of f(x) and g(x). Certainly if 2^64 is bigger than the >>> Mignotte bound times the gcd of b(2^64) and d(2^64) then we are ok. >>> The polynomial gcd that we are after will just end up being multiplied >>> by some integer constant which we can remove by computing the content. >>> >>> But there is no way that I can see to bound gcd of b(2^64) and d(2^64). >>> >>> Thus it seems essential to me that one actually check that the >>> purported poly gcd actually divides f(x) and g(x). >>> >>> But the time it takes to do this "divides" test is too great. I don't >>> have a solution to this problem at present. >>> >>> Bill. >>> >>> 2008/11/9 Bill Hart <goo...@go...>: >>>> I coded up the packing code to pack into 32 bits instead of 64. There >>>> were a couple of problems with this. When the Mignotte bound is bigger >>>> than 32 bits then a divides gets carried out, which for small gcd's >>>> dominates the time. Thus sometimes it is better to pack into 64 bits >>>> instead. >>>> >>>> Once that is sorted out one soon realises that just computing the >>>> Mignotte bound takes too long. So I implemented an approximation which >>>> is sometimes slightly too large but isn't too bad. >>>> >>>> Once these issues are fixed, here is the result: >>>> >>>> http://sage.math.washington.edu/home/wbhart/flint-trunk/graphing/gcd11.png >>>> >>>> Now the brown in the bottom left must be where packing is done to the >>>> bit. Here the Minotte bound is usually less than the number of bits >>>> that one is packing to. >>>> >>>> Bill. >>>> >>>> 2008/11/8 Bill Hart <goo...@go...>: >>>>> I finally realised that if the number of bits one is packing into is >>>>> bigger than the Mignotte bound, then one does not need to do a divides >>>>> test at all. This is all that Magma is doing. >>>>> >>>>> Here is why: >>>>> >>>>> To do a divides, if one is under the Mignotte bound, one would pack A >>>>> and B and pack the purported GCD and see if the packed A and B divide >>>>> the packed GCD. This would be done by simply doing a multiprecision >>>>> division and seeing if there is a remainder. If not, they divide >>>>> precisely and the quotient is calculated exactly. >>>>> >>>>> But in computing the heuristic gcd, we already have the packed >>>>> versions of A and B, as the gcd of these divides both of them by >>>>> definition. But the packed version of G is precisely the gcd of the >>>>> packed version of A and B. >>>>> >>>>> I've now inserted code to compute the Mignotte bound and omit the >>>>> divides test in this case. The function passes the test code and I am >>>>> now doing a profile run against Magma again: >>>>> >>>>> http://sage.math.washington.edu/home/wbhart/flint-trunk/graphing/gcd10.png >>>>> >>>>> The light blue layer at the bottom left is due to the fact that FLINT >>>>> packs into 64 bits not 30 or 32 as Magma surely does. This means that >>>>> we more easily exceed the Mignotte bound so that the divides test is >>>>> not required. >>>>> >>>>> The browny red region at the bottom left is surely due to the fact >>>>> that Magma is doing a heuristic gcd packing into 32 bits (or perhaps >>>>> 30) instead of 64. >>>>> >>>>> Bill. >>>>> >>>>> 2008/11/7 William Hart <ha...@ya...>: >>>>>> I decided to remove the test for correctness in the Heuristic GCD (which it actually can't do without) and see if that sped things up. It did: >>>>>> >>>>>> http://sage.math.washington.edu/home/wbhart/flint-trunk/graphing/gcd10.png >>>>>> >>>>>> The red on the bottom left turned blue. So it looks like it is the function fmpz_poly_divides which checks if the purported GCD actually divides the original polynomials is taking most of the time (like 80%). >>>>>> >>>>>> One way of doing a quick check is to pack the coefficients of the polynomials into large integers and divide. If the integer divides exactly one gets the quotient, which one can unpack into a polynomial and then if that quotient times the divisor is the dividend, one can prove that the polynomials divide exactly. >>>>>> >>>>>> Clearly the brown at the bottom of the graph on the left is because Magma is packing (both in the GCD and the divides function) into 32 bits instead of 64. >>>>>> >>>>>> So I think it is clear now what has to be implemented. >>>>>> >>>>>> I'm quite surprised that the fmpz_poly_divides_modular function was not sufficiently fast. It seems one can do a lot better if exact division is expected. >>>>>> >>>>>> Bill. >>>>>> >>>>>> --- On Fri, 11/7/08, Bill Hart <goo...@go...> wrote: >>>>>> >>>>>>> From: Bill Hart <goo...@go...> >>>>>>> Subject: [flint-devel] Heuristic gcd >>>>>>> To: "Development list for FLINT" <fli...@li...> >>>>>>> Date: Friday, November 7, 2008, 2:07 PM >>>>>>> I implemented the heuristic gcd for polynomials whose >>>>>>> bitsize is up to >>>>>>> 63 bits. The result looks like this: >>>>>>> >>>>>>> http://sage.math.washington.edu/home/wbhart/flint-trunk/graphing/gcd9.png >>>>>>> >>>>>>> Interestingly although it makes the bottom of the graph >>>>>>> nice and blue, >>>>>>> the big red patch on the left is still there. >>>>>>> >>>>>>> I think perhaps this is just a highly efficient euclidean >>>>>>> algorithm, >>>>>>> or perhaps a Lehmer GCD algorithm. >>>>>>> >>>>>>> The red on the right is nor due to FLINT. That is due to a >>>>>>> lack of >>>>>>> integer half-gcd in GMP. We're still working on fixing >>>>>>> that. >>>>>>> >>>>>>> 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 >>>>>> >>>>>> >>>>>> >>>>>> >>>>>> ------------------------------------------------------------------------- >>>>>> 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 >>>>>> >>>>> >>>> >>> >> > |
|
From: Bill H. <goo...@go...> - 2008-11-21 13:48:38
|
I've started rewriting the F_mpz_poly module based on the new F_mpz module (which now has quite a few functions again). I've coded up F_mpz_poly_add. I ran my benchmarks to see how fast it is. The results are mixed. The benchmark involves constructing poly's A, B, C with signed coefficients and the same number of bits per coefficient and the same length as given in the tables below (except in the final case where C has a different number of bits per coefficient to A and B, as indicated). There are 3 benchmarks run for each problem: 1) Clear and initialise a polynomial R of the given length 2) Do 1 and compute R = A + B 3) Do 2 and compute R = R + C I've included times for the existing fmpz_poly from the FLINT 1.x series, times for FLINT 1.x mpz_poly and the new times for the new FLINT 2.x F_mpz_poly. The good thing is it is always faster than mpz_poly. The most interesting thing is that it is faster to add polynomials with 100 bits per coefficient than with 64 bits. Moreover, FLINT 1.x is much faster at the 64 bit problem. In fact, FLINT 2.x is nearly as slow as mpz_poly. I really don't have an explanation for this. In some senses 64 bits is a worst case for FLINT 2.x. It has to sometimes use 65 bits which means 2 limbs in the mpz_t's (though these are reused and should stabilise at 2 limbs very quickly) and sometimes it has to drop back to 62 bits when there is cancellation, which means repeatedly promoting and demoting coefficients. However when doing R = A + B over and over, I have it set up to do 10000 iterations with the same A and B before changing A and B, so there really ought not be many promotions and demotions between 62 bit coefficients and mpz_t's. Besides all this, it is nearly as slow for 70 bits for which both problems are avoided. The only explanation I have is that GMP is itself very slow at adding 70 bit signed integers. However not even this makes sense. After all, why would it be faster to add 100 bit integers!? It's certainly not a caching issue. So I'm a bit stumped. If anyone has any guesses, please let me know.... Anyhow, here are the times (which I am not unhappy with): length 100, bits = 50, iter = 1000000 ===================================== mpz_poly: clear-init2-add-add : 24.342s clear-init2-add : 17.957s clear-init2 : 9.227s flint fmpz_poly: clear-init-add-add : 10.522s clear-init-add : 5.059s clear-init : 0.183s flint-2 F_mpz_poly: clear-init-add-add : 2.360s clear-init2-add : 1.320s clear-init2 : 0.400s length 100, bits = 64, iter = 1000000 ===================================== mpz_poly: clear-init2-add-add : 24.608s clear-init2-add : 17.758s clear-init2 : 9.245s flint fmpz_poly: clear-init-add-add : 9.861s clear-init-add : 3.804s clear-init : 0.183s flint-2 F_mpz_poly: clear-init-add-add : 23.290s clear-init2-add : 17.259s clear-init2 : 0.450s length 100, bits = 1000, iter = 1000000 ===================================== mpz_poly: clear-init2-add-add : 39.639s clear-init2-add : 28.760s clear-init2 : 10.445s flint fmpz_poly: clear-init-add-add : 19.321s clear-init-add : 8.697s clear-init : 0.235s flint-2 F_mpz_poly: clear-init-add-add : 19.320s clear-init2-add : 9.780s clear-init2 : 0.500s length 100, bits = 100, 200, iter = 1000000 ===================================== mpz_poly: clear-init2-add-add : 33.716s clear-init2-add : 18.921s clear-init2 : 10.237s flint fmpz_poly: clear-init-add-add : 14.611s clear-init-add : 6.425s clear-init : 0.190s flint-2 F_mpz_poly: clear-init-add-add : 13.800s clear-init2-add : 7.220s clear-init2 : 0.450s Bill. |
|
From: Bill H. <goo...@go...> - 2008-11-21 15:43:18
|
I found the problem with F_mpz_poly_add. Basically when coefficients were being demoted as a result of cancellations, the mpz_t's were being lost instead of returned to the memory manager. The times are now much more what I would expect, and close to the FLINT 1.x series times. So now FLINT has a sensible integer format and writing polynomial functions is quite easy. The other advantage of the new model is that attaching polynomials using F_mpz_poly_attach will still be possible. It wasn't under the previous model. length 100, bits = 50, iter = 1000000 ===================================== mpz_poly: clear-init2-add-add : 24.342s clear-init2-add : 17.957s clear-init2 : 9.227s flint fmpz_poly: clear-init-add-add : 10.522s clear-init-add : 5.059s clear-init : 0.183s flint-2 F_mpz_poly: clear-init-add-add: 2.360s clear-init2-add: 1.320s clear-init2: 0.400s length 100, bits = 64, iter = 1000000 ===================================== mpz_poly: clear-init2-add-add : 24.608s clear-init2-add : 17.758s clear-init2 : 9.245s flint fmpz_poly: clear-init-add-add : 9.861s clear-init-add : 3.804s clear-init : 0.183s flint-2 F_mpz_poly: clear-init-add-add: 11.740s clear-init2-add: 7.240s clear-init2: 0.450s length 100, bits = 1000, iter = 1000000 ===================================== mpz_poly: clear-init2-add-add : 39.639s clear-init2-add : 28.760s clear-init2 : 10.445s flint fmpz_poly: clear-init-add-add : 19.321s clear-init-add : 8.697s clear-init : 0.235s flint-2 F_mpz_poly: clear-init-add-add: 18.940s clear-init2-add: 9.240s clear-init2: 0.500s length 100, bits = 100, 200, iter = 1000000 ===================================== mpz_poly: clear-init2-add-add : 33.716s clear-init2-add : 18.921s clear-init2 : 10.237s flint fmpz_poly: clear-init-add-add : 14.611s clear-init-add : 6.425s clear-init : 0.190s flint-2 F_mpz_poly: clear-init-add-add: 12.160s clear-init2-add: 6.520s clear-init2: 0.450s Bill. |
|
From: Bill H. <goo...@go...> - 2008-11-21 13:49:35
|
Apparently the message below never made it to the list, so here it is again. ---------- Forwarded message ---------- From: William Hart <ha...@ya...> Date: 2008/11/18 Subject: New F_mpz_t integer arithmetic modular To: Bill Hart <goo...@go...> I've again been implementing a FLINT integer type (F_mpz_t). This time an F_mpz is either a long (with 62 bits and a sign) or an index into a FLINT wide array of mpz_t's. Then an F_mpz_t is just an array of length 1 of F_mpz's so that call by reference can be implemented (and checking for aliasing is easy). The one downside of this is that if one wants to make this stuff multi-threaded, one needs an array of mpz_t's for each thread. But this shouldn't be too hard to manage. Originally I tried letting an F_mpz_t be directly a long or index into an array of mpz_t's. But then you have to have an interface of functions that look like F_mpz_set(F_mpz_t * f, const F_mpz_t g) and checking for aliasing is then more difficult. Of course the array of length 1 business really ends up causing an extra layer of indirection (especially bad when an integer is already being represented internally as an mpz_t), but the timings don't seem to be affected, as presumably the compiler adapts this extra layer away most of the time. Here are the timings before this extra indirection: Testing F_mpz_getset_ui()... Cpu = 1220 ms Wall = 1225 ms ... ok Testing F_mpz_getset_si()... Cpu = 1240 ms Wall = 1232 ms ... ok Testing F_mpz_getset_mpz()... Cpu = 3600 ms Wall = 3600 ms ... ok Testing F_mpz_set()... Cpu = 180 ms Wall = 189 ms ... ok Here are the timings with the extra indirection: Testing F_mpz_getset_ui()... Cpu = 1220 ms Wall = 1224 ms ... ok Testing F_mpz_getset_si()... Cpu = 1280 ms Wall = 1279 ms ... ok Testing F_mpz_getset_mpz()... Cpu = 3590 ms Wall = 3596 ms ... ok Testing F_mpz_set()... Cpu = 190 ms Wall = 192 ms ... ok That's despite amortising away the cost of random generation and setup cost in the timings and working with very small integers, of say no more than 200 bits. So I have decided to stick with the nice consistent interface that allows easy checks for aliasing but has the extra level of indirection. I've implemented all the memory management stuff and a handful of basic functions and a few test functions. It is implemented in F_mpz.c, F_mpz.h and F_mpz-test.c. Just do make F_mpz-test to make the test module. So now FLINT has it's own integer format. Yay!! Bill. |