robpike / ivy

ivy, an APL-like calculator
Other
1.32k stars 103 forks source link

ivy: implement 'floor' and 'ceil' for complex numbers #148

Closed fzipp closed 1 year ago

fzipp commented 1 year ago

The implementation uses the definition of 'floor' and 'ceil' for complex numbers set by

McDonnell, E. E. "Complex Floor, APL Congress 73." (1973).
https://www.jsoftware.com/papers/eem/complexfloor.htm
robpike commented 1 year ago

Let's check in this code, after editing the comment (or not; your choice). I then have some refactoring I could do that would help smooth it out. At the very least, floor should be factored out of this.

fzipp commented 1 year ago

Let's check in this code, after editing the comment (or not; your choice). I then have some refactoring I could do that would help smooth it out. At the very least, floor should be factored out of this.

Thanks

fzipp commented 1 year ago

The next step would be to make mod work for complex numbers, as defined in the paper:

op z mod w = z - w * floor z / w + w == 0

This would be a straight forward translation:

func (c Complex) mod(ctx Context, w Complex) Complex {
    zeroFlag := complexZero
    if isZero(w) {
        zeroFlag = complexOne
    }
    return c.sub(ctx, w.mul(ctx, (c.div(ctx, w.add(ctx, zeroFlag))).floor(ctx)))
}

It involves several costly Eval calls as well.

robpike commented 1 year ago

I leave this to you.

Dyalog, which claims to implement the same algorithm, says

⌊2.9j3.9 2J4

While ivy says

floor 2.9j3.9 3j3

One of them is wrong.

robpike commented 1 year ago

It's possibly a rounding error because 0.9 is not representable exactly in a float64, while ivy uses rationals. However,

floor 2.9j2.9 3j2

  ⌊2.9j2.9

2J3

which is peculiar.

fzipp commented 1 year ago

I noticed this as well. Dyalog uses > instead of >= in the following condition, deviating from the paper:

    if isTrue("floor", ctx.EvalBinary(x, ">=", y)) {

Our current implementation follows the paper. I have not found a publicly stated reason for why Dyalog does it differently.

robpike commented 1 year ago

So when it's on the diagonal, it could go either way. That makes sense. Peculiar though...

fzipp commented 1 year ago

J (see playground https://www.jdoodle.com/ia/KFQ) and GNU APL both return 3j2, in accordance with the paper.

robpike commented 1 year ago

Thanks for investigating. I'll add this test case.

fzipp commented 1 year ago

It's indirectly covered by these test cases for ceil, but more test cases are a good idea, especially one explicitly for floor.

ceil -0.5j0.5
    -1j1

ceil -0.5j-0.5
    -1