/* mathx.c (emx+gcc) -  additional math function eXtensions.
   Parts: Cephes Math Library Release 2.3:  March, 1995
Copyright 1984, 1995 by Stephen L. Moshier
   Documentation is included on the distribution media as
   Unix-style manual pages that describe the functions and their
   invocation.  The primary documentation for the library functions
   is the book by Moshier, Methods and Programs for Mathematical Functions, 
   Prentice-Hall, 1989.
 */

#include <math.h>
#include <mathx.h>
#include "math_private.h"

int signgam = 0; /* global */

#undef MATH_INL
#define MATH_INL extern /* these are external library functions */

/*							drand.c
 *
 *	Pseudorandom number generator
 *
 *
 *
 * SYNOPSIS:
 *
 * double y, drand();
 *
 * drand( &y );
 *
 *
 *
 * DESCRIPTION:
 *
 * Yields a random number 1.0 <= y < 2.0.
 *
 * The three-generator congruential algorithm by Brian
 * Wichmann and David Hill (BYTE magazine, March, 1987,
 * pp 127-8) is used. The period, given by them, is
 * 6953607871644.
 *
 * Versions invoked by the different arithmetic compile
 * time options DEC, IBMPC, and MIEEE, produce
 * approximately the same sequences, differing only in the
 * least significant bits of the numbers. The UNK option
 * implements the algorithm as recommended in the BYTE
 * article.  It may be used on all computers. However,
 * the low order bits of a double precision number may
 * not be adequately random, and may vary due to arithmetic
 * implementation details on different computers.
 *
 * The other compile options generate an additional random
 * integer that overwrites the low order bits of the double
 * precision number.  This reduces the period by a factor of
 * two but tends to overcome the problems mentioned.
 *
 */

/*  Three-generator random number algorithm
 * of Brian Wichmann and David Hill
 * BYTE magazine, March, 1987 pp 127-8
 *
 * The period, given by them, is (p-1)(q-1)(r-1)/4 = 6.95e12.
 */

static int sx = 1;
static int sy = 10000;
static int sz = 3000;

static union {
 double d;
 unsigned short s[4];
} unkans;

/* This function implements the three
 * congruential generators.
 */
 
	static int 
ranwh(void)
{
int r, s;
/*  sx = sx * 171 mod 30269 */
r = sx/177;
s = sx - 177 * r;
sx = 171 * s - 2 * r;
if( sx < 0 )
	sx += 30269;


/* sy = sy * 172 mod 30307 */
r = sy/176;
s = sy - 176 * r;
sy = 172 * s - 35 * r;
if( sy < 0 )
	sy += 30307;

/* sz = 170 * sz mod 30323 */
r = sz/178;
s = sz - 178 * r;
sz = 170 * s - 63 * r;
if( sz < 0 )
	sz += 30323;
/* The results are in static sx, sy, sz. */
return 0;
}

/*	drand.c
 *
 * Random double precision floating point number between 1 and 2.
 *
 * C callable:
 *	drand( &x );
 */

int 
drand(double *a)
{
unsigned short r;

/* This algorithm of Wichmann and Hill computes a floating point
 * result:
 */
ranwh();
unkans.d = sx/30269.0  +  sy/30307.0  +  sz/30323.0;
r = unkans.d;
unkans.d -= r;
unkans.d += 1.0;

/* make a random 16 bit integer
 * to overwrite the least significant word
 * of unkans.
 */
ranwh();
r = sx * sy + sz;
unkans.s[0] = r;
*a = unkans.d;
return 0;
}

  void
pvec(double x)
{
  union
  {
    double d;
    unsigned short s[4];
  } u;
  int i;

  u.d = x;
  for (i = 0; i < 4; i++)
    printf ("0x%04x ", u.s[i]);
  printf ("\n");
}


/* Parts: 
   Copyright (C) 1991, 92, 93, 95, 96, 97, 98 Free Software 
   Foundation, Inc. This file is part of the GNU C Library.

   The GNU C Library is free software; you can redistribute it and/or
   modify it under the terms of the GNU Library General Public License as
   published by the Free Software Foundation; either version 2 of the
   License, or (at your option) any later version.

   The GNU C Library is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
   Library General Public License for more details.

   You should have received a copy of the GNU Library General Public
   License along with the GNU C Library; see the file COPYING.LIB.  If not,
   write to the Free Software Foundation, Inc., 59 Temple Place - Suite 
   330, Boston, MA 02111-1307, USA.
 *
 *	Extensions to ISO C Standard: 4.5 MATHEMATICS	<math.h>
 */

/* BSD Functions */
		MATH_INL long double
asinhl(long double x) {
 register __const__ long double y = fabsl(x);
	if (__isnanl(x)) return x;
	if (x == 0.0L) return x;
	if (!(__isfinitel(x))) return x;
	if (y > 1.0e10L) return logl(x) + _LOGE2L;
#ifdef GNUFNS  
	return log1pl( (y * y /(sqrtl(y * y + 1.0L) + 1.0L) + y) )
			* _sgn1l(x); 	/* This is corrected */
#else
	return logl( x + sqrtl(x * x + 1.0L));
#endif
}
	MATH_INL long double
acoshl(long double x) {
	if (__isnanl(x)) return x;
	if (x > 1.0e10L) {
		if (x == _INFINITYL) return _INFINITYL; 
		return logl(x) + _LOGE2L; }
        return logl(x + sqrtl(x - 1.0L) * sqrtl(x + 1.0L));
}
	MATH_INL long double
atanhl(long double x) {
  register __const__ long double y = fabsl(x);
	if (x == 0.0L) return x;
	if (y >= 1.0L) {
		if (x == 1.0L) return _INFINITYL;
		if (x == -1.0L) return(- _INFINITYL);
	}
	if (y < 1.0e-8L) return x;
#ifdef GNUFNS
	return -0.5L * log1pl(-(y + y)/(1.0L + y)) * _sgn1(x);
#else
	return  0.5L *   logl( (1.0L + x)/(1.0L - x) );
#endif
}
	MATH_INL double
coshm1(double x) {
  register __const__ double __exm1 = expm1 (fabs (x));
  return 0.5 * (__exm1 / (__exm1 + 1.0)) * __exm1;
}
	MATH_INL double
acosh1p(double x) {
  return log1p(x + sqrt (x) * sqrt (x + 2.0));
}
	MATH_INL void
sincos(double x, double *__sinx, double *__cosx) {
 register double __cosr, __sinr;
  __asm__ __volatile__ ("fsincos"
     : "=t" (__cosr), "=u" (__sinr) : "0" (x));
  *__sinx = __sinr;
  *__cosx = __cosr;
}
	MATH_INL double
pow2(double x) {
  register double value, exponent;
  int p = (int)x;

  if (x ==(double)p)
    return ldexp (1.0, p);

  __asm__ __volatile__
    ("fldl	%%st(0)\n\t"
     "frndint			# int(x)\n\t"
     "fxch\n\t"
     "fsub	%%st(1)		# fract(x)\n\t"
     "f2xm1			# 2^(fract(x)) - 1\n\t"
     : "=t" (value), "=u" (exponent) : "0" (x));
  value += 1.0;
  __asm__ __volatile__
    ("fscale"
     : "=t" (value) : "0" (value), "u" (exponent));

  return value;
}

