richb@sunaus.oz (Rich Burridge) (02/02/90)
**IMPORTANT NOTE** This patch introduces a (hopefully) portable
maths library. Please make sure you read the
details listed below carefully.
This is to try to make calctool portable to all
Unix systems. If you have to make any changes to
get it working on your machine, please let me
know, and I'll adjust the official source code
accordingly.
This patch consists of three parts. As well as patch #6, there are
two new files; mathlib.c and mathlib.h.
---- Cut Here and unpack ----
#!/bin/sh
# shar: Shell Archiver (v1.22)
# Packed Fri Feb 2 13:11:50 EST 1990 by stard!richb
# from directory /home/tmp
#
# This is part 1 of a multipart archive
# do not concatenate these parts, unpack them in order with /bin/sh
#
# Run the following text with /bin/sh to create:
# patch.6
# mathlib.h
# mathlib.c
#
if test -r s2_seq_.tmp
then echo "Must unpack archives in sequence!"
next=`cat s2_seq_.tmp`; echo "Please unpack part $next next"
exit 1; fi
echo "x - extracting patch.6 (Text)"
sed 's/^X//' << 'SHAR_EOF' > patch.6 &&
XThis is official patch #6 for calctool v2.4; please apply it.
X
X**IMPORTANT NOTE** This patch introduces a (hopefully) portable
X maths library. Please make sure you read the
X details listed below carefully.
X
X This is to try to make calctool portable to all
X Unix systems. If you have to make any changes to
X get it working on your machine, please let me
X know, and I'll adjust the official source code
X accordingly.
X
XIt makes the following changes.
X
X 1/ The Makefile did not include the CHANGES file in the OTHERS
X macro definition.
X
X 2/ If you are not using IEEE floating point format, and you need to
X use some of the routines from the portable maths library which is
X now included with this distribution, then you need to uncomment
X the new NOTIEEE definition in the Makefile.
X
X 3/ From Tom Friedel <tpf@jdyx.atlanta.ga.us>
X In the tty driver, in the drawtext routine, there is a character
X array of five elements. There was an invalid assignment to key[5].
X
X 4/ From Stephen Frede <stephenf@softway.oz>
X The definitions in the Makefile for HELPNAME, NEWSFILE and RCNAME
X do not need extra ...GIVEN definitions.
X
X 5/ Calctool.h defined NEWSNAME, when it should have defined NEWSFILE.
X
X 6/ With the X11 driver, the function XIconifyWindow is new to R4.
X If you are running an X11 release which is not R4, then you need
X to uncomment the new definition in the Makefile called NOT_R4.
X This nullifys the function of the calctool OFF key, but hey, you
X can also you the window manager to iconify calctool (assuming
X your window manager can do that!).
X
X 7/ The factorial function has been adjusted to only work with
X positive integers. I couldn't find a public domain or freely
X distributable lgamma function written in C.
X
X 8/ fillbox in graphics.c no longer needs to be passed a window
X parameter.
X
X 9/ From Stephen Frede <stephenf@softway.oz>
X Not all linkers have the -L option. The CALCLIB definition in the
X Makefile has been adjusted, and moved to the section where
X LIBNAME is defined. By default, you shouldn't have to do anything,
X but if you are on a Sun and want to use the shared library facility,
X then you need to uncomment three definitions.
X
X 10/ From Stephen Frede <stephenf@softway.oz>
X getcwd is not present on all systems. I've replaced it with getwd.
X
X
XIt add the following files:
X
X mathlib.h
X mathlib.c - These are an attempt to provide a portable maths
X library for those systems that don't have all the
X functions that calctool needs.
X
X You should look near the end of mathlib.h for a set
X of definitions, then set appropriately.
X
XUse Larry Walls patch program to apply these changes, then recompile.
X
X------CUT HERE------patch.6------CUT HERE------
X
X------- calctool.c -------
X*** /tmp/da5214 Fri Feb 2 13:01:42 1990
X--- calctool.c Tue Jan 30 12:08:45 1990
X***************
X*** 17,22 ****
X--- 17,25 ----
X * reported to me then an attempt will be made to fix them.
X */
X
X+ #include <stdio.h>
X+ #include <strings.h>
X+ #include <math.h>
X #include "patchlevel.h"
X #include "color.h"
X #include "calctool.h"
X***************
X*** 361,372 ****
X }
X
X
X! matherr(x) /* Calctools' math library error-handling routine. */
X! struct exception *x ;
X {
X! SPRINTF(display, "Error in %s", x->name) ;
X set_item(DISPLAYITEM, display) ;
X error = 1 ;
X set_item(OPITEM, "CLR") ;
X- return(1) ;
X }
X--- 364,390 ----
X }
X
X
X! /* Calctools' customised math library error-handling routine. */
X!
X! doerr(funcname, etype, eval)
X! char *funcname, *etype ;
X! int eval ;
X {
X! SPRINTF(display, "%s: %s error", funcname, etype) ;
X set_item(DISPLAYITEM, display) ;
X+ error = eval ;
X+ set_item(OPITEM, "CLR") ;
X+ }
X+
X+
X+ /* Default math library exception handling routine. */
X+
X+ /*ARGSUSED*/
X+ matherr(exc)
X+ struct exception *exc ;
X+ {
X+ STRCPY(display, "Error") ;
X+ set_item(DISPLAYITEM, display) ;
X error = 1 ;
X set_item(OPITEM, "CLR") ;
X }
X
X------- functions.c -------
X*** /tmp/da5217 Fri Feb 2 13:01:43 1990
X--- functions.c Fri Feb 2 12:25:38 1990
X***************
X*** 18,27 ****
X--- 18,34 ----
X * reported to me then an attempt will be made to fix them.
X */
X
X+ #include <stdio.h>
X+ #include <strings.h>
X+ #include <math.h>
X #include "calctool.h"
X #include "color.h"
X #include "extern.h"
X
X+ extern double acos(), acosh(), asin(), asinh(), atan(), atanh() ;
X+ extern double cos(), cosh(), exp(), fabs(), log(), log10(), pow() ;
X+ extern double sin(), sinh(), sqrt(), tan(), tanh() ;
X+
X BOOLEAN ibool() ;
X double setbool() ;
X
X***************
X*** 185,202 ****
X double a ;
X int i ;
X
X! if (val == (int) val)
X {
X i = val ;
X a = 1.0 ;
X while ((i > 0) && (a != HUGE)) a *= (float) i-- ;
X }
X- else
X- {
X- a = gamma(val+1) ;
X- a = exp(a) ;
X- if (signgam) a = -a ;
X- }
X return (a) ;
X }
X
X--- 192,204 ----
X double a ;
X int i ;
X
X! a = val ;
X! if (val == (int) val && val > 0.0) /* Only for positive integers. */
X {
X i = val ;
X a = 1.0 ;
X while ((i > 0) && (a != HUGE)) a *= (float) i-- ;
X }
X return (a) ;
X }
X
X***************
X*** 284,290 ****
X break ;
X case '{' : disp_val = exp(disp_val) ; /* e^x */
X break ;
X! case '}' : disp_val = exp(M_LN10*disp_val) ; /* 10^x */
X break ;
X case 'N' : disp_val = log(disp_val) ; /* ln */
X break ;
X--- 286,292 ----
X break ;
X case '{' : disp_val = exp(disp_val) ; /* e^x */
X break ;
X! case '}' : disp_val = exp(LN10*disp_val) ; /* 10^x */
X break ;
X case 'N' : disp_val = log(disp_val) ; /* ln */
X break ;
X***************
X*** 506,513 ****
X
X if (!inverse)
X {
X! if (ttype == DEG) temp = disp_val * M_PI / 180.0 ;
X! else if (ttype == GRAD) temp = disp_val * M_PI / 200.0 ;
X else temp = disp_val ;
X
X if (!hyperbolic)
X--- 508,515 ----
X
X if (!inverse)
X {
X! if (ttype == DEG) temp = disp_val * PI / 180.0 ;
X! else if (ttype == GRAD) temp = disp_val * PI / 200.0 ;
X else temp = disp_val ;
X
X if (!hyperbolic)
X***************
X*** 553,560 ****
X case CCTRL('t') : disp_val = atanh(disp_val) ; /* atanh */
X }
X
X! tresults[(int) DEG] = disp_val * 180.0 / M_PI ;
X! tresults[(int) GRAD] = disp_val * 200.0 / M_PI ;
X tresults[(int) RAD] = disp_val ;
X }
X
X--- 555,562 ----
X case CCTRL('t') : disp_val = atanh(disp_val) ; /* atanh */
X }
X
X! tresults[(int) DEG] = disp_val * 180.0 / PI ;
X! tresults[(int) GRAD] = disp_val * 200.0 / PI ;
X tresults[(int) RAD] = disp_val ;
X }
X
X
X------- graphics.c -------
X*** /tmp/da5220 Fri Feb 2 13:01:45 1990
X--- graphics.c Tue Jan 30 12:09:19 1990
X***************
X*** 14,19 ****
X--- 14,21 ----
X * reported to me then an attempt will be made to fix them.
X */
X
X+ #include <stdio.h>
X+ #include <strings.h>
X #include "calctool.h"
X #include "color.h"
X #include "extern.h"
X***************
X*** 73,80 ****
X drawbox(column*(BWIDTH+BGAP)+BBORDER,
X DISPLAY+row*(BHEIGHT+BGAP)+BBORDER, BWIDTH, BHEIGHT) ;
X fillbox(column*(BWIDTH+BGAP)+BBORDER+1,
X! DISPLAY+row*(BHEIGHT+BGAP)+BBORDER+1, KEYCANVAS,
X! 42, 50, 1, color) ;
X }
X else
X {
X--- 75,81 ----
X drawbox(column*(BWIDTH+BGAP)+BBORDER,
X DISPLAY+row*(BHEIGHT+BGAP)+BBORDER, BWIDTH, BHEIGHT) ;
X fillbox(column*(BWIDTH+BGAP)+BBORDER+1,
X! DISPLAY+row*(BHEIGHT+BGAP)+BBORDER+1, 42, 50, 1, color) ;
X }
X else
X {
X***************
X*** 82,96 ****
X DISPLAY+row*(BHEIGHT+BGAP)+BBORDER+26, 34, 21) ;
X color = (iscolor) ? buttons[n].color : BLACK ;
X fillbox(column*(BWIDTH+BGAP)+BBORDER+6,
X! DISPLAY+row*(BHEIGHT+BGAP)+BBORDER+27, KEYCANVAS,
X! 32, 19, 1, color) ;
X }
X but_text(row, column, portion, state) ;
X }
X
X
X! fillbox(x, y, window, width, height, boundry, color)
X! enum can_type window ;
X int x, y, width, height, boundry, color ;
X {
X if (boundry)
X--- 83,95 ----
X DISPLAY+row*(BHEIGHT+BGAP)+BBORDER+26, 34, 21) ;
X color = (iscolor) ? buttons[n].color : BLACK ;
X fillbox(column*(BWIDTH+BGAP)+BBORDER+6,
X! DISPLAY+row*(BHEIGHT+BGAP)+BBORDER+27, 32, 19, 1, color) ;
X }
X but_text(row, column, portion, state) ;
X }
X
X
X! fillbox(x, y, width, height, boundry, color)
X int x, y, width, height, boundry, color ;
X {
X if (boundry)
X***************
X*** 199,205 ****
X color = (portion) ? WHITE : BLACK ;
X fillbox(column*(BWIDTH+BGAP)+BBORDER+6,
X DISPLAY+row*(BHEIGHT+BGAP)+BBORDER+5+(portion*22),
X! KEYCANVAS, 32, 19, portion, color) ;
X but_text(row, column, portion, state) ;
X }
X }
X--- 198,204 ----
X color = (portion) ? WHITE : BLACK ;
X fillbox(column*(BWIDTH+BGAP)+BBORDER+6,
X DISPLAY+row*(BHEIGHT+BGAP)+BBORDER+5+(portion*22),
X! 32, 19, portion, color) ;
X but_text(row, column, portion, state) ;
X }
X }
X
X------- news.c -------
X*** /tmp/da5223 Fri Feb 2 13:01:46 1990
X--- news.c Tue Jan 30 12:47:57 1990
X***************
X*** 14,19 ****
X--- 14,21 ----
X * reported to me then an attempt will be made to fix them.
X */
X
X+ #include <stdio.h>
X+ #include <strings.h>
X #include "calctool.h"
X #include "color.h"
X #include "extern.h"
X
X------- sunview.c -------
X*** /tmp/da5226 Fri Feb 2 13:01:47 1990
X--- sunview.c Tue Jan 30 12:30:27 1990
X***************
X*** 14,19 ****
X--- 14,21 ----
X * reported to me then an attempt will be made to fix them.
X */
X
X+ #include <stdio.h>
X+ #include <strings.h>
X #include "calctool.h"
X #include "color.h"
X #include "extern.h"
X***************
X*** 27,32 ****
X--- 29,35 ----
X #define NOTIFY_INTERPOSE_DESTROY_FUNC (void) notify_interpose_destroy_func
X #define PW_SETCMSNAME (void) pw_setcmsname
X #define PW_PUTCOLORMAP (void) pw_putcolormap
X+ #define PW_REVERSEVIDEO (void) pw_reversevideo
X #define PW_TTEXT (void) pw_ttext
X #define PW_VECTOR (void) pw_vector
X #define PW_WRITEBACKGROUND (void) pw_writebackground
X***************
X*** 457,463 ****
X
X calc_colorsetup(red, green, blue) ;
X PW_PUTCOLORMAP(cpw, 0, CALC_COLORSIZE, red, green, blue) ;
X! if (inv_video) pw_reversevideo(cpw, 0, CALC_COLORSIZE) ;
X
X if (iscolor)
X {
X--- 460,466 ----
X
X calc_colorsetup(red, green, blue) ;
X PW_PUTCOLORMAP(cpw, 0, CALC_COLORSIZE, red, green, blue) ;
X! if (inv_video) PW_REVERSEVIDEO(cpw, 0, CALC_COLORSIZE) ;
X
X if (iscolor)
X {
X
X------- x11.c -------
X*** /tmp/da5229 Fri Feb 2 13:01:48 1990
X--- x11.c Mon Jan 29 23:10:33 1990
X***************
X*** 14,19 ****
X--- 14,20 ----
X * reported to me then an attempt will be made to fix them.
X */
X
X+ #include <stdio.h>
X #include "calctool.h"
X #include "color.h"
X #include "extern.h"
X***************
X*** 91,117 ****
X
X close_frame()
X {
X! Atom a ;
X! XClientMessageEvent ev ;
X!
X! iconic = 1 ;
X! rstate = 0 ;
X!
X! if (dpy->atoms->wm_change_state == None)
X! {
X! a = XInternAtom (dpy, "WM_CHANGE_STATE", False) ;
X! if (a == None) return ;
X! dpy->atoms->wm_change_state = a ;
X! }
X!
X! ev.type = ClientMessage ;
X! ev.window = frame ;
X! ev.message_type = dpy->atoms->wm_change_state ;
X! ev.format = 32 ;
X! ev.data.l[0] = IconicState ;
X! XSendEvent(dpy, RootWindow(dpy, screen), False,
X! SubstructureRedirectMask | SubstructureNotifyMask,
X! (XEvent *) &ev) ;
X }
X
X
X--- 92,100 ----
X
X close_frame()
X {
X! #ifndef NOT_R4
X! XIconifyWindow(dpy, frame, screen) ;
X! #endif /*NOT_R4*/
X }
X
X
X
X------- calctool.h -------
X*** /tmp/da5232 Fri Feb 2 13:01:49 1990
X--- calctool.h Fri Feb 2 12:54:28 1990
X***************
X*** 14,31 ****
X * reported to me then an attempt will be made to fix them.
X */
X
X! #include <stdio.h>
X! #include <strings.h>
X! #include <ctype.h>
X! #include <signal.h>
X! #include <sys/types.h>
X! #include <sys/ioctl.h>
X! #include <sys/time.h>
X! #include <pwd.h>
X! #include <math.h>
X
X- char *getenv(), *sprintf() ;
X-
X #define CLOSE (void) close /* To make lint happy. */
X #define FCLOSE (void) fclose
X #define FFLUSH (void) fflush
X--- 14,21 ----
X * reported to me then an attempt will be made to fix them.
X */
X
X! char *getenv(), *getwd(), *sprintf() ;
X
X #define CLOSE (void) close /* To make lint happy. */
X #define FCLOSE (void) fclose
X #define FFLUSH (void) fflush
X***************
X*** 91,96 ****
X--- 81,102 ----
X
X enum trig_type { DEG, GRAD, RAD } ; /* Trigonometric types. */
X
X+ /* Mathematical constants used by the routines in functions.c
X+ * These should be declared in math.h, but just in case....
X+ */
X+
X+ #ifndef HUGE
X+ #define HUGE 1.701411733192644270e38
X+ #endif /*HUGE*/
X+
X+ #ifndef LN10
X+ #define LN10 2.30258509299404568402
X+ #endif /*LN10*/
X+
X+ #ifndef PI
X+ #define PI 3.14159265358979323846
X+ #endif /*PI*/
X+
X #define BBORDER 10 /* No of pixels in border. */
X #define BCOLS 6 /* No of columns of buttons. */
X #define BGAP 5 /* No of pixels between buttons. */
X***************
X*** 104,112 ****
X #define EQUAL !strcmp /* For character comparisons. */
X #define EXTRA 5 /* Extra useful character definitions. */
X
X! #ifndef HELPGIVEN
X #define HELPNAME "calctool.help"
X! #endif /*HELPGIVEN*/
X
X #define ICONHEIGHT 64 /* Height of calctool icon. */
X #define ICONWIDTH 42 /* Width of calctool icon. */
X--- 110,118 ----
X #define EQUAL !strcmp /* For character comparisons. */
X #define EXTRA 5 /* Extra useful character definitions. */
X
X! #ifndef HELPNAME
X #define HELPNAME "calctool.help"
X! #endif /*HELPNAME*/
X
X #define ICONHEIGHT 64 /* Height of calctool icon. */
X #define ICONWIDTH 42 /* Width of calctool icon. */
X***************
X*** 119,131 ****
X #endif /*MAXLINE*/
X
X #define MAXMENUS 8 /* Maximum number of popup menus. */
X #define MAXREGS 10 /* Maximum number of memory registers. */
X #define MAXVKEYS 7 /* Number of valid keys after an error. */
X #define MIN(x,y) ((x) < (y) ? (x) : (y))
X
X! #ifndef NEWSGIVEN
X! #define NEWSNAME "calctool.ps"
X! #endif /*NEWSGIVEN*/
X
X #define NOBUTTONS BROWS * BCOLS
X
X--- 125,145 ----
X #endif /*MAXLINE*/
X
X #define MAXMENUS 8 /* Maximum number of popup menus. */
X+
X+ #ifndef MAXPATHLEN
X+ #define MAXPATHLEN 1024 /* Longest possible path length. */
X+ #endif /*MAXPATHLEN*/
X+
X #define MAXREGS 10 /* Maximum number of memory registers. */
X #define MAXVKEYS 7 /* Number of valid keys after an error. */
X+
X+ #ifndef MIN
X #define MIN(x,y) ((x) < (y) ? (x) : (y))
X+ #endif /*MIN*/
X
X! #ifndef NEWSFILE
X! #define NEWSFILE "calctool.ps"
X! #endif /*NEWSFILE*/
X
X #define NOBUTTONS BROWS * BCOLS
X
X***************
X*** 133,141 ****
X #define index strchr
X #endif /*NOINDEX*/
X
X! #ifndef RCGIVEN
X #define RCNAME ".calctoolrc"
X! #endif /*RCGIVEN*/
X
X #ifndef NO_4_3SIGNAL
X #define SIGRET void
X--- 147,155 ----
X #define index strchr
X #endif /*NOINDEX*/
X
X! #ifndef RCNAME
X #define RCNAME ".calctoolrc"
X! #endif /*RCNAME*/
X
X #ifndef NO_4_3SIGNAL
X #define SIGRET void
X
X------- Makefile -------
X*** /tmp/da5235 Fri Feb 2 13:01:50 1990
X--- Makefile Fri Feb 2 12:24:19 1990
X***************
X*** 28,35 ****
X # It can also be specified here by uncommenting the following
X # macro definition and setting appropriately.
X #
X! #HELPNAME = -DHELPGIVEN -DHELPNAME=\"$(LIBDIR)/calctool.help\"
X #-----------------------------------------------------------------------
X # Not all machines have the index() string library function. If you
X # don't have this function then you should uncomment the NOINDEX
X # definition below.
X--- 28,44 ----
X # It can also be specified here by uncommenting the following
X # macro definition and setting appropriately.
X #
X! #HELPNAME = -DHELPNAME=\"$(LIBDIR)/calctool.help\"
X #-----------------------------------------------------------------------
X+ # An attempt is now made at using a portable math library. You should
X+ # carefully setup the following three definitions.
X+ #
X+ # One of the things this is not the same on all machines is the
X+ # floating point format. If your machine does *not* use IEEE floating
X+ # point format, then the following definition needs to be uncommented.
X+ #
X+ #NOIEEE = -DIEEE
X+ #-----------------------------------------------------------------------
X # Not all machines have the index() string library function. If you
X # don't have this function then you should uncomment the NOINDEX
X # definition below.
X***************
X*** 47,53 ****
X # should be uncommented and set to the default location for the calctool
X # PostScript file, which is downloaded to the NeWS server.
X #
X! #NEWSFILE = -DNEWSGIVEN -DNEWSFILE=\"$(LIBDIR)/calctool.ps\"
X #-----------------------------------------------------------------------
X # Calctool loads a run control file when it starts. By default
X # this file is called ".calctoolrc", and is taken from the users
X--- 56,62 ----
X # should be uncommented and set to the default location for the calctool
X # PostScript file, which is downloaded to the NeWS server.
X #
X! #NEWSFILE = -DNEWSFILE=\"$(LIBDIR)/calctool.ps\"
X #-----------------------------------------------------------------------
X # Calctool loads a run control file when it starts. By default
X # this file is called ".calctoolrc", and is taken from the users
X***************
X*** 56,62 ****
X # can also be specified here by uncommenting the following macro
X # definition and setting appropriately.
X #
X! #RCNAME = -DRCGIVEN -DRCNAME=\".calctoolrc\"
X #-----------------------------------------------------------------------
X # If you not running under a BSD4.3 derived system, the parameters
X # to the select call are different, and this definition should be
X--- 65,71 ----
X # can also be specified here by uncommenting the following macro
X # definition and setting appropriately.
X #
X! #RCNAME = -DRCNAME=\".calctoolrc\"
X #-----------------------------------------------------------------------
X # If you not running under a BSD4.3 derived system, the parameters
X # to the select call are different, and this definition should be
X***************
X*** 66,77 ****
X #-----------------------------------------------------------------------
X # Note that with SunOS v4.x it is possible to create a shared library
X # for the calctool routines. If you wish to do this, then you should
X! # change LIBNAME below to libcalctool.so.1.1, and uncomment the SHLIB
X! # line.
X #
X LIBNAME = libcalctool.a
X #
X #LIBNAME = libcalctool.so.1.1
X #SHLIB = -pic
X #-----------------------------------------------------------------------
X # If you are not running under a BSD4.3 derived system, then the
X--- 75,89 ----
X #-----------------------------------------------------------------------
X # Note that with SunOS v4.x it is possible to create a shared library
X # for the calctool routines. If you wish to do this, then you should
X! # uncomment the other three definitions below, commented out the initial
X! # two.
X #
X LIBNAME = libcalctool.a
X+ CALCLIB = $(LIBNAME)
X #
X+ # Shared library definitions. Uncomment if required.
X #LIBNAME = libcalctool.so.1.1
X+ #CALCLIB = -L. -lcalctool
X #SHLIB = -pic
X #-----------------------------------------------------------------------
X # If you are not running under a BSD4.3 derived system, then the
X***************
X*** 99,104 ****
X--- 111,122 ----
X #
X #X11INCDIR = -I$(OPENWINHOME)/include
X #X11LIBDIR = -L$(OPENWINHOME)/lib
X+ #
X+ # If you are not running X11R4, then you should uncomment the following
X+ # definition. This is needed to disable the call to XIconifyWindow in
X+ # the routine close_frame. This routines is new in X11R4.
X+ #
X+ #NOT_R4 = -DNOT_R4
X #-------------------------------------------------------------------------
X # If you are compiling the XView version, then the following two lines
X # should be uncommented.
X***************
X*** 118,124 ****
X #
X # Compilation flags and standard macro definitions.
X #
X! CFLAGS = -g $(HELPNAME) $(NEWSFILE) $(NOINDEX) $(RCNAME) \
X $(SELTYPE) $(SHLIB) $(SIGRET) $(SUN4_KEYBOARD) \
X $(TTEXT) $(MGRPARAM) $(MGRINCDIR) $(X11INCDIR) \
X $(XVIEWINCDIR)
X--- 136,142 ----
X #
X # Compilation flags and standard macro definitions.
X #
X! CFLAGS = -g $(HELPNAME) $(NEWSFILE) $(NOIEEE) $(NOINDEX) $(RCNAME) \
X $(SELTYPE) $(SHLIB) $(SIGRET) $(SUN4_KEYBOARD) \
X $(TTEXT) $(MGRPARAM) $(MGRINCDIR) $(X11INCDIR) \
X $(XVIEWINCDIR)
X***************
X*** 130,155 ****
X
X LIBSRCS = graphics.c display.c functions.c get.c
X LIBOBJS = graphics.o display.o functions.o get.o
X! STDSRCS = calctool.c
X! STDOBJS = calctool.o
X OBJS = $(STDOBJS) $(LIBNAME)
X
X GSRCS = mgr.c news.c sunview.c tty.c x11.c xview.c
X! HDRS = calctool.h color.h extern.h
X IMAGES = calctool.icon calctool.color.icon help.cursor
X OTHERS = Makefile README TODO calctool.help calctool.1 \
X! calctool.ps patchlevel.h .calctoolrc
X
X! SFILES1 = $(LIBSRCS) $(STDSRCS)
X! SFILES2 = $(GSRCS)
X! SFILES3 = $(HDRS) $(IMAGES)
X! SFILES4 = $(OTHERS)
X
X- CALCLIB = -L. -lcalctool
X MGRLIBS = $(CALCLIB) $(MGRHOME)/lib/libmgr.a -lm
X! NEWSLIBS = $(CALCLIB) $$NEWSHOME/lib/libcps.a -lm
X SVIEWLIBS = $(CALCLIB) -lsuntool -lsunwindow -lpixrect -lm
X! TTYLIBS = $(CALCLIB) -ltermcap -lm
X X11LIBS = $(CALCLIB) -lX11 -lm
X XVIEWLIBS = $(CALCLIB) -lxview -lolgx -lX11 -lm
X
X--- 148,173 ----
X
X LIBSRCS = graphics.c display.c functions.c get.c
X LIBOBJS = graphics.o display.o functions.o get.o
X! STDSRCS = calctool.c mathlib.c
X! STDOBJS = calctool.o mathlib.o
X OBJS = $(STDOBJS) $(LIBNAME)
X
X GSRCS = mgr.c news.c sunview.c tty.c x11.c xview.c
X! HDRS = calctool.h color.h extern.h mathlib.h
X IMAGES = calctool.icon calctool.color.icon help.cursor
X OTHERS = Makefile README TODO calctool.help calctool.1 \
X! CHANGES calctool.ps patchlevel.h .calctoolrc
X
X! SFILES1 = $(STDSRCS)
X! SFILES2 = $(LIBSRCS)
X! SFILES3 = $(GSRCS)
X! SFILES4 = $(HDRS) $(IMAGES)
X! SFILES5 = $(OTHERS)
X
X MGRLIBS = $(CALCLIB) $(MGRHOME)/lib/libmgr.a -lm
X! NEWSLIBS = $(CALCLIB) $$OPENWINHOME/lib/libcps.a -lm
X SVIEWLIBS = $(CALCLIB) -lsuntool -lsunwindow -lpixrect -lm
X! TTYLIBS = $(CALCLIB) -ltermcap -lm
X X11LIBS = $(CALCLIB) -lX11 -lm
X XVIEWLIBS = $(CALCLIB) -lxview -lolgx -lX11 -lm
X
X***************
X*** 230,235 ****
X--- 248,254 ----
X shar.script $(SFILES2) > archive.2
X shar.script $(SFILES3) > archive.3
X shar.script $(SFILES4) > archive.4
X+ shar.script $(SFILES5) > archive.5
X
X create: SCCS
X -sccs create $(STDSRCS) $(GSRCS) $(HDRS) $(IMAGES) $(OTHERS)
X***************
X*** 236,238 ****
X--- 255,270 ----
X SCCS:
X -mkdir SCCS
X chmod 755 SCCS
X+
X+ calctool.o: calctool.c patchlevel.h color.h calctool.h
X+ display.o: display.c calctool.h color.h extern.h
X+ functions.o: functions.c calctool.h color.h extern.h
X+ get.o: get.c patchlevel.h calctool.h color.h extern.h
X+ graphics.o: graphics.c calctool.h color.h extern.h
X+ mathlib.o: mathlib.c mathlib.h calctool.h
X+ mgr.o: mgr.c calctool.h color.h extern.h
X+ news.o: news.c calctool.h color.h extern.h
X+ sunview.o: sunview.c calctool.h color.h extern.h
X+ tty.o: tty.c calctool.h color.h extern.h
X+ x11.o: x11.c calctool.h color.h extern.h
X+ xview.o: xview.c calctool.h color.h extern.h
X
X------- README -------
X*** /tmp/da5238 Fri Feb 2 13:01:51 1990
X--- README Tue Jan 30 00:07:24 1990
X***************
X*** 39,44 ****
X--- 39,62 ----
X defined by the Makefile macro definition LIBDIR, so you might need root
X permission to successfully do this.
X
X+
X+ Portable maths routines.
X+ ------------------------
X+
X+ It is hoped that your system supplies all the mathematical functions
X+ required by calctool. If not then, it is possible to use the needed
X+ ones from the portable math library routines that comes with these
X+ sources.
X+
X+ In mathlib.h, there is one definition for each routine used by calctool.
X+ These are commented out by default to signify that this system has that
X+ routine in it's maths library. If you are missing any, then uncomment
X+ the appropriate definitions.
X+
X+
X+ Todo.
X+ -----
X+
X There is a TODO file included, which lists the current known bugs, and
X a set of possible enhancements, some relatively trivial, others that
X are fairly major enhancements which would make interesting future
X***************
X*** 46,51 ****
X--- 64,70 ----
X
X
X Acknowledgements.
X+ -----------------
X
X Thanks to Ed Falk at Sun Microsystems (Mountain View) for most of the
X basic arithmetical algorithms used, to Andrew Nicholson for revising the
X***************
X*** 55,65 ****
X
X Thanks go also to James Buster, David Weaver, Steve Damron, Mike Bender,
X Charles Tierney, Trevor Watson, Marla Berg, David Hough, Jeff Donsbach,
X! Mel Melchner, Peter Allott, Skip Gilbrech and Keith McNeill for bug
X! reports and/or bug fixes plus sugggested enhancements.
X
X Suggestions for furthur improvement would be most welcome, plus bugs,
X comments and flames.
X
X! Rich Burridge, DOMAIN: richb@sunaus.oz.au
X! PHONE: +61 2 413 2666 UUCP: {uunet,mcvax,ukc}!munnari!sunaus.oz!richb
X--- 74,86 ----
X
X Thanks go also to James Buster, David Weaver, Steve Damron, Mike Bender,
X Charles Tierney, Trevor Watson, Marla Berg, David Hough, Jeff Donsbach,
X! Mel Melchner, Peter Allott, Skip Gilbrech, Tom Friedel, Keith McNeill
X! and Stephen Frede for bug reports and/or bug fixes plus sugggested
X! enhancements.
X
X Suggestions for furthur improvement would be most welcome, plus bugs,
X comments and flames.
X
X! Rich Burridge, DOMAIN: richb@Aus.Sun.COM
X! PHONE: +61 2 413 2666 ACSNET: richb@sunaus.sun.oz
X! UUCP: {uunet,mcvax,ukc}!munnari!sunaus.oz!richb
X
X------- calctool.1 -------
X*** /tmp/da5241 Fri Feb 2 13:01:52 1990
X--- calctool.1 Mon Jan 29 23:55:44 1990
X***************
X*** 259,265 ****
X .IP "\fB1/x [ R ]\fP" 18
X Return the value of 1 divided by the current entry.
X .IP "\fBx! [ ! ]\fP" 18
X! Return the factorial of the current entry.
X .IP "\fBx^2 [ @ ]\fP" 18
X Return the square of the current entry.
X .SS Number Manipulation Operators.
X--- 259,266 ----
X .IP "\fB1/x [ R ]\fP" 18
X Return the value of 1 divided by the current entry.
X .IP "\fBx! [ ! ]\fP" 18
X! Return the factorial of the current entry. Note that this only works
X! for positive integers.
X .IP "\fBx^2 [ @ ]\fP" 18
X Return the square of the current entry.
X .SS Number Manipulation Operators.
X
X------- patchlevel.h -------
X*** /tmp/da5244 Fri Feb 2 13:01:53 1990
X--- patchlevel.h Tue Jan 23 10:00:38 1990
X***************
X*** 14,17 ****
X * reported to me then an attempt will be made to fix them.
X */
X
X! #define PATCHLEVEL 5
X--- 14,17 ----
X * reported to me then an attempt will be made to fix them.
X */
X
X! #define PATCHLEVEL 6
X
X------- mgr.c -------
X*** /tmp/da5247 Fri Feb 2 13:01:53 1990
X--- mgr.c Tue Jan 30 00:21:25 1990
X***************
X*** 14,19 ****
X--- 14,24 ----
X * reported to me then an attempt will be made to fix them.
X */
X
X+ #include <stdio.h>
X+ #include <signal.h>
X+ #include <sys/types.h>
X+ #include <sys/ioctl.h>
X+ #include <sys/time.h>
X #include "dump.h"
X #include "term.h"
X #include "calctool.h"
X
X------- tty.c -------
X*** /tmp/da5250 Fri Feb 2 13:01:54 1990
X--- tty.c Tue Jan 30 12:23:19 1990
X***************
X*** 1,3 ****
X--- 1,4 ----
X+ /*LINTLIBRARY*/
X
X /* %Z%%M% %I% %E%
X *
X***************
X*** 14,19 ****
X--- 15,26 ----
X * reported to me then an attempt will be made to fix them.
X */
X
X+ #include <stdio.h>
X+ #include <strings.h>
X+ #include <signal.h>
X+ #include <sys/types.h>
X+ #include <sys/ioctl.h>
X+ #include <sys/time.h>
X #include "calctool.h"
X #include "color.h"
X #include "extern.h"
X***************
X*** 105,111 ****
X /*ARGSUSED*/
X do_menu(mtype) /* Popup appropriate menu (null routine). */
X enum menu_type mtype ;
X! {}
X
X
X do_move(x, y) /* Move to character position (x, y). */
X--- 112,120 ----
X /*ARGSUSED*/
X do_menu(mtype) /* Popup appropriate menu (null routine). */
X enum menu_type mtype ;
X! {
X! return 0 ; /* To make lint happy. */
X! }
X
X
X do_move(x, y) /* Move to character position (x, y). */
X***************
X*** 234,240 ****
X case 3 : key[0] = ' ' ;
X STRNCPY(&key[1], str, 3) ;
X }
X! key[5] = '\0' ;
X STRCPY(str, key) ;
X if (ty % 3) invert = 1 ;
X }
X--- 243,249 ----
X case 3 : key[0] = ' ' ;
X STRNCPY(&key[1], str, 3) ;
X }
X! key[4] = '\0' ;
X STRCPY(str, key) ;
X if (ty % 3) invert = 1 ;
X }
X
X------- xview.c -------
X*** /tmp/da5253 Fri Feb 2 13:01:55 1990
X--- xview.c Mon Jan 29 23:11:47 1990
X***************
X*** 14,19 ****
X--- 14,20 ----
X * reported to me then an attempt will be made to fix them.
X */
X
X+ #include <stdio.h>
X #include "calctool.h"
X #include "color.h"
X #include "extern.h"
X
X------- get.c -------
X*** /tmp/da5256 Fri Feb 2 13:01:56 1990
X--- get.c Fri Feb 2 12:54:34 1990
X***************
X*** 14,19 ****
X--- 14,23 ----
X * reported to me then an attempt will be made to fix them.
X */
X
X+ #include <stdio.h>
X+ #include <strings.h>
X+ #include <sys/param.h>
X+ #include <pwd.h>
X #include "patchlevel.h"
X #include "calctool.h"
X #include "color.h"
X***************
X*** 281,286 ****
X--- 285,291 ----
X {
X char *home ; /* Pathname for users home directory. */
X char name[MAXLINE] ; /* Full name of users .calctoolrc file. */
X+ char pathname[MAXPATHLEN] ; /* Current working directory. */
X int n ;
X struct passwd *entry ;
X
X***************
X*** 299,305 ****
X SPRINTF(name, "%s/%s", home, RCNAME) ;
X get_rcfile(name) ; /* Read .calctoolrc from users home directory. */
X
X! SPRINTF(name, "%s/%s", getcwd((char *) NULL, MAXLINE), RCNAME) ;
X get_rcfile(name) ; /* Read .calctoolrc file from current directory. */
X }
X
X--- 304,310 ----
X SPRINTF(name, "%s/%s", home, RCNAME) ;
X get_rcfile(name) ; /* Read .calctoolrc from users home directory. */
X
X! SPRINTF(name, "%s/%s", getwd(pathname), RCNAME) ;
X get_rcfile(name) ; /* Read .calctoolrc file from current directory. */
X }
X
X
X------- display.c -------
X*** /tmp/da5259 Fri Feb 2 13:01:57 1990
X--- display.c Tue Jan 30 12:09:53 1990
X***************
X*** 17,22 ****
X--- 17,25 ----
X * reported to me then an attempt will be made to fix them.
X */
X
X+ #include <stdio.h>
X+ #include <strings.h>
X+ #include <math.h>
X #include "calctool.h"
X #include "color.h"
X #include "extern.h"
X
X------- CHANGES -------
X*** /tmp/da5262 Fri Feb 2 13:01:58 1990
X--- CHANGES Fri Feb 2 13:01:01 1990
X***************
X*** 93,95 ****
X--- 93,152 ----
X
X * CHANGES - documented history of the changes made with each new
X calctool patch.
X+
X+ v2.4 - patchlevel 6. - Posted to comp.sources.bugs (January 1990).
X+
X+ Changes:
X+
X+ * The Makefile did not include the CHANGES file in the OTHERS
X+ macro definition.
X+
X+ * If you are not using IEEE floating point format, and you need to
X+ use some of the routines from the portable maths library which is
X+ now included with this distribution, then you need to uncomment
X+ the new NOTIEEE definition in the Makefile.
X+
X+ * From Tom Friedel <tpf@jdyx.atlanta.ga.us>
X+ In the tty driver, in the drawtext routine, there is a character
X+ array of five elements. There was an invalid assignment to key[5].
X+
X+ * From Stephen Frede <stephenf@softway.oz>
X+ The definitions in the Makefile for HELPNAME, NEWSFILE and RCNAME
X+ do not need extra ...GIVEN definitions.
X+
X+ * Calctool.h defined NEWSNAME, when it should have defined NEWSFILE.
X+
X+ * With the X11 driver, the function XIconifyWindow is new to R4.
X+ If you are running an X11 release which is not R4, then you need
X+ to uncomment the new definition in the Makefile called NOT_R4.
X+ This nullifys the function of the calctool OFF key, but hey, you
X+ can also you the window manager to iconify calctool (assuming
X+ your window manager can do that!).
X+
X+ * The factorial function has been adjusted to only work with
X+ positive integers. I couldn't find a public domain or freely
X+ distributable lgamma function written in C.
X+
X+ * fillbox in graphics.c no longer needs to be passed a window
X+ parameter
X+
X+ * From Stephen Frede <stephenf@softway.oz>
X+ Not all linkers have the -L option. The CALCLIB definition in the
X+ Makefile has been adjusted, and moved to the section where
X+ LIBNAME is defined. By default, you shouldn't have to do anything,
X+ but if you are on a Sun and want to use the shared library facility,
X+ then you need to uncomment three definitions.
X+
X+ * From Stephen Frede <stephenf@softway.oz>
X+ getcwd is not present on all systems. I've replaced it with getwd.
X+
X+
X+ New files:
X+
X+ * mathlib.h
X+ mathlib.c - These are an attempt to provide a portable maths
X+ library for those systems that don't have all the
X+ functions that calctool needs.
X+
X+ You should look near the end of mathlib.h for a set
X+ of definitions, then set appropriately.
SHAR_EOF
chmod 0644 patch.6 || echo "restore of patch.6 fails"
set `wc -c patch.6`;Sum=$1
if test "$Sum" != "35130"
then echo original size 35130, current size $Sum;fi
echo "x - extracting mathlib.h (Text)"
sed 's/^X//' << 'SHAR_EOF' > mathlib.h &&
X
X/* @(#)mathlib.h 1.2 90/02/02
X *
X * Definitions used with the portable math library.
X *
X * This is being done because libm.a doesn't appear to be as portable
X * as originally assumed.
X *
X * These routines are taken from two sources:
X *
X * 1/ Fred Fishs' portable maths library which was posted to the
X * comp.sources.unix newsgroup on April 1987.
X *
X * acos, acosh, asin, asinh, atan, atanh, cos, cosh, dabss,
X * exp, log, log10, mod, poly, sin, sinh, sqrt, tan, tanh.
X *
X * 2/ The BSD4.3 maths library found on uunet in
X * bsd_sources/src/usr.lib/libm.
X *
X * pow, pow_p, scalb, logb, copysign, finite, drem, exp__E,
X * log__L.
X *
X * Customised and condensed by Rich Burridge,
X * Sun Microsystems, Australia.
X *
X * Permission is given to distribute these sources, as long as the
X * copyright messages are not removed, and no monies are exchanged.
X *
X * No responsibility is taken for any errors or inaccuracies inherent
X * either to the comments or the code of this program, but if
X * reported to me then an attempt will be made to fix them.
X */
X
X/************************************************************************
X * *
X * N O T I C E *
X * *
X * Copyright Abandoned, 1987, Fred Fish *
X * *
X * This previously copyrighted work has been placed into the *
X * public domain by the author (Fred Fish) and may be freely used *
X * for any purpose, private or commercial. I would appreciate *
X * it, as a courtesy, if this notice is left in all copies and *
X * derivative works. Thank you, and enjoy... *
X * *
X * The author makes no warranty of any kind with respect to this *
X * product and explicitly disclaims any implied warranties of *
X * merchantability or fitness for any particular purpose. *
X * *
X ************************************************************************
X */
X
X/* This file gets included with all of the floating point math
X * library routines when they are compiled. Note that
X * this is the proper place to put machine dependencies
X * whenever possible.
X *
X * It should be pointed out that for simplicity's sake, the
X * environment parameters are defined as floating point constants,
X * rather than octal or hexadecimal initializations of allocated
SHAR_EOF
echo "End of part 1"
echo "File mathlib.h is continued in part 2"
echo "2" > s2_seq_.tmp
exit 0richb@sunaus.oz (Rich Burridge) (02/02/90)
---- Cut Here and unpack ----
#!/bin/sh
# this is part 2 of a multipart archive
# do not concatenate these parts, unpack them in order with /bin/sh
# file mathlib.h continued
#
CurArch=2
if test ! -r s2_seq_.tmp
then echo "Please unpack part 1 first!"
exit 1; fi
( read Scheck
if test "$Scheck" != $CurArch
then echo "Please unpack part $Scheck next!"
exit 1;
else exit 0; fi
) < s2_seq_.tmp || exit 1
echo "x - Continuing file mathlib.h"
sed 's/^X//' << 'SHAR_EOF' >> mathlib.h
X * storage areas. This means that the range of allowed numbers
X * may not exactly match the hardware's capabilities. For example,
X * if the maximum positive double precision floating point number
X * is EXACTLY 1.11...E100 and the constant "MAXDOUBLE is
X * defined to be 1.11E100 then the numbers between 1.11E100 and
X * 1.11...E100 are considered to be undefined. For most
X * applications, this will cause no problems.
X *
X * An alternate method is to allocate a global static "double" variable,
X * say "maxdouble", and use a union declaration and initialization
X * to initialize it with the proper bits for the EXACT maximum value.
X * This was not done because the only compilers available to the
X * author did not fully support union initialization features.
X */
X
Xextern double acos(), acosh(), asin(), asinh(), atan(), atanh() ;
Xextern double cos(), cosh(), exp(), fabs(), log(), log10(), pow() ;
Xextern double sin(), sinh(), sqrt(), tan(), tanh() ;
X
X/*START============start of definitions from <values.h>============START
X *
X * If your system has a /usr/include/values.h, or has another include
X * file which defines:
X * MAXDOUBLE => Maximum double precision number
X * MINDOUBLE => Minimum double precision number
X * DMAXEXP => Maximum exponent of a double
X * DMINEXP => Minimum exponent of a double
X *
X * you can comment out these definitions down to the END line below.
X */
X
X#ifndef BITSPERBYTE
X/* These values work with any binary representation of integers
X * where the high-order bit contains the sign. */
X
X/* a number used normally for size of a shift */
X#if gcos
X#define BITSPERBYTE 9
X#else
X#define BITSPERBYTE 8
X#endif
X#define BITS(type) (BITSPERBYTE * (int)sizeof(type))
X
X/* Various values that describe the binary floating-point representation
X * MAXDOUBLE - the largest double
X * ((_EXPBASE ** DMAXEXP) * (1 - (_EXPBASE ** -DSIGNIF)))
X * MINDOUBLE - the smallest double (_EXPBASE ** (DMINEXP - 1))
X * DMAXEXP - the maximum exponent of a double (as returned by frexp())
X * DMINEXP - the minimum exponent of a double (as returned by frexp())
X * DSIGNIF - the number of significant bits in a double
X * _IEEE - 1 if IEEE standard representation is used
X * _DEXPLEN - the number of bits for the exponent of a double
X * _HIDDENBIT - 1 if high-significance bit of mantissa is implicit
X */
X
X#if u3b || u3b5 || sun
X#define MAXDOUBLE 1.79769313486231470e+308
X#define MINDOUBLE 4.94065645841246544e-324
X#define _IEEE 1
X#define _DEXPLEN 11
X#define _HIDDENBIT 1
X#define DMINEXP (-(DMAXEXP + DSIGNIF - _HIDDENBIT - 3))
X#endif
X
X#if pdp11 || vax
X#define MAXDOUBLE 1.701411834604692293e+38
X#define MINDOUBLE (0.01 * 2.938735877055718770e-37)
X#define _IEEE 0
X#define _DEXPLEN 8
X#define _HIDDENBIT 1
X#define DMINEXP (-DMAXEXP)
X#endif
X
X#if gcos
X#define MAXDOUBLE 1.7014118346046923171e+38
X#define MINDOUBLE 2.9387358770557187699e-39
X#define _IEEE 0
X#define _DEXPLEN 8
X#define _HIDDENBIT 0
X#define DMINEXP (-(DMAXEXP + 1))
X#endif
X
X#define DSIGNIF (BITS(double) - _DEXPLEN + _HIDDENBIT - 1)
X#define DMAXEXP ((1 << _DEXPLEN - 1) - 1 + _IEEE)
X
X#endif /*BITSPERBYTE*/
X
X/*END==============end of definitions from <values.h>==================END*/
X
X
X#define LOG2_MAXDOUBLE (DMAXEXP + 1)
X#define LOG2_MINDOUBLE (DMINEXP - 1)
X#define LOGE_MAXDOUBLE (LOG2_MAXDOUBLE / LOG2E)
X#define LOGE_MINDOUBLE (LOG2_MINDOUBLE / LOG2E)
X
X/*
X * The following are hacks which should be fixed when I understand all
X * the issues a little better. |tanh(TANH_MAXARG)| = 1.0
X */
X
X#define TANH_MAXARG 16
X#define SQRT_MAXDOUBLE 1.304380e19
X
X#define TWOPI (2.0 * PI)
X#define HALFPI (PI / 2.0)
X#define FOURTHPI (PI / 4.0)
X#define SIXTHPI (PI / 6.0)
X#define LOG2E 1.4426950408889634074 /* Log to base 2 of e */
X#define LOG10E 0.4342944819032518276
X#define SQRT2 1.4142135623730950488
X#define SQRT3 1.7320508075688772935
X#define LN2 0.6931471805599453094
X#define LNSQRT2 0.3465735902799726547
X
X/* MC68000 HARDWARE DEPENDENCIES
X *
X * cc -DIEEE => uses IEEE floating point format
X *
X * Apologies for the double negative. I want as few definitions
X * needed in the default case as possible.
X */
X
X#ifndef NOIEEE
X#define X6_UNDERFLOWS (4.209340e-52) /* X**6 almost underflows */
X#define X16_UNDERFLOWS (5.421010e-20) /* X**16 almost underflows */
X#endif /*NOIEEE*/
X
X/* It is hoped that your system supplies all the mathematical functions
X * required by calctool. If not then, it is possible to use the needed
X * ones from the portable math library routines that comes with these
X * sources.
X *
X * There is one definition for each routine used by calctool. These are
X * commented out by default to signify that this system has that routine.
X * If you are missing any, then uncomment the appropriate definitions.
X */
X
X/*#define NEED_ACOS */
X/*#define NEED_ACOSH */
X/*#define NEED_ASIN */
X/*#define NEED_ASINH */
X/*#define NEED_ATAN */
X/*#define NEED_ATANH */
X/*#define NEED_COS */
X/*#define NEED_COSH */
X/*#define NEED_EXP */
X/*#define NEED_LOG */
X/*#define NEED_LOG10 */
X/*#define NEED_POW */
X/*#define NEED_SIN */
X/*#define NEED_SINH */
X/*#define NEED_SQRT */
X/*#define NEED_TAN */
X/*#define NEED_TANH */
SHAR_EOF
echo "File mathlib.h is complete"
chmod 0444 mathlib.h || echo "restore of mathlib.h fails"
set `wc -c mathlib.h`;Sum=$1
if test "$Sum" != "8385"
then echo original size 8385, current size $Sum;fi
echo "x - extracting mathlib.c (Text)"
sed 's/^X//' << 'SHAR_EOF' > mathlib.c &&
X
X/* @(#)mathlib.c 1.2 90/02/02
X *
X * These are the mathematical routines used by calctool.
X *
X * This is being done because libm.a doesn't appear to be as portable
X * as originally assumed.
X *
X * It is hoped that your system supplies all the mathematical functions
X * required by calctool. If not then, it is possible to use the needed
X * ones from this library. Look in mathlib.h for a set of definitions,
X * and uncomment and set appropriately.
X *
X * These routines are taken from two sources:
X *
X * 1/ Fred Fishs' portable maths library which was posted to the
X * comp.sources.unix newsgroup on April 1987.
X *
X * acos, acosh, asin, asinh, atan, atanh, cos, cosh, dabs,
X * exp, log, log10, mod, poly, sin, sinh, sqrt, tan, tanh.
X *
X * 2/ The BSD4.3 maths library found on uunet in
X * bsd_sources/src/usr.lib/libm.
X *
X * pow, pow_p, scalb, logb, copysign, finite, drem, exp__E,
X * log__L.
X *
X * Customised and condensed by Rich Burridge,
X * Sun Microsystems, Australia.
X *
X * Permission is given to distribute these sources, as long as the
X * copyright messages are not removed, and no monies are exchanged.
X *
X * No responsibility is taken for any errors or inaccuracies inherent
X * either to the comments or the code of this program, but if
X * reported to me then an attempt will be made to fix them.
X */
X
X/************************************************************************
X * *
X * N O T I C E *
X * *
X * Copyright Abandoned, 1987, Fred Fish *
X * *
X * This previously copyrighted work has been placed into the *
X * public domain by the author (Fred Fish) and may be freely used *
X * for any purpose, private or commercial. I would appreciate *
X * it, as a courtesy, if this notice is left in all copies and *
X * derivative works. Thank you, and enjoy... *
X * *
X * The author makes no warranty of any kind with respect to this *
X * product and explicitly disclaims any implied warranties of *
X * merchantability or fitness for any particular purpose. *
X * *
X ************************************************************************
X */
X
X#include <stdio.h>
X#include <errno.h>
X#include "mathlib.h"
X#include "calctool.h"
X
Xdouble atan(), cos(), cosh(), dabs(), exp(), frexp() ;
Xdouble ldexp(), log(), mod(), modf(), poly(), sin() ;
Xdouble sinh(), sqrt() ;
X
Xextern double frexp(), ldexp(), modf() ;
X
X
X/* FUNCTION
X *
X * acos double precision arc cosine
X *
X * DESCRIPTION
X *
X * Returns double precision arc cosine of double precision
X * floating point argument.
X *
X * USAGE
X *
X * double acos(x)
X * double x ;
X *
X * REFERENCES
X *
X * Fortran IV-plus user's guide, Digital Equipment Corp. pp B-1.
X *
X * RESTRICTIONS
X *
X * The maximum relative error for the approximating polynomial
X * in atan is 10**(-16.82). However, this assumes exact arithmetic
X * in the polynomial evaluation. Additional rounding and
X * truncation errors may occur as the argument is reduced
X * to the range over which the polynomial approximation
X * is valid, and as the polynomial is evaluated using
X * finite-precision arithmetic.
X *
X * INTERNALS
X *
X * Computes arccosine(x) from:
X *
X * (1) If x < -1.0 or x > +1.0 then call
X * doerr and return 0.0 by default.
X *
X * (2) If x = 0.0 then acos(x) = PI/2.
X *
X * (3) If x = 1.0 then acos(x) = 0.0
X *
X * (4) If x = -1.0 then acos(x) = PI.
X *
X * (5) If 0.0 < x < 1.0 then
X * acos(x) = atan(Y)
X * Y = sqrt[1-(x**2)] / x
X *
X * (6) If -1.0 < x < 0.0 then
X * acos(x) = atan(Y) + PI
X * Y = sqrt[1-(x**2)] / x
X */
X
X#ifdef NEED_ACOS
Xdouble
Xacos(x)
Xdouble x ;
X{
X double y ;
X auto double retval ;
X
X if (x > 1.0 || x < -1.0)
X {
X doerr("acos", "DOMAIN", EDOM) ;
X retval = 0.0 ;
X }
X else if (x == 0.0) retval = HALFPI ;
X else if (x == 1.0) retval = 0.0 ;
X else if (x == -1.0) retval = PI ;
X else
X {
X y = atan(sqrt(1.0 - (x * x)) / x) ;
X if (x > 0.0) retval = y ;
X else retval = y + PI ;
X }
X return(retval) ;
X}
X#endif /*NEED_ACOS*/
X
X
X/* FUNCTION
X *
X * acosh double precision hyperbolic arc cosine
X *
X * DESCRIPTION
X *
X * Returns double precision hyperbolic arc cosine of double precision
X * floating point number.
X *
X * USAGE
X *
X * double acosh(x)
X * double x ;
X *
X * RESTRICTIONS
X *
X * The range of the ACOSH function is all real numbers greater
X * than or equal to 1.0 however large arguments may cause
X * overflow in the x squared portion of the function evaluation.
X *
X * For precision information refer to documentation of the
X * floating point library primatives called.
X *
X * INTERNALS
X *
X * Computes acosh(x) from:
X *
X * 1. If x < 1.0 then report illegal
X * argument and return zero.
X *
X * 2. If x > sqrt(MAXDOUBLE) then
X * set x = sqrt(MAXDOUBLE and
X * continue after reporting overflow.
X *
X * 3. acosh(x) = log [x+sqrt(x**2 - 1)]
X */
X
X#ifdef NEED_ACOSH
Xdouble
Xacosh(x)
Xdouble x ;
X{
X auto double retval ;
X
X if (x < 1.0)
X {
X doerr("acosh", "DOMAIN", ERANGE) ;
X retval = 0.0 ;
X }
X else if (x > SQRT_MAXDOUBLE)
X {
X doerr("acosh", "OVERFLOW", ERANGE) ;
X x = SQRT_MAXDOUBLE ;
X retval = log(2* SQRT_MAXDOUBLE) ;
X }
X else retval = log(x + sqrt(x*x - 1.0)) ;
X return(retval) ;
X}
X#endif /*NEED_ACOSH*/
X
X
X/*
X * FUNCTION
X *
X * asin double precision arc sine
X *
X * DESCRIPTION
X *
X * Returns double precision arc sine of double precision
X * floating point argument.
X *
X * If argument is less than -1.0 or greater than +1.0, calls
X * doerr with a DOMAIN error. If doerr does not handle
X * the error then prints error message and returns 0.
X *
X * USAGE
X *
X * double asin(x)
X * double x ;
X *
X * REFERENCES
X *
X * Fortran IV-plus user's guide, Digital Equipment Corp. pp B-2.
X *
X * RESTRICTIONS
X *
X * For precision information refer to documentation of the floating
X * point library primatives called.
X *
X * INTERNALS
X *
X * Computes arcsine(x) from:
X *
X * (1) If x < -1.0 or x > +1.0 then
X * call doerr and return 0.0 by default.
X *
X * (2) If x = 0.0 then asin(x) = 0.0
X *
X * (3) If x = 1.0 then asin(x) = PI/2.
X *
X * (4) If x = -1.0 then asin(x) = -PI/2
X *
X * (5) If -1.0 < x < 1.0 then
X * asin(x) = atan(y)
X * y = x / sqrt[1-(x**2)]
X */
X
X#ifdef NEED_ASIN
Xdouble
Xasin(x)
Xdouble x ;
X{
X auto double retval ;
X
X if ( x > 1.0 || x < -1.0)
X {
X doerr("asin", "DOMAIN", EDOM) ;
X retval = 0.0 ;
X }
X else if (x == 0.0) retval = 0.0 ;
X else if (x == 1.0) retval = HALFPI ;
X else if (x == -1.0) retval = -HALFPI ;
X else retval = atan(x / sqrt (1.0 - (x * x))) ;
X return(retval) ;
X}
X#endif /*NEED_ASIN*/
X
X
X/* FUNCTION
X *
X * asinh double precision hyperbolic arc sine
X *
X * DESCRIPTION
X *
X * Returns double precision hyperbolic arc sine of double precision
X * floating point number.
X *
X * USAGE
X *
X * double asinh(x)
X * double x ;
X *
X * RESTRICTIONS
X *
X * The domain of the ASINH function is the entire real axis
X * however the evaluation of x squared may cause overflow
X * for large magnitudes.
X *
X * For precision information refer to documentation of the
X * floating point library routines called.
X *
X * INTERNALS
X *
X * Computes asinh(x) from:
X *
X * 1. Let xmax = sqrt(MAXDOUBLE - 1)
X *
X * 2. If x < -xmax or xmax < x then
X * let x = xmax and flag overflow.
X *
X * 3. asinh(x) = log [x+sqrt(x**2 + 1)]
X */
X
X#ifdef NEED_ASINH
Xdouble
Xasinh(x)
Xdouble x ;
X{
X auto double retval ;
X
X if (x < -SQRT_MAXDOUBLE || x > SQRT_MAXDOUBLE)
X {
X doerr("asinh", "OVERFLOW", ERANGE) ;
X retval = log(2 * SQRT_MAXDOUBLE) ;
X }
X else retval = log(x + sqrt(x*x + 1.0)) ;
X return(retval) ;
X}
X#endif /*NEED_ASINH*/
X
X
X/* FUNCTION
X *
X * atan double precision arc tangent
X *
X * DESCRIPTION
X *
X * Returns double precision arc tangent of double precision
X * floating point argument.
X *
X * USAGE
X *
X * double atan(x)
X * double x ;
X *
X * REFERENCES
X *
X * Fortran 77 user's guide, Digital Equipment Corp. pp B-3
X *
X * Computer Approximations, J.F. Hart et al, John Wiley & Sons,
X * 1968, pp. 120-130.
X *
X * RESTRICTIONS
X *
X * The maximum relative error for the approximating polynomial
X * is 10**(-16.82). However, this assumes exact arithmetic
X * in the polynomial evaluation. Additional rounding and
X * truncation errors may occur as the argument is reduced
X * to the range over which the polynomial approximation
X * is valid, and as the polynomial is evaluated using
X * finite-precision arithmetic.
X *
X * INTERNALS
X *
X * Computes arctangent(x) from:
X *
X * (1) If x < 0 then negate x, perform steps
X * 2, 3, and 4, and negate the returned
X * result. This makes use of the identity
X * atan(-x) = -atan(x).
X *
X * (2) If argument x > 1 at this point then
X * test to be sure that x can be inverted
X * without underflowing. If not, reduce
X * x to largest possible number that can
X * be inverted, issue warning, and continue.
X * Perform steps 3 and 4 with arg = 1/x
X * and subtract returned result from
X * pi/2. This makes use of the identity
X * atan(x) = pi/2 - atan(1/x) for x>0.
X *
X * (3) At this point 0 <= x <= 1. If
X * x > tan(pi/12) then perform step 4
X * with arg = (x*sqrt(3)-1)/(sqrt(3)+x)
X * and add pi/6 to returned result. This
X * last transformation maps arguments
X * greater than tan(pi/12) to arguments
X * in range 0...tan(pi/12).
X *
X * (4) At this point the argument has been
X * transformed so that it lies in the
X * range -tan(pi/12)...tan(pi/12).
X * Since very small arguments may cause
X * underflow in the polynomial evaluation,
X * a final check is performed. If the
X * argument is less than the underflow
X * bound, the function returns x, since
X * atan(x) approaches asin(x) which
X * approaches x, as x goes to zero.
X *
X * (5) atan(x) is now computed by one of the
X * approximations given in the cited
X * reference (Hart). That is:
X *
X * atan(x) = x * SUM [ C[i] * x**(2*i) ]
X * over i = {0,1,...8}.
X *
X * Where:
X *
X * C[0] = .9999999999999999849899
X * C[1] = -.333333333333299308717
X * C[2] = .1999999999872944792
X * C[3] = -.142857141028255452
X * C[4] = .11111097898051048
X * C[5] = -.0909037114191074
X * C[6] = .0767936869066
X * C[7] = -.06483193510303
X * C[8] = .0443895157187
X * (coefficients from HART table #4945 pg 267)
X */
X
X#ifdef NEED_ATAN
X
X#define LAST_BOUND 0.2679491924311227074725 /* tan (PI/12) */
X
Xstatic double atan_coeffs[] = {
X .9999999999999999849899, /* P0 must be first */
X -.333333333333299308717,
X .1999999999872944792,
X -.142857141028255452,
X .11111097898051048,
X -.0909037114191074,
X .0767936869066,
X -.06483193510303,
X .0443895157187 /* Pn must be last */
X} ;
X
Xdouble
Xatan(x)
Xdouble x ;
X{
X register int order ;
X double t1, t2, xt2 ;
X auto double retval ;
X
X if (x < 0.0) retval = -(atan(-x)) ;
X else if (x > 1.0)
X {
X if (x < MAXDOUBLE && x > -MAXDOUBLE)
X retval = HALFPI - atan(1.0 / x) ;
X else
X {
X doerr("atan", "UNDERFLOW", EDOM) ;
X retval = 0.0 ;
X }
X }
X else if (x > LAST_BOUND)
X {
X t1 = x * SQRT3 - 1.0 ;
X t2 = SQRT3 + x ;
X retval = SIXTHPI + atan(t1 / t2) ;
X }
X else if (x < X16_UNDERFLOWS)
X {
X doerr("atan", "PLOSS", EDOM) ;
X retval = x ;
X }
X else
X {
X xt2 = x * x ;
X order = sizeof(atan_coeffs) / sizeof(double) ;
X order -= 1 ;
X retval = x * poly(order, atan_coeffs, xt2) ;
X }
X return(retval) ;
X}
X#endif /*NEED_ATAN*/
X
X
X/* FUNCTION
X *
X * atanh double precision hyperbolic arc tangent
X *
X * DESCRIPTION
X *
X * Returns double precision hyperbolic arc tangent of double precision
X * floating point number.
X *
X * USAGE
X *
X * double atanh(x)
X * double x ;
X *
X * RESTRICTIONS
X *
X * The range of the atanh function is -1.0 to +1.0 exclusive.
X * Certain pathological cases near 1.0 may cause overflow
X * in evaluation of 1+x / 1-x, depending upon machine exponent
X * range and mantissa precision.
X *
X * For precision information refer to documentation of the
X * other floating point library routines called.
X *
X * INTERNALS
X *
X * Computes atanh(x) from:
X *
X * 1. If x <= -1.0 or x >= 1.0
X * then report argument out of range and return 0.0
X *
X * 2. atanh(x) = 0.5 * log((1+x)/(1-x))
X */
X
X#ifdef NEED_ATANH
Xdouble
Xatanh(x)
Xdouble x ;
X{
X auto double retval ;
X
X if (x <= -1.0 || x >= 1.0)
X {
X doerr("atan", "DOMAIN", ERANGE) ;
X retval = 0.0 ;
X }
X else retval = 0.5 * log((1+x) / (1-x)) ;
X return(retval) ;
X}
X#endif /*NEED_ATANH*/
X
X
X/* FUNCTION
X *
X * cos - double precision cosine
X *
X * DESCRIPTION
X *
X * Returns double precision cosine of double precision
X * floating point argument.
X *
X * USAGE
X *
X * double cos(x)
X * double x ;
X *
X * REFERENCES
X *
X * Computer Approximations, J.F. Hart et al, John Wiley & Sons,
X * 1968, pp. 112-120.
X *
X * RESTRICTIONS
X *
X * The sin and cos routines are interactive in the sense that
X * in the process of reducing the argument to the range -PI/4
X * to PI/4, each may call the other. Ultimately one or the
X * other uses a polynomial approximation on the reduced
X * argument. The sin approximation has a maximum relative error
X * of 10**(-17.59) and the cos approximation has a maximum
X * relative error of 10**(-16.18).
X *
X * These error bounds assume exact arithmetic
X * in the polynomial evaluation. Additional rounding and
X * truncation errors may occur as the argument is reduced
X * to the range over which the polynomial approximation
X * is valid, and as the polynomial is evaluated using
X * finite-precision arithmetic.
X *
X * INTERNALS
X *
X * Computes cos(x) from:
X *
X * (1) Reduce argument x to range -PI to PI.
X *
X * (2) If x > PI/2 then call cos recursively
X * using relation cos(x) = -cos(x - PI).
X *
X * (3) If x < -PI/2 then call cos recursively
X * using relation cos(x) = -cos(x + PI).
X *
X * (4) If x > PI/4 then call sin using
X * relation cos(x) = sin(PI/2 - x).
X *
X * (5) If x < -PI/4 then call cos using
X * relation cos(x) = sin(PI/2 + x).
X *
X * (6) If x would cause underflow in approx
X * evaluation arithmetic then return
X * sqrt(1.0 - x**2).
X *
X * (7) By now x has been reduced to range
X * -PI/4 to PI/4 and the approximation
X * from HART pg. 119 can be used:
X *
X * cos(x) = ( p(y) / q(y) )
X * Where:
X *
X * y = x * (4/PI)
X *
X * p(y) = SUM [ Pj * (y**(2*j)) ]
X * over j = {0,1,2,3}
X *
X * q(y) = SUM [ Qj * (y**(2*j)) ]
X * over j = {0,1,2,3}
X *
X * P0 = 0.12905394659037374438571854e+7
X * P1 = -0.3745670391572320471032359e+6
X * P2 = 0.134323009865390842853673e+5
X * P3 = -0.112314508233409330923e+3
X * Q0 = 0.12905394659037373590295914e+7
X * Q1 = 0.234677731072458350524124e+5
X * Q2 = 0.2096951819672630628621e+3
X * Q3 = 1.0000...
X * (coefficients from HART table #3843 pg 244)
X *
X *
X * **** NOTE **** The range reduction relations used in
X * this routine depend on the final approximation being valid
X * over the negative argument range in addition to the positive
X * argument range. The particular approximation chosen from
X * HART satisfies this requirement, although not explicitly
X * stated in the text. This may not be true of other
X * approximations given in the reference.
X */
X
X#ifdef NEED_COS
X
Xstatic double cos_pcoeffs[] = {
X 0.12905394659037374438e7,
X -0.37456703915723204710e6,
X 0.13432300986539084285e5,
X -0.11231450823340933092e3
X} ;
X
Xstatic double cos_qcoeffs[] = {
X 0.12905394659037373590e7,
X 0.23467773107245835052e5,
X 0.20969518196726306286e3,
X 1.0
X} ;
X
Xdouble
Xcos(x)
Xdouble x ;
X{
X auto double y, yt2 ;
X auto double retval ;
X
X if (x < -PI || x > PI)
X {
X x = mod(x, TWOPI) ;
X if (x > PI) x = x - TWOPI ;
X else if (x < -PI) x = x + TWOPI ;
X }
X if (x > HALFPI) retval = -(cos(x - PI)) ;
X else if (x < -HALFPI) retval = -(cos(x + PI)) ;
X else if (x > FOURTHPI) retval = sin(HALFPI - x) ;
X else if (x < -FOURTHPI) retval = sin(HALFPI + x) ;
X else if (x < X6_UNDERFLOWS && x > -X6_UNDERFLOWS)
X retval = sqrt(1.0 - (x * x)) ;
X else
X {
X y = x / FOURTHPI ;
X yt2 = y * y ;
X retval = poly(3, cos_pcoeffs, yt2) / poly(3, cos_qcoeffs, yt2) ;
X }
X return(retval) ;
X}
X#endif /*NEED_COS*/
X
X
X/* FUNCTION
X *
X * cosh double precision hyperbolic cosine
X *
X * DESCRIPTION
X *
X * Returns double precision hyperbolic cosine of double precision
X * floating point number.
X *
X * USAGE
X *
X * double cosh(x)
X * double x ;
X *
X * REFERENCES
X *
X * Fortran IV plus user's guide, Digital Equipment Corp. pp B-4
X *
X * RESTRICTIONS
X *
X * Inputs greater than log(MAXDOUBLE) result in overflow.
X * Inputs less than log(MINDOUBLE) result in underflow.
X *
X * For precision information refer to documentation of the
X * floating point library routines called.
X *
X * INTERNALS
X *
X * Computes hyperbolic cosine from:
X *
X * cosh(X) = 0.5 * (exp(X) + exp(-X))
X */
X
X
X#ifdef NEED_COSH
Xdouble
Xcosh(x)
Xdouble x ;
X{
X auto double retval ;
X
X if (x > LOGE_MAXDOUBLE)
X {
X doerr("cosh", "OVERFLOW", ERANGE) ;
X retval = MAXDOUBLE ;
X }
X else if (x < LOGE_MINDOUBLE)
X {
X doerr("cosh", "UNDERFLOW", ERANGE) ;
X retval = MINDOUBLE ;
X }
X else
X {
X x = exp(x) ;
X retval = 0.5 * (x + 1.0 / x) ;
X }
X return(retval) ;
X}
X#endif /*NEED_COSH*/
X
X
X/* FUNCTION
X *
X * dabs double precision absolute value
X *
X * DESCRIPTION
X *
X * Returns absolute value of double precision number.
X *
X * USAGE
X *
X * double dabs(x)
X * double x ;
X */
X
X#ifdef NEED_EXP
Xdouble
Xdabs(x)
Xdouble x ;
X{
X if (x < 0.0) x = -x ;
X return(x) ;
X}
X#endif /*NEED_EXP*/
X
X
X/* FUNCTION
X *
X * exp double precision exponential
X *
X * DESCRIPTION
X *
X * Returns double precision exponential of double precision
X * floating point number.
X *
X * USAGE
X *
X * double exp(x)
X * double x ;
X *
X * REFERENCES
X *
X * Fortran IV plus user's guide, Digital Equipment Corp. pp B-3
X *
X * Computer Approximations, J.F. Hart et al, John Wiley & Sons,
X * 1968, pp. 96-104.
X *
X * RESTRICTIONS
X *
X * Inputs greater than log(MAXDOUBLE) result in overflow.
X * Inputs less than log(MINDOUBLE) result in underflow.
X *
X * The maximum relative error for the approximating polynomial
X * is 10**(-16.4). However, this assumes exact arithmetic
X * in the polynomial evaluation. Additional rounding and
X * truncation errors may occur as the argument is reduced
X * to the range over which the polynomial approximation
X * is valid, and as the polynomial is evaluated using
X * finite precision arithmetic.
X *
X *
X * INTERNALS
X *
X * Computes exponential from:
X *
X * exp(x) = 2**y * 2**z * 2**w
X *
X * Where:
X *
X * y = int ( x * log2(e) )
X *
X * v = 16 * frac ( x * log2(e))
X *
X * z = (1/16) * int (v)
X *
X * w = (1/16) * frac (v)
X *
X * Note that:
X *
X * 0 =< v < 16
X *
X * z = {0, 1/16, 2/16, ...15/16}
X *
X * 0 =< w < 1/16
X *
X * Then:
X *
X * 2**z is looked up in a table of 2**0, 2**1/16, ...
X *
X * 2**w is computed from an approximation:
X *
X * 2**w = (A + B) / (A - B)
X *
X * A = C + (D * w * w)
X *
X * B = w * (E + (F * w * w))
X *
X * C = 20.8137711965230361973
X *
X * D = 1.0
X *
X * E = 7.2135034108448192083
X *
X * F = 0.057761135831801928
X *
X * Coefficients are from HART, table #1121, pg 206.
X *
X * Effective multiplication by 2**y is done by a
X * floating point scale with y as scale argument.
X */
X
X#ifdef NEED_EXP
X
X#define C 20.8137711965230361973 /* Polynomial approx coeff. */
X#define D 1.0 /* Polynomial approx coeff. */
X#define E 7.2135034108448192083 /* Polynomial approx coeff. */
X#define F 0.057761135831801928 /* Polynomial approx coeff. */
X
X/************************************************************************
X * *
X * This table is fixed in size and reasonably hardware *
X * independent. The given constants were generated on a *
X * DECSYSTEM 20 using double precision FORTRAN. *
X * *
X ************************************************************************
X */
X
Xstatic double fpof2tbl[] = {
X 1.00000000000000000000, /* 2 ** 0/16 */
X 1.04427378242741384020, /* 2 ** 1/16 */
X 1.09050773266525765930, /* 2 ** 2/16 */
X 1.13878863475669165390, /* 2 ** 3/16 */
X 1.18920711500272106640, /* 2 ** 4/16 */
X 1.24185781207348404890, /* 2 ** 5/16 */
X 1.29683955465100966610, /* 2 ** 6/16 */
X 1.35425554693689272850, /* 2 ** 7/16 */
X 1.41421356237309504880, /* 2 ** 8/16 */
X 1.47682614593949931110, /* 2 ** 9/16 */
X 1.54221082540794082350, /* 2 ** 10/16 */
X 1.61049033194925430820, /* 2 ** 11/16 */
X 1.68179283050742908600, /* 2 ** 12/16 */
X 1.75625216037329948340, /* 2 ** 13/16 */
X 1.83400808640934246360, /* 2 ** 14/16 */
X 1.91520656139714729380 /* 2 ** 15/16 */
X} ;
X
Xdouble
Xexp(x)
Xdouble x ;
X{
X register int ival, y ;
X auto double a, b, rtnval, t, temp, v, w, wpof2, zpof2 ;
X auto double retval ;
X
X if (x > LOGE_MAXDOUBLE)
X {
X doerr("exp", "OVERFLOW", ERANGE) ;
X retval = MAXDOUBLE ;
X }
X else if (x <= LOGE_MINDOUBLE)
X {
X doerr("exp", "UNDERFLOW", ERANGE) ;
X retval = 0.0 ;
X }
X else
X {
X t = x * LOG2E ;
X (void) modf(t, &temp) ;
X y = temp ;
X v = 16.0 * modf(t, &temp) ;
X (void) modf(dabs(v), &temp) ;
X ival = temp ;
X if (x < 0.0) zpof2 = 1.0 / fpof2tbl[ival] ;
X else zpof2 = fpof2tbl[ival] ;
X w = modf(v, &temp) / 16.0 ;
X a = C + (D * w * w) ;
X b = w * (E + (F * w * w)) ;
X wpof2 = (a + b) / (a - b) ;
X retval = ldexp((wpof2 * zpof2), y) ;
X }
X return(retval) ;
X}
X#endif /*NEED_EXP*/
X
X
X/* FUNCTION
X *
X * log double precision natural log
X *
X * DESCRIPTION
X *
X * Returns double precision natural log of double precision
X * floating point argument.
X *
X * USAGE
X *
X * double log(x)
X * double x ;
X *
X * REFERENCES
X *
X * Computer Approximations, J.F. Hart et al, John Wiley & Sons,
X * 1968, pp. 105-111.
X *
X * RESTRICTIONS
X *
X * The absolute error in the approximating polynomial is
X * 10**(-19.38). Note that contrary to DEC's assertion
X * in the F4P user's guide, the error is absolute and not
X * relative.
X *
X * This error bound assumes exact arithmetic
X * in the polynomial evaluation. Additional rounding and
X * truncation errors may occur as the argument is reduced
X * to the range over which the polynomial approximation
X * is valid, and as the polynomial is evaluated using
X * finite-precision arithmetic.
X *
X * INTERNALS
X *
X * Computes log(X) from:
X *
X * (1) If argument is zero then flag an error
X * and return minus infinity (or rather its
X * machine representation).
X *
X * (2) If argument is negative then flag an
X * error and return minus infinity.
X *
X * (3) Given that x = m * 2**k then extract
X * mantissa m and exponent k.
X *
X * (4) Transform m which is in range [0.5,1.0]
X * to range [1/sqrt(2), sqrt(2)] by:
X *
X * s = m * sqrt(2)
X *
X * (5) Compute z = (s - 1) / (s + 1)
X *
X * (6) Now use the approximation from HART
X * page 111 to find log(s):
X *
X * log(s) = z * ( P(z**2) / Q(z**2) )
X *
X * Where:
X *
X * P(z**2) = SUM [ Pj * (z**(2*j)) ]
X * over j = {0,1,2,3}
X *
X * Q(z**2) = SUM [ Qj * (z**(2*j)) ]
X * over j = {0,1,2,3}
X *
X * P0 = -0.240139179559210509868484e2
X * P1 = 0.30957292821537650062264e2
X * P2 = -0.96376909336868659324e1
X * P3 = 0.4210873712179797145
X * Q0 = -0.120069589779605254717525e2
X * Q1 = 0.19480966070088973051623e2
X * Q2 = -0.89111090279378312337e1
X * Q3 = 1.0000
X *
X * (coefficients from HART table #2705 pg 229)
X *
X * (7) Finally, compute log(x) from:
X *
X * log(x) = (k * log(2)) - log(sqrt(2)) + log(s)
X */
X
X#ifdef NEED_LOG
X
Xstatic double log_pcoeffs[] = {
X -0.24013917955921050986e2,
X 0.30957292821537650062e2,
X -0.96376909336868659324e1,
X 0.4210873712179797145
X} ;
X
Xstatic double log_qcoeffs[] = {
X -0.12006958977960525471e2,
X 0.19480966070088973051e2,
X -0.89111090279378312337e1,
X 1.0000
X} ;
X
Xdouble
Xlog(x)
Xdouble x ;
X{
X auto int k ;
X auto double pqofz, s, z, zt2 ;
X auto double retval ;
X
X if (x == 0.0)
X {
X doerr("log", "SINGULARITY", EDOM) ;
X retval = -MAXDOUBLE ;
X }
X else if (x < 0.0)
X {
X doerr("log", "DOMAIN", EDOM) ;
X retval = -MAXDOUBLE ;
X }
X else
X {
X s = SQRT2 * frexp(x, &k) ;
X z = (s - 1.0) / (s + 1.0) ;
X zt2 = z * z ;
X pqofz = z * (poly(3, log_pcoeffs, zt2) / poly(3, log_qcoeffs, zt2)) ;
X x = k * LN2 ;
X x -= LNSQRT2 ;
X x += pqofz ;
X retval = x ;
X }
X return(retval) ;
X}
X#endif /*NEED_LOG*/
X
X
X/* FUNCTION
X *
X * log10 double precision common log
X *
X * DESCRIPTION
X *
X * Returns double precision common log of double precision
X * floating point argument.
X *
X * USAGE
X *
X * double log10(x)
X * double x ;
X *
X * REFERENCES
X *
X * PDP-11 Fortran IV-plus users guide, Digital Equip. Corp.,
X * 1975, pp. B-3.
X *
X * RESTRICTIONS
X *
X * For precision information refer to documentation of the
X * floating point library routines called.
X *
X * INTERNALS
X *
X * Computes log10(x) from:
X *
X * log10(x) = log10(e) * log(x)
X */
X
X#ifdef NEED_LOG10
Xdouble
Xlog10(x)
Xdouble x ;
X{
X x = LOG10E * log(x) ;
X return(x) ;
X}
X#endif /*NEED_LOG10*/
X
X
X/* FUNCTION
X *
X * mod double precision modulo
X *
X * DESCRIPTION
X *
X * Returns double precision modulo of two double
X * precision arguments.
X *
X * USAGE
X *
X * double mod(value, base)
X * double value ;
X * double base ;
X */
X
X#ifdef NEED_COS || NEED_SIN
Xdouble mod(value, base)
Xdouble value ;
Xdouble base ;
X{
X auto double intpart ;
X
X value /= base ;
X value = modf(value, &intpart) ;
X value *= base ;
X return(value) ;
X}
X#endif /*NEED_COS || NEED_SIN*/
X
X
X/* FUNCTION
X *
X * poly double precision polynomial evaluation
X *
X * DESCRIPTION
X *
X * Evaluates a polynomial and returns double precision
X * result. Is passed a the order of the polynomial,
X * a pointer to an array of double precision polynomial
X * coefficients (in ascending order), and the independent
X * variable.
X *
X * USAGE
X *
X * double poly(order, coeffs, x)
X * int order ;
X * double *coeffs ;
X * double x ;
X *
X * INTERNALS
X *
X * Evalates the polynomial using recursion and the form:
X *
X * P(x) = P0 + x(P1 + x(P2 +...x(Pn)))
X */
X
X#ifdef NEED_ATAN || NEED_COS || NEED_LOG || NEED_SIN
Xdouble
Xpoly(order, coeffs, x)
Xregister int order ;
Xdouble *coeffs ;
Xdouble x ;
X{
X auto double curr_coeff, rtn_value ;
X
X if (order <= 0) rtn_value = *coeffs ;
X else
X {
X curr_coeff = *coeffs ; /* Bug in Unisoft's compiler. Does not */
X coeffs++ ; /* generate good code for *coeffs++ */
X rtn_value = curr_coeff + x * poly(--order, coeffs, x) ;
X }
X return(rtn_value) ;
X}
X#endif /*NEED_ATAN || NEED_COS || NEED_LOG || NEED_SIN*/
SHAR_EOF
echo "End of part 2"
echo "File mathlib.c is continued in part 3"
echo "3" > s2_seq_.tmp
exit 0richb@sunaus.oz (Rich Burridge) (02/02/90)
---- Cut Here and unpack ----
#!/bin/sh
# this is part 3 of a multipart archive
# do not concatenate these parts, unpack them in order with /bin/sh
# file mathlib.c continued
#
CurArch=3
if test ! -r s2_seq_.tmp
then echo "Please unpack part 1 first!"
exit 1; fi
( read Scheck
if test "$Scheck" != $CurArch
then echo "Please unpack part $Scheck next!"
exit 1;
else exit 0; fi
) < s2_seq_.tmp || exit 1
echo "x - Continuing file mathlib.c"
sed 's/^X//' << 'SHAR_EOF' >> mathlib.c
X
X
X/* FUNCTION
X *
X * sin double precision sine
X *
X * DESCRIPTION
X *
X * Returns double precision sine of double precision
X * floating point argument.
X *
X * USAGE
X *
X * double sin(x)
X * double x ;
X *
X * REFERENCES
X *
X * Computer Approximations, J.F. Hart et al, John Wiley & Sons,
X * 1968, pp. 112-120.
X *
X * RESTRICTIONS
X *
X * The sin and cos routines are interactive in the sense that
X * in the process of reducing the argument to the range -PI/4
X * to PI/4, each may call the other. Ultimately one or the
X * other uses a polynomial approximation on the reduced
X * argument. The sin approximation has a maximum relative error
X * of 10**(-17.59) and the cos approximation has a maximum
X * relative error of 10**(-16.18).
X *
X * These error bounds assume exact arithmetic
X * in the polynomial evaluation. Additional rounding and
X * truncation errors may occur as the argument is reduced
X * to the range over which the polynomial approximation
X * is valid, and as the polynomial is evaluated using
X * finite-precision arithmetic.
X *
X * INTERNALS
X *
X * Computes sin(x) from:
X *
X * (1) Reduce argument x to range -PI to PI.
X *
X * (2) If x > PI/2 then call sin recursively
X * using relation sin(x) = -sin(x - PI).
X *
X * (3) If x < -PI/2 then call sin recursively
X * using relation sin(x) = -sin(x + PI).
X *
X * (4) If x > PI/4 then call cos using
X * relation sin(x) = cos(PI/2 - x).
X *
X * (5) If x < -PI/4 then call cos using
X * relation sin(x) = -cos(PI/2 + x).
X *
X * (6) If x is small enough that polynomial
X * evaluation would cause underflow
X * then return x, since sin(x)
X * approaches x as x approaches zero.
X *
X * (7) By now x has been reduced to range
X * -PI/4 to PI/4 and the approximation
X * from HART pg. 118 can be used:
X *
X * sin(x) = y * ( p(y) / q(y) )
X * Where:
X *
X * y = x * (4/PI)
X *
X * p(y) = SUM [ Pj * (y**(2*j)) ]
X * over j = {0,1,2,3}
X *
X * q(y) = SUM [ Qj * (y**(2*j)) ]
X * over j = {0,1,2,3}
X *
X * P0 = 0.206643433369958582409167054e+7
X * P1 = -0.18160398797407332550219213e+6
X * P2 = 0.359993069496361883172836e+4
X * P3 = -0.2010748329458861571949e+2
X * Q0 = 0.263106591026476989637710307e+7
X * Q1 = 0.3927024277464900030883986e+5
X * Q2 = 0.27811919481083844087953e+3
X * Q3 = 1.0000...
X * (coefficients from HART table #3063 pg 234)
X *
X *
X * **** NOTE **** The range reduction relations used in
X * this routine depend on the final approximation being valid
X * over the negative argument range in addition to the positive
X * argument range. The particular approximation chosen from
X * HART satisfies this requirement, although not explicitly
X * stated in the text. This may not be true of other
X * approximations given in the reference.
X */
X
X#ifdef NEED_SIN
X
Xstatic double sin_pcoeffs[] = {
X 0.20664343336995858240e7,
X -0.18160398797407332550e6,
X 0.35999306949636188317e4,
X -0.20107483294588615719e2
X} ;
X
Xstatic double sin_qcoeffs[] = {
X 0.26310659102647698963e7,
X 0.39270242774649000308e5,
X 0.27811919481083844087e3,
X 1.0
X} ;
X
Xdouble
Xsin(x)
Xdouble x ;
X{
X double y, yt2 ;
X auto double retval ;
X
X if (x < -PI || x > PI)
X {
X x = mod(x, TWOPI) ;
X if (x > PI) x = x - TWOPI ;
X else if (x < -PI) x = x + TWOPI ;
X }
X if (x > HALFPI) retval = -(sin(x - PI)) ;
X else if (x < -HALFPI) retval = -(sin(x + PI)) ;
X else if (x > FOURTHPI) retval = cos(HALFPI - x) ;
X else if (x < -FOURTHPI) retval = -(cos(HALFPI + x)) ;
X else if (x < X6_UNDERFLOWS && x > -X6_UNDERFLOWS) retval = x ;
X else
X {
X y = x / FOURTHPI ;
X yt2 = y * y ;
X retval = y * (poly(3, sin_pcoeffs, yt2) / poly(3, sin_qcoeffs, yt2)) ;
X }
X return(retval) ;
X}
X#endif /*NEED_SIN*/
X
X
X/* FUNCTION
X *
X * sinh double precision hyperbolic sine
X *
X * DESCRIPTION
X *
X * Returns double precision hyperbolic sine of double precision
X * floating point number.
X *
X * USAGE
X *
X * double sinh(x)
X * double x ;
X *
X * REFERENCES
X *
X * Fortran IV plus user's guide, Digital Equipment Corp. pp B-5
X *
X * RESTRICTIONS
X *
X * Inputs greater than ln(MAXDOUBLE) result in overflow.
X * Inputs less than ln(MINDOUBLE) result in underflow.
X *
X * For precision information refer to documentation of the
X * floating point library routines called.
X *
X * INTERNALS
X *
X * Computes hyperbolic sine from:
X *
X * sinh(x) = 0.5 * (exp(x) - 1.0/exp(x))
X */
X
X#ifdef NEED_SINH
Xdouble
Xsinh(x)
Xdouble x ;
X{
X auto double retval ;
X
X if (x > LOGE_MAXDOUBLE)
X {
X doerr("sinh", "OVERFLOW", ERANGE) ;
X retval = MAXDOUBLE ;
X }
X else if (x < LOGE_MINDOUBLE)
X {
X doerr("sinh", "UNDERFLOW", ERANGE) ;
X retval = MINDOUBLE ;
X }
X else
X {
X x = exp(x) ;
X retval = 0.5 * (x - (1.0 / x)) ;
X }
X return(retval) ;
X}
X#endif /*NEED_SINH*/
X
X
X/* FUNCTION
X *
X * sqrt double precision square root
X *
X * DESCRIPTION
X *
X * Returns double precision square root of double precision
X * floating point argument.
X *
X * USAGE
X *
X * double sqrt(x)
X * double x ;
X *
X * REFERENCES
X *
X * Fortran IV-PLUS user's guide, Digital Equipment Corp. pp B-7.
X *
X * Computer Approximations, J.F. Hart et al, John Wiley & Sons,
X * 1968, pp. 89-96.
X *
X * RESTRICTIONS
X *
X * The relative error is 10**(-30.1) after three applications
X * of Heron's iteration for the square root.
X *
X * However, this assumes exact arithmetic in the iterations
X * and initial approximation. Additional errors may occur
X * due to truncation, rounding, or machine precision limits.
X *
X * INTERNALS
X *
X * Computes square root by:
X *
X * (1) Range reduction of argument to [0.5,1.0]
X * by application of identity:
X *
X * sqrt(x) = 2**(k/2) * sqrt(x * 2**(-k))
X *
X * k is the exponent when x is written as
X * a mantissa times a power of 2 (m * 2**k).
X * It is assumed that the mantissa is
X * already normalized (0.5 =< m < 1.0).
X *
X * (2) An approximation to sqrt(m) is obtained
X * from:
X *
X * u = sqrt(m) = (P0 + P1*m) / (Q0 + Q1*m)
X *
X * P0 = 0.594604482
X * P1 = 2.54164041
X * Q0 = 2.13725758
X * Q1 = 1.0
X *
X * (coefficients from HART table #350 pg 193)
X *
X * (3) Three applications of Heron's iteration are
X * performed using:
X *
X * y[n+1] = 0.5 * (y[n] + (m/y[n]))
X *
X * where y[0] = u = approx sqrt(m)
X *
X * (4) If the value of k was odd then y is either
X * multiplied by the square root of two or
X * divided by the square root of two for k positive
X * or negative respectively. This rescales y
X * by multiplying by 2**frac(k/2).
X *
X * (5) Finally, y is rescaled by int(k/2) which
X * is equivalent to multiplication by 2**int(k/2).
X *
X * The result of steps 4 and 5 is that the value
X * of y between 0.5 and 1.0 has been rescaled by
X * 2**(k/2) which removes the original rescaling
X * done prior to finding the mantissa square root.
X *
X * NOTES
X *
X * The Convergent Technologies compiler optimizes division
X * by powers of two to "arithmetic shift right" instructions.
X * This is ok for positive integers but gives different
X * results than other C compilers for negative integers.
X * For example, (-1)/2 is -1 on the CT box, 0 on every other
X * machine I tried.
X */
X
X#ifdef NEED_SQRT
X
X#define ITERATIONS 3 /* Number of iterations */
X#define P0 0.594604482 /* Approximation coeff */
X#define P1 2.54164041 /* Approximation coeff */
X#define Q0 2.13725758 /* Approximation coeff */
X#define Q1 1.0 /* Approximation coeff */
X
Xdouble
Xsqrt(x)
Xdouble x ;
X{
X register int bugfix, count, kmod2 ;
X auto int exponent, k ;
X auto double m, u, y ;
X auto double retval ;
X
X if (x == 0.0) retval = 0.0 ;
X else if (x < 0.0)
X {
X doerr("sqrt", "DOMAIN", EDOM) ;
X retval = 0.0 ;
X }
X else
X {
X m = frexp(x, &k) ;
X u = (P0 + (P1 * m)) / (Q0 + (Q1 * m)) ;
X for (count = 0, y = u; count < ITERATIONS; count++)
X y = 0.5 * (y + (m / y)) ;
X if ((kmod2 = (k % 2)) < 0) y /= SQRT2 ;
X else if (kmod2 > 0) y *= SQRT2 ;
X bugfix = 2 ;
X retval = ldexp(y, k / bugfix) ;
X }
X return(retval) ;
X}
X#endif /*NEED_SQRT*/
X
X
X/* FUNCTION
X *
X * tan Double precision tangent
X *
X * DESCRIPTION
X *
X * Returns tangent of double precision floating point number.
X *
X * USAGE
X *
X * double tan(x)
X * double x ;
X *
X * INTERNALS
X *
X * Computes the tangent from tan(x) = sin(x) / cos(x).
X *
X * If cos(x) = 0 and sin(x) >= 0, then returns largest
X * floating point number on host machine.
X *
X * If cos(x) = 0 and sin(x) < 0, then returns smallest
X * floating point number on host machine.
X *
X * REFERENCES
X *
X * Fortran IV plus user's guide, Digital Equipment Corp. pp. B-8
X */
X
X#ifdef NEED_TAN
Xdouble
Xtan(x)
Xdouble x ;
X{
X double cosx, sinx ;
X auto double retval ;
X
X sinx = sin(x) ;
X cosx = cos(x) ;
X if (cosx == 0.0)
X {
X doerr("tan", "OVERFLOW", ERANGE) ;
X if (sinx >= 0.0) retval = MAXDOUBLE ;
X else retval = -MAXDOUBLE ;
X }
X else retval = sinx / cosx ;
X return(retval) ;
X}
X#endif /*NEED_TAN*/
X
X
X/* FUNCTION
X *
X * tanh double precision hyperbolic tangent
X *
X * DESCRIPTION
X *
X * Returns double precision hyperbolic tangent of double precision
X * floating point number.
X *
X * USAGE
X *
X * double tanh(x)
X * double x ;
X *
X * REFERENCES
X *
X * Fortran IV plus user's guide, Digital Equipment Corp. pp B-5
X *
X * RESTRICTIONS
X *
X * For precision information refer to documentation of the
X * floating point library routines called.
X *
X * INTERNALS
X *
X * Computes hyperbolic tangent from:
X *
X * tanh(x) = 1.0 for x > TANH_MAXARG
X * = -1.0 for x < -TANH_MAXARG
X * = sinh(x) / cosh(x) otherwise
X */
X
X#ifdef NEED_TANH
Xdouble
Xtanh(x)
Xdouble x ;
X{
X auto double retval ;
X register int positive ;
X
X if (x > TANH_MAXARG || x < -TANH_MAXARG)
X {
X if (x > 0.0) positive = 1 ;
X else positive = 0 ;
X doerr("tanh", "PLOSS", ERANGE) ;
X if (positive) retval = 1.0 ;
X else retval = -1.0 ;
X }
X else retval = sinh(x) / cosh(x) ;
X return(retval) ;
X}
X#endif /*NEED_TANH*/
X
X
X/*
X * Copyright (c) 1985 Regents of the University of California.
X * All rights reserved.
X *
X * Redistribution and use in source and binary forms are permitted
X * provided that the above copyright notice and this paragraph are
X * duplicated in all such forms and that any documentation,
X * advertising materials, and other materials related to such
X * distribution and use acknowledge that the software was developed
X * by the University of California, Berkeley. The name of the
X * University may not be used to endorse or promote products derived
X * from this software without specific prior written permission.
X * THIS SOFTWARE IS PROVIDED ``AS IS'' AND WITHOUT ANY EXPRESS OR
X * IMPLIED WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED
X * WARRANTIES OF MERCHANTIBILITY AND FITNESS FOR A PARTICULAR PURPOSE.
X *
X * All recipients should regard themselves as participants in an ongoing
X * research project and hence should feel obligated to report their
X * experiences (good or bad) with these elementary function codes, using
X * the sendbug(8) program, to the authors.
X */
X
X#ifdef NEED_POW
X/* POW(X,Y)
X * RETURN X**Y
X * DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS)
X * CODED IN C BY K.C. NG, 1/8/85;
X * REVISED BY K.C. NG on 7/10/85.
X *
X * Required system supported functions:
X * scalb(x,n)
X * logb(x)
X * copysign(x,y)
X * finite(x)
X * drem(x,y)
X *
X * Required kernel functions:
X * exp__E(a, c) ...return exp(a+c) - 1 - a*a/2
X * log__L(x) ...return (log(1+x) - 2s)/s, s=x/(2+x)
X * pow_p(x, y) ...return +(anything)**(finite non zero)
X *
X * Method
X * 1. Compute and return log(x) in three pieces:
X * log(x) = n*ln2 + hi + lo,
X * where n is an integer.
X * 2. Perform y * log(x) by simulating muti-precision arithmetic and
X * return the answer in three pieces:
X * y * log(x) = m * ln2 + hi + lo,
X * where m is an integer.
X * 3. Return x ** y = exp(y * log(x))
X * = 2^m * (exp(hi + lo)).
X *
X * Special cases:
X * (anything) ** 0 is 1 ;
X * (anything) ** 1 is itself;
X * (anything) ** NaN is NaN;
X * NaN ** (anything except 0) is NaN;
X * +-(anything > 1) ** +INF is +INF;
X * +-(anything > 1) ** -INF is +0;
X * +-(anything < 1) ** +INF is +0;
X * +-(anything < 1) ** -INF is +INF;
X * +-1 ** +-INF is NaN and signal INVALID;
X * +0 ** +(anything except 0, NaN) is +0;
X * -0 ** +(anything except 0, NaN, odd integer) is +0;
X * +0 ** -(anything except 0, NaN) is +INF and signal DIV-BY-ZERO;
X * -0 ** -(anything except 0, NaN, odd integer) is +INF with signal;
X * -0 ** (odd integer) = -( +0 ** (odd integer) );
X * +INF ** +(anything except 0,NaN) is +INF;
X * +INF ** -(anything except 0,NaN) is +0;
X * -INF ** (odd integer) = -( +INF ** (odd integer) );
X * -INF ** (even integer) = ( +INF ** (even integer) );
X * -INF ** -(anything except integer,NaN) is NaN with signal;
X * -(x=anything) ** (k=integer) is (-1)**k * (x ** k);
X * -(anything except 0) ** (non-integer) is NaN with signal;
X *
X * Accuracy:
X * pow(x, y) returns x**y nearly rounded. In particular, on a SUN, a VAX,
X * and a Zilog Z8000,
X * pow(integer, integer)
X * always returns the correct integer provided it is representable.
X * In a test run with 100,000 random arguments with 0 < x, y < 20.0
X * on a VAX, the maximum observed error was 1.79 ulps (units in the
X * last place).
X *
X * Constants :
X * The hexadecimal values are the intended ones for the following constants.
X * The decimal values may be used, provided that the compiler will convert
X * from decimal to binary accurately enough to produce the hexadecimal values
X * shown.
X */
X
X#if defined(vax) || defined(tahoe) /* VAX D format */
X#include <errno.h>
Xextern double infnan() ;
X#ifdef vax
X#define _0x(A,B) 0x/**/A/**/B
X#else /* vax */
X#define _0x(A,B) 0x/**/B/**/A
X#endif /* vax */
X/* static double */
X/* ln2hi = 6.9314718055829871446E-1 , Hex 2^ 0 * .B17217F7D00000 */
X/* ln2lo = 1.6465949582897081279E-12 , Hex 2^-39 * .E7BCD5E4F1D9CC */
X/* invln2 = 1.4426950408889634148E0 , Hex 2^ 1 * .B8AA3B295C17F1 */
X/* sqrt2 = 1.4142135623730950622E0 ; Hex 2^ 1 * .B504F333F9DE65 */
Xstatic long ln2hix[] = { _0x(7217,4031), _0x(0000,f7d0)} ;
Xstatic long ln2lox[] = { _0x(bcd5,2ce7), _0x(d9cc,e4f1)} ;
Xstatic long invln2x[] = { _0x(aa3b,40b8), _0x(17f1,295c)} ;
Xstatic long sqrt2x[] = { _0x(04f3,40b5), _0x(de65,33f9)} ;
X#define ln2hi (*(double*) ln2hix)
X#define ln2lo (*(double*) ln2lox)
X#define invln2 (*(double*) invln2x)
X#define sqrt2 (*(double*) sqrt2x)
X#else /* defined(vax)||defined(tahoe) */
Xstatic double
Xln2hi = 6.9314718036912381649E-1 , /* Hex 2^ -1 * 1.62E42FEE00000 */
Xln2lo = 1.9082149292705877000E-10 , /* Hex 2^-33 * 1.A39EF35793C76 */
Xinvln2 = 1.4426950408889633870E0 , /* Hex 2^ 0 * 1.71547652B82FE */
Xsqrt2 = 1.4142135623730951455E0 ; /* Hex 2^ 0 * 1.6A09E667F3BCD */
X#endif /* defined(vax) || defined(tahoe) */
X
Xstatic double zero = 0.0 ;
X
Xdouble
Xpow(x, y)
Xdouble x, y ;
X{
X double drem(), pow_p(), copysign(), t ;
X int finite() ;
X
X if (y == 0.0) return(1.0) ;
X#if !defined(vax) && !defined(tahoe)
X else if (y == 1.0 || x != x) return(x) ; /* if x is NaN or y = 1 */
X#else
X else if (y == 1.0) return(x) ; /* if y = 1 */
X#endif /* !defined(vax) && !defined(tahoe) */
X
X#if !defined(vax) && !defined(tahoe)
X else if (y != y) return(y) ; /* if y is NaN */
X#endif /* !defined(vax) && !defined(tahoe) */
X else if (!finite(y)) /* if y is INF */
X if ((t = copysign(x, 1.0)) == 1.0)
X return(0.0 / zero) ;
X else if (t > 1.0) return((y > 0.0) ? y : 0.0) ;
X else return((y < 0.0) ? -y : 0.0) ;
X else if (y == 2.0) return(x * x) ;
X else if (y == -1.0) return(1.0 / x) ;
X
X/* sign(x) = 1 */
X
X else if (copysign(1.0, x) == 1.0) return(pow_p(x, y)) ;
X
X/* sign(x)= -1 */
X/* if y is an even integer */
X
X else if ((t = drem(y, 2.0)) == 0.0) return(pow_p(-x, y)) ;
X
X/* if y is an odd integer */
X
X else if (copysign(t, 1.0) == 1.0) return(-pow_p(-x, y)) ;
X
X/* Henceforth y is not an integer */
X
X else if (x == 0.0) return((y > 0.0) ? -x : 1.0 / (-x)) ; /* x is -0 */
X else /* return NaN */
X {
X#if defined(vax) || defined(tahoe)
X return(infnan(EDOM)) ; /* NaN */
X#else /* defined(vax) || defined(tahoe) */
X return(0.0 / zero) ;
X#endif /* defined(vax) || defined(tahoe) */
X }
X}
X
X
X/* pow_p(x,y) return x**y for x with sign = 1 and finite y */
X
Xstatic double
Xpow_p(x, y)
Xdouble x, y ;
X{
X double logb(), scalb(), copysign(), log__L(), exp__E() ;
X double c, s, t, z, tx, ty ;
X
X#ifdef tahoe
X double tahoe_tmp ;
X#endif /* tahoe */
X float sx, sy ;
X long k = 0 ;
X int n, m ;
X
X if (x == 0.0 || !finite(x)) /* if x is +INF or +0 */
X {
X
X#if defined(vax) || defined(tahoe)
X return((y > 0.0) ? x : infnan(ERANGE)) ; /* if y < 0.0, return +INF */
X#else /* defined(vax) || defined(tahoe) */
X return((y > 0.0) ? x : 1.0 / x) ;
X#endif /* defined(vax) || defined(tahoe) */
X }
X if (x == 1.0) return(x) ; /* if x = 1.0, return 1 since y is finite */
X
X/* reduce x to z in [sqrt(1/2)-1, sqrt(2)-1] */
X
X z = scalb(x, -(n = logb(x))) ;
X
X#if !defined(vax) && !defined(tahoe) /* IEEE double; subnormal number */
X if (n <= -1022)
X {
X n += (m = logb(z)) ;
X z = scalb(z, -m) ;
X }
X#endif /* !defined(vax) && !defined(tahoe) */
X
X if (z >= sqrt2)
X {
X n += 1 ;
X z *= 0.5 ;
X }
X z -= 1.0 ;
X
X/* log(x) = nlog2 + log(1 + z) ~ nlog2 + t + tx */
X
X s = z / (2.0 + z) ;
X c = z * z * 0.5 ;
X tx = s * (c + log__L(s * s)) ;
X t = z - (c - tx) ;
X tx += (z - t) - c ;
X
X/* if y * log(x) is neither too big nor too small */
X
X if ((s = logb(y) + logb(n + t)) < 12.0)
X if (s > -60.0)
X {
X
X/* compute y * log(x) ~ mlog2 + t + c */
X
X s = y * (n + invln2 * t) ;
X m = s + copysign(0.5, s) ; /* m := nint(y * log(x)) */
X k = y ;
X if ((double) k == y) /* if y is an integer */
X {
X k = m - k * n ;
X sx = t ;
X tx += (t - sx) ;
X }
X else /* if y is not an integer */
X {
X k = m ;
X tx += n * ln2lo ;
X sx = (c = n * ln2hi) + t ;
X tx += (c - sx) + t ;
X }
X
X/* end of checking whether k == y */
X
X sy = y ;
X ty = y - sy ; /* y ~ sy + ty */
X
X#ifdef tahoe
X s = (tahoe_tmp = sx) * sy - k * ln2hi ;
X#else /* tahoe */
X s = (double) sx * sy - k * ln2hi ; /* (sy + ty) * (sx + tx) - kln2 */
X#endif /* tahoe */
X
X z = (tx * ty - k * ln2lo) ;
X tx = tx * sy ; ty = sx * ty ;
X t = ty + z ; t += tx ; t += s ;
X c = -((((t-s)-tx)-ty)-z) ;
X
X/* return exp(y * log(x)) */
X
X t += exp__E(t,c) ;
X return(scalb(1.0 + t, m)) ;
X }
X
X/* end of if log(y * log(x)) > -60.0 */
X
X else /* exp(+- tiny) = 1 with inexact flag */
X {
X ln2hi + ln2lo ;
X return(1.0) ;
X }
X else if (copysign(1.0, y) * (n + invln2 * t) < 0.0)
X return(scalb(1.0, -5000)) ; /* exp(-(big#)) underflows to zero */
X else
X return(scalb(1.0, 5000)) ; /* exp(+(big#)) overflows to INF */
X}
X
X
X/* Some IEEE standard 754 recommended functions and remainder and sqrt for
X * supporting the C elementary functions.
X ******************************************************************************
X * WARNING:
X * These codes are developed (in double) to support the C elementary
X * functions temporarily. They are not universal, and some of them are very
X * slow (in particular, drem and sqrt is extremely inefficient). Each
X * computer system should have its implementation of these functions using
X * its own assembler.
X ******************************************************************************
X *
X * IEEE 754 required operations:
X * drem(x, p)
X * returns x REM y = x - [x/y]*y , where [x/y] is the integer
X * nearest x/y; in half way case, choose the even one.
X * sqrt(x)
X * returns the square root of x correctly rounded according to
X * the rounding mod.
X *
X * IEEE 754 recommended functions:
X * (a) copysign(x, y)
X * returns x with the sign of y.
X * (b) scalb(x, N)
X * returns x * (2 ** N), for integer values N.
X * (c) logb(x)
X * returns the unbiased exponent of x, a signed integer in
X * double precision, except that logb(0) is -INF, logb(INF)
X * is +INF, and logb(NAN) is that NAN.
X * (d) finite(x)
X * returns the value TRUE if -INF < x < +INF and returns
X * FALSE otherwise.
X *
X * CODED IN C BY K.C. NG, 11/25/84;
X * REVISED BY K.C. NG on 1/22/85, 2/13/85, 3/24/85.
X */
X
X#if defined(vax) || defined(tahoe) /* VAX D format */
X static unsigned short msign = 0x7fff, mexp = 0x7f80 ;
X static short prep1 = 57, gap = 7, bias = 129 ;
X static double novf = 1.7E38, nunf = 3.0E-39 ;
X#else /* defined(vax) || defined(tahoe) */
X static unsigned short msign = 0x7fff, mexp = 0x7ff0 ;
X static short prep1 = 54, gap = 4, bias = 1023 ;
X static double novf = 1.7E308, nunf = 3.0E-308 ;
X#endif /* defined(vax) || defined(tahoe) */
X
Xdouble
Xscalb(x, N)
Xdouble x ;
Xint N ;
X{
X int k ;
X double scalb() ;
X
X#ifdef national
X unsigned short *px = (unsigned short *) &x + 3 ;
X#else /* national */
X unsigned short *px = (unsigned short *) &x ;
X#endif /* national */
X
X if (x == 0.0) return(x) ;
X
X#if defined(vax) || defined(tahoe)
X
X if ((k = *px & mexp) != ~msign)
X {
X if (N < -260) return(nunf * nunf) ;
X else if (N > 260)
X {
X extern double infnan(), copysign() ;
X return(copysign(infnan(ERANGE), x)) ;
X }
X#else /* defined(vax) || defined(tahoe) */
X
X if ((k = *px & mexp) != mexp)
X {
X if (N < -2100) return(nunf * nunf) ;
X else if (N > 2100) return(novf + novf) ;
X if (k == 0)
X {
X x *= scalb(1.0, (int) prep1) ;
X N -= prep1 ;
X return(scalb(x, N)) ;
X }
X#endif /* defined(vax) || defined(tahoe) */
X
X if ((k = (k >> gap) + N) > 0)
X if (k < (mexp >> gap)) *px = (*px & ~mexp) | (k << gap) ;
X else x = novf + novf ; /* overflow */
X else
X if (k > -prep1)
X { /* gradual underflow */
X *px = (*px & ~mexp) | (short) (1 << gap) ;
X x *= scalb(1.0, k-1) ;
X }
X else return(nunf * nunf) ;
X }
X return(x) ;
X}
X
X
Xdouble
Xcopysign(x, y)
Xdouble x, y ;
X{
X#ifdef national
X unsigned short *px = (unsigned short *) &x + 3,
X *py = (unsigned short *) &y + 3 ;
X#else /* national */
X unsigned short *px = (unsigned short *) &x,
X *py = (unsigned short *) &y ;
X#endif /* national */
X
X#if defined(vax) || defined(tahoe)
X if ((*px & mexp) == 0) return(x) ;
X#endif /* defined(vax) || defined(tahoe) */
X
X *px = (*px & msign) | (*py & ~msign) ;
X return(x) ;
X}
X
X
Xdouble
Xlogb(x)
Xdouble x ;
X{
X#ifdef national
X short *px = (short *) &x + 3, k ;
X#else /* national */
X short *px = (short *) &x, k ;
X#endif /* national */
X
X#if defined(vax) || defined(tahoe)
X return (int) (((*px & mexp) >> gap) - bias) ;
X#else /* defined(vax) || defined(tahoe) */
X
X if ((k = *px & mexp) != mexp)
X if (k != 0) return((k >> gap) - bias) ;
X else if (x != 0.0) return(-1022.0) ;
X else return(-(1.0 / zero)) ;
X else if (x != x) return(x) ;
X else
X {
X *px &= msign ;
X return(x) ;
X }
X#endif /* defined(vax) || defined(tahoe) */
X}
X
X
Xfinite(x)
Xdouble x ;
X{
X#if defined(vax) || defined(tahoe)
X return(1) ;
X#else /* defined(vax) || defined(tahoe) */
X#ifdef national
X return((*((short *) &x + 3) & mexp) != mexp) ;
X#else /* national */
X return((*((short *) &x) & mexp) != mexp) ;
X#endif /* national */
X#endif /* defined(vax) || defined(tahoe) */
X}
X
X
Xdouble
Xdrem(x, p)
Xdouble x, p ;
X{
X short sign ;
X double hp, dp, tmp, drem(), scalb() ;
X unsigned short k ;
X
X#ifdef national
X unsigned short *px = (unsigned short *) &x + 3,
X *pp = (unsigned short *) &p + 3,
X *pd = (unsigned short *) &dp + 3,
X *pt = (unsigned short *) &tmp + 3 ;
X#else /* national */
X unsigned short *px = (unsigned short *) &x,
X *pp = (unsigned short *) &p ,
X *pd = (unsigned short *) &dp ,
X *pt = (unsigned short *) &tmp ;
X#endif /* national */
X
X *pp &= msign ;
X
X#if defined(vax) || defined(tahoe)
X if ((*px & mexp) == ~msign) /* is x a reserved operand? */
X#else /* defined(vax) || defined(tahoe) */
X if ((*px & mexp) == mexp)
X#endif /* defined(vax) || defined(tahoe) */
X return (x-p) - (x-p) ; /* create nan if x is inf */
X
X if (p == 0.0)
X {
X doerr("drem", "SINGULARITY", EDOM) ;
X#if defined(vax) || defined(tahoe)
X extern double infnan() ;
X return(infnan(EDOM)) ;
X#else /* defined(vax) || defined(tahoe) */
X return(0.0 / zero) ;
X#endif /* defined(vax) || defined(tahoe) */
X }
X
X#if defined(vax) || defined(tahoe)
X if ((*pp & mexp) == ~msign) /* is p a reserved operand? */
X#else /* defined(vax) || defined(tahoe) */
X if ((*pp & mexp) == mexp)
X#endif /* defined(vax)||defined(tahoe) */
X {
X if (p != p) return p ;
X else return x ;
X }
X else if (((*pp & mexp) >> gap) <= 1)
X
X/* subnormal p, or almost subnormal p */
X
X {
X double b ;
X b = scalb(1.0, (int) prep1) ;
X p *= b ;
X x = drem(x, p) ;
X x *= b ;
X return(drem(x, p) / b) ;
X }
X else if (p >= novf / 2)
X {
X p /= 2 ;
X x /= 2 ;
X return(drem(x, p) * 2) ;
X }
X else
X {
X dp = p + p ;
X hp = p / 2 ;
X sign = *px & ~msign ;
X *px &= msign ;
X while (x > dp)
X {
X k = (*px & mexp) - (*pd & mexp) ;
X tmp = dp ;
X *pt += k ;
X
X#if defined(vax) || defined(tahoe)
X if (x < tmp) *pt -= 128 ;
X#else /* defined(vax) || defined(tahoe) */
X if (x < tmp) *pt -= 16 ;
X#endif /* defined(vax) || defined(tahoe) */
X
X x -= tmp ;
X }
X if (x > hp)
X {
X x -= p ;
X if (x >= hp) x -= p ;
X }
X
X#if defined(vax) || defined(tahoe)
X if (x)
X#endif /* defined(vax) || defined(tahoe) */
X *px ^= sign ;
X return(x) ;
X }
X}
X
X
X/* exp__E(x,c)
X * ASSUMPTION: c << x SO THAT fl(x+c)=x.
X * (c is the correction term for x)
X * exp__E RETURNS
X *
X * / exp(x+c) - 1 - x , 1E-19 < |x| < .3465736
X * exp__E(x,c) = |
X * \ 0 , |x| < 1E-19.
X *
X * DOUBLE PRECISION (IEEE 53 bits, VAX D FORMAT 56 BITS)
X * KERNEL FUNCTION OF EXP, EXPM1, POW FUNCTIONS
X * CODED IN C BY K.C. NG, 1/31/85;
X * REVISED BY K.C. NG on 3/16/85, 4/16/85.
X *
X * Required system supported function:
X * copysign(x,y)
X *
X * Method:
X * 1. Rational approximation. Let r=x+c.
X * Based on
X * 2 * sinh(r/2)
X * exp(r) - 1 = ---------------------- ,
X * cosh(r/2) - sinh(r/2)
X * exp__E(r) is computed using
X * x*x (x/2)*W - ( Q - ( 2*P + x*P ) )
X * --- + (c + x*[---------------------------------- + c ])
X * 2 1 - W
X * where P := p1*x^2 + p2*x^4,
X * Q := q1*x^2 + q2*x^4 (for 56 bits precision, add q3*x^6)
X * W := x/2-(Q-x*P),
X *
X * (See the listing below for the values of p1,p2,q1,q2,q3. The poly-
X * nomials P and Q may be regarded as the approximations to sinh
X * and cosh :
X * sinh(r/2) = r/2 + r * P , cosh(r/2) = 1 + Q . )
X *
X * The coefficients were obtained by a special Remez algorithm.
X *
X * Approximation error:
X *
X * | exp(x) - 1 | 2**(-57), (IEEE double)
X * | ------------ - (exp__E(x,0)+x)/x | <=
X * | x | 2**(-69). (VAX D)
X *
X * Constants:
X * The hexadecimal values are the intended ones for the following constants.
X * The decimal values may be used, provided that the compiler will convert
X * from decimal to binary accurately enough to produce the hexadecimal values
X * shown.
X */
X
X#if defined(vax) || defined(tahoe) /* VAX D format */
X#ifdef vax
X#define _0x(A,B) 0x/**/A/**/B
X#else /* vax */
X#define _0x(A,B) 0x/**/B/**/A
X#endif /* vax */
X
X/* static double */
X/* p1 = 1.5150724356786683059E-2 , Hex 2^ -6 * .F83ABE67E1066A */
X/* p2 = 6.3112487873718332688E-5 , Hex 2^-13 * .845B4248CD0173 */
X/* q1 = 1.1363478204690669916E-1 , Hex 2^ -3 * .E8B95A44A2EC45 */
X/* q2 = 1.2624568129896839182E-3 , Hex 2^ -9 * .A5790572E4F5E7 */
X/* q3 = 1.5021856115869022674E-6 ; Hex 2^-19 * .C99EB4604AC395 */
X
Xstatic long p1x[] = { _0x(3abe,3d78), _0x(066a,67e1) } ;
Xstatic long p2x[] = { _0x(5b42,3984), _0x(0173,48cd) } ;
Xstatic long q1x[] = { _0x(b95a,3ee8), _0x(ec45,44a2) } ;
Xstatic long q2x[] = { _0x(7905,3ba5), _0x(f5e7,72e4) } ;
Xstatic long q3x[] = { _0x(9eb4,36c9), _0x(c395,604a) } ;
X#define p1 (*(double*) p1x)
X#define p2 (*(double*) p2x)
X#define q1 (*(double*) q1x)
X#define q2 (*(double*) q2x)
X#define q3 (*(double*) q3x)
X#else /* defined(vax) || defined(tahoe) */
X
Xstatic double
Xp1 = 1.3887401997267371720E-2, /* Hex 2^ -7 * 1.C70FF8B3CC2CF */
Xp2 = 3.3044019718331897649E-5, /* Hex 2^-15 * 1.15317DF4526C4 */
Xq1 = 1.1110813732786649355E-1, /* Hex 2^ -4 * 1.C719538248597 */
Xq2 = 9.9176615021572857300E-4 ; /* Hex 2^-10 * 1.03FC4CB8C98E8 */
X#endif /* defined(vax) || defined(tahoe) */
X
Xdouble
Xexp__E(x, c)
Xdouble x, c ;
X{
X static double small = 1.0E-19 ;
X double copysign(), z, p, q, xp, xh, w ;
X
X if (copysign(x, 1.0) > small)
X {
X z = x * x ;
X p = z * (p1 + z * p2) ;
X
X#if defined(vax) || defined(tahoe)
X q = z * (q1 + z * (q2 + z * q3)) ;
X#else /* defined(vax) || defined(tahoe) */
X q = z * (q1 + z * q2) ;
X#endif /* defined(vax) || defined(tahoe) */
X xp = x * p ;
X xh = x * 0.5 ;
X w = xh - (q - xp) ;
X p = p + p ;
X c += x * ((xh * w - (q - (p + xp))) / (1.0 - w) + c) ;
X return(z * 0.5 + c) ;
X }
X
X/* end of |x| > small */
X
X else
X {
X if (x != 0.0) 1.0 + small ; /* raise the inexact flag */
X return(copysign(0.0, x)) ;
X }
X}
X
X
X/* log__L(Z)
X * LOG(1+X) - 2S X
X * RETURN --------------- WHERE Z = S*S, S = ------- , 0 <= Z <= .0294... * S 2 + X
X *
X * DOUBLE PRECISION (VAX D FORMAT 56 bits or IEEE DOUBLE 53 BITS)
X * KERNEL FUNCTION FOR LOG; TO BE USED IN LOG1P, LOG, AND POW FUNCTIONS
X * CODED IN C BY K.C. NG, 1/19/85;
X * REVISED BY K.C. Ng, 2/3/85, 4/16/85.
X *
X * Method :
X * 1. Polynomial approximation: let s = x/(2+x).
X * Based on log(1+x) = log(1+s) - log(1-s)
X * = 2s + 2/3 s**3 + 2/5 s**5 + .....,
X *
X * (log(1+x) - 2s)/s is computed by
X *
X * z*(L1 + z*(L2 + z*(... (L7 + z*L8)...)))
X *
X * where z=s*s. (See the listing below for Lk's values.) The
X * coefficients are obtained by a special Remez algorithm.
X *
X * Accuracy:
X * Assuming no rounding error, the maximum magnitude of the approximation
X * error (absolute) is 2**(-58.49) for IEEE double, and 2**(-63.63)
X * for VAX D format.
X *
X * Constants:
X * The hexadecimal values are the intended ones for the following constants.
X * The decimal values may be used, provided that the compiler will convert
X * from decimal to binary accurately enough to produce the hexadecimal values
X * shown.
X */
X
X#if defined(vax) || defined(tahoe) /* VAX D format (56 bits) */
X#ifdef vax
X#define _0x(A,B) 0x/**/A/**/B
X#else /* vax */
X#define _0x(A,B) 0x/**/B/**/A
X#endif /* vax */
X
X/* static double */
X/* L1 = 6.6666666666666703212E-1 , Hex 2^ 0 * .AAAAAAAAAAAAC5 */
X/* L2 = 3.9999999999970461961E-1 , Hex 2^ -1 * .CCCCCCCCCC2684 */
X/* L3 = 2.8571428579395698188E-1 , Hex 2^ -1 * .92492492F85782 */
X/* L4 = 2.2222221233634724402E-1 , Hex 2^ -2 * .E38E3839B7AF2C */
X/* L5 = 1.8181879517064680057E-1 , Hex 2^ -2 * .BA2EB4CC39655E */
X/* L6 = 1.5382888777946145467E-1 , Hex 2^ -2 * .9D8551E8C5781D */
X/* L7 = 1.3338356561139403517E-1 , Hex 2^ -2 * .8895B3907FCD92 */
X/* L8 = 1.2500000000000000000E-1 , Hex 2^ -2 * .80000000000000 */
X
Xstatic long L1x[] = { _0x(aaaa,402a), _0x(aac5,aaaa) } ;
Xstatic long L2x[] = { _0x(cccc,3fcc), _0x(2684,cccc) } ;
Xstatic long L3x[] = { _0x(4924,3f92), _0x(5782,92f8) } ;
Xstatic long L4x[] = { _0x(8e38,3f63), _0x(af2c,39b7) } ;
Xstatic long L5x[] = { _0x(2eb4,3f3a), _0x(655e,cc39) } ;
Xstatic long L6x[] = { _0x(8551,3f1d), _0x(781d,e8c5) } ;
Xstatic long L7x[] = { _0x(95b3,3f08), _0x(cd92,907f) } ;
Xstatic long L8x[] = { _0x(0000,3f00), _0x(0000,0000) } ;
X
X#define L1 (*(double*) L1x)
X#define L2 (*(double*) L2x)
X#define L3 (*(double*) L3x)
X#define L4 (*(double*) L4x)
X#define L5 (*(double*) L5x)
X#define L6 (*(double*) L6x)
X#define L7 (*(double*) L7x)
X#define L8 (*(double*) L8x)
X#else /* defined(vax) || defined(tahoe) */
X
Xstatic double
XL1 = 6.6666666666667340202E-1, /* Hex 2^ -1 * 1.5555555555592 */
XL2 = 3.9999999999416702146E-1, /* Hex 2^ -2 * 1.999999997FF24 */
XL3 = 2.8571428742008753154E-1, /* Hex 2^ -2 * 1.24924941E07B4 */
XL4 = 2.2222198607186277597E-1, /* Hex 2^ -3 * 1.C71C52150BEA6 */
XL5 = 1.8183562745289935658E-1, /* Hex 2^ -3 * 1.74663CC94342F */
XL6 = 1.5314087275331442206E-1, /* Hex 2^ -3 * 1.39A1EC014045B */
XL7 = 1.4795612545334174692E-1 ; /* Hex 2^ -3 * 1.2F039F0085122 */
X#endif /* defined(vax) || defined(tahoe) */
X
Xdouble
Xlog__L(z)
Xdouble z ;
X{
X#if defined(vax) || defined(tahoe)
X return(z * (L1 + z * (L2 + z * (L3 + z * (L4 + z *
X (L5 + z * (L6 + z * (L7 + z * L8)))))))) ;
X#else /* defined(vax) || defined(tahoe) */
X return(z * (L1 + z * (L2 + z * (L3 + z * (L4 + z *
X (L5 + z * (L6 + z * L7))))))) ;
X#endif /* defined(vax) || defined(tahoe) */
X}
X#endif /*NEED_POW*/
SHAR_EOF
echo "File mathlib.c is complete"
chmod 0444 mathlib.c || echo "restore of mathlib.c fails"
set `wc -c mathlib.c`;Sum=$1
if test "$Sum" != "66590"
then echo original size 66590, current size $Sum;fi
rm -f s2_seq_.tmp
echo "You have unpacked the last part"
exit 0