Improve formatting of comments in cal/lucas.cal

Fixed a man page warning about ./myfile where the leading dot
was mistook for an nroff macro.  Thanks goes to David Haller
<dnh at opensuse dot org> for providing the patch.

Improved gen_v1(h,n) in lucas.cal for cases where h is not a
multiple of 3. Optimized the search for v(1) when h is a
multiple of 3.

Fixed a Makefile problem, reported by Doug Hays <doughays6 at gmail
dot com>, where if a macOS user set BINDIR, LIBDIR, CALC_SHAREDIR
or INCDIR in the top section, their values will be overwritten by
the Darwin specific section.

CHANGES

@@ -1,4 +1,29 @@
-The following are the changes from calc version 2.12.6.0 to date:
+The following are the changes from calc version 2.12.6.1 to date:
+
+    Improved gen_v1(h,n) in lucas.cal to use an even faster search method.
+
+    Improved are checking in lucas.cal.  In particular both h and n must be
+    integers >= 1.  In the case of both rodseth_xhn(x, h, n) and gen_v1(h, n)
+    h must be odd.
+
+    Fixed an C code indenting issue that was reported by Thomas Walter
+    <th dot walter42 at gmx dot de> in zfunc.c.
+
+    Fixed a man page warning about ./myfile where the leading dot
+    was mistook for an nroff macro.  Thanks goes to David Haller
+    <dnh at opensuse dot org> for providing the patch.
+
+    Improved gen_v1(h,n) in lucas.cal for cases where h is not a
+    multiple of 3. Optimized the search for v(1) when h is a
+    multiple of 3.
+
+    Fixed a Makefile problem, reported by Doug Hays <doughays6 at gmail
+    dot com>, where if a macOS user set BINDIR, LIBDIR, CALC_SHAREDIR
+    or INCDIR in the top section, their values will be overwritten by
+    the Darwin specific section.
+
+
+The following are the changes from calc version 2.12.6.0 to 2.12.6.0:
 
     Added the makefile variable ${COMMON_ADD} that will add flags
     to all compile and link commands. The ${COMMON_ADD} flags are
@@ -68,7 +93,6 @@ The following are the changes from calc version 2.12.6.0 to date:
 
     Fixed a number of typos in the CHANGES file.
 
-
 The following are the changes from calc version 2.12.5.4 to 2.12.5.6:
 
     Recompile to match current RHEL7.2 libc and friends.

Makefile.ship

@@ -23,7 +23,7 @@
 #	 READLINE_LIB= -lreadline
 #	 READLINE_EXTRAS= -lhistory -lncurses
 #
-# Copyright (C) 1999-2017  Landon Curt Noll
+# Copyright (C) 1999-2018  Landon Curt Noll
 #
 # Calc is open software; you can redistribute it and/or modify it under
 # the terms of version 2.1 of the GNU Lesser General Public License
@@ -564,14 +564,30 @@ HAVE_UNUSED=
 #
 #	INCDIR= /dev/env/DJDIR/include
 #
-# If in doubt, set:
+# If in doubt, for non-macOS hosts set:
 #
 #	INCDIR= /usr/include
 #
+# However, if you are on macOS then set:
+#
+#	INCDIR= /usr/local/include
+#if 0	/* start of skip for non-Gnu makefiles */
+ifeq ($(target),Darwin)
 
+# default INCDIR for macOS
+INCDIR= /usr/local/include
+
+else
+#endif	/* end of skip for non-Gnu makefiles */
+
+# default INCDIR for non-macOS
+INCDIR= /usr/include
 #INCDIR= /usr/local/include
 #INCDIR= /dev/env/DJDIR/include
-INCDIR= /usr/include
+
+#if 0	/* start of skip for non-Gnu makefiles */
+endif
+#endif	/* end of skip for non-Gnu makefiles */
 
 # Where to install calc related things
 #
@@ -602,23 +618,71 @@ INCDIR= /usr/include
 #	LIBDIR= /dev/env/DJDIR/lib
 #	CALC_SHAREDIR= /dev/env/DJDIR/share/calc
 #
-# If in doubt, set:
+# If in doubt, for non-macOS hosts set:
 #
 #	BINDIR= /usr/bin
 #	LIBDIR= /usr/lib
 #	CALC_SHAREDIR= /usr/share/calc
 #
+# However, if you are on macOS then set:
+#
+#	BINDIR= /usr/local/bin
+#	LIBDIR= /usr/local/lib
+#	CALC_SHAREDIR= /usr/local/share/calc
+#
+# 	NOTE: Starting with macOS El Capitan OS X 10.11, root by default
+# 	      could not mkdir under system locations, so macOS must now
+# 	      use the /usr/local tree.
+
+#if 0	/* start of skip for non-Gnu makefiles */
+ifeq ($(target),Darwin)
+
+# default BINDIR for macOS
+BINDIR= /usr/local/bin
+
+else
+#endif	/* end of skip for non-Gnu makefiles */
+
+# default BINDIR for non-macOS
+BINDIR= /usr/bin
 #BINDIR= /usr/local/bin
 #BINDIR= /dev/env/DJDIR/bin
-BINDIR= /usr/bin
 
+#if 0	/* start of skip for non-Gnu makefiles */
+endif
+
+ifeq ($(target),Darwin)
+
+# default LIBDIR for macOS
+LIBDIR= /usr/local/lib
+
+else
+#endif	/* end of skip for non-Gnu makefiles */
+
+# default LIBDIR for non-macOS
+LIBDIR= /usr/lib
 #LIBDIR= /usr/local/lib
 #LIBDIR= /dev/env/DJDIR/lib
-LIBDIR= /usr/lib
 
+#if 0	/* start of skip for non-Gnu makefiles */
+endif
+
+ifeq ($(target),Darwin)
+
+# default CALC_SHAREDIR for macOS
+CALC_SHAREDIR= /usr/local/share/calc
+
+else
+#endif	/* end of skip for non-Gnu makefiles */
+
+# default CALC_SHAREDIR for non-macOS
+CALC_SHAREDIR= /usr/share/calc
 #CALC_SHAREDIR= /usr/local/lib/calc
 #CALC_SHAREDIR= /dev/env/DJDIR/share/calc
-CALC_SHAREDIR= /usr/share/calc
+
+#if 0	/* start of skip for non-Gnu makefiles */
+endif
+#endif	/* end of skip for non-Gnu makefiles */
 
 # NOTE: Do not set CALC_INCDIR to /usr/include or /usr/local/include!!!
 #	Always be sure that the CALC_INCDIR path ends in /calc to avoid
@@ -990,7 +1054,7 @@ EXT=
 
 # The default calc versions
 #
-VERSION= 2.12.6.0
+VERSION= 2.12.6.4
 
 # Names of shared libraries with versions
 #
@@ -1263,16 +1327,6 @@ LDCONFIG:=
 # DARWIN_ARCH= -arch ppc		# PPC binary
 # DARWIN_ARCH= -arch x86_64		# native 64-bit binary
 DARWIN_ARCH=				# native binary
-#
-# Starting with El Capitan OS X 10.11, root by default could not
-# mkdir under system locations, so we now use the /usr/local tree.
-#
-OPTDIR:= /usr/local
-BINDIR:= /${OPTDIR}/bin
-LIBDIR:= /${OPTDIR}/lib
-CALC_SHAREDIR:= /${OPTDIR}/share
-CALC_INCDIR:= /${OPTDIR}/include
-SCRIPTDIR:= ${BINDIR}/cscript
 endif
 
 ##################

cal/lucas.cal

@@ -28,6 +28,10 @@
  * 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
+ *
+ * Also see the reference code, available both C and go:
+ *
+ *	https://github.com/arcetri/goprime
  */
 
 /*
@@ -154,6 +158,12 @@
  *
  * Finally, one should eliminate all values of 'h*2^n-1' where
  * 'h*2^n+1' is divisible by a small primes.
+ *
+ * NOTE: Today, for world record sized h*2^n-1 primes, one might
+ *	 search for factors < 2^46 or more.  By excluding h*2^n-1
+ *	 with prime factors < 2^46, where h*2^n-1 is a bit larger
+ *	 than the largest known prime, one may exclude about 96.5%
+ *	 of candidates that have "small" prime factors.
  */
 
 pprod256 = 0;	/* product of  "primes up to 256" / "primes up to 46" */
@@ -243,8 +253,8 @@ pprod256 = 0;	/* product of  "primes up to 256" / "primes up to 46" */
  * do make this so.
  *
  * input:
- *	h    the h as in h*2^n-1
- *	n    the n as in h*2^n-1
+ *	h    h as in h*2^n-1 (must be >= 1)
+ *	n    n as in h*2^n-1 (must be >= 1)
  *
  * returns:
  *	1 => h*2^n-1 is prime
@@ -267,13 +277,17 @@ lucas(h, n)
 	 * check arg types
 	 */
 	if (!isint(h)) {
-		ldebug("lucas", "h is non-int");
 		quit "FATAL: bad args: h must be an integer";
 	}
+	if (h < 1) {
+		quit "FATAL: bad args: h must be an integer >= 1";
+	}
 	if (!isint(n)) {
-		ldebug("lucas", "n is non-int");
 		quit "FATAL: bad args: n must be an integer";
 	}
+	if (n < 1) {
+		quit "FATAL: bad args: n must be an integer >= 1";
+	}
 
 	/*
 	 * reduce h if even
@@ -484,9 +498,9 @@ lucas(h, n)
  * See the function gen_v1() for details on the value of v(1).
  *
  * 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)
+ *	h	- h as in h*2^n-1	(must be >= 1)
+ *	n	- n as in h*2^n-1	(must be >= 1)
+ *	v1	- gen_v1(h,n)		(must be >= 3) (see function below)
  *
  * returns:
  *	u(2)	- initial value for Lucas test on h*2^n-1
@@ -499,6 +513,8 @@ gen_u2(h, n, v1)
 	local r;		/* low value: v(n) */
 	local s;		/* high value: v(n+1) */
 	local hbits;		/* highest bit set in h */
+	local oldh;		/* pre-reduced h */
+	local oldn;		/* pre-reduced n */
 	local i;
 
 	/*
@@ -507,24 +523,45 @@ gen_u2(h, n, v1)
 	if (!isint(h)) {
 		quit "bad args: h must be an integer";
 	}
+	if (h < 0) {
+		quit "bad args: h must be an integer >= 1";
+	}
 	if (!isint(n)) {
 		quit "bad args: n must be an integer";
 	}
+	if (n < 1) {
+		quit "bad args: n must be an integer >= 1";
+	}
 	if (!isint(v1)) {
 		quit "bad args: v1 must be an integer";
 	}
-	if (v1 <= 0) {
-		quit "bogus arg: v1 is <= 0";
+	if (v1 < 3) {
+		quit "bogus arg: v1 must be an integer >= 3";
+	}
+
+	/*
+	 * reduce h if even
+	 *
+	 * we will force h to be odd by moving powers of two over to 2^n
+	 */
+	oldh = h;
+	oldn = n;
+	shiftdown = fcnt(h,2);	/* h % 2^shiftdown == 0, max shiftdown */
+	if (shiftdown > 0) {
+		h >>= shiftdown;
+		n += shiftdown;
 	}
 
 	/*
 	 * enforce the h > 0 and n >= 2 rules
 	 */
 	if (h <= 0 || n < 1) {
+		print "	   ERROR: h=":oldh, "n=":oldn, "reduced h=":h, "n=":n;
 		quit "reduced args violate the rule: 0 < h < 2^n";
 	}
 	hbits = highbit(h);
 	if (hbits >= n) {
+		print "	   ERROR: h=":oldh, "n=":oldn, "reduced h=":h, "n=":n;
 		quit "reduced args violate the rule: 0 < h < 2^n";
 	}
 
@@ -599,8 +636,8 @@ gen_u2(h, n, v1)
  * See the function gen_u2() for details.
  *
  * input:
- *	h	- h as in h*2^n-1
- *	n	- n as in h*2^n-1
+ *	h	- h as in h*2^n-1	(must be >= 1)
+ *	n	- n as in h*2^n-1	(must be >= 1)
  *	v1	- gen_v1(h,n)		(see function below)
  *
  * returns:
@@ -638,9 +675,9 @@ gen_u0(h, n, v1)
  *	x > 2
  *
  * input:
- *	x	- potential v(1) value
- *	h	- h as in h*2^n-1
- *	n	- n as in h*2^n-1
+ *	x	potential v(1) value
+ *	h	h as in h*2^n-1 (h must be odd >= 1)
+ *	n	n as in h*2^n-1	(must be >= 1)
  *
  * returns:
  *	1	if v(1) == x for h*2^n-1
@@ -657,9 +694,18 @@ rodseth_xhn(x, h, n)
 	if (!isint(h)) {
 		quit "bad args: h must be an integer";
 	}
+	if (iseven(h)) {
+		quit "bad args: h must be an odd integer";
+	}
+	if (h < 1) {
+		quit "bad args: h must be an integer >= 1";
+	}
 	if (!isint(n)) {
 		quit "bad args: n must be an integer";
 	}
+	if (n < 1) {
+		quit "bad args: n must be an integer >= 1";
+	}
 	if (!isint(x)) {
 		quit "bad args: x must be an integer";
 	}
@@ -703,9 +749,9 @@ rodseth_xhn(x, h, n)
  * 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:
+ * of h*2^n-1, h odd 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
  *    ----------
@@ -763,45 +809,161 @@ rodseth_xhn(x, h, n)
  *           1 161
  *           1 155
  *
+ * It is important that we select the smallest possible v(1).  While testing
+ * various values of X for V(1) is fast, using larger than necessary values
+ * of V(1) of can slow down calculating V(h).
+ *
  * 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.
+ * odd h that are also a multiple of 3 and for many values of n where h < 2^n.
+ *
+ * For example for in a sample size of 1000000 numbers of the form h*2^n-1
+ * where h is an odd multiple of 3, 12996351 <= h <= 13002351,
+ * 4331116 <= n <= 4332116, these are the smallest v(1) values that were found:
+ *
+ *  smallest  percentage
+ *	v(1)     used
+ *  --------  ---------
+ *	  3   40.0000 %
+ *	  5   25.6833 %
+ *	  9   11.6924 %
+ *	 11   10.4528 %
+ *	 15    4.8048 %
+ *	 17    2.3458 %
+ *	 21    1.3734 %
+ *	 29    1.0527 %
+ *	 20    0.8595 %
+ *	 27    0.5758 %
+ *	 35    0.4420 %
+ *	 36    0.2433 %
+ *	 39    0.1779 %
+ *	 41    0.0885 %
+ *	 45    0.0571 %
+ *	 32    0.0337 %
+ *	 51    0.0289 %
+ *	 44    0.0205 %
+ *	 49    0.0176 %
+ *	 56    0.0137 %
+ *	 59    0.0108 %
+ *	 57    0.0053 %
+ *	 65    0.0047 %
+ *	 55    0.0045 %
+ *	 69    0.0031 %
+ *	 71    0.0024 %
+ *	 66    0.0011 %
+ *	 95    0.0008 %
+ *	 81    0.0008 %
+ *	 77    0.0006 %
+ *	 72    0.0005 %
+ *	 99    0.0004 %
+ *	 80    0.0003 %
+ *	 74    0.0003 %
+ *	 84    0.0002 %
+ *	 67    0.0002 %
+ *	 87    0.0001 %
+ *	104    0.0001 %
+ *	129    0.0001 %
+ *
+ * However, a case can be made for considering only odd values for v(1)
+ * candidates.  When h * 2^n-1 is prime and h is an odd multiple of 3,
+ * a smallest v(1) that is even is extremely rate.  Of the list of 146553
+ * known primes of the form h*2^n-1 when h is an odd a multiple of 3,
+ * none has an smallest v(1) that was even.
+ *
+ * See:
+ *
+ *	https://github.com/arcetri/verified-prime
+ *
+ * for that list of 146553 known primes of the form h*2^n-1.
+ *
+ * That same example for in a sample size of 1000000 numbers of the
+ * form h*2^n-1 where h is an odd multiple of 3, 12996351 <= h <= 13002351,
+ * 4331116 <= n <= 4332116, these are the smallest odd v(1) values that were
+ * found:
+ *
+ *  smallest    percentage
+ *  odd v(1)       used
+ *  --------    ---------
+ *	  3	40.0000 %
+ *	  5	25.6833 %
+ *	  9	11.6924 %
+ *	 11	10.4528 %
+ *	 15	 4.8048 %
+ *	 17	 2.3458 %
+ *	 21	 1.6568 %
+ *	 29	 1.6174 %
+ *	 35	 0.4529 %
+ *	 27	 0.3546 %
+ *	 39	 0.3470 %
+ *	 41	 0.2159 %
+ *	 45	 0.1173 %
+ *	 31	 0.0661 %
+ *	 51	 0.0619 %
+ *	 55	 0.0419 %
+ *	 59	 0.0250 %
+ *	 49	 0.0170 %
+ *	 69	 0.0110 %
+ *	 65	 0.0098 %
+ *	 71	 0.0078 %
+ *	 85	 0.0048 %
+ *	 81	 0.0044 %
+ *	 95	 0.0038 %
+ *	 99	 0.0021 %
+ *	125	 0.0009 %
+ *	 57	 0.0007 %
+ *	111	 0.0005 %
+ *	 77	 0.0003 %
+ *	165	 0.0003 %
+ *	155	 0.0002 %
+ *	129	 0.0002 %
+ *	101	 0.0002 %
+ *	 53	 0.0001 %
+ *
+ * Moreover when evaluating odd candidates for v(1), one may cache Jacobi
+ * symbol evaluations to reduce the number of Jacobi symbol evaluations to
+ * a minimum. For example, if one tests 5 and finds that the 2nd case fails:
+ *
+ *	jacobi(5+2, h*2^n-1) != -1
+ *
+ * Then if one is later testing 9, the Jacobi symbol value for the first
+ * 1st case:
+ *
+ *	jacobi(7-2, h*2^n-1)
+ *
+ * is already known.
+ *
+ * Without Jacobi symbol value caching, it requires on average
+ * 4.851377 Jacobi symbol evaluations.  With Jacobi symbol value caching
+ * cacheing, an averare of 4.348820 Jacobi symbol evaluations is needed.
  *
  * Given this information, when odd h is a multiple of 3 we try, in order,
- * these values of X:
+ * these odd 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
+ *	3, 5, 9, 11, 15, 17, 21, 29, 27, 35, 39, 41, 31, 45, 51, 55, 49, 59,
+ *	69, 65, 71, 57, 85, 81, 95, 99, 77, 53, 67, 125, 111, 105, 87, 129,
+ *	101, 83, 165, 155, 149, 141, 121, 109
  *
  * 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.
+ * Less than 1 case out of 1000000 will not be satisifed by the above list.
+ * If no value in that list works, we start simple search starting with X = 167
+ * and incrementing by 2 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.
+ * The x_tbl[] matrix contains those values of X to try in order.
+ * If all x_tbl_len fail to satisfy Ref4 condition 1 (this happens less than
+ * 1 in 1000000 cases), then we begin a linear search of odd values starting at
+ * next_x until we find a proper X value.
  */
-x_tbl_len = 45;
+x_tbl_len = 42;
 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
+    3, 5, 9, 11, 15, 17, 21, 29, 27, 35, 39, 41, 31, 45, 51, 55, 49, 59,
+    69, 65, 71, 57, 85, 81, 95, 99, 77, 53, 67, 125, 111, 105, 87, 129,
+    101, 83, 165, 155, 149, 141, 121, 109
 };
-next_x = 120;
+next_x = 167;	/* must be 2 more than the largest value in x_tbl[] */
 
 /*
  * gen_v1 - compute the v(1) for a given h*2^n-1 if we can
@@ -859,12 +1021,22 @@ next_x = 120;
  *
  *	u(2) = v(h)				 (NOTE: some call this u(2))
  *
- * so we simply return
+ * so we can always return
  *
  *	v(1) = alpha^1 + alpha^(-1)
  *	     = (2+sqrt(3)) + (2-sqrt(3))
  *	     = 4
  *
+ * In 40% of the cases when h is not a multiple of 3, 3 is a valid value
+ * for v(1).  We can test if 3 is a valid value for v(1) in this case:
+ *
+ *	if jacobi(1, h*2^n-1) == 1 and jacobi(5, h*2^n-1) == -1, then
+ *	    v(1) = 3
+ *	else
+ *	    v(1) = 4
+ *
+ * HOTE: The above "if then else" works only of h is not a multiple of 3.
+ *
  ***
  *
  * Case 2:	(h mod 3 == 0)
@@ -956,17 +1128,22 @@ next_x = 120;
  ***
  *
  * input:
- *	h	h as in h*2^n-1
- *	n	n as in h*2^n-1
+ *	h	h as in h*2^n-1 (h must be odd >= 1)
+ *	n	n as in h*2^n-1	(must be >= 1)
  *
  * output:
- *	returns v(1), or -1 is there is no quick way
+ *	returns v(1), or
+ *		-1 when h*2^n-1 is a multiple of 3
  */
 define
 gen_v1(h, n)
 {
 	local x;	/* potential v(1) to test */
 	local i;	/* x_tbl index */
+	local v1m2;	/* X-2 1st case */
+	local v1p2;	/* X+2 2nd case */
+	local testval;	/* h*2^n-1 - value we are testing if prime */
+	local mat cached_v1[next_x];	/* cached Jacobi symbol values or 0 */
 
 	/*
 	 * check arg types
@@ -974,58 +1151,111 @@ gen_v1(h, n)
 	if (!isint(h)) {
 		quit "bad args: h must be an integer";
 	}
+	if (iseven(h)) {
+		quit "bad args: h must be an odd integer";
+	}
+	if (h < 1) {
+		quit "bad args: h must be an integer >= 1";
+	}
 	if (!isint(n)) {
 		quit "bad args: n must be an integer";
 	}
+	if (n < 1) {
+		quit "bad args: n must be an integer >= 1";
+	}
+
+	/*
+	 * pretest: Verify that h*2^n-1 is not a multiple of 3
+	 */
+	if (((h % 3 == 1) && (n % 2 == 0)) || ((h % 3 == 2) && (n % 2 == 1))) {
+		/* no need to test h*2^n-1, it is not prime */
+		return -1;
+	}
 
 	/*
 	 * check for Case 1:      (h mod 3 != 0)
 	 */
 	if (h % 3 != 0) {
-		/* v(1) is easy to compute */
-		return 4;
+		if (rodseth_xhn(3, h, n) == 1) {
+			/* 40% of the time, 3 works when h mod 3 != 0 */
+			return 3;
+		} else {
+			/* otherwise 4 always works when h mod 3 != 0 */
+			return 4;
+		}
 	}
 
 	/*
 	 * What follow is Case 2:      (h mod 3 == 0)
 	 */
 
+	/*
+	 * clear cache
+	 */
+	matfill(cached_v1, 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
+	 *
+	 * NOTE: If we wanted to be super optimial, we would cache
+	 *	 jacobi(X+2, h*2^n-1) that that when we increment X
+	 *	 to the next odd value, the now jacobi(X-2, h*2^n-1)
+	 *	 does not need to be re-evaluted.
 	 */
+	testval = h*2^n-1;
 	for (i=0; i < x_tbl_len; ++i) {
 
 		/*
-		 * test Ref4 condition 1:
+		 * obtain the next test candidate
 		 */
 		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;
+		/*
+		 * Check x for condition 1 part 1
+		 *
+		 *	jacobi(x-2, h*2^n-1) == 1
+		 */
+		v1m2 = x-2;
+		if (cached_v1[v1m2] == 0) {
+			cached_v1[v1m2] = jacobi(v1m2, testval);
 		}
+		if (cached_v1[v1m2] != 1) {
+			continue;
+		}
+
+		/*
+		 * Check x for condition 1 part 2
+		 *
+		 *	jacobi(x+2, h*2^n-1) == -1
+		 */
+		v1p2 = x+2;
+		if (cached_v1[v1p2] == 0) {
+			cached_v1[v1p2] = jacobi(v1p2, testval);
+		}
+		if (cached_v1[v1p2] != -1) {
+			continue;
+		}
+
+		/*
+		 * 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
+	 * We are in that rare case (less than 1 in 1 000 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.
+	 * of odd vules at next_x from here on.
 	 */
 	x = next_x;
 	while (rodseth_xhn(x, h, n) != 1) {
-		++x;
+		x += 2;
 	}
 	/* finally found a v(1) value */
 	ldebug("gen_v1", "h= " + str(h) + " n= " + str(n) +
@@ -1457,8 +1687,8 @@ legacy_d_qval[7] = 19;	legacy_v1_qval[7] = 74;	/* a=38	 b=1  r=2  */
  ***
  *
  * input:
- *	h	h as in h*2^n-1
- *	n	n as in h*2^n-1
+ *	h	h as in h*2^n-1	(must be >= 1)
+ *	n	n as in h*2^n-1	(must be >= 1)
  *
  * output:
  *	returns v(1), or -1 is there is no quick way
@@ -1470,6 +1700,22 @@ legacy_gen_v1(h, n)
 	local val_mod;	/* h*2^n-1 mod 'D' */
 	local i;
 
+	/*
+	 * check arg types
+	 */
+	if (!isint(h)) {
+		quit "bad args: h must be an integer";
+	}
+	if (h < 1) {
+		quit "bad args: h must be an integer >= 1";
+	}
+	if (!isint(n)) {
+		quit "bad args: n must be an integer";
+	}
+	if (n < 1) {
+		quit "bad args: n must be an integer >= 1";
+	}
+
 	/*
 	 * check for case 1
 	 */

calc.man

@@ -657,8 +657,8 @@ searches in succession:
 .sp 1
 .in +5n
 .nf
-./myfile
-./myfile.cal
+\a./myfile
+\a./myfile.cal
 ${LIBDIR}/myfile
 ${LIBDIR}/myfile.cal
 ${CUSTOMCALDIR}/myfile

custom/Makefile

@@ -348,7 +348,7 @@ EXT=
 
 # The default calc versions
 #
-VERSION= 2.12.6.0
+VERSION= 2.12.6.4
 
 # Names of shared libraries with versions
 #
@@ -628,16 +628,6 @@ LDCONFIG:=
 # DARWIN_ARCH= -arch ppc		# PPC binary
 # DARWIN_ARCH= -arch x86_64		# native 64-bit binary
 DARWIN_ARCH=				# native binary
-#
-# Starting with El Capitan OS X 10.11, root by default could not
-# mkdir under system locations, so we now use the /usr/local tree.
-#
-OPTDIR:= /usr/local
-BINDIR:= /${OPTDIR}/bin
-LIBDIR:= /${OPTDIR}/lib
-CALC_SHAREDIR:= /${OPTDIR}/share
-CALC_INCDIR:= /${OPTDIR}/include
-SCRIPTDIR:= ${BINDIR}/cscript
 endif
 
 ##################

custom/Makefile.head

@@ -348,7 +348,7 @@ EXT=
 
 # The default calc versions
 #
-VERSION= 2.12.6.0
+VERSION= 2.12.6.4
 
 # Names of shared libraries with versions
 #

version.c

@@ -45,7 +45,7 @@ static char *program;
 #define MAJOR_VER	2	/* major library version */
 #define MINOR_VER	12	/* minor library version */
 #define MAJOR_PATCH	6	/* major software level under library version */
-#define MINOR_PATCH	0	/* minor software level or 0 if not patched */
+#define MINOR_PATCH	4	/* minor software level or 0 if not patched */
 
 
 /*

zfunc.c

@@ -1029,93 +1029,98 @@ zgcd(ZVALUE z1, ZVALUE z2, ZVALUE *res)
 				needw = FALSE;
 			}
 			g = *a0 * w;
-			if (h < BASEB) g &= (1 << h) - 1;
-			else g &= BASE1;
+			if (h < BASEB) {
+			    g &= (1 << h) - 1;
+			} else {
+			    g &= BASE1;
+			}
+		} else {
+			g = 1;
 		}
-		else g = 1;
-	   } else
+	    } else {
 		g = (HALF) *a0 * w;
-		a = a0;
-		b = b0;
-		i = n;
-		if (g > 1) {			/* a - g * b case */
-			f = 0;
-			while (i--) {
-				f = (FULL) *a - g * *b++ - f;
-				*a++ = (HALF) f;
-				f >>= BASEB;
-				f = -f & BASE1;
-			}
-			if (f) {
-				i = m - n;
-				while (i-- && f) {
-					f = *a - f;
-					*a++ = (HALF) f;
-					f >>= BASEB;
-					f = -f & BASE1;
-				}
-			}
-			while (m && !*a0) { /* Removing trailing zeros */
-				m--;
-				a0++;
-			}
-			if (f) {		/* a - g * b < 0 */
-				while (m > 1 && a0[m-1] == BASE1) m--;
-				*a0 = - *a0;
-				a = a0;
-				i = m;
-				while (--i) {
-					a++;
-					*a = ~*a;
-				}
-			}
-		} else {				/* abs(a - b) case */
-			while (i && *a++ == *b++) i--;
-			q = n - i;
-			if (m == n) {		/* a and b same length */
-				if (i) {	 /* a not equal to b */
-					while (m && a0[m-1] == b0[m-1]) m--;
-					if (a0[m-1] < b0[m-1]) {
-						/* Swapping since a < b */
-						a = a0;
-						a0 = b0;
-						b0 = a;
-						k = j;
-					}
-					a = a0 + q;
-					b = b0 + q;
-					i = m - q;
-					f = 0;
-					while (i--) {
-						f = (FULL) *a - *b++ - f;
-						*a++ = (HALF) f;
-						f >>= BASEB;
-						f = -f & BASE1;
-					}
-				}
-			} else {		/* a has more digits than b */
-				a = a0 + q;
-				b = b0 + q;
-				i = n - q;
-				f = 0;
-				while (i--) {
-					f = (FULL) *a - *b++ - f;
-					*a++ = (HALF) f;
-					f >>= BASEB;
-					f = -f & BASE1;
-				}
-				if (f) { while (!*a) *a++ = BASE1;
-					(*a)--;
-				}
-			}
-			a0 += q;
-			m -= q;
-			while (m && !*a0) { /* Removing trailing zeros */
-				m--;
-				a0++;
-			}
-		}
-		while (m && !a0[m-1]) m--;	/* Removing leading zeros */
+	    }
+	    a = a0;
+	    b = b0;
+	    i = n;
+	    if (g > 1) {			/* a - g * b case */
+		    f = 0;
+		    while (i--) {
+			    f = (FULL) *a - g * *b++ - f;
+			    *a++ = (HALF) f;
+			    f >>= BASEB;
+			    f = -f & BASE1;
+		    }
+		    if (f) {
+			    i = m - n;
+			    while (i-- && f) {
+				    f = *a - f;
+				    *a++ = (HALF) f;
+				    f >>= BASEB;
+				    f = -f & BASE1;
+			    }
+		    }
+		    while (m && !*a0) { /* Removing trailing zeros */
+			    m--;
+			    a0++;
+		    }
+		    if (f) {		/* a - g * b < 0 */
+			    while (m > 1 && a0[m-1] == BASE1) m--;
+			    *a0 = - *a0;
+			    a = a0;
+			    i = m;
+			    while (--i) {
+				    a++;
+				    *a = ~*a;
+			    }
+		    }
+	    } else {				/* abs(a - b) case */
+		    while (i && *a++ == *b++) i--;
+		    q = n - i;
+		    if (m == n) {		/* a and b same length */
+			    if (i) {	 /* a not equal to b */
+				    while (m && a0[m-1] == b0[m-1]) m--;
+				    if (a0[m-1] < b0[m-1]) {
+					    /* Swapping since a < b */
+					    a = a0;
+					    a0 = b0;
+					    b0 = a;
+					    k = j;
+				    }
+				    a = a0 + q;
+				    b = b0 + q;
+				    i = m - q;
+				    f = 0;
+				    while (i--) {
+					    f = (FULL) *a - *b++ - f;
+					    *a++ = (HALF) f;
+					    f >>= BASEB;
+					    f = -f & BASE1;
+				    }
+			    }
+		    } else {		/* a has more digits than b */
+			    a = a0 + q;
+			    b = b0 + q;
+			    i = n - q;
+			    f = 0;
+			    while (i--) {
+				    f = (FULL) *a - *b++ - f;
+				    *a++ = (HALF) f;
+				    f >>= BASEB;
+				    f = -f & BASE1;
+			    }
+			    if (f) { while (!*a) *a++ = BASE1;
+				    (*a)--;
+			    }
+		    }
+		    a0 += q;
+		    m -= q;
+		    while (m && !*a0) { /* Removing trailing zeros */
+			    m--;
+			    a0++;
+		    }
+	    }
+	    while (m && !a0[m-1]) m--;	/* Removing leading zeros */
 	}
 	if (m == 1) {				/* a has one digit */
 		v = *a0;