I've researched this further, and it should be called signum(), because it can be defined for complex numbers, whereas sign() cannot.
signum() returns a number with the same phase as z, but with unit magnitude, except that if z is a complex zero, it returns z.
Again, it obeys

signum(z) * abs(z) ===>  z

The code below is not correct for nans, but gives the idea.
Further work is required (need to do complex abs() at the same time).
It's amazing how something so apparently trivial becomes complicated.

creal signum(creal z)
{
   return nonzero(z) ? z/abs(z) : z;
}

// returns false if z is zero or has a nan component
bit nonzero(creal z)
{
   return !isnan(z.im) && !isnan(z.re) && z != 0 + 0i;
}


Don Clugston wrote:
> Don Clugston wrote:
> 
> Actually, more accurate and simpler is to return +-NAN when presented with a NAN.
> This also satisfies the identity, because -NAN * NAN = -NAN.
> So actually there are six different possible return values from this function!
> It can be done with only a single branch.
> How about it, Walter?
> -----------------------------------
> 
> /** The mathematical sign() or signum() function.
>     Returns 1 if x is positive, -1 if x is negative,
>             0 if x is zero, real.nan if x is a nan.
>     Preserves sign of +-0 and +-nan.
> */
> real sign(real x)
> {
>    return x<>0 ? copysign(1, x) : x;
> }
> 
> unittest {
> assert( sign(-2.5) == -1 );
> assert( isNegZero(sign(-0.0)) );
> assert( isPosZero(sign(0.0)) );
> assert( sign(real.infinity)==1 );
> assert( isPosNan(sign(real.nan)) );
> assert( isNegNan(sign(-real.nan)) );
> }
> 
> // internal only
> 
> bit isPosNan(real x)
> {
>    return isnan(x) && !signbit(x);
> }
> 
> bit isNegNan(real x)
> {
>    return isnan(x) && signbit(x);
> }
> 
> // for completeness
> 
> int sign(int x)
> {
>   return x>0 ? 1 : x<0 ? -1 : x;
> }
> 
> long sign(long x)
> {
>   return x>0 ? 1 : x<0 ? -1 : x;
> }

"Don Clugston" <dac@nospam.com.au> wrote in message news:dg8jvj$a5e$1@digitaldaemon.com...
> behaves intuitively, and satisfies
> sign(x)*abs(x) is the same as x, for all x,
> even when x=-0.0, +-NAN.

NaN math follows its own logic, which is why most math functions will first check for NaN operands and do something special with them separately from the main logic.

The sign bit for a NaN has no meaning, and any program that relies on it is faulty. abs(NaN) should return a NaN, not +NaN, as there really is no such thing as a + or - NaN even though the bit is there. sign(NaN) should therefore return NaN.

This implies that floating point hardware may give random results for the NaN sign bit. What one must be careful of, in using things like signbit(), is that one does not generate a dependence on the sign of NaN.

(sign(NaN)*abs(NaN)==NaN) must return false, even if the bit patterns match
exactly, as any comparison op with either or both operands being NaN will
return false.

Walter Bright wrote:
> "Don Clugston" <dac@nospam.com.au> wrote in message
> news:dg8jvj$a5e$1@digitaldaemon.com...
> 
>>behaves intuitively, and satisfies
>>sign(x)*abs(x) is the same as x, for all x,
>>even when x=-0.0, +-NAN.
> 
> 
> NaN math follows its own logic, which is why most math functions will first
> check for NaN operands and do something special with them separately from
> the main logic.
> 
> The sign bit for a NaN has no meaning, and any program that relies on it is
> faulty. abs(NaN) should return a NaN, not +NaN, as there really is no such
> thing as a + or - NaN even though the bit is there. sign(NaN) should
> therefore return NaN.
> 
> This implies that floating point hardware may give random results for the
> NaN sign bit. What one must be careful of, in using things like signbit(),
> is that one does not generate a dependence on the sign of NaN.
> 
> (sign(NaN)*abs(NaN)==NaN) must return false, even if the bit patterns match
> exactly, as any comparison op with either or both operands being NaN will
> return false.

Right. I realized the problem when I considered the creal case.
If a creal z contains a nan and a non-nan, there is no way that
signum(z)*abs(z) can be the same as z, because abs(z) must be nan.

I was thinking there are some cases where a NaN is produced where the sign is still known,  (eg +infinity/+infinity) but since NaNs are not comparable it doesn't really make any sense.

Does D have a policy on dealing with NaNs? If a function recieves a NaN,
it is required to pass on the same NaN, or is it OK to simply return real.nan? (thereby destroying any information stored in the spare NaN bits, which it seems that D is not using). The former would cause issues if a function recieves two nans -- which should be passed on?

In this case the non-signed-ness of NaN doesn't make any difference to the code, just to the comments and the unit tests.
I'll do a version for complex numbers and resubmit, together with abs.
---------------------------

/** The mathematical sign (positive, negative, or zero)
    Defined such that signum(x)*abs(x) == x
    Returns: 1 if x is positive,
            -1 if x is negative,
             0 if x is zero.
    Preserves sign of +-0. Returns NaN if x is a NaN.
*/
real signum(real x)
{
   return x<>0 ? copysign(1, x) : x;
}

unittest {
assert( signum(-2.5) == -1 );
assert( isNegZero(signum(-0.0)) );
assert( isPosZero(signum(0.0)) );
assert( signum(real.infinity)==1 );
assert( isnan(signum(real.nan)) );
}

Don Clugston wrote:
> 
> Walter Bright wrote:
> 
>> "Don Clugston" <dac@nospam.com.au> wrote in message news:dg8jvj$a5e$1@digitaldaemon.com...
>> 
>>> behaves intuitively, and satisfies sign(x)*abs(x) is the same as
>>> x, for all x, even when x=-0.0, +-NAN.
>> 
>> NaN math follows its own logic, which is why most math functions
>> will first check for NaN operands and do something special with
>> them separately from the main logic.
>> 
>> The sign bit for a NaN has no meaning, and any program that relies
>> on it is faulty. abs(NaN) should return a NaN, not +NaN, as there
>> really is no such thing as a + or - NaN even though the bit is
>> there. sign(NaN) should therefore return NaN.
>> 
>> This implies that floating point hardware may give random results
>> for the NaN sign bit. What one must be careful of, in using things
>> like signbit(), is that one does not generate a dependence on the
>> sign of NaN.
>> 
>> (sign(NaN)*abs(NaN)==NaN) must return false, even if the bit
>> patterns match exactly, as any comparison op with either or both
>> operands being NaN will return false.


> I was thinking there are some cases where a NaN is produced where the
>  sign is still known,  (eg +infinity/+infinity) but since NaNs are
> not comparable it doesn't really make any sense.
> 
> Does D have a policy on dealing with NaNs? If a function recieves a
> NaN, it is required to pass on the same NaN, or is it OK to simply
> return real.nan? (thereby destroying any information stored in the
> spare NaN bits, which it seems that D is not using). The former would
> cause issues if a function recieves two nans -- which should be
> passed on?
> 
> In this case the non-signed-ness of NaN doesn't make any difference
> to the code, just to the comments and the unit tests. I'll do a
> version for complex numbers and resubmit, together with abs.

The standard leaves so many bit patterns for NaN that theoretically we could do serious math without ever resorting to other that NaN numbers. ;-)

However, looking at this from the other side, this just means that no bits in the entire NaN are significant, except for those that explicitly denote NaN.

NaN exists for the sole purpose of propagating "the improperness" of a calculation result up the chain to where it can be intelligently handled, somewhat similar to exceptions in D, Java, and other languages.

Therefore, the *Politically Correct* way to handle a NaN (as one or more of the arguments to the function) is to return the relevant compiler constant for NaN. In so doing the library writer explicitly acknowledges that we are propagating *nothing more* than the knowledge that this time we could not produce a numeric result -- thus invalidating also all dependent numeric calculations.

Since the standard does not specify anything about the other bits, it is reasonable to assume that different brands of math processors (and most probably even different versions from the same manufacturer) return different kinds of garbage in a NaN. Similarly, passing one or more NaN values to the math processor, and of course getting a NaN as a result, may preserve the original bit pattern(s), or then again it may not. (The fancy word for this is Undefined.)

> If a function recieves a NaN, it is required to pass on the same NaN?

Let me rephrase:

"If a function receives something that is not a number" ... "it is required to report this misconduct." And only that.

-----

One might of course argue that it is costlier to fetch and return a constant NaN when we already have the offending one ready and at hand. IMHO that is not relevant because that only makes a difference in a tight loop (or whatever piece of code that gets run monstrous times in a short while). I reality, if that is the case, then the programmer ought to rethink his program logic anyway.

"Don Clugston" <dac@nospam.com.au> wrote in message news:dgbkav$ruj$1@digitaldaemon.com...
> Does D have a policy on dealing with NaNs?

Yes. It follows the IEEE 754 spec and the general convention on them.

> If a function recieves a NaN,
> it is required to pass on the same NaN, or is it OK to simply return
> real.nan? (thereby destroying any information stored in the spare NaN
> bits, which it seems that D is not using).

The function should do its best to propagate the original NaN.

> The former would cause issues
> if a function recieves two nans -- which should be passed on?

Either.

> In this case the non-signed-ness of NaN doesn't make any difference to
> the code, just to the comments and the unit tests.
> I'll do a version for complex numbers and resubmit, together with abs.
> ---------------------------
>
> /** The mathematical sign (positive, negative, or zero)
>      Defined such that signum(x)*abs(x) == x
>      Returns: 1 if x is positive,
>              -1 if x is negative,
>               0 if x is zero.
>      Preserves sign of +-0. Returns NaN if x is a NaN.
> */
> real signum(real x)
> {
>     return x<>0 ? copysign(1, x) : x;
> }
>
> unittest {
> assert( signum(-2.5) == -1 );
> assert( isNegZero(signum(-0.0)) );
> assert( isPosZero(signum(0.0)) );
> assert( signum(real.infinity)==1 );
> assert( isnan(signum(real.nan)) );
> }

I don't understand why this function is necessary.

Walter Bright wrote:
> "Don Clugston" <dac@nospam.com.au> wrote in message
>>If a function recieves a NaN,
>>it is required to pass on the same NaN, or is it OK to simply return
>>real.nan? (thereby destroying any information stored in the spare NaN
>>bits, which it seems that D is not using).
> 
> The function should do its best to propagate the original NaN.
> 
>>The former would cause issues
>>if a function recieves two nans -- which should be passed on?
> 
> Either.

Great. That clarification is very helpful.

>>real signum(real x)
> 
> I don't understand why this function is necessary.

Nor do I, really. copysign() achieves the same thing in every case I've been able to think of. There seemed to be a lot of objection to removing it from std.math2, though. I figured that if it was included, it should be implemented correctly. If you've unmoved by the clamour, just throw it away.



It seems as though the only remaining significant _mathematical_ function in std.math2 is trunc() -- or is it already covered by floor() and ceil() ?

The other copyright question surrounds hypot(). In the cephes documentation, this is the only license info I've found:

"  Some software in this archive may be from the book _Methods and
Programs for Mathematical Functions_ (Prentice-Hall or Simon & Schuster
International, 1989) or from the Cephes Mathematical Library, a
commercial product. In either event, it is copyrighted by the author.
What you see here may be used freely but it comes with no support or
guarantee. " -- Cephes readme.

Does this mean a BSD-style license?
If so, then Cephes could also be used as a reference for some of the trickier math and statistic functions.  (Starting with lgamma and tgamma, which are declared in std.c.math, but which cause linker errors if actually used :-) ).
If not, then hypot() might be a problem.

Opinion?

On second thoughts, cis should be renamed to exp(ireal).
Start to take advantage of the inbuilt support for imaginary numbers!
BTW, I was delighted to find that the .im and .re properties are
supported for real and ireal, despite it not being mentioned anywhere
in the docs.
Corollary to The First Rule of D: It's always better than you think.

--------------------

//  exp(i theta) = cos(theta) + sin(theta)i.
version (X86)
{
    // On x86, this is almost as fast as sin(x).
    creal exp(ireal x)
    {
     asm {
         fld x;
         fsincos;
         fxch st(1), st(0);
      }
    }
} else {
    creal exp(ireal x)
    {
      return cos(x.im) + sin(x.im)*1i;
    }
}

unittest {
    assert(exp(0i)==1+0i);
    assert(exp(1.3e5i)==cos(1.3e5)+sin(1.3e5)*1i);
    creal c = exp(ireal.nan);
    assert(isnan(c.re) && isnan(c.im));
    c = exp(ireal.infinity);
    assert(isnan(c.re) && isnan(c.im));
}

"Don Clugston" <dac@nospam.com.au> wrote in message news:dgdq0d$30tl$1@digitaldaemon.com...
> Nor do I, really. copysign() achieves the same thing in every case I've been able to think of. There seemed to be a lot of objection to removing it from std.math2, though. I figured that if it was included, it should be implemented correctly. If you've unmoved by the clamour, just throw it away.

copysign() is a standard function and needs to be there. It just doesn't
need to be in std.math2.

> It seems as though the only remaining significant _mathematical_
> function in std.math2 is trunc() -- or is it already covered by floor()
> and ceil() ?

trunc() is in the DMC standard library. I just need to trannsfer it.

> The other copyright question surrounds hypot(). In the cephes documentation, this is the only license info I've found:
>
> "  Some software in this archive may be from the book _Methods and Programs for Mathematical Functions_ (Prentice-Hall or Simon & Schuster International, 1989) or from the Cephes Mathematical Library, a commercial product. In either event, it is copyrighted by the author. What you see here may be used freely but it comes with no support or guarantee. " -- Cephes readme.
>
> Does this mean a BSD-style license?
> If so, then Cephes could also be used as a reference for some of the
> trickier math and statistic functions.  (Starting with lgamma and
> tgamma, which are declared in std.c.math, but which cause linker errors
> if actually used :-) ).
> If not, then hypot() might be a problem.
>
> Opinion?

I think the Cephes libraries are the best out there. A couple years ago, I emailed Steve Moshier, the author, about the license with the same question, and here was his response:

---------------------------------------------------------------------
Thank you for writing. I have not put any restrictions on the use of
material posted to the net. Bell Labs has suggested this for the files
in netlib/cephes:
http://www.netlib.org/research/boilerplate
---------------------------------------------------------------------

Since then, I have filled in many gaps in the DMC runtime library with Steve's Cephes code. The Cephes code seems based on the old Cody&Waite algorithms, which are based on the earlier Hart&Cheney book.

Cephes does have functions missing from DMC like asinh, and the only reason I hadn't incorporated them yet is because it's a fair amount of work to adapt them, and I'm lazy. The thicket of #ifdef's need to be removed, the symbols need to be converted to DMC symbols, the tables need to be converted from hex kludges to the modern hex float format, the poly calculation needs to be rewritten to use DMC's poly function, and the special cases need to be rewritten to use the x87's features.

There's an interesting article on slashdot today about a new method for solving trig. problems without using sin/cos/tan.  A link to info on the book is here:

http://web.maths.unsw.edu.au.nyud.net:8090/~norman/book.htm

I'll be curious to see the difference in speed and complexity of this new method compared to the classic approach.  I don't suppose anyone has heard of this guy before?


Sean

In article <dghg6n$h8d$1@digitaldaemon.com>, Sean Kelly says...
>
>There's an interesting article on slashdot today about a new method for solving trig. problems without using sin/cos/tan.  A link to info on the book is here:
>
>http://web.maths.unsw.edu.au.nyud.net:8090/~norman/book.htm
>
>I'll be curious to see the difference in speed and complexity of this new method compared to the classic approach.  I don't suppose anyone has heard of this guy before?

Upon further reading, while the author talks up his ideas quite a lot, I don't see much that's actually revolutionary.  It may turn out that the best aspect of this approach is its use as a teaching aid.


Sean