[Jprogramming] Round to nearest integer: harder than it seems

Henry Rich henryhrich at gmail.com
Fri Jan 15 01:53:04 UTC 2016


Comparison tolerance.

     z=.0.5+2^50

     x: z

1125899906842625

     x:!.0 z

2251799813685249r2


Henry Rich

On 1/14/2016 8:47 PM, Don Guinn wrote:
> I tried using x: on a floating point number to see what happens. The
> dictionary said that it produces a rational for floating point numbers.
> Well, it is a rational, but it rounded to the nearest integer. I expected
> it to give a rational with a denominator of 2. This was done on a 64 bit J.
>
>     z=.0.5+2^50
>     z
> 1.1259e15
>     0j16":z
> 1125899906842624.5000000000000000
>     x:z
> 1125899906842625
>     x:^:_1 x:z
> 1125899906842625
>     3!:0 x:^:_1 x:z
> 4
>     3!:0 x:z
> 128
>     z=.1r2+2x^50
>     z
> 2251799813685249r2
>     x:^:_1 z
> 1.1259e15
>     0j16":x:^:_1 z
> 1125899906842624.5000000000000000
>
> On Thu, Jan 14, 2016 at 6:30 PM, Marshall Lochbaum <mwlochbaum at gmail.com>
> wrote:
>
>> There are definitely cases where you absolutely need your numbers to be
>> integers, and I was wondering how to address that. Thankfully, J
>> functions that take integers will take floats which are close enough,
>> but if you want to pass numbers to C, or write them to a binary file,
>> they need to be the type you have set for them.
>>
>> Interestingly, my function (fl) can replace the final (<.) in roundNums
>> to guarantee the input is an integer if TO is (integer) 1. It's
>> compatible with any rounding direction.
>>
>> Marshall
>>
>> On Thu, Jan 14, 2016 at 07:33:11PM -0500, Devon McCormick wrote:
>>> Personally, I avoid paying attention to internal data types whenever
>>> possible as this leads us down an implementation-dependent rabbit hole
>> with
>>> little upside.
>>>
>>> My own rounding utility "roundNums" - found here:
>>> http://code.jsoftware.com/wiki/User:Devon_McCormick/Utilities - allows
>>> options for different kinds of rounding - up, down, and banker's - for
>> the
>>> cases where we're on an exact halfway point.  Banker's rounding - also
>>> referenced here http://code.jsoftware.com/wiki/Addons/general/misc -
>>> attempts to avoid systematic bias by rounding either up or down.  This
>>> would be important, e.g., in the case of a series resulting from the
>>> average of two series of whole numbers where rounding either up or down
>>> would introduce directional bias.
>>>
>>> Also, I like my way of specifying the rounding precision with a minimum
>>> increment, i.e. "1" for whole numbers or "0.05" for nickels, as this is
>>> intuitive and flexible.
>>>
>>> On Thu, Jan 14, 2016 at 7:01 PM, Marshall Lochbaum <mwlochbaum at gmail.com
>>>
>>> wrote:
>>>
>>>> Actually, you're off by one on the low end--a quirk of two's complement
>>>> notation is that the minimum possible value is actually one less than
>>>> minus the maximum. But the bigger problem is that doing things with
>> rank
>>>> zero is very slow. I was working with audio data where there are
>>>> typically almost a hundred thousand values per second (and, admittedly,
>>>> where any more than 24 bits in the output is overkill), so this
>> wouldn't
>>>> be fast enough.
>>>>
>>>> Marshall
>>>>
>>>> On Thu, Jan 14, 2016 at 09:39:20PM +0000, 'Pascal Jasmin' via
>> Programming
>>>> wrote:
>>>>> anything wrong with this?
>>>>>
>>>>> (0.5 <.@:+ ])`(9223372036854775807 * *)@.(9223372036854775807 <
>> |)("0) _
>>>> 234234.3 __
>>>>> works as an adverb where you provide your rounding function.
>>>>>
>>>>>
>>>>> roundA =: `(9223372036854775807 * *)(@.(9223372036854775807 < |))("0)
>>>>>
>>>>> (0.5 <.@:+ ]) roundA _ 234234.3 __
>>>>>> . roundA _ 234234.3 __
>>>>> ----- Original Message -----
>>>>> From: Marshall Lochbaum <mwlochbaum at gmail.com>
>>>>> To: programming at jsoftware.com
>>>>> Sent: Thursday, January 14, 2016 2:56 PM
>>>>> Subject: [Jprogramming] Round to nearest integer: harder than it
>> seems
>>>>> Here's something I spent far too long on, and consequently thought
>> was
>>>>> worth sharing. I can turn it into an essay on the J wiki if people
>> want
>>>>> that.
>>>>>
>>>>> Recently I ran into the problem of rounding a J floating point
>> number to
>>>>> an integer, and forcing the result to have integer type. This seems
>> like
>>>>> a simple task: using a standard rounding function, we have
>>>>>     0.5 <.@:+ 2.3 5.1 7.6 3.9
>>>>> 2 5 8 4
>>>>> But with numbers that are too large, the result still contains
>>>>> floating-point numbers, and has type 8 (floating point) rather than 4
>>>>> (integer).
>>>>>     0.5 <.@:+ _1e50 2e30 _
>>>>> _1e50 2e30 _
>>>>>     3!:0 ]0.5 <.@:+ _1e50 2e30 _
>>>>> 8
>>>>> When applied to a float, (<.) applies the C floor function, which
>> yields
>>>>> another float, and than casts the results to integers if all of them
>> are
>>>>> exactly representable as integers. They're not here, so they are
>> left as
>>>>> floating-point numbers.
>>>>>
>>>>> To give a more accurate problem statement, I want the 64-bit signed
>>>>> integer which is closest to the function input. Thus numbers above
>> the
>>>>> maximum representable integer should round to that integer, and
>> likewise
>>>>> for numbers below the minimum representable integer. We define these
>> two
>>>>> bounds now.
>>>>>     MAX =: ->: MIN =: _2 <.@^ 63
>>>>> Note that since MAX is one less than 2^63, trying to take (2<.@^63)
>>>>> would give us a float, and subtracting one would still leave us with
>> a
>>>>> floating point number, which is not actually equal to MAX since
>> (>:MAX)
>>>>> is representable as a float, while nearby integers are not. MIN on
>> the
>>>>> other hand is safely computed as an exponent. Note the negative base,
>>>>> which works because 63 is odd.
>>>>>
>>>>> With these bounds, our problem should be easy: clamp to the integer
>>>>> range, then use (<.).
>>>>>     ([: <. MIN>.MAX<.]) 0.5 + __ _1e30 _1e10 _100.3
>>>>> _9223372036854775808 _9223372036854775808 _10000000000 _100
>>>>> So far, so good...
>>>>>     ([: <. MIN>.MAX<.]) 0.5 + 50.4 2e10 1e30 _
>>>>> 50 2e10 9.22337e18 9.22337e18
>>>>> Oops. What happened?
>>>>>     MAX <. _
>>>>> 9.22337e18
>>>>> Since one of the arguments is a float, (<.) casts both to floats, and
>>>>> takes the minimum. But the closest floating-point number to MAX is
>>>>> (MAX+1), and that number's floor (MAX+1) isn't representable as an
>>>>> integer--it's one too big. We didn't have this problem with MIN,
>> since
>>>>> it is exactly a negative power of two.
>>>>>     <. MAX+1
>>>>> 9.22337e18
>>>>>
>>>>> We'll make a test case that contains numbers close to both bounds.
>> I've
>>>>> included the addition of 0.5 in t so I can focus on the floor
>> function
>>>>> from now on. The type of our result is a float, so we failed.
>>>>>     t =. 0.5 + __,_,~ (MIN,0,MAX) +/(,@:) i:1e5
>>>>>     3!:0 ([: <. MIN>.MAX<.]) t
>>>>> 8
>>>>>
>>>>> We can fix the problem by using exact integers, but it's extremely
>> slow.
>>>>> However, it serves as a good answer key. The ("0) is there for a
>>>>> reason--otherwise the big array of exact integers tends to flood RAM.
>>>>>     fl_e =: (MIN>.MAX<.<.)&.:x:"0  NB. exact floor
>>>>>     3!:0 key =. fl_e t
>>>>> 4
>>>>>     10 (6!:2) 'fl_e t'
>>>>> 3.62018
>>>>>
>>>>> If we use a number small enough that its floating-point
>> representation
>>>>> is equal to a 64-bit integer, then we can force our answer to be a
>>>>> float, but it's not correct since the results are sometimes too
>> small.
>>>>> If that doesn't matter and speed is critical, this is the right
>> method
>>>>> to use.
>>>>>     MAX1 =. MAX - 512
>>>>>     fl_f =: [: <. MIN>.MAX1<.]  NB. fast floor
>>>>>     3!:0 fl_f t
>>>>> 4
>>>>>     key -: fl_f t
>>>>> 0
>>>>>     10 (6!:2) 'fl_f t'
>>>>> 0.0179205
>>>>>
>>>>> Finally, my solution. It's not particularly elegant, but it is
>> correct
>>>>> and has good performance. We reduce all the values larger than MAX to
>>>>> zero, then clamp on the minimum side and take the floor. For the
>> values
>>>>> that we removed, we add MAX back in. The comparison (<:&MAX) is only
>>>>> computed once to save a little time.
>>>>>     fl =: ((MAX*-.@]) + [: <. MIN>.*) <:&MAX
>>>>>     key -: fl t
>>>>> 1
>>>>>     10 (6!:2) 'fl t'
>>>>> 0.0426181
>>>>> It's critical to use (<:) rather than (<) to test whether numbers are
>>>>> acceptable even though it fails MAX, which wouldn't break (<.).
>> That's
>>>>> because comparisons cast their arguments to floats before comparing,
>> so
>>>>>     MAX < MAX+1
>>>>> 0
>>>>>
>>>>> Maybe there's a quicker solution to be found. The following rounds
>>>>> towards zero quickly by negating all the positive numbers, and
>> restoring
>>>>> their signs later. However, adding in the cases to make it equal to
>> (<.)
>>>>> on small numbers removes its advantage.
>>>>>     fl_o =: (] * MIN <.@:>. *) -@:*  NB. floor towards zero
>>>>>     fl_o _4.6 _3 _2.8 _1.2 3.4 5.8 9
>>>>> _5 _3 _3 _2 4 6 9
>>>>>     10 (6!:2) 'fl_o t'
>>>>> 0.0324293
>>>>>
>>>>> Any takers?
>>>>>
>>>>> Marshall
>>>>>
>> ----------------------------------------------------------------------
>>>>> For information about J forums see
>> http://www.jsoftware.com/forums.htm
>> ----------------------------------------------------------------------
>>>>> For information about J forums see
>> http://www.jsoftware.com/forums.htm
>>>> ----------------------------------------------------------------------
>>>> For information about J forums see http://www.jsoftware.com/forums.htm
>>>>
>>>
>>>
>>> --
>>>
>>> Devon McCormick, CFA
>>>
>>> Quantitative Consultant
>>> ----------------------------------------------------------------------
>>> For information about J forums see http://www.jsoftware.com/forums.htm
>> ----------------------------------------------------------------------
>> For information about J forums see http://www.jsoftware.com/forums.htm
>>
> ----------------------------------------------------------------------
> For information about J forums see http://www.jsoftware.com/forums.htm



More information about the Programming mailing list