mirror of
https://github.com/lcn2/calc.git
synced 2025-08-19 01:13:27 +03:00
Compare commits
19 Commits
Author | SHA1 | Date | |
---|---|---|---|
|
e555a718c0 | ||
|
b29fcf2dd5 | ||
|
4f86703843 | ||
|
07d8bf0f3e | ||
|
b4cd692bae | ||
|
83c898cc2b | ||
|
c585d7aa78 | ||
|
f42a003d04 | ||
|
8da0471f07 | ||
|
1c20261b93 | ||
|
aeb9a9d473 | ||
|
66883b390d | ||
|
ea533659ce | ||
|
9e81971f25 | ||
|
cbbd866535 | ||
|
bf23f82c29 | ||
|
ec5c584785 | ||
|
6bbb8c0e42 | ||
|
438554b0ed |
58
CHANGES
58
CHANGES
@@ -1,4 +1,44 @@
|
||||
The following are the changes from calc version 2.12.6.1 to date:
|
||||
The following are the changes from calc version 2.12.6.6 to date:
|
||||
|
||||
For historical purposes, in lucas.cal, gen_v1(1, n) always returns 4.
|
||||
|
||||
Fixed some compiler warnings, thanks to a report by Mike
|
||||
<michael dot d dot ince at gmail dot com>.
|
||||
|
||||
Added work around for a gcc warning bug, thanks to a report by Mike
|
||||
<michael dot d dot ince at gmail dot com>.
|
||||
|
||||
Fixed errors in various help files such as:
|
||||
|
||||
mat randbit seed srandom types
|
||||
|
||||
Removed the MAXSTRING symbol because it was no longer used by calc.
|
||||
|
||||
Increased HIST_SIZE (depth of the history stack) from 10k to 32k.
|
||||
|
||||
Increased TTYSIZE (reallocation size for terminal buffers) from 100 to 8191.
|
||||
|
||||
Increased MAXDEPTH (maximum depth of input stack) from 10 to 255.
|
||||
|
||||
Increased interactive input buffer size from 1024 to 256k. This has the
|
||||
effect of increasing the maximum length of an input line from a tty.
|
||||
This helps with an interactive bug that was reported by Ruslan Kabatsayev
|
||||
(b7 dot 10110111 at gmail dot com).
|
||||
|
||||
|
||||
The following are the changes from calc version 2.12.6.4 to 2.12.6.5:
|
||||
|
||||
Fixed warning about undefined operations involving the qlink(q)
|
||||
macro by replacing that macro with an inline-function. Thanks goes
|
||||
to David Haller <dnh at opensuse dot org> for this fix.
|
||||
|
||||
NOTE for Windows 10 users: Pavel Nemec <pane at seznam dot cz>
|
||||
reported that calc version 2.12.6.4 has been successfully
|
||||
compiled, installed and running on Windows 10. See README.WINDOWS
|
||||
for more details.
|
||||
|
||||
|
||||
The following are the changes from calc version 2.12.6.1 to 2.12.6.3:
|
||||
|
||||
Improved gen_v1(h,n) in lucas.cal to use an even faster search method.
|
||||
|
||||
@@ -6,6 +46,22 @@ The following are the changes from calc version 2.12.6.1 to date:
|
||||
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:
|
||||
|
||||
|
@@ -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.1
|
||||
VERSION= 2.12.6.7
|
||||
|
||||
# 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
|
||||
|
||||
##################
|
||||
|
@@ -10,6 +10,21 @@ NOTE: The main developers do not have access to a Windoz based platform.
|
||||
Of course you are welcome to send us any patches that fix your
|
||||
Windoz build environment.
|
||||
|
||||
|
||||
=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
|
||||
=-= compiling with Windows Subsystem for Linux (WSL) =-Cygwin =-=
|
||||
=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
|
||||
|
||||
It has been reported that calc version 2.12.6.4 has been successfully
|
||||
compiled, installed and running on Windows 10 on 2018 Jan 21.
|
||||
|
||||
We were told:
|
||||
|
||||
"The Windows Subsystem for Linux (WSL) is a new Windows 10 feature that
|
||||
enables you to run native Linux command-line tools directly on Windows"
|
||||
|
||||
https://docs.microsoft.com/cs-cz/windows/wsl/about
|
||||
|
||||
=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
|
||||
=-= compiling with Cygwin =-=
|
||||
=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
|
||||
|
357
cal/lucas.cal
357
cal/lucas.cal
@@ -1,7 +1,7 @@
|
||||
/*
|
||||
* lucas - perform a Lucas primality test on h*2^n-1
|
||||
*
|
||||
* Copyright (C) 1999,2017 Landon Curt Noll
|
||||
* Copyright (C) 1999,2017,2018 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
|
||||
@@ -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" */
|
||||
@@ -363,41 +373,53 @@ lucas(h, n)
|
||||
return 1; /* 239 is prime */
|
||||
}
|
||||
|
||||
/*
|
||||
* Verify that h*2^n-1 is not a multiple of 3
|
||||
*
|
||||
* The case for h*2^n-1 == 3 is handled above.
|
||||
*/
|
||||
if (((h % 3 == 1) && (n % 2 == 0)) || ((h % 3 == 2) && (n % 2 == 1))) {
|
||||
/* no need to test h*2^n-1, it is a multiple of 3 */
|
||||
ldebug("lucas","not-prime: != 3 and is a multiple of 3");
|
||||
return 0;
|
||||
}
|
||||
|
||||
/*
|
||||
* Avoid any numbers divisible by small primes
|
||||
*/
|
||||
/*
|
||||
* check for 3 <= prime factors < 29
|
||||
* pfact(28)/2 = 111546435
|
||||
* check for 5 <= prime factors < 31
|
||||
* pfact(30)/6 = 1078282205
|
||||
*/
|
||||
testval = h*2^n - 1;
|
||||
if (gcd(testval, 111546435) > 1) {
|
||||
/* a small 3 <= prime < 29 divides h*2^n-1 */
|
||||
ldebug("lucas","not-prime: 3<=prime<29 divides h*2^n-1");
|
||||
if (gcd(testval, 1078282205) > 1) {
|
||||
/* a small 5 <= prime < 31 divides h*2^n-1 */
|
||||
ldebug("lucas",\
|
||||
"not-prime: a small 5<=prime<31 divides h*2^n-1");
|
||||
return 0;
|
||||
}
|
||||
/*
|
||||
* check for 29 <= prime factors < 47
|
||||
* pfact(46)/pfact(28) = 5864229
|
||||
* check for 31 <= prime factors < 53
|
||||
* pfact(52)/pfact(30) = 95041567
|
||||
*/
|
||||
if (gcd(testval, 58642669) > 1) {
|
||||
/* a small 29 <= prime < 47 divides h*2^n-1 */
|
||||
ldebug("lucas","not-prime: 29<=prime<47 divides h*2^n-1");
|
||||
if (gcd(testval, 95041567) > 1) {
|
||||
/* a small 31 <= prime < 53 divides h*2^n-1 */
|
||||
ldebug("lucas","not-prime: 31<=prime<53 divides h*2^n-1");
|
||||
return 0;
|
||||
}
|
||||
/*
|
||||
* check for prime 47 <= factors < 257, if h*2^n-1 is large
|
||||
* 2^282 > pfact(256)/pfact(46) > 2^281
|
||||
* check for prime 53 <= factors < 257, if h*2^n-1 is large
|
||||
* 2^276 > pfact(256)/pfact(52) > 2^275
|
||||
*/
|
||||
bits = highbit(testval);
|
||||
if (bits >= 281) {
|
||||
if (bits >= 275) {
|
||||
if (pprod256 <= 0) {
|
||||
pprod256 = pfact(256)/pfact(46);
|
||||
pprod256 = pfact(256)/pfact(52);
|
||||
}
|
||||
if (gcd(testval, pprod256) > 1) {
|
||||
/* a small 47 <= prime < 257 divides h*2^n-1 */
|
||||
/* a small 53 <= prime < 257 divides h*2^n-1 */
|
||||
ldebug("lucas",\
|
||||
"not-prime: 47<=prime<257 divides h*2^n-1");
|
||||
"not-prime: 53<=prime<257 divides h*2^n-1");
|
||||
return 0;
|
||||
}
|
||||
}
|
||||
@@ -413,7 +435,9 @@ lucas(h, n)
|
||||
* 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.
|
||||
* will never return -1 if h*2^n-1 is not a multiple of 3.
|
||||
* Because the "multiple of 3" case is handled above, the
|
||||
* call below to gen_v1() will never return -1.
|
||||
*/
|
||||
v1 = gen_v1(h, n);
|
||||
if (v1 < 0) {
|
||||
@@ -739,7 +763,7 @@ 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 odd multiple of 3 and large n (some around 1e4, some near 1e6,
|
||||
* 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:
|
||||
*
|
||||
@@ -805,82 +829,155 @@ rodseth_xhn(x, h, n)
|
||||
*
|
||||
* The above distribution was found to hold fairly well over many values of
|
||||
* 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 835823 numbers of the form h*2^n-1
|
||||
* where odd h >= 12996351 is a multiple of 3, n >= 12996351, these are the
|
||||
* smallest v(1) values that were found:
|
||||
*
|
||||
* smallest percentage
|
||||
* v(1) used
|
||||
* -------------------
|
||||
* 3 40.000%
|
||||
* 5 25.683%
|
||||
* 9 11.693%
|
||||
* 11 10.452%
|
||||
* 15 4.806%
|
||||
* 17 2.348%
|
||||
* 21 1.656%
|
||||
* 29 1.281%
|
||||
* 27 0.6881%
|
||||
* 35 0.4536%
|
||||
* 39 0.3121%
|
||||
* 41 0.1760%
|
||||
* 31 0.1414%
|
||||
* 45 0.1173%
|
||||
* 51 0.05576%
|
||||
* 55 0.03300%
|
||||
* 49 0.03185%
|
||||
* 59 0.02090%
|
||||
* 69 0.00980%
|
||||
* 65 0.009367%
|
||||
* 71 0.007205%
|
||||
* 57 0.006341%
|
||||
* 85 0.004611%
|
||||
* 81 0.004179%
|
||||
* 95 0.002882%
|
||||
* 99 0.001873%
|
||||
* 77 0.001153%
|
||||
* 53 0.0007205%
|
||||
* 67 0.0005764%
|
||||
* 125 0.0005764%
|
||||
* 105 0.0005764%
|
||||
* 87 0.0004323%
|
||||
* 111 0.0004323%
|
||||
* 101 0.0002882%
|
||||
* 83 0.0001441%
|
||||
* 127 0.0001196%
|
||||
* 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:
|
||||
*
|
||||
* 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 127287 known primes of the form
|
||||
* h*2^n-1 when h was a multiple of 3, none has an smallest v(1) that was even.
|
||||
* 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 %
|
||||
*
|
||||
* About 1 out of 835000 cases when h is a multiple of 3 use v(1) > 127 as the
|
||||
* smallest value of v(1).
|
||||
* 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 sorted values of X:
|
||||
* these odd values of X:
|
||||
*
|
||||
* 3, 5, 9, 11, 15, 17, 21, 27, 29, 31, 35, 39, 41, 45, 49, 51, 53, 55,
|
||||
* 57, 59, 65, 67, 69, 71, 77, 81, 83, 85, 87, 95, 99, 101, 105, 111, 125
|
||||
* 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 = 12/
|
||||
* 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.
|
||||
* 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 = 35;
|
||||
x_tbl_len = 42;
|
||||
mat x_tbl[x_tbl_len];
|
||||
x_tbl = {
|
||||
3, 5, 9, 11, 15, 17, 21, 27, 29, 31, 35, 39, 41, 45, 49, 51, 53, 55,
|
||||
57, 59, 65, 67, 69, 71, 77, 81, 83, 85, 87, 95, 99, 101, 105, 111, 125
|
||||
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 = 127;
|
||||
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
|
||||
@@ -938,12 +1035,22 @@ next_x = 127;
|
||||
*
|
||||
* 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)
|
||||
@@ -1039,13 +1146,18 @@ next_x = 127;
|
||||
* 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
|
||||
@@ -1066,18 +1178,56 @@ gen_v1(h, n)
|
||||
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;
|
||||
}
|
||||
|
||||
/*
|
||||
* Common Mersenne number case:
|
||||
*
|
||||
* For Mersenne numbers:
|
||||
*
|
||||
* 2^n-1
|
||||
*
|
||||
* we can use, 40% of the time, v(1) == 3. However nearly all code that
|
||||
* implements the Lucas-Lehmer test uses v(1) == 4. Whenever for
|
||||
* h != 0 mod 3, and particular the Mersenne number case of when h == 1:
|
||||
*
|
||||
* 1*2^n-1
|
||||
*
|
||||
* v(1) == 4 always works. For this reason, we return 4 when h == 1.
|
||||
*/
|
||||
if (h == 1) {
|
||||
/* v(1) == 4 always works for the Mersenne number case */
|
||||
return 4;
|
||||
}
|
||||
|
||||
/*
|
||||
* 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:
|
||||
*
|
||||
@@ -1089,26 +1239,51 @@ gen_v1(h, n)
|
||||
* 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 835 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
|
||||
* of odd vules at next_x from here on.
|
||||
*/
|
||||
|
17
calc.c
17
calc.c
@@ -104,6 +104,8 @@ main(int argc, char **argv)
|
||||
int c; /* option */
|
||||
int index;
|
||||
int maxindex;
|
||||
/* fix gcc warning bug */
|
||||
int unusedint = 0;
|
||||
char *cp;
|
||||
char *endcp;
|
||||
char *bp;
|
||||
@@ -278,7 +280,9 @@ main(int argc, char **argv)
|
||||
exit(6);
|
||||
}
|
||||
calc_debug = cp;
|
||||
(void) strtol(cp, &endcp, 10);
|
||||
/* fix gcc warning bug */
|
||||
unusedint =
|
||||
strtol(cp, &endcp, 10);
|
||||
cp = endcp;
|
||||
if (*cp != '\0' &&
|
||||
*cp != ' ' && *cp != ':') {
|
||||
@@ -310,7 +314,9 @@ main(int argc, char **argv)
|
||||
exit(9);
|
||||
}
|
||||
resource_debug = cp;
|
||||
(void) strtol(cp, &endcp, 10);
|
||||
/* fix gcc warning bug */
|
||||
unusedint =
|
||||
strtol(cp, &endcp, 10);
|
||||
cp = endcp;
|
||||
if (*cp != '\0' &&
|
||||
*cp != ' ' && *cp != ':') {
|
||||
@@ -340,7 +346,8 @@ main(int argc, char **argv)
|
||||
exit(12);
|
||||
}
|
||||
user_debug = cp;
|
||||
(void) strtol(cp, &endcp, 10);
|
||||
/* unusedint avoids gcc warning bug */
|
||||
unusedint = strtol(cp, &endcp, 10);
|
||||
cp = endcp;
|
||||
if (*cp != '\0' && *cp != ' ') {
|
||||
fprintf(stderr, "Bad syntax in"
|
||||
@@ -721,8 +728,10 @@ main(int argc, char **argv)
|
||||
printf("main: run_state = %s\n", run_state_name(run_state));
|
||||
|
||||
/*
|
||||
* all done
|
||||
* All done! - Jessica Noll, Age 2
|
||||
*/
|
||||
/* fix gcc warning bug */
|
||||
unusedint++;
|
||||
libcalc_call_me_last();
|
||||
return (run_state == RUN_EXIT_WITH_ERROR ||
|
||||
run_state == RUN_ZERO) ? 1 : 0;
|
||||
|
1
calc.h
1
calc.h
@@ -65,7 +65,6 @@
|
||||
|
||||
#define SYMBOLSIZE 256 /* maximum symbol name size */
|
||||
#define MAXLABELS 100 /* maximum number of user labels in function */
|
||||
#define MAXSTRING 1024 /* maximum size of string constant */
|
||||
#define MAXSTACK 2048 /* maximum depth of evaluation stack */
|
||||
#define MAXFILES 20 /* maximum number of opened files */
|
||||
#define PROMPT1 "> " /* default normal prompt*/
|
||||
|
4
calc.man
4
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
|
||||
|
@@ -348,7 +348,7 @@ EXT=
|
||||
|
||||
# The default calc versions
|
||||
#
|
||||
VERSION= 2.12.6.1
|
||||
VERSION= 2.12.6.7
|
||||
|
||||
# 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
|
||||
|
||||
##################
|
||||
|
@@ -348,7 +348,7 @@ EXT=
|
||||
|
||||
# The default calc versions
|
||||
#
|
||||
VERSION= 2.12.6.1
|
||||
VERSION= 2.12.6.7
|
||||
|
||||
# Names of shared libraries with versions
|
||||
#
|
||||
|
@@ -173,8 +173,6 @@ STATIC struct infoname sys_info[] = {
|
||||
(FULL)MAXSCANCOUNT},
|
||||
{"MAXSTACK", "max depth of evaluation stack", NULL,
|
||||
(FULL)MAXSTACK},
|
||||
{"MAXSTRING", "max size of string constant", NULL,
|
||||
(FULL)MAXSTRING},
|
||||
{"MAXUFULL", "largest FULL value", NULL,
|
||||
(FULL)MAXUFULL},
|
||||
{"MAXULONG", "largest unsigned long val", NULL,
|
||||
|
8
help/mat
8
help/mat
@@ -54,13 +54,13 @@ DESCRIPTION
|
||||
The elements of the matrix are stored internally as a linear array
|
||||
in which locations are arranged in order of increasing indices.
|
||||
For example, in order of location, the six element of A = mat [2,3]
|
||||
are
|
||||
are:
|
||||
|
||||
A[0,0], A[0,1], A[0,2], A[1,0], A[1,,1], A[1,2].
|
||||
A[0,0], A[0,1], A[0,2], A[1,0], A[1,1], and A[1,2].
|
||||
|
||||
These elements may also be specified using the double-bracket operator
|
||||
with a single integer index as in A[[0]], A[[1]], ..., A[[5]].
|
||||
If p is assigned the value &A[0.0], the address of A[[i]] for 0 <= i < 6
|
||||
If p is assigned the value &A[0,0], the address of A[[i]] for 0 <= i < 6
|
||||
is p + i as long as A exists and a new value is not assigned to A.
|
||||
|
||||
When a matrix is created, each element is initially assigned the
|
||||
@@ -414,7 +414,7 @@ SEE ALSO
|
||||
det, inverse, isident, test, config, search, rsearch, reverse, copy,
|
||||
blkcpy, dp, cp, randperm, sort
|
||||
|
||||
## Copyright (C) 1999-2006 Landon Curt Noll
|
||||
## Copyright (C) 1999-2006,2018 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
|
||||
|
@@ -1,5 +1,5 @@
|
||||
NAME
|
||||
randbit - additive 55 shuffle pseudo-random number generator
|
||||
randbit - subtractive 100 shuffle pseudo-random number generator
|
||||
|
||||
SYNOPSIS
|
||||
randbit([x])
|
||||
@@ -42,7 +42,7 @@ LINK LIBRARY
|
||||
SEE ALSO
|
||||
seed, srand, randbit, isrand, random, srandom, israndom
|
||||
|
||||
## Copyright (C) 1999 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 the version 2.1 of the GNU Lesser General Public License
|
||||
|
11
help/seed
11
help/seed
@@ -18,14 +18,7 @@ DESCRIPTION
|
||||
It should be pointed out that the information collected by seed
|
||||
is almost certainly non-chaotic. This function is likely not
|
||||
suitable for applications (such as cryptographic applications)
|
||||
where the unpredictability of seeds is critical. For such critical
|
||||
applications, LavaRnd should be used. See the URL:
|
||||
|
||||
http://www.LavaRnd.org/
|
||||
|
||||
for information about seeding a pseudo-random number generator
|
||||
(such as rand() or random()) with the cryptographic hash of the
|
||||
digitization of chaotic system.
|
||||
where the unpredictability of seeds is critical.
|
||||
|
||||
Given the above warning, this builtin function produces a seed that is
|
||||
suitable for most applications that desire a different pseudo-random
|
||||
@@ -51,7 +44,7 @@ LINK LIBRARY
|
||||
SEE ALSO
|
||||
seed, srand, randbit, isrand, rand, random, srandom, israndom
|
||||
|
||||
## Copyright (C) 1999 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 the version 2.1 of the GNU Lesser General Public License
|
||||
|
@@ -332,13 +332,13 @@ LIMITS
|
||||
iq >= 2^16
|
||||
|
||||
LINK LIBRARY
|
||||
RAND *zsrandom(ZVALUE *pseed, MATRIX *pmat55)
|
||||
RAND *zsetrandom(RAND *state)
|
||||
RANDOM *zsrandom(ZVALUE *pseed, MATRIX *pmat55)
|
||||
RANDOM *zsetrandom(RAND *state)
|
||||
|
||||
SEE ALSO
|
||||
seed, srand, randbit, isrand, random, srandom, israndom
|
||||
|
||||
## Copyright (C) 1999 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 the version 2.1 of the GNU Lesser General Public License
|
||||
|
17
help/types
17
help/types
@@ -1,6 +1,6 @@
|
||||
Builtin types
|
||||
|
||||
The calculator has the following built-in types.
|
||||
The calculator has the following built-in types:
|
||||
|
||||
null value
|
||||
This is the undefined value type. The function 'null'
|
||||
@@ -101,7 +101,20 @@ Builtin types
|
||||
using the result of the 'files' function. Such copies are
|
||||
indistinguishable from each other.
|
||||
|
||||
## Copyright (C) 1999 Landon Curt Noll
|
||||
The calculator also has a number of special types that as used
|
||||
by some special builtin functions:
|
||||
|
||||
rand
|
||||
A subtractive 100 shuffle pseudo-random number generator
|
||||
state. This state is used by functions such as isrand()
|
||||
and srand().
|
||||
|
||||
random
|
||||
A Blum-Blum-Shub pseudo-random number generator state.
|
||||
This state is used by functions such as israndom() and
|
||||
srandom().
|
||||
|
||||
## Copyright (C) 1999,2018 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
|
||||
|
2
hist.h
2
hist.h
@@ -36,7 +36,7 @@
|
||||
#endif
|
||||
|
||||
#ifndef HIST_SIZE
|
||||
#define HIST_SIZE (1024*10)
|
||||
#define HIST_SIZE (1024*32)
|
||||
#endif
|
||||
|
||||
|
||||
|
6
input.c
6
input.c
@@ -59,8 +59,8 @@ E_FUNC FILE *f_open(char *name, char *mode);
|
||||
E_FUNC FILE *curstream(void);
|
||||
|
||||
|
||||
#define TTYSIZE 100 /* reallocation size for terminal buffers */
|
||||
#define MAXDEPTH 10 /* maximum depth of input */
|
||||
#define TTYSIZE 8191 /* reallocation size for terminal buffers */
|
||||
#define MAXDEPTH 255 /* maximum depth of input */
|
||||
#define IS_READ 1 /* reading normally */
|
||||
#define IS_REREAD 2 /* reread current character */
|
||||
#define chartoint(ch) ((ch) & 0xff) /* make sure char is not negative */
|
||||
@@ -789,7 +789,7 @@ ttychar(void)
|
||||
{
|
||||
int ch; /* current char */
|
||||
int len; /* length of current command */
|
||||
STATIC char charbuf[1024];
|
||||
STATIC char charbuf[256*1024];
|
||||
|
||||
/*
|
||||
* If we have more to read from the saved command line, then do that.
|
||||
|
4
qmath.h
4
qmath.h
@@ -254,7 +254,9 @@ E_FUNC NUMBER *swap_HALF_in_NUMBER(NUMBER *dest, NUMBER *src, BOOL all);
|
||||
#define qhighbit(q) (zhighbit((q)->num))
|
||||
#define qlowbit(q) (zlowbit((q)->num))
|
||||
#define qdivcount(q1, q2) (zdivcount((q1)->num, (q2)->num))
|
||||
#define qlink(q) ((q)->links++, (q))
|
||||
/* operation on #q may be undefined, so replace with an inline-function */
|
||||
/* was: #define qlink(q) ((q)->links++, (q)) */
|
||||
static inline NUMBER* qlink(NUMBER* q) { if(q) { (q)->links++; } return q; }
|
||||
|
||||
#define qfree(q) {if (--((q)->links) <= 0) qfreenum(q);}
|
||||
|
||||
|
2
value.c
2
value.c
@@ -2944,7 +2944,7 @@ printestr(VALUE *vp)
|
||||
bp = vp->v_nblock->blk;
|
||||
}
|
||||
i = bp->datalen;
|
||||
math_fmt("%ld,%d)", i, bp->blkchunk);
|
||||
math_fmt("%ld,%d)", i, (int) bp->blkchunk);
|
||||
cp = bp->data;
|
||||
if (i > 0) {
|
||||
math_str("={");
|
||||
|
@@ -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 1 /* minor software level or 0 if not patched */
|
||||
#define MINOR_PATCH 7 /* minor software level or 0 if not patched */
|
||||
|
||||
|
||||
/*
|
||||
|
175
zfunc.c
175
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;
|
||||
|
Reference in New Issue
Block a user