Huge Factorials
Ryan Duchesne
mod p than this
for i from 1 to n do
r:=r*i mod p
od:
then r will be n! mod p the big problem that I am having is that I am
working with numbers 1 to 10^10 which this algorithim is just way too
slow it takes roughly 8 hours to compute (10^9)! mod p I was wondering
if there is a better way to do this I have tried splitting it up into
different loops and calculating it that way, but it still takes just
as long. Any help would be great.
Cheers,
Ryan Duchesne
Steve Bellenot
How big is p? If p <= n then n! = 0 mod p. Is p a prime? A product of
small primes?
A better quesion is what are you doing that requires n! mod p?
--
http://www.math.fsu.edu/~bellenot
bellenot <At/> math.fsu.edu
+1.850.644.7189 (4053fax)
Ryan Duchesne
> In article <7e12f50.03113...@posting.google.com>,
> Ryan Duchesne <catsan...@hotmail.com> wrote:
> >Hello I was wondering if there is a faster way to compute n factorial
> >mod p than this
> >
> > for i from 1 to n do
> > r:=r*i mod p
> > od:
> >
> >then r will be n! mod p the big problem that I am having is that I am
> >working with numbers 1 to 10^10 which this algorithim is just way too
> >slow it takes roughly 8 hours to compute (10^9)! mod p I was wondering
> >if there is a better way to do this I have tried splitting it up into
> >different loops and calculating it that way, but it still takes just
> >as long. Any help would be great.
>
> How big is p? If p <= n then n! = 0 mod p. Is p a prime? A product of
> small primes?
>
> A better quesion is what are you doing that requires n! mod p?
I am trying to find n! mod p for p's which are greater than n. I am
doing large legredre symbols. and I am finding it near impossible to
do this.
Cheers
Alec Mihailovs
news:7e12f50.03120...@posting.google.com...
> I am trying to find n! mod p for p's which are greater than n. I am
> doing large legredre symbols. and I am finding it near impossible to
> do this.
If you meant Legendre symbols, Maple calculates them pretty fast,
> with(numtheory);
> legendre(5*10^19+1,10^20+39);
-1
In general, such calculations as huge factorials mod p should be programmed
directly in assembly. It is pretty easy and doesn't require almost any
assembly knowledge - just few instructions. It can be compiled as a library
(dll) and accessed from Maple through external calling.
Alec Mihailovs
http://webpages.shepherd.edu/amihailo/
Alec Mihailovs
news:dZLyb.457$Zu5.6...@news3.news.adelphia.net...
> In general, such calculations as huge factorials mod p should be
programmed
> directly in assembly. It is pretty easy and doesn't require almost any
> assembly knowledge - just few instructions. It can be compiled as a
library
> (dll) and accessed from Maple through external calling.
Since a few people expressed an interest in that, I am posting here the
assembly code and other instructions. To produce a dll, one needs an
assembler. I used MASM32 v.8 which can be downloaded from
http://www.masm32.com .
To produce a dll, put 3 files, Facmodp.asm, Facmodp.def, and Makeit.bat in
one directory, and run Makeit.bat (by clicking on it in Windows Explorer, or
opening Facmodp.asm in qeditor and running Makeit.bat from the Project
menu). Here are the files:
----------Start of Facmodp.asm----------------
.486
.model flat,stdcall
option casemap:none
include \masm32\include\windows.inc
.code
LibMain proc hInstDLL:HINSTANCE, reason:DWORD, unused:DWORD
mov eax,TRUE
ret
LibMain Endp
Facmodp proc n:DWORD, p:DWORD
push ebx
mov ebx, p
mov ecx, n
mov eax, 1
L:
mul ecx
div ebx
mov eax, edx
dec ecx
jnz L
pop ebx
ret
Facmodp endp
End LibMain
------------End of Facmodp.asm---------------------
------------Start of Facmodp.def-------------------
LIBRARY Facmodp
EXPORTS Facmodp
------------End of Facmodp.def---------------------
------------Start of Makeit.bat--------------------
@echo off
if exist Facmodp.obj del Facmodp.obj
if exist Facmodp.dll del Facmodp.dll
\masm32\bin\ml /c /coff Facmodp.asm
\masm32\bin\Link /SUBSYSTEM:WINDOWS /DLL /DEF:Facmodp.def Facmodp.obj
pause
------------End of Makeit.bat----------------------
That will produce 4 files, Facmodp.dll, Facmodp.lib, Facmodp.exp, and
Facmodp.obj.
In Maple,
> Facmodp:=define_external(
> 'Facmodp',
> 'n'::integer[4],
> 'p'::integer[4],
> 'RETURN'::integer[4],
> 'LIB'="F:/MyProjects/Assembly/Facmodp.dll"
> ):
with changing the LIB to the location of Facmodp.dll creates a function
calculating n! mod p.
For example,
> Facmodp(10^9,10^9+7);
698611116
which took just a few seconds on my computer instead of 8 hours. The
procedure Facmodp does just calculations, without any checking, so it
requires a sane input. Entering p=0 would crash Maple's kernel (so Maple
would have to be restarted after that).
Files Facmodp.def and Makeit.bat are self-explanatory. I'll add here
detailed comments on Facmodp. asm.
.486
.model flat,stdcall
option casemap:none
include \masm32\include\windows.inc
; This is a standart beginning of asm files for MASM32.
; They tell assembler that the program will be used in
; computers with 486 processor or better, with Windows,
; and use stdcall calling convention. casemap:none means
; that the code is case sensitive - that is a requirement
; for using windows.inc include file supplied with MASM32.
.code
; The beginning of code section.
LibMain proc hInstDLL:HINSTANCE, reason:DWORD, unused:DWORD
mov eax,TRUE
ret
LibMain Endp
; Creating the dll entry point. If eax is False,
; the dll won't start, so we put True there.
Facmodp proc n:DWORD, p:DWORD
; The beginning of our function.
; It has two 4byte=32bit parameters, n and p.
push ebx
; In Windows, registers eax, ecx, and edx can be arbitrarily
; modified by programs, but ebx should be preserved, so we are
; pushing it on the stack at the beginning of the program and will
; restore it back before the end of the program.
mov ebx, p
; put p in the ebx register
mov ecx, n
; put n in the (counter) register ecx
mov eax, 1
; put 1 in the eax register
L:
; a label
mul ecx
; multiply eax*ecx and put result in edx:eax, i.e. the beginning of it
; in edx and the last 32 bits in eax
div ebx
; divide edx:eax (the result of multiplication) by ebx=p, and put the
; quotient in eax and the remainder in edx.
mov eax, edx
; move the remainder from edx to eax
dec ecx
; decrease the counter (ecx) by 1. We start from n there. It becomes
; n-1 after first loop, n-2 after second etc.
jnz L
; if counter is not 0, return to label L, otherwise continue below.
pop ebx
; Restore the ebx value from the stack.
ret
; Return.
Facmodp endp
; The end of our function.
End LibMain
;The end of our DLL.
The return value of the function is located in eax register (that is a
standard thing in Windows).
Note that in Maple we used integer[4] data descriptor for n and p. It takes
4 bytes, the same as DWORD. So, we can use values less than 2^31 =
2,147,483,648 for n and p. Actually, keeping in mind that negative integers
have highest bit F, we can use negative integers for representing unsigned
integers up to 2^32-1 = 4,294,967,295.
For the sake of simplicity, I didn't do any optimization of the calculating
procedure, or adding at the beginning that if n>p, then the result is 0, or
checking whether p is 0.
Alec Mihailovs
http://webpages.shepherd.edu/amihailo/
Ryan Duchesne
Cheers
"Alec Mihailovs" <al...@mihailovs.com> wrote in message news:<oHNzb.7942$Zu5.1...@news3.news.adelphia.net>...
Ryan Duchesne
quick question since we use the fact DWORD for n and p which are 4bit
and leave us to a max possible value of 2,000,000,000 is it possible
to increase this number? If it is how would I modify the sample code
in order to achieve this I have tried just a straight multiplication,
but of course this did not work out for me. I promise this is the
last question on this topic. Thank you very much for your help thus
far
Cheers,
Ryan Duchesne
catsan...@hotmail.com (Ryan Duchesne) wrote in message news:<7e12f50.03120...@posting.google.com>...
Carl Devore
> Hello Alec I thank you very much for your time I just have one more
> quick question since we use the fact DWORD for n and p which are 4bit
> and leave us to a max possible value of 2,000,000,000 is it possible
> to increase this number? If it is how would I modify the sample code
> in order to achieve this I have tried just a straight multiplication,
> but of course this did not work out for me.
The following idea deserves to be much more widely known:
It is possible to do exact integer computation with "double" variables
for a much wider range than is possible with 4-byte integer (so-called
"long") variables. I think the range is 2^52 compared to 2^31.
Furthermore, these computations are sometimes faster (depending on
processor) than the 4-byte operations. This is simply because the
double operations have been much more thoroughly engineered because
they are perceived as being more important.
This idea is the key to the amazing speed of Maple's
LinearAlgebra:-Modular package (which works in a 4-byte mode, an
8-byte mode, or unlimited precision).
Ryan Duchesne
Cheers
dev...@math.udel.edu (Carl Devore) wrote in message news:<3a64b911.0312...@posting.google.com>...
Ryan Duchesne
or 8 bit in a Intel surrounding system (32 bit) and the problem that I
am having is that when I use n::integer[4] it gives an error because
the primes that I am using in the provided program (on previous
thread) are larger than 2,000,000,000 I am using numbers higher than
10^10 and so I tried n::float[8] and this gave a float(undefined)
answer where I know that this is impossible I am trying to get this
program to work and I really enjoy Maple but it is really giving me a
problem with the ability to use these large numbers. Any help from
any one would be fantastic thank you all very much for your time.
Cheers,
Ryan
dev...@math.udel.edu (Carl Devore) wrote in message news:<3a64b911.0312...@posting.google.com>...
Stephen Forrest
> I understand that indeed Maple has the amazing ability to computer 4
> or 8 bit in a Intel surrounding system (32 bit) and the problem that I
> am having is that when I use n::integer[4] it gives an error because
> the primes that I am using in the provided program (on previous
> thread) are larger than 2,000,000,000 I am using numbers higher than
> 10^10 and so I tried n::float[8] and this gave a float(undefined)
> answer where I know that this is impossible I am trying to get this
> program to work and I really enjoy Maple but it is really giving me a
> problem with the ability to use these large numbers. Any help from
> any one would be fantastic thank you all very much for your time.
Maple can only deal with 8-byte or 4-byte integers in a specialized
context: when they appear in Matrices, Vectors, or Arrays.
Thus I can do
> a := Array([1,2,3,4], datatype=integer[4]);
and manipulate the Array, etc.
However, you can't simply create 4-bit integers as assign them to
variables, etc. You mentioned using n::integer[8]. Presumably you did
this inside a procedure declaration like this:
> proc(a::integer[4]) 2*a end proc:
Note however, that this syntax does not _declare_ n as an integer, but
merely performs a type check. It's equivalent to
> type(a, integer[4]);
Reading ?type,integer, you can see that an integer a is of type
integer[n] 'if, and only if, expr is an integer representable as an n-bit
hardware integer'. The key bit is 'representable': an integer a may be
representable as a 4-bit integer, but that doesn't mean it's stored as
one.
The situation is slightly different for floats: you cannot store hardware
floats outside of table structures, but you can use fast hardware float
arithmetic using the evalhf command, if the precision offered by your CPU
is sufficient for your needs.
Perhaps if you elaborate more on the purpose of your program, someone
here could offer an opinion on whether the special efficiency of
LinearAlgebra:-Modular could be used to your advantage.
Steve
Carl Devore
> I understand that indeed Maple has the amazing ability to computer 4
> or 8 bit in a Intel surrounding system (32 bit) and the problem that I
> am having is that when I use n::integer[4] it gives an error because
> the primes that I am using in the provided program (on previous
> thread) are larger than 2,000,000,000 I am using numbers higher than
> 10^10 and so I tried n::float[8] and this gave a float(undefined)
> answer where I know that this is impossible I am trying to get this
> program to work and I really enjoy Maple but it is really giving me a
> problem with the ability to use these large numbers. Any help from
> any one would be fantastic thank you all very much for your time.
Sorry, I did not mean to give you the impression that what you want
can be done in Maple and hardware arithmetic at the same time. I was
just pointing out that 8-byte exact arithmetic is possible in
assembly, and that the idea was used by Maple. The modulus for 8-byte
exact arithmetic in LinearAlgebra:-Modular is limited to 2^25 because
it needs to represent exactly the product of two such numbers, and it
needs to be able to add a block of a certain number of these products
together as a dot product.
You want to use moduli greater than 2^31, so your intermediate
products will be up to 64 bits Maybe this is possible using the
guard digits on the Pentium (8-byte arithmetic is done internally in
10 bytes). I don't know enough assembler to know.
Maybe it is time for you to start using an arbitrary-precision
arithmetic package like GMP?
Ryan Duchesne
lagendre symbol for n from 10^9..10^10 for the next fourty primes
which are great for each n. Therefore I have created a program which
finds the next fourty primes for each n and then it computes n! mod
p[i] where i runs from 1 to 40 and then I compute the lagendre of each
of these calculations. With the great skills of Alec who helped me
with a program to compute n! mod p rather fast and although the
program I have created it still takes a long time to run the program,
but that is fine unfortunately when n > 2,000,000,000 the program
quits because it cannot compute n! mod p because of this 4-bit
problem. The mod p program is on an attached string in this group and
I feel that it is as good as this is going to get. I am wondering if
someone knows of a faster way to compute the legendre symbol for n!
mod p? or could someone please tell me what to do to the existing
program that I am using for n mod p and please look at the part that I
am supposed to type into Maple and tell me how I can change this in
order to computer n! greater than 2,000,000,000? I thank you very
much for all of your time.
Cheers,
Ryan Duchesne
Stephen Forrest <safo...@alumni.uwaterloo.ca> wrote in message news:<Pine.LNX.4.43.031210...@perpugilliam.csclub.uwaterloo.ca>...
Edwin Clark
> OK here is the whole idea that I want to do. I would like to find the
> lagendre symbol for n from 10^9..10^10 for the next fourty primes
> which are great for each n.
Have you tried numtheory[legendre]? For example it does the following
computations in the blink of an eye.
> n:=rand(10^99..10^100-1)();
> p:=nextprime(n);
n := 85380320622220857229741217686043056139217455800374092598119\
52655310075487163797179490457039169594160
p := 85380320622220857229741217686043056139217455800374092598119\
52655310075487163797179490457039169594177
> legendre(n,p);
1
But, of course taking n from 10^9 to 10^10 for 40 primes each is going to
take some time no matter how quick your computations are.
Edwin Clark
I see that you want legendre(n!,p) not legendre(n,p).
Someone on the number theory mailing list:
http://listserv.nodak.edu/archives/nmbrthry.html
may know how to compute legendre(n!,p) rapidly.
----Edwin Clark
"Edwin Clark" <ecl...@math.usf.edu> wrote in message
news:PZQBb.2255$Dt6....@twister.tampabay.rr.com...
Jack Fearnley
> OK here is the whole idea that I want to do. I would like to find the
> lagendre symbol for n from 10^9..10^10 for the next fourty primes
> which are great for each n. Therefore I have created a program which
> finds the next fourty primes for each n and then it computes n! mod
> p[i] where i runs from 1 to 40 and then I compute the lagendre of each
> of these calculations. With the great skills of Alec who helped me
> with a program to compute n! mod p rather fast and although the
> program I have created it still takes a long time to run the program,
> but that is fine unfortunately when n > 2,000,000,000 the program
> quits because it cannot compute n! mod p because of this 4-bit
> problem. The mod p program is on an attached string in this group and
> I feel that it is as good as this is going to get. I am wondering if
> someone knows of a faster way to compute the legendre symbol for n!
> mod p? or could someone please tell me what to do to the existing
> program that I am using for n mod p and please look at the part that I
> am supposed to type into Maple and tell me how I can change this in
> order to computer n! greater than 2,000,000,000? I thank you very
> much for all of your time.
>
> Cheers,
>
>
> Ryan Duchesne
>
Here is a way to calculate (n!/q) without overflow. I don't know if it will
be fast enough for you.
Note that the Legendre symbol is fully multiplicative, that is
(mn/q)=(m/q)(n/q).
We use the Legendre formula for factorial:
n! = prod(p^s(n,p),p=2..n prime)
where
s(n,p)=sum([n/p^k],k=1..infinity)
The brackets[] mean "integer part of" and the sum is not really infinite
since [n/p^k] is zero for p^k greater than n.
So (n!/q) = prod((p/q)^s(n,p),p=2..n and prime)
Note that since the Legendre symbol is +1 or -1 we only need to compute the
Legendre symbol when s(n,p) is odd and then only to the first power. This
gives
(n!/q) = prod((p/q),p=2..n and prime and s(n,p) odd)
I hope my notation is clear and that this formulation helps.
Best Regards,
Jack Fearnley
Alec Mihailovs
> OK here is the whole idea that I want to do. I would like to find the
> lagendre symbol for n from 10^9..10^10 for the next fourty primes
> which are great for each n. Therefore I have created a program which
> finds the next fourty primes for each n and then it computes n! mod
> p[i] where i runs from 1 to 40 and then I compute the lagendre of each
> of these calculations. With the great skills of Alec who helped me
> with a program to compute n! mod p rather fast and although the
> program I have created it still takes a long time to run the program,
> but that is fine unfortunately when n > 2,000,000,000 the program
> quits because it cannot compute n! mod p because of this 4-bit
> problem. The mod p program is on an attached string in this group and
> I feel that it is as good as this is going to get. I am wondering if
> someone knows of a faster way to compute the legendre symbol for n!
> mod p? or could someone please tell me what to do to the existing
> program that I am using for n mod p and please look at the part that I
> am supposed to type into Maple and tell me how I can change this in
> order to computer n! greater than 2,000,000,000? I thank you very
> much for all of your time.
First, as I mentioned earlier, negative integers can be used for
representing n and p from 2^31 to 2^32-1, by subtracting 2^32 from them. For
example, to find (2^31+1)! mod (2^31+11), one can do
Facmodp(2^31+1-2^32,2^31+11-2^32);
1455688321
However, it can be done much faster using Wilson's theorem (p-1)! = -1 mod p
for prime p, which gives n! = -1/((n+1)*(n+2)*...*(p-1)) mod p =
(-1)^(p-n)*(p-n-1)! mod p. Both these ideas can be written in Maple as
follows:
f := proc(n, p)
local P;
if n = p - 1 then return p - 1 end if;
P := `if`(2147483648 <= p, p - 4294967296, p);
(-1)^(p - n)/Facmodp(p - n - 1, P) mod p
end proc:
For example,
f(2^31+1,2^31+11);
1455688321
The same result, but obtained much faster. We can calculate 40 Legendre
symbols of n! mod p using following procedure,
F := proc(n)
local A, i, p;
global f;
A := Array(1 .. 40, datatype = integer[1]);
p := n;
for i to 40 do
p := nextprime(p);
A[i] := numtheory[legendre](f(n, p), p)
end do;
eval(A)
end proc:
For example,
A := F(2^31+1):
printf("%d",A);
-1 1 1 -1 1 -1 -1 -1 1 -1 1 1 -1 1 1 1 -1 1 1 1 -1 -1 1 1 1 1 1 -1 -1
1 -1 1 1 1 -1 1 -1 -1 1 -1
It works for n < 2^32 = 4294967296. For larger n, one can use similar
procedures calculating the products of Legendre symbols instead of
factorials. Something like this,
f1:=(n,p)->mul(numtheory[legendre](-i,p),i=1..p-n-1)*numtheory[legendre](-1,
p):
F1 := proc(n)
local A, i, p;
global f1;
A := Array(1 .. 40, datatype = integer[1]);
p := n;
for i to 40 do p := nextprime(p); A[i] := f1(n, p) end do
;
eval(A)
end proc:
Certainly, it is slower than F, but still gives the answer rather fast. For
example,
A1:=F1(10^10):
> printf("%d",A1);
1 1 -1 1 1 1 1 -1 1 -1 -1 1 1 -1 1 -1 -1 1 -1 -1 1 -1 -1 1 -1 1 -1 1 1
1 -1 -1 1 1 1 1 1 -1 -1 -1
Alec Mihailovs
http://webpages.shepherd.edu/amihailo/
Alec Mihailovs
>
> First, as I mentioned earlier, negative integers can be used for
> representing n and p from 2^31 to 2^32-1, by subtracting 2^32 from them.
For
> example, to find (2^31+1)! mod (2^31+11), one can do
>
> Facmodp(2^31+1-2^32,2^31+11-2^32);
>
> 1455688321
I just realized that return also can be negative for p > 2^31. For example,
Facmodp(17,-5);
-288108170
In such cases, the correct answer is greater by 2^32,
%+2^32;
4006859126
17! mod (2^32-5);
4006859126
That means that f should me modified from
> f := proc(n, p)
> local P;
> if n = p - 1 then return p - 1 end if;
> P := `if`(2147483648 <= p, p - 4294967296, p);
> (-1)^(p - n)/Facmodp(p - n - 1, P) mod p
> end proc:
to
f := proc(n, p)
local r;
if n = p - 1 then return p - 1 end if;
if p < 2147483648 then r := Facmodp(p - n - 1, p)
else
r := Facmodp(p - n - 1, p - 4294967296);
r := `if`(r < 0, 4294967296 + r, r)
end if;
(-1)^(p - n)/r mod p
end proc:
> It works for n < 2^32 = 4294967296.
Actually, f works for n < p < 2^32, but since in F 40 primes greater than n
should also be < 2^32, F works for n < 2^32 - 1079 = 4294966217.
Alec Mihailovs
http://webpages.shepherd.edu/amihailo/
Ryan Duchesne
granted great knowledge and I thank you for thank. All the best in
the new year.
Cheers,
Ryan Duchesne
"Alec Mihailovs" <al...@mihailovs.com> wrote in message news:<zyECb.914$9s.2...@news3.news.adelphia.net>...