drygdryg/calc
1
0
Fork 0
mirror of https://github.com/lcn2/calc.git synced 2025-08-16 01:03:29 +03:00
Code Issues Packages Projects Releases Wiki Activity

Release calc version 2.12.4.12

Browse Source
This commit is contained in:
Landon Curt Noll
2013-09-01 20:07:26 -07:00
parent 85bfa30897
commit 57a22a6f39
346 changed files with 2444 additions and 502 deletions
Download Patch File Download Diff File Expand all files Collapse all files

58
cal/Makefile
View File

@@ -18,8 +18,8 @@
# 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.7 $
# @(#) $Id: Makefile,v 30.7 2013/08/11 09:07:26 chongo Exp $
# @(#) $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
@@ -179,19 +179,32 @@ MV= mv
CO= co
TRUE= true
TOUCH= touch
SED= sed
SORT= sort
FMT= fmt
# The calc files to install
#
CALC_FILES= 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 factorial2.cal factorial.cal \
gvec.cal hello.cal hms.cal intfile.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 README \
regress.cal repeat.cal screen.cal seedrandom.cal set8700.cal \
set8700.line solve.cal specialfunctions.cal statistics.cal sumsq.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 \
@@ -244,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

28
cal/README
View File

@@ -1382,25 +1382,17 @@ statistics.cal
strings.cal
toupper(s)
tolower(s)
strcasecmp(s1,s2)
strncasecmp(s1,s2,length)
isascii(c)
isalnum(c)
isalpha(c)
iscntrl(c)
isdigit(c)
isgraph(c)
islower(c)
isprint(c)
ispunct(c)
isspace(c)
isupper(c)
isblank(c)
isxdigit(c)
Implements most of the functions of libc's ctype.h and strings.h.
Implements some of the functions of libc's ctype.h and strings.h.
NOTE: A number of the ctype.h and strings.h functions are now builtin
functions in calc.
WARNING: If the remaining functions in this calc resource file become
calc builtin functions, then strings.cal may be removed in
a future release.
sumsq.cal
@@ -1792,8 +1784,8 @@ zeta2.cal
## 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.6 $
## @(#) $Id: README,v 30.6 2013/08/18 20:01:53 chongo Exp $
## @(#) $Revision: 30.8 $
## @(#) $Id: README,v 30.8 2013/09/02 01:46:05 chongo Exp $
## @(#) $Source: /usr/local/src/bin/calc/cal/RCS/README,v $
##
## Under source code control: 1990/02/15 01:50:32

88
cal/infinities.cal Normal file
View File

@@ -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()";
}

728
cal/intnum.cal Normal file
View File

@@ -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,[...])";
}

214
cal/regress.cal
View File

@@ -17,8 +17,8 @@
* 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.10 $
* @(#) $Id: regress.cal,v 30.10 2013/09/01 22:08:44 chongo Exp $
* @(#) $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
@@ -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()
{
@@ -1445,91 +1447,10 @@ define test_functions()
vrfy(jacobi(-1,-1) == 0, '1236: jacobi(-1,-1) == 0');
vrfy(jacobi(0,-1) == 0, '1237: jacobi(0,-1) == 0');
/* ctype function tests */
vrfy(isalnum("A") == 1, '1238: isalnum("A") == 1');
vrfy(isalnum("a") == 1, '1239: isalnum("a") == 1');
vrfy(isalnum("2") == 1, '1239: isalnum("2") == 1');
vrfy(isalnum("\t") == 0, '1240: isalnum("\t") == 0');
vrfy(isalpha("A") == 1, '1241: isalpha("A") == 1');
vrfy(isalpha("a") == 1, '1242: isalpha("a") == 1');
vrfy(isalpha("2") == 0, '1243: isalpha("2") == 0');
vrfy(isalpha("\t") == 0, '1244: isalpha("\t") == 0');
vrfy(iscntrl("A") == 0, '1245: iscntrl("A") == 0');
vrfy(iscntrl("a") == 0, '1246: iscntrl("a") == 0');
vrfy(iscntrl("2") == 0, '1247: iscntrl("2") == 0');
vrfy(iscntrl("\t") == 1, '1248: iscntrl("\t") == 1');
vrfy(isdigit("A") == 0, '1249: isdigit("A") == 0');
vrfy(isdigit("a") == 0, '1250: isdigit("a") == 0');
vrfy(isdigit("2") == 1, '1251: isdigit("2") == 1');
vrfy(isdigit("\t") == 0, '1252: isdigit("\t") == 0');
vrfy(isgraph("A") == 1, '1253: isgraph("A") == 1');
vrfy(isgraph("a") == 1, '1254: isgraph("a") == 1');
vrfy(isgraph("2") == 1, '1255: isgraph("2") == 1');
vrfy(isgraph("\t") == 0, '1255: isgraph("\t") == 0');
vrfy(islower("A") == 0, '1256: islower("A") == 0');
vrfy(islower("a") == 1, '1257: islower("a") == 1');
vrfy(islower("1") == 0, '1258: islower("1") == 0');
vrfy(isprint("A") == 1, '1259: isprint("A") == 1');
vrfy(isprint("a") == 1, '1260: isprint("a") == 1');
vrfy(isprint(" ") == 1, '1261: isprint(" ") == 1');
vrfy(isprint("\t") == 0, '1262: isprint("\t") == 0');
vrfy(ispunct("A") == 0, '1263: ispunct("A") == 0');
vrfy(ispunct("a") == 0, '1264: ispunct("a") == 0');
vrfy(ispunct(" ") == 0, '1265: ispunct(" ") == 0');
vrfy(ispunct("?") == 1, '1266: ispunct("?") == 1');
vrfy(isspace("A") == 0, '1267: isspace("A") == 0');
vrfy(isspace("Krik") == 0, '1268: isspace("Krik") == 0');
vrfy(isspace(" ") == 1, '1269: isspace(" ") == 1');
vrfy(isspace("?") == 0, '1270: isspace("?") == 0');
vrfy(isupper("A") == 1, '1271: isupper("A") == 1');
vrfy(isupper("a") == 0, '1272: isupper("a") == 0');
vrfy(isupper("1") == 0, '1273: isupper("1") == 0');
vrfy(isxdigit("A") == 1, '1274: isxdigit("A") == 1');
vrfy(isxdigit("f") == 1, '1275: isxdigit("f") == 1');
vrfy(isxdigit("2") == 1, '1276: isxdigit("2") == 1');
vrfy(isxdigit("x") == 0, '1277: isxdigit("x") == 0');
vrfy(strcasecmp("ab", "aBc") == -1,
'1278: strcasecmp("ab", "aBc") == -1');
vrfy(strcasecmp("abc", "aBb") == 1,
'1279: strcasecmp("abc", "aBb") == 1');
vrfy(strcasecmp("abc", "abc") == 0,
'1280: strcasecmp("abc", "abc") == 0');
vrfy(strcasecmp("abc", "aBc") == 0,
'1281: strcasecmp("abc", "aBc") == 0');
vrfy(strcasecmp("abc", "aBd") == -1,
'1282: strcasecmp("abc", "aBd") == -1');
vrfy(strcasecmp("abc\0", "aBc") == 1,
'1283: strcasecmp("abc\0", "aBc") == 1');
vrfy(strcasecmp("a\0b", "A\0c") == -1,
'1284: strcasecmp("a\0b", "A\0c") == -1');
vrfy(strncasecmp("abc", "xyz", 0) == 0,
'1285: strncasecmp("abc", "xyz", 0) == 0');
vrfy(strncasecmp("abc", "xyz", 1) == -1,
'1286: strncasecmp("abc", "xyz", 1) == -1');
vrfy(strncasecmp("abc", "", 1) == 1,
'1287: strncasecmp("abc", "", 1) == 1');
vrfy(strncasecmp("a", "b", 2) == -1,
'1288: strncasecmp("a", "b", 2) == -1');
vrfy(strncasecmp("ab", "Ac", 2) == -1,
'1289: strncasecmp("ab", "Ac", 2) == -1');
vrfy(strncasecmp("\0ac", "\0b", 2) == -1,
'1290: strncasecmp("\0ac", "\0b", 2) == -1');
vrfy(strncasecmp("ab", "aBc", 2) == 0,
'1291: strncasecmp("ab", "aBc", 2) == 0');
vrfy(strncasecmp("abc", "abd", 2) == 0,
'1292: strncasecmp("abc", "abd", 2) == 0');
/*
* NOTE: Function tests are continued in test_functionss()
* starting at test 9000.
*/
print '1293: Ending test_functions';
}
@@ -8015,9 +7936,128 @@ 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
*

71
cal/smallfactors.cal Normal file
View File

@@ -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)";
}

41
cal/strings.cal Normal file
View File

@@ -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)";
}