csound / manual

Csound Reference Manual (English)
Other
46 stars 29 forks source link

What's correct behavior of round()? Example results no longer correct! #219

Closed tgrey1 closed 4 years ago

tgrey1 commented 4 years ago

Looking for feedback about what is correct here... the manual page for round() states:

Returns the integer value nearest to x ; if the fractional part of x is exactly 0.5, the direction of rounding is undefined.

But this no longer appears to be the case. I tested manually with a few more values, and it seems like it's consistently rounding down when the decimal portion equals .5

Running the manual example yields different results than what is listed as output on the manual page. (trimmed output to the one relevant line) Manual: instr 1: inumber = 9.000 idiv = 2.000 ifl = 5.000 Actual: instr 1: inumber = 9.000 idiv = 2.000 iro = 4.000

This was tested with Csound version 6.14 (double samples) Apr 21 2020 [commit: 393bfeb2f358c2baa3e302b47337ab479b298b64]

I'll also update the output for the minor change of ifl vs. iro for the variable name, but that's secondary to making sure the opcode works correctly and is described accurately.

Thanks!

jpffitch commented 4 years ago

Does it round the same on all platforms? I assume that the comment was to allow for different compilers a f libraries.

⁣Get TypeApp for Android ​

On 7 Aug 2020, 09:38, at 09:38, tgrey1 notifications@github.com wrote:

Looking for feedback about what is correct here... the manual page for round() states:

Returns the integer value nearest to x ; if the fractional part of x is exactly 0.5, the direction of rounding is undefined.

But this no longer appears to be the case. I tested manually with a few more values, and it seems like it's consistently rounding down when the decimal portion equals .5

Running the manual example yields different results than what is listed as output on the manual page. (trimmed output to the one relevant line) Manual: instr 1: inumber = 9.000 idiv = 2.000 ifl = 5.000 Actual: instr 1: inumber = 9.000 idiv = 2.000 iro = 4.000

This was tested with Csound version 6.14 (double samples) Apr 21 2020 [commit: 393bfeb2f358c2baa3e302b47337ab479b298b64]

I'll also update the output for the minor change of ifl vs. iro for the variable name, but that's secondary to making sure the opcode works correctly and is described accurately.

Thanks!

-- You are receiving this because you are subscribed to this thread. Reply to this email directly or view it on GitHub: https://github.com/csound/manual/issues/219

vlazzarini commented 4 years ago

Since round uses the C library, we may need to see what the documentation says. On MacOs

     The round() functions return the integral value nearest to x rounding half-way
     cases away from zero, regardless of the current rounding direction.

     The lround() and llround() functions return the integral value nearest to x
     (rounding half-way cases away from zero, regardless of the current rounding direc-
     tion) in the return formats specified.  If the rounded value is outside the range
     of the return type, the numeric result is unspecified and the "invalid" floating-
     point exception is raised. A range error may occur if the magnitude of x is too
     large.

which appears to be the opposite of always rounding down.

tgrey1 commented 4 years ago

Sorry I should have specified, my testing was run on Mac OS 10.13.6.

This also gets weirder... I was wrong saying it consistently rounded down. It looks like maybe even numbers (including 0) round down, whereas odd numbers round up?

ival1 = round(.5)
ival2 = round(1.5)
ival3 = round(2.5)
ival4 = round(3.5)
print .5, ival1
print 1.5, ival2
print 2.5, ival3
print 3.5, ival4

yields:

instr 1:  .5 = 0.500  ival1 = 0.000
instr 1:  1.5 = 1.500  ival2 = 2.000
instr 1:  2.5 = 2.500  ival3 = 2.000
instr 1:  3.5 = 3.500  ival4 = 4.000
vlazzarini commented 4 years ago

I forgot that we actually use lrint in linux/MacOs (doubles), so

    The rint() functions return the integral value nearest to x (according to the pre-
     vailing rounding mode) in floating-point format.

     The lrint() and llrint() functions return the integral value nearest to x (accord-
     ing to the prevailing rounding mode) in the return formats specified.  If the
     rounded value is outside the range of the return type, the numeric result is
     unspecified and the "invalid" floating-point exception is raised. A range error
     may occur if the magnitude of x is too large.

Now on other systems/float-precision combinations, other code is used. I think what's in the manual is the best description.

tgrey1 commented 4 years ago

Ok, so differing results in this case is expected? I find that odd, but if that's the case perhaps the statement if the fractional part of x is exactly 0.5, the direction of rounding is undefined. could be made a little more specific, to something like: if the fractional part of x is exactly 0.5, the direction of rounding is dependent on the host operating system.

Ideally, maybe we could even collect descriptions of the expected behaviors on specific different common operating systems? I feel like that would be helpful and avoid future confusion for others.

vlazzarini commented 4 years ago

That could be, but in the case of lrint() that is not even explicitly defined, so I am not sure, but undefined appears to be more appropriate.

tgrey1 commented 4 years ago

Ok, I agree... not specifying each expected behavior makes sense, since it may be unclear in some cases and might even change over time unbeknownst to us.

But I still feel it's important to explicitly mention in the manual that results may differ between operating systems with this opcode, especially since it's observable (for some) with the example. If it's ok, I'd still like to add that suggestion of: if the fractional part of x is exactly 0.5, the direction of rounding is undefined and dependent on the host operating system.

Tho I'm also open to any similar wording if anyone has better suggestions. I feel like stating it as undefined doesn't make clear to users that the output may change.

As an aside, is there a similar function to round that can be depended on to behave the same on all systems?

vlazzarini commented 4 years ago

Sounds OK, "and may depend" could be used, since say with lrint, it's the same function used in both MacOs and Linux and I suspect it has the same (undefined) behaviour in both.

tgrey1 commented 4 years ago

Out of curiosity, I'll test further with @tjingboem tomorrow and figure out for sure. Between the two of us we can test it on OSX/Lin/Win.

For the manual, I'll use: if the fractional part of x is exactly 0.5, the direction of rounding is undefined and may depend on the host operating system.

vlazzarini commented 4 years ago

sounds good. I think we are using the fastest functions in each system that have it (lrint etc).

tgrey1 commented 4 years ago

We've tested a little more, and the results were same across platforms: OSX 10.13.6, Windows 10, and Linux Mint 20.

Each time, if the integral portion is even it rounds down... and if the integral portion is odd it rounds up. It seems like it's using what this site calls "Banker's Rounding", aka round to even: https://www.mathsisfun.com/numbers/rounding-methods.html

This doesn't seem intuitive to me as a default behavior... but maybe there's a reason it's being done this way? Obviously it wasn't always this way, since at one point the example yielded the results that are currently represented in the manual.

Tomorrow I'll try grabbing the oldest csound build I can find and testing with that.

Test conducted:

<CsoundSynthesizer>
<CsOptions>
; Select audio/midi flags here according to platform
; Audio out   Audio in
-odac           -iadc    ;;;RT audio I/O
; For Non-realtime ouput leave only the line below:
; -o subinstr_named.wav -W ;;; for file output any platform
</CsOptions>
<CsInstruments>
; Initialize the global variables.
sr = 44100
ksmps = 10
nchnls = 2
0dbfs  = 1
instr 1
iCount init 0
while iCount < 10 do
    iVal = iCount+.5
    iRound = round(iVal)
    print iVal, iRound
    iCount +=1
od
endin
</CsInstruments>
<CsScore>
i1 0 .1
e
</CsScore>
</CsoundSynthesizer>

Output on all tested systems so far:

instr 1:  iVal = 0.500  iRound = 0.000
instr 1:  iVal = 1.500  iRound = 2.000
instr 1:  iVal = 2.500  iRound = 2.000
instr 1:  iVal = 3.500  iRound = 4.000
instr 1:  iVal = 4.500  iRound = 4.000
instr 1:  iVal = 5.500  iRound = 6.000
instr 1:  iVal = 6.500  iRound = 6.000
instr 1:  iVal = 7.500  iRound = 8.000
instr 1:  iVal = 8.500  iRound = 8.000
instr 1:  iVal = 9.500  iRound = 10.000
tgrey1 commented 4 years ago

Ok, I have confirmed this is a regression/behavior change. I installed csound 5.19 from: http://dream.cs.bath.ac.uk/Csound-archive/OldReleases/csound5/csound5.19/

I had to alter the test slightly since the "while" opcode didn't exist yet.

New test:

<CsoundSynthesizer>
<CsOptions>
; Select audio/midi flags here according to platform
; Audio out   Audio in
-odac           -iadc    ;;;RT audio I/O
; For Non-realtime ouput leave only the line below:
; -o subinstr_named.wav -W ;;; for file output any platform
</CsOptions>
<CsInstruments>
; Initialize the global variables.
sr = 44100
ksmps = 10
nchnls = 2
0dbfs  = 1

instr 1
iCount init 0

loop:
    iVal = iCount+.5
    iRound = round(iVal)
    print iVal, iRound
    iCount = iCount + 1
if (iCount < 10) igoto loop

endin

</CsInstruments>
<CsScore>
i1 0 .1
e
</CsScore>
</CsoundSynthesizer>

Output from 5.19:

instr 1:  iVal = 0.500  iRound = 1.000
instr 1:  iVal = 1.500  iRound = 2.000
instr 1:  iVal = 2.500  iRound = 3.000
instr 1:  iVal = 3.500  iRound = 4.000
instr 1:  iVal = 4.500  iRound = 5.000
instr 1:  iVal = 5.500  iRound = 6.000
instr 1:  iVal = 6.500  iRound = 7.000
instr 1:  iVal = 7.500  iRound = 8.000
instr 1:  iVal = 8.500  iRound = 9.000
instr 1:  iVal = 9.500  iRound = 10.000

Output from 6.14:

instr 1:  iVal = 0.500  iRound = 0.000
instr 1:  iVal = 1.500  iRound = 2.000
instr 1:  iVal = 2.500  iRound = 2.000
instr 1:  iVal = 3.500  iRound = 4.000
instr 1:  iVal = 4.500  iRound = 4.000
instr 1:  iVal = 5.500  iRound = 6.000
instr 1:  iVal = 6.500  iRound = 6.000
instr 1:  iVal = 7.500  iRound = 8.000
instr 1:  iVal = 8.500  iRound = 8.000
instr 1:  iVal = 9.500  iRound = 10.000

Should I go ahead and open this as a bug against csound core?

vlazzarini commented 4 years ago

I don't think it is a regression, it may be that a different C function was used. That does not mean a regression. It is still undefined for that value. Unless you can demonstrate two things

(1) that the rounding behaviour is that one and the same for all platforms; and

(2) that the computation load is less or equal to what is now;

it is not an issue in the code. The documentation continues to be correct as well.

joachimheintz commented 4 years ago

but i think the result from 6.14 is correct, whereas the 5.19 output is strange, not?

On 09/08/2020 12:12, tgrey1 wrote:

Ok, I have confirmed this is a regression/behavior change. I installed csound 5.19 from: http://dream.cs.bath.ac.uk/Csound-archive/OldReleases/csound5/csound5.19/

I had to alter the test slightly since the "while" opcode didn't exist yet.

New test:

; Select audio/midi flags here according to platform ; Audio out Audio in -odac -iadc;;;RT audio I/O ; For Non-realtime ouput leave only the line below: ; -o subinstr_named.wav -W ;;; for file output any platform ; Initialize the global variables. sr = 44100 ksmps = 10 nchnls = 2 0dbfs = 1 instr 1 iCount init 0 loop: iVal = iCount+.5 iRound = round(iVal) print iVal, iRound iCount = iCount + 1 if (iCount < 10) igoto loop endin i1 0 .1 e

Output from 5.19:

|instr 1: iVal = 0.500 iRound = 1.000 instr 1: iVal = 1.500 iRound = 2.000 instr 1: iVal = 2.500 iRound = 3.000 instr 1: iVal = 3.500 iRound = 4.000 instr 1: iVal = 4.500 iRound = 5.000 instr 1: iVal = 5.500 iRound = 6.000 instr 1: iVal = 6.500 iRound = 7.000 instr 1: iVal = 7.500 iRound = 8.000 instr 1: iVal = 8.500 iRound = 9.000 instr 1: iVal = 9.500 iRound = 10.000 |

Output from 6.14:

|instr 1: iVal = 0.500 iRound = 0.000 instr 1: iVal = 1.500 iRound = 2.000 instr 1: iVal = 2.500 iRound = 2.000 instr 1: iVal = 3.500 iRound = 4.000 instr 1: iVal = 4.500 iRound = 4.000 instr 1: iVal = 5.500 iRound = 6.000 instr 1: iVal = 6.500 iRound = 6.000 instr 1: iVal = 7.500 iRound = 8.000 instr 1: iVal = 8.500 iRound = 8.000 instr 1: iVal = 9.500 iRound = 10.000 |

Should I go ahead and open this as a bug against csound core?

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub https://github.com/csound/manual/issues/219#issuecomment-671033344, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAQYHKREUU4PHQSG6U2MBPDR7ZZCXANCNFSM4PXM3YEA.

vlazzarini commented 4 years ago

The fact is that rounding of mid-way values is undefined.

tgrey1 commented 4 years ago

I can not demonstrate those two things, so I guess I'll stick with the previously agreed to changes to the manual .

Is it reasonable to put in a request for a new rounding function that can have predictable defined results across OS and csound versions in the future? This difference could have a very significant impact on instruments or entire compositions across versions if they depend on round(), which seems like a pretty common and critical function to be able to rely on.

Maybe, if possible, this improved version could even take an additional optional i-rate argument to specify a rounding method, so a user could explicitly request whichever is appropriate to their use case (half round up, half round down, bankers, random, etc)

vlazzarini commented 4 years ago

There is no way to specify rounding behaviour using the C library functions. However, you can implement this by hand; an UDO seems to suffice. Something like

opcode Round,i,i ival xin xout (frac(ival) < 0.5 ? int(ival) : int(ival) + 1) endop

kunstmusik commented 4 years ago

I find the current behavior a bit odd myself. I think the Round UDO doesn't work for negative numbers?

For comparison sake, we have this in CS6:

https://github.com/csound/csound/blob/develop/include/sysdep.h#L353-L412

and this in Csound 5.17:

https://github.com/csound/csound/blob/csound5_17_3/H/sysdep.h#L314-L373

jpffitch commented 4 years ago

My thoughts entirely

⁣Sent from TypeApp ​

On Aug 9, 2020, 12:33, at 12:33, vlazzarini notifications@github.com wrote:

I don't think it is a regression, it may be that a different C function was used. That does not mean a regression. It is still undefined for that value. Unless you can demonstrate two things

(1) that the rounding behaviour is that one and the same for all platforms; and

(2) that the computation load is less or equal to what is now;

it is not an issue in the code. The documentation continues to be correct as well.

-- You are receiving this because you commented. Reply to this email directly or view it on GitHub: https://github.com/csound/manual/issues/219#issuecomment-671040636

vlazzarini commented 4 years ago

I don't see why it should be confusing. It's just undefined as elsewhere, if you are using the C library. The UDO I showed it's just an example. You can use floor() or ceil() instead of int(). What I meant is that we use the fastest code (lrint) in 64 bit unices, I don't think we should change that.

opcode Round,i,i
  ival xin
  xout (ival < 0 ? (frac(ival) < 0.5 ? int(ival) - 1 : int(ival)) : (frac(ival) < 0.5 ? int(ival) : int(ival) + 1))
endop
kunstmusik commented 4 years ago

Took a look and the behavior we have now is apparently the default for IEEE-754 rounding operations:

https://en.wikipedia.org/wiki/Rounding#Round_half_to_even

I suppose we could expose calling of setting rounding mode via opcode if desired but leave what we have as default.

vlazzarini commented 4 years ago

To be fair, I think this is not needed, since it can be easily solved with an UDO.

tgrey1 commented 4 years ago

@vlazzarini - I'd already written something similar, but yours is more elegant and works better as one-liner! I used it as a basis for a new version with an optional switch for rounding halves up, down, or to continue using the undefined round() behavior.

opcode Round,i,io
    iVal, iRoundType xin
    if iRoundType>0 then
        iVal = (frac(iVal) < 0.5 ? int(iVal) : int(iVal) + 1)
    elseif iRoundType<0 then
        iVal = (frac(iVal) <= 0.5 ? int(iVal) : int(iVal) + 1)
    else
        iVal = round(iVal)
    endif
    xout iVal
endop

@kunstmusik

I suppose we could expose calling of setting rounding mode via opcode if desired but leave what we have as default.

That's essentially what I had been hoping for, either with the original or with a new opcode. The problem I see is that if a user calls round() with data input that will become audible (particularly turning streams of data into midi notes), the composition could be noticeably different between systems and version compilations if half numbers are regularly encountered in the input. It seems to me like a pretty important functionality to have built in... but if I'm in a minority on that opinion, maybe it's not really worth it and I'm imagining this usefulness (a real possibility).

To be fair, I think this is not needed, since it can be easily solved with an UDO.

Most of the new opcodes being added could have been UDOs, in fact many of them were before their inclusion!

But regardless, that discussion is outside of the scope of this particular issue, which was originally just meant to help make sure we update the manual to be as accurate as possible... but the results of testing and the responses went down a deeper rabbit hole than expected. Perhaps I'll open a new issue and refer back here for context so that this one can return it's focus on the manual page clarification.

Please leave this open until I update the entry, it's already assigned to me. I plan on taking a day or two soon and knocking out a bunch of these small updates. Thanks.