prime?(-5) false
Martin Rubey
prime?(-5) gives false.
Martin
(Apart from this, I was a little surprised that we don't even have
prime? for polynomials, etc.)
Ralf Hemmecke
> Do we really want to declare all integers less than 2 non-prime? Eg.,
> prime?(-5) gives false.
I'd say, the following line in the "add" part of the category
IntegerNumberSystem is a bug.
prime? x == prime?(x)$IntegerPrimesPackage(%)
Why?
Just compare the documentation...
In IntegerPrimesPackage we have:
prime?: I -> Boolean
++ \spad{prime?(n)} returns true if n is prime and false if not.
++ The algorithm used is Rabin's probabilistic primality test
++ (reference: Knuth Volume 2 Semi Numerical Algorithms).
++ If \spad{prime? n} returns false, n is proven composite.
++ If \spad{prime? n} returns true, prime? may be in error
++ however, the probability of error is very low.
++ and is zero below 25*10^9 (due to a result of Pomerance et al),
++ below 10^12 and 10^13 due to results of Pinch,
++ and below 341550071728321 due to a result of Jaeschke.
++ Specifically, this implementation does at least 10 pseudo prime
++ tests and so the probability of error is \spad{< 4^(-10)}.
++ The running time of this method is cubic in the length
++ of the input n, that is \spad{O( (log n)^3 )}, for n<10^20.
++ beyond that, the algorithm is quartic, \spad{O( (log n)^4 )}.
++ Two improvements due to Davenport have been incorporated
++ which catches some trivial strong pseudo-primes, such as
++ [Jaeschke, 1991] 1377161253229053 * 413148375987157, which
++ the original algorithm regards as prime
whereas the exported function of Integer should match what is given by
UniqueFactorizationDomain, i.e.:
prime?: % -> Boolean
++ prime?(x) tests if x can never be written as the product
++ of two non-units of the ring,
++ i.e., x is an irreducible element.
So better file a bug report.
Ralf
Martin Rubey
> On 01/07/2010 04:34 PM, Martin Rubey wrote:
>> Do we really want to declare all integers less than 2 non-prime? Eg.,
>> prime?(-5) gives false.
>
> I'd say, the following line in the "add" part of the category
> IntegerNumberSystem is a bug.
>
> prime? x == prime?(x)$IntegerPrimesPackage(%)
>
> Why?
>
> Just compare the documentation...
are you saying we should use a deterministic algorithm? I cannot quite
believe this - I'd guess it would make FriCAS unusably slow.
Or did I misunderstand?
Martin
Waldek Hebisch
In UniqueFactorizationDomain we have:
prime?: % -> Boolean
++ prime?(x) tests if x can never be written as the product of two
++ non-units of the ring,
++ i.e., x is an irreducible element.
So the problem is that apropriate polynomial domains should have
UniqueFactorizationDomain.
And for consistency integers should use the same definition, that
is apply absolute value before proper primality test.
--
Waldek Hebisch
heb...@math.uni.wroc.pl
Martin Rubey
> Martin Rubey wrote:
>> Do we really want to declare all integers less than 2 non-prime? Eg.,
>> prime?(-5) gives false.
>>
>> Martin
>>
>> (Apart from this, I was a little surprised that we don't even have
>> prime? for polynomials, etc.)
>
> In UniqueFactorizationDomain we have:
>
> prime?: % -> Boolean
> ++ prime?(x) tests if x can never be written as the product of two
> ++ non-units of the ring,
> ++ i.e., x is an irreducible element.
>
>
> So the problem is that apropriate polynomial domains should have
> UniqueFactorizationDomain.
OK, I see there is a default definition. I guess this goes to the
discussion whether PolynomialFactorizationExplicit is a bad idea...
> And for consistency integers should use the same definition, that
> is apply absolute value before proper primality test.
Yes. Should I apply the following patch?:
prime? n ==
n < 0 => n := -n
n < two => false
Martin
Andrey G. Grozin
I've just committed fricas-1.0.8 to Gentoo portage tree (i.e., it is now
an official package in Gentoo Linux). It can be compiled by 5 different
lisps: sbcl, clisp, gcl, ecl, clozurecl (checked with versions:
sbcl-1.0.34, clisp-2.48, gcl-2.6.8 cvs snapshot, ecl-9.12.3,
clozurecl-1.4). ./configure is called with --disable-aldor --with-x
--with-lisp=$LISP . Strangely enough, I had no success with cmucl (2010-01
snapshot, sse2). When I configure with --with-lisp=lisp (the name of cmucl
executable is lisp), make aborts with
....
; Loading
#P"/var/tmp/portage/sci-mathematics/fricas-1.0.8/work/fricas-1.0.8/src/lisp/fricas-lisp.
sse2f".
T
*
NIL
*
#<The FRICAS-LISP package, 76/190 internal, 1004/1574 external>
* [Doing purification: Done.]
[Undoing binding stack... done]
[Saving current lisp image as executable into
"/var/tmp/portage/sci-mathematics/fricas-1.0.8/work/
fricas-1.0.8/build/i686-pc-linux-gnu/bin/lisp":
[Writing core objects... read-only... static... dynamic... done]
Linking executable...
[/usr/bin/./../lib/cmucl/lib/linker.sh: linking
/var/tmp/portage/sci-mathematics/fricas-1.
0.8/work/fricas-1.0.8/build/i686-pc-linux-gnu/bin/lisp...
+ '[' 3 -ne 3 ']'
+ CCOMPILER=gcc
+ shift
+ '[' gcc = cc ']'
+ PATH=/bin:/usr/bin:/usr/local/bin
++ which gcc
+ GCC=/usr/bin/gcc
+ '[' -z /usr/bin/gcc ']'
++ /usr/bin/gcc -print-libgcc-file-name
+ CRTPATH=/usr/lib/gcc/i686-pc-linux-gnu/4.4.2/libgcc.a
++ dirname /usr/lib/gcc/i686-pc-linux-gnu/4.4.2/libgcc.a
+ LIBROOT=/usr/lib/gcc/i686-pc-linux-gnu/4.4.2
+ echo 'LIBROOT is /usr/lib/gcc/i686-pc-linux-gnu/4.4.2'
LIBROOT is /usr/lib/gcc/i686-pc-linux-gnu/4.4.2
++ uname
+ OPSYS=Linux
+ VER=
+ OUTPUT='-o
/var/tmp/portage/sci-mathematics/fricas-1.0.8/work/fricas-1.0.8/build/i686-pc-linux-gnu/bin/lisp'
++ dirname
/var/tmp/portage/sci-mathematics/fricas-1.0.8/work/fricas-1.0.8/build/i686-pc-linux-gnu/bin/lisp
+
OUTDIR=/var/tmp/portage/sci-mathematics/fricas-1.0.8/work/fricas-1.0.8/build/i686-pc-linux-gnu/bin
++ pwd
+
CURDIR=/var/tmp/portage/sci-mathematics/fricas-1.0.8/work/fricas-1.0.8/src/lisp
+ LINKER=/usr/bin/ld
++ dirname /usr/bin/./../lib/cmucl/lib/linker.sh
+ CMUCLLIB=/usr/bin/./../lib/cmucl/lib
+ OBJS='--whole-archive /usr/bin/./../lib/cmucl/lib/lisp.a
--no-whole-archive'
+ FLAGS=-export-dynamic
+ SCRIPT='-T /usr/bin/./../lib/cmucl/lib/Linux-cmucl-linker-script'
+ BIFLAG='--defsym builtin_image_flag=0x08048000'
+ IFADDR='--defsym initial_function_addr=0x58100111'
+ case "$OPSYS" in
++ uname -m
+ ARCH=i686
+ '[' i686 = x86_64 ']'
+ STARTCRT='/usr/lib/crt1.o /usr/lib/crti.o
/usr/lib/gcc/i686-pc-linux-gnu/4.4.2/crtbegin.o'
+ ENDCRT='/usr/lib/gcc/i686-pc-linux-gnu/4.4.2/crtend.o /usr/lib/crtn.o'
+ DLINKER='-dynamic-linker /lib/ld-linux.so.2'
+ LIBS='-L/usr/lib/gcc/i686-pc-linux-gnu/4.4.2 -ldl -lm -lgcc -lc -lgcc'
+ cd
/var/tmp/portage/sci-mathematics/fricas-1.0.8/work/fricas-1.0.8/build/i686-pc-linux-gnu/bin
+ /usr/bin/ld -T /usr/bin/./../lib/cmucl/lib/Linux-cmucl-linker-script
-dynamic-linker /lib/ld-linux.so.2 -o
/var/tmp/portage/sci-mathematics/fricas-1.0.8/work/fricas-1.0.8/build/i686-pc-linux-gnu/bin/lisp
/usr/lib/crt1.o /usr/lib/crti.o
/usr/lib/gcc/i686-pc-linux-gnu/4.4.2/crtbegin.o -export-dynamic --defsym
builtin_image_flag=0x08048000 --defsym initial_function_addr=0x58100111
--whole-archive /usr/bin/./../lib/cmucl/lib/lisp.a --no-whole-archive
-L/usr/lib/gcc/i686-pc-linux-gnu/4.4.2 -ldl -lm -lgcc -lc -lgcc
/usr/lib/gcc/i686-pc-linux-gnu/4.4.2/crtend.o /usr/lib/crtn.o
/usr/bin/ld:
/var/tmp/portage/sci-mathematics/fricas-1.0.8/work/fricas-1.0.8/build/i686-pc-linux-gnu/bin/lisp:
warning: allocated section `CORRO' not in segment
/usr/bin/ld:
/var/tmp/portage/sci-mathematics/fricas-1.0.8/work/fricas-1.0.8/build/i686-pc-linux-gnu/bin/lisp:
warning: allocated section `CORSTA' not in segment
/usr/bin/ld:
/var/tmp/portage/sci-mathematics/fricas-1.0.8/work/fricas-1.0.8/build/i686-pc-linux-gnu/bin/lisp:
warning: allocated section `CORDYN' not in segment
+ cd
/var/tmp/portage/sci-mathematics/fricas-1.0.8/work/fricas-1.0.8/src/lisp
+ exit 0
done]
done.
echo timestamp > do_it.cmucl
make[2]: Leaving directory
`/var/tmp/portage/sci-mathematics/fricas-1.0.8/work/fricas-1.0.8/src/lisp'
make[2]: Entering directory
`/var/tmp/portage/sci-mathematics/fricas-1.0.8/work/fricas-1.0.8/src/boot'
Building stage 0
mkdir -p -- stage0
rm -rf prev-stage
rm -f stage0/ptyout.sse2f stage0/btincl2.sse2f stage0/btscan2.sse2f
stage0/typrops.sse2f stage0/btpile2.sse2f stage0/typars.sse2f
stage0/tytree1.sse2f
rm -f stage0/ptyout.clisp stage0/btincl2.clisp stage0/btscan2.clisp
stage0/typrops.clisp stage0/btpile2.clisp stage0/typars.clisp
stage0/tytree1.clisp
make OBJECTS="stage0/ptyout.sse2f stage0/btincl2.sse2f
stage0/btscan2.sse2f stage0/typrops.sse2f stage0/btpile2.sse2f
stage0/typars.sse2f stage0/tytree1.sse2f" stage0/bootsys
make[3]: Entering directory
`/var/tmp/portage/sci-mathematics/fricas-1.0.8/work/fricas-1.0.8/src/boot'
/var/tmp/portage/sci-mathematics/fricas-1.0.8/work/fricas-1.0.8/build/scripts/document
--tag=lisp --mode=compile --output=initial-env.sse2f
--use=/var/tmp/portage/sci-mathematics/fricas-1.0.8/work/fricas-1.0.8/build/i686-pc-linux-gnu/bin/lisp
initial-env.lisp
Couldn't map all core sections! Exiting!
make[3]: *** [initial-env.sse2f] Error 255
make[3]: Leaving directory
`/var/tmp/portage/sci-mathematics/fricas-1.0.8/work/fricas-1.0.8/src/boot'
make[2]: *** [stage0/stamp_bootsys] Error 2
make[2]: Leaving directory
`/var/tmp/portage/sci-mathematics/fricas-1.0.8/work/fricas-1.0.8/src/boot'
make[1]: *** [all-boot] Error 2
make[1]: Leaving directory
`/var/tmp/portage/sci-mathematics/fricas-1.0.8/work/fricas-1.0.8/src'
make: *** [all-src] Error 2
Am I doing something stupid? Also, what's the difference between
--with-lisp and --with-lisp-flavor ? Which of them should I use?
Andrey
Waldek Hebisch
> Waldek Hebisch <heb...@math.uni.wroc.pl> writes:
>
> > Martin Rubey wrote:
> >> Do we really want to declare all integers less than 2 non-prime? Eg.,
> >> prime?(-5) gives false.
> >>
> >> Martin
> >>
> >> (Apart from this, I was a little surprised that we don't even have
> >> prime? for polynomials, etc.)
> >
> > In UniqueFactorizationDomain we have:
> >
> > prime?: % -> Boolean
> > ++ prime?(x) tests if x can never be written as the product of two
> > ++ non-units of the ring,
> > ++ i.e., x is an irreducible element.
> >
> >
> > So the problem is that apropriate polynomial domains should have
> > UniqueFactorizationDomain.
>
> OK, I see there is a default definition. I guess this goes to the
> discussion whether PolynomialFactorizationExplicit is a bad idea...
>
IMHO code to support PolynomialFactorizationExplicit is bad, and
we should do it in a different way. We should definitely have
some way to propagate information about factorization so that
domains for which we can factorize have UniqueFactorizationDomain.
ATM it is not clear for me if PolynomialFactorizationExplicit is
good for this purpose.
> > And for consistency integers should use the same definition, that
> > is apply absolute value before proper primality test.
>
> Yes. Should I apply the following patch?:
>
> prime? n ==
> n < 0 => n := -n
> n < two => false
>
Yes.
--
Waldek Hebisch
heb...@math.uni.wroc.pl
Bertfried Fauser
[caution possibly stupid remark]
>> > prime?: % -> Boolean
>> > ++ prime?(x) tests if x can never be written as the product of two
>> > ++ non-units of the ring,
>> > ++ i.e., x is an irreducible element.
* Given this definition, why should -5=(-1)*5 be a prime at all?
* given prime? uses a statistical method, the definition is also not
met and should possibly be stated as:
++ prime?(x) tests if x can not be written as the product of two
++ non-units of the ring (with very small probability)
Cheers and all the best wishes for 2010
BF.
--
% PD Dr Bertfried Fauser
% Research Fellow, School of Computer Science, Univ. of Birmingham
% Honorary Associate, University of Tasmania
% Privat Docent: University of Konstanz, Physics Dept
<http://www.uni-konstanz.de>
% contact |-> URL : http://www.cs.bham.ac.uk/~fauserb/
% Phone : +44-121-41-42795 and +49 1520 9874517
Waldek Hebisch
>
> >> > ? ? ? prime?: % -> Boolean
> >> > ? ? ? ? ? ? ++ prime?(x) tests if x can never be written as the product of two
> >> > ? ? ? ? ? ? ++ non-units of the ring,
> >> > ? ? ? ? ? ? ++ i.e., x is an irreducible element.
>
> * Given this definition, why should -5=(-1)*5 be a prime at all?
-1 is a unit (= invertible element)...
> * given prime? uses a statistical method, the definition is also not
> met and should possibly be stated as:
> ++ prime?(x) tests if x can not be written as the product of two
> ++ non-units of the ring (with very small probability)
>
This is similar to '+' and '*' in Float: both give approximate
result but we still claim that Float is a Field, and do not put
disclaimer in AbelianGroup (and SemiGroup)...
--
Waldek Hebisch
heb...@math.uni.wroc.pl
Ralf Hemmecke
>> met and should possibly be stated as:
>> ++ prime?(x) tests if x can not be written as the product of two
>> ++ non-units of the ring (with very small probability)
> This is similar to '+' and '*' in Float: both give approximate
> result but we still claim that Float is a Field, and do not put
> disclaimer in AbelianGroup (and SemiGroup)...
I wouldn't say this is similar. For prime?$Integer there is a
deterministic correct algorithm. It's only that for efficiency reasons
FriCAS (and most other systems) rather implement a probabilistic test.
Addition in Float is not associative, so Float cannot be a Semigroup.
In the future I'll push stronger for removing the Field property from
Float. (Think about computing a Gröbnerbasis with polynomials over
Float. On identical input, one can end up with completely different
results just by using different critical pair selection strategies in
Buchbergers algorithm.)
Ralf
Ralf Hemmecke
> prime? n ==
> n< 0 => n := -n
> n< two => false
Why not just write
prime? x == prime?(abs x)$IntegerPrimesPackage(%)
?
Ralf
Waldek Hebisch
> Hello *,
>
> I've just committed fricas-1.0.8 to Gentoo portage tree (i.e., it is now
> an official package in Gentoo Linux). It can be compiled by 5 different
> lisps: sbcl, clisp, gcl, ecl, clozurecl (checked with versions:
> sbcl-1.0.34, clisp-2.48, gcl-2.6.8 cvs snapshot, ecl-9.12.3,
> clozurecl-1.4). ./configure is called with --disable-aldor --with-x
> --with-lisp=$LISP .
This is good news, thanks.
> Strangely enough, I had no success with cmucl (2010-01
> snapshot, sse2). When I configure with --with-lisp=lisp (the name of cmucl
> executable is lisp), make aborts with
>
<snip>
> ....
> make[3]: Entering directory
> `/var/tmp/portage/sci-mathematics/fricas-1.0.8/work/fricas-1.0.8/src/boot'
> /var/tmp/portage/sci-mathematics/fricas-1.0.8/work/fricas-1.0.8/build/scripts/document
> --tag=lisp --mode=compile --output=initial-env.sse2f
> --use=/var/tmp/portage/sci-mathematics/fricas-1.0.8/work/fricas-1.0.8/build/i686-pc-linux-gnu/bin/lisp
> initial-env.lisp
> Couldn't map all core sections! Exiting!
> make[3]: *** [initial-env.sse2f] Error 255
>
> Am I doing something stupid?
ATM cmucl port works on some machines and fails on others. From INTALL:
ATM CMU CL port should be considered experimental, it received only
little testing. Also CMU CL seem to have problems on some machines.
The 'seem to have problems on some machines' means that build may
crash like above. Port was done by Anatoly Raportirenko and
for him it worked on a few different machines (both 32-bit and
64-bit). For me it failed on 64-bit machine, but worked in 32-bit
chroot. He also wrote that unicode cmucl have bug, so one have to
patch and rebuild it (I tried non-unicode version).
> Also, what's the difference between
> --with-lisp and --with-lisp-flavor ? Which of them should I use?
>
You should use '--with-lisp', it specifies command to used to run
Lisp compiler. Normally FriCAS should automatically detect which
of supported Lisp is in use. In the (very unlikely) case that
FriCAS can not correctly detect which implementation is in use
one can use '--with-lisp-flavor', but given that automatic detection
it will probably go away in the future.
--
Waldek Hebisch
heb...@math.uni.wroc.pl
Waldek Hebisch
Your way would make 'prime?$Integer' and 'prime?$IntegerPrimesPackage'
incompatible, which may be confusing.
--
Waldek Hebisch
heb...@math.uni.wroc.pl
Andrey G. Grozin
> The 'seem to have problems on some machines' means that build may
> crash like above. Port was done by Anatoly Raportirenko and
> for him it worked on a few different machines (both 32-bit and
> 64-bit). For me it failed on 64-bit machine, but worked in 32-bit
> chroot. He also wrote that unicode cmucl have bug, so one have to
> patch and rebuild it (I tried non-unicode version).
bug in cmucl-20a was fixed soon after the release, and later snapshots
should be OK.)
Is there a testsuit for fricas?
make check
says that there is no target check. I'd like to run the testsuit with all
lisps. First, to see if any of them has some problems. Second, this would
be an interesting benchmark.
Andrey
Martin Rubey
-- probability n prime is > 3/4 for each iteration
-- for most n this probability is much greater than 3/4
t := powmod(p, q, n)
-- neither of these cases tells us anything
if not (one? t or t = nm1) then
for j in 1..k-1 repeat
oldt := t
t := mulmod(t, t, n)
one? t => return true
-- we have squared someting not -1 and got 1
t = nm1 =>
leave
not (t = nm1) => return true
false
shouldn't the penultimate line read
not (oldt = nm1) => return true
??? We know that n is composite if
p^(2^(j+1)*q) is different from -1 and +1 but
p^(2^k*q) = -1
Martin
PS: p is confusing: it's not necessarily a prime, I vote for changing it
to b. Also, the loop would be better to understand if it read
for j in 2..k repeat
Ralf Hemmecke
> Ralf Hemmecke wrote:
>>> Should I apply the following patch?:
>>> prime? n ==
>>> n< 0 => n := -n
>>> n< two => false
>>
>> Why not just write
>>
>> prime? x == prime?(abs x)$IntegerPrimesPackage(%)
> Your way would make 'prime?$Integer' and 'prime?$IntegerPrimesPackage'
> incompatible, which may be confusing.
Then you have misunderstood my note in
http://groups.google.com/group/fricas-devel/msg/8101eaff595fb5f8?hl=en
These two function must be different. By definition, i.e. documentation.
Yes, there is confusion. One of the confusing things is that the
function in IntegerPrimesPackage is called prime?. It should better be
called
rabinPrime?
or something that makes a bit more clear that it is a probabilistic test
on numbers > 0.
Note that IntegerPrimesPackage is just a package and thus functions need
not have names that one would expect in a mathematical context.
A different name would make it also more obvious that
prime? x == rabinPrime?(x)$IntegerPrimesPackage(%)
is a bug because of different documentation of
prime?$UniqueFactorizationDomain and rabinPrime?$IntegerPrimesPackage.
I even disregard here that rabinPrime? is probabilistic. Just the lines
++ \spad{prime?(n)} returns true if n is prime and false if not.
++ (reference: Knuth Volume 2 Semi Numerical Algorithms).
. But how is "n is prime" in I: IntegerNumberSystem actually defined?
I'm not sure that every mathematician in the world would immediately say
that -5 is a prime integer.
If Martin's patch goes in, also the documentation in
IntegerPrimesPackage should be changed appropriately. Currently it says
nothing about negative numbers, but an appropriate note would clarify
the situation.
Side remark: Note that according to
http://mathworld.wolfram.com/PrimeElement.html
the documentation for prime?$UniqueFactorizationDomain is actually
testing for irreducibility.
That
x irreducible implies x prime
in unique factorization domains is a theorem.
Ralf
Waldek Hebisch
make all-input
runs test. Note that currently a few tests prints summary at the
end, for other one have to use textual comparison between runs:
fricas-1.0.8/src/scripts/norm-out build-tree1 build-tree2
The norm-out script is an utility which gets rid of most inessential
differences. Remaining legal differences are due to floating point
(rounding) noise and Lisp diagnostics.
Concerning using tests as a benchmark: I have noticed that after
building from release tarball 'make all-input' is spending
surprisinly large amount of time in activities other than
running tests, so the result has little in common with actual
use of FriCAS (when you build from SVN checkout time outside
tests is reasonably small). Also, most tests are short
running so for them we measure startup time of Lisp executables.
For benchmarking one should probably take on few long running
tests: mapleok, marcbench, r21bugsbig and easter.
BTW: More than year ago I posted my results of
speed comparison:
http://groups.google.com/group/fricas-devel/browse_thread/thread/fbd1c7c9796a7173?hl=en
--
Waldek Hebisch
heb...@math.uni.wroc.pl
Waldek Hebisch
> Andrey G. Grozin wrote:
> This is good news, thanks.
>
> > Strangely enough, I had no success with cmucl (2010-01
> > snapshot, sse2).
>
>
> ATM CMU CL port should be considered experimental, it received only
> little testing. Also CMU CL seem to have problems on some machines.
>
> The 'seem to have problems on some machines' means that build may
> crash like above. Port was done by Anatoly Raportirenko and
> for him it worked on a few different machines (both 32-bit and
> 64-bit). For me it failed on 64-bit machine, but worked in 32-bit
> chroot. He also wrote that unicode cmucl have bug, so one have to
> patch and rebuild it (I tried non-unicode version).
>
part). Namely given three files below I still get crash:
----------------------<Makefile>--<cut here>-----------------
FRICAS_LISP = /mnt/a5/fricas/usr/bin/cmucl
.PHONY: all
all: do_it.cmucl
echo '(quit)' | ./lisp
do_it.cmucl: fricas-lisp.lisp fricas-package.lisp
echo '(load "fricas-package.lisp")' \
'(load (compile-file "fricas-lisp.lisp"))' \
'(in-package "FRICAS-LISP") (save-core "lisp")' \
| $(FRICAS_LISP)
clean:
rm fricas-lisp.sse2f lisp
---------------------<cut here>-------------------------------
----------------<fricas-package.lisp>----<cut here>-----------
;;; We put this in separate file to avoid problems with compilation.
(eval-when (:execute :compile-toplevel :load-toplevel)
(if (not (find-package "FRICAS-LISP"))
(make-package "FRICAS-LISP"
:use (list (or (find-package "COMMON-LISP")
"LISP")))))
---------------------<cut here>-------------------------------
----------------<fricas-lisp.lisp>----<cut here>--------------
(in-package "FRICAS-LISP")
;; Save current image on disk as executable and quit.
(defun save-core (core-image)
(let ((top-fun #'(lambda ()
(lisp::%top-level))))
(ext::save-lisp
(unix::unix-maybe-prepend-current-directory core-image)
:init-function top-fun :executable t :print-herald nil))
)
---------------------<cut here>-------------------------------
--
Waldek Hebisch
heb...@math.uni.wroc.pl
Martin Rubey
> We know that n is composite if
>
> p^(2^(j+1)*q) is different from -1 and +1 but
>
> p^(2^k*q) = -1
Actually, I can't find an example for this behaviour. Can this happen
at all?
Martin
Martin Rubey
rabinProvesComposite and prime?
There are at least two things unclear. For example, we have (globally)
rootsMinus1:Set I := empty()
-- used to check whether we detect too many roots of -1
count2Order:Vector NonNegativeInteger := new(1,0)
-- used to check whether we observe an element of maximal two-order
rootsMinus1 is the set {b^(q 2^(j-1)) : b^(q 2^j) = -1 mod n}
Clearly, when rootsMinus1 contains more than two elements, then n cannot
be a prime. However, I could not find *any* number such that this
condition would be triggered in rabinProvesComposite.
count2Order(j) = #{b: b^(q 2^(j-1)) = -1 mod n}
When count2Order(k)>0 (remember q 2^k = n-1) then n must be prime, since
then b has order n-1. However, count2Order(j) for smaller j seems to be
unused, and I can't think of a good use for them. Moreover, the first
test count2Order(k)>0 appears very late. Testing earlier seems to be
*very* effective.
Can I commit? (after a little more testing
rabinProvesCompositeSmall(p,n,nm1,q,k) ==
-- probability n prime is > 3/4 for each iteration
-- for most n this probability is much greater than 3/4
t := powmod(p, q, n)
-- neither of these cases tells us anything
if one? t or t = nm1 then
false
else
for j in 1..k-1 repeat
t := mulmod(t, t, n)
-- we have squared someting not -1 and got 1
one? t => return true
-- we encountered a -1 before the very end
t = nm1 => return false
-- all but possibly b^(2^k*q) are different from +1 and -1
true
rabinProvesComposite(p,n,nm1,q,k) ==
-- probability n prime is > 3/4 for each iteration
-- for most n this probability is much greater than 3/4
t := powmod(p, q, n)
if t = nm1 then count2Order(1) := count2Order(1)+1
-- neither of these cases tells us anything
if one? t or t = nm1 then
false
else
for j in 1..k-1 repeat
oldt := t
t := mulmod(t, t, n)
-- we have squared someting not -1 and got 1
one? t => return true
t = nm1 =>
rootsMinus1 := union(rootsMinus1, oldt)
# rootsMinus1 > 2 =>
return true -- Z/nZ can't be a field
count2Order(j+1) := count2Order(j+1)+1
return false
-- all but possibly b^(2^k*q) are different from +1 and -1
true
prime? n ==
n < two => false
n < nextSmallPrime => member?(n, smallPrimes)
-- not one? gcd(n, productSmallPrimes) => false
not (gcd(n, productSmallPrimes) = 1) => false
n < nextSmallPrimeSquared => true
nm1 := n-1
q := (nm1) quo two
for k in 1.. while not odd? q repeat q := q quo two
-- q = (n-1) quo 2^k for largest possible k
n < JaeschkeLimit =>
rabinProvesCompositeSmall(2::I,n,nm1,q,k) => return false
rabinProvesCompositeSmall(3::I,n,nm1,q,k) => return false
n < PomeranceLimit =>
rabinProvesCompositeSmall(5::I,n,nm1,q,k) => return false
member?(n,PomeranceList) => return false
true
rabinProvesCompositeSmall(7::I,n,nm1,q,k) => return false
n < PinchLimit =>
rabinProvesCompositeSmall(10::I,n,nm1,q,k) => return false
member?(n,PinchList) => return false
true
rabinProvesCompositeSmall(5::I,n,nm1,q,k) => return false
rabinProvesCompositeSmall(11::I,n,nm1,q,k) => return false
n < PinchLimit2 =>
member?(n,PinchList2) => return false
true
rabinProvesCompositeSmall(13::I,n,nm1,q,k) => return false
rabinProvesCompositeSmall(17::I,n,nm1,q,k) => return false
true
rootsMinus1:= empty()
count2Order := new(k,0) -- vector of k zeroes
mn := minIndex smallPrimes
for i in mn+1..mn+10 repeat
rabinProvesComposite(smallPrimes i,n,nm1,q,k) => return false
(count2Order(k) > 0) => return true
import IntegerRoots(I)
q > 1 and perfectSquare?(3*n+1) => false
((n9:=n rem (9::I))=1 or n9 = -1) and perfectSquare?(8*n+1) => false
-- Both previous tests from Damgard & Landrock
currPrime:=smallPrimes(mn+10)
probablySafe:=tenPowerTwenty
while count2Order(k) = 0 or n > probablySafe repeat
currPrime := nextPrime currPrime
probablySafe:=probablySafe*(100::I)
rabinProvesComposite(currPrime,n,nm1,q,k) => return false
if n <= probablySafe then
print("count2Order"::OutputForm)
print(count2Order::OutputForm)
true
Waldek Hebisch
> I played a little more with PRIMES, more precisely with
> rabinProvesComposite and prime?
>
> There are at least two things unclear. For example, we have (globally)
>
> rootsMinus1:Set I := empty()
> -- used to check whether we detect too many roots of -1
> count2Order:Vector NonNegativeInteger := new(1,0)
> -- used to check whether we observe an element of maximal two-order
>
> rootsMinus1 is the set {b^(q 2^(j-1)) : b^(q 2^j) = -1 mod n}
>
> Clearly, when rootsMinus1 contains more than two elements, then n cannot
> be a prime. However, I could not find *any* number such that this
> condition would be triggered in rabinProvesComposite.
>
AFAIK such cases are extremally rare, it is almost impossible to find
such case via random search and systemtic methods to find them
are quite heavy.
> count2Order(j) = #{b: b^(q 2^(j-1)) = -1 mod n}
>
> When count2Order(k)>0 (remember q 2^k = n-1) then n must be prime, since
> then b has order n-1. However, count2Order(j) for smaller j seems to be
> unused, and I can't think of a good use for them. Moreover, the first
> test count2Order(k)>0 appears very late. Testing earlier seems to be
> *very* effective.
>
> Can I commit? (after a little more testing
Please wait a bit: I looked at rabinProvesCompositeSmall in the past
and it looked correct and well optimized. It is quite easy to
make it worse -- I would like to have some time to look more
carefully at your version.
--
Waldek Hebisch
heb...@math.uni.wroc.pl
Martin Rubey
> Martin Rubey wrote:
>> I played a little more with PRIMES, more precisely with
>> rabinProvesComposite and prime?
>>
>> There are at least two things unclear. For example, we have (globally)
>>
>> rootsMinus1:Set I := empty()
>> -- used to check whether we detect too many roots of -1
>> count2Order:Vector NonNegativeInteger := new(1,0)
>> -- used to check whether we observe an element of maximal two-order
>>
>> rootsMinus1 is the set {b^(q 2^(j-1)) : b^(q 2^j) = -1 mod n}
>>
>> Clearly, when rootsMinus1 contains more than two elements, then n cannot
>> be a prime. However, I could not find *any* number such that this
>> condition would be triggered in rabinProvesComposite.
>>
>
> AFAIK such cases are extremally rare, it is almost impossible to find
> such case via random search and systemtic methods to find them
> are quite heavy.
right. I just found Davenports report, and
56897193526942024370326972321
is such a number. Moreover, I found that
1195068768795265792518361315725116351898245581
is considered prime after my modification, which is not true. I'll put
these numbers (and any others I can find in Davenport's paper) into the
testfile.
Martin
Martin Rubey
>> AFAIK such cases are extremally rare, it is almost impossible to find
>> such case via random search and systemtic methods to find them
>> are quite heavy.
>
> right. I just found Davenports report, and
>
> 56897193526942024370326972321
>
> is such a number. Moreover, I found that
>
> 1195068768795265792518361315725116351898245581
>
> is considered prime after my modification, which is not true. I'll put
> these numbers (and any others I can find in Davenport's paper) into the
> testfile.
I just asked a colleague to explain the two “maximal 2-part” test from
Primality Testing Revisited by Davenport
http://portal.acm.org/citation.cfm?id=143290
to me. I think I should put it into the documentation. Short
explanation:
count2Order(k) counts the number of b such that
b^((n-1)/2) = -1 (mod n).
If n is prime, then b^((n-1)/2) is just the Legendre symbol
legendre(b, n), and -1 for exactly half of the b's. Of course, if n is
not prime, we might still get some -1, which explains why we should
first do some checks disregarding count2Order.
We could not see why one tests count2Order(k) = 0, and does not require
count2Order(k) to exceed a proportion of the number of tests, but that's
certainly only tuning.
It is also unclear why Davenport collects more information than
needed...
Well, that's it for the moment.
Martin
Martin Rubey
I was just looking at your modifications of the primality test in axiom,
especially th two “maximal 2-part” test from Primality Testing Revisited.
After I while I found the following short explanation:
count2Order(k) counts the number of b such that
b^((n-1)/2) = -1 (mod n).
If n is prime, then b^((n-1)/2) is just the Legendre symbol
legendre(b, n), and -1 for exactly half of the b's. Of course, if n is
not prime, we might still get some -1, which explains why we should
first do some checks disregarding count2Order.
I could not see however why one tests count2Order(k) = 0, and does not
require count2Order(k) to exceed a proportion of the number of tests.
It is also unclear to my why you collect more information than, namely
the number of b^(q*2^j) congruent to -1 mod n for all j.
Did you have plans for a more accurate check? For your convenience, I
attach the source.
Besides that, I also noticed the failure prime?(149491*747451*34233211)
giving true. Zhenxiang Zhang, Two kinds of strong pseudoprimes up to
10^36 conjectures that this number is actually psi_11, i.e. the smallest
strong pseudoprime to all the first 11 prime bases. Maybe we need to
update the algorithm, it's just one more step :-)
(49) -> rabinProvesComposite(37, n, n-1, q, k)
(49) true
I hope you are doing well,
Martin Rubey
rabinProvesCompositeSmall(p,n,nm1,q,k) ==
-- probability n prime is > 3/4 for each iteration
-- for most n this probability is much greater than 3/4
t := powmod(p, q, n)
-- neither of these cases tells us anything
if not (one? t or t = nm1) then
for j in 1..k-1 repeat
oldt := t
t := mulmod(t, t, n)
one? t => return true
-- we have squared someting not -1 and got 1
t = nm1 =>
leave
not (t = nm1) => return true
false
rabinProvesComposite(p,n,nm1,q,k) ==
-- probability n prime is > 3/4 for each iteration
-- for most n this probability is much greater than 3/4
t := powmod(p, q, n)
-- neither of these cases tells us anything
if t=nm1 then count2Order(1):=count2Order(1)+1
if not (one? t or t = nm1) then
for j in 1..k-1 repeat
oldt := t
t := mulmod(t, t, n)
one? t => return true
-- we have squared someting not -1 and got 1
t = nm1 =>
rootsMinus1:=union(rootsMinus1,oldt)
count2Order(j+1):=count2Order(j+1)+1
leave
not (t = nm1) => return true
# rootsMinus1 > 2 => true -- Z/nZ can't be a field
false
prime? n ==
n < two => false
n < nextSmallPrime => member?(n, smallPrimes)
not one? gcd(n, productSmallPrimes) => false
n < nextSmallPrimeSquared => true
import IntegerRoots(I)
q > 1 and perfectSquare?(3*n+1) => false
((n9:=n rem (9::I))=1 or n9 = -1) and perfectSquare?(8*n+1) => false
-- Both previous tests from Damgard & Landrock
currPrime:=smallPrimes(mn+10)
probablySafe:=tenPowerTwenty
while count2Order(k) = 0 or n > probablySafe repeat
currPrime := nextPrime currPrime
probablySafe:=probablySafe*(100::I)
rabinProvesComposite(currPrime,n,nm1,q,k) => return false
true