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

Release calc version 2.12.5.5

Browse Source
This commit is contained in:
Landon Curt Noll
2017-05-20 19:24:06 -07:00
parent 2726ae9d23
commit 8dd7a3cd2a
23 changed files with 1126 additions and 422 deletions
Download Patch File Download Diff File Expand all files Collapse all files

57
cal/Makefile
View File

@@ -2,7 +2,7 @@
#
# cal - makefile for calc standard resource files
#
# Copyright (C) 1999-2006 Landon Curt Noll
# Copyright (C) 1999-2006,2017 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
@@ -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.11 $
# @(#) $Id: Makefile,v 30.11 2014/09/02 07:14:49 chongo Exp $
# @(#) $Source: /usr/local/src/bin/calc/cal/RCS/Makefile,v $
# @(#) $Revision: 30.12 $
# @(#) $Id: Makefile,v 30.12 2017/05/19 16:09:14 chongo Exp $
# @(#) $Source: /usr/local/src/bin/calc-RHEL7/cal/RCS/Makefile,v $
#
# Under source code control: 1991/07/21 05:00:54
# File existed as early as: 1991
@@ -199,7 +199,7 @@ CALC_FILES= README alg_config.cal beer.cal bernoulli.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 \
lucas.cal lucas_chk.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 \
@@ -211,6 +211,10 @@ CALC_FILES= README alg_config.cal beer.cal bernoulli.cal \
test8400.cal test8500.cal test8600.cal test8900.cal toomcook.cal \
unitfrac.cal varargs.cal xx_print.cal zeta2.cal
# These calc files are now obsolete and are removed by the install rule.
#
DEAD_CALC_FILES= lucas_tbl.cal
# These files are found (but not built) in the distribution
#
DISTLIST= ${CALC_FILES} ${MAKE_FILE}
@@ -301,6 +305,23 @@ clean:
clobber: clean
${RM} -f .all
-${Q} if [ -e .DS_Store ]; then \
echo ${RM} -rf .DS_Store; \
${RM} -rf .DS_Store; \
fi
-${Q} for i in ${DEAD_CALC_FILES} /dev/null; do \
if [ "$$i" = "/dev/null" ]; then \
continue; \
fi; \
if [ -e "${T}${CALC_SHAREDIR}/$$i" ]; then \
echo "${RM} -f ${T}${CALC_SHAREDIR}/$$i"; \
${RM} -f ${T}${CALC_SHAREDIR}/$$i; \
fi; \
if [ -e "./$$i" ]; then \
echo "${RM} -f ./$$i"; \
${RM} -f ./$$i; \
fi; \
done
# install everything
#
@@ -333,12 +354,38 @@ install: all
echo "installed ${T}${CALC_SHAREDIR}/$$i"; \
fi; \
done
${Q} for i in ${DEAD_CALC_FILES} /dev/null; do \
if [ "$$i" = "/dev/null" ]; then \
continue; \
fi; \
if [ -e "${T}${CALC_SHAREDIR}/$$i" ]; then \
echo "${RM} -f ${T}${CALC_SHAREDIR}/$$i"; \
${RM} -f ${T}${CALC_SHAREDIR}/$$i; \
fi; \
if [ -e "./$$i" ]; then \
echo "${RM} -f ./$$i"; \
${RM} -f ./$$i; \
fi; \
done
# Try to remove everything that was installed
#
# NOTE: Keep the uninstall rule in reverse order to the install rule
#
uninstall:
- ${Q} for i in ${DEAD_CALC_FILES} /dev/null; do \
if [ "$$i" = "/dev/null" ]; then \
continue; \
fi; \
if [ -e "${T}${CALC_SHAREDIR}/$$i" ]; then \
echo "${RM} -f ${T}${CALC_SHAREDIR}/$$i"; \
${RM} -f ${T}${CALC_SHAREDIR}/$$i; \
fi; \
if [ -e "./$$i" ]; then \
echo "${RM} -f ./$$i"; \
${RM} -f ./$$i; \
fi; \
done
-${Q} for i in ${CALC_FILES} /dev/null; do \
if [ "$$i" = "/dev/null" ]; then \
continue; \

38
cal/README
View File

@@ -14,9 +14,12 @@ For example:
; read lucas
lucas(h,n) defined
gen_u2(h,n,v1) defined
gen_u0(h,n,v1) defined
rodseth_xhn(x,h,n) defined
gen_v1(h,n) defined
ldebug(funct,str) defined
legacy_gen_v1(h,n) defined
will cause calc to load and execute the 'lucas.cal' resource file.
Executing the resource file will cause several functions to be defined.
@@ -792,7 +795,27 @@ lucas.cal
lucas(h, n)
Perform a primality test of h*2^n-1, with 1<=h<2*n.
Perform a primality test of h*2^n-1.
gen_u2(h, n, v1)
Generate u(2) for h*2^n-1. This function is used by lucas(h, n),
as the first term in the lucas sequence that is needed to
prove that h*2^n-1 is prime or not prime.
NOTE: Some call this term u(0). The function gen_u0(h, n, v1)
simply calls gen_u2(h, n, v1) for such people. :-)
gen_v1(h, v)
Generate v(1) for h*2^n-1. This function is used by lucas(h, n),
via the gen_u2(h, n, v1), to supply the 3rd argument to gen_u2.
legacy_gen_v1(h, n)
Generate v(1) for h*2^n-1 using the legacy Amdahl 6 method.
This function sometimes returns -1 for a few cases when
h is a multiple of 3. This function is NOT used by lucas(h, n).
lucas_chk.cal
@@ -805,11 +828,6 @@ lucas_chk.cal
Used by regress.cal during the 2100 test set.
lucas_tbl.cal
Lucasian criteria for primality tables.
mersenne.cal
mersenne(p)
@@ -1798,7 +1816,7 @@ zeta2.cal
for information on this special zeta function.
## Copyright (C) 2000,2014 David I. Bell and Landon Curt Noll
## Copyright (C) 2000,2014,2017 David I. Bell and Landon Curt Noll
##
## Primary author: Landon Curt Noll
##
@@ -1816,9 +1834,9 @@ 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.9 $
## @(#) $Id: README,v 30.9 2014/10/06 08:44:18 chongo Exp $
## @(#) $Source: /usr/local/src/bin/calc/cal/RCS/README,v $
## @(#) $Revision: 30.10 $
## @(#) $Id: README,v 30.10 2017/05/19 16:09:14 chongo Exp $
## @(#) $Source: /usr/local/src/bin/calc-RHEL7/cal/RCS/README,v $
##
## Under source code control: 1990/02/15 01:50:32
## File existed as early as: before 1990

666
cal/lucas.cal
View File

@@ -1,7 +1,7 @@
/*
* lucas - perform a Lucas primality test on h*2^n-1
*
* Copyright (C) 1999 Landon Curt Noll
* Copyright (C) 1999,2017 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
@@ -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.2 $
* @(#) $Id: lucas.cal,v 30.2 2013/09/27 08:58:46 chongo Exp $
* @(#) $Source: /usr/local/src/bin/calc/cal/RCS/lucas.cal,v $
* @(#) $Revision: 30.4 $
* @(#) $Id: lucas.cal,v 30.4 2017/05/20 21:54:16 chongo Exp $
* @(#) $Source: /usr/local/src/bin/calc-RHEL7/cal/RCS/lucas.cal,v $
*
* Under source code control: 1990/05/03 16:49:51
* File existed as early as: 1990
@@ -28,6 +28,12 @@
* Share and enjoy! :-) http://www.isthe.com/chongo/tech/comp/calc/
*/
/*
* For a general tutorial on how to find a new largest known prime, see:
*
* http://www.isthe.com/chongo/tech/math/prime/prime-tutorial.pdf
*/
/*
* NOTE: This is a standard calc resource file. For information on calc see:
*
@@ -71,10 +77,15 @@
* NOTE: Both largest known and largest known twin prime records have been
* broken. Rather than update this file each time, I'll just
* congratulate the finders and encourage others to try for
* larger finds. Records were made to be broken afterall!
* larger finds. Records were made to be broken after all!
*/
/* ON GAINING A WORLD RECORD:
/*
* ON GAINING A WORLD RECORD:
*
* For a general tutorial on how to find a new largest known prime, see:
*
* http://www.isthe.com/chongo/tech/math/prime/prime-tutorial.pdf
*
* The routines in calc were designed to be portable, and to work on
* numbers of 'sane' size. The Amdahl 6 team used a 'ultra-high speed
@@ -83,6 +94,13 @@
* The heart of the package was a multiplication and square routine that
* was based on the PFA Fast Fourier Transform and on Winograd's radix FFTs.
*
* NOTE: While the PFA Fast Fourier Transform and Winograd's radix FFTs
* might have been optimal for the Amdahl 6 team at the time,
* they might not be optimal for your CPU architecture. See
* the above mentioned tutorial for information on better
* methods of performing multiplications and squares of very
* large numbers.
*
* Having a fast computer, and a good multi-precision package are
* critical, but one also needs to know where to look in order to have
* a good chance at a record. Knowing what to test is beyond the scope
@@ -139,12 +157,10 @@
* be the factors of another candidate.
*
* Finally, one should eliminate all values of 'h*2^n-1' where
* 'h*2^n+1' is divisible by a small primes. The ideas behind this
* point is beyond the scope of this program.
* 'h*2^n+1' is divisible by a small primes.
*/
global pprod256; /* product of "primes up to 256" / "primes up to 46" */
pprod256 = 0; /* product of "primes up to 256" / "primes up to 46" */
/*
* lucas - lucas primality test on h*2^n-1
@@ -171,15 +187,33 @@ global pprod256; /* product of "primes up to 256" / "primes up to 46" */
* "Introduction to Analytic Number Theory", by Tom A. Apostol,
* Springer-Verlag, 1984, p 188.
*
* An excellent 5-page paper by Oystein J. Rodseth (we apologize that the
* ASCII character set does not allow us to spell his name with the
* umlaut marks on the O's):
*
* NOTE: The original Amdahl 6 method predates the publication of Ref4.
* The gen_v1() function used by lucas() uses the Ref4 method.
* See the 'Amdahl 6 legacy code' section below for the original
* method of generating v(1).
*
* Ref4:
*
* "A note on primality tests for N = h*2^n-1", by Oystein J. Rodseth,
* Department of Mathematics, University of Bergen, BIT Numerical
* Mathematics. 34 (3): pp 451-454.
*
* http://folk.uib.no/nmaoy/papers/luc.pdf
*
* This test is performed as follows: (see Ref1, Theorem 5)
*
* a) generate u(0) (see the function gen_u0() below)
* a) generate u(2) (see the function gen_u2() below)
* (NOTE: some call this u(0))
*
* b) generate u(n-2) according to the rule:
* b) generate u(n) according to the rule:
*
* u(i+1) = u(i)^2-2 mod h*2^n-1
*
* c) h*2^n-1 is prime if and only if u(n-2) == 0 Q.E.D. :-)
* c) h*2^n-1 is prime if and only if u(n) == 0 Q.E.D. :-)
*
* Now the following conditions must be true for the test to work:
*
@@ -188,7 +222,7 @@ global pprod256; /* product of "primes up to 256" / "primes up to 46" */
* h < 2^n
* h mod 2 == 1
*
* A few misc notes:
* A few miscellaneous notes:
*
* In order to reduce the number of tests, as attempt to eliminate
* any number that is divisible by a prime less than 257. Valid prime
@@ -222,7 +256,7 @@ lucas(h, n)
local testval; /* h*2^n-1 */
local shiftdown; /* the power of 2 that divides h */
local u; /* the u(i) sequence value */
local v1; /* the v(1) generator of u(0) */
local v1; /* the v(1) generator of u(2) */
local i; /* u sequence cycle number */
local oldh; /* pre-reduced h */
local oldn; /* pre-reduced n */
@@ -364,18 +398,17 @@ lucas(h, n)
}
/*
* try to compute u(0)
* try to compute u(2) (NOTE: some call this u(0))
*
* We will use gen_v1() to give us a v(1) using the values
* of 'h' and 'n'. We will then use gen_u0() to convert
* the v(1) into u(0).
* of 'h' and 'n'. We will then use gen_u2() to convert
* the v(1) into u(2).
*
* If gen_v1() returns a negative value, then we failed to
* generate a test for h*2^n-1. This is because h mod 3 == 0
* is hard to do, and in rare cases, exceed the tables found
* in this program. We will generate an message and assume
* the number is not prime, even though if we had a larger
* table, we might have been able to show that it is prime.
* generate a test for h*2^n-1. The legacy function,
* legacy_gen_v1() used by the Amdahl 6 could have returned
* -1. The new gen_v1() based on the method outlined in Ref4
* will never return -1.
*/
v1 = gen_v1(h, n);
if (v1 < 0) {
@@ -384,10 +417,10 @@ lucas(h, n)
ldebug("lucas", "unknown: no v(1)");
return -1;
}
u = gen_u0(h, n, v1);
u = gen_u2(h, n, v1);
/*
* compute u(n-2)
* compute u(n) (NOTE: some call this u(n-2))
*/
for (i=3; i <= n; ++i) {
/* u = (u^2 - 2) % testval; */
@@ -407,11 +440,19 @@ lucas(h, n)
}
/*
* gen_u0 - determine the initial Lucas sequence for h*2^n-1
* gen_u2 - determine the initial Lucas sequence for h*2^n-1
*
* Historically many start the Lucas sequence with u(0).
* Some, like the author of this code, prefer to start
* with U(2). This is so one may say:
*
* 2^p-1 is prime if u(p) = 0 mod 2^p-1
* or:
* h*2^p-1 is prime if u(p) = 0 mod h*2^p-1
*
* According to Ref1, Theorem 5:
*
* u(0) = alpha^h + alpha^(-h)
* u(2) = alpha^h + alpha^(-h) (NOTE: Ref1 calls it u(0))
*
* Now:
*
@@ -419,7 +460,7 @@ lucas(h, n)
*
* Therefore:
*
* u(0) = v(h)
* u(2) = v(h) (NOTE: Ref1 calls it u(0))
*
* We calculate v(h) as follows: (Ref1, top of page 873)
*
@@ -447,11 +488,11 @@ lucas(h, n)
* v1 - gen_v1(h,n) (see function below)
*
* returns:
* u(0) - initial value for Lucas test on h*2^n-1
* -1 - failed to generate u(0)
* u(2) - initial value for Lucas test on h*2^n-1
* -1 - failed to generate u(2)
*/
define
gen_u0(h, n, v1)
gen_u2(h, n, v1)
{
local shiftdown; /* the power of 2 that divides h */
local r; /* low value: v(n) */
@@ -500,7 +541,7 @@ gen_u0(h, n, v1)
* at least 2 bits long for the loop below to work.
*/
if (h == 1) {
ldebug("gen_u0", "quick h == 1 case");
ldebug("gen_u2", "quick h == 1 case");
/* return r%(h*2^n-1); */
return hnrmod(r, h, n, -1);
}
@@ -540,21 +581,502 @@ gen_u0(h, n, v1)
return r;
}
/*
* gen_u0 - determine the initial Lucas sequence for h*2^n-1
*
* Historically many start the Lucas sequence with u(0).
* Some, like the author of this code, prefer to start
* with u(2). This is so one may say:
*
* 2^p-1 is prime if u(p) = 0 mod 2^p-1
* or:
* h*2^n-1 is prime if U(n) = 0 mod h*2^n-1
*
* For those using the old code with gen_u0(), we
* simply call gen_u2() instead.
*
* See the function gen_u2() for details.
*
* input:
* h - h as in h*2^n-1
* n - n as in h*2^n-1
* v1 - gen_v1(h,n) (see function below)
*
* returns:
* u(2) - initial value for Lucas test on h*2^n-1
* -1 - failed to generate u(2)
*/
define
gen_u0(h, n, v1)
{
return gen_u2(h, n, v1);
}
/*
* rodseth_xhn - determine if v(1) == x for h*2^n-1
*
* For a given h*2^n-1, v(1) == x if:
*
* jacobi(x-2, h*2^n-1) == 1 (Ref4, condition 1) part 1
* jacobi(x+2, h*2^n-1) == -1 (Ref4, condition 1) part 2
*
* Now when x-2 <= 0:
*
* jacobi(x-2, h*2^n-1) == 0
*
* because:
*
* jacobi(x,y) == 0 if x <= 0
*
* So for (Ref4, condition 1) part 1 to be true:
*
* x-2 > 0
*
* And therefore:
*
* x > 2
*
* input:
* x - potential v(1) value
* h - h as in h*2^n-1
* n - n as in h*2^n-1
*
* returns:
* 1 if v(1) == x for h*2^n-1
* 0 otherwise
*/
define
rodseth_xhn(x, h, n)
{
local testval; /* h*2^n-1 */
/*
* check arg types
*/
if (!isint(h)) {
quit "bad args: h must be an integer";
}
if (!isint(n)) {
quit "bad args: n must be an integer";
}
if (!isint(x)) {
quit "bad args: x must be an integer";
}
/*
* firewall
*/
if (x <= 2) {
return 0;
}
/*
* Check for jacobi(x-2, h*2^n-1) == 1 (Ref4, condition 1) part 1
*/
testval = h*2^n-1;
if (jacobi(x-2, testval) != 1) {
return 0;
}
/*
* Check for jacobi(x+2, h*2^n-1) == -1 (Ref4, condition 1) part 2
*/
if (jacobi(x+2, testval) != -1) {
return 0;
}
/*
* v(1) == x for this h*2^n-1
*/
return 1;
}
/*
* Trial tables used by gen_v1()
*
* When h mod 3 == 0, one needs particular values of D, a and b (see gen_v1
* documentation) in order to find a value of v(1).
* When h mod 3 == 0, according to Ref4 we need to find the first value X where:
*
* This table defines 'quickmax' possible tests to be taken in ascending
* order. The v1_qval[x] refers to a v(1) value from Ref1, Table 1. A
* related D value is found in d_qval[x]. All D values expect d_qval[1]
* are also taken from Ref1, Table 1. The case of D == 21 as listed in
* Ref1, Table 1 can be changed to D == 7 for the sake of the test because
* of {note 6}.
* jacobi(X-2, h*2^n-1) == 1 (Ref4, condition 1) part 1
* jacobi(X+2, h*2^n-1) == -1 (Ref4, condition 1) part 2
*
* We can show that X > 2. See the comments in the rodseth_xhn(x,h,n) above.
*
* Some values of X satisfy more often than others. For example a large sample
* of odd h, h multiple of 3 and large n (some around 1e4, some near 1e6, others
* near 3e7) where the sample size was 66 973 365, here is the count of the
* smallest value of X that satisfies conditions in Ref4, condition 1:
*
* count X
* ----------
* 26791345 3
* 17223016 5
* 7829600 9
* 6988774 11
* 3301093 15
* 1517149 17
* 910346 21
* 711791 29
* 573403 20
* 390395 27
* 288637 35
* 149751 36
* 107733 39
* 58743 41
* 35619 45
* 25052 32
* 17775 51
* 13031 44
* 7563 56
* 7540 49
* 7060 59
* 4407 57
* 2948 65
* 2502 55
* 2388 69
* 2094 71
* 689 77
* 626 81
* 491 66
* 426 95
* 219 80
* 203 67
* 185 84
* 152 99
* 127 72
* 102 74
* 98 87
* 67 90
* 55 104
* 48 101
* 32 105
* 17 109
* 16 116
* 15 111
* 13 92
* 12 125
* 7 129
* 3 146
* 2 140
* 2 120
* 1 165
* 1 161
* 1 155
*
* The above distribution was found to hold fairly well over many values of
* odd h that are a multiple of 3 and for many values of n where h < 2^n.
*
* Given this information, when odd h is a multiple of 3 we try, in order,
* these values of X:
*
* 3, 5, 9, 11, 15, 17, 21, 29, 20, 27, 35, 36, 39, 41, 45, 32, 51, 44,
* 56, 49, 59, 57, 65, 55, 69, 71, 77, 81, 66, 95, 80, 67, 84, 99, 72,
* 74, 87, 90, 104, 101, 105, 109, 116, 111, 92
*
* And stop on the first value of X where:
*
* jacobi(X-2, h*2^n-1) == 1
* jacobi(X+2, h*2^n-1) == -1
*
* If no value in that list works, we start simple search starting with X = 120
* and incrementing by 1 until a value of X is found.
*
* The x_tbl[] matrix contains those common values of X to try in order.
* If all x_tbl_len fail to satisfy Ref4 condition 1, then we begin a
* linear search at next_x until we find a proper X value.
*
* IMPORTANT NOTE: Using this table will not find the smallest possible v(1)
* for a given h and n. This is not a problem because using
* a larger value of v(1) does not impact the primality test.
* Furthermore after lucas(h, n) generates a few u(n) terms,
* the values will wrap (due to computing mod h*2^n-1).
* Finally on average, about 1/4 of the values of X work as
* v(1) for a given n when h is a multiple of 3. Skipping
* rarely used v(1) will not doom gen_v1() to a long search.
*/
x_tbl_len = 45;
mat x_tbl[x_tbl_len];
x_tbl = {
3, 5, 9, 11, 15, 17, 21, 29, 20, 27, 35, 36, 39, 41, 45, 32, 51, 44,
56, 49, 59, 57, 65, 55, 69, 71, 77, 81, 66, 95, 80, 67, 84, 99, 72,
74, 87, 90, 104, 101, 105, 109, 116, 111, 92
};
next_x = 120;
/*
* gen_v1 - compute the v(1) for a given h*2^n-1 if we can
*
* This function assumes:
*
* n > 2 (n==2 has already been eliminated)
* h mod 2 == 1
* h < 2^n
* h*2^n-1 mod 3 != 0 (h*2^n-1 has no small factors, such as 3)
*
* The generation of v(1) depends on the value of h. There are two cases
* to consider, h mod 3 != 0, and h mod 3 == 0.
*
***
*
* Case 1: (h mod 3 != 0)
*
* This case is easy.
*
* In Ref1, page 869, one finds that if: (or see Ref2, page 131-132)
*
* h mod 6 == +/-1
* h*2^n-1 mod 3 != 0
*
* which translates, gives the functions assumptions, into the condition:
*
* h mod 3 != 0
*
* If this case condition is true, then:
*
* u(2) = (2+sqrt(3))^h + (2-sqrt(3))^h (see Ref1, page 869)
* = (2+sqrt(3))^h + (2+sqrt(3))^(-h) (NOTE: some call this u(2))
*
* and since Ref1, Theorem 5 states:
*
* u(2) = alpha^h + alpha^(-h) (NOTE: some call this u(2))
* r = abs(2^2 - 1^2*3) = 1
*
* and the bottom of Ref1, page 872 states:
*
* v(x) = alpha^x + alpha^(-x)
*
* If we let:
*
* alpha = (2+sqrt(3))
*
* then
*
* u(2) = v(h) (NOTE: some call this u(2))
*
* so we simply return
*
* v(1) = alpha^1 + alpha^(-1)
* = (2+sqrt(3)) + (2-sqrt(3))
* = 4
*
***
*
* Case 2: (h mod 3 == 0)
*
* For the case where h is a multiple of 3, we turn to Ref4.
*
* The central theorem on page 3 of that paper states that
* we may set v(1) to the first value X that satisfies:
*
* jacobi(X-2, h*2^n-1) == 1 (Ref4, condition 1)
* jacobi(X+2, h*2^n-1) == -1 (Ref4, condition 1)
*
* NOTE: Ref4 uses P, which we shall refer to as X.
* Ref4 uses N, which we shall refer to as h*2^n-1.
*
* NOTE: Ref4 uses the term Legendre-Jacobi symbol, which
* we shall refer to as the Jacobi symbol.
*
* Before we address the two conditions, we need some background information
* on two symbols, Legendre and Jacobi. In Ref 2, pp 278, 284-285, we find
* the following definitions of jacobi(a,b) and L(a,p):
*
* The Legendre symbol L(a,p) takes the value:
*
* L(a,p) == 1 => a is a quadratic residue of p
* L(a,p) == -1 => a is NOT a quadratic residue of p
*
* when:
*
* p is prime
* p mod 2 == 1
* gcd(a,p) == 1
*
* The value a is a quadratic residue of b if there exists some integer z
* such that:
*
* z^2 mod b == a
*
* The Jacobi symbol jacobi(a,b) takes the value:
*
* jacobi(a,b) == 1 => b is not prime,
* or a is a quadratic residue of b
* jacobi(a,b) == -1 => a is NOT a quadratic residue of b
*
* when
*
* b mod 2 == 1
* gcd(a,b) == 1
*
* It is worth noting for the Legendre symbol, in order for L(X+/-2,
* h*2^n-1) to be defined, we must ensure that neither X-2 nor X+2 are
* factors of h*2^n-1. This is done by pre-screening h*2^n-1 to not
* have small factors and keeping X+2 less than that small factor
* limit. It is worth noting that in lucas(h, n), we first verify
* that h*2^n-1 does not have a factor < 257 before performing the
* primality test. So while X+/-2 < 257, we know that
* gcd(X+/-2, h*2^n-1) == 1.
*
* Returning to the testing of conditions in Ref4, condition 1:
*
* jacobi(X-2, h*2^n-1) == 1
* jacobi(X+2, h*2^n-1) == -1
*
* When such an X is found, we set:
*
* v(1) = X
*
***
*
* In conclusion, we can compute v,(1) by attempting to do the following:
*
* h mod 3 != 0
*
* we return:
*
* v(1) == 4
*
* h mod 3 == 0
*
* we return:
*
* v(1) = X
*
* where X > 2 in a integer such that:
*
* jacobi(X-2, h*2^n-1) == 1
* jacobi(X+2, h*2^n-1) == -1
*
***
*
* input:
* h h as in h*2^n-1
* n n as in h*2^n-1
*
* output:
* returns v(1), or -1 is there is no quick way
*/
define
gen_v1(h, n)
{
local x; /* potential v(1) to test */
local i; /* x_tbl index */
/*
* check arg types
*/
if (!isint(h)) {
quit "bad args: h must be an integer";
}
if (!isint(n)) {
quit "bad args: n must be an integer";
}
/*
* check for Case 1: (h mod 3 != 0)
*/
if (h % 3 != 0) {
/* v(1) is easy to compute */
return 4;
}
/*
* What follow is Case 2: (h mod 3 == 0)
*/
/*
* We will look for x that satisfies conditions in Ref4, condition 1:
*
* jacobi(X-2, h*2^n-1) == 1 part 1
* jacobi(X+2, h*2^n-1) == -1 part 2
*/
for (i=0; i < x_tbl_len; ++i) {
/*
* test Ref4 condition 1:
*/
x = x_tbl[i];
if (rodseth_xhn(x, h, n) == 1) {
/*
* found a x that satisfies Ref4 condition 1
*/
ldebug("gen_v1", "h= " + str(h) + " n= " + str(n) +
" v1= " + str(x) + " using tbl[ " +
str(i) + " ]");
return x;
}
}
/*
* We are in that rare case (about 1 in 2 300 000) where none of the
* common X values satisfy Ref4 condition 1. We start a linear search
* at next_x from here on.
*
* However, we also need to keep in mind that when x+2 >= 257, we
* need to verify that gcd(x-2, h*2^n-1) == 1 and
* and to verify that gcd(x+2, h*2^n-1) == 1.
*/
x = next_x;
while (rodseth_xhn(x, h, n) != 1) {
++x;
}
/* finally found a v(1) value */
ldebug("gen_v1", "h= " + str(h) + " n= " + str(n) +
" v1= " + str(x) + " beyond tbl");
return x;
}
/*
* ldebug - print a debug statement
*
* input:
* funct name of calling function
* str string to print
*/
define
ldebug(funct, str)
{
if (config("resource_debug") & 8) {
print "DEBUG:", funct:":", str;
}
return;
}
/*
************************
* Amdahl 6 legacy code *
************************
*
* NOTE: What follows is legacy code based on the method used by the
* Amdahl 6 group:
*
* John Brown, Landon Curt Noll, Bodo Parady, Gene Smith,
* Joel Smith and Sergio Zarantonello
*
* This method generated v(1) for nearly all values, except for a
* few rare cases when h mod 3 == 0. The code is NOT used by lucas.cal
* above. The gen_v1() function above is based on an improved method
* outlined in Ref4. That method generated v(1) for all h.
*
* The code below is kept for historical purposes only. The functions
* and global variables of the Amdahl 6 legacy code all begin with legacy_.
*/
/*
* Trial tables used by legacy_gen_v1()
*
* When h mod 3 == 0, one needs particular values of D, a and b (see
* legacy_gen_v1 documentation) in order to find a value of v(1).
*
* This table defines 'legacy_quickmax' possible tests to be taken in ascending
* order. The legacy_v1_qval[x] refers to a v(1) value from Ref1, Table 1. A
* related D value is found in legacy_d_qval[x]. All D values expect
* legacy_d_qval[1] are also taken from Ref1, Table 1. The case of D == 21 as
* listed in Ref1, Table 1 can be changed to D == 7 for the sake of the test
* because of {note 6}.
*
* It should be noted that the D values all satisfy the selection values
* as outlined in the gen_v1() function comments. That is:
* as outlined in the legacy_gen_v1() function comments. That is:
*
* D == P*(2^f)*(3^g)
*
@@ -571,20 +1093,20 @@ gen_u0(h, n, v1)
* where Q == 1. No further processing is needed to compute v(1) when r
* is of this form.
*/
quickmax = 8;
mat d_qval[quickmax];
mat v1_qval[quickmax];
d_qval[0] = 5; v1_qval[0] = 3; /* a=1 b=1 r=4 */
d_qval[1] = 7; v1_qval[1] = 5; /* a=3 b=1 r=12 D=21 */
d_qval[2] = 13; v1_qval[2] = 11; /* a=3 b=1 r=4 */
d_qval[3] = 11; v1_qval[3] = 20; /* a=3 b=1 r=2 */
d_qval[4] = 29; v1_qval[4] = 27; /* a=5 b=1 r=4 */
d_qval[5] = 53; v1_qval[5] = 51; /* a=53 b=1 r=4 */
d_qval[6] = 17; v1_qval[6] = 66; /* a=17 b=1 r=1 */
d_qval[7] = 19; v1_qval[7] = 74; /* a=38 b=1 r=2 */
legacy_quickmax = 8;
mat legacy_d_qval[legacy_quickmax];
mat legacy_v1_qval[legacy_quickmax];
legacy_d_qval[0] = 5; legacy_v1_qval[0] = 3; /* a=1 b=1 r=4 */
legacy_d_qval[1] = 7; legacy_v1_qval[1] = 5; /* a=3 b=1 r=12 D=21 */
legacy_d_qval[2] = 13; legacy_v1_qval[2] = 11; /* a=3 b=1 r=4 */
legacy_d_qval[3] = 11; legacy_v1_qval[3] = 20; /* a=3 b=1 r=2 */
legacy_d_qval[4] = 29; legacy_v1_qval[4] = 27; /* a=5 b=1 r=4 */
legacy_d_qval[5] = 53; legacy_v1_qval[5] = 51; /* a=53 b=1 r=4 */
legacy_d_qval[6] = 17; legacy_v1_qval[6] = 66; /* a=17 b=1 r=1 */
legacy_d_qval[7] = 19; legacy_v1_qval[7] = 74; /* a=38 b=1 r=2 */
/*
* gen_v1 - compute the v(1) for a given h*2^n-1 if we can
* legacy_gen_v1 - compute the v(1) for a given h*2^n-1 if we can
*
* This function assumes:
*
@@ -613,12 +1135,12 @@ d_qval[7] = 19; v1_qval[7] = 74; /* a=38 b=1 r=2 */
*
* If this case condition is true, then:
*
* u(0) = (2+sqrt(3))^h + (2-sqrt(3))^h (see Ref1, page 869)
* = (2+sqrt(3))^h + (2+sqrt(3))^(-h)
* u(2) = (2+sqrt(3))^h + (2-sqrt(3))^h (see Ref1, page 869)
* = (2+sqrt(3))^h + (2+sqrt(3))^(-h) (some call this u(0))
*
* and since Ref1, Theorem 5 states:
*
* u(0) = alpha^h + alpha^(-h)
* u(2) = alpha^h + alpha^(-h)
* r = abs(2^2 - 1^2*3) = 1
*
* and the bottom of Ref1, page 872 states:
@@ -631,7 +1153,7 @@ d_qval[7] = 19; v1_qval[7] = 74; /* a=38 b=1 r=2 */
*
* then
*
* u(0) = v(h)
* u(2) = v(h)
*
* so we simply return
*
@@ -666,7 +1188,7 @@ d_qval[7] = 19; v1_qval[7] = 74; /* a=38 b=1 r=2 */
*
* where L(x,y) is the Legendre symbol (see below), then:
*
* u(0) = alpha^h + alpha^(-h)
* u(2) = alpha^h + alpha^(-h)
*
* The bottom of Ref1, page 872 states:
*
@@ -674,7 +1196,7 @@ d_qval[7] = 19; v1_qval[7] = 74; /* a=38 b=1 r=2 */
*
* thus since:
*
* u(0) = v(h)
* u(2) = v(h)
*
* so we want to return:
*
@@ -929,7 +1451,7 @@ d_qval[7] = 19; v1_qval[7] = 74; /* a=38 b=1 r=2 */
* returns v(1), or -1 is there is no quick way
*/
define
gen_v1(h, n)
legacy_gen_v1(h, n)
{
local d; /* the 'D' value to try */
local val_mod; /* h*2^n-1 mod 'D' */
@@ -947,10 +1469,10 @@ gen_v1(h, n)
* We will try all 'D' values until we find a proper v(1)
* or run out of 'D' values.
*/
for (i=0; i < quickmax; ++i) {
for (i=0; i < legacy_quickmax; ++i) {
/* grab our 'D' value */
d = d_qval[i];
d = legacy_d_qval[i];
/* compute h*2^n-1 mod 'D' quickly */
val_mod = (h*pmod(2,n%(d-1),d)-1) % d;
@@ -965,13 +1487,13 @@ gen_v1(h, n)
/* D mod 4 == 1, so check for J(D, h*2^n-1) == -1 */
if (jacobi(val_mod, d) == -1) {
/* it worked, return the related v(1) value */
return v1_qval[i];
return legacy_v1_qval[i];
}
} else {
/* D mod 4 == -1, so check for J(D, h*2^n-1) == 1 */
if (jacobi(val_mod, d) == 1) {
/* it worked, return the related v(1) value */
return v1_qval[i];
return legacy_v1_qval[i];
}
}
}
@@ -1029,19 +1551,3 @@ gen_v1(h, n)
/* no quick and dirty v(1), so return -1 */
return -1;
}
/*
* ldebug - print a debug statement
*
* input:
* funct name of calling function
* str string to print
*/
define
ldebug(funct, str)
{
if (config("resource_debug") & 8) {
print "DEBUG:", funct:":", str;
}
return;
}

165
cal/lucas_tbl.cal
View File

@@ -1,165 +0,0 @@
/*
* lucas_tbl - lucasian criteria for primality tables
*
* Copyright (C) 1999 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.1 $
* @(#) $Id: lucas_tbl.cal,v 30.1 2007/03/16 11:09:54 chongo Exp $
* @(#) $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
*
* chongo <was here> /\oo/\ http://www.isthe.com/chongo/
* Share and enjoy! :-) http://www.isthe.com/chongo/tech/comp/calc/
*/
/*
* Lucasian criteria for primality
*
* The following table is taken from:
*
* "Lucasian Criteria for the Primality of N=h*2^n-1", by Hans Riesel,
* Mathematics of Computation, Vol 23 #108, p 872.
*
* The index of the *_val[] arrays correspond to the v(1) values found
* in the table. That is, for v(1) == x:
*
* D == d_val[x]
* a == a_val[x]
* b == b_val[x]
* r == r_val[x] (r == abs(a^2 - b^2*D))
*
*
* Note that when *_val[i] is not a number, the related v(1) value
* is not found in Table 1.
*/
trymax = 100;
mat d_val[trymax+1];
mat a_val[trymax+1];
mat b_val[trymax+1];
mat r_val[trymax+1];
/* v1= 0 INVALID */
/* v1= 1 INVALID */
/* v1= 2 INVALID */
d_val[ 3]= 5; a_val[ 3]= 1; b_val[ 3]=1; r_val[ 3]=4;
d_val[ 4]= 3; a_val[ 4]= 1; b_val[ 4]=1; r_val[ 4]=2;
d_val[ 5]= 21; a_val[ 5]= 3; b_val[ 5]=1; r_val[ 5]=12;
d_val[ 6]= 2; a_val[ 6]= 1; b_val[ 6]=1; r_val[ 6]=1;
/* v1= 7 INVALID */
d_val[ 8]= 15; a_val[ 8]= 3; b_val[ 8]=1; r_val[ 8]=6;
d_val[ 9]= 77; a_val[ 9]= 7; b_val[ 9]=1; r_val[ 9]=28;
d_val[10]= 6; a_val[10]= 2; b_val[10]=1; r_val[10]=2;
d_val[11]= 13; a_val[11]= 3; b_val[11]=1; r_val[11]=4;
d_val[12]= 35; a_val[12]= 5; b_val[12]=1; r_val[12]=10;
d_val[13]= 165; a_val[13]=11; b_val[13]=1; r_val[13]=44;
/* v1=14 INVALID */
d_val[15]= 221; a_val[15]=13; b_val[15]=1; r_val[15]=52;
d_val[16]= 7; a_val[16]= 3; b_val[16]=1; r_val[16]=2;
d_val[17]= 285; a_val[17]=15; b_val[17]=1; r_val[17]=60;
/* v1=18 INVALID */
d_val[19]= 357; a_val[19]=17; b_val[19]=1; r_val[19]=68;
d_val[20]= 11; a_val[20]= 3; b_val[20]=1; r_val[20]=2;
d_val[21]= 437; a_val[21]=19; b_val[21]=1; r_val[21]=76;
d_val[22]= 30; a_val[22]= 5; b_val[22]=1; r_val[22]=5;
/* v1=23 INVALID */
d_val[24]= 143; a_val[24]=11; b_val[24]=1; r_val[24]=22;
d_val[25]= 69; a_val[25]= 9; b_val[25]=1; r_val[25]=12;
d_val[26]= 42; a_val[26]= 6; b_val[26]=1; r_val[26]=6;
d_val[27]= 29; a_val[27]= 5; b_val[27]=1; r_val[27]=4;
d_val[28]= 195; a_val[28]=13; b_val[28]=1; r_val[28]=26;
d_val[29]= 93; a_val[29]= 9; b_val[29]=1; r_val[29]=12;
d_val[30]= 14; a_val[30]= 4; b_val[30]=1; r_val[30]=2;
d_val[31]= 957; a_val[31]=29; b_val[31]=1; r_val[31]=116;
d_val[32]= 255; a_val[32]=15; b_val[32]=1; r_val[32]=30;
d_val[33]=1085; a_val[33]=31; b_val[33]=1; r_val[33]=124;
/* v1=34 INVALID */
d_val[35]=1221; a_val[35]=33; b_val[35]=1; r_val[35]=132;
d_val[36]= 323; a_val[36]=17; b_val[36]=1; r_val[36]=34;
d_val[37]=1365; a_val[37]=35; b_val[37]=1; r_val[37]=140;
d_val[38]= 10; a_val[38]= 3; b_val[38]=1; r_val[38]=1;
d_val[39]=1517; a_val[39]=37; b_val[39]=1; r_val[39]=148;
d_val[40]= 399; a_val[40]=19; b_val[40]=1; r_val[40]=38;
d_val[41]=1677; a_val[41]=39; b_val[41]=1; r_val[41]=156;
d_val[42]= 110; a_val[42]=10; b_val[42]=1; r_val[42]=10;
d_val[43]= 205; a_val[43]=15; b_val[43]=1; r_val[43]=20;
d_val[44]= 483; a_val[44]=21; b_val[44]=1; r_val[44]=42;
d_val[45]=2021; a_val[45]=43; b_val[45]=1; r_val[45]=172;
d_val[46]= 33; a_val[46]= 6; b_val[46]=1; r_val[46]=3;
/* v1=47 INVALID */
d_val[48]= 23; a_val[48]= 5; b_val[48]=1; r_val[48]=2;
d_val[49]=2397; a_val[49]=47; b_val[49]=1; r_val[49]=188;
d_val[50]= 39; a_val[50]= 6; b_val[50]=1; r_val[50]=3;
d_val[51]= 53; a_val[51]= 7; b_val[51]=1; r_val[51]=4;
/* v1=52 INVALID */
d_val[53]=2805; a_val[53]=51; b_val[53]=1; r_val[53]=204;
d_val[54]= 182; a_val[54]=13; b_val[54]=1; r_val[54]=13;
d_val[55]=3021; a_val[55]=53; b_val[55]=1; r_val[55]=212;
d_val[56]= 87; a_val[56]= 9; b_val[56]=1; r_val[56]=6;
d_val[57]=3245; a_val[57]=55; b_val[57]=1; r_val[57]=220;
d_val[58]= 210; a_val[58]=14; b_val[58]=1; r_val[58]=14;
d_val[59]=3477; a_val[59]=57; b_val[59]=1; r_val[59]=228;
d_val[60]= 899; a_val[60]=29; b_val[60]=1; r_val[60]=58;
d_val[61]= 413; a_val[61]=21; b_val[61]=1; r_val[61]=28;
/* v1=62 INVALID */
d_val[63]=3965; a_val[63]=61; b_val[63]=1; r_val[63]=244;
d_val[64]=1023; a_val[64]=31; b_val[64]=1; r_val[64]=62;
d_val[65]= 469; a_val[65]=21; b_val[65]=1; r_val[65]=28;
d_val[66]= 17; a_val[66]= 4; b_val[66]=1; r_val[66]=1;
d_val[67]=4485; a_val[67]=65; b_val[67]=1; r_val[67]=260;
d_val[68]=1155; a_val[68]=33; b_val[68]=1; r_val[68]=66;
d_val[69]=4757; a_val[69]=67; b_val[69]=1; r_val[69]=268;
d_val[70]= 34; a_val[70]= 6; b_val[70]=1; r_val[70]=2;
d_val[71]=5037; a_val[71]=69; b_val[71]=1; r_val[71]=276;
d_val[72]=1295; a_val[72]=35; b_val[72]=1; r_val[72]=70;
d_val[73]= 213; a_val[73]=15; b_val[73]=1; r_val[73]=12;
d_val[74]= 38; a_val[74]= 6; b_val[74]=1; r_val[74]=2;
d_val[75]=5621; a_val[75]=73; b_val[75]=1; r_val[75]=292;
d_val[76]=1443; a_val[76]=37; b_val[76]=1; r_val[76]=74;
d_val[77]= 237; a_val[77]=15; b_val[77]=1; r_val[77]=12;
d_val[78]= 95; a_val[78]=10; b_val[78]=1; r_val[78]=5;
/* v1=79 INVALID */
d_val[80]=1599; a_val[80]=39; b_val[80]=1; r_val[80]=78;
d_val[81]=6557; a_val[81]=79; b_val[81]=1; r_val[81]=316;
d_val[82]= 105; a_val[82]=10; b_val[82]=1; r_val[82]=5;
d_val[83]= 85; a_val[83]= 9; b_val[83]=1; r_val[83]=4;
d_val[84]=1763; a_val[84]=41; b_val[84]=1; r_val[84]=82;
d_val[85]=7221; a_val[85]=83; b_val[85]=1; r_val[85]=332;
d_val[86]= 462; a_val[86]=21; b_val[86]=1; r_val[86]=21;
d_val[87]=7565; a_val[87]=85; b_val[87]=1; r_val[87]=340;
d_val[88]= 215; a_val[88]=15; b_val[88]=1; r_val[88]=10;
d_val[89]=7917; a_val[89]=87; b_val[89]=1; r_val[89]=348;
d_val[90]= 506; a_val[90]=22; b_val[90]=1; r_val[90]=22;
d_val[91]=8277; a_val[91]=89; b_val[91]=1; r_val[91]=356;
d_val[92]= 235; a_val[92]=15; b_val[92]=1; r_val[92]=10;
d_val[93]=8645; a_val[93]=91; b_val[93]=1; r_val[93]=364;
d_val[94]= 138; a_val[94]=12; b_val[94]=1; r_val[94]=6;
d_val[95]=9021; a_val[95]=93; b_val[95]=1; r_val[95]=372;
d_val[96]= 47; a_val[96]= 7; b_val[96]=1; r_val[96]=2;
d_val[97]=1045; a_val[97]=33; b_val[97]=1; r_val[97]=44;
/* v1=98 INVALID */
d_val[99]=9797; a_val[99]=97; b_val[99]=1; r_val[99]=388;
d_val[100]= 51; a_val[100]= 7; b_val[100]=1; r_val[100]=2;
if (config("resource_debug") & 3) {
print "d_val[100] defined";
print "a_val[100] defined";
print "b_val[100] defined";
print "r_val[100] defined";
}

11
cal/regress.cal
View File

@@ -1,7 +1,7 @@
/*
* regress - calc regression and correctness test suite
*
* Copyright (C) 1999-2006 David I. Bell and Landon Curt Noll
* Copyright (C) 1999-2017 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
@@ -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.13 $
* @(#) $Id: regress.cal,v 30.13 2016/02/06 08:39:35 chongo Exp $
* @(#) $Source: /usr/local/src/bin/calc/cal/RCS/regress.cal,v $
* @(#) $Revision: 30.14 $
* @(#) $Id: regress.cal,v 30.14 2017/05/19 16:09:14 chongo Exp $
* @(#) $Source: /usr/local/src/bin/calc-RHEL7/cal/RCS/regress.cal,v $
*
* Under source code control: 1990/02/15 01:50:36
* File existed as early as: before 1990
@@ -8068,7 +8068,6 @@ return test_functions2();
* hello.cal - designed to go into an infinite loop
* lucal.cal - already read by this file
* lucas_chk.cal - already read by this file
* lucas_tbl.cal - duplicatres code already read by another file
* regress.cal - this file
* surd.cal - already read by this file
* test9999.cal - files of this form are already read by this file
@@ -8141,8 +8140,6 @@ read -once intfile;
print '9830: read -once intfile';
read -once lucas;
print '9831: read -once lucas';
read -once lucas_tbl;
print '9832: read -once lucas_tbl';
read -once natnumset;
print '9833: read -once natnumset';
read -once repeat;