Numerical problems

[Erich Neuwirth]

The following behavior is rather surprising:

julia> 100^10
7766279631452241920

No warning, no error for the integer overflow.
Are there any plans to change this?


With floating point arithmetic, things work.
julia> 100.0^10
1.0e20

But in teaching, no warning with integer overflow is rather dangerous.

[John Myles White]

As far as I recall, I don't think this is going to change anytime soon. If you want to run integer ops that might overflow, use BigInt's:

julia> load("bigint.jl")
julia> BigInt(100)^10

-- John

> --
>
>
>

[Gabor]

When trying to use Float64 values, I was surprised to observe this:

> 100.0^10
1.0e20

This works as expected.

> 100.^10
1661992960

Here the dot becomes part of the .^ elementwise operator.

Is other cases the number formats 100.0 and 100. are identical.
Do you consider this a problem, or it's no big deal?

[John Myles White]

Buf. Looks like the . operator introduces a parsing ambiguity. But it also seems like .^ should perhaps not even be defined in the example you've given, since it's operating on a single value rather than a vector.

 -- John

--
 
 
 

[Ohad]

Does the same apply to:

julia> 9223372036854775807+1

-9223372036854775808

or is this an issue?

Why isn't it going to change, btw?

[Viral Shah]

There is an issue open to track this:

https://github.com/JuliaLang/julia/issues/855

-viral

[Stefan Karpinski]

I don't understand why this behavior is surprising. We all know that integers, as represented by computers, aren't Platonic mathematical integers and don't behave like them. Checking for overflow everywhere by default is way too slow. In my mind this is the same as complaining about this:

julia> 0.1 + 0.2 == 0.3

false

Float64's aren't real numbers. Int64's aren't integers. Julia isn't in the business of lying to you and trying to make it seem like these things are true even though they aren't. It is in the business of giving you easy access to efficient operations on these common hardware-supported imperfect representations of numbers.

If you're teaching, this is an excellent opportunity to explain to students how computers represent integers. Here's a very instructive loop:

julia> for k=0:10

           println(bits(100^k))

       end

0000000000000000000000000000000000000000000000000000000000000001

0000000000000000000000000000000000000000000000000000000001100100

0000000000000000000000000000000000000000000000000010011100010000

0000000000000000000000000000000000000000000011110100001001000000

0000000000000000000000000000000000000101111101011110000100000000

0000000000000000000000000000001001010100000010111110010000000000

0000000000000000000000001110100011010100101001010001000000000000

0000000000000000010110101111001100010000011110100100000000000000

0000000000100011100001101111001001101111110000010000000000000000

0000110111100000101101101011001110100111011001000000000000000000

0110101111000111010111100010110101100011000100000000000000000000

I suspect that you'd be doing your students a huge favor by teaching them early on that ints aren't integers and floats aren't real numbers. It will save them from a lot of confusion and bad computational thinking.

--
 
 
 

[Keno Fischer]

I still think that we should have optional overflow exceptions.

--
 
 
 

[Stefan Karpinski]

Absolutely — there are lots of situations where you need a way to check for overflow. The parse_int function is a good example. That's what the issue linked above is about.

--
 
 
 

[John Myles White]

I also think optional exceptions could be useful.

I think this is another spot where Julia's core philosophy differs a lot from other dynamic languages like Python or Matlab, which are largely designed to conceal the fact that you're using a real machine with substantial limitations in its ability to represent real numbers and not a machine as defined by a traditional theoretician. Driving that home somewhere in the manual is probably worth doing. Something like, "Julia is able to produce code that runs faster in part because Julia insists on representing problems using units of computation like integer and floating point arithmetic that are hardcoded by modern chips. Lots of other dynamic languages are less efficient because they invest a lot of computational resources in hiding the underlying machine from the user."

 -- John

--
 
 
 

[Jasper]

For the record there seems to be a bug in
power_by_squaring at intfuncs.jl:80(e519c91d571a7caf9dc96e8ca355d6e5d958bf83)
that can be seen with `2^-2` or `^(Integer,NegativeInteger)` generally.

[Stefan Karpinski]

You mean that if throws a DomainError? No, that's intentional:

https://github.com/JuliaLang/julia/issues/745

I'm not entirely convinced this is the best behavior, but it was the upshot of that issue discussion. The error message certainly ought to be better.

[Stefan Karpinski]

Definitely. I think we are in need of a Philosophy page in the manual explicitly stating things like this. It is a very different approach than what other dynamic languages take, and is much closer to C.

--
 
 
 

[John Myles White]

If you have a place where you'd like to put that, I'll throw in these tidbits as I collect them.

 -- John

--
 
 
 

[Jasper]

On Tue, 21 Aug 2012 13:49:03 -0400
Stefan Karpinski <ste...@karpinski.org> wrote:

> You mean that if throws a DomainError? No, that's intentional:
>
> https://github.com/JuliaLang/julia/issues/745
>
> I'm not entirely convinced this is the best behavior, but it was the
> upshot of that issue discussion. The error message certainly ought to
> be better.

I see, probably should have searched the issues for the function
name&error before mentioning it..

[Stefan Karpinski]

No, no. That's good. It's a lousy error since it gives no indication of why there's a problem.

--

[Neuwirth Erich]

But even  the old 8086 had flags for overflow. So the argument that julia behaves this way because it is "close to the hardware" does not seem valid for me.

--
 
 
 

[Elliot Saba]

The 8086 has flags for overflow, but C, Assembly, etc. didn't check that register for overflow after every operation.  Stefan's point is that it is the checking after every operation that people expect and is unreasonable.  Manual checking after certain operations the programmer, not the language deems "suspect" is what is being promoted here, which is exactly how many languages (such as Assembly and C) have done things for decades.

-E

--
 
 
 

[Jeff Bezanson]

We are planning to do at least 2 things about this: providing an easy
way to add an overflow check to a given operation, and adding a switch
that will turn on overflow detection globally at a performance cost.

When I was 10 or 11 I learned that maxint = 32767, and that doesn't
seem to have done me much visible harm.

> --
>
>
>

[Neuwirth Erich]

BTW: Java has very similar problems!
And therefore Scala and Clojure also.

> --
>
>
>

[Jeffrey Sarnoff]

thanks for that, Jeff.

[dslate]

For what it's worth, here is a perspective on this issue from the point of

view of an experienced computer programmer (nearly 50 years) but a relative

newcomer to Julia (last few months):

On one hand, I agree with Stefan Karpinski, John Myles White, et al, that part

of the process of growing up involves abandoning the cherished illusions of

childhood and accepting some hard truths: there is no Santa Claus, no free

lunch, and 0.1 + 0.2 does not equal 0.3.

However, since Julia aspires to be "a high-level, high-performance dynamic

programming language for technical computing, with syntax that is familiar to

users of other technical computing environments", it would be desirable that

not only the syntax but also the semantics of Julia, in particular the basic

arithmetic operations, not only resemble those of other technical computing

environments, but also conform as much as possible, within some limitations

imposed by the need for high performance, to the reasonable expectations of

Julia's prospective users, many of whom may be "end-user" scientists and

engineers first and computer programmers second.  That is presumably why

Julia's ordinary division operator, '/', converts its args to Float64, so

that, unlike C (and even some dynamic languages, such as Python < version

3.0), 5/2 yields 2.5 and not 2.

Perhaps the best default representation in a technical computing environment

for real numbers, including literals without decimal points, is Float64.

Speaking as someone who studied physics before switching to computers, that

choice arguably produces the fewest surprises and "gotchas" in the results of

arithmetic expressions.  Note that this is what Lua, which tries to pack lots

of expressive power and performance into a very small language, does; in

describing Lua's few built-in data types, the book "Programming in Lua"

explains why separate integer types are unnecessary.

However, as someone whose first programming experience was with assembly

language for an IBM 7090 mainframe, and who later wrote lots of assembly code

for 8-bit microprocessors such as MOS Technology's 6502, I certainly

appreciate the utility of low-level bit and byte fiddling.  Using small

integer types when appropriate saves memory space over double precision

floats.  In the "old days", back when "maxint = 32767" was first engraved on

Jeff Bezanson's grey matter, it also saved computing time, although with

modern hardware the speedup is less evident.

Regarding neuwirthe's example of 100^10 overflowing silently, I think a

legitimate case could be made that, like the '/' operator, the '^' operator

ought to convert its args to Float64 and produce a Float64 result, both

because of the likelihood of an integer overflow "gotcha", and also because

such an expression just "feels" (to me anyway) that it should take place in

the real rather than integer domain.  Perhaps a problem with this point of

view is that unlike 100^10, 2^7 does look somewhat like a synonym for bit-

shifting (1 << 7).

I'm not arguing here to actually change Julia's specs for '^', just throwing

in my 2 cents worth of philosophy on the subject.

-- Dave Slate

[Jeffrey Sarnoff]

worth at least 3 cents

[dslate]

Thanks Jeffrey, although that might just be due to inflation :)

-- Dave Slate

[Steve Jaffe]

I too cannot see a reason, unrelated to “implementation details”, that all numerical values would not be floating-point numbers. End-users who are scientists rather than programmers will definitely expect this and be extremely upset when their expectations are contradicted.

 

I agree 110% with Jeffrey Sarnoff that anything else completely contradicts Julia’s stated aspirations.

 

Steve

--
 
 
 

[John Cowan]

On Thu, Aug 23, 2012 at 10:01 PM, Steve Jaffe <sja...@riskspan.com> wrote:

> I too cannot see a reason, unrelated to “implementation details”, that all
> numerical values would not be floating-point numbers. End-users who are
> scientists rather than programmers will definitely expect this and be
> extremely upset when their expectations are contradicted.

They're doomed to be upset anyway when they discover that floats don't
behave like real numbers.

In any case, integers are still needed for indices and discrete
mathematics (e.g. discrete statistics).

--
GMail doesn't have rotating .sigs, but you can see mine at
http://www.ccil.org/~cowan/signatures

[dslate]

Actually, integers aren't really needed for indexing, since  IEEE double precision floats (Float64) represent integers perfectly well up to about 10^15.

-- Dave Slate

[Alessandro "Jake" Andrioni]

On 23 August 2012 23:01, Steve Jaffe <sja...@riskspan.com> wrote:
> I too cannot see a reason, unrelated to “implementation details”, that all
> numerical values would not be floating-point numbers. End-users who are
> scientists rather than programmers will definitely expect this and be
> extremely upset when their expectations are contradicted.

As someone who ocasionally writes code for friends who work in
combinatorics, having a good and efficient way to represent integers
is essential. (FP multiplication is still slower in the hardware
level, AFAIK) (and probably other operations are also a bit slower)

[John Cowan]

On Thu, Aug 23, 2012 at 11:04 PM, dslate <rusti...@gmail.com> wrote:

> Actually, integers aren't really needed for indexing, since IEEE double
> precision floats (Float64) represent integers perfectly well up to about
> 10^15.

True. However, hardware needs an integer to actually determine the
address of an array element, and float-to-int conversion is not
efficient *at all*.

Specifically, it needs to set the FPU mode to "truncate", which has
the effect of flushing the FP pipeline. Bad news. If you just leave
the FP mode alone, you get much more performance, but you get very
unpredictable (i.e. undebuggable) results if the value happens not to
be an integer.

Since Lua isn't all about performance, being a bytecode interpreter,
it doesn't have to worry about this cost.

[dslate]

I agree, integers can be useful, and I'm not actually advocating getting rid of them.  The main point of my post was that scientific and technical users may expect expressions involving numbers to behave as if the numbers were treated like floats more often than Julia actually does.  For example,

in R, "10^100" evaluates to "1e+100", and there is no issue of an integer overflow  Of course Julia doesn't want to emulate R in general, but if it aims to eventually replace R these points may confuse users. 

As to Lua, it has a pretty good JIT called LuaJIT2, although the nature of the language probably limits its performance to rather less than what Julia is aiming for.

[Stefan Karpinski]

On Thu, Aug 23, 2012 at 10:01 PM, Steve Jaffe <sja...@riskspan.com> wrote:

 

I too cannot see a reason, unrelated to “implementation details”, that all numerical values would not be floating-point numbers. End-users who are scientists rather than programmers will definitely expect this and be extremely upset when their expectations are contradicted.

 

Julia is all about the implementation details and not only exposing them, but also giving you sufficient control over them. In this, Julia is much more like C than any other high-level language, especially the traditional dynamic ones. That is why we have all sorts of numerical types: 8, 16, 32, 64, 128 bits, signed, unsigned, float. These implementation details are often very important when you're tackling sufficiently hard problems. We did make some judgement calls here: you don't really need to do 8-bit integer arithmetic — you really just need Int8s for compact storage. (That thinking took a *long* time coming too and is probably still contentious, but I think in practice its a good simplification.) And integers may not "in principle" be needed, but in practice, you really do need them. As John says, for indexing at least, but also for efficient loops — modern floating point ops may be fast, but they're still not as fast as integer ops. Languages like JavaScript and Lua that don't have and integer type end up having to guess when numbers are actually integers so that they can run fast. We're not into guessing games, though — I'd much rather let the programmer tell the compiler that what they want is an actual integer.

I agree 110% with Jeffrey Sarnoff that anything else completely contradicts Julia’s stated aspirations.

 

Since I'm the one who probably stated those aspirations (the provenance of that particular bit of copy is fairly shrouded in the mists of time now), I feel fairly comfortable saying that having integer ops overflow doesn't, at least in spirit, contradict those aspirations. If it does, maybe we should reword the aspirations ;-)

On Thu, Aug 23, 2012 at 8:23 PM, dslate <rusti...@gmail.com> wrote:

Regarding neuwirthe's example of 100^10 overflowing silently, I think a legitimate case could be made that, like the '/' operator, the '^' operator ought to convert its args to Float64 and produce a Float64 result, both because of the likelihood of an integer overflow "gotcha", and also because such an expression just "feels" (to me anyway) that it should take place in the real rather than integer domain.  Perhaps a problem with this point of view is that unlike 100^10, 2^7 does look somewhat like a synonym for bit-shifting (1 << 7).

It is definitely fairly compelling to make x^y return a float for integer x and y, and that was discussed in issue #745. The problem is that far more often than one wants to write 10^100, one wants to write x^2 or x^3 and have it just mean x*x or x*x*x. This dictates that when x is an integer, x^2 and x^3 are also integers. To that end, we decided to just say that if you want to take big powers and want to have them come out accurately, you should convert your values to floating point. Likewise, if you want to take negative powers and not get an error, use floating point. To be practical about it, if you're expecting to take really big powers and get really big answers, you should be working with floating point numbers in the first place. Writing 10^100 and getting a huge, correct answer is mostly a gimmick: it's nice to see it come out right in the repl, but if you're doing actual big computations, you should either be using floats or bigints, not 64-bit integers.

[John Cowan]

On Fri, Aug 24, 2012 at 1:30 AM, Stefan Karpinski <ste...@karpinski.org> wrote:

> It is definitely fairly compelling to make x^y return a float for integer x
> and y, and that was discussed in issue #745. The problem is that far more
> often than one wants to write 10^100, one wants to write x^2 or x^3 and have
> it just mean x*x or x*x*x.

I agree, so perhaps we should take advantage of having two plausible
operators, ^ and ** (right now ** just prints an error telling you to
use ^ instead). One of them would accept arguments of any numeric
type and return a Float64 or Complex128 result. The other would
accept only an Int on the right side, and return the same type as the
left.

So you would write 10^100 to get 1e100, but x**2 to get x * x. Or
vice versa, I'm not picky.

[John Cowan]

On Fri, Aug 24, 2012 at 1:44 AM, John Cowan <johnw...@gmail.com> wrote:

> The other would
> accept only an Int on the right side, and return the same type as the
> left.

That is, in accordance with the regular numeric type rules (so an Int8
** Int8 is an Int, not an Int8, just as with +.

[Stefan Karpinski]

On Fri, Aug 24, 2012 at 1:44 AM, John Cowan <johnw...@gmail.com> wrote:

On Fri, Aug 24, 2012 at 1:30 AM, Stefan Karpinski <ste...@karpinski.org> wrote:

> It is definitely fairly compelling to make x^y return a float for integer x
> and y, and that was discussed in issue #745. The problem is that far more
> often than one wants to write 10^100, one wants to write x^2 or x^3 and have
> it just mean x*x or x*x*x.

I agree, so perhaps we should take advantage of having two plausible
operators, ^ and ** (right now ** just prints an error telling you to
use ^ instead).  One of them would accept arguments of any numeric
type and return a Float64 or Complex128 result.  The other would
accept only an Int on the right side, and return the same type as the
left.

So you would write 10^100 to get  1e100, but x**2 to get x * x.  Or
vice versa, I'm not picky.

Eh. This seems unobvious and hard to remember — I'm sure I would forget which to use. I'd rather write float(x)^y if I want to do floating point exponentiation. Easy to remember and completely obvious to read.

[John Cowan]

On Fri, Aug 24, 2012 at 1:47 AM, Stefan Karpinski <ste...@karpinski.org> wrote:

>> So you would write 10^100 to get 1e100, but x**2 to get x * x. Or
>> vice versa, I'm not picky.
>
> Eh. This seems unobvious and hard to remember — I'm sure I would forget
> which to use.

The conceit is that ^ is mathematically correct exponentiation,
whereas ** is (logically, not necessarily in implementation) repeated
multiplication, so it uses the multiplication sign doubled.

[Jeffrey Sarnoff]

with a_posint, b_posint of type Union(Unsigned, Integer) and b_posint >= 0

     answer = a_posint^b_posint # using ^ not ** because
                                #   "brevity is the soul of wit"
                                #    -- William Shakespeare, "Hamlet"

     where the infinitely precise result happens to occur within
           typemin( promote_type(typeof(a_posint), typeof(b_posint)) ) ) .. typemax( same )
     there find answer's value as that value and answer's type as that type
     otherwise if the precise_result is negative and promote_type() is an unsigned type
           and precise_result >= typemin( signed integer type of same bitwidth )
     there find answer's value as that value and answer's type as that type
     otherwise     
     where the infinitely precise result happens to occur within
           typemin( morebits( promote_type(typeof(a_posint), typeof(b_posint)) ) ) .. typemax( same )
     there find answer's value as that value and answer's type as that type
     otherwise -- if precise_result is negative and morebits() is an unsigned type etc
     otherwise -- try morebits(morebits( _ )) while wider integer bitstypes are available etc

     still too large -- ok float it

     same for multiplication

because
  ints do not act as if they want floats to be ints
  but floats do act as if they want ints to be floats
  so float types tend to proliferate through ints present

anthropomorphic because 
  having taken care to specify a float free pair
    to be multiplied or non-negatively powered
  Julia, ever well-mannered, would reflect that

nitty-gritty because
  all the date and time information and all temporally aware constuents
  are integer valued and must remain so to ensure alacrity and accuracy
  (please, no more autofloat conversion parties .. pick a better theme)