interesting, so what is the type of pi/100 in matlab? is it a rational or
something? or do they do some kind of hacky approximation in their
(in)equality tests?
andrew
On Thursday, 24 April 2014 10:28:16 UTC-3, Peter Simon wrote:
>
> *Matlab behavior:*
>
> >> format long
> >> x = 0:pi/100:pi;
> >> length(x)
>
> ans =
>
> 101
>
> >> x(1) - 0
>
> ans =
>
> 0
>
> >> x(end) - pi
>
> ans =
>
> 0
>
> >> diff([max(diff(x)), min(diff(x))])
>
> ans =
>
> -4.440892098500626e-16
>
>
> *Julia behavior:*
>
> julia> x = 0:pi/100:pi;
>
> julia> length(x)
> 100
>
> julia> x[1]-0
> 0.0
>
> julia> x[end]-pi
> -0.031415926535897754
>
> julia> diff([maximum(diff(x)), minimum(diff(x))])
> 1-element Array{Float64,1}:
> -4.44089e-16
>
>
> I highlighted the important differences in red. IMO the Matlab behavior
> is more intuitive. If you choose an increment that "very nearly"
> (apparently within a few ulp, as stated by Stefan) evenly divides the
> difference between the first and last entries in the colon notation, Matlab
> assumes that you want the resulting vector to begin and end with exactly
> these values. This makes sense in most situations, IMO. E.g. when using a
> range to specify sample points in a numerical integration scheme, or when
> using a range to specify x-coordinates in an interpolation scheme. In both
> of these cases you don't want the last point to be adjusted away from what
> you specified using the colon notation.
>
> --Pete
>
> --Peter
>
> On Thursday, April 24, 2014 5:28:19 AM UTC-7, andrew cooke wrote:
>>
>>
>> how does matlab differ from julia here? (i though the original problem
>> was related to "fundamental" issues of how you represent numbers on a
>> computer, not language specific).
>>
>>
>> On Thursday, 24 April 2014 02:24:59 UTC-3, Peter Simon wrote:
>>>
>>> Thanks, Simon, that construct works nicely to solve the problem I posed.
>>>
>>>
>>> I have to say, though, that I find Matlab's colon range behavior more
>>> intuitive and generally useful, even if it isn't as "exact" as Julia's.
>>>
>>> --Peter
>>>
>>> On Wednesday, April 23, 2014 7:17:23 PM UTC-7, Simon Kornblith wrote:
>>>>
>>>> pi*(0:0.01:1) or similar should work.
>>>>
>>>> On Wednesday, April 23, 2014 7:12:58 PM UTC-4, Peter Simon wrote:
>>>>>
>>>>> Thanks for the explanation--it makes sense now. This question arose
>>>>> for me because of the example presented in
>>>>> https://groups.google.com/d/msg/julia-users/CNYaDUYog8w/QH9L_Q9Su9YJ :
>>>>>
>>>>> x = [0:0.01:pi]
>>>>>
>>>>> used as the set of x-coordinates for sampling a function to be
>>>>> integrated (ideally over the interval (0,pi)). But the range defined in
>>>>> x
>>>>> has a last entry of 3.14, which will contribute to the inaccuracy of the
>>>>> integral being sought in that example. As an exercise, I was trying to
>>>>> implement the interpolation solution described later in that thread by
>>>>> Cameron McBride: "BTW, another possibility is to use a spline
>>>>> interpolation on the original data and integrate the spline evaluation
>>>>> with quadgk()". It seems that one cannot use e.g. linspace(0,pi,200)
>>>>> for
>>>>> the x values, because CoordInterpGrid will not accept an array as its
>>>>> first
>>>>> argument, so you have to use a range object. But the range object has a
>>>>> built-in error for the last point because of the present issue. Any
>>>>> suggestions?
>>>>>
>>>>> Thanks,
>>>>>
>>>>> --Peter
>>>>>
>>>>> On Wednesday, April 23, 2014 3:24:10 PM UTC-7, Stefan Karpinski wrote:
>>>>>>
>>>>>> The issue is that float(pi) < 100*(pi/100). The fact that pi is not
>>>>>> rational – or rather, that float64(pi) cannot be expressed as the
>>>>>> division
>>>>>> of two 24-bit integers as a 64-bit float – prevents rational lifting of
>>>>>> the
>>>>>> range from kicking in. I worried about this kind of issue when I was
>>>>>> working on FloatRanges, but I'm not sure what you can really do, aside
>>>>>> from
>>>>>> hacks where you just decide that things are "close enough" based on some
>>>>>> ad
>>>>>> hoc notion of close enough (Matlab uses 3 ulps). For example, you can't
>>>>>> notice that pi/(pi/100) is an integer – because it isn't:
>>>>>>
>>>>>> julia> pi/(pi/100)
>>>>>> 99.99999999999999
>>>>>>
>>>>>>
>>>>>> One approach is to try to find a real value x such that
>>>>>> float64(x/100) == float64(pi)/100 and float64(x) == float64(pi). If any
>>>>>> such value exists, it makes sense to do a lifted FloatRange instead of
>>>>>> the
>>>>>> default naive stepping seen here. In this case there obviously exists
>>>>>> such
>>>>>> a real number – π itself is one such value. However, that doesn't quite
>>>>>> solve the problem since many such values exist and they don't
>>>>>> necessarily
>>>>>> all produce the same range values – which one should be used? In this
>>>>>> case,
>>>>>> π is a good guess, but only because we know that's a special and
>>>>>> important
>>>>>> number. Adding in ad hoc special values isn't really satisfying or
>>>>>> acceptable. It would be nice to give the right behavior in cases where
>>>>>> there is only one possible range that could have been intended (despite
>>>>>> there being many values of x), but I haven't figured out how determine
>>>>>> if
>>>>>> that is the case or not. The current code handles the relatively
>>>>>> straightforward case where the start, step and stop values are all
>>>>>> rational.
>>>>>>
>>>>>>
>>>>>> On Wed, Apr 23, 2014 at 5:59 PM, Peter Simon <[email protected]>wrote:
>>>>>>
>>>>>>> The first three results below are what I expected. The fourth
>>>>>>> result surprised me:
>>>>>>>
>>>>>>> julia> (0:pi:pi)[end]
>>>>>>> 3.141592653589793
>>>>>>>
>>>>>>> julia> (0:pi/2:pi)[end]
>>>>>>> 3.141592653589793
>>>>>>>
>>>>>>> julia> (0:pi/3:pi)[end]
>>>>>>> 3.141592653589793
>>>>>>>
>>>>>>> julia> (0:pi/100:pi)[end]
>>>>>>> 3.1101767270538954
>>>>>>>
>>>>>>> Is this behavior correct?
>>>>>>>
>>>>>>> Version info:
>>>>>>> julia> versioninfo()
>>>>>>> Julia Version 0.3.0-prerelease+2703
>>>>>>> Commit 942ae42* (2014-04-22 18:57 UTC)
>>>>>>> Platform Info:
>>>>>>> System: Windows (x86_64-w64-mingw32)
>>>>>>> CPU: Intel(R) Core(TM) i7 CPU 860 @ 2.80GHz
>>>>>>> WORD_SIZE: 64
>>>>>>> BLAS: libopenblas (USE64BITINT DYNAMIC_ARCH NO_AFFINITY)
>>>>>>> LAPACK: libopenblas
>>>>>>> LIBM: libopenlibm
>>>>>>>
>>>>>>>
>>>>>>> --Peter
>>>>>>>
>>>>>>>
>>>>>>
