Release calc version 2.12.4.12

BUGS

@@ -84,7 +84,7 @@ Problems that have known work-a-rounds:
 
 mis-features in calc:
 
-    Some problems are not bugs but rarther mis-features / things that could
+    Some problems are not bugs but rather mis-features / things that could
     work better.  The following is a list of mis-features that should be
     addressed and improved someday.
 
@@ -129,7 +129,7 @@ mis-features in calc:
 
       will not.
 
-## Copyright (C) 1999-2007  Landon Curt Noll
+## Copyright (C) 1999-2013  Landon Curt Noll
 ##
 ## Calc is open software; you can redistribute it and/or modify it under
 ## the terms of the version 2.1 of the GNU Lesser General Public License
@@ -145,9 +145,9 @@ mis-features in calc:
 ## received a copy with calc; if not, write to Free Software Foundation, Inc.
 ## 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
 ##
-## @(#) $Revision: 30.1 $
-## @(#) $Id: BUGS,v 30.1 2007/03/16 11:09:46 chongo Exp $
-## @(#) $Source: /usr/local/src/cmd/calc/RCS/BUGS,v $
+## @(#) $Revision: 30.2 $
+## @(#) $Id: BUGS,v 30.2 2013/08/11 01:09:27 chongo Exp $
+## @(#) $Source: /usr/local/src/bin/calc/RCS/BUGS,v $
 ##
 ## Under source code control:	1994/03/18 14:06:13
 ## File existed as early as:	1994

CHANGES

@@ -1,4 +1,368 @@
-The following are the changes from calc version 2.12.1.10 to date:
+The following are the changes from calc version 2.12.4.11 to date:
+
+    Fixed many typos in comments of the Makefile thanks to the review
+    work of Michael Somos.
+
+    Fixed typo in "help sysinfo".
+
+    The Makefile rule, debug, is now more verbose and prints more information
+    about the calc compiled constants.
+
+    Added a more of calc resource files by
+    Christoph Zurnieden <czurnieden at gmx dot de> including:
+
+    infinities.cal - handle infinities symbolically, a little helper file
+    intnum.cal - implementation of tanh sinh and Gauss-Legendre quadrature
+    smallfactors.cal - find the factors of a number < 2^32
+    strings.cal - implementation of isascii() and isblank()
+
+    Reformatted some calc resource files.  Cleanup in comment the headers
+    of some calc resource files.
+
+    Minor formatting changes to a few help files.
+
+    No need to be special picky about the test8900.cal calc resource file.
+
+    Added a number of ctype-like builtins:
+
+    isalnum - whether character is alpha-numeric
+    isalpha - whether character is alphabetic
+    iscntrl - whether character is a control character
+    isdigit - whether character is a digit character
+    isgraph - whether character is a graphical character
+    islower - whether character is lower case
+    isprint - whether character is a printable
+    ispunct - whether character is a punctuation
+    isspace - whether character is a space character
+    isupper - whether character is upper case
+    isxdigit - whether character a hexadecimal digit
+    strcasecmp - compare two strings, case independent
+    strncasecmp - compare two strings up to n characters, case independent
+    strtolower - transform an ASCII string to lower case
+    strtoupper - transform an ASCII string to upper case
+
+    For details on these new builtins, see their help messages.
+    Thanks goes to Inge Zurnieden <inge dot zurnieden at gmx dot de> for
+    these new builtins.
+
+    Calc source code is now picky v2.3 clean using:
+
+	picky -s -v file file2 ..
+
+    With the exception of:
+
+	help/errorcodes.sed
+	cal/set8700.line
+
+    Due to the long lines in those files, we use:
+
+	picky -w -s -v help/errorcodes.sed cal/set8700.line
+
+    For more information about the picky tool, see:
+
+	http://cis.csuohio.edu/~somos/picky.html
+
+    Removed functions from strings.cal that have been replaced by
+    the new ctype-like builtin functions.
+
+    Fixed cal/Makefile to include missing intnum.cal file.
+
+    Added detail_help_list make target to cal/Makefile.
+
+    The detaillist make target in help/Makefile is now
+    called detail_help_list.
+
+
+The following are the changes from calc version 2.12.4.6 to version 2.12.4.10:
+
+    Updated RPM build process to remove use of deprecated flags.
+
+    Applied a number of fixes to calc.spec and rpm.mk file.
+    See calc.spec.in for details.  Changed rpm release to 2.1.
+
+    Set MACOSX_DEPLOYMENT_TARGET=10.8 as we upgraded all of
+    our development Mac OS X to 10.8.
+
+    Libraries are chmodded as 0644 to allow for building rpms
+    without root.
+
+    Silenced annoying warning about unused variable 'intp'
+    while compiline endian.c under some circumstances.
+
+    Fixed typo in redeclaration warnings.  Thanks to
+    Christoph Zurnieden <czurnieden at gmx dot de> for this report.
+
+    Added a number of calc resource files by
+    Christoph Zurnieden <czurnieden at gmx dot de> including:
+
+    bernpoly.cal - Computes the nth Bernoulli polynomial at z for any n,z
+    brentsolve.cal - root-finder implementwed with the Brent-Dekker trick
+    factorial.cal - product of the positive integers
+    factorial2.cal - variety of integer functions quasi-related to factoral
+    lambertw.cal - Computes Lambert's W-function at "z" at branch "branch"
+    lnseries.cal - Calculates a series of natural logarithms at 1,2,3,4...n
+    specialfunctions.cal - Calculates the value of the beta function
+    statistics.cal - a wide vareity of stastical functions
+    toomcook.cal - Multiply by way of the Toom-Cook algorithm
+    zeta2.cal - Calculate the value of the Hurwitz Zeta function
+
+    Fixed a makefile bug that prevented the those new calc resource
+    files from being installed.
+
+    Improved the formatting of the output from:
+
+	help resource
+
+    We replaced COPYING-LGPL with the version that is found at
+    http://www.gnu.org/licenses/lgpl-2.1.txt because that version
+    contans some whitespace formatting cleanup.  Otherwise the
+    license is the same.
+
+    We fixed a number of places where "the the" was used
+    when just "the" should be used.
+
+	NOTE: Fixes to grammar, spelling and minor formatting
+	      problems are welcome.  Please send us your patches!
+
+    With the exception of 3 source files, we became "picky" about
+    line lengths and other issues reported by the picky tool:
+
+    	cal/test8900.cal
+	cal/set8700.line
+	help/errorcodes.sed
+
+    The above 3 files now pass picky -w (OK except for line length).
+    For more information about the picky tool, see:
+
+	http://cis.csuohio.edu/~somos/picky.html
+
+
+The following are the changes from calc version 2.12.4.3 to 2.12.4.5:
+
+    Added gvec.cal resource script.
+
+    Added calc-symlink make rule to setup symlinks from stardard locations
+    into a tree specified by a non-empty ${T} makefile variable.  Added
+    calc-unsymlink to remove any symlinks that may have been created by
+    the calc-symlink rule.
+
+    If is OK for the calc-symlink make rule to pre-remove a symlink.
+
+    Fixed bug was uncovered in calc that caused script failures when calc
+    is called within a while loop in BASH if the while loop is fed from
+    stdin due to calc's redirection/inheritance of stdin and no option
+    to change this behavior.  Thanks gores to David C. Rankin
+    <drankinatty at gmail dot com> for the bug fix and to David Haller
+    <dnh at opensuse dot org> for helping debug this problem.
+
+
+The following are the changes from calc version 2.12.4.0 to 2.12.4.2:
+
+    Fixed a documentation bug for the sgn() builtin.
+
+    Added the 1<<8/2 evaluation example to "help unexpected".  That
+    expression evalutes to 128, not 16 as some C programmers might expect.
+
+    Fixed a bug in solve.cal where high was not returned in some situations.
+
+    Fixed a bug reported by Paul & Karen Tomlinson (paulnkaz at pktomlinson
+    dot fsnet dot co dot uk) where calling log multiple times with different
+    values of epsilon resulted in an incorrect value.
+
+    Removed cvd rule from Makefiles.
+
+    The Makefile used in the source rpm (calc-*.src.rpm) no longer uses
+    the -Werror compile flag.  This is to help those distributions with
+    compilers that make produce (hopefully) complination warnings.
+    NOTE: For testing and calc build purposes will recommend and will
+    continue to use the -Werror flag.
+
+    Fixed a typo in the Makefile where the make variable ${SAMPLE_OBJ}
+    was misspelled as ${SAMPLE_OBJS}.
+
+    Added prep makefile rule to make is easier to compile calc without
+    an optimizer.  By doing:
+
+    	make clobber prep
+
+    one may build a calc binary that is easier to debug.
+
+    Fixed a bug where an certains typos (e.g., calling an unknown
+    function) would previously cause calc to exit.
+
+    Updated the COPYING file to reflect the new filenames associated
+    with the SHA1 hash function, and removed mention of files related
+    to the SHA (SHA0, not SHA1) and the MD5 hash functions (which is
+    no longer supported in calc).
+
+    Fixed a bug where a calling vsnprintf() twice created problems.
+    The thanks for this fix goes to Matthew Miller (mattdm at mattdm
+    dot org) for this patch.
+
+    Michael Penk (mpenk at wuska dot com) reported success in installs
+    under windoz via Cygwin by making a change to the Cygwin target.
+    These changes have been folded into the main calc Makefile.
+    The old recommendation of using 'make win32_hsrc' is no longer
+    required for Cygwin.  See the README.WINDOWS file for details.
+
+    Added dms.cal and hms.cal resource files.  The dms.cal is a more
+    functional version of deg.cal.  It is a superset except that increment
+    and decrement is on the arc second level.  The hms.cal is for
+    24-hour cycle instread of the 360 degree cycle of dms.cal.
+
+    Changed deg.cal object name from dms to deg so that the more functional
+    dms.cal can own the dms object name.
+
+    Updated 'help obj' to reflect changes to 'show objfunctions' and
+    resource file example list since 1999.
+
+    Fixed problem where CALC_BYTE_ORDER refering to CALC_BIG_ENDIAN
+    and CALC_LITTLE_ENDIAN instead of BIG_ENDIAN and LITTLE_ENDIAN.
+
+
+The following are the changes from calc version 2.12.3.0 to 2.12.3.3:
+
+    Fixed the Jacobi function where it returned 1 when it should have
+    returned 0.  Thanks goes to Kevin Sopp (baraclese at googlemail dot com)
+    for discovering the problem and suggesting the nature if the fix.
+
+    Calc versions will always be of the form x.y.z.w even when the
+    MINOR_PATCH (w) is 0.  Thus, 2.12.3.0 will be printed as 2.12.3.0
+    instread of just 2.12.3.
+
+    Added MINGW32_NT-5.0 compile target based on a patch from
+    Brian L. Angus (angus at eng dot utah dot edu).
+
+    Removed the use of rpm.release in the Makefile.
+
+    Mac OS Darwin targets no longer attempt to use ldconfig.  Under the
+    Darwin target, the LDCONFIG make variable is redefined to be
+    an empty value.  Thanks goes to Ralf Trinler (art at infra dot de)
+    for reporting this problem.
+
+    The ${CALC_INCDIR}/custom is no longer being removed at install time
+    if it is empty.  Now when ${ALLOW_CUSTOM} make variable is empty,
+    an empty ${CALC_INCDIR}/custom may be left hehind.
+
+    Fixed a problem where a "make clobber" would remove custom/Makefile
+    and fail to rebuilt it.
+
+
+The following are the changes from calc version 2.12.2.3 to 2.12.2.4:
+
+    Added OpenBSD target.
+
+    Using the -r test instead of the -e test in Makefiles because some
+    out of date shells still do not have the -e test.
+
+    The Makefile now avoids the use of if ! command because some out of
+    date shells to not support the ! construct.
+
+
+The following are the changes from calc version 2.12.1.1 to 2.12.2.2:
+
+    Added an explicit Solaris target.
+
+    Fixed confusion in Makefile where some uses of ${EXT} were misnamed ${EXE}.
+
+    Added a "make strip" rule, per suggestion from Igor Furlan <primorec
+    at sbcglobal dot net>, to allow one to strip previously built binary
+    executables and libraries.
+
+    Under the Darwin / OS X target, ${DARWIN_ARCH} is left empty meaning
+    that calc is compiled for the native CPU type instead of Universal
+    Binary (Intel and PPC).
+
+    By default, the calc binary that is built for the rpm forces
+    ${LD_SHARE} to be empty.  An empty ${LD_SHARE} means that the calc
+    from the rpm does not set rpath.  This in turn causes the default
+    system path to be searched when looking for libcalc and libcustcalc.
+
+    The Makefile shipped with calc still sets ${LD_SHARE} for host targets.
+    By default, the dynamic shared library search path for all targets
+    starts with the source directory.  Starting the search in the source
+    directory is convenient for testing and debugging but is not appropriate
+    for installation on a production system.  To get the same effect
+    as the calc binary in the calc rpm, try:
+
+	make clobber
+	make calc-dynamic-only BLD_TYPE=calc-dynamic-only LD_SHARE=
+	make install
+
+    The libcalc and libcustcalc shared libraries are now tied to
+    the 4 level calc version instead of just 3 levels.  For example,
+    under Linux calc version 2.12.2.1 uses /usr/lib/libcalc.so.2.12.2.1
+    instead of just the /usr/lib/libcalc.so.2.12.2 file.  This change
+    was made so that calc produced by 'make clobber; make all install'
+    is consistent with the calc rpm.
+
+    Calc is now releasing the calc-debuginfo rpm for those RPM users who
+    which to use non-stripped libraries and binaries for debugging
+    purposes.  By default, the calc rpm installed stripped binaries
+    and libraries.
+
+    Added this high priority item to the calc help/todo list:
+
+	It is overkill to have nearly everything wind up in libcalc.
+	Form a libcalcmath and a libcalclang so that an application
+	that just wants to link with the calc math libs can use them
+	without dragging in all of the other calc language, I/O,
+	and builtin functions.
+
+    Fixed the wording for the -i flag in the calc man page.
+
+    Added some notes to the help/unexpected file regarding calc
+    and interactice shells.
+
+    Fixed bug where a FILEPOS was copied FPOS_POS_BITS octets instead of
+    FPOS_POS_LEN octets.
+
+    Split out ${READLINE_EXTRAS} Makefile variables from ${READLINE_LIB}
+    to better deal with Fedora rpm requirements.
+
+    Bit 8 (0x80) of calc_debug is reserved for custom debugging.
+    See help/config and custom/HOW_TO_ADD for details.
+
+    When the Makefile variable ${ALLOW_CUSTOM} is not defined or empty,
+    the libcustcalc library is not built or linked against, certain make
+    rules skip going into the custom sub-directory, the install
+    rule skips certain custom installation actions, and the common
+    C flags (${COMMON_CFLAGS}) is given -UCUSTOM.  Other make rules such
+    as "make clean" and "make clobber" still work as before.  Also
+    the Makefile.simple assumes that the Makefile variable ${ALLOW_CUSTOM}
+    is -DCUSTOM.
+
+    Clarified that the calc builtin functions rand() and random()
+    operate over a half closed interval.  The help/rand and help/random
+    refer to  the top of the interval as "beyond" instead of "max".
+
+    Releaseing source tar balls using bzip2 instead of with gzip.  So
+    what was calc-something.tar.gz is now calc-something.tar.bz2.
+    To "uncompress" use:
+
+    	bunzip2 calc-something.tar.bz2
+
+    On some systems, one may untar directly by:
+
+	tar -jxvf calc-something.tar.bz2
+
+   The Makefile variable ${BYTE_ORDER} was replaced by ${CALC_BYTE_ORDER}.
+
+   Changed the way the Makefile can force the calc byte order.  If you set
+   the Makefile variable ${CALC_BYTE_ORDER} to be -DCALC_BIG_ENDIAN then
+   endian.h will force the CPP symbol CALC_BYTE_ORDER to be BIG_ENDIAN.
+   If you set ${CALC_BYTE_ORDER} to be -DCALC_LITTLE_ENDIAN then endian.h
+   will force the CPP symbol CALC_BYTE_ORDER to be LITTLE_ENDIAN.
+   If the Makefile variable ${CALC_BYTE_ORDER} is empty, then the CPP
+   symbol CALC_BYTE_ORDER will set to the CPP symbol BYTE_ORDER as
+   defined by some system include file (if the Makefile can find such
+   an include file), or the Makefile compiling endian.c and hopefully
+   using that result to set CPP symbol CALC_BYTE_ORDER.  Regardless of
+   how it happens, the CPP symbol CALC_BYTE_ORDER should end up set in
+   endian_calc.h include file.
+
+
+The following are the changes from calc version 2.12.1.10 to 2.12.2:
 
     Put back the missing -s flags on the cscripts:  mersenne, 4dsphere,
     fprodcut, plus, and powerterm.  Thanks goes to Bradley Reed
@@ -927,8 +1291,7 @@ The following are the changes from calc version 2.11.10 to 2.11.10:
 The following are the changes from calc version 2.11.9 to 2.11.9.3:
 
     Fixed calc man page examples to move -f to the end of the line.
-    Thanks goes to Michael Somos <somos at grail dot cba dot csuohio
-    dot edu> for pointing this out.
+    Thanks goes to Michael Somos for pointing this out.
 
     Linux and gcc now compiled with -Wall -W -Wno-comment.
 
@@ -1883,7 +2246,7 @@ The following are the changes from calc version 2.11.1 to 2.11.1t2.2:
     Cleaned up help display system.  Help file lines that begin with
     '##' are not displayed.
 
-    Calc source and documentation now uses the the these terms:
+    Calc source and documentation now uses these terms:
 
 	*.cal files	calc resource file
 	*.a files	calc binary link library
@@ -1921,7 +2284,7 @@ The following are the changes from calc version 2.11.0t10 to 2.11.0t11:
     Fixed whitespace to make the best use of 8 character tabs.
 
     Fixed some bugs relating to '// and %' in combination with some
-    of the the rounding modes based on a patch from Ernest Bowen
+    of the rounding modes based on a patch from Ernest Bowen
     <ernie at turing dot une dot edu dot au>.
 
     A patch from Klaus Alexander Seistrup <klaus at seistrup dot dk>, when
@@ -2115,7 +2478,7 @@ The following are the changes from calc version 2.11.0t10 to 2.11.0t11:
 
 	Where ``/usr/local/bin/calc'' is the path to the calc binary.
 	Additional -options may be added to the line, but it MUST
-	start with -S.	For example, the the executable file ``plus''
+	start with -S.	For example, the executable file ``plus''
 	contain the following:
 
 	    #!/usr/local/bin/calc -S -e
@@ -4578,7 +4941,7 @@ The following are the changes from calc version 2.10.3t0 to 2.10.3t2:
 
 	global variable: to the end of the session with calc
 
-	static within a function definition: to the the first of:
+	static within a function definition: to the first of:
 
 	    an end of a global, static or local declaration (including
 	    initialization code) with the same identifier
@@ -4676,7 +5039,7 @@ The following are the changes from calc version 2.10.3t0 to 2.10.3t2:
 
 	global a, b, mat A, B[2] = {3,4}, C[2] = {4,5}, obj point P = {5,6}, Q
 
-    Fixed some bugs related to global and static scoping.  See the the
+    Fixed some bugs related to global and static scoping.  See the
     5200 regress test and lib/test5200.cal for details.
 
     Optimized opcode generator so that functions defined using '=' do not
@@ -5173,7 +5536,7 @@ The following are the changes from calc version 2.10.2t4 to 2.10.2t24:
 	 ILDFLAGS are flags given to ${CC} for linking .o files
 		  for intermediate progs
 
-	 CC is how the the C compiler is invoked
+	 CC is how the C compiler is invoked
 
     Added more tests to regress.cal.
 
@@ -5398,7 +5761,7 @@ The following are the changes from calc version 2.10.2t1 to 2.10.2t3:
     Some systems has a libc symbolic qadd() that conflicted with calc's
     qadd function.  To avoid this, qadd() has been renamed to qqadd().
 
-    The calc error codes are produced from the the calcerr.tbl file.
+    The calc error codes are produced from the calcerr.tbl file.
     Instead of changing #defines in value.h, one can not edit calcerr.tbl.
     The Makefile builds calcerr.h from this file.
 
@@ -5941,7 +6304,7 @@ The following are the changes from calc version 2.9.3t11 to 2.10.0t12:
 	LDFLAGS are flags given to ${CC} for linking .o files
 	ILDFLAGS are given to ${CC} for linking .o's for intermediate progs
 
-	CC is how the the C compiler is invoked
+	CC is how the C compiler is invoked
 
     The syntax error:
 
@@ -6592,9 +6955,9 @@ Following is a list of visible changes to calc from version 1.24.7 to 1.26.1:
 ## received a copy with calc; if not, write to Free Software Foundation, Inc.
 ## 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
 ##
-## @(#) $Revision: 30.8 $
-## @(#) $Id: CHANGES,v 30.8 2007/09/02 06:26:38 chongo Exp $
-## @(#) $Source: /usr/local/src/cmd/calc/RCS/CHANGES,v $
+## @(#) $Revision: 30.37 $
+## @(#) $Id: CHANGES,v 30.37 2013/09/02 02:31:57 chongo Exp $
+## @(#) $Source: /usr/local/src/bin/calc/RCS/CHANGES,v $
 ##
 ## Under source code control:	1993/06/02 18:12:57
 ## File existed as early as:	1989

COPYING

@@ -6,17 +6,17 @@ This file is Copyrighted
 
     This file is covered under the following Copyright:
 
-	Copyright (C) 1999  Landon Curt Noll
+	Copyright (C) 1999-2008  Landon Curt Noll
 	All rights reserved.
 
 	Everyone is permitted to copy and distribute verbatim copies
 	of this license document, but changing it is not allowed.
 
-    # @(#) $Revision: 30.1 $
-    # @(#) $Id: COPYING,v 30.1 2007/03/16 11:09:46 chongo Exp $
-    # @(#) $Source: /usr/local/src/cmd/calc/RCS/COPYING,v $
+    # @(#) $Revision: 30.4 $
+    # @(#) $Id: COPYING,v 30.4 2013/09/01 20:14:30 chongo Exp $
+    # @(#) $Source: /usr/local/src/bin/calc/RCS/COPYING,v $
 
-=-=
+-=-
 
 Calc is covered by the GNU Lesser General Public License
 --------------------------------------------------------
@@ -78,7 +78,7 @@ Calc is covered by the GNU Lesser General Public License
 
     Feel free to follow the name line with additional EMail text as desired.
 
-=-=
+-=-
 
     Calc bug reports and calc bug fixes should be sent to:
 
@@ -93,7 +93,7 @@ Calc is covered by the GNU Lesser General Public License
 
     You may have additional words in your subject line.
 
-=-=
+-=-
 
 Calc's relationship to the GNU Lesser General Public License
 ------------------------------------------------------------
@@ -149,7 +149,7 @@ Calc's relationship to the GNU Lesser General Public License
      except for the exception files explicitly listed in the ``Calc
      copyrights and exception files'' section below.
 
-=-=
+-=-
 
 Calc copyrights and exception files
 -----------------------------------
@@ -165,12 +165,12 @@ Calc copyrights and exception files
 	Copyright (C) year  Ernest Bowen and Landon Curt Noll
 	Copyright (C) year  Ernest Bowen
 	Copyright (C) year  Petteri Kettunen and Landon Curt Noll
+	Copyright (C) year  Christoph Zurnieden
 
     These files are not covered under one of the Copyrights listed above:
 
-	    shs1.c	    shs1.h	 shs.c		shs.h
-	    md5.c	    md5.h	 COPYING	COPYING-LGPL
-	    cal/qtime.cal   cal/screen.cal
+	    sha1.c		sha1.h		COPYING
+	    COPYING-LGPL	cal/qtime.cal	cal/screen.cal
 
     The file COPYING-LGPL, which contains a copy of the version 2.1
     GNU Lesser General Public License, is itself Copyrighted by the
@@ -182,17 +182,14 @@ Calc copyrights and exception files
     top of this file.  It is important to note that you may distribute
     verbatim copies of this file but you may not modify this file.
 
-    Some of these exception files are in the public domain.  The md5.c
-    and md5.h files were "derived from the RSA Data Security, Inc. MD5
-    Message-Digest Algorithm" and are under a copyright that allows these
-    two files to be freely used and distributed.  Other files are under
-    the LGPL but have different authors that those listed above.
+    Some of these exception files are in the public domain.  Other files
+    are under the LGPL but have different authors that those listed above.
 
     In all cases one may use and distribute these exception files freely.
     And because one may freely distribute the LGPL covered files, the
     entire calc source may be freely used and distributed.
 
-=-=
+-=-
 
 General Copyleft and License info
 ---------------------------------
@@ -206,7 +203,7 @@ General Copyleft and License info
 	http://www.gnu.org/copyleft/lesser.html
 	http://www.gnu.org/copyleft/lesser.txt
 
-=-=
+-=-
 
 Why calc did not use the GNU General Public License
 ---------------------------------------------------

COPYING-LGPL

@@ -1,8 +1,8 @@
-		  GNU LESSER GENERAL PUBLIC LICENSE
-		       Version 2.1, February 1999
+                  GNU LESSER GENERAL PUBLIC LICENSE
+                       Version 2.1, February 1999
 
  Copyright (C) 1991, 1999 Free Software Foundation, Inc.
-     51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA
+ 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
  Everyone is permitted to copy and distribute verbatim copies
  of this license document, but changing it is not allowed.
 
@@ -10,7 +10,7 @@
  as the successor of the GNU Library Public License, version 2, hence
  the version number 2.1.]
 
-			    Preamble
+                            Preamble
 
   The licenses for most software are designed to take away your
 freedom to share and change it.  By contrast, the GNU General Public
@@ -112,7 +112,7 @@ modification follow.  Pay close attention to the difference between a
 former contains code derived from the library, whereas the latter must
 be combined with the library in order to run.
 
-		  GNU LESSER GENERAL PUBLIC LICENSE
+                  GNU LESSER GENERAL PUBLIC LICENSE
    TERMS AND CONDITIONS FOR COPYING, DISTRIBUTION AND MODIFICATION
 
   0. This License Agreement applies to any software library or other
@@ -146,7 +146,7 @@ such a program is covered only if its contents constitute a work based
 on the Library (independent of the use of the Library in a tool for
 writing it).  Whether that is true depends on what the Library does
 and what the program that uses the Library does.
-  
+
   1. You may copy and distribute verbatim copies of the Library's
 complete source code as you receive it, in any medium, provided that
 you conspicuously and appropriately publish on each copy an
@@ -432,7 +432,7 @@ decision will be guided by the two goals of preserving the free status
 of all derivatives of our free software and of promoting the sharing
 and reuse of software generally.
 
-			    NO WARRANTY
+                            NO WARRANTY
 
   15. BECAUSE THE LIBRARY IS LICENSED FREE OF CHARGE, THERE IS NO
 WARRANTY FOR THE LIBRARY, TO THE EXTENT PERMITTED BY APPLICABLE LAW.
@@ -455,7 +455,7 @@ FAILURE OF THE LIBRARY TO OPERATE WITH ANY OTHER SOFTWARE), EVEN IF
 SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE POSSIBILITY OF SUCH
 DAMAGES.
 
-		     END OF TERMS AND CONDITIONS
+                     END OF TERMS AND CONDITIONS
 
            How to Apply These Terms to Your New Libraries
 
@@ -485,7 +485,8 @@ convey the exclusion of warranty; and each file should have at least the
 
     You should have received a copy of the GNU Lesser General Public
     License along with this library; if not, write to the Free Software
-    Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA
+    Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301
+    USA
 
 Also add information on how to contact you by electronic and paper mail.
 
@@ -500,5 +501,3 @@ necessary.  Here is a sample; alter the names:
   Ty Coon, President of Vice
 
 That's all there is to it!
-
-

HOWTO.INSTALL

@@ -1,4 +1,4 @@
-Installing calc from the gziped tarball in 4 easy steps:
+Installing calc from the bzip2-ed tarball in 4 easy steps:
 
     0)  If your platform supports i686 RPMs, you may want to go to:
 
@@ -17,13 +17,13 @@ Installing calc from the gziped tarball in 4 easy steps:
 
 	The following 4 steps apply to calc source tree that comes from either:
 
-		gunzip -c calc-*.tar.gz | tar -xvf -
+		bunzip2 -c calc-*.tar.bz2 | tar -xvf -
 
 	or from:
 
 		rpm -ivh calc-*.src.rpm
 		cd /var/tmp
-		gunzip -c /usr/src/redhat/SOURCES/calc-*.tar.gz | tar -xvf -
+		bunzip2 -c /usr/src/redhat/SOURCES/calc-*.tar.bz2 | tar -xvf -
 
     1)	Look at the makefile, and adjust it to suit your needs.
 
@@ -204,9 +204,9 @@ the calc help subsystem.  See the README file for details.
 ## received a copy with calc; if not, write to Free Software Foundation, Inc.
 ## 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
 ##
-## @(#) $Revision: 30.5 $
-## @(#) $Id: HOWTO.INSTALL,v 30.5 2007/09/01 19:54:03 chongo Exp $
-## @(#) $Source: /usr/local/src/cmd/calc/RCS/HOWTO.INSTALL,v $
+## @(#) $Revision: 30.6 $
+## @(#) $Id: HOWTO.INSTALL,v 30.6 2007/10/16 12:22:22 chongo Exp $
+## @(#) $Source: /usr/local/src/bin/calc/RCS/HOWTO.INSTALL,v $
 ##
 ## Under source code control:	1999/09/27 20:48:44
 ## File existed as early as:	1999

LIBRARY

@@ -640,7 +640,7 @@ need call libcalc_call_me_last() only once.
 ##
 ## @(#) $Revision: 30.1 $
 ## @(#) $Id: LIBRARY,v 30.1 2007/03/16 11:09:46 chongo Exp $
-## @(#) $Source: /usr/local/src/cmd/calc/RCS/LIBRARY,v $
+## @(#) $Source: /usr/local/src/bin/calc/RCS/LIBRARY,v $
 ##
 ## Under source code control:	1993/07/30 19:44:49
 ## File existed as early as:	1993

README

@@ -136,7 +136,7 @@ The calc web site is located at:
 ##
 ## @(#) $Revision: 30.1 $
 ## @(#) $Id: README,v 30.1 2007/03/16 11:09:46 chongo Exp $
-## @(#) $Source: /usr/local/src/cmd/calc/RCS/README,v $
+## @(#) $Source: /usr/local/src/bin/calc/RCS/README,v $
 ##
 ## Under source code control:	1995/10/25 05:27:59
 ## File existed as early as:	1995

README.WINDOWS

@@ -10,6 +10,63 @@ NOTE: The main developers do not have access to a Windoz based platform.
       Of course you are welcome to send us any patches that fix your
       Windoz build environment.
 
+=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
+=-= compiling with Cygwin =-=
+=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
+
+An effort is being made to allow windows users to compile calc using the
+Cygwin project (http://sources.redhat.com/cygwin/) with the GCC compiler
+and Un*x tools for Windows.
+
+The major porting work for Cygwin was performed by Thomas Jones-Low
+(tjoneslo at softstart dot com).
+
+In March 2009, Michael Penk (mpenk at wuska dot com) reported success in
+installs under Cygwin:
+
+    On my fairly complete Cygwin installs, everything compiles,
+    checks, and installs correctly.  My Cygwin is configured
+    in a very standard way (out of the box, using all of Cygwin's
+    defaults).  The install worked on 5 different machines with
+    Cygwin on them: one XP home, one XP professional, and three
+    Vista professionals.
+
+Using the calc Makefile, he did the following:
+
+    make all target=Cygwin
+    make check
+    make install
+
+He also reports:
+
+    Of course, one should be logged in as an Administrator when
+    one builds and installs calc.
+
+He was compiling calc 2.12.4.0 with Cygwin version 1.5.25-15.
+
+=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
+=-= If all else fails, for Cygwin =-=
+=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
+
+Much earlier (2001?) Thomas Jones-Low (tjoneslo at softstart dot com)
+recommended that you generate by hand all of the header files that
+by the Makefile.  This has been done for you via the makefile rule:
+
+	make win32_hsrc
+
+which uses the Makefile variables in win32.mkdef to form these header
+files under win32 directory.
+
+You will find generated versions of these files located in the win32
+sub-directory.  These files may be appropriate for your Cygwin building
+needs.
+
+In particular:
+
+	Just copy the win32/*.[ch] files up into the top level calc
+	source directory, edit them (if needed) and build using the
+	Cygwin GCC compiler and Cygwin build environment.
+
 =-=-=-=-=-=-=-=-=-=-=-=-=-=-=
 =-= compiling under DJGPP =-=
 =-=-=-=-=-=-=-=-=-=-=-=-=-=-=
@@ -68,52 +125,8 @@ Look for Makefile comments of the form:
 Follow those recommendations.  In cases where they conflict with
 the above Makefile list, follow the recommendation in the Makefile.
 
-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
-=-= compiling with Cygwin =-=
-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
 
-An effort is being made to allow windows users to compile calc using the
-Cygwin project (http://sources.redhat.com/cygwin/) with the GCC compiler
-and Un*x tools for Windows.
-
-The major porting work was performed by Thomas Jones-Low
-(tjoneslo at softstart dot com).  He said:
-
-    I had previous stated to this group that I have successfully managed
-    to port a version of Calc to Windows, and promised some point to
-    post what was required, so here it is.
-
-    One obvious manner of doing this port is to get the latest version
-    of the Cygwin project (http://sources.redhat.com/cygwin/) with the
-    GCC compiler and Un*x tools for Windows and recompile.
-
-    I built my working version using Calc ... I am using Visual C++
-    version 7.0, which is an older version of the Microsoft development
-    tools. The make file provided with Calc is not compatible with
-    NMAKE, so I used the Visual Studio tools to generate another one
-    (not included). Calc is built in two parts, calc.dll, which is the
-    library, and calc.exe which is the command line interface.
-
-He recommended that you generate by hand all of the header files that
-by the Makefile.  This has been done for you via the makefile rule:
-
-	make win32_hsrc
-
-which uses the Makefile variables in win32.mkdef to form these header
-files under win32 directory.
-
-You will find generated versions of these files located in the win32
-sub-directory.  These files may be appropriate for your Cygwin building
-needs.
-
-In particular:
-
-	Just copy the win32/*.[ch] files up into the top level calc
-	source directory, edit them (if needed) and build using the
-	Cygwin GCC compiler and Cygwin build environment.
-
-
-## Copyright (C) 2002-2007  Landon Curt Noll and Thomas Jones-Low
+## Copyright (C) 2002-2009  Landon Curt Noll and Thomas Jones-Low
 ##
 ## Calc is open software; you can redistribute it and/or modify it under
 ## the terms of the version 2.1 of the GNU Lesser General Public License
@@ -129,9 +142,9 @@ In particular:
 ## received a copy with calc; if not, write to Free Software Foundation, Inc.
 ## 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
 ##
-## @(#) $Revision: 30.1 $
-## @(#) $Id: README.WINDOWS,v 30.1 2007/03/16 11:09:46 chongo Exp $
-## @(#) $Source: /usr/local/src/cmd/calc/RCS/README.WINDOWS,v $
+## @(#) $Revision: 30.2 $
+## @(#) $Id: README.WINDOWS,v 30.2 2009/03/14 02:29:31 chongo Exp $
+## @(#) $Source: /usr/local/src/bin/calc/RCS/README.WINDOWS,v $
 ##
 ## Under source code control:	2001/02/25 14:00:05
 ## File existed as early as:	2001

addop.c

@@ -21,7 +21,7 @@
  *
  * @(#) $Revision: 30.1 $
  * @(#) $Id: addop.c,v 30.1 2007/03/16 11:09:46 chongo Exp $
- * @(#) $Source: /usr/local/src/cmd/calc/RCS/addop.c,v $
+ * @(#) $Source: /usr/local/src/bin/calc/RCS/addop.c,v $
  *
  * Under source code control:	1990/02/15 01:48:10
  * File existed as early as:	before 1990

align32.c

@@ -19,7 +19,7 @@
  *
  * @(#) $Revision: 30.1 $
  * @(#) $Id: align32.c,v 30.1 2007/03/16 11:09:46 chongo Exp $
- * @(#) $Source: /usr/local/src/cmd/calc/RCS/align32.c,v $
+ * @(#) $Source: /usr/local/src/bin/calc/RCS/align32.c,v $
  *
  * Under source code control:	1995/11/23 05:18:06
  * File existed as early as:	1995

alloc.h

@@ -17,9 +17,9 @@
  * received a copy with calc; if not, write to Free Software Foundation, Inc.
  * 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
  *
- * @(#) $Revision: 30.1 $
- * @(#) $Id: alloc.h,v 30.1 2007/03/16 11:09:46 chongo Exp $
- * @(#) $Source: /usr/local/src/cmd/calc/RCS/alloc.h,v $
+ * @(#) $Revision: 30.3 $
+ * @(#) $Id: alloc.h,v 30.3 2013/08/11 08:41:38 chongo Exp $
+ * @(#) $Source: /usr/local/src/bin/calc/RCS/alloc.h,v $
  *
  * Under source code control:	1990/02/15 01:48:29
  * File existed as early as:	before 1990
@@ -46,11 +46,15 @@
 # include <string.h>
 
 #else
+#if defined(_WIN32) && defined(NOTCYGWIN)
+#include <stdio.h>
+#endif
 
 # if defined(HAVE_NEWSTR)
 E_FUNC void *memcpy();
 E_FUNC void *memset();
-#if defined(FORCE_STDC) || (defined(__STDC__) && __STDC__ != 0) || defined(__cplusplus)
+#if defined(FORCE_STDC) || \
+   (defined(__STDC__) && __STDC__ != 0) || defined(__cplusplus)
 E_FUNC size_t strlen();
 #  else
 E_FUNC long strlen();
@@ -79,7 +83,8 @@ E_FUNC int strcmp();
 
 #if !defined(HAVE_MEMMOVE)
 # undef MEMMOVE_SIZE_T
-#if defined(FORCE_STDC) || (defined(__STDC__) && __STDC__ != 0) || defined(__cplusplus)
+#if defined(FORCE_STDC) || \
+   (defined(__STDC__) && __STDC__ != 0) || defined(__cplusplus)
 #  define MEMMOVE_SIZE_T size_t
 # else
 #  define MEMMOVE_SIZE_T long

assocfunc.c

@@ -17,9 +17,9 @@
  * received a copy with calc; if not, write to Free Software Foundation, Inc.
  * 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
  *
- * @(#) $Revision: 30.1 $
- * @(#) $Id: assocfunc.c,v 30.1 2007/03/16 11:09:46 chongo Exp $
- * @(#) $Source: /usr/local/src/cmd/calc/RCS/assocfunc.c,v $
+ * @(#) $Revision: 30.2 $
+ * @(#) $Id: assocfunc.c,v 30.2 2013/08/11 08:41:38 chongo Exp $
+ * @(#) $Source: /usr/local/src/bin/calc/RCS/assocfunc.c,v $
  *
  * Under source code control:	1993/07/20 23:04:27
  * File existed as early as:	1993
@@ -332,7 +332,8 @@ assoccopy(ASSOC *oldap)
 			oldep = oldep->e_next) {
 			ep = (ASSOCELEM *) malloc(ELEMSIZE(oldep->e_dim));
 			if (ep == NULL) {
-				math_error("Cannot allocate association element");
+				math_error("Cannot allocate "
+					   "association element");
 				/*NOTREACHED*/
 			}
 			ep->e_dim = oldep->e_dim;
@@ -340,7 +341,8 @@ assoccopy(ASSOC *oldap)
 			ep->e_value.v_type = V_NULL;
 			ep->e_value.v_subtype = V_NOSUBTYPE;
 			for (i = 0; i < ep->e_dim; i++)
-				copyvalue(&oldep->e_indices[i], &ep->e_indices[i]);
+				copyvalue(&oldep->e_indices[i],
+					  &ep->e_indices[i]);
 			copyvalue(&oldep->e_value, &ep->e_value);
 			listhead = &ap->a_table[ep->e_hash % ap->a_size];
 			ep->e_next = *listhead;

blkcpy.c

@@ -19,9 +19,9 @@
  * received a copy with calc; if not, write to Free Software Foundation, Inc.
  * 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
  *
- * @(#) $Revision: 30.1 $
- * @(#) $Id: blkcpy.c,v 30.1 2007/03/16 11:09:46 chongo Exp $
- * @(#) $Source: /usr/local/src/cmd/calc/RCS/blkcpy.c,v $
+ * @(#) $Revision: 30.2 $
+ * @(#) $Id: blkcpy.c,v 30.2 2013/08/11 08:41:38 chongo Exp $
+ * @(#) $Source: /usr/local/src/bin/calc/RCS/blkcpy.c,v $
  *
  * Under source code control:	1997/04/18 20:41:26
  * File existed as early as:	1997
@@ -374,7 +374,8 @@ copyblk2mat(BLOCK *blk, long ssi, long num, MATRIX *dmat, long dsi)
  * copymat2blk - copy matrix to block
  */
 int
-copymat2blk(MATRIX *smat, long ssi, long num, BLOCK *dblk, long dsi, BOOL noreloc)
+copymat2blk(MATRIX *smat, long ssi, long num, BLOCK *dblk, long dsi,
+	    BOOL noreloc)
 {
 	long i;
 	long newlen;
@@ -720,7 +721,8 @@ copystr2file(STRING *str, long ssi, long num, FILEID id, long dsi)
  * copyblk2blk - copy block to block
  */
 int
-copyblk2blk(BLOCK *sblk, long ssi, long num, BLOCK *dblk, long dsi, BOOL noreloc)
+copyblk2blk(BLOCK *sblk, long ssi, long num, BLOCK *dblk, long dsi,
+	    BOOL noreloc)
 {
 	long	newlen;
 	long	newsize;
@@ -762,7 +764,8 @@ copyblk2blk(BLOCK *sblk, long ssi, long num, BLOCK *dblk, long dsi, BOOL noreloc
  * copystr2blk - copy string to block
  */
 int
-copystr2blk(STRING *str, long ssi, long num, BLOCK *dblk, long dsi, BOOL noreloc)
+copystr2blk(STRING *str, long ssi, long num, BLOCK *dblk, long dsi,
+	    BOOL noreloc)
 {
 	long	len;
 	long	newlen;
@@ -982,7 +985,8 @@ memmove(void *s1, CONST void *s2, MEMMOVE_SIZE_T n)
  * copynum2blk - copy number numerator to block
  */
 int
-copynum2blk(NUMBER *snum, long ssi, long num, BLOCK *dblk, long dsi, BOOL noreloc)
+copynum2blk(NUMBER *snum, long ssi, long num, BLOCK *dblk, long dsi,
+	    BOOL noreloc)
 {
 	size_t	newlen;
 	size_t	newsize;
@@ -1033,7 +1037,8 @@ copynum2blk(NUMBER *snum, long ssi, long num, BLOCK *dblk, long dsi, BOOL norelo
  * copyblk2num - copy block to number
  */
 int
-copyblk2num(BLOCK *sblk, long ssi, long num, NUMBER *dnum, long dsi, NUMBER **res)
+copyblk2num(BLOCK *sblk, long ssi, long num, NUMBER *dnum, long dsi,
+	    NUMBER **res)
 {
 	size_t	newlen;
 	NUMBER *ret;		/* cloned and modified numerator */

blkcpy.h

@@ -21,7 +21,7 @@
  *
  * @(#) $Revision: 30.1 $
  * @(#) $Id: blkcpy.h,v 30.1 2007/03/16 11:09:46 chongo Exp $
- * @(#) $Source: /usr/local/src/cmd/calc/RCS/blkcpy.h,v $
+ * @(#) $Source: /usr/local/src/bin/calc/RCS/blkcpy.h,v $
  *
  * Under source code control:	1997/04/18 20:41:25
  * File existed as early as:	1997

block.c

@@ -21,7 +21,7 @@
  *
  * @(#) $Revision: 30.1 $
  * @(#) $Id: block.c,v 30.1 2007/03/16 11:09:46 chongo Exp $
- * @(#) $Source: /usr/local/src/cmd/calc/RCS/block.c,v $
+ * @(#) $Source: /usr/local/src/bin/calc/RCS/block.c,v $
  *
  * Under source code control:	1997/02/27 00:29:40
  * File existed as early as:	1997

block.h

@@ -21,7 +21,7 @@
  *
  * @(#) $Revision: 30.1 $
  * @(#) $Id: block.h,v 30.1 2007/03/16 11:09:46 chongo Exp $
- * @(#) $Source: /usr/local/src/cmd/calc/RCS/block.h,v $
+ * @(#) $Source: /usr/local/src/bin/calc/RCS/block.h,v $
  *
  * Under source code control:	1997/02/21 05:03:39
  * File existed as early as:	1997

byteswap.c

@@ -17,9 +17,9 @@
  * received a copy with calc; if not, write to Free Software Foundation, Inc.
  * 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
  *
- * @(#) $Revision: 30.1 $
- * @(#) $Id: byteswap.c,v 30.1 2007/03/16 11:09:46 chongo Exp $
- * @(#) $Source: /usr/local/src/cmd/calc/RCS/byteswap.c,v $
+ * @(#) $Revision: 30.2 $
+ * @(#) $Id: byteswap.c,v 30.2 2013/08/11 08:41:38 chongo Exp $
+ * @(#) $Source: /usr/local/src/bin/calc/RCS/byteswap.c,v $
  *
  * Under source code control:	1995/10/11 04:44:01
  * File existed as early as:	1995
@@ -102,7 +102,8 @@ swap_b8_in_ZVALUE(ZVALUE *dest, ZVALUE *src, BOOL all)
 		 */
 		dest = malloc(sizeof(ZVALUE));
 		if (dest == NULL) {
-			math_error("swap_b8_in_ZVALUE: swap_b8_in_ZVALUE: Not enough memory");
+			math_error("swap_b8_in_ZVALUE: swap_b8_in_ZVALUE: "
+				    "Not enough memory");
 			/*NOTREACHED*/
 		}
 

byteswap.h

@@ -19,7 +19,7 @@
  *
  * @(#) $Revision: 30.1 $
  * @(#) $Id: byteswap.h,v 30.1 2007/03/16 11:09:46 chongo Exp $
- * @(#) $Source: /usr/local/src/cmd/calc/RCS/byteswap.h,v $
+ * @(#) $Source: /usr/local/src/bin/calc/RCS/byteswap.h,v $
  *
  * Under source code control:	1995/10/11 04:44:01
  * File existed as early as:	1995

cal/Makefile

@@ -18,9 +18,9 @@
 # received a copy with calc; if not, write to Free Software Foundation, Inc.
 # 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
 #
-# @(#) $Revision: 30.2 $
-# @(#) $Id: Makefile,v 30.2 2007/03/22 07:34:47 chongo Exp $
-# @(#) $Source: /usr/local/src/cmd/calc/cal/RCS/Makefile,v $
+# @(#) $Revision: 30.10 $
+# @(#) $Id: Makefile,v 30.10 2013/09/02 03:02:00 chongo Exp $
+# @(#) $Source: /usr/local/src/bin/calc/cal/RCS/Makefile,v $
 #
 # Under source code control:	1991/07/21 05:00:54
 # File existed as early as:	1991
@@ -111,9 +111,6 @@ CALC_SHAREDIR= /usr/share/calc
 # ---------------------------------------------------------------
 # ${HELPDIR}		where the help directory is installed
 # ${CALC_INCDIR}	where the calc include files are installed
-# ${CUSTOMCALDIR}	where custom *.cal files are installed
-# ${CUSTOMHELPDIR}	where custom help files are installed
-# ${CUSTOMINCDIR}	where custom .h files are installed
 # ${SCRIPTDIR}		where calc shell scripts are installed
 #
 # NOTE: The install rule prepends installation paths with ${T}, which
@@ -124,16 +121,10 @@ CALC_SHAREDIR= /usr/share/calc
 #
 #	HELPDIR= ${CALC_SHAREDIR}/help
 #	CALC_INCDIR= ${INCDIR}/calc
-#	CUSTOMCALDIR= ${CALC_SHAREDIR}/custom
-#	CUSTOMHELPDIR= ${CALC_SHAREDIR}/custhelp
-#	CUSTOMINCDIR= ${CALC_INCDIR}/custom
 #	SCRIPTDIR= ${BINDIR}/cscript
 #
 HELPDIR= ${CALC_SHAREDIR}/help
 CALC_INCDIR= ${INCDIR}/calc
-CUSTOMCALDIR= ${CALC_SHAREDIR}/custom
-CUSTOMHELPDIR= ${CALC_SHAREDIR}/custhelp
-CUSTOMINCDIR= ${CALC_INCDIR}/custom
 SCRIPTDIR= ${BINDIR}/cscript
 
 # T - top level directory under which calc will be installed
@@ -188,21 +179,37 @@ MV= mv
 CO= co
 TRUE= true
 TOUCH= touch
+SED= sed
+SORT= sort
+FMT= fmt
 
 # The calc files to install
 #
-CALC_FILES= README bigprime.cal deg.cal ellip.cal lucas.cal lucas_chk.cal \
-	lucas_tbl.cal mersenne.cal mod.cal pell.cal pi.cal pix.cal \
-	pollard.cal poly.cal psqrt.cal quat.cal regress.cal solve.cal \
-	sumsq.cal surd.cal unitfrac.cal varargs.cal chrem.cal mfactor.cal \
-	bindings randmprime.cal test1700.cal randrun.cal linear.cal \
-	randbitrun.cal bernoulli.cal test2300.cal test2600.cal \
-	test2700.cal test3100.cal test3300.cal test3400.cal prompt.cal \
-	test3500.cal seedrandom.cal test4000.cal test4100.cal test4600.cal \
-	beer.cal hello.cal test5100.cal test5200.cal randombitrun.cal \
-	randomrun.cal repeat.cal xx_print.cal natnumset.cal qtime.cal \
-	test8400.cal test8500.cal test8600.cal chi.cal intfile.cal screen.cal \
-	dotest.cal set8700.cal set8700.line alg_config.cal sumtimes.cal
+# This list is prodiced by the detaillist rule when no WARNINGS are detected.
+#
+# Please use:
+#
+#	make calc_files_list
+#
+# to keep this list in nice sorted order and to check that these
+# deailed help files are under RCS control.
+#
+CALC_FILES= README alg_config.cal beer.cal bernoulli.cal \
+	bernpoly.cal bigprime.cal bindings brentsolve.cal chi.cal chrem.cal \
+	constants.cal deg.cal dms.cal dotest.cal ellip.cal factorial.cal \
+	factorial2.cal gvec.cal hello.cal hms.cal infinities.cal \
+	intfile.cal intnum.cal lambertw.cal linear.cal lnseries.cal \
+	lucas.cal lucas_chk.cal lucas_tbl.cal mersenne.cal mfactor.cal \
+	mod.cal natnumset.cal pell.cal pi.cal pix.cal pollard.cal poly.cal \
+	prompt.cal psqrt.cal qtime.cal quat.cal randbitrun.cal randmprime.cal \
+	randombitrun.cal randomrun.cal randrun.cal regress.cal repeat.cal \
+	screen.cal seedrandom.cal set8700.cal set8700.line smallfactors.cal \
+	solve.cal specialfunctions.cal statistics.cal strings.cal sumsq.cal \
+	sumtimes.cal surd.cal test1700.cal test2300.cal test2600.cal \
+	test2700.cal test3100.cal test3300.cal test3400.cal test3500.cal \
+	test4000.cal test4100.cal test4600.cal test5100.cal test5200.cal \
+	test8400.cal test8500.cal test8600.cal test8900.cal toomcook.cal \
+	unitfrac.cal varargs.cal xx_print.cal zeta2.cal
 
 # These files are found (but not built) in the distribution
 #
@@ -250,6 +257,27 @@ calcliblist:
 	    fi; \
 	done
 
+# These next rule help form the ${CALC_FILES} makefile variables above.
+#
+calc_files_list:
+	${Q} -(find .  -mindepth 1 -maxdepth 1 -type f -name '*.cal' -print | \
+	   while read i; do \
+		if [ X"$$i" != X"/dev/null" ]; then \
+		    if [ ! -f RCS/$$i,v ]; then \
+			echo "WARNING: $$i not under RCS control" 1>&2; \
+		    else \
+			echo $$i; \
+		    fi; \
+		fi; \
+	done; \
+	echo '--first_line--'; \
+	echo README; \
+	echo set8700.line; \
+	echo bindings) | \
+	    ${SED} -e 's:^\./::' | LANG=C ${SORT} | ${FMT} -70 | \
+	    ${SED} -e '1s/--first_line--/CALC_FILES=/' -e '2,$$s/^/	/' \
+		   -e 's/$$/ \\/' -e '$$s/ \\$$//'
+
 ##
 #
 # rpm rules
@@ -301,7 +329,7 @@ install: all
 		${RM} -f ${T}${CALC_SHAREDIR}/$$i.new; \
 		${CP} -f $$i ${T}${CALC_SHAREDIR}/$$i.new; \
 		${CHMOD} 0444 ${T}${CALC_SHAREDIR}/$$i.new; \
-		${MV} -f ${T}${CALC_SHAREDIR}/$$i.new ${T}${CALC_SHAREDIR}/$$i; \
+		${MV} -f ${T}${CALC_SHAREDIR}/$$i.new ${T}${CALC_SHAREDIR}/$$i;\
 		echo "installed ${T}${CALC_SHAREDIR}/$$i"; \
 	    fi; \
 	done

cal/alg_config.cal

@@ -17,9 +17,9 @@
  * received a copy with calc; if not, write to Free Software Foundation, Inc.
  * 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
  *
- * @(#) $Revision: 30.1 $
- * @(#) $Id: alg_config.cal,v 30.1 2007/03/16 11:09:54 chongo Exp $
- * @(#) $Source: /usr/local/src/cmd/calc/cal/RCS/alg_config.cal,v $
+ * @(#) $Revision: 30.2 $
+ * @(#) $Id: alg_config.cal,v 30.2 2013/08/11 08:41:38 chongo Exp $
+ * @(#) $Source: /usr/local/src/bin/calc/cal/RCS/alg_config.cal,v $
  *
  * Under source code control:	2006/06/07 14:10:11
  * File existed as early as:	2006
@@ -80,7 +80,8 @@ define mul_loop(repeat, x)
     len = sizeof((*x)[0]) / baseb_bytes;
     for (i=1; i < 4; ++i) {
 	if ((sizeof((*x)[i]) / baseb_bytes) != len) {
-	    quit "mul_loop: 2nd arg matrix elements are not of equal BASEB-bit word length";
+	    quit "mul_loop: 2nd arg matrix elements are not of "
+	         "equal BASEB-bit word length";
 	}
     }
 
@@ -378,7 +379,8 @@ define best_mul2()
 	}
     } while (ratio >= 1.0);
     if (config("user_debug") > 0) {
-	printf("alg1/alg2 ratio now < 1.0, starting binary search between %d and %d\n",
+	printf("alg1/alg2 ratio now < 1.0, starting binary search "
+	       "between %d and %d\n",
 	    low, high);
     }
 
@@ -472,7 +474,8 @@ define sq_loop(repeat, x)
     len = sizeof((*x)[0]) / baseb_bytes;
     for (i=1; i < 4; ++i) {
 	if ((sizeof((*x)[i]) / baseb_bytes) != len) {
-	    quit "sq_loop: 2nd arg matrix elements are not of equal BASEB-bit word length";
+	    quit "sq_loop: 2nd arg matrix elements are not of equal "
+	         "BASEB-bit word length";
 	}
     }
 
@@ -769,7 +772,8 @@ define best_sq2()
 	}
     } while (ratio >= 1.0);
     if (config("user_debug") > 0) {
-	printf("alg1/alg2 ratio now < 1.0, starting binary search between %d and %d\n",
+	printf("alg1/alg2 ratio now < 1.0, starting binary search "
+	        "between %d and %d\n",
 	    low, high);
     }
 
@@ -866,7 +870,8 @@ define pow_loop(repeat, x, ex)
     len = sizeof((*x)[0]) / baseb_bytes;
     for (i=1; i < 4; ++i) {
 	if ((sizeof((*x)[i]) / baseb_bytes) != len) {
-	    quit "pow_loop: 2nd arg matrix elements are not of equal BASEB-bit word length";
+	    quit "pow_loop: 2nd arg matrix elements are not of "
+	         "equal BASEB-bit word length";
 	}
     }
     if (!isint(ex) || ex < 3) {
@@ -1151,7 +1156,8 @@ define best_pow2()
 	if (config("user_debug") > 1) {
 	    printf("    pmod alg1/alg2 ratio = %.3f\n", ratio);
 	    if (ratio > 1.0 && ratio <= 1.02) {
-	    printf("    while alg1 is slightly better than alg2, it is not clearly better\n");
+	    printf("    while alg1 is slightly better than alg2, "
+	    	   "it is not clearly better\n");
 	    }
 	}
     } while (ratio <= 1.02);
@@ -1205,8 +1211,8 @@ define best_pow2()
 	looped = 1;
     } while (ratio >= 1.0);
     if (config("user_debug") > 0) {
-	printf("alg1/alg2 ratio now < 1.0, starting binary search between %d and %d\n",
-	    low, high);
+	printf("alg1/alg2 ratio now < 1.0, starting binary search "
+	       "between %d and %d\n", low, high);
     }
 
     /*

cal/beer.cal

@@ -19,7 +19,7 @@
  *
  * @(#) $Revision: 30.1 $
  * @(#) $Id: beer.cal,v 30.1 2007/03/16 11:09:54 chongo Exp $
- * @(#) $Source: /usr/local/src/cmd/calc/cal/RCS/beer.cal,v $
+ * @(#) $Source: /usr/local/src/bin/calc/cal/RCS/beer.cal,v $
  *
  * Under source code control:	1996/11/13 13:21:05
  * File existed as early as:	1996

cal/bernoulli.cal

@@ -19,7 +19,7 @@
  *
  * @(#) $Revision: 30.1 $
  * @(#) $Id: bernoulli.cal,v 30.1 2007/03/16 11:09:54 chongo Exp $
- * @(#) $Source: /usr/local/src/cmd/calc/cal/RCS/bernoulli.cal,v $
+ * @(#) $Source: /usr/local/src/bin/calc/cal/RCS/bernoulli.cal,v $
  *
  * Under source code control:	1991/09/30 11:18:41
  * File existed as early as:	1991

cal/bernpoly.cal

@@ -0,0 +1,59 @@
+/*
+ * bernpoly - Bernoully polynomials B_n(z) for arbitrary n,z..
+ *
+ * Copyright (C) 2013 Christoph Zurnieden
+ *
+ * Calc is open software; you can redistribute it and/or modify it under
+ * the terms of the version 2.1 of the GNU Lesser General Public License
+ * as published by the Free Software Foundation.
+ *
+ * Calc is distributed in the hope that it will be useful, but WITHOUT
+ * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
+ * or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General
+ * Public License for more details.
+ *
+ * A copy of version 2.1 of the GNU Lesser General Public License is
+ * distributed with calc under the filename COPYING-LGPL.  You should have
+ * received a copy with calc; if not, write to Free Software Foundation, Inc.
+ * 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
+ *
+ * @(#) $Revision: 30.4 $
+ * @(#) $Id: bernpoly.cal,v 30.4 2013/08/18 20:01:53 chongo Exp $
+ * @(#) $Source: /usr/local/src/bin/calc/cal/RCS/bernpoly.cal,v $
+ *
+ * Under source code control:	2013/08/11 01:31:28
+ * File existed as early as:	2013
+ */
+
+
+static resource_debug_level;
+resource_debug_level = config("resource_debug", 0);
+
+
+read -once zeta2
+
+
+/* Idea by Don Zagier */
+define bernpoly(n,z){
+  local h s c k;
+  if(isint(n) && n>=0){
+    h=0;s=0;c=-1;
+    for(k=1;k<=n+1;k++){
+      c*=1-(n+2)/k;
+      s+=z^n;
+      z++;
+      h+=c*s/k;
+    }
+    return h;
+  }
+  else return -n*hurwitzzeta(1-n,z);
+}
+
+
+/*
+ * restore internal function from resource debugging
+ */
+config("resource_debug", resource_debug_level),;
+if (config("resource_debug") & 3) {
+    print "bernpoly(n,z)";
+}

cal/bigprime.cal

@@ -19,7 +19,7 @@
  *
  * @(#) $Revision: 30.1 $
  * @(#) $Id: bigprime.cal,v 30.1 2007/03/16 11:09:54 chongo Exp $
- * @(#) $Source: /usr/local/src/cmd/calc/cal/RCS/bigprime.cal,v $
+ * @(#) $Source: /usr/local/src/bin/calc/cal/RCS/bigprime.cal,v $
  *
  * Under source code control:	1991/05/22 21:56:32
  * File existed as early as:	1991

cal/bindings

@@ -18,7 +18,7 @@
 #
 # @(#) $Revision: 30.1 $
 # @(#) $Id: bindings,v 30.1 2007/03/16 11:09:54 chongo Exp $
-# @(#) $Source: /usr/local/src/cmd/calc/cal/RCS/bindings,v $
+# @(#) $Source: /usr/local/src/bin/calc/cal/RCS/bindings,v $
 #
 # Under source code control:	1993/05/02 20:09:19
 # File existed as early as:	1993

cal/brentsolve.cal

@@ -0,0 +1,258 @@
+/*
+ * brentsolve - Root finding with the Brent-Dekker trick
+ *
+ * Copyright (C) 2013 Christoph Zurnieden
+ *
+ * Calc is open software; you can redistribute it and/or modify it under
+ * the terms of the version 2.1 of the GNU Lesser General Public License
+ * as published by the Free Software Foundation.
+ *
+ * Calc is distributed in the hope that it will be useful, but WITHOUT
+ * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
+ * or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General
+ * Public License for more details.
+ *
+ * A copy of version 2.1 of the GNU Lesser General Public License is
+ * distributed with calc under the filename COPYING-LGPL.  You should have
+ * received a copy with calc; if not, write to Free Software Foundation, Inc.
+ * 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
+ *
+ * @(#) $Revision: 30.4 $
+ * @(#) $Id: brentsolve.cal,v 30.4 2013/08/18 20:01:53 chongo Exp $
+ * @(#) $Source: /usr/local/src/bin/calc/cal/RCS/brentsolve.cal,v $
+ *
+ * Under source code control:	2013/08/11 01:31:28
+ * File existed as early as:	2013
+ */
+
+
+static resource_debug_level;
+resource_debug_level = config("resource_debug", 0);
+
+
+/*
+  A short explanation is at http://en.wikipedia.org/wiki/Brent%27s_method
+  I tried to follow the description at wikipedia as much as possible to make
+  the slight changes I did more visible.
+  You may give http://people.sc.fsu.edu/~jburkardt/cpp_src/brent/brent.html a
+  short glimpse (Brent's originl Fortran77 versions and some translations of
+  it).
+*/
+
+static true = 1;
+static false = 0;
+define brentsolve(low, high,eps){
+  local a b c d fa fb fc fa2 fb2 fc2 s fs  tmp tmp2  mflag i places;
+  a = low;
+  b = high;
+  c = 0;
+
+  if(isnull(eps))
+    eps = epsilon(epsilon()*1e-3);
+  places = highbit(1 + int( 1/epsilon() ) ) + 1;
+
+  d = 1/eps;
+
+  fa = f(a);
+  fb = f(b);
+
+  fc = 0;
+  s = 0;
+  fs = 0;
+
+  if(fa * fb >= 0){
+    if(fa < fb){
+      epsilon(eps);
+      return a;
+    }
+    else{
+      epsilon(eps);
+      return b;
+    }
+  }
+
+  if(abs(fa) < abs(fb)){
+    tmp = a; a = b; b = tmp;
+    tmp = fa; fa = fb; fb = tmp;
+  }
+
+  c = a;
+  fc = fa;
+  mflag = 1;
+  i = 0;
+
+  while(!(fb==0) && (abs(a-b) > eps)){
+    if((fa != fc) && (fb != fc)){
+      /* Inverse quadratic interpolation*/
+      fc2 = fc^2;
+      fa2 = fa^2;
+      s = bround(((fb^2*((fc*a)-(c*fa)))+(fb*((c*fa2)-(fc2*a)))+(b*((fc2*fa)
+                    -(fc*fa2))))/((fc - fb)*(fa - fb)*(fc - fa)),places++);
+    }
+    else{
+      /* Secant Rule*/
+      s =bround( b - fb * (b - a) / (fb - fa),places++);
+    }
+    tmp2 = (3 * a + b) / 4;
+    if( (!( ((s > tmp2) && (s < b))||((s < tmp2) && (s > b))))
+        || (mflag && (abs(s - b) >= (abs(b - c) / 2)))
+        || (!mflag && (abs(s - b) >= (abs(c - d) / 2)))) {
+      s = (a + b) / 2;
+      mflag = true;
+    }
+    else{
+      if( (mflag && (abs(b - c) < eps))
+          || (!mflag && (abs(c - d) < eps))) {
+        s = (a + b) / 2;
+        mflag = true;
+      }
+      else
+        mflag = false;
+    }
+    fs = f(s);
+    c = b;
+    fc = fb;
+    if (fa * fs < 0){
+      b = s;
+      fb = fs;
+    }
+    else {
+      a = s;
+      fa = fs;
+    }
+
+    if (abs(fa) < abs(fb)){
+      tmp = a; a = b; b = tmp;
+      tmp = fa; fa = fb; fb = tmp;
+    }
+    i++;
+    if (i > 1000){
+      epsilon(eps);
+      return newerror("brentsolve: does not converge");
+    }
+  }
+  epsilon(eps);
+  return b;
+}
+
+/*
+  A variation of the solver to accept functions named differently from "f". The
+  code should explain it.
+*/
+define brentsolve2(low, high,which,eps){
+  local a b c d fa fb fc fa2 fb2 fc2 s fs  tmp tmp2  mflag i places;
+  a = low;
+  b = high;
+  c = 0;
+
+  switch(param(0)){
+    case 0:
+    case 1: return newerror("brentsolve2: not enough argments");
+    case 2: eps = epsilon(epsilon()*1e-2);
+            which = 0;break;
+    case 3: eps = epsilon(epsilon()*1e-2);break;
+    default: break;
+  };
+  places = highbit(1 + int(1/epsilon())) + 1;
+
+  d = 1/eps;
+
+  switch(which){
+    case 1:  fa = __CZ__invbeta(a);
+             fb = __CZ__invbeta(b); break;
+    case 2:  fa = __CZ__invincgamma(a);
+             fb = __CZ__invincgamma(b); break;
+
+    default: fa = f(a);fb = f(b);   break;
+  };
+
+  fc = 0;
+  s = 0;
+  fs = 0;
+
+  if(fa * fb >= 0){
+    if(fa < fb)
+      return a;
+    else
+      return b;
+  }
+
+  if(abs(fa) < abs(fb)){
+    tmp = a; a = b; b = tmp;
+    tmp = fa; fa = fb; fb = tmp;
+  }
+
+  c = a;
+  fc = fa;
+  mflag = 1;
+  i = 0;
+
+  while(!(fb==0) && (abs(a-b) > eps)){
+
+    if((fa != fc) && (fb != fc)){
+      /* Inverse quadratic interpolation*/
+      fc2 = fc^2;
+      fa2 = fa^2;
+      s = bround(((fb^2*((fc*a)-(c*fa)))+(fb*((c*fa2)-(fc2*a)))+(b*((fc2*fa)
+                    -(fc*fa2))))/((fc - fb)*(fa - fb)*(fc - fa)),places);
+      places++;
+    }
+    else{
+      /* Secant Rule*/
+      s =bround( b - fb * (b - a) / (fb - fa),places);
+      places++;
+    }
+    tmp2 = (3 * a + b) / 4;
+    if( (!( ((s > tmp2) && (s < b))||((s < tmp2) && (s > b))))
+        || (mflag && (abs(s - b) >= (abs(b - c) / 2)))
+        || (!mflag && (abs(s - b) >= (abs(c - d) / 2)))) {
+      s = (a + b) / 2;
+      mflag = true;
+    }
+    else{
+      if( (mflag && (abs(b - c) < eps))
+          || (!mflag && (abs(c - d) < eps))) {
+        s = (a + b) / 2;
+        mflag = true;
+      }
+      else
+        mflag = false;
+    }
+    switch(which){
+      case 1:  fs = __CZ__invbeta(s); break;
+      case 2:  fs = __CZ__invincgamma(s); break;
+
+      default: fs = f(s);             break;
+    };
+    c = b;
+    fc = fb;
+    if (fa * fs < 0){
+      b = s;
+      fb = fs;
+    }
+    else {
+      a = s;
+      fa = fs;
+    }
+
+    if (abs(fa) < abs(fb)){
+      tmp = a; a = b; b = tmp;
+      tmp = fa; fa = fb; fb = tmp;
+    }
+    i++;
+    if (i > 1000){
+      return newerror("brentsolve2: does not converge");
+    }
+  }
+  return b;
+}
+
+
+/*
+ * restore internal function from resource debugging
+ */
+config("resource_debug", resource_debug_level),;
+if (config("resource_debug") & 3) {
+    print "brentsolve(low, high,eps)";
+    print "brentsolve2(low, high,which,eps)";
+}

cal/chi.cal

@@ -19,7 +19,7 @@
  *
  * @(#) $Revision: 30.1 $
  * @(#) $Id: chi.cal,v 30.1 2007/03/16 11:09:54 chongo Exp $
- * @(#) $Source: /usr/local/src/cmd/calc/cal/RCS/chi.cal,v $
+ * @(#) $Source: /usr/local/src/bin/calc/cal/RCS/chi.cal,v $
  *
  * Under source code control:	2001/03/27 14:10:11
  * File existed as early as:	2001

cal/chrem.cal

@@ -21,7 +21,7 @@
  *
  * @(#) $Revision: 30.1 $
  * @(#) $Id: chrem.cal,v 30.1 2007/03/16 11:09:54 chongo Exp $
- * @(#) $Source: /usr/local/src/cmd/calc/cal/RCS/chrem.cal,v $
+ * @(#) $Source: /usr/local/src/bin/calc/cal/RCS/chrem.cal,v $
  *
  * Under source code control:	1992/09/26 01:00:47
  * File existed as early as:	1992

cal/constants.cal

@@ -0,0 +1,104 @@
+/*
+ * constants - implementation of different constants to arbitrary precision
+ *
+ * Copyright (C) 2013 Christoph Zurnieden
+ *
+ * Calc is open software; you can redistribute it and/or modify it under
+ * the terms of the version 2.1 of the GNU Lesser General Public License
+ * as published by the Free Software Foundation.
+ *
+ * Calc is distributed in the hope that it will be useful, but WITHOUT
+ * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
+ * or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General
+ * Public License for more details.
+ *
+ * A copy of version 2.1 of the GNU Lesser General Public License is
+ * distributed with calc under the filename COPYING-LGPL.  You should have
+ * received a copy with calc; if not, write to Free Software Foundation, Inc.
+ * 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
+ *
+ * @(#) $Revision: 30.4 $
+ * @(#) $Id: constants.cal,v 30.4 2013/08/18 20:01:53 chongo Exp $
+ * @(#) $Source: /usr/local/src/bin/calc/cal/RCS/constants.cal,v $
+ *
+ * Under source code control:	2013/08/11 01:31:28
+ * File existed as early as:	2013
+ */
+
+
+static resource_debug_level;
+resource_debug_level = config("resource_debug", 0);
+
+
+static __CZ__euler_mascheroni = 0;
+static __CZ__euler_mascheroni_prec = 0;
+
+
+define e(){
+  local k temp1 temp2 ret eps factor upperlimit prec;
+
+  prec = digits(1/epsilon());
+  if(__CZ__euler_mascheroni != 0 &&  __CZ__euler_mascheroni_prec >= prec)
+    return __CZ__euler_mascheroni;
+  if(prec<=20) return 2.718281828459045235360287471;
+  if(prec<=1800){
+    __CZ__euler_mascheroni = exp(1);
+    __CZ__euler_mascheroni_prec = prec;
+  }
+
+  eps=epsilon(1e-20);
+  factor = 1;
+  k = 0;
+  upperlimit = prec * ln(10);
+  while(k<upperlimit){
+    k += ln(factor);
+    factor++;
+  }
+  epsilon(eps);
+  temp1 = 0;
+  ret   = 1;
+  for(k=3;k<=factor;k++){
+    temp2 = temp1;
+    temp1 = ret;
+    ret =  (k-1) *(temp1 + temp2);
+  }
+
+  ret = inverse( ret * inverse(factorial(factor) ) ) ;
+  __CZ__euler_mascheroni = ret;
+  __CZ__euler_mascheroni_prec = prec;
+  return ret;
+}
+
+
+/* Lupas' series */
+static __CZ__catalan = 0;
+static __CZ__catalan_prec = 0;
+define G(){
+  local eps a s t n;
+  eps = epsilon(epsilon()*1e-10);
+  if(__CZ__catalan != 0 &&  __CZ__catalan >= log(1/eps))
+    return __CZ__catalan;
+  a = 1;
+  s = 0;
+  t = 1;
+  n = 1;
+  while(abs(t)> eps){
+    a *= 32 * n^3 * (2*n-1);
+    a /=((3-16*n+16*n^2)^2);
+    t = a * (-1)^(n-1) * (40*n^2-24*n+3) / (n^3 * (2*n-1));
+    s += t;
+    n += 1;
+  }
+  s = s/64;
+  __CZ__catalan = s;
+  __CZ__catalan_prec = log(1/eps);
+  epsilon(eps);
+  return s;
+}
+
+
+config("resource_debug", resource_debug_level),;
+if (config("resource_debug") & 3) {
+    print "e()";
+    print "G()";
+}

cal/deg.cal

@@ -17,9 +17,9 @@
  * received a copy with calc; if not, write to Free Software Foundation, Inc.
  * 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
  *
- * @(#) $Revision: 30.1 $
- * @(#) $Id: deg.cal,v 30.1 2007/03/16 11:09:54 chongo Exp $
- * @(#) $Source: /usr/local/src/cmd/calc/cal/RCS/deg.cal,v $
+ * @(#) $Revision: 30.2 $
+ * @(#) $Id: deg.cal,v 30.2 2010/09/02 06:01:14 chongo Exp $
+ * @(#) $Source: /usr/local/src/bin/calc/cal/RCS/deg.cal,v $
  *
  * Under source code control:	1990/02/15 01:50:33
  * File existed as early as:	before 1990
@@ -28,9 +28,9 @@
  */
 
 
-obj dms {deg, min, sec};
+obj deg {deg, min, sec};
 
-define dms(deg, min, sec)
+define deg(deg, min, sec)
 {
 	local ans;
 
@@ -38,18 +38,18 @@ define dms(deg, min, sec)
 		sec = 0;
 	if (isnull(min))
 		min = 0;
-	obj dms ans;
+	obj deg ans;
 	ans.deg = deg;
 	ans.min = min;
 	ans.sec = sec;
-	fixdms(ans);
+	fixdeg(ans);
 	return ans;
 }
 
 
-define dms_add(a, b)
+define deg_add(a, b)
 {
-	local obj dms	ans;
+	local obj deg	ans;
 
 	ans.deg = 0;
 	ans.min = 0;
@@ -66,14 +66,14 @@ define dms_add(a, b)
 		ans.sec += b.sec;
 	} else
 		ans.deg += b;
-	fixdms(ans);
+	fixdeg(ans);
 	return ans;
 }
 
 
-define dms_neg(a)
+define deg_neg(a)
 {
-	local obj dms	ans;
+	local obj deg	ans;
 
 	ans.deg = -a.deg;
 	ans.min = -a.min;
@@ -82,15 +82,15 @@ define dms_neg(a)
 }
 
 
-define dms_sub(a, b)
+define deg_sub(a, b)
 {
 	return a - b;
 }
 
 
-define dms_mul(a, b)
+define deg_mul(a, b)
 {
-	local obj dms	ans;
+	local obj deg	ans;
 
 	if (istype(a, ans) && istype(b, ans))
 		quit "Cannot multiply degrees together";
@@ -103,24 +103,24 @@ define dms_mul(a, b)
 		ans.min = b.min * a;
 		ans.sec = b.sec * a;
 	}
-	fixdms(ans);
+	fixdeg(ans);
 	return ans;
 }
 
 
-define dms_print(a)
+define deg_print(a)
 {
 	print a.deg : 'd' : a.min : 'm' : a.sec : 's' :;
 }
 
 
-define dms_abs(a)
+define deg_abs(a)
 {
 	return a.deg + a.min / 60 + a.sec / 3600;
 }
 
 
-define fixdms(a)
+define fixdeg(a)
 {
 	a.min += frac(a.deg) * 60;
 	a.deg = int(a.deg);
@@ -134,5 +134,5 @@ define fixdms(a)
 }
 
 if (config("resource_debug") & 3) {
-    print "obj dms {deg, min, sec} defined";
+    print "obj deg {deg, min, sec} defined";
 }

cal/dms.cal

@@ -0,0 +1,384 @@
+/*
+ * dms - calculate in degrees, minutes, and seconds (based on deg)
+ *
+ * Copyright (C) 1999,2010  David I. Bell and Landon Curt Noll
+ *
+ * Calc is open software; you can redistribute it and/or modify it under
+ * the terms of the version 2.1 of the GNU Lesser General Public License
+ * as published by the Free Software Foundation.
+ *
+ * Calc is distributed in the hope that it will be useful, but WITHOUT
+ * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
+ * or FITNESS FOR A PARTICULAR PURPOSE.	 See the GNU Lesser General
+ * Public License for more details.
+ *
+ * A copy of version 2.1 of the GNU Lesser General Public License is
+ * distributed with calc under the filename COPYING-LGPL.  You should have
+ * received a copy with calc; if not, write to Free Software Foundation, Inc.
+ * 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
+ *
+ * @(#) $Revision: 30.2 $
+ * @(#) $Id: dms.cal,v 30.2 2010/09/02 06:14:16 chongo Exp $
+ * @(#) $Source: /usr/local/src/bin/calc/cal/RCS/dms.cal,v $
+ *
+ * Under source code control:	1990/02/15 01:50:33
+ * File existed as early as:	before 1990
+ *
+ * Share and enjoy!  :-)	http://www.isthe.com/chongo/tech/comp/calc/
+ */
+
+
+obj dms {deg, min, sec};
+
+define dms(deg, min, sec)
+{
+    local obj dms ans;		/* return value */
+
+    /* default missing args to 0 */
+    if (isnull(sec)) {
+	sec = 0;
+    }
+    if (isnull(min)) {
+	min = 0;
+    }
+
+    /* load object */
+    ans.deg = deg;
+    ans.min = min;
+    ans.sec = sec;
+
+    /* return properly formed object */
+    ans = fixdms(ans);
+    return ans;
+}
+
+
+define dms_add(a, b)
+{
+    local obj dms ans;		/* return value */
+
+    /* initalize value to 1st arg */
+    if (istype(a, ans)) {
+	/* 1st arg is dms object, load it */
+	ans.deg = a.deg;
+	ans.min = a.min;
+	ans.sec = a.sec;
+    } else {
+	/* 1st arg is not dms, assume scalar degrees */
+	ans.deg = a;
+	ans.min = 0;
+	ans.sec = 0;
+    }
+
+    /* add value of 2nd arg */
+    if (istype(b, ans)) {
+	/* 2nd arg is dms object, add it */
+	ans.deg += b.deg;
+	ans.min += b.min;
+	ans.sec += b.sec;
+    } else {
+	/* 2nd arg is not dms, add scalar degrees */
+	ans.deg += b;
+    }
+
+    /* return normalized result */
+    ans = fixdms(ans);
+    return ans;
+}
+
+
+define dms_neg(a)
+{
+    local obj dms ans;		/* return value */
+
+    /* negate argument */
+    if (istype(a, ans)) {
+	/* 1st arg is dms object, load it */
+	ans.deg = -a.deg;
+	ans.min = -a.min;
+	ans.sec = -a.sec;
+    } else {
+	/* 2nd arg is not dms, negate scalar degrees */
+	ans.deg = -a;
+	ans.min = 0;
+	ans.sec = 0;
+    }
+
+    /* return normalized result */
+    ans = fixdms(ans);
+    return ans;
+}
+
+
+define dms_sub(a, b)
+{
+    local obj dms ans;		/* return value */
+
+    /* initalize value to 1st arg */
+    if (istype(a, ans)) {
+	/* 1st arg is dms object, load it */
+	ans.deg = a.deg;
+	ans.min = a.min;
+	ans.sec = a.sec;
+    } else {
+	/* 1st arg is not dms, assume scalar degrees */
+	ans.deg = a;
+	ans.min = 0;
+	ans.sec = 0;
+    }
+
+    /* subtract value of 2nd arg */
+    if (istype(b, ans)) {
+	/* 2nd arg is dms object, subtract it */
+	ans.deg -= b.deg;
+	ans.min -= b.min;
+	ans.sec -= b.sec;
+    } else {
+	/* 2nd arg is not dms, subtract scalar degrees */
+	ans.deg -= b;
+    }
+
+    /* return normalized result */
+    ans = fixdms(ans);
+    return ans;
+}
+
+
+define dms_mul(a, b)
+{
+    local obj dms ans;		/* return value */
+
+    /* dms object multiplication */
+    if (istype(a, ans) && istype(b, ans)) {
+	ans.deg = dms_abs(a) * dms_abs(b);
+	ans.min = 0;
+	ans.sec = 0;
+
+    /* scalar multiplication */
+    } else if (istype(a, ans)) {
+	ans.deg = a.deg * b;
+	ans.min = a.min * b;
+	ans.sec = a.sec * b;
+    } else {
+	ans.deg = b.deg * a;
+	ans.min = b.min * a;
+	ans.sec = b.sec * a;
+    }
+
+    /* return normalized result */
+    ans = fixdms(ans);
+    return ans;
+}
+
+
+define dms_print(a)
+{
+    local obj dms ans;		/* temp object for dms type testing */
+
+    /* firewall - arg must be a dms object */
+    if (! istype(a, ans)) {
+	quit "dms_print called with non dms object";
+    }
+
+    /* print in dms form */
+    print a.deg : 'd' : a.min : 'm' : a.sec : 's' :;
+}
+
+
+define dms_abs(a)
+{
+    local obj dms ans;		/* temp object for dms type testing */
+    local deg;			/* return scalar value */
+
+    /* firewall - just absolute value non dms objects */
+    if (! istype(a, ans)) {
+    	return abs(a);
+    }
+
+    /* compute degrees */
+    deg = a.deg + a.min / 60 + a.sec / 3600;
+
+    /* return degrees */
+    return deg;
+}
+
+
+define dms_norm(a)
+{
+    local obj dms ans;		/* temp object for dms type testing */
+    local deg;			/* degrees */
+
+    /* firewall - arg must be a dms object */
+    if (! istype(a, ans)) {
+	quit "dms_norm called with non dms object";
+    }
+
+    /* square degrees (norm is the square of absolute value */
+    deg = dms_abs(a);
+
+    /* return degrees */
+    return deg*deg;
+}
+
+
+define dms_test(a)
+{
+    local obj dms ans;		/* temp value */
+
+    /* firewall - arg must be a dms object */
+    if (! istype(a, ans)) {
+	quit "dms_test called with non dms object";
+    }
+
+    /* return false of non-zero */
+    ans = fixdms(a);
+    if (ans.deg == 0 && ans.min == 0 && ans.sec == 0) {
+	/* false */
+	return 0;
+    }
+    /* true */
+    return 1;
+}
+
+
+define dms_int(a)
+{
+    local obj dms ans;		/* return value */
+
+    /* firewall - arg must be a dms object */
+    if (! istype(a, ans)) {
+	quit "dms_int called with non dms object";
+    }
+
+    /* normalize the argument */
+    ans = fixdms(a);
+
+    /* truncate to the nearest second */
+    ans.sec = int(ans.sec);
+
+    /* return value to the nearest second */
+    return ans;
+}
+
+
+define dms_frac(a)
+{
+    local obj dms ans;		/* return value */
+
+    /* firewall - arg must be a dms object */
+    if (! istype(a, ans)) {
+	quit "dms_frac called with non dms object";
+    }
+
+    /* normalize the argument */
+    ans = fixdms(a);
+
+    /* remove all but fractional seconds */
+    ans.deg = 0;
+    ans.min = 0;
+    ans.sec = frac(ans.sec);
+
+    /* return value to the second fraction */
+    return ans;
+}
+
+
+define dms_rel(a,b)
+{
+    local abs_a, abs_b;		/* scalars of the arguments */
+
+    /* compute scalars of the arguments */
+    abs_a = dms_abs(a);
+    abs_b = dms_abs(b);
+
+    /* return the comparison */
+    return cmp(abs_a, abs_b);
+}
+
+
+define dms_cmp(a,b)
+{
+    local abs_a, abs_b;		/* scalars of the arguments */
+
+    /* compute scalars of the arguments */
+    abs_a = dms_abs(a);
+    abs_b = dms_abs(b);
+
+    /* return the equality comparison */
+    return (abs_a == abs_b);
+}
+
+
+define dms_inc(a)
+{
+    local obj dms ans;		/* return value */
+
+    /* increment a dms object */
+    if (istype(a, ans)) {
+	ans = a;
+	++ans.sec;
+
+	/* return normalized result */
+	ans = fixdms(ans);
+	return ans;
+    }
+
+    /* increment a scalar */
+    return a+1;
+}
+
+
+define dms_dec(a)
+{
+    local obj dms ans;		/* return value */
+
+    /* decrement a dms object */
+    if (istype(a, ans)) {
+	ans = a;
+	--ans.sec;
+
+	/* return normalized result */
+	ans = fixdms(ans);
+	return ans;
+    }
+
+    /* decrement a scalar */
+    return a-1;
+}
+
+
+define fixdms(a)
+{
+    local obj dms ans;		/* temp value */
+
+    /* firewall */
+    if (! istype(a, ans)) {
+	quit "attempt to fix a non dms object";
+    }
+
+    /* force minutes to be intergral */
+    a.min += frac(a.deg) * 60;
+    a.deg = int(a.deg);
+
+    /* force degrees to be intergral */
+    a.sec += frac(a.min) * 60;
+    a.min = int(a.min);
+
+    /* carry excess seconds into minutes */
+    a.min += a.sec // 60;
+    a.sec %= 60;
+
+    /* carry excess minutes into degrees */
+    a.deg += a.min // 60;
+    a.min %= 60;
+
+    /* round degrees :-) */
+    a.deg %= 360;
+
+    /* return normalized result */
+    return a;
+}
+
+if (config("resource_debug") & 3) {
+    print "obj dms {deg, min, sec} defined";
+}

cal/dotest.cal

@@ -18,7 +18,7 @@
  *
  * @(#) $Revision: 30.2 $
  * @(#) $Id: dotest.cal,v 30.2 2007/03/16 11:09:54 chongo Exp $
- * @(#) $Source: /usr/local/src/cmd/calc/cal/RCS/dotest.cal,v $
+ * @(#) $Source: /usr/local/src/bin/calc/cal/RCS/dotest.cal,v $
  *
  * This file is not covered under version 2.1 of the GNU LGPL.
  *

cal/ellip.cal

@@ -19,7 +19,7 @@
  *
  * @(#) $Revision: 30.1 $
  * @(#) $Id: ellip.cal,v 30.1 2007/03/16 11:09:54 chongo Exp $
- * @(#) $Source: /usr/local/src/cmd/calc/cal/RCS/ellip.cal,v $
+ * @(#) $Source: /usr/local/src/bin/calc/cal/RCS/ellip.cal,v $
  *
  * Under source code control:	1990/02/15 01:50:33
  * File existed as early as:	before 1990

cal/factorial.cal

@@ -0,0 +1,204 @@
+/*
+ * factorial - implementation of different algorithms for the factorial
+ *
+ * Copyright (C) 2013 Christoph Zurnieden
+ *
+ * Calc is open software; you can redistribute it and/or modify it under
+ * the terms of the version 2.1 of the GNU Lesser General Public License
+ * as published by the Free Software Foundation.
+ *
+ * Calc is distributed in the hope that it will be useful, but WITHOUT
+ * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
+ * or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General
+ * Public License for more details.
+ *
+ * A copy of version 2.1 of the GNU Lesser General Public License is
+ * distributed with calc under the filename COPYING-LGPL.  You should have
+ * received a copy with calc; if not, write to Free Software Foundation, Inc.
+ * 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
+ *
+ * @(#) $Revision: 30.4 $
+ * @(#) $Id: factorial.cal,v 30.4 2013/08/18 20:01:53 chongo Exp $
+ * @(#) $Source: /usr/local/src/bin/calc/cal/RCS/factorial.cal,v $
+ *
+ * Under source code control:	2013/08/11 01:31:28
+ * File existed as early as:	2013
+ */
+
+
+/*
+ * hide internal function from resource debugging
+ */
+static resource_debug_level;
+resource_debug_level = config("resource_debug", 0);
+
+
+/*
+  get dependencies
+*/
+read -once toomcook;
+
+
+/* A simple list to keep things...uhm...simple?*/
+static __CZ__primelist = list();
+
+/* Helper for primorial: fill list with primes in range a,b */
+define __CZ__fill_prime_list(a,b)
+{
+  local k;
+  k=a;
+  if(isprime(k))k--;
+  while(1){
+    k = nextprime(k);
+    if(k > b) break;
+    append(__CZ__primelist,k );
+  }
+}
+
+/* Helper for factorial: how often prime p divides the factorial of n */
+define __CZ__prime_divisors(n,p)
+{
+    local q,m;
+    q = n;
+    m = 0;
+    if (p > n) return 0;
+    if (p > n/2) return 1;
+    while (q >= p) {
+        q  = q//p;
+        m += q;
+    }
+    return m;
+}
+
+/*
+  Wrapper. Please set cut-offs to own taste and hardware.
+*/
+define factorial(n){
+  local prime result shift prime_list k k1 k2 expo_list pix cut primorial;
+
+  result = 1;
+  prime  = 2;
+
+  if(!isint(n))  {
+    return newerror("factorial(n): n is not an integer"); ## or gamma(n)?
+  }
+  if(n < 0)      return newerror("factorial(n): n < 0");
+  if(n < 9000 && !isdefined("test8900"))   {
+    ## builtin is implemented with splitting but only with
+    ## Toom-Cook 2 (by Karatsuba (the father))
+    return n!;
+  }
+
+  shift = __CZ__prime_divisors(n,prime);
+  prime = 3;
+  cut = n//2;
+  pix   = pix(cut);
+  prime_list = mat[pix];
+  expo_list  = mat[pix];
+
+  k = 0;
+/*
+  Peter Borwein's algorithm
+
+  @Article{journals/jal/Borwein85,
+  author = {Borwein, Peter B.},
+  title  = {On the Complexity of Calculating Factorials.},
+  journal = {J. Algorithms},
+  year   = {1985},
+  number = {3},
+  url    = {http://dblp.uni-trier.de/db/journals/jal/jal6.html#Borwein85}
+*/
+
+  do {
+    prime_list[k]  = prime;
+    expo_list[k++] = __CZ__prime_divisors(n,prime);
+    prime          = nextprime(prime);
+  }while(prime <= cut);
+
+  /* size of the largest exponent in bits  */
+  k1 = highbit(expo_list[0]);
+  k2 = size(prime_list)-1;
+
+  for(;k1>=0;k1--){
+    /*
+       the cut-off for T-C-4 ist still to low, using T-C-3 here
+       TODO: check cutoffs
+     */
+    result = toomcook3square(result);
+    /*
+      almost all time is spend in this loop, so cutting of the
+      upper half of the primes makes sense
+    */
+    for(k=0; k<=k2; k++) {
+      if((expo_list[k] & (1 << k1)) != 0) {
+        result *= prime_list[k];
+      }
+    }
+
+  }
+  primorial = primorial( cut, n);
+  result *= primorial;
+  result <<= shift;
+  return result;
+}
+
+/*
+  Helper for primorial: do the product with binary splitting
+  TODO: do it without the intermediate list
+*/
+define __CZ__primorial__lowlevel( a, b ,p)
+{
+  local  c;
+  if( b == a) return p ;
+  if( b-a > 1){
+    c= (b + a) >> 1;
+    return  __CZ__primorial__lowlevel( a   , c , __CZ__primelist[a] )
+         *  __CZ__primorial__lowlevel( c+1 , b , __CZ__primelist[b] ) ;
+  }
+  return  __CZ__primelist[a] * __CZ__primelist[b];
+}
+
+/*
+   Primorial, Product of consecutive primes in range a,b
+
+  Originally meant to do primorials with a start different from 2, but
+  found out that this is faster at about a=1,b>=10^5 than the builtin
+  function pfact().  With the moderately small list a=1,b=10^6 (78498
+  primes) it is 3 times faster.  A quick look-up showed what was already
+  guessed: pfact() does it linearly.  (BTW: what is the time complexity
+  of the primorial with the naive algorithm?)
+*/
+define primorial(a,b)
+{
+  local C1 C2;
+  if(!isint(a)) return newerror("primorial(a,b): a is not an integer");
+  else if(!isint(b)) return newerror("primorial(a,b): b is not an integer");
+  else if(a < 0) return newerror("primorial(a,b): a < 0");
+  else if( b < 2 ) return newerror("primorial(a,b): b < 2");
+  else if( b < a) return newerror("primorial(a,b): b < a");
+  else{
+    /* last prime < 2^32 is also  max. prime for nextprime()*/
+    if(b >= 4294967291) return newerror("primorial(a,b): max. prime exceeded");
+    if(b ==  2) return 2;
+    /*
+      Can be extended by way of pfact(b)/pfact(floor(a-1/2)) for small a
+     */
+    if(a<=2 && b < 10^5) return pfact(b);
+    /* TODO: use pix() and a simple array (mat[])instead*/
+     __CZ__primelist = list();
+     __CZ__fill_prime_list(a,b);
+    C1 = size(__CZ__primelist)-1;
+    return __CZ__primorial__lowlevel( 0, C1,1)
+  }
+}
+
+
+/*
+ * restore internal function from resource debugging
+ * report important interface functions
+ */
+config("resource_debug", resource_debug_level),;
+if (config("resource_debug") & 3) {
+    print "factorial(n)";
+    print "primorial(a, b)";
+}

cal/factorial2.cal

@@ -0,0 +1,723 @@
+/*
+ * factorial2 - implementation of different factorial related functions
+ *
+ * Copyright (C) 2013 Christoph Zurnieden
+ *
+ * Calc is open software; you can redistribute it and/or modify it under
+ * the terms of the version 2.1 of the GNU Lesser General Public License
+ * as published by the Free Software Foundation.
+ *
+ * Calc is distributed in the hope that it will be useful, but WITHOUT
+ * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
+ * or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General
+ * Public License for more details.
+ *
+ * A copy of version 2.1 of the GNU Lesser General Public License is
+ * distributed with calc under the filename COPYING-LGPL.  You should have
+ * received a copy with calc; if not, write to Free Software Foundation, Inc.
+ * 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
+ *
+ * @(#) $Revision: 30.4 $
+ * @(#) $Id: factorial2.cal,v 30.4 2013/08/18 20:01:53 chongo Exp $
+ * @(#) $Source: /usr/local/src/bin/calc/cal/RCS/factorial2.cal,v $
+ *
+ * Under source code control:	2013/08/11 01:31:28
+ * File existed as early as:	2013
+ */
+
+
+/*
+ * hide internal function from resource debugging
+ */
+static resource_debug_level;
+resource_debug_level = config("resource_debug", 0);
+
+
+/*
+  get dependencies
+*/
+read -once factorial toomcook specialfunctions;
+
+
+/*
+  Factorize a factorial and put the result in a 2-column matrix with pi(n) rows
+  mat[ primes , exponent ]
+  Result can be restricted to start at a prime different from 2 with the second
+  argument "start". That arguments gets taken at face value if it prime and
+  smaller than n, otherwise the next larger prime is taken if that prime is
+  smaller than n.
+*/
+
+define __CZ__factor_factorial(n,start){
+  local prime prime_list k pix stop;
+
+
+  if(!isint(n)) return
+    newerror("__CZ__factor_factorial(n,start): n is not integer");
+  if(n <  0)    return newerror("__CZ__factor_factorial(n,start): n < 0");
+  if(n == 1)    return newerror("__CZ__factor_factorial(n,start): n == 1");
+
+  if(start){
+    if(!isint(start) && start < 0 && start > n)
+      return newerror("__CZ__factor_factorial(n,start): value of "
+      		      "parameter 'start' out of range");
+    if(start == n && isprime(n)){
+      prime_list = mat[1 , 2];
+      prime_list[0,0] = n;
+      prime_list[0,1] = 1;
+    }
+    else if(!isprime(start) && nextprime(start) >n)
+      return newerror("__CZ__factor_factorial(n,start): value of parameter "
+      		      "'start' out of range");
+    else{
+      if(!isprime(start))  prime = nextprime(start);
+      else prime = start;
+    }
+  }
+  else
+    prime = 2;
+
+  pix   = pix(n);
+  if(start){
+    pix -= pix(prime) -1;
+  }
+  prime_list = mat[pix , 2];
+
+  k = 0;
+
+  do {
+    prime_list[k  ,0] = prime;
+    prime_list[k++,1] = __CZ__prime_divisors(n,prime);
+    prime          = nextprime(prime);
+  }while(prime <= n);
+
+  return prime_list;
+}
+
+/*
+
+  subtracts exponents of n_1! from exponents of n_2! with n_1<=n_2
+
+  Does not check for size or consecutiveness of the primes or a carry
+*/
+
+define __CZ__subtract_factored_factorials(matrix_2n,matrix_n){
+  local k ret len1,len2,tmp count p e;
+  len1 = size(matrix_n)/2;
+  len2 = size(matrix_2n)/2;
+  if(len2<len1){
+
+    swap(len1,len2);
+    tmp = matrix_n;
+    matrix_n = matrix_2n;
+    matrix_2n = tmp;
+  }
+  tmp = mat[len1,2];
+  k   = 0;
+
+  for(;k<len1;k++){
+     p = matrix_2n[k,0];
+     e = matrix_2n[k,1] - matrix_n[k,1];
+     if(e!=0){
+       tmp[count  ,0] = p;
+       tmp[count++,1] = e;
+     }
+  }
+  ret = mat[count + (len2-len1),2];
+  for(k=0;k<count;k++){
+     ret[k,0] = tmp[k,0];
+     ret[k,1] = tmp[k,1];
+  }
+
+  free(tmp);
+  for(k=len1;k<len2;k++){
+     ret[count,0] = matrix_2n[k,0];
+     ret[count++,1] = matrix_2n[k,1];
+  }
+  return ret;
+}
+
+/*
+
+  adds exponents of n_1! to exponents of n_2! with n_1<=n_2
+
+  Does not check for size or consecutiveness of the primes or a carry
+*/
+
+define __CZ__add_factored_factorials(matrix_2n,matrix_n){
+  local k ret len1,len2,tmp;
+  len1 = size(matrix_n)/2;
+  len2 = size(matrix_2n)/2;
+  if(len2<len1){
+    swap(len1,len2);
+    tmp = matrix_n;
+    matrix_n = matrix_2n;
+    matrix_2n = tmp;
+  }
+  ret = mat[len2,2];
+  k   = 0;
+  for(;k<len1;k++){
+     ret[k,0] = matrix_2n[k,0];
+     ret[k,1] = matrix_2n[k,1] + matrix_n[k,1];
+  }
+  for(;k<len2;k++){
+     ret[k,0] = matrix_2n[k,0];
+     ret[k,1] = matrix_2n[k,1];
+  }
+  return ret;
+}
+
+/*
+   Does not check if all exponents are positive
+
+
+                   timings
+                   this		  comb	    comb-this     rel. k/n
+;  benchmark_binomial(10,13)
+n=2^13 k=2^10 	 0.064004 	 0.016001	+ 	0.76923076923076923077
+n=2^13 k=2^11 	 0.064004 	 0.048003	+ 	0.84615384615384615385
+n=2^13 k=2^12 	 0.068004 	 0.124008	- 	0.92307692307692307692
+;  benchmark_binomial(10,15)
+n=2^15 k=2^10 	 0.216014 	 0.024001	+ 	0.66666666666666666667
+n=2^15 k=2^11 	 0.220014 	 0.064004	+ 	0.73333333333333333333
+n=2^15 k=2^12 	 0.228014 	 0.212014	+ 	0.8
+n=2^15 k=2^13 	 0.216013 	 0.664042	- 	0.86666666666666666667
+n=2^15 k=2^14 	 0.240015 	 1.868117	- 	0.93333333333333333333
+;  benchmark_binomial(11,15)
+n=2^15 k=2^11 	 0.216014 	 0.068004	+ 	0.73333333333333333333
+n=2^15 k=2^12 	 0.236015 	 0.212013	+ 	0.8
+n=2^15 k=2^13 	 0.216013 	 0.656041	- 	0.86666666666666666667
+n=2^15 k=2^14 	 0.244016 	 1.872117	- 	0.93333333333333333333
+;  benchmark_binomial(11,18)
+n=2^18 k=2^11 	 1.652103 	 0.100006	+ 	0.61111111111111111111
+n=2^18 k=2^12 	 1.608101 	 0.336021	+ 	0.66666666666666666667
+n=2^18 k=2^13 	 1.700106 	 1.140071	+ 	0.72222222222222222222
+n=2^18 k=2^14 	 1.756109 	 3.924245	- 	0.77777777777777777778
+n=2^18 k=2^15 	 2.036127 	 13.156822	- 	0.83333333333333333333
+n=2^18 k=2^16 	 2.172135 	 41.974624	- 	0.88888888888888888889
+n=2^18 k=2^17 	 2.528158 	 121.523594	- 	0.94444444444444444444
+;  benchmark_binomial(15,25)
+n=2^25 k=2^15 	 303.790985 	 38.266392	+ 	0.6
+;  benchmark_binomial(17,25)
+n=2^25 k=2^17 	 319.127944 	 529.025062	- 	0.68
+*/
+
+define benchmark_binomial(s,limit){
+  local ret k A B T1 T2 start end N K;
+    N = 2^(limit);
+  for(k=s;k<limit;k++){
+    K = 2^k;
+   start=usertime();A=binomial(N,K);end=usertime();
+   T1 = end-start;
+   start=usertime();B=comb(N,K);end=usertime();
+   T2 = end-start;
+   print "n=2^"limit,"k=2^"k,"	",T1,"	",T2,T1<T2?"-":"+","	"k/limit;
+   if(A!=B){
+     print "false";
+    break;
+   }
+  }
+}
+
+define __CZ__multiply_factored_factorial(matrix,stop){
+  local prime result shift prime_list k k1 k2 expo_list pix count start;
+  local hb flag;
+
+  result = 1;
+  shift  = 0;
+
+
+  if(!ismat(matrix))
+   return newerror("__CZ__multiply_factored_factorial(matrix): "
+   		   "argument matrix not a matrix ");
+  if(!matrix[0,0])
+    return
+      newerror("__CZ__multiply_factored_factorial(matrix): "
+      	       "matrix[0,0] is null/0");
+
+  if(!isnull(stop))
+    pix = stop;
+  else
+    pix   = size(matrix)/2-1;
+
+  if(matrix[0,0] == 2 && matrix[0,1] > 0){
+    shift = matrix[0,1];
+    if(pix-1 == 0)
+      return 2^matrix[0,1];
+  }
+
+  /*
+    This is a more general way to do the multiplication, so any optimization
+    must have been done by the caller.
+  */
+  k = 0;
+  /*
+     The size of the largest exponent in bits is calculated dynamically.
+     Can be done more elegantly and saves one run over the whole array if done
+     inside the main loop.
+   */
+  hb =0;
+  for(k=0;k<pix;k++){
+    k1=highbit(matrix[k,1]);
+    if(hb < k1)hb=k1;
+  }
+
+  k2 = pix;
+  start = 0;
+  if(shift) start++;
+
+  for(k1=hb;k1>=0;k1--){
+    /*
+       the cut-off for T-C-4 ist still too low, using T-C-3 here
+       TODO: check cutoffs
+     */
+    result = toomcook3square(result);
+
+    for(k=start; k<=k2; k++) {
+      if((matrix[k,1] & (1 << k1)) != 0) {
+        result *= matrix[k,0];
+      }
+    }
+  }
+
+  result <<= shift;
+  return result;
+}
+
+/*
+    Compute binomial coeficients n!/(k!(n-k)!)
+
+    One of the rare cases where a formula once meant to ease manual computation
+    is actually the (aymptotically) fastest way to do it (in July 2013) for
+    the extreme case binomial(2N,N) but for a high price, the memory
+    needed is pi(N)--theoretically.
+*/
+define binomial(n,k){
+  local ret factored_n factored_k factored_nk denom num quot K prime_list prime;
+  local pix diff;
+
+  if(!isint(n) || !isint(k))
+    return newerror("binomial(n,k): input is not integer");
+  if(n<0 || k<0)
+    return newerror("binomial(n,k): input is not >= 0"); ;
+  if(n<k )   return 0;
+  if(n==k)   return 1;
+  if(k==0)   return 1;
+  if(k==1)   return n;
+  if(n-k==1) return n;
+  /*
+      cut-off depends on real size of n,k and size of n/k
+      The current cut-off is to small for large n, e.g.:
+      for 2n=2^23, k=n-n/2 the quotient is q=2n/k=0.25. Empirical tests showed
+      that 2n=2^23 and k=2^16  with q=0.0078125 are still faster than the
+      builtin function.
+
+      The symmetry (n,k) = (n,n-k) is of not much advantage here. One way
+      might be to get closer to k=n/2 if k<n-k but only if the difference
+      is small and n very large.
+   */
+  if(n<2e4 && !isdefined("test8900"))  return comb(n,k);
+  if(n<2e4 && k< n-n/2 && !isdefined("test8900")) return comb(n,k);
+  /*
+    This should be done in parallel to save some memory, e.g. no temporary
+    arrays are needed, all can be done inline.
+    The theoretical memory needed is pi(k).
+    Which is still a lot.
+  */
+
+  prime = 2;
+  pix   = pix(n);
+  prime_list = mat[pix , 2];
+  K = 0;
+  do {
+    prime_list[K  ,0] = prime;
+    diff = __CZ__prime_divisors(n,prime)-
+	   ( __CZ__prime_divisors(n-k,prime)+__CZ__prime_divisors(k,prime));
+    if(diff != 0)
+      prime_list[K++,1] = diff;
+    prime          = nextprime(prime);
+  }while(prime <= k);
+
+  do {
+    prime_list[K  ,0] = prime;
+    diff = __CZ__prime_divisors(n,prime)-__CZ__prime_divisors(n-k,prime);
+    if(diff != 0)
+      prime_list[K++,1] = diff;
+    prime          = nextprime(prime);
+  }while(prime <= n-k);
+
+  do {
+    prime_list[K  ,0] = prime;
+    prime_list[K++,1] = __CZ__prime_divisors(n,prime);
+    prime          = nextprime(prime);
+  }while(prime <= n);
+  ##print K,pix(k),pix(n-k),pix(n);
+  ##factored_k  = __CZ__factor_factorial(k,1);
+  ##factored_nk = __CZ__factor_factorial(n-k,1);
+
+  ##denom = __CZ__add_factored_factorials(factored_k,factored_nk);
+    ##free(factored_k,factored_nk);
+  ##num  = __CZ__factor_factorial(n,1);
+  ##quot  = __CZ__subtract_factored_factorials( num , denom );
+    ##free(num,denom);
+
+  ret = __CZ__multiply_factored_factorial(`prime_list,K-1);
+
+  return ret;
+}
+
+/*
+    Compute large catalan numbers  C(n) = binomial(2n,n)/(n+1) with
+    cut-off: (n>5e4)
+    Needs a lot of memory.
+*/
+define bigcatalan(n){
+  if(!isint(n) )return newerror("bigcatalan(n): n is not integer");
+  if( n<0) return newerror("bigcatalan(n): n < 0");
+  if( n<5e4 && !isdefined("test8900") ) return catalan(n);
+  return binomial(2*n,n)/(n+1);
+}
+
+/*
+  df(-111) = -1/3472059605858239446587523014902616804783337112829102414124928
+  		7753332469144201839599609375
+
+ df(-3+1i) = 0.12532538977287649201-0.0502372106177184607i
+ df(2n + 1) = (2*n)!/(n!*2^n)
+*/
+define __CZ__double_factorial(n){
+  local n1 n2 diff prime pix K prime_list k;
+  prime = 3;
+  pix   = pix(2*n)+1;
+  prime_list = mat[pix , 2];
+  K = 0;
+  do {
+    prime_list[K  ,0] = prime;
+    diff = __CZ__prime_divisors(2*n,prime)-( __CZ__prime_divisors(n,prime));
+    if(diff != 0)
+      prime_list[K++,1] = diff;
+    prime          = nextprime(prime);
+  }while(prime <= n);
+  do {
+    prime_list[K  ,0] = prime;
+    prime_list[K++,1] = __CZ__prime_divisors(2*n,prime);
+    prime          = nextprime(prime);
+  }while(prime <= 2*n);
+  return __CZ__multiply_factored_factorial(prime_list,K);
+/*
+        n1=__CZ__factor_factorial(2*n,1);
+        n1[0,1] = n1[0,1]-n;
+        n2=__CZ__factor_factorial(n,1);
+        diff=__CZ__subtract_factored_factorials( n1 , n2 );
+        return __CZ__multiply_factored_factorial(diff);
+*/
+}
+
+##1, 1, 3, 15, 105, 945, 10395, 135135, 2027025, 34459425, 654729075,
+##13749310575, 316234143225, 7905853580625, 213458046676875,
+##6190283353629375, 191898783962510625, 6332659870762850625,
+##221643095476699771875, 8200794532637891559375
+
+## 1, 2, 8, 48, 384, 3840, 46080, 645120, 10321920, 185794560,
+##3715891200, 81749606400, 1961990553600, 51011754393600,
+##1428329123020800, 42849873690624000, 1371195958099968000,
+##46620662575398912000, 1678343852714360832000, 63777066403145711616000
+define doublefactorial(n){
+    local n1 n2 diff eps ret;
+    if(!isint(n) ){
+       /*
+           Probably one of the not-so-good ideas. See result of
+            http://www.wolframalpha.com/input/?i=doublefactorial%28a%2Bbi%29
+        */
+       eps=epsilon(epsilon()*1e-2);
+       ret =  2^(n/2-1/4 * cos(pi()* n)+1/4) * pi()^(1/4 *
+       	      cos(pi()* n)-1/4)* gamma(n/2+1);
+       epsilon(eps);
+       return ret;
+    }
+    if(n==2) return 2;
+    if(n==3) return 3;
+    switch(n){
+      case -1:
+      case 0 : return 1;break;
+      case 2 : return 2;break;
+      case 3 : return 3;break;
+      case 4 : return 8;break;
+      default: break;
+    }
+    if(isodd(n)){
+      /*
+        TODO: find reasonable cutoff
+        df(2n + 1) = (2*n)!/(n!*2^n)
+      */
+      if(n>0){
+         n = (n+1)//2;
+         return __CZ__double_factorial(n);
+      }
+      else{
+        if(n == -3 ) return -1;
+        n = ((-n)-1)/2;
+        return ((-1)^-n)/__CZ__double_factorial(n);
+       }
+    }
+    else{
+      /*
+         I'm undecided here. The formula for complex n is valid for the negative
+         integers, too.
+      */
+      n = n>>1;
+      if(n>0){
+        if(!isdefined("test8900"))
+          return factorial(n)<<n;
+        else
+          return n!<<n;
+      }
+      else
+        return newerror("doublefactorial(n): even(n) < 0");
+   }
+}
+
+/*
+    Algorithm 3.17,
+    Donald Kreher and Douglas Simpson,
+    Combinatorial Algorithms,
+    CRC Press, 1998, page 89.
+*/
+static __CZ__stirling1;
+static __CZ__stirling1_n = -1;
+static __CZ__stirling1_m = -1;
+
+define stirling1(n,m){
+  local i j k;
+  if(n<0)return newerror("stirling1(n,m): n <= 0");
+  if(m<0)return newerror("stirling1(n,m): m < 0");
+  if(n<m) return 0;
+  if(n==m) return 1;
+  if(m==0 || n==0) return 0;
+  /* We always use the list */
+  /*
+  if(m=1){
+    if(iseven(n)) return -factorial(n-1);
+    else return factorial(n-1);
+  }
+  if(m == n-1){
+    if(iseven(n)) return -binomial(n,2);
+    else return -binomial(n,2);
+  }
+  */
+  if(__CZ__stirling1_n >= n && __CZ__stirling1_m  >= m){
+    return __CZ__stirling1[n,m];
+  }
+  else{
+    __CZ__stirling1      = mat[n+1,m+1];
+    __CZ__stirling1[0,0] = 1;
+    for(i=1;i<=n;i++)
+      __CZ__stirling1[i,0] = 0;
+    for(i=1;i<=n;i++){
+      for(j=1;j<=m;j++){
+        if(j<=i){
+          __CZ__stirling1[i, j] =   __CZ__stirling1[i - 1, j - 1] - (i - 1)\
+                                  * __CZ__stirling1[i - 1, j];
+        }
+        else{
+         __CZ__stirling1[i, j] = 0;
+        }
+      }
+    }
+    __CZ__stirling1_n = n;
+    __CZ__stirling1_m = m;
+    return __CZ__stirling1[n,m];
+  }
+}
+
+define stirling2(n,m){
+  local k sum;
+  if(n<0)return newerror("stirling2(n,m): n < 0");
+  if(m<0)return newerror("stirling2(n,m): m < 0");
+  if(n<m) return 0;
+  if(n==0 && n!=m) return 0;
+  if(n==m) return 1;
+  if(m==0 )return 0;
+  if(m==1) return 1;
+  if(m==2) return 2^(n-1)-1;
+  /*
+    There are different methods to speed up alternating sums.
+    This one doesn't.
+   */
+  if(isdefined("test8900")){
+    for(k=0;k<=m;k++){
+      sum += (-1)^(m-k)*comb(m,k)*k^n;
+    }
+  return sum/(m!);
+  }
+  else{
+    for(k=0;k<=m;k++){
+      sum += (-1)^(m-k)*binomial(m,k)*k^n;
+    }
+  return sum/factorial(m);
+  }
+}
+
+static __CZ__stirling2;
+static __CZ__stirling2_n = -1;
+static __CZ__stirling2_m = -1;
+define stirling2caching(n,m){
+  local nm i j ;
+  if(n<0)return newerror("stirling2iter(n,m): n < 0");
+  if(m<0)return newerror("stirling2iter(n,m): m < 0");
+  /* no shortcuts here */
+
+  if(n<m) return 0;
+  if(n==0 && n!=m) return 0;
+  if(n==m) return 1;
+  if(m==0 )return 0;
+  if(m==1) return 1;
+  if(m==2) return 2^(n-1)-1;
+
+  nm = n-m;
+  if(__CZ__stirling2_n >= n && __CZ__stirling2_m >= m){
+    return __CZ__stirling2[n,m];
+  }
+  else{
+    __CZ__stirling2 = mat[n+1,m+1];
+    __CZ__stirling2[0,0] = 1;
+    for(i=1;i<=n;i++){
+      __CZ__stirling2[i,0] = 0;
+      for(j=1;j<=m;j++){
+        if(j<=i){
+          __CZ__stirling2[i, j] =   __CZ__stirling2[i -1, j -1] + (j )\
+                                  * __CZ__stirling2[i - 1, j];
+        }
+        else{
+         __CZ__stirling2[i, j] = 0;
+        }
+      }
+    }
+  }
+  __CZ__stirling2_n = (n);
+  __CZ__stirling2_m = (m);
+  return __CZ__stirling2[n,m];
+}
+
+define bell(n){
+  local sum s2list k A;
+
+  if(!isint(n)) return newerror("bell(n): n is not integer");
+  if(n < 0) return newerror("bell(n): n is not positive");
+  /* place some more shortcuts here?*/
+  if(n<=15){
+    mat A[16] = {
+                  1, 1, 2, 5, 15, 52, 203, 877, 4140, 21147, 115975, 678570,
+                  4213597, 27644437, 190899322, 1382958545
+                };
+    return A[n];
+  }
+  /* Start by generating the list of stirling numbers of the second kind */
+  s2list = stirling2caching(n,n//2);
+  if(iserror(s2list))
+    return newerror("bell(n): could not build stirling num. list");
+  sum = 0;
+  for(k=1;k<=n;k++){
+      sum += stirling2caching(n,k);
+  }
+  return sum;
+}
+
+define subfactorialrecursive(n){
+  if(n==0) return 1;
+  if(n==1) return 0;
+  if(n==2) return 1;
+  return n * subfactorialrecursive(n-1) + (-1)^n;
+}
+
+/* This is, quite amusingely, faster than the very same algorithm in
+   PARI/GP + GMP*/
+define subfactorialiterative(n){
+  local k temp1 temp2 ret;
+  if(n==0) return 1;
+  if(n==1) return 0;
+  if(n==2) return 1;
+  temp1 = 0;
+  ret   = 1;
+  for(k=3;k<=n;k++){
+    temp2 = temp1;
+    temp1 = ret;
+    ret =  (k-1) *(temp1 + temp2);
+  }
+  return ret;
+}
+
+define subfactorial(n){
+  local epsilon eps ret lnfact;
+  if(!isint(n))return  newerror("subfactorial(n): n is not integer.");
+  if(n < 0)return  newerror("subfactorial(n): n < 0");
+  return subfactorialiterative(n);
+}
+
+define risingfactorial(x,n){
+  local num denom quot ret;
+  if(n == 1) return x;
+  if(x==0) return newerror("risingfactorial(x,n): x == 0");
+  if(!isint(x) || !isint(n)){
+    return gamma(x+n)/gamma(x);
+  }
+  if(x<1)return newerror("risingfactorial(x,n): integer x and x < 1");
+  if(x+n < 1)return newerror("risingfactorial(x,n): integer x+n and x+n < 1");
+  if(x<9000&&n<9000){
+    return (x+n-1)!/(x-1)!;
+  }
+  else{
+    num   = __CZ__factor_factorial(x+n-1,1);
+    denom = __CZ__factor_factorial(x-1,1);
+    quot  = __CZ__subtract_factored_factorials( num , denom );
+      free(num,denom);
+    ret = __CZ__multiply_factored_factorial(quot);
+    return ret;
+  }
+}
+
+define fallingfactorial(x,n){
+  local num denom quot ret;
+  if(n == 0) return 1;
+
+  if(!isint(x) || !isint(n)){
+    if(x == n) return gamma(x+1);
+    return gamma(x+1)/gamma(x-n+1);
+  }
+  else{
+    if(x<0 || x-n < 0)
+     return newerror("fallingfactorial(x,n): integer x<0 or x-n < 0");
+    if(x == n) return factorial(x);
+    if(x<9000&&n<9000){
+      return (x)!/(x-n)!;
+    }
+    else{
+      num   = __CZ__factor_factorial(x,1);
+      denom = __CZ__factor_factorial(x-n,1);
+      quot  = __CZ__subtract_factored_factorials( num , denom );
+        free(num,denom);
+      ret = __CZ__multiply_factored_factorial(quot);
+      return ret;
+    }
+  }
+}
+
+
+/*
+ * restore internal function from resource debugging
+ * report important interface functions
+ */
+config("resource_debug", resource_debug_level),;
+if (config("resource_debug") & 3) {
+    print "binomial(n,k)";
+    print "bigcatalan(n)";
+    print "doublefactorial(n)";
+    print "subfactorial(n)";
+    print "stirling1(n,m)";
+    print "stirling2(n,m)";
+    print "stirling2caching(n,m)";
+    print "bell(n)";
+    print "subfactorial(n)";
+    print "risingfactorial(x,n)";
+    print "fallingfactorial(x,n)";
+}

cal/gvec.cal

@@ -0,0 +1,107 @@
+/*
+ * gvec - vectorize any single-input function or trailing operator
+ *
+ * This version accepts arbitrary number of arguments, but of course
+ * they must all be same length vectors.
+ *
+ * The gvec function is for use in either a two-arg function or a two-arg
+ * operation "function" must be first; calc doesn't care how many more
+ * arguments there actually are.
+ *
+ * @(#) $Revision: 30.3 $
+ * @(#) $Id: gvec.cal,v 30.3 2011/05/23 23:00:55 chongo Exp $
+ * @(#) $Source: /usr/local/src/bin/calc/cal/RCS/gvec.cal,v $
+ *
+ * Under source code control:   2011/03/31 17:54:55
+ * File existed as early as:    2010
+ *
+ * By Carl Witthoft carl at witthoft dot com
+ */
+
+define gvec(function, vector)
+{
+    local xlen,y,foo;
+    local precx = 1e-50;	/* default for now */
+    local argc = param(0)-1;
+    local old_tilde;		/* previous config("tilde") */
+
+    /*
+     * parse args
+     */
+    local plist = mat[argc];
+    if (config("resource_debug") & 8) {
+	print "plist=", plist;
+	print "argc=", argc;
+    }
+    for(local i = 0; i< argc; i++) {
+	local ii = i + 2;
+	if (config("resource_debug") & 8) {
+	    print "ii=", ii;
+	    print "param(" : ii : "}=", param(ii);
+	    print "size(param(" : ii : ")=", size(param(ii));
+	}
+	plist[i] = size(param(ii));
+    }
+    local slist=sort(plist);
+    if (config("resource_debug") & 8) {
+	print "plist=", plist;
+    }
+    local argm = argc-1;
+    if (config("resource_debug") & 8) {
+	print "argm=", argm;
+    }
+    if (slist[0] != slist[argm]) {
+	quit "lengths don't match";
+    }
+    xlen = size(vector);
+    y = mat[xlen];
+
+    /*
+     * We can't do str(vector[j]) outside loop, eval() petulantly refuses to
+     * look at local variables.
+     *
+     * Also we need to config("tilde",0) to turn off lead tilde
+     * (so str(vector[j]) looks like a number.
+     */
+    old_tilde = config("tilde",0);
+
+    /*
+     * Ok, now check to see if "function" is a function.  If not, it's an
+     * operation and it's up to user to make it valid
+     */
+    if (isdefined(function)) {
+
+	/* yep, it's a function, either builtin or user-defined */
+	for (local j=0; j<xlen; j++) {
+
+	    /* build the function call */
+	    foo = strcat(function, "(");
+	    for (local jj = 0; jj<argc; jj++) {
+		foo = strcat(foo , str(param(jj+2)[j]), ",");
+	    }
+	    foo = strcat(foo, str(precx), ")");
+	    if (config("resource_debug") & 8) {
+		print "foo=", foo;
+	    }
+	    y[j] = eval(foo);
+	}
+
+    /*
+     * it is an operator --  multi-argument operator makes no sense
+     */
+    } else {
+	if (argc > 1) {
+	    quit "Error: operator can accept only one argument";
+	}
+	for (j=0; j<xlen; j++) {
+	    foo = strcat(str(vector[j]), function);
+	    y[j] = eval(foo);
+	}
+    }
+
+    /* restore tilde mode if needed */
+    config("tilde", old_tilde);
+
+    /* return result */
+    return y;
+}

cal/hello.cal

@@ -19,7 +19,7 @@
  *
  * @(#) $Revision: 30.1 $
  * @(#) $Id: hello.cal,v 30.1 2007/03/16 11:09:54 chongo Exp $
- * @(#) $Source: /usr/local/src/cmd/calc/cal/RCS/hello.cal,v $
+ * @(#) $Source: /usr/local/src/bin/calc/cal/RCS/hello.cal,v $
  *
  * Under source code control:	1996/11/13 13:25:43
  * File existed as early as:	1996

cal/hms.cal

@@ -0,0 +1,384 @@
+/*
+ * hms - calculate in hours, minutes, and seconds
+ *
+ * Copyright (C) 2010  Landon Curt Noll
+ *
+ * Calc is open software; you can redistribute it and/or modify it under
+ * the terms of the version 2.1 of the GNU Lesser General Public License
+ * as published by the Free Software Foundation.
+ *
+ * Calc is distributed in the hope that it will be useful, but WITHOUT
+ * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
+ * or FITNESS FOR A PARTICULAR PURPOSE.	 See the GNU Lesser General
+ * Public License for more details.
+ *
+ * A copy of version 2.1 of the GNU Lesser General Public License is
+ * distributed with calc under the filename COPYING-LGPL.  You should have
+ * received a copy with calc; if not, write to Free Software Foundation, Inc.
+ * 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
+ *
+ * @(#) $Revision: 30.2 $
+ * @(#) $Id: hms.cal,v 30.2 2010/09/02 06:14:16 chongo Exp $
+ * @(#) $Source: /usr/local/src/bin/calc/cal/RCS/hms.cal,v $
+ *
+ * Under source code control:	2010/09/01 17:14:55
+ * File existed as early as:	2010
+ *
+ * Share and enjoy!  :-)	http://www.isthe.com/chongo/tech/comp/calc/
+ */
+
+
+obj hms {hour, min, sec};
+
+define hms(hour, min, sec)
+{
+    local obj hms ans;		/* return value */
+
+    /* default missing args to 0 */
+    if (isnull(sec)) {
+	sec = 0;
+    }
+    if (isnull(min)) {
+	min = 0;
+    }
+
+    /* load object */
+    ans.hour = hour;
+    ans.min = min;
+    ans.sec = sec;
+
+    /* return properly formed object */
+    ans = fixhms(ans);
+    return ans;
+}
+
+
+define hms_add(a, b)
+{
+    local obj hms ans;		/* return value */
+
+    /* initalize value to 1st arg */
+    if (istype(a, ans)) {
+	/* 1st arg is hms object, load it */
+	ans.hour = a.hour;
+	ans.min = a.min;
+	ans.sec = a.sec;
+    } else {
+	/* 1st arg is not hms, assume scalar hours */
+	ans.hour = a;
+	ans.min = 0;
+	ans.sec = 0;
+    }
+
+    /* add value of 2nd arg */
+    if (istype(b, ans)) {
+	/* 2nd arg is hms object, add it */
+	ans.hour += b.hour;
+	ans.min += b.min;
+	ans.sec += b.sec;
+    } else {
+	/* 2nd arg is not hms, add scalar hours */
+	ans.hour += b;
+    }
+
+    /* return normalized result */
+    ans = fixhms(ans);
+    return ans;
+}
+
+
+define hms_neg(a)
+{
+    local obj hms ans;		/* return value */
+
+    /* negate argument */
+    if (istype(a, ans)) {
+	/* 1st arg is hms object, load it */
+	ans.hour = -a.hour;
+	ans.min = -a.min;
+	ans.sec = -a.sec;
+    } else {
+	/* 2nd arg is not hms, negate scalar hours */
+	ans.hour = -a;
+	ans.min = 0;
+	ans.sec = 0;
+    }
+
+    /* return normalized result */
+    ans = fixhms(ans);
+    return ans;
+}
+
+
+define hms_sub(a, b)
+{
+    local obj hms ans;		/* return value */
+
+    /* initalize value to 1st arg */
+    if (istype(a, ans)) {
+	/* 1st arg is hms object, load it */
+	ans.hour = a.hour;
+	ans.min = a.min;
+	ans.sec = a.sec;
+    } else {
+	/* 1st arg is not hms, assume scalar hours */
+	ans.hour = a;
+	ans.min = 0;
+	ans.sec = 0;
+    }
+
+    /* subtract value of 2nd arg */
+    if (istype(b, ans)) {
+	/* 2nd arg is hms object, subtract it */
+	ans.hour -= b.hour;
+	ans.min -= b.min;
+	ans.sec -= b.sec;
+    } else {
+	/* 2nd arg is not hms, subtract scalar hours */
+	ans.hour -= b;
+    }
+
+    /* return normalized result */
+    ans = fixhms(ans);
+    return ans;
+}
+
+
+define hms_mul(a, b)
+{
+    local obj hms ans;		/* return value */
+
+    /* hms object multiplication */
+    if (istype(a, ans) && istype(b, ans)) {
+	ans.hour = hms_abs(a) * hms_abs(b);
+	ans.min = 0;
+	ans.sec = 0;
+
+    /* scalar multiplication */
+    } else if (istype(a, ans)) {
+	ans.hour = a.hour * b;
+	ans.min = a.min * b;
+	ans.sec = a.sec * b;
+    } else {
+	ans.hour = b.hour * a;
+	ans.min = b.min * a;
+	ans.sec = b.sec * a;
+    }
+
+    /* return normalized result */
+    ans = fixhms(ans);
+    return ans;
+}
+
+
+define hms_print(a)
+{
+    local obj hms ans;		/* temp object for hms type testing */
+
+    /* firewall - arg must be a hms object */
+    if (! istype(a, ans)) {
+	quit "hms_print called with non hms object";
+    }
+
+    /* print in hms form */
+    print a.hour : ':' : a.min : ':' : a.sec :;
+}
+
+
+define hms_abs(a)
+{
+    local obj hms ans;		/* temp object for hms type testing */
+    local hour;			/* return scalar value */
+
+    /* firewall - just absolute value non hms objects */
+    if (! istype(a, ans)) {
+    	return abs(a);
+    }
+
+    /* compute hours */
+    hour = a.hour + a.min / 60 + a.sec / 3600;
+
+    /* return hours */
+    return hour;
+}
+
+
+define hms_norm(a)
+{
+    local obj hms ans;		/* temp object for hms type testing */
+    local hour;			/* hours */
+
+    /* firewall - arg must be a hms object */
+    if (! istype(a, ans)) {
+	quit "hms_norm called with non hms object";
+    }
+
+    /* square hours (norm is the square of absolute value */
+    hour = hms_abs(a);
+
+    /* return hours */
+    return hour*hour;
+}
+
+
+define hms_test(a)
+{
+    local obj hms ans;		/* temp value */
+
+    /* firewall - arg must be a hms object */
+    if (! istype(a, ans)) {
+	quit "hms_test called with non hms object";
+    }
+
+    /* return false of non-zero */
+    ans = fixhms(a);
+    if (ans.hour == 0 && ans.min == 0 && ans.sec == 0) {
+	/* false */
+	return 0;
+    }
+    /* true */
+    return 1;
+}
+
+
+define hms_int(a)
+{
+    local obj hms ans;		/* return value */
+
+    /* firewall - arg must be a hms object */
+    if (! istype(a, ans)) {
+	quit "hms_int called with non hms object";
+    }
+
+    /* normalize the argument */
+    ans = fixhms(a);
+
+    /* truncate to the nearest second */
+    ans.sec = int(ans.sec);
+
+    /* return value to the nearest second */
+    return ans;
+}
+
+
+define hms_frac(a)
+{
+    local obj hms ans;		/* return value */
+
+    /* firewall - arg must be a hms object */
+    if (! istype(a, ans)) {
+	quit "hms_frac called with non hms object";
+    }
+
+    /* normalize the argument */
+    ans = fixhms(a);
+
+    /* remove all but fractional seconds */
+    ans.hour = 0;
+    ans.min = 0;
+    ans.sec = frac(ans.sec);
+
+    /* return value to the second fraction */
+    return ans;
+}
+
+
+define hms_rel(a,b)
+{
+    local abs_a, abs_b;		/* scalars of the arguments */
+
+    /* compute scalars of the arguments */
+    abs_a = hms_abs(a);
+    abs_b = hms_abs(b);
+
+    /* return the comparison */
+    return cmp(abs_a, abs_b);
+}
+
+
+define hms_cmp(a,b)
+{
+    local abs_a, abs_b;		/* scalars of the arguments */
+
+    /* compute scalars of the arguments */
+    abs_a = hms_abs(a);
+    abs_b = hms_abs(b);
+
+    /* return the equality comparison */
+    return (abs_a == abs_b);
+}
+
+
+define hms_inc(a)
+{
+    local obj hms ans;		/* return value */
+
+    /* increment a hms object */
+    if (istype(a, ans)) {
+	ans = a;
+	++ans.sec;
+
+	/* return normalized result */
+	ans = fixhms(ans);
+	return ans;
+    }
+
+    /* increment a scalar */
+    return a+1;
+}
+
+
+define hms_dec(a)
+{
+    local obj hms ans;		/* return value */
+
+    /* decrement a hms object */
+    if (istype(a, ans)) {
+	ans = a;
+	--ans.sec;
+
+	/* return normalized result */
+	ans = fixhms(ans);
+	return ans;
+    }
+
+    /* decrement a scalar */
+    return a-1;
+}
+
+
+define fixhms(a)
+{
+    local obj hms ans;		/* temp value */
+
+    /* firewall */
+    if (! istype(a, ans)) {
+	quit "attempt to fix a non hms object";
+    }
+
+    /* force minutes to be intergral */
+    a.min += frac(a.hour) * 60;
+    a.hour = int(a.hour);
+
+    /* force hours to be intergral */
+    a.sec += frac(a.min) * 60;
+    a.min = int(a.min);
+
+    /* carry excess seconds into minutes */
+    a.min += a.sec // 60;
+    a.sec %= 60;
+
+    /* carry excess minutes into hours */
+    a.hour += a.min // 60;
+    a.min %= 60;
+
+    /* round hours by day */
+    a.hour %= 24;
+
+    /* return normalized result */
+    return a;
+}
+
+if (config("resource_debug") & 3) {
+    print "obj hms {hour, min, sec} defined";
+}

cal/infinities.cal

@@ -0,0 +1,88 @@
+/*
+ * infinities - handle infinities symbolically, a little helper file
+ *
+ * Copyright (C) 2013 Christoph Zurnieden
+ *
+ * Calc is open software; you can redistribute it and/or modify it under
+ * the terms of the version 2.1 of the GNU Lesser General Public License
+ * as published by the Free Software Foundation.
+ *
+ * Calc is distributed in the hope that it will be useful, but WITHOUT
+ * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
+ * or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General
+ * Public License for more details.
+ *
+ * A copy of version 2.1 of the GNU Lesser General Public License is
+ * distributed with calc under the filename COPYING-LGPL.  You should have
+ * received a copy with calc; if not, write to Free Software Foundation, Inc.
+ * 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
+ */
+
+
+static resource_debug_level;
+resource_debug_level = config("resource_debug", 0);
+
+
+define isinfinite(x)
+{
+    if (isstr(x)) {
+	if (strncmp(x, "cinf", 4) == 0
+	    || strncmp(x, "pinf", 4) == 0 || strncmp(x, "ninf", 4) == 0)
+	    return 1;
+    }
+    return 0;
+}
+
+define iscinf(x)
+{
+    if (isstr(x)) {
+	if (strncmp(x, "cinf", 4) == 0)
+	    return 1;
+    }
+    return 0;
+}
+
+define ispinf(x)
+{
+    if (isstr(x)) {
+	if (strncmp(x, "pinf", 4) == 0)
+	    return 1;
+    }
+    return 0;
+}
+
+define isninf(x)
+{
+    if (isstr(x)) {
+	if (strncmp(x, "ninf", 4) == 0)
+	    return 1;
+    }
+    return 0;
+}
+
+define cinf()
+{
+    return "cinf";
+}
+
+define ninf()
+{
+    return "ninf";
+}
+
+define pinf()
+{
+    return "pinf";
+}
+
+
+config("resource_debug", resource_debug_level),;
+if (config("resource_debug") & 3) {
+    print "isinfinite(x)";
+    print "iscinf(x)";
+    print "ispinf(x)";
+    print "isninf(x)";
+    print "cinf()";
+    print "ninf()";
+    print "pinf()";
+}

cal/intfile.cal

@@ -19,7 +19,7 @@
  *
  * @(#) $Revision: 30.1 $
  * @(#) $Id: intfile.cal,v 30.1 2007/03/16 11:09:54 chongo Exp $
- * @(#) $Source: /usr/local/src/cmd/calc/cal/RCS/intfile.cal,v $
+ * @(#) $Source: /usr/local/src/bin/calc/cal/RCS/intfile.cal,v $
  *
  * Under source code control:   2001/03/31 08:13:11
  * File existed as early as:    2001

cal/intnum.cal

@@ -0,0 +1,728 @@
+/*
+ * intnum - implementation of tanhsinh- and Gauss-Legendre quadrature
+ *
+ * Copyright (C) 2013 Christoph Zurnieden
+ *
+ * Calc is open software; you can redistribute it and/or modify it under
+ * the terms of the version 2.1 of the GNU Lesser General Public License
+ * as published by the Free Software Foundation.
+ *
+ * Calc is distributed in the hope that it will be useful, but WITHOUT
+ * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
+ * or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General
+ * Public License for more details.
+ *
+ * A copy of version 2.1 of the GNU Lesser General Public License is
+ * distributed with calc under the filename COPYING-LGPL.  You should have
+ * received a copy with calc; if not, write to Free Software Foundation, Inc.
+ * 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
+ */
+
+
+static resource_debug_level;
+resource_debug_level = config("resource_debug", 0);
+
+
+read -once infinities;
+
+static __CZ__tanhsinh_x;
+static __CZ__tanhsinh_w;
+static __CZ__tanhsinh_order;
+static __CZ__tanhsinh_prec;
+
+define quadtsdeletenodes()
+{
+    free(__CZ__tanhsinh_x);
+    free(__CZ__tanhsinh_w);
+    free(__CZ__tanhsinh_order);
+    free(__CZ__tanhsinh_prec);
+}
+
+define quadtscomputenodes(order, expo, eps)
+{
+    local t cht sht chp sum k PI places;
+    local h t0 x w;
+    if (__CZ__tanhsinh_order == order && __CZ__tanhsinh_prec == eps)
+	return 1;
+    __CZ__tanhsinh_order = order;
+    __CZ__tanhsinh_prec = eps;
+    __CZ__tanhsinh_x = list();
+    __CZ__tanhsinh_w = list();
+    /* The tanhsinh algorithm needs a slightly higher precision than G-L */
+    eps = epsilon(eps * 1e-2);
+    places = highbit(1 + int (1 / epsilon())) +1;
+    PI = pi();
+    sum = 0;
+    t0 = 2 ^ (-expo);
+    h = 2 * t0;
+    /*
+     * The author wanted to use the mpmath trick here which was
+     * advertised---and reasonably so!---to be faster. Didn't work out
+     * so well with calc.
+     * PI4 = PI/4;
+     * expt0 = bround(exp(t0),places);
+     * a = bround( PI4 * expt0,places);
+     * b = bround(PI4 / expt0,places);
+     * udelta = bround(exp(h),places);
+     * urdelta = bround(1/udelta,places);
+     */
+    /* make use of x(-t) = -x(t), w(-t) = w(t)  */
+    for (k = 0; k < 20 * order + 1; k++) {
+	/*
+	 * x = tanh(pi/2 * sinh(t))
+	 * w = pi/2 * cosh(t) / cosh(pi/2 * sinh(t))^2
+	 */
+	t = bround(t0 + k * h, places);
+
+	cht = bround(cosh(t), places);
+	sht = bround(sinh(t), places);
+	chp = bround(cosh(0.5 * PI * sht), places);
+	x = bround(tanh(0.5 * PI * sht), places);
+	w = bround((PI * h * cht) / (2 * chp ^ 2), places);
+	/*
+	 * c = bround(exp(a-b),places);
+	 * d = bround(1/c,places);
+	 * co =bround( (c+d)/2,places);
+	 * si =bround( (c-d)/2,places);
+	 * x = bround(si / co,places);
+	 * w = bround((a+b) / co^2,places);
+	 */
+	if (abs(x - 1) <= eps)
+	    break;
+
+	append(__CZ__tanhsinh_x, x);
+	append(__CZ__tanhsinh_w, w);
+	/*
+	 * a *= udelta;
+	 * b *= urdelta;
+	 */
+    }
+
+    /* Normalize the weights to make them add up to 2 (two) */
+    /*
+     * for(k=0;k < size(__CZ__tanhsinh_w);k++)
+     * sum = bround(sum + __CZ__tanhsinh_w[k],places);
+     * sum *= 2;
+     * for(k=0;k < size(__CZ__tanhsinh_w);k++)
+     * __CZ__tanhsinh_w[k] = bround(2.0 * __CZ__tanhsinh_w[k] / sum,places);
+     */
+
+    epsilon(eps);
+    return 1;
+}
+
+define quadtscore(a, b, n)
+{
+    local k c d order eps places sum ret x x1 x2 xm w w1 w2 m sizel;
+
+    eps = epsilon(epsilon() * 1e-2);
+    places = highbit(1 + int (1 / epsilon())) +1;
+    m = int (4 + max(0, ln(places / 30.0) / ln(2))) + 2;
+    if (!isnull(n)) {
+	order = n;
+	m = ilog(order / 3, 2) + 1;
+    } else
+	order = 3 * 2 ^ (m - 1);
+
+    quadtscomputenodes(order, m, epsilon());
+    sizel = size(__CZ__tanhsinh_w);
+
+    if (isinfinite(a) || isinfinite(b)) {
+	/*
+	 *           x
+	 * t =  ------------
+	 *                2
+	 *      sqrt(1 - y )
+	 */
+	if (isninf(a) && ispinf(b)) {
+	    for (k = 0; k < sizel; k++) {
+		x1 = __CZ__tanhsinh_x[k];
+		x2 = -__CZ__tanhsinh_x[k];
+		w1 = __CZ__tanhsinh_w[k];
+
+		x = bround(x1 * (1 - x1 ^ 2) ^ (-1 / 2), places);
+		xm = bround(x2 * (1 - x2 ^ 2) ^ (-1 / 2), places);
+		w = bround(w1 * (((1 - x1 ^ 2) ^ (-1 / 2)) / (1 - x1 ^ 2)),
+			   places);
+		w2 = bround(w1 * (((1 - x2 ^ 2) ^ (-1 / 2)) / (1 - x2 ^ 2)),
+			    places);
+		sum += bround(w * f(x), places);
+		sum += bround(w2 * f(xm), places);
+	    }
+	}
+	/*
+	 *        1
+	 * t =  - - + b + 1
+	 *        x
+	 */
+	else if (isninf(a) && !iscinf(b)) {
+	    for (k = 0; k < sizel; k++) {
+		x1 = __CZ__tanhsinh_x[k];
+		x2 = -__CZ__tanhsinh_x[k];
+		w1 = __CZ__tanhsinh_w[k];
+
+		x = bround((b + 1) - (2 / (x1 + 1)), places);
+		xm = bround((b + 1) - (2 / (x2 + 1)), places);
+		w = bround(w1 * (1 / 2 * (2 / (x1 + 1)) ^ 2), places);
+		w2 = bround(w1 * (1 / 2 * (2 / (x2 + 1)) ^ 2), places);
+		sum += bround(w * f(x), places);
+		sum += bround(w2 * f(xm), places);
+	    }
+	}
+	/*
+	 *      1
+	 * t =  - + a - 1
+	 *      x
+	 */
+	else if (!iscinf(a) && ispinf(b)) {
+	    for (k = 0; k < sizel; k++) {
+		x1 = __CZ__tanhsinh_x[k];
+		x2 = -__CZ__tanhsinh_x[k];
+		w1 = __CZ__tanhsinh_w[k];
+		x = bround((a - 1) + (2 / (x1 + 1)), places);
+		xm = bround((a - 1) + (2 / (x2 + 1)), places);
+		w = bround(w1 * (((1 / 2) * (2 / (x1 + 1)) ^ 2)), places);
+		w2 = bround(w1 * (((1 / 2) * (2 / (x2 + 1)) ^ 2)), places);
+		sum += bround(w * f(x), places);
+		sum += bround(w2 * f(xm), places);
+	    }
+	} else if (isninf(a) || isninf(b)) {
+	    /*TODO: swap(a,b) and negate(w)? Lookup! */
+	    return newerror("quadtscore: reverse limits?");
+	} else {
+	    return
+		newerror("quadtscore: complex infinity not yet implemented");
+        }
+	ret = sum;
+    } else {
+	/* Avoid rounding errors */
+	if (a == -1 && b == 1) {
+	    c = 1;
+	    d = 0;
+	} else {
+	    c = (b - a) / 2;
+	    d = (b + a) / 2;
+	}
+	sum = 0;
+	for (k = 0; k < sizel; k++) {
+	    sum +=
+		bround(__CZ__tanhsinh_w[k] * f(c * __CZ__tanhsinh_x[k] + d),
+		       places);
+	    sum +=
+		bround(__CZ__tanhsinh_w[k] * f(c * -__CZ__tanhsinh_x[k] + d),
+		       places);
+	}
+	ret = c * sum;
+    }
+    epsilon(eps);
+    return ret;
+}
+
+static __CZ__quadts_error;
+
+define quadts(a, b, points)
+{
+    local k sp results epsbits nsect interval length segment slope C ;
+    local x1 x2 y1 y2  sum D1 D2 D3 D4;
+    if (param(0) < 2)
+	return newerror("quadts: not enough arguments");
+    epsbits = highbit(1 + int (1 / epsilon())) +1;
+    if (param(0) < 3 || isnull(points)) {
+	/* return as given */
+	return quadtscore(a, b);
+    } else {
+	if ((isinfinite(a) || isinfinite(b))
+	    && (!ismat(points) && !islist(points)))
+	    return
+		newerror(strcat
+			 ("quadts: segments of infinite length ",
+			  "are not yet supported"));
+	if (ismat(points) || islist(points)) {
+	    sp = size(points);
+	    if (sp == 0)
+		return
+		    newerror(strcat
+			     ("quadts: variable 'points` must be a list or ",
+			      "1d-matrix of a length > 0"));
+	    /* check if all points are numbers */
+	    for (k = 0; k < sp; k++) {
+		if (!isnum(points[k]))
+		    return
+			newerror(strcat
+				 ("quadts: elements of 'points` must be",
+				  " numbers only"));
+	    }
+	    /* We have n-1 intervals and a and b, hence n-1 + 2 results */
+	    results = mat[sp + 1];
+            if (a != points[0]) {
+	        results[0] = quadtscore(a, points[0]);
+            } else {
+                results[0] = 0;
+            }
+	    if (sp == 1) {
+                if (b != points[0]) {
+		    results[1] = quadtscore(points[0], b);
+                } else {
+                    results[1] = 0;
+                }
+	    } else {
+		for (k = 1; k < sp; k++) {
+		    results[k] = quadtscore(points[k - 1], points[k]);
+		}
+                if (b != points[k - 1]) {
+		    results[k] = quadtscore(points[k - 1], b);
+                } else {
+                    results[k] = 0;
+                }
+	    }
+	} else {
+	    if (!isint(points) || points <= 0)
+		return newerror(strcat("quadts: variable 'points` must be a ",
+				       "list or a positive integer"));
+	    /* Taking "points" as the number of equally spaced intervals */
+	    results = mat[points + 1];
+	    /* It is easy if a,b lie on the real line */
+	    if (isreal(a) && isreal(b)) {
+		length = abs(a - b);
+		segment = length / points;
+
+		for (k = 1; k <= points; k++) {
+		    results[k - 1] =
+			quadtscore(a + (k - 1) * segment, a + k * segment);
+		}
+	    } else {
+		/* We have at least one complex limit but treat "points" still
+                 * as the number of equally spaced intervals on a straight line
+                 * connecting a and b. Computing the segments here is a bit
+                 * more complicated but not much, it should have been taught in
+                 * highschool.
+		 * Other contours by way of a list of points */
+		slope = (im(b) - im(a)) / (re(b) - re(a));
+		C = (im(a) + slope) * re(a);
+		length = abs(re(a) - re(b));
+		segment = length / points;
+
+		/* y = mx+C where m is the slope, x is the real part and y the
+		 * imaginary part  */
+                if(re(a)>re(b))swap(a,b);
+		for (k = re(a); k <= (re(b)); k+=segment) {
+		    x1 = slope*(k) +  C;
+		    results[k] = quadtscore(k + x1 * 1i);
+		}
+	    }			/* else of isreal */
+	}			/* else of ismat|islist */
+    }				/* else of isnull(points) */
+    /* With a bit of undeserved luck we have a result by now. */
+    sp = size(results);
+    for (k = 0; k < sp; k++) {
+	sum += results[k];
+    }
+    return sum;
+}
+
+static __CZ__gl_x;
+static __CZ__gl_w;
+static __CZ__gl_order;
+static __CZ__gl_prec;
+
+define quadglcomputenodes(N)
+{
+    local places k l x w t1 t2 t3 t4 t5 r tmp;
+
+    if (__CZ__gl_order == N && __CZ__gl_prec == epsilon())
+	return;
+
+    __CZ__gl_x = mat[N];
+    __CZ__gl_w = mat[N];
+    __CZ__gl_order = N;
+    __CZ__gl_prec = epsilon();
+
+    places = highbit(1 + int (1 / epsilon())) +1;
+
+    /*
+     * Compute roots and weights (doing it inline seems to be fastest)
+     * Trick shamelessly stolen from D. Bailey et .al (program "arprec")
+     */
+    for (k = 1; k <= N//2; k++) {
+	r = bround(cos(pi() * (k - .25) / (N + .5)), places);
+	while (1) {
+	    t1 = 1, t2 = 0;
+	    for (l = 1; l <= N; l++) {
+		t3 = t2;
+		t2 = t1;
+		t1 = bround(((2 * l - 1) * r * t2 - (l - 1) * t3) / l, places);
+	    }
+	    t4 = bround(N * (r * t1 - t2) / ((r ^ 2) - 1), places);
+	    t5 = r;
+	    tmp = t1 / t4;
+	    r = r - tmp;
+	    if (abs(tmp) <= epsilon())
+		break;
+	}
+	x = r;
+	w = bround(2 / ((1 - r ^ 2) * t4 ^ 2), places);
+
+	__CZ__gl_x[k - 1] = x;
+	__CZ__gl_w[k - 1] = w;
+	__CZ__gl_x[N - k] = -__CZ__gl_x[k - 1];
+	__CZ__gl_w[N - k] = __CZ__gl_w[k - 1];
+    }
+    return;
+}
+
+define quadgldeletenodes()
+{
+    free(__CZ__gl_x);
+    free(__CZ__gl_w);
+    free(__CZ__gl_order);
+    free(__CZ__gl_prec);
+}
+
+define quadglcore(a, b, n)
+{
+    local k c d digs order eps places sum ret err x x1 w w1 m;
+    local phalf x2 px1 spx1 u b1 a1 half;
+
+    eps = epsilon(epsilon() * 1e-2);
+    places = highbit(1 + int (1 / epsilon())) +1;
+    if (!isnull(n))
+	order = n;
+    else {
+	m = int (4 + max(0, ln(places / 30.0) / ln(2))) + 2;
+	order = 3 * 2 ^ (m - 1);
+    }
+
+
+    quadglcomputenodes(order, 1);
+
+    if (isinfinite(a) || isinfinite(b)) {
+	if (isninf(a) && ispinf(b)) {
+	    for (k = 0; k < order; k++) {
+		x1 = __CZ__gl_x[k];
+		w1 = __CZ__gl_w[k];
+
+		x = bround(x1 * (1 - x1 ^ 2) ^ (-1 / 2), places);
+		w = bround(w1 * (((1 - x1 ^ 2) ^ (-1 / 2)) / (1 - x1 ^ 2)),
+			   places);
+		sum += bround(w * f(x), places);
+	    }
+	} else if (isninf(a) && !iscinf(b)) {
+	    for (k = 0; k < order; k++) {
+		x1 = __CZ__gl_x[k];
+		w1 = __CZ__gl_w[k];
+
+		x = bround((b + 1) - (2 / (x1 + 1)), places);
+		w = bround(w1 * (1 / 2 * (2 / (x1 + 1)) ^ 2), places);
+		sum += bround(w * f(x), places);
+	    }
+	} else if (!iscinf(a) && ispinf(b)) {
+	    for (k = 0; k < order; k++) {
+		x1 = __CZ__gl_x[k];
+		w1 = __CZ__gl_w[k];
+		x = bround((a - 1) + (2 / (x1 + 1)), places);
+		w = bround(w1 * (((1 / 2) * (2 / (x1 + 1)) ^ 2)), places);
+		sum += bround(w * f(x), places);
+	    }
+	} else if (isninf(a) || isninf(b)) {
+	    /*TODO: swap(a,b) and negate(w)? Lookup! */
+	    return newerror("quadglcore: reverse limits?");
+	} else
+	    return
+		newerror("quadglcore: complex infinity not yet implemented");
+	ret = sum;
+    } else {
+	/* Avoid rounding errors */
+	if (a == -1 && b == 1) {
+	    c = 1;
+	    d = 0;
+	} else {
+	    c = (b - a) / 2;
+	    d = (b + a) / 2;
+	}
+	sum = 0;
+	for (k = 0; k < order; k++) {
+	    sum += bround(__CZ__gl_w[k] * f(c * __CZ__gl_x[k] + d), places);
+	}
+	ret = c * sum;
+    }
+    epsilon(eps);
+    return ret;
+}
+
+define quadgl(a, b, points)
+{
+    local k sp results epsbits nsect interval length segment slope C x1 y1 x2
+	y2;
+    local sum D1 D2 D3 D4;
+    if (param(0) < 2)
+	return newerror("quadgl: not enough arguments");
+    epsbits = highbit(1 + int (1 / epsilon())) +1;
+    if (isnull(points)) {
+	/* return as given */
+	return quadglcore(a, b);
+    } else {
+	/* But if we could half the time needed to execute a single operation
+	 * we could do all of it in just twice that time. */
+	if (isinfinite(a) || isinfinite(b)
+	    && (!ismat(points) && !islist(points)))
+	    return
+		newerror(strcat
+			 ("quadgl: multiple segments of infinite length ",
+			  "are not yet supported"));
+	if (ismat(points) || islist(points)) {
+	    sp = size(points);
+	    if (sp == 0)
+		return
+		    newerror(strcat
+			     ("quadgl: variable 'points` must be a list or ",
+			      "1d-matrix of a length > 0"));
+	    /* check if all points are numbers */
+	    for (k = 0; k < sp; k++) {
+		if (!isnum(points[k]))
+		    return
+			newerror(strcat
+				 ("quadgl: elements of 'points` must be ",
+				  "numbers only"));
+	    }
+	    /* We have n-1 intervals and a and b, hence n-1 + 2 results */
+	    results = mat[sp + 1];
+            if (a != points[0]) {
+	        results[0] = quadglcore(a, points[0]);
+            } else {
+                results[0] = 0;
+            }
+	    if (sp == 1) {
+		if (b != points[0]) {
+		    results[1] = quadglcore(points[0], b);
+                } else {
+                    results[1] = 0;
+                }
+	    } else {
+		for (k = 1; k < sp; k++) {
+		    results[k] = quadglcore(points[k - 1], points[k]);
+		}
+                if (b != points[k - 1]) {
+		    results[k] = quadglcore(points[k - 1], b);
+                } else {
+                    results[k] = 0;
+                }
+	    }
+	} else {
+	    if (!isint(points) || points <= 0)
+		return newerror(strcat("quadgl: variable 'points` must be a ",
+				       "list or a positive integer"));
+	    /* Taking "points" as the number of equally spaced intervals */
+	    results = mat[points + 1];
+	    /* It is easy if a,b lie on the real line */
+	    if (isreal(a) && isreal(b)) {
+		length = abs(a - b);
+		segment = length / points;
+
+		for (k = 1; k <= points; k++) {
+		    results[k - 1] =
+			quadglcore(a + (k - 1) * segment, a + k * segment);
+		}
+	    } else {
+		/* Other contours by way of a list of points */
+		slope = (im(b) - im(a)) / (re(b) - re(a));
+		C = (im(a) + slope) * re(a);
+		length = abs(re(a) - re(b));
+		segment = length / points;
+
+		/* y = mx+C where m is the slope, x is the real part and y the
+		 * imaginary part  */
+                if(re(a)>re(b))swap(a,b);
+		for (k = re(a); k <= (re(b)); k+=segment) {
+		    x1 = slope*(k) +  C;
+		    results[k] = quadglcore(k + x1 * 1i);
+		}
+	    }			/* else of isreal */
+	}			/* else of ismat|islist */
+    }				/* else of isnull(points) */
+    /* With a bit of undeserved luck we have a result by now. */
+    sp = size(results);
+    for (k = 0; k < sp; k++) {
+	sum += results[k];
+    }
+    return sum;
+}
+
+define quad(a, b, points = -1, method = "tanhsinh")
+{
+    if (isnull(a) || isnull(b) || param(0) < 2)
+	return newerror("quad: both limits must be given");
+    if (isstr(a)) {
+	if (strncmp(a, "cinf", 1) == 0)
+	    return
+		newerror(strcat
+			 ("quad: complex infinity not yet supported, use",
+			  " 'pinf' or 'ninf' respectively"));
+    }
+    if (isstr(b)) {
+	if (strncmp(b, "cinf", 1) == 0)
+	    return
+		newerror(strcat
+			 ("quad: complex infinity not yet supported, use",
+			  " 'pinf' or 'ninf' respectively"));
+    }
+
+    if (param(0) == 3) {
+	if (isstr(points))
+	    method = points;
+    }
+
+    if (strncmp(method, "tanhsinh", 1) == 0) {
+	if (!isstr(points)) {
+	    if (points == -1) {
+		return quadts(a, b);
+	    } else {
+		return quadts(a, b, points);
+	    }
+	} else {
+	    return quadts(a, b);
+	}
+    }
+
+    if (strncmp(method, "gausslegendre", 1) == 0) {
+	if (!isstr(points)) {
+	    if (points == -1) {
+		return quadgl(a, b);
+	    } else {
+		return quadgl(a, b, points);
+	    }
+	} else {
+	    return quadgl(a, b);
+	}
+    }
+}
+
+define makerange(start, end, steps)
+{
+    local ret k l step C length slope x1 x2 y1 y2;
+    local segment;
+    steps = int (steps);
+    if (steps < 1) {
+	return newerror("makerange: number of steps must be > 0");
+    }
+    if (!isnum(start) || !isnum(end)) {
+	return newerror("makerange: only numbers are supported yet");
+    }
+    if (isreal(start) && isreal(end)) {
+	step = (end - start) / (steps);
+	print step;
+	ret = mat[steps + 1];
+	for (k = 0; k <= steps; k++) {
+	    ret[k] = k * step + start;
+	}
+    } else {
+	ret = mat[steps + 1];
+	if (re(start) > re(end)) {
+	    swap(start, end);
+	}
+
+	slope = (im(end) - im(start)) / (re(end) - re(start));
+	C = im(start) - slope * re(start);
+	length = abs(re(start) - re(end));
+	segment = length / (steps);
+
+	for (k = re(start), l = 0; k <= (re(end)); k += segment, l++) {
+	    x1 = slope * (k) + C;
+	    ret[l] = k + x1 * 1i;
+	}
+
+    }
+    return ret;
+}
+
+define makecircle(radius, center, points)
+{
+    local ret k a b twopi centerx centery;
+    if (!isint(points) || points < 2) {
+	return
+	    newerror("makecircle: number of points is not a positive integer");
+    }
+    if (!isnum(center)) {
+	return newerror("makecircle: center does not lie on the complex plane");
+    }
+    if (!isreal(radius) || radius <= 0) {
+	return newerror("makecircle: radius is not a real > 0");
+    }
+    ret = mat[points];
+    twopi = 2 * pi();
+    centerx = re(center);
+    centery = im(center);
+    for (k = 0; k < points; k++) {
+	a = centerx + radius * cos(twopi * k / points);
+	b = centery + radius * sin(twopi * k / points);
+	ret[k] = a + b * 1i;
+    }
+    return ret;
+}
+
+define makeellipse(angle, a, b, center, points)
+{
+    local ret k x y twopi centerx centery;
+    if (!isint(points) || points < 2) {
+	return
+	    newerror("makeellipse: number of points is not a positive integer");
+    }
+    if (!isnum(center)) {
+	return
+	    newerror("makeellipse: center does not lie on the complex plane");
+    }
+    if (!isreal(a) || a <= 0) {
+	return newerror("makecircle: a is not a real > 0");
+    }
+    if (!isreal(b) || b <= 0) {
+	return newerror("makecircle: b is not a real > 0");
+    }
+    if (!isreal(angle)) {
+	return newerror("makecircle: angle is not a real");
+    }
+    ret = mat[points];
+    twopi = 2 * pi();
+    centerx = re(center);
+    centery = im(center);
+    for (k = 0; k < points; k++) {
+	x = centerx + a * cos(twopi * k / points) * cos(angle)
+	    - b * sin(twopi * k / points) * sin(angle);
+	y = centerx + a * cos(twopi * k / points) * sin(angle)
+	    + b * sin(twopi * k / points) * cos(angle);
+	ret[k] = x + y * 1i;
+    }
+    return ret;
+}
+
+define makepoints()
+{
+    local ret k;
+    ret = mat[param(0)];
+    for (k = 0; k < param(0); k++) {
+	if (!isnum(param(k + 1))) {
+	    return
+		newerror(strcat
+			 ("makepoints: parameter number \"", str(k + 1),
+			  "\" is not a number"));
+	}
+	ret[k] = param(k + 1);
+    }
+    return ret;
+}
+
+
+config("resource_debug", resource_debug_level),;
+if (config("resource_debug") & 3) {
+    print "quadtsdeletenodes()";
+    print "quadtscomputenodes(order, expo, eps)";
+    print "quadtscore(a,b,n)";
+    print "quadts(a,b,points)";
+    print "quadglcomputenodes(N)";
+    print "quadgldeletenodes()";
+    print "quadglcore(a,b,n)";
+    print "quadgl(a,b,points)";
+    print "quad(a,b,points=-1,method=\"tanhsinh\")";
+    print "makerange(start, end, steps)";
+    print "makecircle(radius, center, points)";
+    print "makeellipse(angle, a, b, center, points)";
+    print "makepoints(a1,[...])";
+}

cal/lambertw.cal

@@ -0,0 +1,288 @@
+/*
+ * lambertw - Lambert's W-function
+ *
+ * Copyright (C) 2013 Christoph Zurnieden
+ *
+ * Calc is open software; you can redistribute it and/or modify it under
+ * the terms of the version 2.1 of the GNU Lesser General Public License
+ * as published by the Free Software Foundation.
+ *
+ * Calc is distributed in the hope that it will be useful, but WITHOUT
+ * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
+ * or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General
+ * Public License for more details.
+ *
+ * A copy of version 2.1 of the GNU Lesser General Public License is
+ * distributed with calc under the filename COPYING-LGPL.  You should have
+ * received a copy with calc; if not, write to Free Software Foundation, Inc.
+ * 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
+ *
+ * @(#) $Revision: 30.4 $
+ * @(#) $Id: lambertw.cal,v 30.4 2013/08/18 20:01:53 chongo Exp $
+ * @(#) $Source: /usr/local/src/bin/calc/cal/RCS/lambertw.cal,v $
+ *
+ * Under source code control:	2013/08/11 01:31:28
+ * File existed as early as:	2013
+ */
+
+
+static resource_debug_level;
+resource_debug_level = config("resource_debug", 0);
+
+
+/*
+
+     R. M. Corless and G. H. Gonnet and D. E. G. Hare and D. J. Jeffrey and
+     D. E. Knuth, "On the Lambert W Function", Advances n Computational
+     Mathematics, 329--359, (1996)
+     http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.112.6117
+
+     D. J. Jeffrey, D. E. G. Hare, R. M. Corless, "Unwinding the branches of the
+     Lambert W function", The Mathematical Scientist, 21, pp 1-7, (1996)
+     http://www.apmaths.uwo.ca/~djeffrey/Offprints/wbranch.pdf
+
+     Darko Verebic, "Having Fun with Lambert W(x) Function"
+     arXiv:1003.1628v1, March 2010, http://arxiv.org/abs/1003.1628
+
+     Winitzki, S. "Uniform Approximations for Transcendental Functions",
+     In Part 1 of Computational Science and its Applications - ICCSA 2003,
+     Lecture Notes in Computer Science, Vol. 2667, Springer-Verlag,
+     Berlin, 2003, 780-789. DOI 10.1007/3-540-44839-X_82
+     A copy may be found by Google.
+
+
+*/
+static true = 1;
+static false = 0;
+
+/* Branch 0, Winitzki (2003) , the well known Taylor series*/
+define __CZ__lambertw_0(z,eps){
+  local a=2.344e0, b=0.8842e0, c=0.9294e0, d=0.5106e0, e=-1.213e0;
+  local y=sqrt(2*exp(1)*z+2);
+  return (2*ln(1+b*y)-ln(1+c*ln(1+d*y))+e)/(1+1/(2*ln(1+b*y)+2*a));
+}
+/* branch -1 */
+define __CZ__lambertw_m1(z,eps){
+  local wn k;
+  /* Cut-off found in Maxima */
+  if(z < 0.3) return __CZ__lambertw_app(z,eps);
+  wn = z;
+  /*  Verebic (2010) eqs. 16-18*/
+  for(k=0;k<10;k++){
+    wn = ln(-z)-ln(-wn);
+  }
+  return wn;
+}
+
+/*
+  generic approximation
+
+  series for 1+W((z-2)/(2 e))
+
+  Corless et al (1996) (4.22)
+  Verebic (2010) eqs. 35-37; more coefficients given at the end of sect. 3.1
+  or online
+     http://www.wolframalpha.com/input/?
+       i=taylor+%28+1%2Bproductlog%28+%28z-2%29%2F%282*e%29+%29+%29
+  or by using the function lambertw_series_print() after running
+  lambertw_series(z,eps,branch,terms) at least once with the wanted number of
+  terms and z = 1 (which might throw an error because the series will not
+  converge in anybodies lifetime for something that far from the branchpoint).
+
+
+*/
+define __CZ__lambertw_app(z,eps){
+  local  b0=-1, b1=1, b2=-1/3, b3=11/72;
+  local  y=sqrt(2*exp(1)*z+2);
+  return b0 + ( y * (b1 + (y * (b2 + (b3 * y)))));
+}
+
+static __CZ__Ws_a;
+static __CZ__Ws_c;
+static __CZ__Ws_len=0;
+
+define lambertw_series_print(){
+  local k;
+  for(k=0;k<__CZ__Ws_len;k++){
+    print num(__CZ__Ws_c[k]):"/":den(__CZ__Ws_c[k]):"*p^":k;
+  }
+}
+
+/*
+   The series is fast but only if _very_ close to the branchpoint
+   The exact branch must be given explicitly, e.g.:
+
+; lambertw(-exp(-1)+.001)-lambertw_series(-exp(-1)+.001,epsilon()*1e-10,0)
+	-0.14758879113205794065490184399030194122136720202792-
+                  0.00000000000000000000000000000000000000000000000000i
+; lambertw(-exp(-1)+.001)-lambertw_series(-exp(-1)+.001,epsilon()*1e-10,1)
+	0.00000000000000000000000000000000000000000000000000-
+                  0.00000000000000000000000000000000000000000000000000i
+ */
+define lambertw_series(z,eps,branch,terms){
+  local k l limit tmp sum A C P PP epslocal;
+  if(!isnull(terms))
+    limit = terms;
+  else
+    limit = 100;
+
+  if(isnull(eps))
+    eps = epsilon(epsilon()*1e-10);
+  epslocal = epsilon(eps);
+
+  P = sqrt(2*(exp(1)*z+1));
+  if(branch != 0) P = -P;
+  tmp=0;sum=0;PP=P;
+
+    __CZ__Ws_a   = mat[limit+1];
+    __CZ__Ws_c   = mat[limit+1];
+    __CZ__Ws_len = limit;
+   /*
+      c0 = -1; c1 = 1
+      a0 =  2; a1 =-1
+    */
+    __CZ__Ws_c[0] = -1; __CZ__Ws_c[1] =  1;
+    __CZ__Ws_a[0] =  2; __CZ__Ws_a[1] = -1;
+    sum += __CZ__Ws_c[0];
+    sum += __CZ__Ws_c[1] * P;
+    PP *= P;
+    for(k=2;k<limit;k++){
+      for(l=2;l<k;l++){
+         __CZ__Ws_a[k] += __CZ__Ws_c[l]*__CZ__Ws_c[k+1-l];
+      }
+
+      __CZ__Ws_c[k] = (k-1) * (  __CZ__Ws_c[k-2]/2
+                            +__CZ__Ws_a[k-2]/4)/
+                      (k+1)-__CZ__Ws_a[k]/2-__CZ__Ws_c[k-1]/(k+1);
+      tmp = __CZ__Ws_c[k] * PP;
+      sum += tmp;
+      if(abs(tmp) <= eps){
+        epsilon(epslocal);
+        return sum;
+      }
+      PP *= P;
+    }
+    epsilon(epslocal);
+    return
+      newerror(strcat("lambertw_series: does not converge in ",
+                        str(limit)," terms" ));
+}
+
+/* */
+define lambertw(z,branch){
+  local eps epslarge ret branchpoint bparea w we ew w1e wn k places m1e;
+  local closeness;
+
+  eps = epsilon();
+  if(branch == 0){
+    if(!im(z)){
+      if(abs(z) <= eps) return 0;
+      if(abs(z-exp(1)) <= eps) return 1;
+      if(abs(z - (-ln(2)/2)) <= eps ) return -ln(2);
+      if(abs(z - (-pi()/2)) <= eps ) return 1i*pi()/2;
+    }
+  }
+
+  branchpoint = -exp(-1);
+  bparea = .2;
+  if(branch == 0){
+    if(!im(z) && abs(z-branchpoint) == 0) return -1;
+    ret = __CZ__lambertw_0(z,eps);
+   /* Yeah, C&P, I know, sorry */
+   ##ret = ln(z) + 2*pi()*1i*branch - ln(ln(z)+2*pi()*1i*branch);
+  }
+  else if(branch == 1){
+    if(im(z)<0 && abs(z-branchpoint) <= bparea)
+      ret = __CZ__lambertw_app(z,eps);
+   /* Does calc have a goto? Oh, it does!  */
+    ret =ln(z) + 2*pi()*1i*branch - ln(ln(z)+2*pi()*1i*branch);
+  }
+  else if(branch == -1){##print "-1";
+    if(!im(z) && abs(z-branchpoint) == 0) return -1;
+    if(!im(z) && z>branchpoint && z < 0){##print "0";
+      ret = __CZ__lambertw_m1(z,eps);}
+    if(im(z)>=0 && abs(z-branchpoint) <= bparea){##print "1";
+      ret = __CZ__lambertw_app(z,eps);}
+    ret =ln(z) + 2*pi()*1i*branch - ln(ln(z)+2*pi()*1i*branch);
+  }
+  else
+    ret = ln(z) + 2*pi()*1i*branch - ln(ln(z)+2*pi()*1i*branch);
+
+  /*
+     Such a high precision is only needed _very_ close to the branchpoint
+     and might even be insufficient if z has not been computed with
+     sufficient precision itself (M below was calculated by Mathematica and also
+     with the series above with epsilon(1e-200)):
+    ; epsilon(1e-50)
+	0.00000000000000000001
+    ; display(50)
+	20
+    ; M=-0.9999999999999999999999997668356018402875796636464119050387
+    ; lambertw(-exp(-1)+1e-50,0)-M
+	-0.00000000000000000000000002678416515423276355643684
+    ; epsilon(1e-60)
+	0.0000000000000000000000000000000000000000000000000
+    ; A=-exp(-1)+1e-50
+    ; epsilon(1e-50)
+	0.00000000000000000000000000000000000000000000000000
+    ; lambertw(A,0)-M
+	-0.00000000000000000000000000000000000231185460220585
+    ; lambertw_series(A,epsilon(),0)-M
+	-0.00000000000000000000000000000000000132145133161626
+    ; epsilon(1e-100)
+	0.00000000000000000000000000000000000000000000000001
+    ; A=-exp(-1)+1e-50
+    ; epsilon(1e-65)
+	0.00000000000000000000000000000000000000000000000000
+    ; lambertw_series(A,epsilon(),0)-M
+	0.00000000000000000000000000000000000000000000000000
+    ; lambertw_series(-exp(-1)+1e-50,epsilon(),0)-M
+	-0.00000000000000000000000000000000000000002959444084
+    ; epsilon(1e-74)
+	0.00000000000000000000000000000000000000000000000000
+    ; lambertw_series(-exp(-1)+1e-50,epsilon(),0)-M
+	-0.00000000000000000000000000000000000000000000000006
+   */
+  closeness = abs(z-branchpoint);
+  if( closeness< 1){
+    if(closeness != 0)
+      eps = epsilon(epsilon()*( closeness));
+    else
+      eps = epsilon(epsilon()^2);
+  }
+  else
+    eps = epsilon(epsilon()*1e-2);
+
+
+  epslarge =epsilon();
+
+  places = highbit(1 + int(1/epslarge)) + 1;
+  w = ret;
+  for(k=0;k<100;k++){
+    ew = exp(w);
+    we = w*ew;
+    if(abs(we-z)<= 4*epslarge*abs(z))break;
+    w1e = (1+w)*ew;
+    wn = bround(w- ((we - z) / ( w1e - ( (w+2)*(we-z) )/(2*w+2)  ) ),places++) ;
+    if( abs(wn - w) <= epslarge*abs(wn)) break;
+    else  w = wn;
+  }
+
+  if(k==100){
+    epsilon(eps);
+    return newerror("lambertw: Halley iteration does not converge");
+  }
+  /* The Maxima coders added a check if the iteration converged to the correct
+     branch. This coder deems it superfluous. */
+
+  epsilon(eps);
+  return wn;
+}
+
+
+config("resource_debug", resource_debug_level),;
+if (config("resource_debug") & 3) {
+    print "lambertw(z,branch)";
+    print "lambertw_series(z,eps,branch,terms)";
+    print "lambertw_series_print()";
+}

cal/linear.cal

@@ -19,7 +19,7 @@
  *
  * @(#) $Revision: 30.2 $
  * @(#) $Id: linear.cal,v 30.2 2007/03/17 05:57:42 chongo Exp $
- * @(#) $Source: /usr/local/src/cmd/calc/cal/RCS/linear.cal,v $
+ * @(#) $Source: /usr/local/src/bin/calc/cal/RCS/linear.cal,v $
  *
  * Under source code control:	2005/12/12 06:41:50
  * File existed as early as:	2005

cal/lnseries.cal

@@ -0,0 +1,112 @@
+/*
+ * lnseries - special functions (e.g.: gamma, zeta, psi)
+ *
+ * Copyright (C) 2013 Christoph Zurnieden
+ *
+ * Calc is open software; you can redistribute it and/or modify it under
+ * the terms of the version 2.1 of the GNU Lesser General Public License
+ * as published by the Free Software Foundation.
+ *
+ * Calc is distributed in the hope that it will be useful, but WITHOUT
+ * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
+ * or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General
+ * Public License for more details.
+ *
+ * A copy of version 2.1 of the GNU Lesser General Public License is
+ * distributed with calc under the filename COPYING-LGPL.  You should have
+ * received a copy with calc; if not, write to Free Software Foundation, Inc.
+ * 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
+ *
+ * @(#) $Revision: 30.4 $
+ * @(#) $Id: lnseries.cal,v 30.4 2013/08/18 20:01:53 chongo Exp $
+ * @(#) $Source: /usr/local/src/bin/calc/cal/RCS/lnseries.cal,v $
+ *
+ * Under source code control:	2013/08/11 01:31:28
+ * File existed as early as:	2013
+ */
+
+
+/*
+ * hide internal function from resource debugging
+ */
+static resource_debug_level;
+resource_debug_level = config("resource_debug", 0);
+
+
+static __CZ__int_logs;
+static __CZ__int_logs_limit;
+static __CZ__int_logs_prec;
+
+
+define deletelnseries(){
+  free(__CZ__int_logs,__CZ__int_logs_limit,__CZ__int_logs_prec);
+}
+
+define lnfromseries(n){
+  if(    isnull(__CZ__int_logs)
+      || __CZ__int_logs_limit < n
+      || __CZ__int_logs_prec < log(1/epsilon())){
+
+    lnseries(n+1);
+  }
+  return __CZ__int_logs[n,0];
+}
+
+define lnseries(limit){
+  local k j eps ;
+  if(    isnull(__CZ__int_logs)
+      || __CZ__int_logs_limit < limit
+      || __CZ__int_logs_prec < log(1/epsilon())){
+    __CZ__int_logs = mat[limit+1,2];
+    __CZ__int_logs_limit = limit;
+    __CZ__int_logs_prec = log(1/epsilon());
+
+    /* probably still too much */
+    eps = epsilon(epsilon()*10^(-(5+log(limit))));
+    k =2;
+    while(1){
+      /* the prime itself, compute logarithm */
+      __CZ__int_logs[k,0] = ln(k);
+      __CZ__int_logs[k,1] = k;
+
+      for(j = 2*k;j<=limit;j+=k){
+       /* multiples of prime k, add logarithm of k computed earlier */
+       __CZ__int_logs[j,0] += __CZ__int_logs[k,0];
+       /* First hit, set counter to number */
+       if(__CZ__int_logs[j,1] ==0)
+         __CZ__int_logs[j,1]=j;
+       /* reduce counter by prime added */
+       __CZ__int_logs[j,1] //= __CZ__int_logs[k,1];
+      }
+
+      k++;
+      if(k>=limit) break;
+      /* Erastothenes-sieve: look for next prime. */
+      while(__CZ__int_logs[k,0]!=0){
+        k++;
+        if(k>=limit) break;
+      }
+    }
+    /* Second run to include the last factor */
+    for(k=1;k<=limit;k++){
+      if(__CZ__int_logs[k,1] != k){
+        __CZ__int_logs[k,0] +=__CZ__int_logs[ __CZ__int_logs[k,1],0];
+        __CZ__int_logs[k,1] = 0;
+      }
+    }
+
+    epsilon(eps);
+  }
+  return 1;
+}
+
+
+/*
+ * restore internal function from resource debugging
+ */
+config("resource_debug", resource_debug_level),;
+if (config("resource_debug") & 3) {
+    print "lnseries(limit)";
+    print "lnfromseries(n)";
+    print "deletelnseries()";
+}

cal/lucas.cal

@@ -19,7 +19,7 @@
  *
  * @(#) $Revision: 30.1 $
  * @(#) $Id: lucas.cal,v 30.1 2007/03/16 11:09:54 chongo Exp $
- * @(#) $Source: /usr/local/src/cmd/calc/cal/RCS/lucas.cal,v $
+ * @(#) $Source: /usr/local/src/bin/calc/cal/RCS/lucas.cal,v $
  *
  * Under source code control:	1990/05/03 16:49:51
  * File existed as early as:	1990

cal/lucas_chk.cal

@@ -19,7 +19,7 @@
  *
  * @(#) $Revision: 30.1 $
  * @(#) $Id: lucas_chk.cal,v 30.1 2007/03/16 11:09:54 chongo Exp $
- * @(#) $Source: /usr/local/src/cmd/calc/cal/RCS/lucas_chk.cal,v $
+ * @(#) $Source: /usr/local/src/bin/calc/cal/RCS/lucas_chk.cal,v $
  *
  * Under source code control:	1991/01/11 05:41:43
  * File existed as early as:	1991

cal/lucas_tbl.cal

@@ -19,7 +19,7 @@
  *
  * @(#) $Revision: 30.1 $
  * @(#) $Id: lucas_tbl.cal,v 30.1 2007/03/16 11:09:54 chongo Exp $
- * @(#) $Source: /usr/local/src/cmd/calc/cal/RCS/lucas_tbl.cal,v $
+ * @(#) $Source: /usr/local/src/bin/calc/cal/RCS/lucas_tbl.cal,v $
  *
  * Under source code control:	1991/01/26 02:43:43
  * File existed as early as:	1991

cal/mersenne.cal

@@ -21,7 +21,7 @@
  *
  * @(#) $Revision: 30.1 $
  * @(#) $Id: mersenne.cal,v 30.1 2007/03/16 11:09:54 chongo Exp $
- * @(#) $Source: /usr/local/src/cmd/calc/cal/RCS/mersenne.cal,v $
+ * @(#) $Source: /usr/local/src/bin/calc/cal/RCS/mersenne.cal,v $
  *
  * Under source code control:	1991/05/22 21:56:36
  * File existed as early as:	1991

cal/mfactor.cal

@@ -19,7 +19,7 @@
  *
  * @(#) $Revision: 30.1 $
  * @(#) $Id: mfactor.cal,v 30.1 2007/03/16 11:09:54 chongo Exp $
- * @(#) $Source: /usr/local/src/cmd/calc/cal/RCS/mfactor.cal,v $
+ * @(#) $Source: /usr/local/src/bin/calc/cal/RCS/mfactor.cal,v $
  *
  * Under source code control:	1996/07/06 06:09:40
  * File existed as early as:	1996

cal/mod.cal

@@ -19,7 +19,7 @@
  *
  * @(#) $Revision: 30.1 $
  * @(#) $Id: mod.cal,v 30.1 2007/03/16 11:09:54 chongo Exp $
- * @(#) $Source: /usr/local/src/cmd/calc/cal/RCS/mod.cal,v $
+ * @(#) $Source: /usr/local/src/bin/calc/cal/RCS/mod.cal,v $
  *
  * Under source code control:	1990/02/15 01:50:34
  * File existed as early as:	before 1990

cal/natnumset.cal

@@ -19,7 +19,7 @@
  *
  * @(#) $Revision: 30.1 $
  * @(#) $Id: natnumset.cal,v 30.1 2007/03/16 11:09:54 chongo Exp $
- * @(#) $Source: /usr/local/src/cmd/calc/cal/RCS/natnumset.cal,v $
+ * @(#) $Source: /usr/local/src/bin/calc/cal/RCS/natnumset.cal,v $
  *
  * Under source code control:	1997/09/07 23:53:51
  * File existed as early as:	1997

cal/pell.cal

@@ -19,7 +19,7 @@
  *
  * @(#) $Revision: 30.1 $
  * @(#) $Id: pell.cal,v 30.1 2007/03/16 11:09:54 chongo Exp $
- * @(#) $Source: /usr/local/src/cmd/calc/cal/RCS/pell.cal,v $
+ * @(#) $Source: /usr/local/src/bin/calc/cal/RCS/pell.cal,v $
  *
  * Under source code control:	1990/02/15 01:50:34
  * File existed as early as:	before 1990

cal/pi.cal

@@ -19,7 +19,7 @@
  *
  * @(#) $Revision: 30.1 $
  * @(#) $Id: pi.cal,v 30.1 2007/03/16 11:09:54 chongo Exp $
- * @(#) $Source: /usr/local/src/cmd/calc/cal/RCS/pi.cal,v $
+ * @(#) $Source: /usr/local/src/bin/calc/cal/RCS/pi.cal,v $
  *
  * Under source code control:	1991/05/22 21:56:37
  * File existed as early as:	1991

cal/pix.cal

@@ -19,7 +19,7 @@
  *
  * @(#) $Revision: 30.1 $
  * @(#) $Id: pix.cal,v 30.1 2007/03/16 11:09:54 chongo Exp $
- * @(#) $Source: /usr/local/src/cmd/calc/cal/RCS/pix.cal,v $
+ * @(#) $Source: /usr/local/src/bin/calc/cal/RCS/pix.cal,v $
  *
  * Under source code control:	1996/07/09 03:14:14
  * File existed as early as:	1996

cal/pollard.cal

@@ -19,7 +19,7 @@
  *
  * @(#) $Revision: 30.1 $
  * @(#) $Id: pollard.cal,v 30.1 2007/03/16 11:09:54 chongo Exp $
- * @(#) $Source: /usr/local/src/cmd/calc/cal/RCS/pollard.cal,v $
+ * @(#) $Source: /usr/local/src/bin/calc/cal/RCS/pollard.cal,v $
  *
  * Under source code control:	1991/05/22 21:56:37
  * File existed as early as:	1991

cal/poly.cal

@@ -19,7 +19,7 @@
  *
  * @(#) $Revision: 30.1 $
  * @(#) $Id: poly.cal,v 30.1 2007/03/16 11:09:54 chongo Exp $
- * @(#) $Source: /usr/local/src/cmd/calc/cal/RCS/poly.cal,v $
+ * @(#) $Source: /usr/local/src/bin/calc/cal/RCS/poly.cal,v $
  *
  * Under source code control:	1990/02/15 01:50:35
  * File existed as early as:	before 1990

cal/prompt.cal

@@ -19,7 +19,7 @@
  *
  * @(#) $Revision: 30.1 $
  * @(#) $Id: prompt.cal,v 30.1 2007/03/16 11:09:54 chongo Exp $
- * @(#) $Source: /usr/local/src/cmd/calc/cal/RCS/prompt.cal,v $
+ * @(#) $Source: /usr/local/src/bin/calc/cal/RCS/prompt.cal,v $
  *
  * Under source code control:	1995/12/18 04:43:25
  * File existed as early as:	1995

cal/psqrt.cal

@@ -19,7 +19,7 @@
  *
  * @(#) $Revision: 30.1 $
  * @(#) $Id: psqrt.cal,v 30.1 2007/03/16 11:09:54 chongo Exp $
- * @(#) $Source: /usr/local/src/cmd/calc/cal/RCS/psqrt.cal,v $
+ * @(#) $Source: /usr/local/src/bin/calc/cal/RCS/psqrt.cal,v $
  *
  * Under source code control:	1990/02/15 01:50:35
  * File existed as early as:	before 1990

cal/qtime.cal

@@ -22,7 +22,7 @@
  *
  * @(#) $Revision: 30.1 $
  * @(#) $Id: qtime.cal,v 30.1 2007/03/16 11:09:54 chongo Exp $
- * @(#) $Source: /usr/local/src/cmd/calc/cal/RCS/qtime.cal,v $
+ * @(#) $Source: /usr/local/src/bin/calc/cal/RCS/qtime.cal,v $
  *
  * Under source code control:	1999/10/13 04:10:33
  * File existed as early as:	1999

cal/quat.cal

@@ -17,9 +17,9 @@
  * received a copy with calc; if not, write to Free Software Foundation, Inc.
  * 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
  *
- * @(#) $Revision: 30.1 $
- * @(#) $Id: quat.cal,v 30.1 2007/03/16 11:09:54 chongo Exp $
- * @(#) $Source: /usr/local/src/cmd/calc/cal/RCS/quat.cal,v $
+ * @(#) $Revision: 30.2 $
+ * @(#) $Id: quat.cal,v 30.2 2013/08/11 08:41:38 chongo Exp $
+ * @(#) $Source: /usr/local/src/bin/calc/cal/RCS/quat.cal,v $
  *
  * Under source code control:	1990/02/15 01:50:35
  * File existed as early as:	before 1990
@@ -55,7 +55,8 @@ define quat(a,b,c,d)
 
 define quat_print(a)
 {
-	print "quat(" : a.s : ", " : a.v[0] : ", " : a.v[1] : ", " : a.v[2] : ")" :;
+	print "quat(" : a.s : ", " : a.v[0] : ", " :
+	      a.v[1] : ", " : a.v[2] : ")" :;
 }
 
 

cal/randbitrun.cal

@@ -19,7 +19,7 @@
  *
  * @(#) $Revision: 30.1 $
  * @(#) $Id: randbitrun.cal,v 30.1 2007/03/16 11:09:54 chongo Exp $
- * @(#) $Source: /usr/local/src/cmd/calc/cal/RCS/randbitrun.cal,v $
+ * @(#) $Source: /usr/local/src/bin/calc/cal/RCS/randbitrun.cal,v $
  *
  * Under source code control:	1995/02/13 03:43:11
  * File existed as early as:	1995

cal/randmprime.cal

@@ -19,7 +19,7 @@
  *
  * @(#) $Revision: 30.1 $
  * @(#) $Id: randmprime.cal,v 30.1 2007/03/16 11:09:54 chongo Exp $
- * @(#) $Source: /usr/local/src/cmd/calc/cal/RCS/randmprime.cal,v $
+ * @(#) $Source: /usr/local/src/bin/calc/cal/RCS/randmprime.cal,v $
  *
  * Under source code control:	1994/03/14 23:11:21
  * File existed as early as:	1994

cal/randombitrun.cal

@@ -19,7 +19,7 @@
  *
  * @(#) $Revision: 30.1 $
  * @(#) $Id: randombitrun.cal,v 30.1 2007/03/16 11:09:54 chongo Exp $
- * @(#) $Source: /usr/local/src/cmd/calc/cal/RCS/randombitrun.cal,v $
+ * @(#) $Source: /usr/local/src/bin/calc/cal/RCS/randombitrun.cal,v $
  *
  * Under source code control:	1995/02/13 03:43:11
  * File existed as early as:	1995

cal/randomrun.cal

@@ -19,7 +19,7 @@
  *
  * @(#) $Revision: 30.1 $
  * @(#) $Id: randomrun.cal,v 30.1 2007/03/16 11:09:54 chongo Exp $
- * @(#) $Source: /usr/local/src/cmd/calc/cal/RCS/randomrun.cal,v $
+ * @(#) $Source: /usr/local/src/bin/calc/cal/RCS/randomrun.cal,v $
  *
  * Under source code control:	1997/02/19 03:35:59
  * File existed as early as:	1997

cal/randrun.cal

@@ -19,7 +19,7 @@
  *
  * @(#) $Revision: 30.1 $
  * @(#) $Id: randrun.cal,v 30.1 2007/03/16 11:09:54 chongo Exp $
- * @(#) $Source: /usr/local/src/cmd/calc/cal/RCS/randrun.cal,v $
+ * @(#) $Source: /usr/local/src/bin/calc/cal/RCS/randrun.cal,v $
  *
  * Under source code control:	1995/02/12 20:00:06
  * File existed as early as:	1995

cal/regress.cal

@@ -17,9 +17,9 @@
  * received a copy with calc; if not, write to Free Software Foundation, Inc.
  * 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
  *
- * @(#) $Revision: 30.3 $
- * @(#) $Id: regress.cal,v 30.3 2007/07/11 22:58:26 chongo Exp $
- * @(#) $Source: /usr/local/src/cmd/calc/cal/RCS/regress.cal,v $
+ * @(#) $Revision: 30.12 $
+ * @(#) $Id: regress.cal,v 30.12 2013/09/02 02:32:55 chongo Exp $
+ * @(#) $Source: /usr/local/src/bin/calc/cal/RCS/regress.cal,v $
  *
  * Under source code control:	1990/02/15 01:50:36
  * File existed as early as:	before 1990
@@ -768,6 +768,8 @@ print '016: parsed test_bignums()';
 
 /*
  * Test many of the built-in functions.
+ *
+ * See test_functionss() starting at test 9000 for more built-in function tests.
  */
 define test_functions()
 {
@@ -1397,7 +1399,8 @@ define test_functions()
 	vrfy(quomod(10,-3,a,b,12) == 1, '1193: vrfy(quomod(10,-3,a,b,12) == 1');
 	vrfy(a == -3,			'1194: a == -3');
 	vrfy(b == 1,			'1195: b == 1');
-	vrfy(quomod(-10,-3,a,b,13) == 1,'1196: vrfy(quomod(-10,-3,a,b,13) == 1');
+	vrfy(quomod(-10,-3,a,b,13) == 1,
+	    '1196: vrfy(quomod(-10,-3,a,b,13) == 1');
 	vrfy(a == 4,			'1197: a == 4');
 	vrfy(b == 2,			'1198: b == 2');
 	vrfy(quomod(10,3,a,b,14) == 1,	'1199: vrfy(quomod(10,3,a,b,14) == 1');
@@ -1412,7 +1415,44 @@ define test_functions()
 	vrfy(isnum(systime()),		'1206: isnum(systime())');
 	vrfy(isnum(usertime()),		'1207: isnum(usertime())');
 
-	print '1208: Ending test_functions';
+	/* more jacobi tests */
+	vrfy(jacobi(987897,987) == 0,	'1208: jacobi(987897,987) == 0');
+	vrfy(jacobi(897,987) == 0,	'1209: jacobi(897,987) == 0');
+	vrfy(jacobi(987,897) == 0,	'1210: jacobi(987,897) == 0');
+	vrfy(jacobi(90,897) == 0,	'1211: jacobi(90,897) == 0');
+	vrfy(jacobi(45,897) == 0,	'1212: jacobi(45,897) == 0');
+	vrfy(jacobi(897,45) == 0,	'1213: jacobi(897,45) == 0');
+	vrfy(jacobi(42,45) == 0,	'1214: jacobi(42,45) == 0');
+	vrfy(jacobi(21,45) == 0,	'1215: jacobi(21,45) == 0');
+	vrfy(jacobi(45,21) == 0,	'1216: jacobi(45,21) == 0');
+	vrfy(jacobi(3,21) == 0,		'1217: jacobi(3,21) == 0');
+	vrfy(jacobi(0,21) == 0,		'1218: jacobi(0,21) == 0');
+	vrfy(jacobi(0,20003) == 0,	'1219: jacobi(0,20003) == 0');
+	vrfy(jacobi(1,20003) == 1,	'1220: jacobi(1,20003) == 1');
+	vrfy(jacobi(1236,20003) == 1,	'1221: jacobi(1236,20003) == 1');
+	vrfy(jacobi(618,20003) == -1,	'1222: jacobi(618,20003) == -1');
+	vrfy(jacobi(309,20003) == 1,	'1223: jacobi(309,20003) == 1');
+	vrfy(jacobi(227,309) == 1,	'1224: jacobi(227,309) == 1');
+	vrfy(jacobi(82,227) == 1,	'1225: jacobi(82,227) == 1');
+	vrfy(jacobi(41,227) == -1,	'1226: jacobi(41,227) == -1');
+	vrfy(jacobi(22,41) == -1,	'1227: jacobi(22,41) == -1');
+	vrfy(jacobi(11,41) == -1,	'1228: jacobi(11,41) == -1');
+	vrfy(jacobi(8,11) == -1,	'1229: jacobi(8,11) == -1');
+	vrfy(jacobi(4,11) == 1,		'1230: jacobi(4,11) == 1');
+	vrfy(jacobi(2,11) == -1,	'1231: jacobi(2,11) == -1');
+	vrfy(jacobi(1,11) == 1,		'1232: jacobi(1,11) == 1');
+	vrfy(jacobi(0,11) == 0,		'1233: jacobi(0,11) == 0');
+	vrfy(jacobi(0,0) == 0,		'1234: jacobi(0,0) == 0');
+	vrfy(jacobi(-1,0) == 0,		'1235: jacobi(-1,0) == 0');
+	vrfy(jacobi(-1,-1) == 0,	'1236: jacobi(-1,-1) == 0');
+	vrfy(jacobi(0,-1) == 0,		'1237: jacobi(0,-1) == 0');
+
+	/*
+	 * NOTE: Function tests are continued in test_functionss()
+	 *	 starting at test 9000.
+	 */
+
+	print '1293: Ending test_functions';
 }
 print '017: parsed test_functions()';
 
@@ -1429,14 +1469,14 @@ define _test_underscore()
 	local _a = 27;
 	local __a = 23209;
 
-	print "1290: Beginning _test_underscore";
+	print "1294: Beginning _test_underscore";
 
-	vrfy(_a == 27,			    '1291: _a == 27');
-	vrfy(_ == 49,			    '1292: _ == 49');
-	vrfy(__ == 63,			    '1293: __ == 63');
-	vrfy(__a == 23209,		    '1294: __a == 23209');
+	vrfy(_a == 27,			    '1295: _a == 27');
+	vrfy(_ == 49,			    '1296: _ == 49');
+	vrfy(__ == 63,			    '1297: __ == 63');
+	vrfy(__a == 23209,		    '1298: __a == 23209');
 
-	print "1295: Ending _test_underscore";
+	print "1299: Ending _test_underscore";
 }
 print '020: parsed _test_underscore';
 
@@ -1477,8 +1517,10 @@ define test_assoc()
 	vrfy(isnull(search(a,16)),	'1312: isnull(search(a,16))');
 	a["curds","whey"] = "spider";
 	print				'1313: a["curds","whey"] = "spider"';
-	vrfy(a["curds","whey"] == "spider", '1314: a["curds","whey"] == "spider"');
-	vrfy(a[[rsearch(a,"spider")]] == "spider", '1315: a[[rsearch(a,"spider")]] == "spider"');
+	vrfy(a["curds","whey"] == "spider",
+		'1314: a["curds","whey"] == "spider"');
+	vrfy(a[[rsearch(a,"spider")]] == "spider",
+		'1315: a[[rsearch(a,"spider")]] == "spider"');
 	b = a;
 	print				'1316: b = a';
 	vrfy(b[17] == 19,		'1317: b[17] == 19');
@@ -2706,11 +2748,24 @@ define test_2600()
 	vrfy(log(1e127) == 127,
 			    strcat(str(tnum++), ': log(1e127)) == 127'));
 	vrfy(round(log(17^47),10) == 57.8310993048,
-			    strcat(str(tnum++),
-			    	   ': round(log(17^47),10) == 57.8310993048'));
+	    strcat(str(tnum++),
+		   ': round(log(17^47),10) == 57.8310993048'));
 	vrfy(round(log(127),10) == 2.103803721,
-			    strcat(str(tnum++),
-			    	   ': round(log(127),10) == 2.103803721'));
+	    strcat(str(tnum++),
+		   ': round(log(127),10) == 2.103803721'));
+	vrfy(round(log(0.25,0.00001),5) == -0.60206,
+	    strcat(str(tnum++),
+		   ': round(log(0.25,0.00001),5) == -0.60206'));
+	vrfy(round(log(0.25,1e-10),10) == -0.6020599913,
+	    strcat(str(tnum++),
+		   ': round(log(0.25,1e-10),10) == -0.6020599913'));
+	vrfy(round( log(1.2+1.2i,1e-5),5) == 0.2297+0.34109i,
+	    strcat(str(tnum++),
+		   ': round(log(1.2+1.2i,1e-5),5) == 0.2297+0.34109i'));
+	vrfy(round( log(1.2+1.2i,1e-10),10) == 0.2296962439+0.3410940885i,
+	    strcat(str(tnum++),
+		   ': round(log(1.2+1.2i,1e-10),10) == ',
+		   '0.2296962439+0.3410940885i'));
         epsilon(i),;
 	print tnum++: ': epsilon(i),;';
 
@@ -4847,7 +4902,8 @@ define test_newsyn()
 		vrfy(s5500 == 55,		  '5510: s5500 == 45');
 		vrfy(i == 11,		  '5511: i == 11');
 	}
-	print "5512: { local i; for (s5500 = 0, i = 0; i < 10; i++) s5500 += i; ... }";
+	print "5512: { local i; for (s5500 = 0, i = 0; i < 10; i++) ":
+	      "s5500 += i; ... }";
 	vrfy(s5500 == 55,			  '5513: s5500 == 55');
 	vrfy(i == 11,			  '5514: i == 11');
 
@@ -4982,7 +5038,7 @@ vrfy(i == 9,		'151: i == 9');
 
 
 /*
- * test_commaeq - test changes to , and =
+ * test_commaeq - test changes to = and ,
  */
 obj xx5600 {} xx5600;
 print '152: obj xx5600 {} xx5600';
@@ -6461,7 +6517,8 @@ define test_blk()
 	/* A second named block */
 
 	B1 = blk("+++6700", 15, 10) = {1,2,3,4,5};
-	print				  '6746: B1 = blk("+++6700", 15 , 10) = {1,2,3,4,5};';
+	print
+	    '6746: B1 = blk("+++6700", 15 , 10) = {1,2,3,4,5};';
 	vrfy(size(B1) == 15,		  '6747: size(B1) == 15');
 	vrfy(sizeof(B1) == 20,		  '6748: sizeof(B1) == 20');
 	vrfy(test(B1) == 1,		  '6749: test(B1) == 1');
@@ -6826,7 +6883,8 @@ define test_sha1()
 
 
 	z = sha1(list(1,2,3), "curds and whey", 2^21701-1, pi(1e-100));
-	print '7210: z = sha1(list(1,2,3), "curds and whey", 2^21701-1, pi(1e-100));';
+	print '7210: z = sha1(list(1,2,3), "curds and whey",',
+	      '2^21701-1, pi(1e-100));';
 	vrfy(sha1(z) == 0x158cc87deeb9dd478ca14e3ab359205b0fb15b83,
 		'7211: sha1(z) == 0x158cc87deeb9dd478ca14e3ab359205b0fb15b83');
 
@@ -7869,6 +7927,137 @@ return test_exponentiation();
 /* 88xx: test exponentiation */
 
 
+/*
+ * calc resource functions by Christoph Zurnieden
+ */
+print;
+print '8900: Starting test of calc resource functions by Christoph Zurnieden';
+print '8901: read -once "test8900"';
+read -once "test8900";
+print '8902: about to run test8900(1,,8903)';
+testnum = test8900(1,,8903);
+print '8999: ecnt = 203;'
+ecnt = 203;
+/* 89xx: test calc resource functions by Christoph Zurnieden */
+
+
+/*
+ * Test more of the built-in functions.
+ *
+ * See test_functions() (test 700 - 1238) for other built-in function tests.
+ */
+define test_functions2()
+{
+	print '9001: Beginning test_functions2';
+
+	/* ctype function tests */
+	vrfy(isalnum("A") == 1,		'9002: isalnum("A") == 1');
+	vrfy(isalnum("a") == 1,		'9003: isalnum("a") == 1');
+	vrfy(isalnum("2") == 1,		'9004: isalnum("2") == 1');
+	vrfy(isalnum("\t") == 0,	'9005: isalnum("\\t") == 0');
+
+	vrfy(isalpha("A") == 1,		'9006: isalpha("A") == 1');
+	vrfy(isalpha("a") == 1,		'9007: isalpha("a") == 1');
+	vrfy(isalpha("2") == 0,		'9008: isalpha("2") == 0');
+	vrfy(isalpha("\t") == 0,	'9009: isalpha("\\t") == 0');
+
+	vrfy(iscntrl("A") == 0,		'9010: iscntrl("A") == 0');
+	vrfy(iscntrl("a") == 0,		'9011: iscntrl("a") == 0');
+	vrfy(iscntrl("2") == 0,		'9012: iscntrl("2") == 0');
+	vrfy(iscntrl("\t") == 1,	'9013: iscntrl("\\t") == 1');
+
+	vrfy(isdigit("A") == 0,		'9014: isdigit("A") == 0');
+	vrfy(isdigit("a") == 0,		'9015: isdigit("a") == 0');
+	vrfy(isdigit("2") == 1,		'9016: isdigit("2") == 1');
+	vrfy(isdigit("\t") == 0,	'9017: isdigit("\\t") == 0');
+
+	vrfy(isgraph("A") == 1,		'9018: isgraph("A") == 1');
+	vrfy(isgraph("a") == 1,		'9019: isgraph("a") == 1');
+	vrfy(isgraph("2") == 1,		'9020: isgraph("2") == 1');
+	vrfy(isgraph("\t") == 0,	'9021: isgraph("\\t") == 0');
+
+	vrfy(islower("A") == 0,		'9022: islower("A") == 0');
+	vrfy(islower("a") == 1,		'9023: islower("a") == 1');
+	vrfy(islower("1") == 0,		'9024: islower("1") == 0');
+
+	vrfy(isprint("A") == 1,		'9025: isprint("A") == 1');
+	vrfy(isprint("a") == 1,		'9026: isprint("a") == 1');
+	vrfy(isprint(" ") == 1,		'9027: isprint(" ") == 1');
+	vrfy(isprint("\t") == 0,	'9028: isprint("\\t") == 0');
+
+	vrfy(ispunct("A") == 0,		'9029: ispunct("A") == 0');
+	vrfy(ispunct("a") == 0,		'9030: ispunct("a") == 0');
+	vrfy(ispunct(" ") == 0,		'9031: ispunct(" ") == 0');
+	vrfy(ispunct("?") == 1,		'9032: ispunct("?") == 1');
+
+	vrfy(isspace("A") == 0,		'9033: isspace("A") == 0');
+	vrfy(isspace("Krik") == 0,	'9034: isspace("Krik") == 0');
+	vrfy(isspace(" ") == 1,		'9035: isspace(" ") == 1');
+	vrfy(isspace("?") == 0,		'9036: isspace("?") == 0');
+
+	vrfy(isupper("A") == 1,		'9037: isupper("A") == 1');
+	vrfy(isupper("a") == 0,		'9038: isupper("a") == 0');
+	vrfy(isupper("1") == 0,		'9039: isupper("1") == 0');
+
+	vrfy(isxdigit("A") == 1,	'9040: isxdigit("A") == 1');
+	vrfy(isxdigit("f") == 1,	'9041: isxdigit("f") == 1');
+	vrfy(isxdigit("2") == 1,	'9042: isxdigit("2") == 1');
+	vrfy(isxdigit("x") == 0,	'9043: isxdigit("x") == 0');
+
+	vrfy(strcasecmp("ab", "aBc") == -1,
+			     '9044: strcasecmp("ab", "aBc") == -1');
+	vrfy(strcasecmp("abc", "aBb") == 1,
+			     '9045: strcasecmp("abc", "aBb") == 1');
+	vrfy(strcasecmp("abc", "abc") == 0,
+			     '9046: strcasecmp("abc", "abc") == 0');
+	vrfy(strcasecmp("abc", "aBc") == 0,
+			     '9047: strcasecmp("abc", "aBc") == 0');
+	vrfy(strcasecmp("abc", "aBd") == -1,
+			     '9048: strcasecmp("abc", "aBd") == -1');
+	vrfy(strcasecmp("abc\0", "aBc") == 1,
+			     '9049: strcasecmp("a8c\\0", "aBc") == 1');
+	vrfy(strcasecmp("a\0b", "A\0c") == -1,
+			     '9050: strcasecmp("a\\0b", "A\\0c") == -1');
+
+	vrfy(strncasecmp("abc", "xyz", 0) == 0,
+			     '9051: strncasecmp("abc", "xyz", 0) == 0');
+	vrfy(strncasecmp("abc", "xyz", 1) == -1,
+			     '9052: strncasecmp("abc", "xyz", 1) == -1');
+	vrfy(strncasecmp("abc", "", 1) == 1,
+			     '9053: strncasecmp("abc", "", 1) == 1');
+	vrfy(strncasecmp("a", "b", 2) == -1,
+			     '9054: strncasecmp("a", "b", 2) == -1');
+	vrfy(strncasecmp("ab", "Ac", 2) == -1,
+			     '9055: strncasecmp("ab", "Ac", 2) == -1');
+	vrfy(strncasecmp("\0ac", "\0b", 2) == -1,
+			     '9056: strncasecmp("\\0ac", "\\0b", 2) == -1');
+	vrfy(strncasecmp("ab", "aBc", 2) == 0,
+			     '9057: strncasecmp("ab", "aBc", 2) == 0');
+	vrfy(strncasecmp("abc", "abd", 2) == 0,
+			     '9058: strncasecmp("abc", "abd", 2) == 0');
+
+	local s1 = " gnu lesser general public license";
+	print '9059: local s1 = " gnu lesser general public license";';
+
+	vrfy(strcmp(strtolower(" GNU Lesser General Public License"), s1) == 0,
+	    '9060: strcmp(strtolower(" GNU Lesser General Public License"),' +
+	    ' s1) == 0');
+
+	local s2 = " GNU LESSER GENERAL PUBLIC LICENSE";
+	print '9061: local s2 = " GNU LESSER GENERAL PUBLIC LICENSE";';
+
+	vrfy(strcmp(strtoupper(" GNU Lesser General Public License"), s2) == 0,
+	    '9062: strcmp(strtoupper(" GNU Lesser General Public License"),' +
+	    ' s2) == 0');
+
+	print '9063: Ending test_functions2';
+}
+print;
+print '9000: parsed test_functions2()';
+print;
+return test_functions2();
+
+
 /*
  * read various calc resource files
  *
@@ -7967,7 +8156,11 @@ print	'9838: skipping read -once hello.cal because it is an infinite loop';
 print	'9839: skipping read -once xx_print.cal because it is a printing demo';
 read -once sumtimes;
 print	'9840: read -once sumtimes';
-print	'9841: Ending read of selected calc resource files';
+read -once dms;
+print	'9841: read -once dms';
+read -once hms;
+print	'9842: read -once hms';
+print	'9843: Ending read of selected calc resource files';
 
 
 /*

cal/repeat.cal

@@ -19,7 +19,7 @@
  *
  * @(#) $Revision: 30.1 $
  * @(#) $Id: repeat.cal,v 30.1 2007/03/16 11:09:54 chongo Exp $
- * @(#) $Source: /usr/local/src/cmd/calc/cal/RCS/repeat.cal,v $
+ * @(#) $Source: /usr/local/src/bin/calc/cal/RCS/repeat.cal,v $
  *
  * Under source code control:	2003/01/05 00:00:01
  * File existed as early as:	2003

cal/screen.cal

@@ -17,7 +17,7 @@
  *
  * @(#) $Revision: 30.2 $
  * @(#) $Id: screen.cal,v 30.2 2007/03/16 11:09:54 chongo Exp $
- * @(#) $Source: /usr/local/src/cmd/calc/cal/RCS/screen.cal,v $
+ * @(#) $Source: /usr/local/src/bin/calc/cal/RCS/screen.cal,v $
  *
  * This file is not covered under version 2.1 of the GNU LGPL.
  *

cal/seedrandom.cal

@@ -19,7 +19,7 @@
  *
  * @(#) $Revision: 30.1 $
  * @(#) $Id: seedrandom.cal,v 30.1 2007/03/16 11:09:54 chongo Exp $
- * @(#) $Source: /usr/local/src/cmd/calc/cal/RCS/seedrandom.cal,v $
+ * @(#) $Source: /usr/local/src/bin/calc/cal/RCS/seedrandom.cal,v $
  *
  * Under source code control:	1996/01/01 08:21:00
  * File existed as early as:	1996

cal/set8700.cal

@@ -19,7 +19,7 @@
  *
  * @(#) $Revision: 30.1 $
  * @(#) $Id: set8700.cal,v 30.1 2007/03/16 11:09:54 chongo Exp $
- * @(#) $Source: /usr/local/src/cmd/calc/cal/RCS/set8700.cal,v $
+ * @(#) $Source: /usr/local/src/bin/calc/cal/RCS/set8700.cal,v $
  *
  * Under source code control:	2006/05/20 14:10:11
  * File existed as early as:	2006

cal/set8700.line

@@ -19,7 +19,7 @@
 ##
 ## @(#) $Revision: 30.1 $
 ## @(#) $Id: set8700.line,v 30.1 2007/03/16 11:09:54 chongo Exp $
-## @(#) $Source: /usr/local/src/cmd/calc/cal/RCS/set8700.line,v $
+## @(#) $Source: /usr/local/src/bin/calc/cal/RCS/set8700.line,v $
 ##
 ## Under source code control:	2006/05/20 14:10:11
 ## File existed as early as:	2006

cal/smallfactors.cal

@@ -0,0 +1,71 @@
+/*
+ * smallfactors - find the factors of a number < 2^32
+ *
+ * Copyright (C) 2013 Christoph Zurnieden
+ *
+ * Calc is open software; you can redistribute it and/or modify it under
+ * the terms of the version 2.1 of the GNU Lesser General Public License
+ * as published by the Free Software Foundation.
+ *
+ * Calc is distributed in the hope that it will be useful, but WITHOUT
+ * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
+ * or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General
+ * Public License for more details.
+ *
+ * A copy of version 2.1 of the GNU Lesser General Public License is
+ * distributed with calc under the filename COPYING-LGPL.  You should have
+ * received a copy with calc; if not, write to Free Software Foundation, Inc.
+ * 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
+ */
+
+
+static resource_debug_level;
+resource_debug_level = config("resource_debug", 0);
+
+
+define smallfactors(x0)
+{
+    local d q x flist tuple w;
+
+    if (x >= (2 ^ 32) - 1)
+	return newerror("smallfactors: number must be < 2^32 -1");
+
+    tuple = mat[2];
+    flist = list();
+    x = x0;
+    d = 2;
+    q = 0;
+    tuple[0] = d;
+    if (x < 2)
+	return 0;
+    do {
+	q = x // d;
+	    while (x == (q * d)) {
+	    tuple[0] = d;
+	    tuple[1]++;
+	    x = floor(q);
+	    q = x // d;
+	}
+	d = nextprime(d);
+	if (tuple[1] > 0)
+	    append(flist, tuple);
+	tuple = mat[2];
+    } while (d <= x);
+    return flist;
+}
+
+define printsmallfactors(flist)
+{
+    local k;
+    for (k = 0; k < size(flist); k++) {
+	print flist[k][0]:"^":flist[k][1];
+    }
+}
+
+
+config("resource_debug", resource_debug_level),;
+if (config("resource_debug") & 3) {
+    print "smallfactors(x0)";
+    print "printsmallfactors(flist)";
+
+}

cal/solve.cal

@@ -17,9 +17,9 @@
  * received a copy with calc; if not, write to Free Software Foundation, Inc.
  * 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
  *
- * @(#) $Revision: 30.1 $
- * @(#) $Id: solve.cal,v 30.1 2007/03/16 11:09:54 chongo Exp $
- * @(#) $Source: /usr/local/src/cmd/calc/cal/RCS/solve.cal,v $
+ * @(#) $Revision: 30.3 $
+ * @(#) $Id: solve.cal,v 30.3 2013/08/11 08:41:38 chongo Exp $
+ * @(#) $Source: /usr/local/src/bin/calc/cal/RCS/solve.cal,v $
  *
  * Under source code control:	1990/02/15 01:50:37
  * File existed as early as:	before 1990
@@ -47,12 +47,13 @@ define solve(low, high, epsilon)
 	if (abs(flow) < epsilon)
 		return low;
 	fhigh = f(high);
-	if (abs(flow) < epsilon)
+	if (abs(fhigh) < epsilon)
 		return high;
 	if (sgn(flow) == sgn(fhigh))
 		quit "Non-opposite signs";
 	while (1) {
-		mid = bround(high - fhigh * (high - low) / (fhigh - flow), places);
+		mid = bround(high - fhigh * (high - low) / (fhigh - flow),
+			     places);
 		if ((mid == low) || (mid == high))
 			places++;
 		fmid = f(mid);

cal/statistics.cal

@@ -0,0 +1,502 @@
+/*
+ * statistics - Some assorted statistics functions.
+ *
+ * Copyright (C) 2013 Christoph Zurnieden
+ *
+ * Calc is open software; you can redistribute it and/or modify it under
+ * the terms of the version 2.1 of the GNU Lesser General Public License
+ * as published by the Free Software Foundation.
+ *
+ * Calc is distributed in the hope that it will be useful, but WITHOUT
+ * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
+ * or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General
+ * Public License for more details.
+ *
+ * A copy of version 2.1 of the GNU Lesser General Public License is
+ * distributed with calc under the filename COPYING-LGPL.  You should have
+ * received a copy with calc; if not, write to Free Software Foundation, Inc.
+ * 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
+ *
+ * @(#) $Revision: 30.4 $
+ * @(#) $Id: statistics.cal,v 30.4 2013/08/18 20:01:53 chongo Exp $
+ * @(#) $Source: /usr/local/src/bin/calc/cal/RCS/statistics.cal,v $
+ *
+ * Under source code control:	2013/08/11 01:31:28
+ * File existed as early as:	2013
+ */
+
+
+static resource_debug_level;
+resource_debug_level = config("resource_debug", 0);
+
+
+/*
+ *  get dependencies
+ */
+read -once factorial2 brentsolve
+
+
+/*******************************************************************************
+ *
+ *
+ *                Continuous distributions
+ *
+ *
+ ******************************************************************************/
+
+/* regularized incomplete gamma function like in Octave, hence the name */
+define gammaincoctave(z,a){
+  local tmp;
+  tmp = gamma(z);
+  return (tmp-gammainc(a,z))/tmp;
+}
+
+/* Inverse incomplete beta function. Old and slow. */
+static __CZ__invbeta_a;
+static __CZ__invbeta_b;
+static __CZ__invbeta_x;
+define __CZ__invbeta(x){
+   return __CZ__invbeta_x-__CZ__ibetaas63(x,__CZ__invbeta_a,__CZ__invbeta_b);
+}
+
+define invbetainc_slow(x,a,b){
+  local flag ret eps;
+  /* place checks and balances here */
+  eps = epsilon();
+  if(.5 < x){
+    __CZ__invbeta_x = 1 - x;
+    __CZ__invbeta_a = b;
+    __CZ__invbeta_b = a;
+    flag = 1;
+  }
+  else{
+    __CZ__invbeta_x = x;
+    __CZ__invbeta_a = a;
+    __CZ__invbeta_b = b;
+    flag = 0;
+  }
+
+  ret =  brentsolve2(0,1,1);
+
+  if(flag == 1)
+    ret = 1-ret;
+  epsilon(eps);
+  return ret;
+}
+
+/* Inverse incomplete beta function. Still old but not as slow as the function
+   above. */
+/*
+  Purpose:
+
+    invbetainc computes inverse of the incomplete Beta function.
+
+  Licensing:
+
+    This code is distributed under the GNU LGPL license.
+
+  Modified:
+
+    10 August 2013
+
+  Author:
+
+    Original FORTRAN77 version by GW Cran, KJ Martin, GE Thomas.
+    C version by John Burkardt.
+    Calc version by Christoph Zurnieden
+
+  Reference:
+
+    GW Cran, KJ Martin, GE Thomas,
+    Remark AS R19 and Algorithm AS 109:
+    A Remark on Algorithms AS 63: The Incomplete Beta Integral
+    and AS 64: Inverse of the Incomplete Beta Integeral,
+    Applied Statistics,
+    Volume 26, Number 1, 1977, pages 111-114.
+
+  Parameters:
+
+    Input,  P, Q, the parameters of the incomplete
+    Beta function.
+
+    Input,  BETA, the logarithm of the value of
+    the complete Beta function.
+
+    Input,  ALPHA, the value of the incomplete Beta
+    function.  0 <= ALPHA <= 1.
+
+    Output,  the argument of the incomplete
+    Beta function which produces the value ALPHA.
+
+  Local Parameters:
+
+    Local, SAE, the most negative decimal exponent
+    which does not cause an underflow.
+*/
+define invbetainc(x,a,b){
+  return __CZ__invbetainc(a,b,lnbeta(a,b),x);
+}
+
+define __CZ__invbetainc(p,q,beta,alpha){
+  local a acu adj fpu g h iex indx pp prev qq r s sae sq t tx value;
+  local w xin y yprev places eps;
+
+  /* Dirty trick, don't try at home */
+  eps= epsilon(epsilon()^2);
+  sae = -((log(1/epsilon())/log(2))//2);
+  fpu = 10.0^sae;
+
+  places = highbit(1 + int(1/epsilon())) + 1;
+  value = alpha;
+  if( p <= 0.0 ){
+    epsilon(eps);
+    return newerror("invbeta: argument p <= 0");
+  }
+  if( q <= 0.0 ){
+    epsilon(eps);
+    return newerror("invbeta: argument q <= 0");
+  }
+
+  if( alpha < 0.0 || 1.0 < alpha ){
+    epsilon(eps);
+    return newerror("invbeta: argument alpha out of domain");
+  }
+  if( alpha == 0.0 ){
+    epsilon(eps);
+    return 0;
+  }
+  if( alpha == 1.0 ){
+    epsilon(eps);
+    return 1;
+  }
+  if ( 0.5 < alpha ){
+    a = 1.0 - alpha;
+    pp = q;
+    qq = p;
+    indx = 1;
+  }
+  else{
+    a = alpha;
+    pp = p;
+    qq = q;
+    indx = 0;
+  }
+  r = sqrt ( - ln ( a * a ) );
+
+  y = r-(2.30753+0.27061*r)/(1.0+(0.99229+0.04481*r)*r);
+
+  if ( 1.0 < pp && 1.0 < qq ){
+    r = ( y * y - 3.0 ) / 6.0;
+    s = 1.0 / ( pp + pp - 1.0 );
+    t = 1.0 / ( qq + qq - 1.0 );
+    h = 2.0 / ( s + t );
+    w = y*sqrt(h+r)/h-(t-s)*(r+5.0/6.0-2.0/(3.0*h));
+    value = pp / ( pp + qq * exp ( w + w ) );
+  }
+  else{
+    r = qq + qq;
+    t = 1.0 / ( 9.0 * qq );
+    t = r *  ( 1.0 - t + y * sqrt ( t )^ 3 );
+
+    if ( t <= 0.0 ){
+      value = 1.0 - exp ( ( ln ( ( 1.0 - a ) * qq ) + beta ) / qq );
+    }
+    else{
+      t = ( 4.0 * pp + r - 2.0 ) / t;
+
+      if ( t <= 1.0 ) {
+        value = exp ( ( ln ( a * pp ) + beta ) / pp );
+      }
+      else{
+        value = 1.0 - 2.0 / ( t + 1.0 );
+      }
+    }
+  }
+  r = 1.0 - pp;
+  t = 1.0 - qq;
+  yprev = 0.0;
+  sq = 1.0;
+  prev = 1.0;
+
+  if ( value < 0.0001 )
+    value = 0.0001;
+
+  if ( 0.9999 < value )
+    value = 0.9999;
+
+  acu = 10^sae;
+
+  for ( ; ; ){
+    y = bround(__CZ__ibetaas63( value, pp, qq, beta),places);
+    xin = value;
+    y = bround(exp(ln(y-a)+(beta+r*ln(xin)+t*ln(1.0- xin ) )),places);
+
+    if ( y * yprev <= 0.0 ) {
+      prev = max ( sq, fpu );
+    }
+
+    g = 1.0;
+
+    for ( ; ; ){
+      for ( ; ; ){
+        adj = g * y;
+        sq = adj * adj;
+        if ( sq < prev ){
+          tx = value - adj;
+          if ( 0.0 <= tx && tx <= 1.0 )  break;
+        }
+        g = g / 3.0;
+      }
+      if ( prev <= acu ){
+        if ( indx )
+          value = 1.0 - value;
+        epsilon(eps);
+        return value;
+      }
+      if ( y * y <= acu ){
+        if ( indx )
+          value = 1.0 - value;
+        epsilon(eps);
+        return value;
+      }
+      if ( tx != 0.0 && tx != 1.0 )
+        break;
+      g = g / 3.0;
+    }
+    if ( tx == value )  break;
+    value = tx;
+    yprev = y;
+  }
+  if ( indx )
+    value = 1.0 - value;
+
+  epsilon(eps);
+  return value;
+}
+
+/*******************************************************************************
+ *
+ *
+ *                Beta distribution
+ *
+ *
+ ******************************************************************************/
+
+define betapdf(x,a,b){
+  if(x<0 || x>1) return newerror("betapdf: parameter x out of domain");
+  if(a<=0) return newerror("betapdf: parameter a out of domain");
+  if(b<=0) return newerror("betapdf: parameter b out of domain");
+
+  return 1/beta(a,b) *x^(a-1)*(1-x)^(b-1);
+}
+
+define betacdf(x,a,b){
+  if(x<0 || x>1) return newerror("betacdf: parameter x out of domain");
+  if(a<=0) return newerror("betacdf: parameter a out of domain");
+  if(b<=0) return newerror("betacdf: parameter b out of domain");
+
+  return betainc(x,a,b);
+}
+
+define betacdfinv(x,a,b){
+  return invbetainc(x,a,b);
+}
+
+define betamedian(a,b){
+  local t106 t104 t103 t105 approx ret;
+  if(a == b) return 1/2;
+  if(a == 1 && b > 0) return 1-(1/2)^(1/b);
+  if(a > 0 && b == 1) return   (1/2)^(1/a);
+  if(a == 3 && b == 2){
+    /* Yes, the author is not ashamed to ask Maxima for the exact solution
+       of a quartic equation. */
+    t103 = ( (2^(3/2))/27 +4/27  )^(1/3);
+    t104 = sqrt( ( 9*t103^2 + 4*t103 + 2  )/(t103)  )/3;
+    t105 = -t103-2/(9*t103) +8/9;
+    t106 = sqrt( (27*t104*t105+16)/(t104)  )/(2*3^(3/2));
+    return -t106+t104/2+1/3;
+  }
+  if(a == 2 && b == 3){
+    t103 = ( (2^(3/2))/27 +4/27  )^(1/3);
+    t104 = sqrt( ( 9*t103^2 + 4*t103 + 2  )/(t103)  )/3;
+    t105 = -t103-2/(9*t103) +8/9;
+    t106 = sqrt( (27*t104*t105+16)/(t104)  )/(2*3^(3/2));
+    return 1-(-t106+t104/2+1/3);
+  }
+  return invbetainc(1/2,a,b);
+}
+
+define betamode(a,b){
+  if(a + b == 2) return newerror("betamod: a + b = 2 = division by zero");
+  return (a-1)/(a+b-2);
+}
+
+define betavariance(a,b){
+  return (a*b)/( (a+b)^2*(a+b+1) );
+}
+
+define betalnvariance(a,b){
+  return polygamma(1,a)-polygamma(a+b);
+}
+
+define betaskewness(a,b){
+  return (2*(b-a)*sqrt(a+b+1))/( (a+b+1)*sqrt(a*b) );
+}
+
+define betakurtosis(a,b){
+  local num denom;
+
+  num = 6*( (a-b)^2*(a+b+1)-a*b*(a+b+2));
+  denom = a*b*(a+b+2)*(a+b+3);
+  return num/denom;
+}
+
+define betaentropy(a,b){
+  return lnbeta(a,b)-(a-1)*psi(a)-(b-1)*psi(b)+(a+b+1)*psi(a+b);
+
+}
+
+/*******************************************************************************
+ *
+ *
+ *                Normal (Gaussian) distribution
+ *
+ *
+ ******************************************************************************/
+
+
+define normalpdf(x,mu,sigma){
+  return 1/(sqrt(2*pi()*sigma^2))*exp( ( (x-mu)^2 )/( 2*sigma^2 ) );
+}
+
+define normalcdf(x,mu,sigma){
+  return 1/2*(1+erf( ( x-mu )/( sqrt(2*sigma^2) )  )  );
+}
+
+define probit(p){
+  if(p<0 || p > 1) return newerror("probit: p out of domain 0<=p<=1");
+  return sqrt(2)*ervinv(2*p-1);
+}
+
+define normalcdfinv(p,mu,sigma){
+  if(p<0 || p > 1) return newerror("normalcdfinv: p out of domain 0<=p<=1");
+  return mu+ sigma*probit(p);
+}
+
+define normalmean(mu,sigma){return mu;}
+
+define normalmedian(mu,sigma){return mu;}
+
+define normalmode(mu,sigma){return mu;}
+
+define normalvariance(mu,sigma){return sigma^2;}
+
+define normalskewness(mu,sigma){return 0;}
+
+define normalkurtosis(mu,sigma){return 0;}
+
+define normalentropy(mu,sigma){
+  return 1/3*ln( 2*pi()*exp(1)*sigma^2 );
+}
+
+/* moment generating f. */
+define normalmgf(mu,sigma,t){
+  return exp(mu*t+1/2*sigma^2*t^2);
+}
+
+/* characteristic f. */
+define normalcf(mu,sigma,t){
+  return exp(mu*t-1/2*sigma^2*t^2);
+}
+
+
+/*******************************************************************************
+ *
+ *
+ *                Chi-squared distribution
+ *
+ *
+ ******************************************************************************/
+
+define chisquaredpdf(x,k){
+  if(!isint(k) || k<0) return newerror("chisquaredpdf: k not in N");
+  if(im(x) || x<0) return newerror("chisquaredpdf: x not in +R");
+  /* The gamma function does not check for half integers, do it here? */
+  return 1/(2^(k/2)*gamma(k/2))*x^((k/2)-1)*exp(-x/2);
+}
+
+define chisquaredpcdf(x,k){
+  if(!isint(k) || k<0) return newerror("chisquaredcdf: k not in N");
+  if(im(x) || x<0) return newerror("chisquaredcdf: x not in +R");
+
+  return 1/(gamma(k/2))*gammainc(k/2,x/2);
+}
+
+define chisquaredmean(x,k){return k;}
+
+define chisquaredmedian(x,k){
+  /* TODO: implement a FAST inverse incomplete gamma-{q,p} function */
+  return k*(1-2/(9*k))^3;
+}
+
+define chisquaredmode(x,k){return max(k-2,0);}
+define chisquaredvariance(x,k){return 2*k;}
+define chisquaredskewness(x,k){return sqrt(8/k);}
+define chisquaredkurtosis(x,k){return 12/k;}
+define chisquaredentropy(x,k){
+  return k/2+ln(2*gamma(k/2)) + (1-k/2)*psi(k/2);
+}
+
+define chisquaredmfg(k,t){
+  if(t>=1/2)return newerror("chisquaredmfg: t >= 1/2");
+  return (1-2*t)^(k/2);
+}
+
+define chisquaredcf(k,t){
+  return (1-2*1i*t)^(k/2);
+}
+
+
+/*
+ * restore internal function from resource debugging
+ */
+config("resource_debug", resource_debug_level),;
+if (config("resource_debug") & 3) {
+    print "gammaincoctave(z,a)";
+    print "invbetainc(x,a,b)";
+    print "betapdf(x,a,b)";
+    print "betacdf(x,a,b)";
+    print "betacdfinv(x,a,b)";
+    print "betamedian(a,b)";
+    print "betamode(a,b)";
+    print "betavariance(a,b)";
+    print "betalnvariance(a,b)";
+    print "betaskewness(a,b)";
+    print "betakurtosis(a,b)";
+    print "betaentropy(a,b)";
+    print "normalpdf(x,mu,sigma)";
+    print "normalcdf(x,mu,sigma)";
+    print "probit(p)";
+    print "normalcdfinv(p,mu,sigma)";
+    print "normalmean(mu,sigma)";
+    print "normalmedian(mu,sigma)";
+    print "normalmode(mu,sigma)";
+    print "normalvariance(mu,sigma)";
+    print "normalskewness(mu,sigma)";
+    print "normalkurtosis(mu,sigma)";
+    print "normalentropy(mu,sigma)";
+    print "normalmgf(mu,sigma,t)";
+    print "normalcf(mu,sigma,t)";
+    print "chisquaredpdf(x,k)";
+    print "chisquaredpcdf(x,k)";
+    print "chisquaredmean(x,k)";
+    print "chisquaredmedian(x,k)";
+    print "chisquaredmode(x,k)";
+    print "chisquaredvariance(x,k)";
+    print "chisquaredskewness(x,k)";
+    print "chisquaredkurtosis(x,k)";
+    print "chisquaredentropy(x,k)";
+    print "chisquaredmfg(k,t)";
+    print "chisquaredcf(k,t)";
+}
+

cal/strings.cal

@@ -0,0 +1,41 @@
+/*
+ * strings - implementation of some of the macros in ctype.h
+ *
+ * Copyright (C) 2013 Christoph Zurnieden
+ *
+ * Calc is open software; you can redistribute it and/or modify it under
+ * the terms of the version 2.1 of the GNU Lesser General Public License
+ * as published by the Free Software Foundation.
+ *
+ * Calc is distributed in the hope that it will be useful, but WITHOUT
+ * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
+ * or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General
+ * Public License for more details.
+ *
+ * A copy of version 2.1 of the GNU Lesser General Public License is
+ * distributed with calc under the filename COPYING-LGPL.  You should have
+ * received a copy with calc; if not, write to Free Software Foundation, Inc.
+ * 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
+ */
+
+
+static resource_debug_level;
+resource_debug_level = config("resource_debug", 0);
+
+define isascii(c){
+  c = ord(c);
+  return (c >= 0 && c< 128);
+}
+
+define isblank(c){
+  c = ord(c);
+  return ( c == 32 || c == 9 );
+}
+
+
+config("resource_debug", resource_debug_level),;
+if (config("resource_debug") & 3) {
+    print "isascii(c)";
+    print "isblank(c)";
+}
+

cal/sumsq.cal

@@ -19,7 +19,7 @@
  *
  * @(#) $Revision: 30.1 $
  * @(#) $Id: sumsq.cal,v 30.1 2007/03/16 11:09:54 chongo Exp $
- * @(#) $Source: /usr/local/src/cmd/calc/cal/RCS/sumsq.cal,v $
+ * @(#) $Source: /usr/local/src/bin/calc/cal/RCS/sumsq.cal,v $
  *
  * Under source code control:	1990/02/15 01:50:37
  * File existed as early as:	before 1990

cal/sumtimes.cal

@@ -19,7 +19,7 @@
  *
  * @(#) $Revision: 30.1 $
  * @(#) $Id: sumtimes.cal,v 30.1 2007/03/16 11:09:54 chongo Exp $
- * @(#) $Source: /usr/local/src/cmd/calc/cal/RCS/sumtimes.cal,v $
+ * @(#) $Source: /usr/local/src/bin/calc/cal/RCS/sumtimes.cal,v $
  *
  * Under source code control:	2006/06/22 17:29
  * File existed as early as:	2006

cal/surd.cal

@@ -19,7 +19,7 @@
  *
  * @(#) $Revision: 30.1 $
  * @(#) $Id: surd.cal,v 30.1 2007/03/16 11:09:54 chongo Exp $
- * @(#) $Source: /usr/local/src/cmd/calc/cal/RCS/surd.cal,v $
+ * @(#) $Source: /usr/local/src/bin/calc/cal/RCS/surd.cal,v $
  *
  * Under source code control:	1990/02/15 01:50:38
  * File existed as early as:	before 1990

cal/test1700.cal

@@ -19,7 +19,7 @@
  *
  * @(#) $Revision: 30.1 $
  * @(#) $Id: test1700.cal,v 30.1 2007/03/16 11:09:54 chongo Exp $
- * @(#) $Source: /usr/local/src/cmd/calc/cal/RCS/test1700.cal,v $
+ * @(#) $Source: /usr/local/src/bin/calc/cal/RCS/test1700.cal,v $
  *
  * Under source code control:	1994/03/14 23:12:51
  * File existed as early as:	1994

cal/test2300.cal

@@ -19,7 +19,7 @@
  *
  * @(#) $Revision: 30.1 $
  * @(#) $Id: test2300.cal,v 30.1 2007/03/16 11:09:54 chongo Exp $
- * @(#) $Source: /usr/local/src/cmd/calc/cal/RCS/test2300.cal,v $
+ * @(#) $Source: /usr/local/src/bin/calc/cal/RCS/test2300.cal,v $
  *
  * Under source code control:	1995/07/09 06:12:13
  * File existed as early as:	1995

cal/test2600.cal

@@ -19,9 +19,9 @@
  * received a copy with calc; if not, write to Free Software Foundation, Inc.
  * 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
  *
- * @(#) $Revision: 30.2 $
- * @(#) $Id: test2600.cal,v 30.2 2007/07/11 22:57:23 chongo Exp $
- * @(#) $Source: /usr/local/src/cmd/calc/cal/RCS/test2600.cal,v $
+ * @(#) $Revision: 30.3 $
+ * @(#) $Id: test2600.cal,v 30.3 2013/08/11 08:41:38 chongo Exp $
+ * @(#) $Source: /usr/local/src/bin/calc/cal/RCS/test2600.cal,v $
  *
  * Under source code control:	1995/10/13 00:13:14
  * File existed as early as:	1995
@@ -91,7 +91,8 @@ define testismult(str, n, verbose)
 		if (!ismult(c,a)) {
 			m++;
 			if (verbose > 1) {
-				printf("*** Failure with:\na = %d\nb = %d\n", a,b);
+				printf("*** Failure with:\na = %d\nb = %d\n",
+				       a,b);
 			}
 		}
 	}
@@ -133,7 +134,8 @@ define testsqrt(str, n, eps, verbose)
 		if (abs(c) > 1) {
 			m++;
 			if (verbose > 1) {
-				printf("*** Failure with:\na = %d\neps = %d\n", a,eps);
+				printf("*** Failure with:\na = %d\neps = %d\n",
+				       a,eps);
 			}
 		}
 	}
@@ -178,7 +180,8 @@ define testexp(str, n, eps, verbose)
 		if (abs(c) > 0.02) {
 			m++;
 			if (verbose > 1) {
-				printf("*** Failure with:\na = %d\neps = %d\n", a,eps);
+				printf("*** Failure with:\na = %d\neps = %d\n",
+				       a,eps);
 			}
 		}
 	}
@@ -235,7 +238,8 @@ define testln(str, n, eps, verbose)
 		if (abs(c) > 0.5) {
 			m++;
 			if (verbose > 1) {
-				printf("*** Failure with:\na = %d\neps = %d\n", a,eps);
+				printf("*** Failure with:\na = %d\neps = %d\n",
+				       a,eps);
 			}
 		}
 	}

cal/test2700.cal

@@ -19,9 +19,9 @@
  * received a copy with calc; if not, write to Free Software Foundation, Inc.
  * 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
  *
- * @(#) $Revision: 30.1 $
- * @(#) $Id: test2700.cal,v 30.1 2007/03/16 11:09:54 chongo Exp $
- * @(#) $Source: /usr/local/src/cmd/calc/cal/RCS/test2700.cal,v $
+ * @(#) $Revision: 30.2 $
+ * @(#) $Id: test2700.cal,v 30.2 2013/08/11 08:41:38 chongo Exp $
+ * @(#) $Source: /usr/local/src/bin/calc/cal/RCS/test2700.cal,v $
  *
  * Under source code control:	1995/11/01 22:52:25
  * File existed as early as:	1995
@@ -127,7 +127,8 @@ define testcsqrt(str, n, verbose)
 			if (p) {
 			if (verbose > 0)
 				printf(
-				 "*** Type %d failure for x = %r, y = %r, z = %d\n",
+				 "*** Type %d failure for x = %r, "
+				 "y = %r, z = %d\n",
 				    p, x, y, z);
 				m++;
 			}

cal/test3100.cal

@@ -21,7 +21,7 @@
  *
  * @(#) $Revision: 30.1 $
  * @(#) $Id: test3100.cal,v 30.1 2007/03/16 11:09:54 chongo Exp $
- * @(#) $Source: /usr/local/src/cmd/calc/cal/RCS/test3100.cal,v $
+ * @(#) $Source: /usr/local/src/bin/calc/cal/RCS/test3100.cal,v $
  *
  * Under source code control:	1995/11/28 11:56:57
  * File existed as early as:	1995

cal/test3300.cal

@@ -21,7 +21,7 @@
  *
  * @(#) $Revision: 30.1 $
  * @(#) $Id: test3300.cal,v 30.1 2007/03/16 11:09:54 chongo Exp $
- * @(#) $Source: /usr/local/src/cmd/calc/cal/RCS/test3300.cal,v $
+ * @(#) $Source: /usr/local/src/bin/calc/cal/RCS/test3300.cal,v $
  *
  * Under source code control:	1995/12/02 04:27:41
  * File existed as early as:	1995

cal/test3400.cal

@@ -21,7 +21,7 @@
  *
  * @(#) $Revision: 30.1 $
  * @(#) $Id: test3400.cal,v 30.1 2007/03/16 11:09:54 chongo Exp $
- * @(#) $Source: /usr/local/src/cmd/calc/cal/RCS/test3400.cal,v $
+ * @(#) $Source: /usr/local/src/bin/calc/cal/RCS/test3400.cal,v $
  *
  * Under source code control:	1995/12/02 05:20:11
  * File existed as early as:	1995

cal/test3500.cal

@@ -21,7 +21,7 @@
  *
  * @(#) $Revision: 30.1 $
  * @(#) $Id: test3500.cal,v 30.1 2007/03/16 11:09:54 chongo Exp $
- * @(#) $Source: /usr/local/src/cmd/calc/cal/RCS/test3500.cal,v $
+ * @(#) $Source: /usr/local/src/bin/calc/cal/RCS/test3500.cal,v $
  *
  * Under source code control:	1995/12/18 22:50:46
  * File existed as early as:	1995

cal/test4000.cal

@@ -19,9 +19,9 @@
  * received a copy with calc; if not, write to Free Software Foundation, Inc.
  * 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
  *
- * @(#) $Revision: 30.1 $
- * @(#) $Id: test4000.cal,v 30.1 2007/03/16 11:09:54 chongo Exp $
- * @(#) $Source: /usr/local/src/cmd/calc/cal/RCS/test4000.cal,v $
+ * @(#) $Revision: 30.2 $
+ * @(#) $Id: test4000.cal,v 30.2 2013/08/11 08:41:38 chongo Exp $
+ * @(#) $Source: /usr/local/src/bin/calc/cal/RCS/test4000.cal,v $
  *
  * Under source code control:	1996/03/13 02:38:45
  * File existed as early as:	1996
@@ -199,7 +199,8 @@ define ctimes(str, N, n, count, skip, verbose)
 		p = ptest(A[i], count, skip);
 		if (p) {
 			if (verbose > 0) {
-				printf("*** Error, what should be rare has occurred for x = %d \n", A[i]);
+				printf("*** Error, what should be rare "
+				       "has occurred for x = %d \n", A[i]);
 				m++;
 			}
 		}
@@ -306,7 +307,8 @@ define ntimes(str, N, n, count, skip, residue, modulus, verbose)
 	}
 	tprev = round(usertime() - t, 4);
 	if (verbose > 0) {
-		printf("%d evaluations, nextcand: %d, prevcand: %d\n", n, tnext, tprev);
+		printf("%d evaluations, nextcand: %d, "
+		       "prevcand: %d\n", n, tnext, tprev);
 	}
 }
 

cal/test4100.cal

@@ -21,7 +21,7 @@
  *
  * @(#) $Revision: 30.1 $
  * @(#) $Id: test4100.cal,v 30.1 2007/03/16 11:09:54 chongo Exp $
- * @(#) $Source: /usr/local/src/cmd/calc/cal/RCS/test4100.cal,v $
+ * @(#) $Source: /usr/local/src/bin/calc/cal/RCS/test4100.cal,v $
  *
  * Under source code control:	1996/03/13 03:53:22
  * File existed as early as:	1996

cal/test4600.cal

@@ -21,7 +21,7 @@
  *
  * @(#) $Revision: 30.1 $
  * @(#) $Id: test4600.cal,v 30.1 2007/03/16 11:09:54 chongo Exp $
- * @(#) $Source: /usr/local/src/cmd/calc/cal/RCS/test4600.cal,v $
+ * @(#) $Source: /usr/local/src/bin/calc/cal/RCS/test4600.cal,v $
  *
  * Under source code control:	1996/07/02 20:04:40
  * File existed as early as:	1996

cal/test5100.cal

@@ -21,7 +21,7 @@
  *
  * @(#) $Revision: 30.1 $
  * @(#) $Id: test5100.cal,v 30.1 2007/03/16 11:09:54 chongo Exp $
- * @(#) $Source: /usr/local/src/cmd/calc/cal/RCS/test5100.cal,v $
+ * @(#) $Source: /usr/local/src/bin/calc/cal/RCS/test5100.cal,v $
  *
  * Under source code control:	1996/12/02 23:57:10
  * File existed as early as:	1996

cal/test5200.cal

@@ -21,7 +21,7 @@
  *
  * @(#) $Revision: 30.1 $
  * @(#) $Id: test5200.cal,v 30.1 2007/03/16 11:09:54 chongo Exp $
- * @(#) $Source: /usr/local/src/cmd/calc/cal/RCS/test5200.cal,v $
+ * @(#) $Source: /usr/local/src/bin/calc/cal/RCS/test5200.cal,v $
  *
  * Under source code control:	1997/02/07 02:48:10
  * File existed as early as:	1997

cal/test8400.cal

@@ -19,7 +19,7 @@
  *
  * @(#) $Revision: 30.1 $
  * @(#) $Id: test8400.cal,v 30.1 2007/03/16 11:09:54 chongo Exp $
- * @(#) $Source: /usr/local/src/cmd/calc/cal/RCS/test8400.cal,v $
+ * @(#) $Source: /usr/local/src/bin/calc/cal/RCS/test8400.cal,v $
  *
  * Under source code control:	1999/10/31 01:00:03
  * File existed as early as:	1999

cal/test8500.cal

@@ -19,9 +19,9 @@
  * received a copy with calc; if not, write to Free Software Foundation, Inc.
  * 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
  *
- * @(#) $Revision: 30.1 $
- * @(#) $Id: test8500.cal,v 30.1 2007/03/16 11:09:54 chongo Exp $
- * @(#) $Source: /usr/local/src/cmd/calc/cal/RCS/test8500.cal,v $
+ * @(#) $Revision: 30.2 $
+ * @(#) $Id: test8500.cal,v 30.2 2013/08/11 08:41:38 chongo Exp $
+ * @(#) $Source: /usr/local/src/bin/calc/cal/RCS/test8500.cal,v $
  *
  * Under source code control:	1999/11/12 20:59:59
  * File existed as early as:	1999
@@ -134,8 +134,8 @@ define onetest_8500(a,b,rnd) {
  * The rounding parameter is randomly chosen.
  *
  * After a run of divmod_8500 the a, b, rnd values which gave failure are
- * stored in the list L_8500.  L_8500[0], L_8500[1], L_8500[2] are a, b, rnd for the first
- * test, etc.
+ * stored in the list L_8500.  L_8500[0], L_8500[1], L_8500[2] are a, b,
+ * rnd for the first* test, etc.
  */
 define divmod_8500(N = 10, M1 = 2^128, M2 = 2^64, testnum = 0)
 {

cal/test8600.cal

@@ -21,7 +21,7 @@
  *
  * @(#) $Revision: 30.1 $
  * @(#) $Id: test8600.cal,v 30.1 2007/03/16 11:09:54 chongo Exp $
- * @(#) $Source: /usr/local/src/cmd/calc/cal/RCS/test8600.cal,v $
+ * @(#) $Source: /usr/local/src/bin/calc/cal/RCS/test8600.cal,v $
  *
  * Under source code control:	2000/12/04 19:57:02
  * File existed as early as:	2000

cal/toomcook.cal

@@ -0,0 +1,362 @@
+/*
+ * toomcook - implementation of Toom-Cook(3,4) multiplication algorithm
+ *
+ * Copyright (C) 2013 Christoph Zurnieden
+ *
+ * Calc is open software; you can redistribute it and/or modify it under
+ * the terms of the version 2.1 of the GNU Lesser General Public License
+ * as published by the Free Software Foundation.
+ *
+ * Calc is distributed in the hope that it will be useful, but WITHOUT
+ * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
+ * or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General
+ * Public License for more details.
+ *
+ * A copy of version 2.1 of the GNU Lesser General Public License is
+ * distributed with calc under the filename COPYING-LGPL.  You should have
+ * received a copy with calc; if not, write to Free Software Foundation, Inc.
+ * 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
+ *
+ * @(#) $Revision: 30.4 $
+ * @(#) $Id: toomcook.cal,v 30.4 2013/08/18 20:01:53 chongo Exp $
+ * @(#) $Source: /usr/local/src/bin/calc/cal/RCS/toomcook.cal,v $
+ *
+ * Under source code control:	2013/08/11 01:31:28
+ * File existed as early as:	2013
+ */
+
+/*
+ * hide internal function from resource debugging
+ */
+static resource_debug_level;
+resource_debug_level = config("resource_debug", 0);
+
+
+/* */
+define toomcook3(a,b){
+  local alen blen  a0 a1 a2  b0 b1 b2 m ret sign mask;
+  local S0 S1 S2 S3 S4 T1 T2;
+
+  if(!isint(a) || !isint(b))
+      return newerror("toomcook3(a,b): a and/or b is not an integer");
+
+  alen = digits(a,2);
+  blen = digits(b,2);
+
+  sign = sgn(a) * sgn(b);
+  /* sgn(x) returns 0 if x = 0 */
+  if(sign == 0) return 0;
+
+  m = min(alen,blen)//3;
+  mask = ~-(1<<m);
+
+  /*
+    Cut-off at about 4,000 dec. digits
+    TODO: check
+  */
+  if(isdefined("test8900")){
+    if(m < 20) return a*b;
+  }
+  else{
+    if(m < 4096 ) return a*b;
+  }
+  a = abs(a);
+  b = abs(b);
+
+  a0 = a & mask;
+  a1 = (a>>m) & mask;
+  a2 = (a>>(2*m));
+
+  b0 =  b         & mask;
+  b1 = (b>>m)     & mask;
+  b2 = (b>>(2*m));
+
+  /*
+    Zimmermann
+  */
+
+  S0 = toomcook3(a0 , b0);
+  S1 = toomcook3((a2+a1+a0) , (b2+b1+b0));
+  S2 = toomcook3(((a2<<2)+(a1<<1)+a0) , ((b2<<2)+(b1<<1)+b0));
+  S3 = toomcook3((a2-a1+a0) , (b2-b1+b0));
+  S4 = toomcook3(a2,b2);
+  T1 = (S3<<1) + S2;
+  T1 /= 3;
+  T1 += S0;
+  T1 >>= 1;
+  T1 -= S4<<1;
+  T2 = (S1 + S3)>>1;
+  S1 -= T1;
+  S2 = T2 - S0 - S4;
+  S3 = T1 - T2;
+
+  ret =  (S4<<(4*m)) + (S3<<(3*m)) + (S2<<(2*m)) + (S1<<(1*m)) + S0;
+
+
+  ret = sign *ret;
+
+  return ret;
+}
+
+define toomcook3square(a){
+  local alen  a0 a1 a2  m tmp tmp2 ret sign S0 S1 S2 S3 S4 T1 mask;
+
+  if(!isint(a))return newerror("toomcook3square(a): a is not integer");
+
+  alen = digits(a,2);
+
+  sign = sgn(a) * sgn(a);
+  if(sign == 0) return 0;
+
+  m = alen//3;
+  mask = ~-(1<<m);
+  /*
+    Cut-off at about 5,000 dec. digits
+    TODO: check
+  */
+
+  if(isdefined("test8900")){
+    if(m < 20) return a^2;
+  }
+  else{
+    if(m < 5000 ) return a^2;
+  }
+
+  a = abs(a);
+
+  a0 = a & mask;
+  a1 = (a>>m) & mask;
+  a2 = (a>>(2*m));
+
+  /*
+   Bodrato/Zanoni
+   */
+  S0 = toomcook3square(a0);
+  S1 = toomcook3square(a2+a1+a0);
+  S2 = toomcook3square(a2-a1+a0);
+  S3 = toomcook3(a1<<1,a2);
+  S4 = toomcook3square(a2);
+
+  T1 = (S1 + S2)>>1;
+  S1 = S1 - T1 - S3;
+  S2 = T1 - S4 -S0;
+
+
+  S1 = S1<<(1*m);
+  S2 = S2<<(2*m);
+  S3 = S3<<(3*m);
+  S4 = S4<<(4*m);
+
+  ret = S0 + S1 + S2 + S3 + S4;
+  ret = sign *ret;
+
+  return ret;
+}
+
+define toomcook4(a,b)
+{
+
+  local  a0 a1 a2 a3 b0 b1 b2 b3 b4 ret tmp tmp2 tmp3 sign;
+  local  m alen blen mask;
+  local w1, w2, w3, w4, w5, w6, w7;
+
+  if(!isint(a) || !isint(b))
+    return newerror("toomcook4(a,b): a and/or b is not integer");
+
+  alen = digits(a,2);
+  blen = digits(b,2);
+
+  sign = sgn(a) * sgn(b);
+
+  if(sign == 0) return 0;
+
+  m = min(alen//4,blen//4);
+  mask = ~-(1<<m);
+
+  if(isdefined("test8900")){
+    if(m < 100) return toomcook3(a,b);
+  }
+  else{
+    if(m < 256*3072) return toomcook3(a,b);
+  }
+
+  a = abs(a);
+  b = abs(b);
+
+
+  a0 = a & mask;
+  a1 = (a>>m) & mask;
+  a2 = (a>>(2*m)) & mask;
+  a3 = (a>>(3*m));
+
+  b0 =  b         & mask;
+  b1 = (b>>m)     & mask;
+  b2 = (b>>(2*m)) & mask;
+  b3 = (b>>(3*m));
+
+  /*
+    Bodrato / Zanoni
+  */
+
+  w3 = a3 + (a1 + (a2 + a0));
+  w7 = b3 + (b1 + (b2 + b0));
+
+  w4 = -a3 + (-a1 + (a2 + a0));
+  w5 = -b3 + (-b1 + (b2 + b0));
+
+  w3 = toomcook4(w3, w7);
+  w4 = toomcook4(w4, w5);
+
+  w5 = a3 + ((a1<<2) + ((a2<<1) + (a0<<3)));
+  w2 = b3 + ((b1<<2) + ((b2<<1) + (b0<<3)));
+
+  w6 = -a3 + (-(a1<<2) + ((a2<<1) + (a0<<3)));
+  w7 = -b3 + (-(b1<<2) + ((b2<<1) + (b0<<3)));
+
+  w5 = toomcook4(w5, w2);
+  w6 = toomcook4(w6, w7);
+
+
+  w2 = (a3<<3) + ((a1<<1) + ((a2<<2) + a0));
+  w7 = (b3<<3) + ((b1<<1) + ((b2<<2) + b0));
+
+
+  w2 = toomcook4(w2, w7);
+
+  w1 = toomcook4(a3, b3);
+  w7 = toomcook4(a0, b0);
+
+  w2 = w2 + w5;
+  w6 = w5 - w6;
+  w4 = w3 - w4;
+  w5 = w5 - w1;
+  w5 -= w7 << 6;
+  w4 = w4>>1;
+  w3 = w3 - w4;
+  w5 = w5<<1;
+  w5 = w5 - w6;
+  w2 -= w3 * 65;
+  w3 = w3 - w7;
+  w3 = w3 - w1;
+  w2 += w3 * 45;
+  w5 -= w3<<3;
+  w5 = w5//24;
+  w6 =  w6 - w2;
+  w2 -= w4<<4;
+  w2 = w2//18;
+  w3 = w3 - w5;
+  w4 = w4 - w2;
+  w6 += w2 * 30;
+  w6 = w6//60;
+  w2 = w2 - w6;
+
+
+  ret = w7 + (w6<<m) + (w5<<(2*m)) + (w4<<(3*m))+ (w3<<(4*m))+
+        (w2<<(5*m))+ (w1<<(6*m));
+
+  ret = sign *ret;
+
+  return ret;
+}
+
+define toomcook4square(a){
+  local  a0 a1 a2 a3 ret S0 S1 S2 S3 S4 S5 S6 S7 tmp tmp2 tmp3;
+  local  sign m alen mask;
+  local T0 T1 T2 T3 T4 T5 T6 T7 T8;
+
+  if(!isint(a) )return newerror("toomcook3square(a): a is not integer");
+
+  alen = digits(a,2);
+
+  sign = sgn(a) * sgn(a);
+  /* sgn(x) returns 0 if x = 0 */
+  if(sign == 0) return 0;
+
+  m = (alen)//4;
+  mask = ~-( 1 << m );
+
+  /*
+    cut-off at about 2 mio. dec. digits
+    TODO: check!
+  */
+
+  if(isdefined("test8900")){
+    if(m < 100) return toomcook3square(a);
+  }
+  else{
+    if(m < 512*3072) return toomcook3square(a);
+  }
+
+  a = abs(a);
+
+  a0 = a          & mask;
+  a1 = (a>>m)     & mask;
+  a2 = (a>>(2*m)) & mask;
+  a3 = (a>>(3*m)) ;
+
+  /*
+    Bodrato / Zanoni
+  */
+
+  S1 = toomcook4square(a0);
+  S2 = toomcook4(a0<<1,a1);
+  S3 = toomcook4((a0 + a1 - a2 - a3 ) , (a0 - a1 - a2 + a3 ));
+  S4 = toomcook4square(a0 + a1 + a2 + a3 );
+  S5 = toomcook4( (a0 - a2 )<<1 , (a1 - a3 ));
+  S6 = toomcook4(a3<<1 , a2);
+  S7 = toomcook4square(a3);
+
+  T1 = S3 + S4;
+  T2 = (T1 + S5 )>>1;
+  T3 = S2 + S6;
+  T4 = T2 - T3;
+  T5 = T3 - S5;
+  T6 = T4 - S3;
+  T7 = T4 - S1;
+  T8 = T6 - S7;
+
+  ret =   (S7<<(6*m)) + (S6<<(5*m)) + (T7<<(4*m))
+        + (T5<<(3*m)) + (T8<<(2*m)) + (S2<<(1*m))  + S1;
+
+  ret = sign *ret;
+
+  return ret;
+}
+
+/*
+  TODO: Implement the asymmetric variations
+*/
+
+/*
+  produce_long_random_number(n) returns large pseudorandom numbers. Really large
+  numbers, e.g.:
+      produce_long_random_number(16)
+  is ca 4,128,561 bits (ca 1,242,821 dec. digits) large. Exact length is not
+  predeterminable because of the chaotic output of the function random().
+*/
+define __CZ__produce_long_random_number(n)
+{
+  local ret k;
+  ret = 1;
+  if(!isint(n) || n<1)
+    return newerror("__CZ__produce_long_random_number(n): "
+    		    "n is not an integer >=1");
+  for(k=0;k<n;k++){
+    ret += random();
+    ret  = toomcook4square(ret);
+  }
+  return ret;
+}
+
+
+/*
+ * restore internal function from resource debugging
+ * report important interface functions
+ */
+config("resource_debug", resource_debug_level),;
+if (config("resource_debug") & 3) {
+    print "toomcook3(a,b)";
+    print "toomcook3square(a)";
+    print "toomcook4(a,b)";
+    print "toomcook4square(a)";
+}