Dear all, Thank you for your very useful suggestions.
I: (i) changed the package argument to Join(Field,RadicalCategory), which simplified the code; (ii) included Johannes' signValue function, which is an elegant solution to AlgebraicNumber not having a sign function; and simplified the householder function by using the function signature Vector(R)->Matrix(R). I did not include the complex case or pivoting in the solution at http://rosettacode.org/wiki/QR_decomposition#Axiom, simply because they were outside the problem statement. As http://axiom-wiki.newsynthesis.org/ was down, a draft solution to the complex case (without pivots) is given below. Kindly, Mark. )abbrev package TESTP TestPackage VR ==> Vector(R) MR ==> Matrix(R) NNI ==> NonNegativeInteger TestPackage(F:Join(Field,RadicalCategory,TranscendentalFunctionCategory),R:ComplexCategory(F)): with transposeH: MR -> MR unitVector: NNI -> VR norm: VR -> F conjugate: VR -> VR "/": (R,MR) -> MR "/": (R,VR) -> VR "/": (VR,R) -> VR "^": (VR,NNI) -> VR solveUpperTriangular: (MR,VR) -> VR householder: VR -> MR qr: MR -> Record(q:MR,r:MR) lsqr: (MR,VR) -> VR polyfit: (VR,VR,NNI) -> VR == add import VectorFunctions2(R,F) transposeH(m:MR):MR == transpose(map(conjugate, m)) conjugate(v:VR):VR == map((vi:R):R +-> conjugate(vi)$R, v)$VR unitVector(dim) == out := new(dim,0@R)$VR out(1) := 1@R out a:R / m:MR == map((mij:R):R +-> a/mij, m)$MR a:R / v:VR == map((vi:R):R +-> a/vi, v)$VR v:VR / a:R == map((vi:R):R +-> vi/a, v)$VR v:VR ^ n:NonNegativeInteger == map((vi:R):R +-> vi^n, v)$VR norm(v:Vector(R)):F == reduce("+",map(norm,v)) solveUpperTriangular(r,b) == n := ncols r x := new(n,0@R)$VR for k in n..1 by -1 repeat index := min(n,k+1) x(k) := (b(k)-reduce("+",subMatrix(r,k,k,index,n)*x.(index..n)))/r(k,k) x householder(a) == m := #a alpha : R := -a(1)*inv(sqrt norm a(1))*sqrt(norm(a)) e : VR := unitVector(m) u : VR := a + alpha*e v : VR := u/complex(sqrt norm(u),0) w : R := dot(conjugate a,v)/dot(conjugate v,a) scalarMatrix(m,1) - (1+w)*transpose(outerProduct(v,conjugate v)) qr(a) == m := nrows a n := ncols a qm := scalarMatrix(m,1) rm : MR := copy a for i in 1..(if m=n then n-1 else n) repeat x := column(subMatrix(rm,i,m,i,n),1) h := scalarMatrix(m,1) setsubMatrix!(h,i,i,householder x) qm := qm*h rm := h*rm [qm,rm] )if false )co qr-complex.spad -- example m := matrix [[12, -51, 4+%i], [6, 167, -68], [-4, 24, -41]] out := qr m -- tests - OK! transpose(map(conjugate,out.q))*out.q out.q*transpose(map(conjugate,out.q)) out.q * out.r - m -- R code -- a <- matrix(complex(real=c(12, -51, 4, 6, 167, -68, -4, 24, -41),im=c(0,0,1,0,0,0,0,0,0)),byrow=T,nrow=3) -- d <- qr(a) -- qr.Q(d) -- qr.R(d) -- NOT the same: pivoting? )end if -----Original Message----- From: [email protected] [mailto:[email protected]] On Behalf Of Prof. Dr. Johannes Grabmeier Sent: Tuesday, 13 March 2012 8:58 AM To: [email protected] Subject: Re: [fricas-devel] Domain that combines Field, RadicalCategory and OrderedRing? Mark, a few years ago I did some similar implementations (householder operations for use to computenumeric Groebner bases) and solved the problems your raised in the style how Ralf suggested it. I did not yet put my code into the system, but we could share our codes and insights, if you want. Here a typical code: signValue(r: R): R == -- todo: make appropriate categories for R R has (sign: R -> R) => sign(r)$R zero? r => r -- hack, c.f. complex numbers R has (sqrt: R -> R) => if sqrt(r*r) = r then 1 else -1 1 Am 13.03.2012 um 08:25 schrieb Ralf Hemmecke: >> Well, sign for AlgebraicNumber is less defined than for complex >> numbers: > > Mark, > > In fact, the only reason why you need AlgebraicNumber is probably because of > the use of the length function to get the length of the vector and that the > length function is only implemented if the R in Vector(R) has RadicalCategory. > > I don't, however, see any reason to use Float anywhere in your code. You > basically use it only for implementing the sign function. > Better would perhaps be to require that R already exports a function > > sign; % -> Integer > > that would also remove the headache made by Waldek above. Or better said, it > moves it from your package to the respective field you are working in. If the > field has no idea what a sign is, then it's anyway hard to introduce one in > your package. > > BTW, is there any reason that you require just R: Ring and not immediately R: > Field in the argument type of TestPackage? Or even Join(Field,OrderedRing)? > > Ralf > > -- > You received this message because you are subscribed to the Google Groups > "FriCAS - computer algebra system" group. > To post to this group, send email to [email protected]. > To unsubscribe from this group, send email to > [email protected]. > For more options, visit this group at > http://groups.google.com/group/fricas-devel?hl=en. > Mit freundlichen Grüßen Johannes Grabmeier Prof. Dr. Johannes Grabmeier Köckstraße 1, D-94469 Deggendorf Tel. +49-(0)-991-2979584, Tel. +49-(0)-171-5503789 Tel. +49-(0)-991-3615-100 (d), Fax: +49-(0)-1803-5518-17745 -- You received this message because you are subscribed to the Google Groups "FriCAS - computer algebra system" group. To post to this group, send email to [email protected]. To unsubscribe from this group, send email to [email protected]. For more options, visit this group at http://groups.google.com/group/fricas-devel?hl=en. -- You received this message because you are subscribed to the Google Groups "FriCAS - computer algebra system" group. To post to this group, send email to [email protected]. To unsubscribe from this group, send email to [email protected]. For more options, visit this group at http://groups.google.com/group/fricas-devel?hl=en.
