This is an Opt In Archive . We would like to hear from you if you want your posts included. For the contact address see About this archive. All posts are copyright (c).
- Contents - Hide Contents - Home - Section 3First Previous Next
2000 2050 2100 2150 2200 2250 2300 2350 2400 2450 2500 2550 2600 2650 2700 2750 2800 2850 2900 2950
2000 - 2025 -
Message: 2000 - Contents - Hide Contents Date: Sun, 18 Nov 2001 06:21:27 Subject: Re: LLL definitions From: genewardsmith@xxxx.xxx --- In tuning-math@y..., "Paul Erlich" <paul@s...> wrote:> --- In tuning-math@y..., genewardsmith@j... wrote:>> In the finite-dimensional case, they define the same topology; > hence>> the identity map is continuous.> So that's promising for being able to use the Tenney metric > eventually, yes?So far as I can see we can't use it unless we dig up an L1 version of LLL, but it does give grounds for hope the answers would often be the same.> OK -- so you should have said "for all g" in the definition . . . > yes? Or is g uniquely defined? It would seem to depend on how you > order the elements of b, since the first vector is made to be the > same.It does depend on the fact that b is ordered, which is why we started with an ordered basis. It is uniquely defined given the ordering.> Can we determine for certain whether we do or don't in any particular > case?I think in such small cases that might be possible; in general this turns into a non-polynomial time mess.>> Existence we do have, >> since LLL will work for this range of d.> I was thinking that two different choices of b might turn out to be > equal contenders for the LLL throne for a particular equivalence > class of b's. Possible or impossible?In theory possible, in practice I think we may be OK. We are not working in the 101 limit, after all.
Message: 2001 - Contents - Hide Contents Date: Sun, 18 Nov 2001 06:24:48 Subject: Re: LLL basis for linear temperaments From: genewardsmith@xxxx.xxx --- In tuning-math@y..., "Paul Erlich" <paul@s...> wrote:> Virtually all of my discussions with you have been geared toward this > one project -- start by motivating periodicity blocks, continue by > motivating temperament of smaller unison vectors, and conclude by > motivating ETs and linear temperaments.Or we could start from ets and detemper to linear and planar temperaments. I think it's good to point out the duality between ets and unison vectors.
Message: 2002 - Contents - Hide Contents Date: Sun, 18 Nov 2001 06:30:08 Subject: Re: LLL definitions From: Paul Erlich --- In tuning-math@y..., genewardsmith@j... wrote:> --- In tuning-math@y..., "Paul Erlich" <paul@s...> wrote:>> --- In tuning-math@y..., genewardsmith@j... wrote: >>>> In the finite-dimensional case, they define the same topology; >> hence>>> the identity map is continuous. >>> So that's promising for being able to use the Tenney metric >> eventually, yes? >> So far as I can see we can't use it unless we dig up an L1 version of > LLL, but it does give grounds for hope the answers would often be the > same.But you said, "if we want another inner product for LLL reduction one way to do it is to transform coordinates and then transform them back." Is there an inner product that is appropriate to the taxicab metric? Am I confused about something?>> OK -- so you should have said "for all g" in the definition . . . >> yes? Or is g uniquely defined? It would seem to depend on how you >> order the elements of b, since the first vector is made to be the >> same. >> It does depend on the fact that b is ordered, which is why we started > with an ordered basis.Hmm . . . shouldn't we then try all permutations of a given basis? Or do the elements have to be in descending order of length? I don't think this is assumed for Gram-Schmidt generally, is it?>> Can we determine for certain whether we do or don't in any > particular >> case? >> I think in such small cases that might be possible; in general this > turns into a non-polynomial time mess.Well, it would be cool to be sure of uniqueness.>>> Existence we do have, >>> since LLL will work for this range of d. >>> I was thinking that two different choices of b might turn out to be >> equal contenders for the LLL throne for a particular equivalence >> class of b's. Possible or impossible? >> In theory possible, in practice I think we may be OK. We are not > working in the 101 limit, after all.Well, having equal contenders us something that, to my intuition, could happen just as easily in few dimensions as in many . . . no?
Message: 2003 - Contents - Hide Contents Date: Sun, 18 Nov 2001 06:35:28 Subject: Re: LLL basis for linear temperaments From: Paul Erlich --- In tuning-math@y..., genewardsmith@j... wrote:> --- In tuning-math@y..., "Paul Erlich" <paul@s...> wrote: >>> Virtually all of my discussions with you have been geared toward > this>> one project -- start by motivating periodicity blocks, continue by >> motivating temperament of smaller unison vectors, and conclude by >> motivating ETs and linear temperaments. >> Or we could start from ets and detemper to linear and planar > temperaments.Many of the readers, coming from a JI background, would be most unhappy if we started from ETs rather than from JI. I'd like the paper to be capable of convincing one or two people as to the value of ETs and LTs. Meanwhile, most modern academic microtonal music theory starts from ETs and MOSs, and do not motivate this starting point correctly or at all. I'd like to provide that foundation. Oh yeah, did I mention the Hypothesis would be part of this paper?> I think it's good to point out the duality between ets > and unison vectors.The duality might be a bit too abstract for musicians to grasp or even care about. And many of the readers, coming from a JI background, would be most unhappy if we started from ETs rather than from JI. I'd like the paper to be capable of convincing one or two people as to the value of ETs. Meanwhile, most modern academic microtonal music theory starts from ETs, and do not motivate this starting point correctly or at all.
Message: 2004 - Contents - Hide Contents Date: Sun, 18 Nov 2001 07:25:44 Subject: LLL reduction of pairs of 7-limit commas From: genewardsmith@xxxx.xxx Paul gave a list previously which was a little more inclusive: 25/24, 28/27, 36/35, 49/48, 50/49, 64/63, 81/80, 126/125, 245/243, 225/224, 1029/1024, 1728/1715, 2401/2400, 3136/3125, 4375/4374, 6144/6125 I did an LLL reduction on all the 16 choose 2 = 120 pairs of elements from this list, and obtained the following 94 reduced LLL-reduced pairs: [36/35, 50/49] [126/125, 1029/1024] [25/24, 48/49] [49/48, 126/125] [64/63, 4375/4374] [81/80, 875/864] [225/224, 4375/4374] [64/63, 686/675] [28/27, 3024/3125] [28/27, 25/24] [3136/3125, 1029/1024] [126/125, 245/243] [126/125, 81/80] [126/125, 64/63] [25/24, 49/45] [49/48, 392/375] [2401/2400, 3136/3125] [126/125, 1715/1728] [1728/1715, 3136/3125] [225/224, 49/48] [28/27, 126/125] [225/224, 1029/1024] [225/224, 1728/1715] [225/224, 245/243] [225/224, 81/80] [225/224, 64/63] [50/49, 128/125] [25/24, 392/405] [4375/4374, 6144/6125] [4375/4374, 3136/3125] [24/25, 49/48] [21/20, 27/28] [64/63, 50/49] [36/35, 200/189] [25/24, 64/63] [126/125, 1024/1029] [81/80, 50/49] [81/80, 49/48] [15/14, 36/35] [4375/4374, 2401/2400] [1728/1715, 4000/3969] [50/49, 4375/4374] [27/28, 49/48] [245/243, 2048/2025] [64/63, 243/245] [245/243, 3136/3125] [25/24, 360/343] [28/27, 625/672] [36/35, 1323/1280] [4375/4374, 1029/1024] [36/35, 225/224] [28/27, 3087/3200] [36/35, 128/135] [28/27, 48/49] [21/20, 24/25] [36/35, 49/50] [50/49, 525/512] [36/35, 16/15] [81/80, 3136/3125] [36/35, 63/64] [21/20, 35/36] [25/24, 81/80] [81/80, 1728/1715] [81/80, 6144/6125] [25/24, 27/28] [28/27, 6144/6125] [49/48, 3136/3125] [50/49, 6144/6125] [81/80, 1715/1728] [25/24, 3087/3200] [36/35, 64/63] [49/48, 80/81] [28/27, 49/48] [64/63, 245/243] [245/243, 3969/4000] [49/48, 5000/5103] [25/24, 224/243] [25/24, 2401/2400] [64/63, 3136/3125] [81/80, 2401/2400] [64/63, 864/875] [6144/6125, 2401/2400] [20/21, 36/35] [36/35, 1029/1000] [25/24, 49/48] [36/35, 5488/5625] [28/27, 1029/1024] [245/243, 1029/1024] [126/125, 1728/1715] [28/27, 50/49] [28/27, 256/245] [50/49, 245/243] [50/49, 864/875] [25/24, 160/147] Some of these have a big difference between the large and small commas, and some introduce commas not on the orginal list.
Message: 2005 - Contents - Hide Contents Date: Sun, 18 Nov 2001 08:26:14 Subject: LLL reduction pairs revised list From: genewardsmith@xxxx.xxx My last list allowed in some duplicates; the following list of 81 pairs is from an improved version of the program: {135/128, 36/35} {36/35, 200/189} {405/392, 25/24} {81/80, 126/125} {4000/3969, 1728/1715} {64/63, 245/243} {25/24, 21/20} {28/27, 50/49} {36/35, 225/224} {160/147, 25/24} {50/49, 81/80} {5103/5000, 49/48} {28/27, 49/48} {64/63, 126/125} {1728/1715, 3136/3125} {128/125, 50/49} {672/625, 28/27} {81/80, 3136/3125} {3200/3087, 28/27} {1029/1024, 4375/4374} {245/243, 2048/2025} {3136/3125, 4375/4374} {225/224, 1029/1024} {50/49, 64/63} {49/48, 392/375} {28/27, 256/245} {126/125, 245/243} {50/49, 245/243} {50/49, 875/864} {50/49, 6144/6125} {25/24, 49/48} {245/243, 225/224} {50/49, 4375/4374} {25/24, 64/63} {25/24, 49/45} {126/125, 1029/1024} {225/224, 4375/4374} {25/24, 2401/2400} {81/80, 6144/6125} {81/80, 1728/1715} {64/63, 686/675} {36/35, 1323/1280} {49/48, 3136/3125} {64/63, 4375/4374} {28/27, 126/125} {49/48, 225/224} {3200/3087, 25/24} {81/80, 225/224} {245/243, 1029/1024} {245/243, 3136/3125} {36/35, 50/49} {28/27, 6144/6125} {25/24, 243/224} {4000/3969, 245/243} {50/49, 525/512} {36/35, 64/63} {28/27, 1029/1024} {1029/1024, 3136/3125} {25/24, 81/80} {49/48, 81/80} {2401/2400, 3136/3125} {81/80, 2401/2400} {49/48, 126/125} {64/63, 875/864} {36/35, 1029/1000} {25/24, 28/27} {28/27, 3125/3024} {2401/2400, 6144/6125} {36/35, 5625/5488} {36/35, 21/20} {81/80, 875/864} {36/35, 16/15} {25/24, 360/343} {126/125, 1728/1715} {36/35, 15/14} {225/224, 1728/1715} {64/63, 225/224} {64/63, 3136/3125} {4375/4374, 6144/6125} {2401/2400, 4375/4374} {28/27, 21/20}
Message: 2006 - Contents - Hide Contents Date: Sun, 18 Nov 2001 09:48:59 Subject: Re: LLL reduction pairs revised list From: genewardsmith@xxxx.xxx --- In tuning-math@y..., genewardsmith@j... wrote:> {126/125, 1728/1715}I selected the above at random from my list, and came to the following conclusions: (1) It should be discussed in conjunction with the 31+27;58 generator in the 58-et, and (2) It should be discussed in conjunction with a corresponding 11-limit planar temperament. Is the world of JI ready for this, I wonder? Paul made the music theory world sound a bit like the big-endians vs the little-endians.
Message: 2007 - Contents - Hide Contents Date: Sun, 18 Nov 2001 09:53:21 Subject: Re: LLL reduction pairs revised list From: genewardsmith@xxxx.xxx --- In tuning-math@y..., genewardsmith@j... wrote:> {126/125, 1728/1715}I selected the above at random from my list, and came to the following conclusions: (1) It should be discussed in conjunction with the 31+27;58 and 27+4;31, and (2) It should be discussed in conjunction with a corresponding 11-limit linear temperament. Is the world of JI ready for this, I wonder? Paul made the music theory world sound a bit like the big-endians vs the little-endians.
Message: 2008 - Contents - Hide Contents Date: Sun, 18 Nov 2001 10:31:30 Subject: Minkowski reduced lattices From: genewardsmith@xxxx.xxx The classic Minkowski reduced lattice definition in the case of a two- dimensional lattice reduces to saying that one basis element is as short as possible, and the other is the shortest such that together with the first it generates the lattice. In small dimensions like this, it is possible to reduce lattices without worrying about exponential time. The definition of Minkowski reduced can also be applied using the Tenney metric.
Message: 2009 - Contents - Hide Contents Date: Mon, 19 Nov 2001 03:35:41 Subject: Re: Minkowski reduced lattices From: Paul Erlich --- In tuning-math@y..., genewardsmith@j... wrote:> The classic Minkowski reduced lattice definition in the case of a two- > dimensional lattice reduces to saying that one basis element is as > short as possible, and the other is the shortest such that together > with the first it generates the lattice. In small dimensions like > this, it is possible to reduce lattices without worrying about > exponential time. The definition of Minkowski reduced can also be > applied using the Tenney metric.Minkowski reduced basis . . . it it unique?
Message: 2010 - Contents - Hide Contents Date: Mon, 19 Nov 2001 21:27:14 Subject: Re: LLL reduction pairs revised list From: Paul Erlich --- In tuning-math@y..., genewardsmith@j... wrote:>Paul gave a list previously which was a little more inclusive: >25/24, 28/27, 36/35, 49/48, 50/49, 64/63, 81/80, 126/125, 245/243, >225/224, 1029/1024, 1728/1715, 2401/2400, 3136/3125, 4375/4374, >6144/6125Well, Gene, the first four entries were intended to be used only as chromatic unison vectors. Actually, we should be taking three at a time, not two. If all three (in the Minkowski reduced basis, hopefully) are commatic UVs, we get an ET. If one is chromatic, we get an MOS. If two are chromatic, a planar temperament. If three are chromatic, a JI block. I think the UVs in the reduced basis have to be in Kees's list or something like it -- if they aren't, then it seems likely that some step in the scale will be smaller than one of the commatic unison vectors -- which we shouldn't allow. Just some imprecise thoughts . . .
Message: 2011 - Contents - Hide Contents Date: Mon, 19 Nov 2001 03:52:31 Subject: Re: Minkowski reduced lattices From: genewardsmith@xxxx.xxx --- In tuning-math@y..., "Paul Erlich" <paul@s...> wrote:> Minkowski reduced basis . . . it it unique?Yes, that's basically why I brought it up. I think usually the LLL- reduced basis in the 7-limit will also be Minkowski reduced, but if it isn't it would be easy enough to fix. Here's some questions for you: how did you pick the commas you wanted to start out from, both this time and the last time? What would be good cut-off criteria for exluding a basis--for instance, numerator times denominator > 2000, let us say, and small basis element not less than a third of the large one?
Message: 2012 - Contents - Hide Contents Date: Mon, 19 Nov 2001 21:39:50 Subject: Re: Minkowski reduced lattices From: Paul Erlich --- In tuning-math@y..., genewardsmith@j... wrote:> --- In tuning-math@y..., genewardsmith@j... wrote: >>> I'll get back to you if I can find this conjecture you speak of. >> I think you either need to find a better keyboardDone. The thinking goes as follows. The "length" of a unison vector n/d, n~=d, in the Tenney lattice with taxicab metric, or van Prooijen lattice with triangular-taxicab metric, is proportional to log(n) + log(d) (hence approx. proportional to log(d)), and also to the "number" (in some weighted sense) of consonant intervals making up that unison vector. Thus, in order to temper this unison vector out (assuming that other UVs being tempered out, if any, are orthogonal to this one), one must temper each consonant interval involved by an "average" amount proportional to w/log(d), where w is the musical width of the unison vector. w=log(n/d) w~=n/d-1 w~=(n-d)/d Hence the amount of tempering implied by the unison vector is approx. proportional to (n-d)/(d*log(d)) Yes?
Message: 2013 - Contents - Hide Contents Date: Mon, 19 Nov 2001 04:36:40 Subject: Re: Minkowski reduced lattices From: Paul Erlich --- In tuning-math@y..., genewardsmith@j... wrote:> --- In tuning-math@y..., "Paul Erlich" <paul@s...> wrote: >>> Minkowski reduced basis . . . it it unique? >> Yes, that's basically why I brought it up. I think usually the LLL- > reduced basis in the 7-limit will also be Minkowski reducedOn the Tenney lattice with a taxicab metric?>, but if > it isn't it would be easy enough to fix. > > Here's some questions for you: how did you pick the commas you wanted > to start out from, both this time and the last time? What would be > good cut-off criteria for exluding a basis--for instance, numerator > times denominator > 2000, let us say, and small basis element not > less than a third of the large one?I'll have to think hard about this. For now, does the conjecture I've posted about the relationship between the numbers in a comma's ratio, and the amount of tempering the comma implied for the constituent consonant intervals, make any sense to you? This keyboard won't allow me to type this right now, but you can search for "conjecture" . . .
Message: 2014 - Contents - Hide Contents Date: Mon, 19 Nov 2001 23:21:42 Subject: "most compact" periodicity block From: Paul Erlich Is it clear what I mean when I say "most compact" periodicity block for a particular equivalence class of bases? Won't it, in general, be delimited by a hexagon in the 5-limit, and a rhombic dodecahedron in the 7-limit? Could we use the resulting three or six unison vectors as a more complete characterization of a given system, rather than an LLL or some such reduced basis, which has some element of arbitrariness because some "second best" reduction might be nearly as good?
Message: 2015 - Contents - Hide Contents Date: Mon, 19 Nov 2001 06:23:42 Subject: Re: Minkowski reduced lattices From: genewardsmith@xxxx.xxx --- In tuning-math@y..., "Paul Erlich" <paul@s...> wrote:>> Yes, that's basically why I brought it up. I think usually the LLL- >> reduced basis in the 7-limit will also be Minkowski reduced> On the Tenney lattice with a taxicab metric?I think so--we are talking about only two vectors, after all, and a Euclidean distance adjusted to be as much like Tenney as possible. I'll get back to you if I can find this conjecture you speak of.
Message: 2016 - Contents - Hide Contents Date: Mon, 19 Nov 2001 23:44:41 Subject: Re: Minkowski reduced lattices From: genewardsmith@xxxx.xxx --- In tuning-math@y..., "Paul Erlich" <paul@s...> wrote:> Yes?Sounds reasonable. Is this a conjecture, and observation or a heuristic argument? If it's a conjecture it needs to be stated more precisely.
Message: 2017 - Contents - Hide Contents Date: Mon, 19 Nov 2001 07:09:28 Subject: Re: Minkowski reduced lattices From: genewardsmith@xxxx.xxx --- In tuning-math@y..., genewardsmith@j... wrote:> I'll get back to you if I can find this conjecture you speak of.I think you either need to find a better keyboard, or give me an article number.
Message: 2018 - Contents - Hide Contents Date: Mon, 19 Nov 2001 23:47:10 Subject: Re: Minkowski reduced lattices From: Paul Erlich --- In tuning-math@y..., genewardsmith@j... wrote:> --- In tuning-math@y..., "Paul Erlich" <paul@s...> wrote: > >> Yes? >> Sounds reasonable. Is this a conjecture, and observation or a > heuristic argument? If it's a conjecture it needs to be stated more > precisely.Heuristic argument. If it can be made more precise I'd be most grateful to whoever does so.
Message: 2019 - Contents - Hide Contents Date: Mon, 19 Nov 2001 23:48:58 Subject: Re: "most compact" periodicity block From: genewardsmith@xxxx.xxx --- In tuning-math@y..., "Paul Erlich" <paul@s...> wrote:> Is it clear what I mean when I say "most compact" periodicity block > for a particular equivalence class of bases? Won't it, in general, be > delimited by a hexagon in the 5-limit, and a rhombic dodecahedron in > the 7-limit?It seems to me that depends on how you define your distance. Incidentally, I suspect that instead of using hexagonal blocks, or rrhombic dodecahedra, ellipsoids would work. I don't think you need a tiling. Could we use the resulting three or six unison vectors> as a more complete characterization of a given system, rather than an > LLL or some such reduced basis, which has some element of > arbitrariness because some "second best" reduction might be nearly as > good?What's the point? I thought you wanted to produce temperaments, not PBs.
Message: 2020 - Contents - Hide Contents Date: Mon, 19 Nov 2001 23:55:03 Subject: Re: LLL reduction pairs revised list From: genewardsmith@xxxx.xxx --- In tuning-math@y..., "Paul Erlich" <paul@s...> wrote:> Actually, we should be taking three at a time, not two. If all three > (in the Minkowski reduced basis, hopefully) are commatic UVs, we get > an ET. If one is chromatic, we get an MOS. If two are chromatic, a > planar temperament. If three are chromatic, a JI block.This will prevent you from considering a lot of interesting temperaments.> I think the UVs in the reduced basis have to be in Kees's list or > something like it -- if they aren't, then it seems likely that some > step in the scale will be smaller than one of the commatic unison > vectors -- which we shouldn't allow.Where is Kees's list? I was thinking of producing a list based on some mathematical conditions, and using that, but I wondered if these lists already represent such an effort.
Message: 2021 - Contents - Hide Contents Date: Tue, 20 Nov 2001 19:56:30 Subject: Re: LLL definitions From: genewardsmith@xxxx.xxx --- In tuning-math@y..., graham@m... wrote:> What's an ordered basis?It's a basis you take in a certain order--first element, second element, etc.>> c[i,j] = <b_i, g_j>/<g_j, g_j>> I've implemented this in Python code: > Please say if I'm going wrong anywhere.It looks right to me, but I don't know Python. One thing that worries me is that> I'm getting floating point results where you stay with integers.If you start with integers, you get rational numbers, and hence you will get floating points in your implementation of Gram-Schmidt. The integers I am staying with are simply the integers in a reduced basis that starts out with integers.> Again, I've implemented that > I'm assuming that repeatedly applying this will give the right result, but > it clearly doesn't.What I gave was a definition of LLL reduced, not the LLL algorithm itself; that has two different kinds of steps in it. Here is some code I took off the web, which comes with a reference. Starting with the pseudocode in a reference book on algorithms, or trying to adapt some code in another implementation, seems like the way to go to me. # This follows the algorithm given in Maurice Mignotte's # "Mathematics for Computer Algebra" LLL := proc() local dim, basis, new_basis, local_iprod, local_norm, local_c, i, j, u, B, k, h, l; basis := args[1]; dim := nops(basis); local_iprod := linalg[innerprod]; local_norm := proc(v::vector) RETURN(linalg[norm](v,2)); end: local_c := 0.75; for i from 2 to nargs do if op(1,args[i]) = iprod then local_iprod := op(2,args[i]); local_norm := proc(v::vector) RETURN(sqrt(local_iprod(v,v))); end; elif op(1,args[i]) = norm then local_norm := op(2,args[i]); elif op(1,args[i]) = c then local_c := op(2,args[i]); fi; od; For one thing, the first ratio can never change.> Also, you haven't said what seed value to use for g. I'm using the > initial value of b.g doesn't need a seed value, it's just the Gram-Schmidt orthogonalization of b.
Message: 2022 - Contents - Hide Contents Date: Tue, 20 Nov 2001 20:39:19 Subject: Another LLL implementation From: genewardsmith@xxxx.xxx This is from Calc; see 404 Not Found * [with cont.] Search for http://www.maths.uq.edu.au/~krm/krm_calc.html#[49 in Wayback Machine] /* The Pohst Algorithm updating only from where necessary */ #include <stdio.h> #include <stdlib.h> #include "integer.h" #include "fun.h" extern MPI *MAXI, *PMAXI; extern unsigned int MLLLVERBOSE; extern unsigned int HERMITEVERBOSE; unsigned int GCDFLAG; MPMATI *BASIS_REDUCTION(MPMATI *Bptr, MPMATI **Eptr, USI rowstage, USI m1, USI n1) /* * Input: *Bptr, a matrix of MPI's, whose first row is not zero. * Output: a pointer to an MPMATI whose rows form a reduced basis for * the lattice spanned by the rows of *Bptr. This basis is reduced in the * sense of the paper "Factoring polynomials with rational coefficients" by * A. K. Lenstra, H. W. Lenstra and L. Lovasz, Math. Ann. 261, 515- 534 (1982) * using the modified version in "Solving exponential Diophantine equations * using lattice basis reduction algorithms" by B. M. M. De Weger, J. No. Theory * 26, 325-367 (1987). A change of basis matrix **Eptr is also returned. * De Weger's algorithm has been changed to cater for arbitrary matrices. The * the rows are now in general linearly dependent. * We use the fact that the Gram Schmidt process detects the first row * which is a linear combination of the preceding rows. We employ a modification * of the LLL algorithm outlined by M. Pohst in J. Symbolic Computation (1987)4, * 123-127. We call this the MLLL algorithm. * The last sigma rows of the matrix **Eptr are relation vectors. * m1 / n1 is usually taken to be 3 / 4. */ { unsigned int i, k, l, n, m, t, flag = 0, Flag = 0; unsigned int flagg, beta, K1 = 0, tau = 2, sigma = 0, rho; MPI **D, *X, *Y, *Z, *H, *Tmp, *R, *M1, *N1; MPMATI *C, *L, *B1ptr; m = Bptr->C; n = Bptr->R; /* We initial Eptr outside the function whenever we call the function. */ /* This is because we have to do so in SMITH(). */ /* MAXI = MAXELTI(Bptr); PMAXI = MAXELTI(*Eptr); */ B1ptr = COPYMATI(Bptr); D = (MPI **)mmalloc((1 + n) * sizeof(MPI *)); D[0] = ONEI(); for (i = 1; i <= n; i++) D[i] = ZEROI(); C = ZEROMNI(n, m); L = ZEROMNI(n, n); found: n = B1ptr->R; i = (K1 == 0) ? 1 : K1; /* K1 = no. of consecutive rows of *B1ptr that don't need updating for the Gram Schmidt process. */ while (i <= B1ptr->R) { BASIS_UPDATE(i, m, &C, &L, B1ptr, D); flag = 1; for (t = 0; t < m; t++) { if (!EQZEROI(C->V[i - 1][t])) { flag = 0; break; } } if (flag) break; X = ZEROI(); for (t = 0; t < m; t++) { H = MULTI(C->V[i - 1][t], C->V[i - 1][t]); Tmp = X; X = ADDI(X, H); FREEMPI(Tmp); FREEMPI(H); } FREEMPI(D[i]); D[i] = INT(X, D[i - 1]); FREEMPI(X); i++; if (MLLLVERBOSE) printf("i = %u\n", i); } beta = (flag) ? i : i - 1; rho = K1 = i - 1; if (MLLLVERBOSE) printf("completed updating the basis\n"); /* Here K1 = no. of LI rows in *B1ptr found by Gram Schmidt process. flag = 0 means all the rho = beta rows of *B1ptr are LI; flag = 1 means that the first rho = beta - 1 rows of *B1ptr are LI, but the beta-th row is a LC of the preceding rows. So beta = number of rows of *B1ptr currently being examined by the LLL algorithm. */ k = tau; if (MLLLVERBOSE) printf("beta = %u\n", beta); M1 = CHANGE(m1); N1 = CHANGE(n1); while (k <= beta) { if (MLLLVERBOSE) printf("beta - k = %u\n", beta -k); l = k - 1; Flag = STEP4(k, l, &L, &B1ptr, Eptr, D, rowstage); if (Flag)/* STEP 9 of POHST. */ { sigma++; if (MLLLVERBOSE) printf("relation vector number %u found\n", sigma); tau = k++; goto found; } X = MULTI(D[k - 2], D[k]); Y = MULTI(D[k - 1], D[k - 1]); Tmp = Y; Y = MULT_I(Y, m1); FREEMPI(Tmp); R = MULTI(L->V[k - 1][k - 2], L->V[k - 1][k - 2]); Z = ADD0I(X, R); Tmp = Z; Z = MULT_I(Z, n1); FREEMPI(Tmp); if (RSV(Y, Z) == 1)/*& STEP 5 of POHST. */ { flagg = 0; if (EQZEROI(D[k]) && EQZEROI(R)) {/* CASE B=0 of STEP 7 of POHST. */ FREEMPI(D[k - 1]); D[k - 1] = ZEROI(); STEP8(k, &B1ptr, &L, Eptr, rowstage); if (k - 1 < K1) K1 = k - 1; /* The swap may have changed 2nd last row */ /* of *B1ptr. */ for (t = 0; t < m; t++) { FREEMPI(C->V[k - 2][t]); C->V[k - 2][t] = ZEROI(); } beta--; flagg = 1; if (k > 2) k--; FREEMPI(X); FREEMPI(Y); FREEMPI(R); FREEMPI(Z); continue; } if (flagg == 0) { for (i = k + 1; i <= beta; i++) STEP7(i, k, &L, D); } STEP8(k, &B1ptr, &L, Eptr, rowstage); if (k - 2 < K1) K1 = k - 2; /* swap will change last two rows of *B1ptr. */ FREEMPI(R); FREEMPI(Y); Y = MULTI(L->V[k - 1][k - 2], L->V[k - 1][k - 2]); Tmp = Y; Y = ADD0I(Y, X); FREEMPI(X); FREEMPI(Tmp); Tmp = D[k - 1]; D[k - 1] = INT0(Y, D[k - 1]); FREEMPI(Tmp); FREEMPI(Y); if (k > 2) k--; } else { /* STEP 6 of POHST. */ FREEMPI(R); FREEMPI(X); FREEMPI(Y); for (l = k - 2; l >= 1; l--) { Flag = STEP4(k, l, &L, &B1ptr, Eptr, D, rowstage); if (Flag) { FREEMPI(Z); sigma++; if (MLLLVERBOSE) printf("relation vector number %u found\n", sigma); tau = k++; goto found; /* STEP 9 of POHST. */ } } k++; } FREEMPI(Z); } FREEMPI(M1); FREEMPI(N1); FREEMATI(C); printf("L = \n"); PRINTMATI(0,L->R-1,0,L->C-1,L); for (i = 0; i <= Bptr->R; i++) { printf("D[%u] = ", i);PRINTI(D[i]);printf(", "); } printf("\n"); FREEMATI(L); for (i = 0; i <= Bptr->R; i++) FREEMPI(D[i]); ffree((char *)D, (1 + Bptr->R) * sizeof(MPI *)); if (MLLLVERBOSE) { printf("number of basis vectors found = %u ;\n", rho); printf("number of relation vectors found = %u .\n", sigma); } return (B1ptr); } unsigned int STEP4(k, l, Lptr, Bptr, Eptr, D, i) /* * updates *Lptr, *Bptr and *Eptr. * returns 1 if row k of *Bptr becomes zero, returns zero otherwise. */ unsigned int k, l, i; MPI *D[]; MPMATI **Lptr, **Bptr, **Eptr; { unsigned int j, flag = 1, t, m, n; MPI *X, *Y, *R, *Tmp; MPMATI *TmpMATI; m = (*Bptr)->C; n = (*Eptr)->R; Y = MULT_I((*Lptr)->V[k - 1][l - 1], 2); if (RSV(Y, D[l]) == 1) { R = NEAREST_INTI((*Lptr)->V[k - 1][l - 1], D[l]); X = MINUSI(R); TmpMATI = *Bptr; *Bptr = ADD_MULT_ROWI(l - 1, k - 1, X, *Bptr); /* MAXI = UPDATEMAXI(MAXI, *Bptr); */ FREEMATI(TmpMATI); TmpMATI = *Eptr; *Eptr = ADD_MULT_ROWI(l + i - 1, k + i - 1, X, *Eptr); /* PMAXI = UPDATEMAXI(PMAXI, *Eptr); */ FREEMATI(TmpMATI); FREEMPI(X); for (j = 1; j < l; j++) { X = MULTI((*Lptr)->V[l - 1][j - 1], R); Tmp = (*Lptr)->V[k - 1][j - 1]; (*Lptr)->V[k - 1][j - 1] = SUBI((*Lptr)->V[k - 1][j - 1], X); FREEMPI(Tmp); FREEMPI(X); } X = MULTI(D[l], R); Tmp = (*Lptr)->V[k - 1][l - 1]; (*Lptr)->V[k - 1][l - 1] = SUBI((*Lptr)->V[k - 1][l - 1], X); FREEMPI(Tmp); FREEMPI(X); FREEMPI(R); } for (t = 0; t < m; t++) { if (!EQZEROI((*Bptr)->V[k - 1][t])) { flag = 0; break; } } if (flag) { TmpMATI = *Bptr; *Bptr = DELETE_ROWI(k, *Bptr); FREEMATI(TmpMATI); for (j = k - 1; j < n - i - 1; j++) *Eptr = SWAP_ROWSI1(j + i, j + i + 1, *Eptr); } FREEMPI(Y); return (flag); } void STEP8(USI k, MPMATI **B1ptr, MPMATI **Lptr, MPMATI **Eptr, USI i) { MPI *T; *B1ptr = SWAP_ROWSI1(k - 2, k - 1, *B1ptr); *Eptr = SWAP_ROWSI1(k + i - 2, k + i - 1, *Eptr); T = COPYI((*Lptr)->V[k - 1][ k - 2]); *Lptr = SWAP_ROWSI1(k - 2, k - 1, *Lptr); FREEMPI((*Lptr)->V[ k - 1][k - 2]); (*Lptr)->V[k - 1][k - 2] = T; FREEMPI((*Lptr)->V[k - 2][k - 2]); (*Lptr)->V[k - 2][k - 2] = ZEROI(); return; } void STEP7(USI i, USI k, MPMATI **Lptr, MPI *D[]) { MPI *X1, *X2, *X3, *Y1, *Y2, *Tmp; X1 = MULTI((*Lptr)->V[i - 1][k - 2], (*Lptr)->V[k - 1][k - 2]); Y1 = MULTI((*Lptr)->V[i - 1][k - 1], D[k - 2]); Tmp = Y1; Y1 = ADDI(Y1, X1); FREEMPI(Tmp); FREEMPI(X1); X2 = MULTI((*Lptr)->V[i - 1][k - 2], D[k]); X3 = MINUSI((*Lptr)->V[k - 1][k - 2]); Y2 = MULTI((*Lptr)->V[i - 1][k - 1], X3); FREEMPI(X3); Tmp = Y2; Y2 = ADDI(Y2, X2); FREEMPI(Tmp); FREEMPI(X2); FREEMPI((*Lptr)->V[i - 1][k - 2]); (*Lptr)->V[i - 1][k - 2] = INT(Y1, D[k - 1]); FREEMPI((*Lptr)->V[i - 1][k - 1]); (*Lptr)->V[i - 1][k - 1] = INT(Y2, D[k - 1]); FREEMPI(Y1); FREEMPI(Y2); return; } void BASIS_UPDATE(USI i, USI m, MPMATI **Cptr, MPMATI **Lptr, MPMATI *B1ptr, MPI *D[]) { unsigned int j, k, t; MPI *X, *Tmp, *X1, *X2, *H; for (k = 1; k <= m; k++) { FREEMPI((*Cptr)->V[i - 1][k - 1]); (*Cptr)->V[i - 1][k - 1] = COPYI(B1ptr->V[i - 1][k - 1]); } for (j = 1; j < i; j++) { X = ZEROI(); for (t = 0; t < m; t++) { if (((*Cptr)->V[j - 1][t])->S != 0 && (B1ptr->V[i - 1][t])->S != 0) {H = MULTI((*Cptr)->V[j - 1][t], B1ptr->V[i - 1][t]);Tmp = X; X = ADDI(X, H); FREEMPI(H); FREEMPI(Tmp); } } FREEMPI((*Lptr)->V[i - 1][j - 1]); (*Lptr)->V[i - 1][j - 1] = X; for (t = 0; t < m; t++) { if (((*Cptr)->V[i - 1][t])->S == 0) X1 = ZEROI(); else X1 = MULTI((*Cptr)->V[i - 1][t], D [j]); if (((*Cptr)->V[j - 1][t])->S != 0 && ((*Lptr)->V[i - 1][j - 1])->S != 0) X2 = MULTI((*Cptr)->V[j - 1][t], (*Lptr)->V[i - 1][j - 1]); else X2 = ZEROI(); Tmp = X1; X1 = SUBI(X1, X2); FREEMPI(Tmp); FREEMPI(X2); FREEMPI((*Cptr)->V[i - 1][t]); (*Cptr)->V[i - 1][t] = INT(X1, D[j - 1]); FREEMPI(X1); } } return; } void CSWAP_UPDATE(USI k, USI m, MPI *S, MPMATI **Cptr, MPI *D[]) { unsigned int t; MPI *Tmp1, *Tmp2, *Tmp3, *Tmp4; for (t = 0; t < m; t++) { Tmp1 = MULTI((*Cptr)->V[k - 1][t], D[k - 2]); Tmp2 = MULTI((*Cptr)->V[k - 2][t], S); Tmp3 = ADDI(Tmp1, Tmp2); FREEMPI(Tmp1); FREEMPI(Tmp2); Tmp1 = MULTI((*Cptr)->V[k - 2][t], D[k]); Tmp2 = MULTI((*Cptr)->V[k - 1][t], S); Tmp4 = SUBI(Tmp1, Tmp2); FREEMPI(Tmp1); FREEMPI(Tmp2); FREEMPI((*Cptr)->V[k - 2][t]); FREEMPI((*Cptr)->V[k - 1][t]); (*Cptr)->V[k - 2][t] = INT(Tmp3, D[k - 1]); (*Cptr)->V[k - 1][t] = INT(Tmp4, D[k - 1]); FREEMPI(Tmp3); FREEMPI(Tmp4); } return; } MPMATI *BASIS_REDUCTION0(MPMATI *Bptr, USI m1, USI n1) /* * Input: *Bptr, a matrix of MPI's, whose first row is not zero. * Output: a pointer to an MPMATI whose rows form a reduced basis for * the lattice spanned by the rows of *Bptr. This basis is reduced in the * sense of the paper "Factoring polynomials with rational coefficients" by * A. K. Lenstra, H. W. Lenstra and L. Lovasz, Math. Ann. 261, 515- 534 (1982) * using the modified version in "Solving exponential Diophantine equations * using lattice basis reduction algorithms" by B. M. M. De Weger, J. No. Theory * 26, 325-367 (1987). No change of basis matrix is returned. * De Weger's algorithm has been changed to cater for arbitrary matrices. The * the rows are now in general linearly dependent. * We use the fact that the Gram Schmidt process detects the first row * which is a linear combination of the preceding rows. We employ a modification * of the LLL algorithm outlined by M. Pohst in J. Symbolic Computation (1987)4, * 123-127. We call this the MLLL algorithm. * If we are using this algorithm to find small multipliers for the extended * gcd problem, GCDFLAG is set in EXTGCD() and gcdflag is set below. * m1 / n1 is usually taken to be 3 / 4. */ { unsigned int i, k, l, n, m, t, flag = 0, Flag = 0, gcdflag = 0; unsigned int flagg, beta, K1 = 0, tau = 2, sigma = 0, rho; MPI **D, *X, *Y, *Z, *H, *Tmp, *R, *M1, *N1; MPMATI *C, *L, *B1ptr; unsigned int norig; Z = NULL; m = Bptr->C; n = Bptr->R; norig = n; B1ptr = COPYMATI(Bptr); D = (MPI **)mmalloc((1 + n) * sizeof(MPI *)); D[0] = ONEI(); for (i = 1; i <= n; i++) D[i] = ZEROI(); C = ZEROMNI(n, m); L = ZEROMNI(n, n); found: n = B1ptr->R; i = (K1 == 0) ? 1 : K1; /* K1 = no. of consecutive rows of *B1ptr that don't need updating for the Gram Schmidt process. */ while (i <= B1ptr->R) { BASIS_UPDATE(i, m, &C, &L, B1ptr, D); flag = 1; for (t = 0; t < m; t++) { if (!EQZEROI(C->V[i - 1][t])) { flag = 0; break; } } if (flag) break; X = ZEROI(); for (t = 0; t < m; t++) { H = MULTI(C->V[i - 1][t], C->V[i - 1][t]); Tmp = X; X = ADDI(X, H); FREEMPI(Tmp); FREEMPI(H); } FREEMPI(D[i]); D[i] = INT(X, D[i - 1]); FREEMPI(X); i++; if (MLLLVERBOSE) printf("i = %u\n", i); } beta = (flag) ? i : i - 1; rho = K1 = i - 1; if (MLLLVERBOSE) printf("BASIS0 completed updating the basis\n"); /* Here K1 = no. of LI rows in *B1ptr found by Gram Schmidt process. flag = 0 means all the rho = beta rows of *B1ptr are LI; flag = 1 means that the first rho = beta - 1 rows of *B1ptr are LI, but the beta-th row is a LC of the preceding rows. So beta = number of rows of *B1ptr currently being examined by the LLL algorithm. */ k = tau; M1 = CHANGE(m1); N1 = CHANGE(n1); while (k <= beta) { if (MLLLVERBOSE) printf("beta - k = %u\n", beta -k); l = k - 1; Flag = STEP40(k, l, &L, &B1ptr, D); if (k >= norig && GCDFLAG) { gcdflag = 1; goto FOUND; } if (Flag)/* STEP 9 of POHST. */ { sigma++; if (MLLLVERBOSE) printf("relation vector number %u found\n", sigma); tau = k++; goto found; } X = MULTI(D[k - 2], D[k]); Y = MULTI(D[k - 1], D[k - 1]); Tmp = Y; Y = MULT_I(Y, m1); FREEMPI(Tmp); R = MULTI(L->V[k - 1][k - 2], L->V[k - 1][k - 2]); Z = ADD0I(X, R); Tmp = Z; Z = MULT_I(Z, n1); FREEMPI(Tmp); if (RSV(Y, Z) == 1)/*& STEP 5 of POHST. */ { flagg = 0; if (EQZEROI(D[k]) && EQZEROI(R)) {/* CASE B=0 of STEP 7 of POHST. */ FREEMPI(D[k - 1]); D[k - 1] = ZEROI(); STEP80(k, &B1ptr, &L); if (k - 1 < K1) K1 = k - 1; /* The swap may have changed 2nd last row */ /* of *B1ptr. */ for (t = 0; t < m; t++) { FREEMPI(C->V[k - 2][t]); C->V[k - 2][t] = ZEROI(); } beta--; flagg = 1; if (k > 2) k--; FREEMPI(X); FREEMPI(Y); FREEMPI(R); FREEMPI(Z); continue; } if (flagg == 0) { for (i = k + 1; i <= beta; i++) STEP7(i, k, &L, D); } STEP80(k, &B1ptr, &L); if (k - 2 < K1) K1 = k - 2; /* swap will change last two rows of *B1ptr. */ FREEMPI(R); FREEMPI(Y); Y = MULTI(L->V[k - 1][k - 2], L->V[k - 1][k - 2]); Tmp = Y; Y = ADD0I(Y, X); FREEMPI(X); FREEMPI(Tmp); Tmp = D[k - 1]; D[k - 1] = INT0(Y, D[k - 1]); FREEMPI(Tmp); FREEMPI(Y); if (k > 2) k--; } else { /* STEP 6 of POHST. */ FREEMPI(R); FREEMPI(X); FREEMPI(Y); FOUND: for (l = k - 2; l >= 1; l--) { Flag = STEP40(k, l, &L, &B1ptr, D); if (Flag) { FREEMPI(Z); sigma++; if (MLLLVERBOSE) printf("relation vector number %u found\n", sigma); tau = k++; goto found; /* STEP 9 of POHST. */ } } k++; } if (!gcdflag) FREEMPI(Z); else break; } FREEMPI(M1); FREEMPI(N1); FREEMATI(C); /* printf("L = \n"); PRINTMATI(0,L->R-1,0,L->C-1,L); for (i = 0; i <= norig; i++) { printf("D[%u] = ", i);PRINTI(D[i]);printf(", "); } printf("\n"); */ FREEMATI(L); for (i = 0; i <= norig; i++) FREEMPI(D[i]); ffree((char *)D, (1 + norig) * sizeof(MPI *)); if (MLLLVERBOSE) { printf("number of basis vectors found = %u ;\n", rho); printf("number of relation vectors found = %u .\n", sigma); } return (B1ptr); } unsigned int STEP40(k, l, Lptr, Bptr, D) /* * updates *Lptr, *Bptr. * returns 1 if row k of *Bptr becomes zero, returns zero otherwise. */ unsigned int k, l; MPI *D[]; MPMATI **Lptr, **Bptr; { unsigned int j, flag = 1, t, m; MPI *X, *Y, *R, *Tmp; MPMATI *TmpMATI; m = (*Bptr)->C; Y = MULT_I((*Lptr)->V[k - 1][l - 1], 2); if (RSV(Y, D[l]) == 1) { R = NEAREST_INTI((*Lptr)->V[k - 1][l - 1], D[l]); X = MINUSI(R); TmpMATI = *Bptr; *Bptr = ADD_MULT_ROWI(l - 1, k - 1, X, *Bptr); if (MLLLVERBOSE) { printf("Row %u -> Row %u + ", k,k);PRINTI (X);printf(" x Row %u\n", l); PRINTMATI(0,(*Bptr)->R-1,0,(*Bptr)->C- 1,*Bptr); GetReturn(); } FREEMATI(TmpMATI); /* MAXI = UPDATEMAXI(MAXI, *Bptr); */ FREEMPI(X); for (j = 1; j < l; j++) { X = MULTI((*Lptr)->V[l - 1][j - 1], R); Tmp = (*Lptr)->V[k - 1][j - 1]; (*Lptr)->V[k - 1][j - 1] = SUBI((*Lptr)->V[k - 1][j - 1], X); FREEMPI(Tmp); FREEMPI(X); } X = MULTI(D[l], R); Tmp = (*Lptr)->V[k - 1][l - 1]; (*Lptr)->V[k - 1][l - 1] = SUBI((*Lptr)->V[k - 1][l - 1], X); FREEMPI(Tmp); FREEMPI(X); FREEMPI(R); } for (t = 0; t < m; t++) { if (!EQZEROI((*Bptr)->V[k - 1][t])) { flag = 0; break; } } if (flag) { TmpMATI = *Bptr; *Bptr = DELETE_ROWI(k, *Bptr); FREEMATI(TmpMATI); } FREEMPI(Y); return (flag); } unsigned int STEP400(k, l, Lptr, Bptr, D, gcdflag, norig) /* * updates *Lptr, *Bptr. * returns 1 if row k of *Bptr becomes zero, returns zero otherwise. */ unsigned int k, l, gcdflag, norig; MPI *D[]; MPMATI **Lptr, **Bptr; { unsigned int j, flag = 1, t, m; MPI *X, *Y, *R, *Tmp; MPMATI *TmpMATI; m = (*Bptr)->C; Y = MULT_I((*Lptr)->V[k - 1][l - 1], 2); if (RSV(Y, D[l]) == 1) { R = NEAREST_INTI((*Lptr)->V[k - 1][l - 1], D[l]); X = MINUSI(R); TmpMATI = *Bptr; *Bptr = ADD_MULT_ROWI(l - 1, k - 1, X, *Bptr); FREEMATI(TmpMATI); /* MAXI = UPDATEMAXI(MAXI, *Bptr); */ FREEMPI(X); for (j = 1; j < l; j++) { X = MULTI((*Lptr)->V[l - 1][j - 1], R); Tmp = (*Lptr)->V[k - 1][j - 1]; (*Lptr)->V[k - 1][j - 1] = SUBI((*Lptr)->V[k - 1][j - 1], X); FREEMPI(Tmp); FREEMPI(X); } X = MULTI(D[l], R); Tmp = (*Lptr)->V[k - 1][l - 1]; (*Lptr)->V[k - 1][l - 1] = SUBI((*Lptr)->V[k - 1][l - 1], X); FREEMPI(Tmp); FREEMPI(X); FREEMPI(R); } for (t = 0; t < m; t++) { if (!EQZEROI((*Bptr)->V[k - 1][t])) { flag = 0; break; } } if (flag) { TmpMATI = *Bptr; *Bptr = DELETE_ROWI(k, *Bptr); FREEMATI(TmpMATI); } FREEMPI(Y); return (flag); } void STEP80(k, B1ptr, Lptr) MPMATI **B1ptr, **Lptr; unsigned int k; { MPI *T; *B1ptr = SWAP_ROWSI1(k - 2, k - 1, *B1ptr); if (MLLLVERBOSE) { printf("Swapping Rows %u and %u\n", k-1,k); PRINTMATI(0,(*B1ptr)->R-1,0,(*B1ptr)->C-1,*B1ptr); GetReturn(); } T = COPYI((*Lptr)->V[k - 1][k - 2]); *Lptr = SWAP_ROWSI1(k - 2, k - 1, *Lptr); FREEMPI((*Lptr)->V[k - 1][k - 2]); (*Lptr)->V[k - 1][k - 2] = T; FREEMPI((*Lptr)->V[k - 2][k - 2]); (*Lptr)->V[k - 2][k - 2] = ZEROI(); return; } MPMATI *EXTGCD(MPMATI *Dptr, MPI **Aptr, MPMATI **Q, USI m1, USI n1) /* * Input: an n x 1 vector of MPI's. * Output: *Aptr = gcd of the vector of MPI's. Also we return a small set of * multipliers using the LLL method of Havas and Matthews. * parameters m1 and n1 were put in at the request of George Havas on 11/7/94. * Normally m1/n1 = 3/4. * matrix Q of the LLL extended gcd paper of Havas-Matthews is returned. */ { MPMATI *Temp, *P, *SM, *R; USI i, m, n, *alpha, nz; MPI **Tptr; n = Dptr->R; m = n - 1; alpha = KB_ROWP(Dptr, &R, &nz); Temp = HERMITE1P(Dptr, R, &P, alpha, nz); FREEMATI(R); ffree((char *)alpha, (Dptr->C) * sizeof(USI)); *Aptr = COPYI(Temp->V[0][0]); FREEMATI(Temp); Tptr = P->V[0]; for (i = 0; i < m; i++) P->V[i] = P->V[i+1]; P->V[n - 1] = Tptr; GCDFLAG = 1; *Q = BASIS_REDUCTION0(P, m1, n1); SM = BUILDMATI(1, n); for (i = 0; i < n; i++) SM->V[0][i] = COPYI((*Q)->V[m][i]); FREEMATI(P); GCDFLAG = 0; return (SM); } MPI *LENGTHSQRI(MPMATI *Mptr, USI i) /* * Returns the square of the length of row i of matrix *Mptr. */ { MPI *SUM, *T, *T1; USI j; SUM = ZEROI(); for (j = 0; j < Mptr->C; j++) { T = SUM; T1 = MULTI(Mptr->V[i][j], Mptr->V[i][j]); SUM = ADD0I(SUM, T1); FREEMPI(T); FREEMPI(T1); } return (SUM); } MPI *LENGTHSQCI(MPMATI *Mptr, USI j) /* * Returns the square of the length of column j of matrix *Mptr. */ { MPI *SUM, *T, *T1; USI i; SUM = ZEROI(); for (i = 0; i < Mptr->R; i++) { T = SUM; T1 = MULTI(Mptr->V[i][j], Mptr->V[i][j]); SUM = ADD0I(SUM, T1); FREEMPI(T); FREEMPI(T1); } return (SUM); } MPMATI *BASIS_REDUCTION00(MPMATI *Bptr, USI m1, USI n1, USI norig) /* * Input: *Bptr, a matrix of MPI's, whose first row is not zero. * Output: a pointer to an MPMATI whose rows form a reduced basis for * the lattice spanned by the rows of *Bptr. This basis is reduced in the * sense of the paper "Factoring polynomials with rational coefficients" by * A. K. Lenstra, H. W. Lenstra and L. Lovasz, Math. Ann. 261, 515- 534 (1982) * using the modified version in "Solving exponential Diophantine equations * using lattice basis reduction algorithms" by B. M. M. De Weger, J. No. Theory * 26, 325-367 (1987). No change of basis matrix is returned. * De Weger's algorithm has been changed to cater for arbitrary matrices. The * the rows are now in general linearly dependent. * We use the fact that the Gram Schmidt process detects the first row * which is a linear combination of the preceding rows. We employ a modification * of the LLL algorithm outlined by M. Pohst in J. Symbolic Computation (1987)4, * 123-127. We call this the MLLL algorithm. * If we are using this algorithm to find small multipliers for the extended * gcd problem, GCDFLAG is set in IMPROVEP() and gcdflag is set below. * m1 / n1 is usually taken to be 3 / 4. * For use in IMPROVEP(). */ { unsigned int i, k, n, m, t, flag = 0, Flag = 0, gcdflag = 0; unsigned int flagg, beta, K1 = 0, tau = 2, sigma = 0, rho, norigg; unsigned int K, j, REPEAT, iterate; int l, tt; MPI **D, *X, *Y, *Z, *H, *Tmp, *R, *M1, *N1; MPI *t1, *t2, *t3, *t4; MPMATI *C, *L, *B1ptr, *temp1; Z = NULL; m = Bptr->C; norigg = n = Bptr->R; B1ptr = COPYMATI(Bptr); D = (MPI **)mmalloc((1 + n) * sizeof(MPI *)); D[0] = ONEI(); for (i = 1; i <= n; i++) D[i] = ZEROI(); C = ZEROMNI(n, m); L = ZEROMNI(n, n); found: n = B1ptr->R; i = (K1 == 0) ? 1 : K1; /* K1 = no. of consecutive rows of *B1ptr that don't need updating for the Gram Schmidt process. */ while (i <= B1ptr->R) { BASIS_UPDATE(i, m, &C, &L, B1ptr, D); flag = 1; for (t = 0; t < m; t++) { if (!EQZEROI(C->V[i - 1][t])) { flag = 0; break; } } if (flag) break; X = ZEROI(); for (t = 0; t < m; t++) { H = MULTI(C->V[i - 1][t], C->V[i - 1][t]); Tmp = X; X = ADDI(X, H); FREEMPI(Tmp); FREEMPI(H); } FREEMPI(D[i]); D[i] = INT(X, D[i - 1]); FREEMPI(X); i++; if (MLLLVERBOSE) printf("i = %u\n", i); } /* strcpy(buff, "L.out"); outfile = fopen(buff, "w"); FPRINTMATI(outfile,0,L->R-1,0,L->C-1,L); fclose(outfile); for (t = 0;t < B1ptr->R; t++){ printf("D[%u] = ", t);PRINTI(D[t]);printf("\n"); } GetReturn(); */ beta = (flag) ? i : i - 1; rho = K1 = i - 1; if (MLLLVERBOSE) printf("completed updating the basis\n"); /* Here K1 = no. of LI rows in *B1ptr found by Gram Schmidt process. flag = 0 means all the rho = beta rows of *B1ptr are LI; flag = 1 means that the first rho = beta - 1 rows of *B1ptr are LI, but the beta-th row is a LC of the preceding rows. So beta = number of rows of *B1ptr currently being examined by the LLL algorithm. */ k = tau; M1 = CHANGE(m1); N1 = CHANGE(n1); while (k <= beta) { if (MLLLVERBOSE) printf("beta - k = %u\n", beta -k); if (gcdflag) l = norig; else l = k - 1; Flag = STEP40(k, l, &L, &B1ptr, D); if (k > norig && GCDFLAG) { gcdflag = 1; if (MLLLVERBOSE) printf("improving row %u\n", k); goto FOUND; } if (Flag)/* STEP 9 of POHST. */ { sigma++; if (MLLLVERBOSE) printf("relation vector number %u found\n", sigma); tau = k++; goto found; } X = MULTI(D[k - 2], D[k]); Y = MULTI(D[k - 1], D[k - 1]); Tmp = Y; Y = MULT_I(Y, m1); FREEMPI(Tmp); R = MULTI(L->V[k - 1][k - 2], L->V[k - 1][k - 2]); Z = ADD0I(X, R); Tmp = Z; Z = MULT_I(Z, n1); FREEMPI(Tmp); if (RSV(Y, Z) == 1)/*& STEP 5 of POHST. */ { flagg = 0; if (EQZEROI(D[k]) && EQZEROI(R)) {/* CASE B=0 of STEP 7 of POHST. */ FREEMPI(D[k - 1]); D[k - 1] = ZEROI(); STEP80(k, &B1ptr, &L); if (k - 1 < K1) K1 = k - 1; /* The swap may have changed 2nd last row */ /* of *B1ptr. */ for (t = 0; t < m; t++) { FREEMPI(C->V[k - 2][t]); C->V[k - 2][t] = ZEROI(); } beta--; flagg = 1; if (k > 2) k--; FREEMPI(X); FREEMPI(Y); FREEMPI(R); FREEMPI(Z); continue; } if (flagg == 0) { for (i = k + 1; i <= beta; i++) STEP7(i, k, &L, D); } STEP80(k, &B1ptr, &L); if (k - 2 < K1) K1 = k - 2; /* swap will change last two rows of *B1ptr. */ FREEMPI(R); FREEMPI(Y); Y = MULTI(L->V[k - 1][k - 2], L->V[k - 1][k - 2]); Tmp = Y; Y = ADD0I(Y, X); FREEMPI(X); FREEMPI(Tmp); Tmp = D[k - 1]; D[k - 1] = INT0(Y, D[k - 1]); FREEMPI(Tmp); FREEMPI(Y); if (k > 2) k--; } else { /* STEP 6 of POHST. */ FREEMPI(R); FREEMPI(X); FREEMPI(Y); FOUND: if (gcdflag) { if (MLLLVERBOSE) printf("gcdflag set\n"); K = norig + 1; } else K = k; for (l = K - 2; l >= 1; l--) { Flag = STEP40(k, l, &L, &B1ptr, D); if (Flag) { FREEMPI(Z); sigma++; if (MLLLVERBOSE) printf("relation vector number %u found\n", sigma); tau = k++; goto found; /* STEP 9 of POHST. */ } } k++; } if (!gcdflag) FREEMPI(Z); } FREEMPI(M1); FREEMPI(N1); FREEMATI(C); FREEMATI(L); for (i = 0; i <= norigg; i++) FREEMPI(D[i]); ffree((char *)D, (1 + norigg) * sizeof(MPI *)); /* Now to improve the last n - norig rows of *B1ptr, using Gauss lattice reduction - pointed out by George Havas in Sims' book - 8/11/94: */ /* strcpy(buff, "progress"); outfile = fopen(buff, "a+"); FPRINTMATI(outfile,0,B1ptr->R-1,0,B1ptr->C-1,B1ptr); */ REPEAT = 1; temp1 = BUILDMATI(1, B1ptr->C); for (k = 0; k < B1ptr->C; k++) temp1->V[0][k] = ZEROI(); iterate = 0; QSORTMATI(B1ptr, 0, norig - 1); while (REPEAT) { REPEAT = 0; for (j = 1; j < norig; j++) { t4 = DOTRI(B1ptr, j, j); for (i = 0; i < j; i++) { t1 = DOTRI(B1ptr, j, i); t2 = DOTRI(B1ptr, i, i); t3 = NEAREST_INTI(t1, t2); FREEMPI(t1); FREEMPI(t2); if (t3->S == 0) { FREEMPI(t3); continue; } t1 = t3; t3 = MINUSI(t3); FREEMPI(t1); for (k = 0; k < B1ptr->C; k++) { FREEMPI(temp1->V[0][k]); X = MULTI(t3, B1ptr->V[i][k]); temp1->V[0][k] = ADDI(X, B1ptr->V[j][k]); FREEMPI(X); } t1 = DOTRI(temp1, 0, 0); tt = RSV(t1, t4); tt = RSV(t1, t4); if (tt == -1) { ADD_MULT_ROWI0(i, j, t3, B1ptr); if (MLLLVERBOSE) { printf("Gauss- improving nullspace row %u: length squared was = ", j + 1); PRINTI(t4); printf("\n"); printf("lengthsquared now = "); PRINTI(t1); printf("\n"); } /* fprintf(outfile, "Gauss- improving nullspace row %u: length squared was = ", j + 1); FPRINTI(outfile, t4); fprintf(outfile, "\n"); fprintf (outfile, "lengthsquared now = "); FPRINTI(outfile, t1); fprintf(outfile, "\n"); fflush(outfile); */ FREEMPI(t4); t4 = t1; REPEAT = 1; } else FREEMPI(t1); FREEMPI(t3); } FREEMPI(t4); } iterate++; QSORTMATI(B1ptr, 0, norig - 1); } REPEAT = 1; iterate = 0; while (REPEAT) { REPEAT = 0; for (j = norig; j < B1ptr->R; j++) { t4 = DOTRI(B1ptr, j, j); for (i = 0; i < norig; i++) { t1 = DOTRI(B1ptr, j, i); t2 = DOTRI(B1ptr, i, i); t3 = NEAREST_INTI(t1, t2); FREEMPI(t1); FREEMPI(t2); if (t3->S == 0) { FREEMPI(t3); continue; } t1 = t3; t3 = MINUSI(t3); FREEMPI(t1); for (k = 0; k < B1ptr->C; k++) { FREEMPI(temp1->V[0][k]); X = MULTI(t3, B1ptr->V[i][k]); temp1->V[0][k] = ADDI(X, B1ptr->V[j][k]); FREEMPI(X); } t1 = DOTRI(temp1, 0, 0); tt = RSV(t1, t4); tt = RSV(t1, t4); if (tt == -1) { ADD_MULT_ROWI0(i, j, t3, B1ptr); if (MLLLVERBOSE) { printf("Gauss- improving row %u: length squared was = ", j + 1); PRINTI(t4); printf("\n"); printf("lengthsquared now = "); PRINTI(t1); printf("\n"); } /* fprintf(outfile, "Gauss- improving row %u: length squared was = ", j + 1); FPRINTI(outfile, t4); fprintf(outfile, "\n"); fprintf (outfile, "lengthsquared now = "); FPRINTI(outfile, t1); fprintf(outfile, "\n"); fflush(outfile); */ FREEMPI(t4); t4 = t1; REPEAT = 1; } else FREEMPI(t1); FREEMPI(t3); } FREEMPI(t4); } iterate++; } FREEMATI(temp1); /* fclose(outfile); */ return (B1ptr); } MPMATI *BASIS_REDUCTION000(MPMATI *Bptr, USI m1, USI n1, MPI *N) /* * Input: *Bptr, a matrix of MPI's, whose first row is not zero. * Output: a pointer to an MPMATI whose rows form a reduced basis for * the lattice spanned by the rows of *Bptr. This basis is reduced in the * sense of the paper "Factoring polynomials with rational coefficients" by * A. K. Lenstra, H. W. Lenstra and L. Lovasz, Math. Ann. 261, 515- 534 (1982) * using the modified version in "Solving exponential Diophantine equations * using lattice basis reduction algorithms" by B. M. M. De Weger, J. No. Theory * 26, 325-367 (1987). No change of basis matrix is returned. * De Weger's algorithm has been changed to cater for arbitrary matrices. The * the rows are now in general linearly dependent. * We use the fact that the Gram Schmidt process detects the first row * which is a linear combination of the preceding rows. We employ a modification * of the LLL algorithm outlined by M. Pohst in J. Symbolic Computation (1987)4, * 123-127. We call this the MLLL algorithm. * If we are using this algorithm to find small multipliers for the extended * gcd problem, GCDFLAG is set in EXTGCD() and gcdflag is set below. * m1 / n1 is usually taken to be 3 / 4. */ { unsigned int i, k, l, n, m, t, flag = 0, Flag = 0, gcdflag = 0; unsigned int flagg, beta, K1 = 0, tau = 2, sigma = 0, rho; MPI **D, *X, *Y, *Z, *H, *Tmp, *R, *M1, *N1; MPMATI *C, *L, *B1ptr; unsigned int norig; Z = NULL; m = Bptr->C; n = Bptr->R; norig = n; B1ptr = COPYMATI(Bptr); D = (MPI **)mmalloc((1 + n) * sizeof(MPI *)); D[0] = ONEI(); for (i = 1; i <= n; i++) D[i] = ZEROI(); C = ZEROMNI(n, m); L = ZEROMNI(n, n); found: n = B1ptr->R; i = (K1 == 0) ? 1 : K1; /* K1 = no. of consecutive rows of *B1ptr that don't need updating for the Gram Schmidt process. */ while (i <= B1ptr->R) { BASIS_UPDATE(i, m, &C, &L, B1ptr, D); flag = 1; for (t = 0; t < m; t++) { if (!EQZEROI(C->V[i - 1][t])) { flag = 0; break; } } if (flag) break; X = ZEROI(); for (t = 0; t < m; t++) { H = MULTI(C->V[i - 1][t], C->V[i - 1][t]); Tmp = X; X = ADDI(X, H); FREEMPI(Tmp); FREEMPI(H); } FREEMPI(D[i]); D[i] = INT(X, D[i - 1]); FREEMPI(X); i++; if (MLLLVERBOSE) printf("i = %u\n", i); } beta = (flag) ? i : i - 1; rho = K1 = i - 1; if (MLLLVERBOSE) printf("completed updating the basis\n"); /* Here K1 = no. of LI rows in *B1ptr found by Gram Schmidt process. flag = 0 means all the rho = beta rows of *B1ptr are LI; flag = 1 means that the first rho = beta - 1 rows of *B1ptr are LI, but the beta-th row is a LC of the preceding rows. So beta = number of rows of *B1ptr currently being examined by the LLL algorithm. */ k = tau; M1 = CHANGE(m1); N1 = CHANGE(n1); while (k <= beta) { if (MLLLVERBOSE) printf("beta - k = %u\n", beta -k); l = k - 1; Flag = STEP4000(k, l, &L, &B1ptr, D, N); if (k >= norig && GCDFLAG) { gcdflag = 1; goto FOUND; } if (Flag)/* STEP 9 of POHST. */ { sigma++; if (MLLLVERBOSE) printf("relation vector number %u found\n", sigma); tau = k++; goto found; } X = MULTI(D[k - 2], D[k]); Y = MULTI(D[k - 1], D[k - 1]); Tmp = Y; Y = MULT_I(Y, m1); FREEMPI(Tmp); R = MULTI(L->V[k - 1][k - 2], L->V[k - 1][k - 2]); Z = ADD0I(X, R); Tmp = Z; Z = MULT_I(Z, n1); FREEMPI(Tmp); if (RSV(Y, Z) == 1)/*& STEP 5 of POHST. */ { flagg = 0; if (EQZEROI(D[k]) && EQZEROI(R)) {/* CASE B=0 of STEP 7 of POHST. */ FREEMPI(D[k - 1]); D[k - 1] = ZEROI(); STEP8000(k, &B1ptr, &L, N); if (k - 1 < K1) K1 = k - 1; /* The swap may have changed 2nd last row */ /* of *B1ptr. */ for (t = 0; t < m; t++) { FREEMPI(C->V[k - 2][t]); C->V[k - 2][t] = ZEROI(); } beta--; flagg = 1; if (k > 2) k--; FREEMPI(X); FREEMPI(Y); FREEMPI(R); FREEMPI(Z); continue; } if (flagg == 0) { for (i = k + 1; i <= beta; i++) STEP7(i, k, &L, D); } STEP8000(k, &B1ptr, &L, N); if (k - 2 < K1) K1 = k - 2; /* swap will change last two rows of *B1ptr. */ FREEMPI(R); FREEMPI(Y); Y = MULTI(L->V[k - 1][k - 2], L->V[k - 1][k - 2]); Tmp = Y; Y = ADD0I(Y, X); FREEMPI(X); FREEMPI(Tmp); Tmp = D[k - 1]; D[k - 1] = INT0(Y, D[k - 1]); FREEMPI(Tmp); FREEMPI(Y); if (k > 2) k--; } else { /* STEP 6 of POHST. */ FREEMPI(R); FREEMPI(X); FREEMPI(Y); FOUND: for (l = k - 2; l >= 1; l--) { Flag = STEP4000(k, l, &L, &B1ptr, D, N); if (Flag) { FREEMPI(Z); sigma++; if (MLLLVERBOSE) printf("relation vector number %u found\n", sigma); tau = k++; goto found; /* STEP 9 of POHST. */ } } k++; } if (!gcdflag) FREEMPI(Z); else break; } FREEMPI(M1); FREEMPI(N1); FREEMATI(C); FREEMATI(L); for (i = 0; i <= norig; i++) FREEMPI(D[i]); ffree((char *)D, (1 + norig) * sizeof(MPI *)); if (MLLLVERBOSE) { printf("number of basis vectors found = %u ;\n", rho); printf("number of relation vectors found = %u .\n", sigma); } return (B1ptr); } unsigned int STEP4000(USI k, USI l, MPMATI **Lptr, MPMATI **Bptr, MPI *D[], MPI *N) /* * updates *Lptr, *Bptr. * returns 1 if row k of *Bptr becomes zero, returns zero otherwise. */ { unsigned int i, j, flag = 1, t, m, n; MPI *X, *Y, *R, *Tmp, *Temp; MPMATI *TmpMATI; m = (*Bptr)->C; Y = MULT_I((*Lptr)->V[k - 1][l - 1], 2); if (RSV(Y, D[l]) == 1) { R = NEAREST_INTI((*Lptr)->V[k - 1][l - 1], D[l]); X = MINUSI(R); TmpMATI = *Bptr; *Bptr = ADD_MULT_ROWI(l - 1, k - 1, X, *Bptr); FREEMATI(TmpMATI); n = (*Bptr)->R; for (i = 0; i < n; i++) { for (j = n; j < m; j++) { Temp = POWERI(N, m - j); (*Bptr)->V[i][j] = INTI((*Bptr)->V[i][j], Temp); FREEMPI(Temp); } } if (HERMITEVERBOSE) { printf("Row %u -> Row %u + ", k,k);PRINTI(X);printf(" x Row %u\n", l); PRINTMATI(0,(*Bptr)->R-1,0,(*Bptr)->C-1,*Bptr); GetReturn(); } for (i = 0; i < n; i++) { for (j = n; j < m; j++) { Temp = POWERI(N, m - j); (*Bptr)->V[i][j] = MULTI((*Bptr)->V[i][j], Temp); FREEMPI(Temp); } } /* MAXI = UPDATEMAXI(MAXI, *Bptr); */ FREEMPI(X); for (j = 1; j < l; j++) { X = MULTI((*Lptr)->V[l - 1][j - 1], R); Tmp = (*Lptr)->V[k - 1][j - 1]; (*Lptr)->V[k - 1][j - 1] = SUBI((*Lptr)->V[k - 1][j - 1], X); FREEMPI(Tmp); FREEMPI(X); } X = MULTI(D[l], R); Tmp = (*Lptr)->V[k - 1][l - 1]; (*Lptr)->V[k - 1][l - 1] = SUBI((*Lptr)->V[k - 1][l - 1], X); FREEMPI(Tmp); FREEMPI(X); FREEMPI(R); } for (t = 0; t < m; t++) { if (!EQZEROI((*Bptr)->V[k - 1][t])) { flag = 0; break; } } if (flag) { TmpMATI = *Bptr; *Bptr = DELETE_ROWI(k, *Bptr); FREEMATI(TmpMATI); } FREEMPI(Y); return (flag); } void STEP8000(USI k, MPMATI **B1ptr, MPMATI **Lptr, MPI *N) { MPI *T, *Temp; USI i, j, n, m; *B1ptr = SWAP_ROWSI1(k - 2, k - 1, *B1ptr); T = COPYI((*Lptr)->V[k - 1][k - 2]); *Lptr = SWAP_ROWSI1(k - 2, k - 1, *Lptr); FREEMPI((*Lptr)->V[k - 1][k - 2]); (*Lptr)->V[k - 1][k - 2] = T; FREEMPI((*Lptr)->V[k - 2][k - 2]); (*Lptr)->V[k - 2][k - 2] = ZEROI(); n = (*B1ptr)->R; m = (*B1ptr)->C; for (i = 0; i < n; i++) { for (j = n; j < m; j++) { Temp = POWERI(N, m - j); (*B1ptr)->V[i][j] = INTI((*B1ptr)->V[i][j], Temp); FREEMPI(Temp); } } if (HERMITEVERBOSE) { printf("Swapping Rows %u and %u\n", k-1,k); PRINTMATI(0,(*B1ptr)->R-1,0,(*B1ptr)->C-1,*B1ptr); GetReturn(); } for (i = 0; i < n; i++) { for (j = n; j < m; j++) { Temp = POWERI(N, m - j); (*B1ptr)->V[i][j] = MULTI((*B1ptr)->V[i][j], Temp); FREEMPI(Temp); } } return; }
Message: 2023 - Contents - Hide Contents Date: Tue, 20 Nov 2001 20:44:05 Subject: Re: LLL basis for linear temperaments From: genewardsmith@xxxx.xxx --- In tuning-math@y..., Graham Breed <graham@m...> wrote:> The confusion was because he was coming up with unique results for > inconsistent temperaments. Firstly, he's said that he takes the nearest > prime approximations, which should clear up the confusion although I'm not > sure our results agree. Secondly, this is exactly the problem he mentions in > that paragraph!The uniqueness was in connection with generators in ets. In my notation, 12+34 uniquely defines a generator, since we add 12 and 34, get 46, and take the nearest primes *for 46*.
Message: 2024 - Contents - Hide Contents Date: Tue, 20 Nov 2001 20:52:41 Subject: Re: LLL definitions From: genewardsmith@xxxx.xxx --- In tuning-math@y..., Graham Breed <graham@m...> wrote:>> Well that kind of slams the door on the "classification" I was hoping >> to acheive, though I suppose a classification by optimal generators >> would be good enough (and would eliminate the need to use LLL, but I >> feel we should still have some lower bound on >> allowable "orthogonality").> The mapping by period and generator is unique for any linear temperament, > which should be relevant here. Some caveats:Minkowski reduction should do what Paul wants; the caveats would concern symmetries which shouldn't happen.> There are always two choices of generator within the period. I used to > always take the smaller, but sometimes it became the larger after > optimization. So you can make the arbitrary choice that the first nonzero > term be negative.The lattice reduction approach doesn't always give us generator/period; it is optimizing steps to get to good intevals, and that isn't always generator/period.
First Previous Next
2000 2050 2100 2150 2200 2250 2300 2350 2400 2450 2500 2550 2600 2650 2700 2750 2800 2850 2900 2950
2000 - 2025 -