Jay,
Before your email, I had read portions of McDonnell's Complex
Floor article. Now I've read it in its entirety. I do not see
anything in the article that contradicts my view the residue function
should be idempotent and capable of reproducing the complete residue
system for a given modulus.
Gaussian integers ARE complex numbers. If the APL residue
function has to work well on all complex numbers, it has to work well
on Gaussian integers which is not happening right now in Gnu APL.
Respectfully,
Fred
On Fri, 2017-06-23 at 08:34 +0100, Jay Foad wrote:
> I urge you to read Eugene McDonnell's Complex Floor, which also
> discusses Residue. I believe the design he comes up with in this
> paper was adopted more or less verbatim in APL. Also bear in mind
> that Floor and Residue in APL have to work well on all complex
> numbers, not just the Gaussian integers.
> Jay.
> On 23 June 2017 at 01:42, Frederick Pitts <[email protected]>
> wrote:
> >
> >
> > Hello Jürgen,
> >
> > Some observations:
> >
> > 1) When performing a residue calculation on positive integers, a
> > straight-forward integer division with remainder calculation
> > suffices. For example, 5 ∣ 13 is computed with 13 / 5 = 2 r 3 and
> > so 5 ∣ 13 = 3 where 3 is in the complete residue system
> > { 0, 1, 2, 3, 4 }. When performing the calculation on negative
> > integers, one has to take advantage of the fact that the
> > integer division quotient and remainder are not unique in order to
> > compute a residue that is in the complete residue system.
> > For 5 ∣ ¯13, ¯13 / 5 = ¯2 r ¯3 where ¯3 is not in the
> > CRS. However, ¯13 / 5 = ¯3 r 2 where 3 is in the CSR. The same
> > concept applies Gaussian integers.
> >
> > 2) I suspect the decision to have the APL2 floor function round
> > toward negative infinity, instead of toward zero,
> > was made based on the desire to save cpu cycles and memory in the
> > residue function code.
> >
> > 3) I read at least one math literature article discussing Gaussian
> > integer Euclidean division algorithms, that recommended
> > rounding down to the nearest real and imaginary part toward
> > negative infinity. Unfortunately I cannot find
> > the article right now. I will continue to look for it. None of
> > the articles discussed using a complex integer floor function.
> >
> > 4) The reason MOD_TEST.apl shows total disagreement MODJ and the
> > builtin residue function is that the complex floor function code
> > change in SVN 965 relocated the CRS's on complex plane. Attached
> > are CRS0-CRS1-6J-6-SVN964.out
> > CRS0-CRS1-6J-6-SVN965.out. The first file contains a CRS map for
> > modulus ¯6J¯6 produced with the residue function
> > followed by a map for the same modulus produced with MODJ using SVN
> > 964. The second file contains the same maps
> > using SVN 965. Observe that for SVN 964 the residue function CRS
> > is in the bottom half of the complex plane, but for SVN 965 it is
> > in the top half. The CRS for the MODJ function is in the bottom
> > half in both SVN cases.
> > 5)The complex floor code change did not help with the issue that
> > the builtin residue function is not idempotent for all possible
> > arguments and consequently generates too many residues. See
> > attached CRSOTST0-SVN965.out. For a grid
> > of Gaussian integers with real and imaginary parts ranging from ¯15
> > to 15, using every value with every other value as modulus and
> > second argument, there were 40 case where the order of CSR exceeded
> > the modulus norm. I think that
> > was the failure count with the previous SVN.
> >
> > Sincerely, I think the complex floor and ceiling functions
> > should not be used by other functions even if IBM and ISO
> > imply they are in their documentations. I'm not seeing them used
> > in the Gaussian integer literature. Again, please correct me if
> > I'm wrong.
> >
> > Regards,
> >
> > Fred
> >
> > On Thu, 2017-06-22 at 18:08 +0200, Juergen Sauermann wrote:
> > > Hi again,
> > >
> > >
> > >
> > > sorry a small typo below. Lines 19/20 should read:
> > >
> > >
> > >
> > > (¯6J¯5 - 0J¯11) ÷ ¯6J¯6
> > >
> > > 0J¯1
> > >
> > >
> > >
> > > /// Jürgen
> > >
> > >
> > >
> > > On 06/22/2017 05:44 PM, Juergen
> > > Sauermann wrote:
> > >
> > >
> > >
> > > >
> > > > Hi Fred at al.,
> > > >
> > > >
> > > >
> > > > I have made another attempt to fix the residue
> > > > function, SVN
> > > > 965.
> > > >
> > > >
> > > >
> > > > For complex m∣b It now rounds down the real() and
> > > > imag()
> > > > parts of the quotient q←b÷m and returns b-q.
> > > >
> > > > Instead of always rounding towards 0 or -infinity, the
> > > > rounding
> > > > direction is now (compared to the previous
> > > >
> > > > attempt) determined by the quadrant in which the
> > > > modulus m
> > > > lies.
> > > >
> > > >
> > > >
> > > > There are still differences to the results displayed by
> > > > MOD_test.apl, but I
> > > > suppose they are
> > > >
> > > > caused by that program. For example, the first line of
> > > > MOD_test.apl, says:
> > > >
> > > >
> > > >
> > > > LA
> > > > RA MODJ |
> > > >
> > > > ¯6J¯6 ¯6J¯5 0J¯11 0J1
> > > >
> > > >
> > > >
> > > > We have:
> > > >
> > > >
> > > >
> > > > (¯6J¯5 - 0J1) ÷ ¯6J¯6
> > > >
> > > > 1
> > > >
> > > > (0J¯11 - 0J1) ÷ ¯6J¯6
> > > >
> > > > 1J1
> > > >
> > > >
> > > >
> > > > so both 0J¯11 and 0J1 are valid remainders
> > > > modulo ¯6J¯6. However, the
> > > >
> > > > magnitude of 0J¯11
> > > > (= 11) is larger than the magnitude of the divisor
> > > > ¯6J¯6 (= around 8.4).
> > > >
> > > > I suppose most people expect the remainder of a
> > > > division to be
> > > > in some sense
> > > >
> > > > smaller than the divisor.
> > > >
> > > >
> > > >
> > > > Regarding Jay's idempotency requirement we now have:
> > > >
> > > >
> > > >
> > > >
> > > > f←{6J6|⍵}
> > > >
> > > > f ¯3 ¯2 ¯3 ¯1 0 1 2 3
> > > >
> > > > 3J6 4J6 3J6 5J6 0 ¯5J6 ¯4J6 ¯3J6
> > > >
> > > > f f ¯3 ¯2 ¯3 ¯1 0 1 2 3
> > > >
> > > > 3J6 4J6 3J6 5J6 0 ¯5J6 ¯4J6 ¯3J6
> > > >
> > > >
> > > >
> > > > f←{5J3|⍵}
> > > >
> > > > f ¯3 ¯2 ¯3 ¯1 0 1 2 3
> > > >
> > > > 2J3 3J3 2J3 4J3 0 ¯2J5 ¯1J5 0J5
> > > >
> > > > f f ¯3 ¯2 ¯3 ¯1 0 1 2 3
> > > >
> > > > 2J3 3J3 2J3 4J3 0 ¯2J5 ¯1J5 0J5
> > > >
> > > >
> > > >
> > > > so at least the first modulus seems to work as well.
> > > > The
> > > > result is still different
> > > >
> > > > from APL2 as reported by Jay, but I can't tell why:
> > > >
> > > >
> > > >
> > > > IBM APL2:
> > > >
> > > >
> > > >
> > > > 5J3 ∣
> > > > 14J5 1J4 ¯4J1
> > > >
> > > > ¯4J1 ¯4J1 ¯4J1
> > > >
> > > >
> > > >
> > > > GNU APL:
> > > >
> > > >
> > > >
> > > > 5J3 ∣ 14J5 1J4 ¯4J1
> > > >
> > > > 1J4 1J4 1J4
> > > >
> > > >
> > > >
> > > > But both 1J4 and ¯4J1 are
> > > > valid remainders. Interestingly Jay's idempotency
> > > > requirement
> > > > seems to
> > > >
> > > > be fulfilled by both the GNU APL and by IBM APL2, so
> > > > that that
> > > > requirement alone does not suffice
> > > >
> > > > to tell which result is correct.
> > > >
> > > >
> > > >
> > > > On the other hand this matter seems to be like
> > > > discussing if
> > > > the square root of 4 is 2 or -2 with
> > > >
> > > > both answers being correct.
> > > >
> > > >
> > > >
> > > > Best Regards,
> > > >
> > > > Jürgen Sauermann
> > > >
> > > >
> > > >
> > > >
> > > >
> > > >
> > > >
> > > >
> > > > On 06/21/2017 10:25 PM, Frederick
> > > > Pitts wrote:
> > > >
> > > >
> > > >
> > > > >
> > > > > Jürgen,
> > > > >
> > > > >
> > > > >
> > > > > The proposed change to DIVJ does not work because
> > > > > 'q1' is
> > > > > a complex number, so the '×' in '× q1' is the
> > > > > complex
> > > > > complement function, not the sign function. I tried
> > > > > the
> > > > > proposed change and every test fails.
> > > > >
> > > > >
> > > > >
> > > > > I will try to hack DIVJ to use a floor towards zero
> > > > > instead of towards minus infinity for the real and
> > > > > imaginary
> > > > > parts of the quotient and see what happens.
> > > > >
> > > > >
> > > > >
> > > > > Correct me if I am wrong, but my mind set is that
> > > > > the APL
> > > > > residue function has to satisfy the following
> > > > > invariants:
> > > > > 1) The order of the complete residue system (residue
> > > > > count)
> > > > > for a given modulo 'n' has to equal the norm of
> > > > > 'n'.
> > > > > 2) And as Jay Foad so succinctly expressed it, the
> > > > > residue
> > > > > function has to be idempotent with respect to its
> > > > > right
> > > > > argument,
> > > > > e.g., ( n | m ) = n | n | m .
> > > > > regardless of the implementation of the residue
> > > > > function.
> > > > >
> > > > >
> > > > >
> > > > > Later,
> > > > >
> > > > >
> > > > >
> > > > > Fred
> > > > >
> > > > >
> > > > >
> > > > >
> > > > >
> > > > >
> > > > >
> > > > >
> > > > >
> > > > >
> > > > >
> > > > >
> > > > > On Wed, 2017-06-21 at 19:46 +0200, Juergen Sauermann
> > > > > wrote:
> > > > >
> > > > > > Hi Fred,
> > > > > >
> > > > > >
> > > > > >
> > > > > > I have a question about the MOD_test.apl that
> > > > > > you
> > > > > > kindly provided.
> > > > > >
> > > > > >
> > > > > >
> > > > > > In function DIVJ on line 57 ff it says:
> > > > > >
> > > > > >
> > > > > >
> > > > > > z ← q , a -
> > > > > > b × q ← CMPLX ⌊ ( 9 11 ) ○ a ÷ b
> > > > > >
> > > > > >
> > > > > >
> > > > > > so the quotient is rounded down towards minus
> > > > > > infinity.
> > > > > >
> > > > > >
> > > > > >
> > > > > > I wonder if that should be something like
> > > > > >
> > > > > >
> > > > > >
> > > > > > z ← q , (×
> > > > > > q1) × a - b × q ← CMPLX ⌊ ∣ 9 11 ○ q1 ← a ÷
> > > > > > b
> > > > > >
> > > > > >
> > > > > >
> > > > > > so that the quotient is rounded towards 0?
> > > > > > Interestingly IBM
> > > > > > and ISO
> > > > > >
> > > > > > give different definitions for the residue in
> > > > > > terms of APL:
> > > > > >
> > > > > >
> > > > > >
> > > > > > IBM (language reference, page 227):
> > > > > >
> > > > > > Z←L∣R
> > > > > >
> > > > > > Z is R-L×⌊ R÷L+L=0
> > > > > >
> > > > > >
> > > > > >
> > > > > > ISO (chapter 7.2.9 Residue):
> > > > > >
> > > > > > R←Q∣P
> > > > > >
> > > > > > R←P-(×P)×|Q×⌊|P÷Q
> > > > > >
> > > > > > and return R if (×R)=×Q, or R+Q
> > > > > >
> > > > > > otherwise.
> > > > > >
> > > > > >
> > > > > >
> > > > > > That suggest that IBM rounds the quotient down
> > > > > > towards minus
> > > > > > infinity while ISO rounds
> > > > > >
> > > > > > towards 0.
> > > > > >
> > > > > >
> > > > > >
> > > > > > My naive view on remainder is that the nearest
> > > > > > integer
> > > > > > quotient shall be smaller in
> > > > > >
> > > > > > magnitude and not smaller in value. Regarding
> > > > > > your proposal
> > > > > > (which is different from
> > > > > >
> > > > > > both IBM and ISO) my concern is that may lead to
> > > > > > different
> > > > > > results for modulo N and
> > > > > >
> > > > > > modulo N×1J0
> > > > > >
> > > > > >
> > > > > >
> > > > > > Best Regards,
> > > > > >
> > > > > > Jürgen Sauermann
> > > > > >
> > > > > >
> > > > > >
> > > > > >
> > > > > >
> > > > > > On 06/21/2017 03:08 AM, Frederick
> > > > > > Pitts wrote:
> > > > > >
> > > > > >
> > > > > >
> > > > > > > Jürgen,
> > > > > > >
> > > > > > >
> > > > > > >
> > > > > > > This message is
> > > > > > > being resent because last minute changes
> > > > > > > I made to
> > > > > > > CRS0.apl and CRS1.apl do not output the
> > > > > > > data I intended.
> > > > > > > This message has corrected versions of
> > > > > > > those files
> > > > > > > attached. Please discard the old
> > > > > > > CRS0.apl and CRS1.apl
> > > > > > > files. The first line of output is the
> > > > > > > modulo basis, the second line is the
> > > > > > > calculated complete residue system values
> > > > > > > and the third
> > > > > > > line is the number of residues in the
> > > > > > > CRS on the previous line.
> > > > > > >
> > > > > > >
> > > > > > >
> > > > > > > CRSOTST0.apl and
> > > > > > > CRSOTST1.apl are unchanged.
> > > > > > >
> > > > > > >
> > > > > > >
> > > > > > > Also please find
> > > > > > > attached MOD_TEST.apl which compares the
> > > > > > > residues
> > > > > > > calculated by MODJ and the builtin
> > > > > > > residue function and
> > > > > > > reports discrepancies. The first column
> > > > > > > of output is
> > > > > > > the modulo basis, the second column the
> > > > > > > right argument
> > > > > > > to the functions, the third column the
> > > > > > > MODJ result and
> > > > > > > the fourth column is the builtin residue
> > > > > > > function
> > > > > > > result.
> > > > > > >
> > > > > > >
> > > > > > >
> > > > > > > Regards
> > > > > > >
> > > > > > >
> > > > > > >
> > > > > > > Fred
> > > > > > >
> > > > > > >
> > > > > > >
> > > > > > > Hello Jürgen,
> > > > > > > SVN 964 moved us in the right
> > > > > > > direction but not completely out of the
> > > > > > > woods. SVN 964 still exhibits errors. For
> > > > > > > instance
> > > > > > > 2J6 | 5J5
> > > > > > > ¯1J7
> > > > > > > 2J6 | ¯1J7
> > > > > > > ¯3J1
> > > > > > > 2J6 | ¯3J1
> > > > > > > ¯3J1
> > > > > > > I found this and previous residue
> > > > > > > function errors using the attached APL
> > > > > > > code files. The files with base name ending
> > > > > > > in '0' use the builtin residue
> > > > > > > function. Those with base name ending in '1'
> > > > > > > use a residue function implemented
> > > > > > > in APL. The files with base name beginning
> > > > > > > with 'CRSOTST' test if the order of
> > > > > > > the complete residue system (CRS) equals the
> > > > > > > norm of the modulo basis. That
> > > > > > > test fails for several modulo bases, 2J6
> > > > > > > being one of them, using the builtin
> > > > > > > residue function. No errors are detected with
> > > > > > > the APL implementation. The other files
> > > > > > > can be used to plot the CRS for a given
> > > > > > > modulo basis where 'a' and 'b' in
> > > > > > > 'a + b * i' are limited to +15 to -15
> > > > > > > range. A full screen terminal window is
> > > > > > > needed to see the plot.
> > > > > > > My APL implementation of the residue
> > > > > > > function is very close to what you
> > > > > > > described in your previous email. Maybe
> > > > > > > comparing the two implementations will
> > > > > > > give insight into why the builtin residue
> > > > > > > function fails for some modulo bases.
> > > > > > > I make no assertion that my
> > > > > > > implementation is correct in all
> > > > > > > aspects.
> > > > > > > Regards,
> > > > > > > Fred
> > > > > > > On Tue, 2017-06-20 at 14:14 +0200, Juergen
> > > > > > > Sauermann
> > > > > > > wrote:
> > > > > > >
> > > > > > > > Hi Frederick,
> > > > > > > >
> > > > > > > >
> > > > > > > >
> > > > > > > > the algorithm for A ∣ B used in GNU APL
> > > > > > > > is this:
> > > > > > > >
> > > > > > > >
> > > > > > > >
> > > > > > > > -
> > > > > > > > compute the quotient Q←B÷A,
> > > > > > > >
> > > > > > > > - "round down" Q to the next
> > > > > > > > (complex) integer
> > > > > > > > Q1,
> > > > > > > >
> > > > > > > > - return B - Q1×A
> > > > > > > >
> > > > > > > >
> > > > > > > >
> > > > > > > > Now the problem seems to be what is
> > > > > > > > meant by "round
> > > > > > > > down". There are two candidates:
> > > > > > > >
> > > > > > > >
> > > > > > > >
> > > > > > > > Q1 ←
> > > > > > > > ⌊
> > > > > > > > Q i.e.
> > > > > > > > use APL floor to round down Q
> > > > > > > >
> > > > > > > > Q1 ← Complex( floor(Q.real(),
> > > > > > > > floor(Q.imag())
> > > > > > > > ) i,e, use C/C++ floor() to round
> > > > > > > > down Q.
> > > > > > > >
> > > > > > > >
> > > > > > > >
> > > > > > > > In your 5J3 ∣ 14J5 example, the
> > > > > > > > quotient is 2.5J¯0.5,
> > > > > > > > which gives different results for the
> > > > > > > > APL floor ⌊
> > > > > > > > and the C/C++ floor().
> > > > > > > >
> > > > > > > >
> > > > > > > >
> > > > > > > > The APL floor ⌊2.5J¯0.5
> > > > > > > > is 3J¯1 (a somewhat dubious
> > > > > > > > invention in the ISO standard on page
> > > > > > > > 19, which I used
> > > > > > > > up to
> > > > > > > >
> > > > > > > > including SVN 963), while the C/C++
> > > > > > > > floor() is 2J¯1.
> > > > > > > > The difference between the APL floor
> > > > > > > > and the C/C++ floor
> > > > > > > > is 1.0 which,
> > > > > > > >
> > > > > > > > multiplied by the divisor, explains the
> > > > > > > > differences that
> > > > > > > > we see.
> > > > > > > >
> > > > > > > >
> > > > > > > >
> > > > > > > > As of SVN 964 I have changed the
> > > > > > > > residue
> > > > > > > > function (∣) to use the C/C++ floor
> > > > > > > > instead of
> > > > > > > > the APL floor. The APL floor and
> > > > > > > >
> > > > > > > > Ceiling functions (⌊ and ⌈) are still
> > > > > > > > using the apparently broken definition
> > > > > > > > in the ISO
> > > > > > > > standard.
> > > > > > > >
> > > > > > > >
> > > > > > > >
> > > > > > > > I hope this works better for you. At
> > > > > > > > least I am getting
> > > > > > > > this in SVN 964:
> > > > > > > >
> > > > > > > >
> > > > > > > >
> > > > > > > >
> > > > > > > > 5J3 | 14J5
> > > > > > > >
> > > > > > > > 1J4
> > > > > > > >
> > > > > > > > 5J3 | 1J4
> > > > > > > >
> > > > > > > > 1J4
> > > > > > > >
> > > > > > > >
> > > > > > > >
> > > > > > > > whereas SVN 963 was giving:
> > > > > > > >
> > > > > > > >
> > > > > > > >
> > > > > > > >
> > > > > > > > 5J3 | 14J5
> > > > > > > >
> > > > > > > > ¯4J1
> > > > > > > >
> > > > > > > > 5J3 | 1J4
> > > > > > > >
> > > > > > > > ¯4J1
> > > > > > > >
> > > > > > > >
> > > > > > > >
> > > > > > > >
> > > > > > > >
> > > > > > > > Best Regards,
> > > > > > > >
> > > > > > > > /// Jürgen
> > > > > > > >
> > > > > > > >
> > > > > > > >
> > > > > > > >
> > > > > > > >
> > > > > > > >
> > > > > > > >
> > > > > > > > On 06/19/2017 07:03 PM,
> > > > > > > > Frederick Pitts wrote:
> > > > > > > >
> > > > > > > >
> > > > > > > >
> > > > > > > > > Jürgen,
> > > > > > > > >
> > > > > > > > > With gnu apl (svn 961 on Fedora 25, Intel(R)
> > > > > > > > > Core(TM) i7-6700
> > > > > > > > > CPU), the residue function (∣) yields the following:
> > > > > > > > >
> > > > > > > > > 5J3 ∣ 14J5
> > > > > > > > > 1J4
> > > > > > > > > 5J3 | 1J4
> > > > > > > > > ¯4J1
> > > > > > > > > 5J3 | ¯4J1
> > > > > > > > > ¯4J1
> > > > > > > > > The above result means that two elements in the
> > > > > > > > > complete residue system
> > > > > > > > > (CSR) for mod 5J3 are equal, i.e. 1J4 = ¯4J1 mod 5J3,
> > > > > > > > > which is not
> > > > > > > > > allowed. None of the elements of a CSR can be equal
> > > > > > > > > modulo the CSR's
> > > > > > > > > basis.
> > > > > > > > >
> > > > > > > > > Regards,
> > > > > > > > >
> > > > > > > > > Fred
> > > > > > > > >
> > > > > > > > >
> > > > > > > > >
> > > > > > > >
> > > > > > > >
> > > > > > > >
> > > > > > > >
> > > > > > >
> > > > > > >
> > > > > >
> > > > > >
> > > > > >
> > > > > >
> > > > >
> > > > >
> > > >
> > > >
> > > >
> > > >
> > >
> > >
> > >
> > >
> > >
