Here are the three least trivial mathematical functions from math2:
cosh, sinh, tanh.
Those functions in math2 had copyright issues, and also had several bugs. The implementations here are more extensively tested, faster, and are public domain.
I also include cis(x) = cos + i sin, which is a fast operation on
x86 CPUs.

Also, productLog() and productLog_1().
This is the inverse of x * exp(x).
I've used the simple Mathematica names rather than more common (but obscure) names, LambertW() and LambertW_1(). This function crops up all over the place but usually goes unrecognised because it was unnamed until quite recently (despite being first described in the 1700s!)
If interested, look it up in Wikipedia.

There's also an internal function for unit testing, which builds on my
feqrel() function.
This is a "second order" assert in that it allows a kind of 'for all'
operation for testing inverses. This function is not intended to be publicly documented, but should be useful for other functions in std.math.

Hopefully this moves us a lot closer to deprecating std.math2.

Note that several of the other trig functions in std.math2 are just plain wrong.
acot(), asec() and acosec() return incorrect results.
They should be deleted. If anyone is using them, their code is wrong!

Also, sign(real x) ignores the sign of negative zero,
so that it violates the identity

sign(x)*abs(x) == x   when x=-0.0.

It should never return 0, since zero has a sign in IEEE.
It should be

real sign(real x)
{
  return signbit(x) ? -1.0 : 1.0;
}

and therefore ints should never return 0 either.

int sign(int x)
{
   return x<0 ? -1 : 1;
}

long sign(long x)
{
   return x<0 ? -1 : 1;
}

-Don.

//------------------------------------------------------------------
// (Private) functions for unit tests

/**
 Test the consistency of a real function which has an inverse.
 Returns the worst (minimum) value of feqrel(x, inversefunc(forwardfunc(x)))
 for all x in the domain.
 Author: Don Clugston. License: Public Domain
 Params:
        forwardfunc   Function to be tested
        inversefunc   Inverse of forwardfunc
        domain        A sequence of pairs of numbers giving the first and last points which are valid for the function.
eg. (-real.infinity, real.infinity) ==> valid everywhere
    (-real.infinity, -real.min, real.min, real.infinity) ==> valid for x!=0.
 Returns:
    number of bits for which x and inversefunc(forwardfunc(x)) are equal for
    every point x in the domain.
    -1 = at least one point is wrong by a factor of 2 or more.
*/
int consistencyRealInverse(real function (real) forwardfunc,
                    real function (real) inversefunc, real [] domain ...)
{
      assert(domain.length>=2); // must have at least one valid range
      assert((domain.length & 1) == 0); // must have even number of endpoints.

      int worsterror=real.mant_dig+1;
      real worstx=domain[0]; // where does the biggest discrepancy occur?

      void testPoint(real x) {
         for (int i=0; i<domain.length; i+=2) {
            if (x>=domain[i] && x<=domain[i+1]) {
				int u=feqrel(x, inversefunc(forwardfunc(x)) );
				if (u<worsterror) { worsterror=u;  worstx=x; }
				return;
		     }
		 }
      }
      // test the edges of the domains
      foreach(real y; domain) testPoint(y);
      real x = 1.01;
	  // first, go from 1 to infinity
	  for (x=1.01; x!=x.infinity; x*=2.83) testPoint(x);
	  // then from 1 to +0
	  for (x=0.98; x!=0; x*=0.401) testPoint(x);
	  // from -1 to -0
	  for (x=-0.93; x!=0; x*=0.403) testPoint(x);
	  // from -1 to -infinity
	  for (x=-1.09; x!=-x.infinity; x*=2.97) testPoint(x);
//	   writefln("Worst error is ", worsterror, " at x= ", worstx);
	  return worsterror;
}

// return true if x is -0.0
bit isNegZero(real x)
{
   return (x==0) && (signbit(x)==1);
}

// return true if x is +0.0
bit isPosZero(real x)
{
   return (x==0) && (signbit(x)==0);
}


//------------------------------------------------------------------

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

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

//------------------------------------------------------------------
//        Inverse hyperbolic trig functions

/** Inverse hyperbolic sine - the inverse of sinh(x)
   Author: Don Clugston   License: Public Domain.
   Definition:
   asinh(x) =  log( x + sqrt( x*x + 1 )) if x >= +0
   asinh(x) = -log(-x + sqrt( x*x + 1 )) if x <= -0
   Preserves sign of -0.0
   Domain: -infinity..+infinity
   Range: -infinity, -log(real.max)..log(real.max), +infinity
   Special values:
   -inf   -inf
   -0     -0
   +0     +0
   +inf   +inf
*/
real asinh(real x)
{
 	if (fabs(x) > 1/real.epsilon)   // beyond this point, x*x + 1 == x*x
 			return copysign(LN2 + log(fabs(x)), x);
	else {
	   // sqrt(x*x + 1) ==  1 +  x * x / ( 1 + sqrt(x*x + 1) )
       return copysign(log1p(fabs(x) + x*x / (1 + sqrt(x*x + 1)) ), x);
	}
}

unittest {
 assert(isPosZero(asinh(0.0)));
 assert(isNegZero(asinh(-0.0)));
 assert(asinh(real.infinity)==real.infinity);
 assert(asinh(-real.infinity)==-real.infinity);
 assert(isnan(asinh(real.nan)));
 version (X86) {
	 assert( consistencyRealInverse(&asinh, &sinh,-double.max, double.max) >= double.mant_dig);
	 assert( consistencyRealInverse(&asinh, &sinh,-1, 1) >= real.mant_dig-2);
	 assert( consistencyRealInverse(&sinh, &asinh,
			 -log(real.max)*(1+real.epsilon), log(real.max)*(1-real.epsilon)
			 )>= double.mant_dig);
 }
}

/** Inverse hyperbolic cosine - the inverse of cosh(x)
   Author: Don Clugston   License: Public Domain.
   Definition:
     acosh(x) = log(x + sqrt( x*x -1))
   Domain:   1..infinity
   Range:    1..log(real.max), infinity
   Special values:
   <1     nan
   1      0
   +inf   +inf
*/
real acosh(real x)
{
	if (x > 1/real.epsilon) return LN2 + log(x);
	else return log(x + sqrt(x*x - 1));
}

unittest {
  assert(isnan(acosh(0.9)));
  assert(isnan(acosh(real.nan)));
  assert(acosh(1)==0.0);
  assert(acosh(real.infinity) == real.infinity);
    version (X86) {
  assert( consistencyRealInverse(&acosh, &cosh, 1, double.max)
	         >= double.mant_dig);
  assert( consistencyRealInverse(&cosh, &acosh, 1, log(real.max)*(1-real.epsilon))
			 >= real.mant_dig-1);
    }
}

/** Inverse hyperbolic tangent - the inverse of tanh(x)
   Author: Don Clugston   License: Public Domain.
   Definition:
	 atanh(x) = log( (1+x)/(1-x) ) / 2
   Domain:   -infinity..infinity
   Range:    -1..1
*/
real atanh(real x)
{
	// log( (1+x)/(1-x) ) == log ( 1 + (2*x)/(1-x) )
	return	 0.5 * log1p( 2 * x / (1 - x) );
}

unittest {
 assert(isPosZero(atanh(0.0)));
 assert(isNegZero(atanh(-0.0)));
 assert(isnan(atanh(real.nan)));
 assert(isNegZero(atanh(-real.infinity)));
    version (X86) {
 assert( consistencyRealInverse(&atanh, &tanh, -1, 1)
         >= real.mant_dig-2);
 assert( consistencyRealInverse(&tanh, &atanh, -1,1)
		 >= real.mant_dig-1);
    }
}

//------------------------------------------------------------------
//     Lambert's W function, also known as ProductLog

/**
   productLog() is the inverse of x*exp(x)
   and is also known as Lambert's W-function.
   productLog(x*exp(x)) == x, where x >= -exp(-1)
   A second inverse, productLog_1(x), can be defined for real x<=-exp(-1).

   Author: Don Clugston. Date: 2005-09-01. License: Public Domain.
   Using an algorithm originally developed by Keith Briggs
   http://keithbriggs.info
*/
real productLog(real x)
{
    real w;
	if (x < 0.5) {
		if (fabs(x) < x.epsilon/2) return x; // exp(x)==1 for small x.
		if (x <= -1/E) {
			if (x== -1/E) return -1; // avoid div by zero
			return x.nan; // no solution
		}
		real p = sqrt( 2.0 * ( E*x + 1.0 ) );
        w=-1.0+p*(1.0+p*(-1/3.0 + p*11.0/72.0));
	} else {
	   if (isnan(x)) return x.nan;
		w = log(x);
		if (x > 3) { w = w-log(w); if (x==x.infinity) return x.infinity; }
	}
	real e,t,w2=w;
	
	int kount=0;
	do { // Halley iterations. Never requires more than 4 iterations.
	       w = w2;
		   e = exp(w);
		   t = w*e - x;
		   t /= e*(w+1.0) - 0.5*(w+2.0)*t / (w+1.0);
		   w2 -= t;
		   kount++;
		   if (kount>20) { writefln(x); return x; }
	} while (fabs(t)>=x.epsilon*(1.0+fabs(w2)));
	return w2;
}

/**
   The inverse of y*exp(y) where y <= -exp(-1)
   Params:
    x     -1/E <= x <=  0
   Returns:
    y such that y*exp(y)==x and y <= -exp(-1).

   Author: Don Clugston. Date: 2005-09-01. License: Public Domain.
   Using an algorithm originally developed by Keith Briggs
   (http://keithbriggs.info)
*/
real productLog_1(real x)
{
  real w;
  if (x<-1/E || x>=0.0 || isnan(x)) return x.nan;
  // initial approx for iteration...
  if (x<-1e-6) { // series about -1/e
    if (x==-1/E) return -1; // avoid div by zero
    real p=-sqrt(2.0*(E*x+1.0));
    w=-1.0+p*(1.0+p*(-1/3.0 + p*11.0/72.0));
  } else { // asymptotic near zero
    real l1=log(-x);
    real l2=log(-l1);
    w=l1-l2+l2/l1;
  }
  if (x>-real.epsilon) return -x.infinity;
  real e, t, w2=w;
	do { // Halley iterations. Never requires more than 10 iterations.
	       w = w2;
		   e = exp(w);
		   t = w*e - x;
		   t /= e*(w+1.0) - 0.5*(w+2.0)*t / (w+1.0);
		   w2 -= t;
	} while (fabs(t)>=x.epsilon*(1.0+fabs(w2)));
	return w2;
}

// x*exp(x). Used only for testing
real productExp(real x)
{
 return x*exp(x);
}

unittest {
 assert(productLog(-1/E)==-1);
 assert(isnan(productLog(real.nan)));
 assert(isNegZero(productLog(-0.0)));
 assert(isPosZero(productLog(0.0)));
 assert(productLog(E)==1);
    version (X86) {
 assert( consistencyRealInverse(&productLog, &productExp,
		-1/E, double.max) >= double.mant_dig);
 assert( consistencyRealInverse(&productExp, &productLog,-1, 11000)
		>= real.mant_dig-3);
    }

  assert(productLog_1(-1/E)==-1);
  assert(productLog_1(-real.infinity)==0);
  assert(isnan(productLog_1(real.nan)));
    version (X86) {
 assert( consistencyRealInverse(&productLog_1, &productExp,
		-1/E, -real.epsilon) >= double.mant_dig);
 assert( consistencyRealInverse(&productExp, &productLog_1,
	log(real.epsilon), -1) >= double.mant_dig);
    }
}

A couple more comments:

* The following functions are defined in std.c.math, but if you use them, you'll get a link time error:

acoshl(), asinhl(), atanhl().

They should be removed from std.c.math and the (commented out) forwarding functions in std.math should also be removed.

* The list of constants in std.math should include

const ireal I = 1.0i;

for compatibility with C99 and for situations like

cos(PI) + I * sin(PI)

where you want to convert a function result from real to imaginary.

-Don.

In article <dg41m1$1o5k$1@digitaldaemon.com>, Don Clugston says...
>
>A couple more comments:
>
>* The following functions are defined in std.c.math, but if you use them, you'll get a link time error:
>
>acoshl(), asinhl(), atanhl().

With DMD?  That's odd.  dm/include/math.h has declarations for these functions.


Sean

Sean Kelly wrote:
>>* The following functions are defined in std.c.math, but if you use them, you'll get a link time error:
>>
>>acoshl(), asinhl(), atanhl().
> 
> 
> With DMD?  That's odd.  dm/include/math.h has declarations for these functions.
> 
> 
> Sean

Yes, it is odd. Test program:
--------------------------------
import std.c.math;

int main()
{
  real x = acoshl(1.0);
  return 0;
}
-------------------------------
c:\dmd\bin\..\..\dm\bin\link.exe bug,bug.exe,,user32+kernel32,bug.def/noi;

OPTLINK (R) for Win32  Release 7.50B1
Copyright (C) Digital Mars 1989 - 2001  All Rights Reserved

bug.obj(bug)
 Error 42: Symbol Undefined _acoshl
--- errorlevel 1
--------------------------------------
Why? Well, they are not standard C routines. Not included in VC++,
for example. So they probably aren't part of gcc, hence they must
be defined for the sake of gdc. They were probably removed from the dmd lib when gdc was born.
I suspect that the functions in dmc are copyright and so can't be
used in phobos.
Just a guess.
- Don.

In article <dg49vr$1vch$1@digitaldaemon.com>, Don Clugston says...
>
>Sean Kelly wrote:
>>>* The following functions are defined in std.c.math, but if you use them, you'll get a link time error:
>>>
>>>acoshl(), asinhl(), atanhl().
>> 
>> 
>> With DMD?  That's odd.  dm/include/math.h has declarations for these functions.
>> 
>> 
>> Sean
>
>Yes, it is odd. Test program:
>--------------------------------
>import std.c.math;
>
>int main()
>{
>   real x = acoshl(1.0);
>   return 0;
>}
>-------------------------------

What about this:

real x = acoshl( cast(real) 1.0 );

>Why? Well, they are not standard C routines.

They are standard C routines as of the C99 spec.  But DMC is more up to date than most compilers with respect to C99 library conformance.  For what it's worth, I have a full implementation of std.c.math here:

http://svn.dsource.org/projects/ares/trunk/src/ares/std/c/math.d

It's correct for DMD as far as I know, but other compilers may not implement all the declared functions.


Sean

"Don Clugston" <dac@nospam.com.au> wrote in message news:dg41m1$1o5k$1@digitaldaemon.com...
> * The following functions are defined in std.c.math, but if you use them, you'll get a link time error:
>
> acoshl(), asinhl(), atanhl().

That's because they aren't implemented in the C runtime library yet. :-(

> * The list of constants in std.math should include
>
> const ireal I = 1.0i;
>
> for compatibility with C99 and for situations like

I disagree with forwarding I from C into D, it's not needed for compatibility.

> cos(PI) + I * sin(PI)
>
> where you want to convert a function result from real to imaginary.

cos(PI) + sin(PI) *1i;

"Don Clugston" <dac@nospam.com.au> wrote in message news:dg40b2$1mud$1@digitaldaemon.com...
>
> Here are the three least trivial mathematical functions from math2:
> cosh, sinh, tanh.
> Those functions in math2 had copyright issues, and also had several
> bugs. The implementations here are more extensively tested, faster, and
> are public domain.
> I also include cis(x) = cos + i sin, which is a fast operation on
> x86 CPUs.

Good!

> Also, productLog() and productLog_1().
> This is the inverse of x * exp(x).
> I've used the simple Mathematica names rather than more common (but
> obscure) names, LambertW() and LambertW_1(). This function crops up all
> over the place but usually goes unrecognised because it was unnamed
> until quite recently (despite being first described in the 1700s!)
> If interested, look it up in Wikipedia.
>
> There's also an internal function for unit testing, which builds on my
> feqrel() function.
> This is a "second order" assert in that it allows a kind of 'for all'
> operation for testing inverses. This function is not intended to be
> publicly documented, but should be useful for other functions in std.math.

Ok.

> Hopefully this moves us a lot closer to deprecating std.math2.

I want to get rid of std.math2.

> Note that several of the other trig functions in std.math2 are just
> plain wrong.
> acot(), asec() and acosec() return incorrect results.
> They should be deleted. If anyone is using them, their code is wrong!

Let's just get rid of std.math2.

> Also, sign(real x) ignores the sign of negative zero,
> so that it violates the identity
>
> sign(x)*abs(x) == x   when x=-0.0.
>
> It should never return 0, since zero has a sign in IEEE.
> It should be
>
> real sign(real x)
> {
>    return signbit(x) ? -1.0 : 1.0;
> }
>
> and therefore ints should never return 0 either.
>
> int sign(int x)
> {
>     return x<0 ? -1 : 1;
> }
>
> long sign(long x)
> {
>     return x<0 ? -1 : 1;
> }

sign() is a redundant function, it's more trouble to document and look up
than it's worth. signbit() is adequate for the job.

On Mon, 12 Sep 2005 13:14:23 -0700, Walter Bright wrote:


[snip]

>> Also, sign(real x) ignores the sign of negative zero,
>> so that it violates the identity
>>
>> sign(x)*abs(x) == x   when x=-0.0.
>>
>> It should never return 0, since zero has a sign in IEEE.
>> It should be
>>
>> real sign(real x)
>> {
>>    return signbit(x) ? -1.0 : 1.0;
>> }
>>
>> and therefore ints should never return 0 either.
>>
>> int sign(int x)
>> {
>>     return x<0 ? -1 : 1;
>> }
>>
>> long sign(long x)
>> {
>>     return x<0 ? -1 : 1;
>> }
> 
> sign() is a redundant function, it's more trouble to document and look up
> than it's worth. signbit() is adequate for the job.

Garbage! Its called 'abstraction'. The term 'signbit' implies that we are asking for a detail of the implementation rather than the mathematical concept of 'sign'.

You could just do an 'alias sign signbit;' if you really didn't want to do a better implementation. However, there isn't a signbit() function for ints and other number types.

-- 
Derek Parnell
Melbourne, Australia
13/09/2005 7:24:32 AM

"Derek Parnell" <derek@psych.ward> wrote in message news:y1h8ymdxkhwn.tg27eopfhs2k.dlg@40tude.net...
> Garbage! Its called 'abstraction'. The term 'signbit' implies that we are asking for a detail of the implementation rather than the mathematical concept of 'sign'.

One always must be careful when thinking about the "mathematical" concept of sign, because floating point doesn't work that way. There's the issue of negative 0, which is not really a mathematical concept, it's more of a kludge to deal with the finite precision of the floating point format. There's also the sign bit of a NaN number, what does that mean?

Any time one wants to test for the signedness of a floating point number, one has to think about "what happens if it is -0" and "what happens if it is a NaN". Making a function called "sign()" can't resolve the problem. It'll pretend to, and that'll lead to bugs. If one wants to go beyond (a<0), using signbit() is the right approach, because it causes one to have to think about the other cases.

The signbit() function is necessary because otherwise getting at the sign bit of the floating point value is a clumsy, nonportable hack. A sign() function does not add sufficient value, and it may even serve to obfuscate things because the reader will have to go look up the doc on it to see what it does with -0 and -NaN. std.math should strive to be an orthogonal set of nontrivial building blocks. Don's work with the arc hyperbolics are a perfect example of great additions to it.

On Mon, 12 Sep 2005 16:59:23 -0700, Walter Bright wrote:

> "Derek Parnell" <derek@psych.ward> wrote in message news:y1h8ymdxkhwn.tg27eopfhs2k.dlg@40tude.net...
>> Garbage! Its called 'abstraction'. The term 'signbit' implies that we are asking for a detail of the implementation rather than the mathematical concept of 'sign'.
> 
> One always must be careful when thinking about the "mathematical" concept of sign, because floating point doesn't work that way. There's the issue of negative 0, which is not really a mathematical concept, it's more of a kludge to deal with the finite precision of the floating point format. There's also the sign bit of a NaN number, what does that mean?

Huh? Okay, it must be me. I'm a simple person with simple ideas.

When somebody asks me what's the sign of such-and-such, I just say "Negative" if its less than zero, "Positive" if greater than zero and "Nothing" if it is zero OR not a number.

What has that got to do with floating point, conceptually or physically?

> Any time one wants to test for the signedness of a floating point number, one has to think about "what happens if it is -0"

is -0 less than zero? If yes, then its negative.

> and "what happens if it is a NaN".

NaN is Not A Number so therefore it doesn't have a sign.

>Making a function called "sign()" can't resolve the problem.

I think it can ...

  enum signedness
  {
    Negative = -1,
    Nothing = 0,
    Positive = 1
  }

  template getsign(TYPE) {
  signedness getsign( TYPE x)
  {
     if (x < 0)
        return signedness.Negative;
     if (x > 0)
        return signedness.Positive;
     return signedness.Nothing;
  }
  }
 alias getsign!(float)  sign;
 alias getsign!(double) sign;
 alias getsign!(real)   sign;
 alias getsign!(int)    sign;
 alias getsign!(long)   sign;
 alias getsign!(short)  sign;




>It'll
> pretend to, and that'll lead to bugs. If one wants to go beyond (a<0), using
> signbit() is the right approach, because it causes one to have to think
> about the other cases.

So does using enums.

> The signbit() function is necessary because otherwise getting at the sign bit of the floating point value is a clumsy, nonportable hack.

Yes, signbit() is required to GET AT THE SIGN BIT. But that is not what I'm talking about. That is an implementation detail of the particular floating point representation system you choose to use. It has nothing to do with the simple mathematical question "what is the sign of X?", which may be a floating point, fixed point, or integer.

> A sign() function does not add sufficient value,

Says who? You? And how come you get to be the Grand Arbiter of "sufficient value"? What empirical measurement did you apply to get to that decision? Define 'sufficient' in this context. I would say that value is in the eye of the beholder.

>and it may even serve to obfuscate
> things because the reader will have to go look up the doc on it to see what it does with -0 and -NaN. std.math should strive to be an orthogonal set of nontrivial building blocks.

Why 'nontrivial'? What may be trivial to one person is not so to others. Some people regard various sort algorithms as trivial. The library should be an enabler.

-- 
Derek
(skype: derek.j.parnell)
Melbourne, Australia
13/09/2005 2:19:42 PM