Julia: csgn() function?

[JaredCrean2]

Currently:

julia> sign(-1.0)
-1.0

Should this return an integer instead?

[yuyichao]

julia> sign(1 + 1im)
0.7071067811865475 + 0.7071067811865475im

[yuyichao]

And I believe there's a few places that relies on this for type stability.

[JaredCrean2]

In that case, should julia implement csgn as well?

[TotalVerb]

Julia's sign _is_ a csgn function. I think it's reasonable for sign to return the type of its input argument. Note that your desired behaviour is already implemented as cmp(x, 0).

[JaredCrean2]

Julia's sign is a csgn function

According to Wikipedia's definition, sgn(1 + 1im) returns (1+ 1im)/sqrt(2), but csgn(1 + 1im) should return 1.

your desired behaviour is already implemented as cmp(x, 0)

That's perfect, thanks.

[TotalVerb]

@JaredCrean2 Ah, my mistake. Perhaps we should have a csgn function then.

[JaredCrean2]

For both real and complex numbers, csgn could return an integer (assuming this holds for all other types of numbers), which would solve my original concern.

[musm]

csign() would be more consistent. for type stability it seems to me that the current behavior is preferable.

[stevengj]

Is there any practical need for csgn? If it is rarely used, anyone who needs it can just implement it themselves.

[JaredCrean2]

The way I noticed this originally was checking sign(a) == sign(b), where a and b are floating point numbers. Having a function that returns integers in all cases is better for this.

[TotalVerb]

                                                                                  Perhaps I am missing something obvious but I don't see why integers are a better return type here.                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                         From: Jared CreanSent: Monday, September 19, 2016 12:02To: JuliaLang/juliaReply To: JuliaLang/juliaCc: Fengyang Wang; CommentSubject: Re: [JuliaLang/julia] csgn() function? (#18549)The way I noticed this originally was checking sign(a) == sign(b), where a and b are floating point numbers.  Having a function that returns integers in all cases is better for this.

—You are receiving this because you commented.Reply to this email directly, view it on GitHub, or mute the thread.

{"api_version":"1.0","publisher":{"api_key":"05dde50f1d1a384dd78767c55493e4bb","name":"GitHub"},"entity":{"external_key":"github/JuliaLang/julia","title":"JuliaLang/julia","subtitle":"GitHub repository","main_image_url":"https://cloud.githubusercontent.com/assets/143418/17495839/a5054eac-5d88-11e6-95fc-7290892c7bb5.png","avatar_image_url":"https://cloud.githubusercontent.com/assets/143418/15842166/7c72db34-2c0b-11e6-9aed-b52498112777.png","action":{"name":"Open in GitHub","url":"https://github.com/JuliaLang/julia"}},"updates":{"snippets":[{"icon":"PERSON","message":"@JaredCrean2 in #18549: The way I noticed this originally was checking sign(a) == sign(b), where a and b are floating point numbers. Having a function that returns integers in all cases is better for this."}],"action":{"name":"View Issue","url":"https://github.com/JuliaLang/julia/issues/18549#issuecomment-248036232"}}}

[JaredCrean2]

Comparing floating point numbers for strict equality is never a good idea.

[musm]

yes but it's either 1.0 or -1.0 so there's no problem

[JaredCrean2]

Because the complex definition is sign = z/abs(z), whether or not the function return exactly 1.0 or approximately 1.0 is an implementation detail.

[yuyichao]

You can't do better for complex anyway. If you want to support complex numbers in that way, you should just do sign(real(a)) == sign(real(b)).

[yuyichao]

yes but it's either 1.0 or -1.0 so there's no problem

FWIW, it can be 0 too.

Basically, if you only want to compare the sign of real numbers, then there's no problem. It's a bug if the implementation doesn't return exactly 1, -1, or 0.

If you want to support complex numbers in general, then you are not dealing with 1, -1, 0 to begin with and you should compare them approximately in your code.

And if you actually want to compare only the sign of the real part, you are not using the right function and you should instead write sign(real(x)). It is true that we can provide csgn for this but as @stevengj mentioned it doesn't seem to be useful enough.

[JaredCrean2]

It's a bug if the implementation doesn't return exactly 1, -1, or 0.

Then it should return an integer. It can either return an integer and be exact or return a floating point and be approximate, but attempting to return an exact floating point is bad practice.

[musm]

1, -1, 0 can be stored exactly in float ...

[JaredCrean2]

For all floating point types that might ever exist in the language? A perfectly valid float type would be one that has similar range as a Float64, but can't represent 1 exactly (perhaps it can represent 1+eps exactly, or something like that) and randomly decides to round up or down to the nearest representable number. If you want sign to return the type of the input, you can't make any assumptions about its representation and still get exactness, or even reliability.

[TotalVerb]

sign is _defined_ to return 1. If some exotic type doesn't have a one value, then it simply can't support sign. It would be invalid for it to randomly round up or down to the nearest representable number.

[yuyichao]

Well, for such types one cannot return exactly 1 either and that's not a good argument to make one(Float64) returning a integer, just that such types are missing many important properties in general and you'd better have a very good reason to use it.

[JaredCrean2]

For a real number, the sign function maps from the value to 1, 0, or -1, _independent of the representation of the number_. That fact that the current implementation fails for some representations is a bug.

[yuyichao]

independent of the representation of the number

I don't see why this is the case and it's pointless to break ppls code without an actual use case.

[simonster]

Are there any examples of number types where multiplicative identity is irrepresentable?

[TotalVerb]

Types with physical units. But in that case, sign still can return the multiplicative identity in the wrapped type, like one.

[stevengj]

@JaredCrean2 the idea that floating-point numbers are "less exact" than integer types is a pernicious myth. (I like this presentation by Kahan, which has much to say on the subject: http://www.cs.berkeley.edu/~wkahan/JAVAhurt.pdf ... start at page 25, and see especially page 33)

Integers (up to the max significand) are exactly represented in floating point, and there is nothing wrong with returning them here.

[stevengj]

And yes, all floating-point types that have ever existed in any language on any hardware, or any sane fp type that could ever conceivably be created, will represent ±1 and 0 exactly. By construction, the set of floating-point numbers are integers × base^exponents for some fixed base and some range of signed integers (significands) and signed exponents. (Floating-point implementations have differed in base, in rounding behaviors for operations, and many other details, but not in that basic construction — something that is not integers × base^exponents is not floating point.)

[stevengj]

I like this quote from Kahan:

Many a textbook asserts that a floating-point number represents the set of all numbers that differ from it by no more than a fraction of the difference between it and its neighbors with the same floating-point format. This figment of the author’s imagination may influence programmers who read it but cannot otherwise affect computers that do not read minds. A number can represent only itself, and does that perfectly.

[JaredCrean2]

I think the link is dead. I get "Page Not Found" saying they recently redesigned their website.

[stevengj]

The link works for me.

[tkelman]

@stevengj you might have it in your browser cache or something? https://people.eecs.berkeley.edu/~wkahan/JAVAhurt.pdf

Couldn't you have some exotic scaled logarithmic number format? And didn't one of the recent Unum proposals include exact representation of some arbitrarily chosen set of points along with their reciprocals? Not sure why you would, but it might be valid to leave out an exact representation of 1 in some of those schemes? Probably wouldn't call those conventional floating point though.

[stevengj]

Those aren't floating-point formats.

[PallHaraldsson]

@tkelman: "And didn't one of the recent Unum proposals include exact representation of some arbitrarily chosen set of points along with their reciprocals?" Yes, Unum 2.0. And yes, they had +1, -1 (+- integers at least up to 10, also exact).

@JaredCrean2: "For all floating point types that might ever exist in the language? A perfectly valid float type would be one that has similar range as a Float64, but can't represent 1 exactly (perhaps it can represent 1+eps exactly, or something like that)"

I can't think of any current [or future] [floating point] number type, where -1 and 1 (but not 0?!) isn't exact. At least not Unums 1.0 or 2.0 and his variants (and wouldn't it be desirable to any type?), [but I'm not familiar with https://github.com/JuliaDiff/DualNumbers.jl (are they something like "1+eps" as they take the form x+y*ɛ? Still y=0 would then do. Then there is https://github.com/JuliaDiff/HyperDualNumbers.jl that I know even less about.. Is it useful to anyone to make a type with only non-exact +1 and -1? Already done?]

[negabinary (base −2), https://en.wikipedia.org/wiki/UMC_(computer) and ternary (base 3) https://en.wikipedia.org/wiki/Setun computers have (all) exact integers; I haven't thought through floating point in those systems, or if they had..]

Probably false alarm: I got same result for sign(1 + 1im), but it took a very long time.. (this "compile"/first time). I'm not sure if this is relevant to anything (e.g. "And I believe there's a few places that relies on this for type stability."). Maybe my machine was just under an unusual load (just after: "load average: 0,66, 0,78, 0,80", then very quick on first use just after restarting Julia with "load average: 0,52, 0,56, 0,68").

[JaredCrean2]

I agree that all the common floating point types do represent -1, 0 and 1 exactly, but I don't understand the motivation for relying in implementation specific behavior. More generally, returning a floating point value signals to the user that the value is the result of some computation and is subject to the usual floating point roundoff issues.

independent of the representation of the number

I don't see why this is the case and it's pointless to break ppls code without an actual use case.

Because that is the definition of the mathematical function. It says to return exactly 1, 0 or -1, not something close. I do see the point you are making with one(Float64), but I see that function as returning the closest representable value to 1, and any error that result from not returning exactly 1 is within the bounds of floating point error for that datatype. In fact, the current behavior could lead to bugs, such as sign(a) == sign(b) returning false even when it should return true when a and b have different representations.

[yuyichao]

Because that is the definition of the mathematical function

It doesn't talk about representation at all.

In fact, the current behavior could lead to bugs, such as sign(a) == sign(b) returning false even when it should return true when a and b have different representations.

Give an example.

[yuyichao]

but I don't understand the motivation for relying in implementation specific behavior.

Because most people are not writing code for a non-existing implementation of numbers. Making it easier to work with existing implementations is much more important.

[martinholters]

Note:

julia> sign(-0.0)
-0.0

So always returning an integer would loose the sign of -0.0, while current sign keeps it. (Although I wonder whether this is relied on anywhere.)

I like the current behavior of returning the input type if possible, choosing something reasonable otherwise (e.g. sign(::Complex{Int64})::Complex{Float64}). Any ever-so-wonky number type could easily follow that, I guess.

If someone really needs a sign function that returns an integer, a definition like sign{T}(::Type{T}, x::Real) = x == zero(x) ? zero(T) : (x < zero(x) ? -one(T) : one(T)) to be able to do sign(Int, x) might make sense.

[PallHaraldsson]

I'm not sure this all matters in this context, but see:

julia> sign(-0.0)
-0.0

@stevengj: "will represent ±1 and 0 exactly." I was sure on the former, but not for 0. Are you sure about 0 (and how does that fit Kahan: "A number can represent only itself, and does that perfectly.", as there is -0.0 and +0.0). There can not be an exact -0.0 that is different from +0.0, mathematically, unless: I always understood -0.0 to be close to 0 from below (a kind of limit); Unum 2.0 does away with separate represenations/unifies +-0 (unlike Unum 1.0 that is a superset of IEEEE floating point).

[stevengj]

@JaredCrean2, saying that small integers are represented exactly is not an implementation detail, it is part of the definition of what floating-point arithmetic means.

returning a floating point value signals to the user that the value is the result of some computation and is subject to the usual floating point roundoff issues.

No it doesn't.

It says to return exactly 1, 0 or -1, not something close

1.0 is exactly 1, not "something close."

In fact, the current behavior could lead to bugs, such as sign(a) == sign(b) returning false even when it should return true when a and b have different representations.

No it can't, because 1.0 == 1 == 1.0f0 == big(1.0) == 1.0+0im: all represent the same value.

@PallHaraldsson, +0.0 == -0.0; they both represent (exactly) the number zero. The sign bit indicates some _additional_ information beyond the value, but does not _change_ the value.

[JaredCrean2]

Do you have a precise definition of floating point in mind? Your earlier definition of integers * base^exponent doesn't seem to require that the integers be the set of all integers and, depending on the exponent values, might not capture the small integers exactly.

btw. the link you posted is working now (not sure what happened before...). Looks like an interesting read.

Edit: does -> doesn't

[stevengj]

@JaredCrean2, floating point numbers are (sign) * (significand with some number of digits in some base) * base ^ (signed exponent with some number of digits), aside from special values like Inf and NaN.

[JaredCrean2]

Ok, I agree to this definition of floating point numbers. Now, from this definition, how do you prove the small integers must be represented exactly?

[stevengj]

By the above definition, a zero exponent is included, hence any ±integer up to the maximum precision of the significand must be exactly representable. The significand has at least one digit, so ±1 and 0 are exactly representable.

[stevengj]

That is, by "signed significand with p digits in some base", I mean any integer of the form ±dₚ⋯d₂d₁ where the d's are arbitrary digits in that base.

[stevengj]

(This isn't just my definition; see e.g. _Numerical Linear Algebra_ by Trefethen and Bau, lecture 13, which provides an essentially equivalent definition except that they make the commonplace idealization that the exponent is an arbitrary integer, i.e. they ignore underflow and overflow, since those effects are somewhat orthogonal to other roundoff-error analyses.)

[stevengj]

I guess you could also define floating-point as something like ±0.dₚ⋯d₂d₁ × βᵉ, where e is the exponent and β is the base, and if the maximum |exponent| is < p, you couldn't get all integers up to the precision (though in practice this is _never_ the case for any real fp format — that would require an insanely small exponent precision). Even then, however, you could still represent ±1 unless the only allowed exponent is zero, in which case it isn't floating point (because the decimal point can't move).

[stevengj]

@JaredCrean2, I can only think that you're trying to imagine a floating-point system in which the minimum exponent is some positive number (so that 1 underflows), or the maximum exponent is some large negative number (so that 1 overflows). Maybe you could argue that this would still be "floating-point" arithmetic, but this has _never_ been the behavior of any real fp system. Nor, I would argue, will it _ever_ be the behavior of any future fp system. People using floating point always want to be able to represent numbers both > 1 and < 1, and it is absurd to suppose otherwise in our standard library.

If we ever come across a wacky number format in which 1 is not representable, we can always define sign(x::WackyNumber) = x < 0 ? -1 : x > 0 ? +1 : 0. But it doesn't make sense to base our treatment of _real_ fp formats on this nonexistent contingency.

[timholy]

If we ever come across a wacky number format in which 1 is not representable

Actually...

julia> using FixedPointNumbers

julia> typemax(Fixed{Int8,7})
0.992Q0f7

julia> sign(typemax(Fixed{Int8,7}))
ERROR: InexactError()
 in sign(::FixedPointNumbers.Fixed{Int8,7}) at ./number.jl:38

See https://en.wikipedia.org/wiki/Fixed-point_arithmetic

[JaredCrean2]

@stevengj I see, you're right. If the significant is required to be an integer in some base, then the small integers will be exactly representable. I was thinking about a more general mapping from the significand bits to their values (something like value(significant) = value_in_binary(significant) + 1/2), but that doesn't fit the definition of floating point

@timholy That's a great example.

[stevengj]

@timholy, I was talking about floating-point formats, but you're quite right that fixed point is different. For fixed point, a sign definition returning ±1 makes sense.

[andreasnoack]

For fixed point, a sign definition returning ±1 makes sense.

It might have surprising consequences

julia> x = Fixed{Int8,7}(0.5)
Fixed{Int8,7}(0.5)

julia> x*1
Fixed{Int8,7}(-0.5)

julia> x*(-1)
Fixed{Int8,7}(-0.5)

[TotalVerb]

I'm starting to become more convinced that a universal One, Zero, and Negative{One} types would be a good idea. Some numeric types can't provide their own equivalents.

[yuyichao]

FWIW, I don't think it's a good idea to make sign type unstable.

[TotalVerb]

The type instability would only occur on incredibly obscure numeric types, which are so exotic as to be practically useless for any reasonable use case. Obviously such sign methods would not be provided by Base, but by packages providing these exotic numeric types.

[yuyichao]

I'm not sure if a One type is useful but it's not really relevant here. If anything, this is better handled by a Sign type which can be used for special number types.

[TotalVerb]

OCInt2, for "one's complement 2-bit integer"? that would support -0, 0, 1, and -1. Though I think this should really be provided in a package.

[timholy]

It might have surprising consequences

On recent FixedPointNumbers that x*1 throws an InexactError rather than weirdly changing sign. Haven't dug down to see why.

[stevengj]

@yuyichao, I don't know why you think this would make the sign function type-unstable. For floating-point types it could continue to return floats, and for fixed-point types it could return Int or whatever. Or it could always return typeof(one(T)). This is type stable.

[TotalVerb]

@stevengj He was referring to my suggestion to make sign return special constants (like Irrational) for types that don't have 1, which is indeed type unstable.

[yuyichao]

Right, I was commenting on returning One(), Zero(), etc. I think it makes sense to keep the current behavior for all Base number types (and most other numbers that can represent 1), and returns whatever it makes sense for FixedPointNumbers (and possibly other special number types that can't represent 1).