Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

More GMP work #1075

Merged
merged 21 commits into from
Jan 23, 2017
Merged

More GMP work #1075

Show file tree
Hide file tree
Changes from 1 commit
Commits
Show all changes
21 commits
Select commit Hold shift + click to select a range
1898ea7
gmpints: introduce fake_mpz_t
fingolfin Jan 21, 2017
8c018e0
gmpints: merge GMPorINTOBJ_INT and ObjInt_Int
fingolfin Jan 21, 2017
e918054
gmpints: rewrite PrintInt using mpz_get_str
fingolfin Jan 21, 2017
1c6d8ee
gmpints: rewrite HexStringInt and STRING_INT using fake_mpz_t
fingolfin Jan 21, 2017
b591686
gmpints: remove GMP_INTOBJ
fingolfin Jan 21, 2017
9a3689f
gmpints: rewrite FuncIntHexString
fingolfin Jan 21, 2017
99418d3
gmpints: replace TypLimb -> mp_limb_t, TypGMPSize -> mp_size_t
fingolfin Jan 21, 2017
b700ea4
gmpints: replace INTEGER_ALLOCATION_SIZE by INTEGER_UNIT_SIZE
fingolfin Jan 21, 2017
af567b1
gmpints: cleanup gmpints.h
fingolfin Jan 21, 2017
6a5815f
gmpints: simplify ModInt
fingolfin Jan 21, 2017
b8a188c
gmpints: remove NEW_INTPOS
fingolfin Jan 21, 2017
8454a6f
gmpints: add IS_INT, IS_POSITIVE, REQUIRE_INT_ARG
fingolfin Jan 21, 2017
1c8ae24
gmpints: change REQUIRE_INT_ARG to use ErrorMayQuit
fingolfin Jan 21, 2017
cdba5d4
gmpints: use REQUIRE_INT_ARG in two more places
fingolfin Jan 21, 2017
3e44492
gmpints: fix RemInt error handling
fingolfin Jan 21, 2017
42c4845
gmpints: add JACOBI_INT and change Jacobi() to use it
fingolfin Jan 21, 2017
0950124
gmpints: add POWERMODINT and INVMODINT
fingolfin Jan 21, 2017
9561e9d
gmpints: add PVALUATION_INT
fingolfin Jan 21, 2017
fbf656a
tests: replace legal header in intarith.tst
fingolfin Jan 21, 2017
9bc7f34
gmpints: explain the buffer in PrintInt
fingolfin Jan 23, 2017
cad52ab
gmpints: improve initial documentation comment
fingolfin Jan 23, 2017
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
Prev Previous commit
Next Next commit
gmpints: add JACOBI_INT and change Jacobi() to use it 
Use GMP's mpz_kronecker to compute the Kronecker-Jacobi symbol J(n,m).
We retain the old GAP implementation for reference, and also to
facilitate testing (we want to make sure that old and new code return
identical results).

Note that the new implementation accepts some more inputs, namely it
allows negative values for m. This is because GMP allows this, and
it seems pointless to restrict this again.

The new implementation is considerably faster. Old code:

gap> for a in [1..1000] do for b in [1..1000] do Jacob(a,b);od;od;time;
1979

New code:

gap> for a in [1..1000] do for b in [1..1000] do Jacob(a,b);od;od;time;
85
  • Loading branch information
@fingolfin
fingolfin committed Jan 21, 2017
commit 42c4845af6a163e41f8986ce0ca9263790190b44
40 changes: 1 addition & 39 deletions lib/numtheor.gi
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -333,45 +333,7 @@ end );
#F Jacobi( <n>, <m> ) . . . . . . . . . . . . . . . . . . . . Jacobi symbol
##
##
InstallGlobalFunction( Jacobi, function ( n, m )
local jac, t;

# check the argument
if m <= 0 then Error("<m> must be positive"); fi;

# compute the Jacobi symbol similar to Euclid's algorithm
jac := 1;
while m <> 1 do

# if the gcd of $n$ and $m$ is $>1$ Jacobi returns $0$
if n = 0 or (n mod 2 = 0 and m mod 2 = 0) then
jac := 0; m := 1;

# $J(n,2*m) = J(n,m) * J(n,2) = J(n,m) * (-1)^{(n^2-1)/8}$
elif m mod 2 = 0 then
if n mod 8 = 3 or n mod 8 = 5 then jac := -jac; fi;
m := m / 2;

# $J(2*n,m) = J(n,m) * J(2,m) = J(n,m) * (-1)^{(m^2-1)/8}$
elif n mod 2 = 0 then
if m mod 8 = 3 or m mod 8 = 5 then jac := -jac; fi;
n := n / 2;

# $J(-n,m) = J(n,m) * J(-1,m) = J(n,m) * (-1)^{(m-1)/2}$
elif n < 0 then
if m mod 4 = 3 then jac := -jac; fi;
n := -n;

# $J(n,m) = J(m,n) * (-1)^{(n-1)*(m-1)/4}$ (quadratic reciprocity)
else
if n mod 4 = 3 and m mod 4 = 3 then jac := -jac; fi;
t := n; n := m mod n; m := t;

fi;
od;

return jac;
end );
InstallGlobalFunction( Jacobi, JACOBI_INT );

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

While its great to be able to do all these in GMP, I would quite like to keep GAP functions around as an alternative. The reason is that while implementing GAP in alternative runtimes, the more code one does not have to write in the runtime the better.

I don't know how to do this cleanly right now though.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

That is exactly how I did it before @ChrisJefferson asked me to remove the GAP implementation sigh. It is still around in the .tst files, by the way.

I can add it back here again, but I don't want to rewrite the code a third time, so I'll wait until at least @ChrisJefferson signs off on it, too.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Caveat: the new Jacobi accepts more inputs than the old (n negative).

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

As discussed in the chat, I am happy with moving the function to tests, and I will attempt to make a structured attempt at detecting which functions are delegated to C code (and collect/provide GAP implementations where appropriate).


#############################################################################
Expand Down
34 changes: 34 additions & 0 deletions src/gmpints.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1962,6 +1962,37 @@ Obj FuncGCD_INT ( Obj self, Obj opL, Obj opR )
}


/****************************************************************************
**
*/
Obj JacobiInt ( Obj opL, Obj opR )
{
fake_mpz_t mpzL, mpzR;
int result;

CHECK_INT(opL);
CHECK_INT(opR);

FAKEMPZ_GMPorINTOBJ( mpzL, opL );
FAKEMPZ_GMPorINTOBJ( mpzR, opR );

result = mpz_kronecker( MPZ_FAKEMPZ(mpzL), MPZ_FAKEMPZ(mpzR) );
CHECK_FAKEMPZ(mpzL);
CHECK_FAKEMPZ(mpzR);

return INTOBJ_INT( result );
}

/****************************************************************************
**
*/
Obj FuncJACOBI_INT ( Obj self, Obj opL, Obj opR )
{
REQUIRE_INT_ARG( "JacobiInt", "left", opL );
REQUIRE_INT_ARG( "JacobiInt", "right", opR );
return JacobiInt( opL, opR );
}

/****************************************************************************
**
** * * * * * * * "Mersenne twister" random numbers * * * * * * * * * * * * *
Expand Down Expand Up @@ -2105,6 +2136,9 @@ static StructGVarFunc GVarFuncs [] = {
{ "POW_OBJ_INT", 2, "obj, gmp",
FuncPOW_OBJ_INT, "src/gmpints.c:POW_OBJ_INT" },

{ "JACOBI_INT", 2, "gmp1, gmp2",
FuncJACOBI_INT, "src/gmpints.c:JACOBI_INT" },

{ "HexStringInt", 1, "gmp",
FuncHexStringInt, "src/gmpints.c:HexStringInt" },

Expand Down
1 change: 1 addition & 0 deletions tst/testinstall/intarith.tst
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1154,6 +1154,7 @@ gap> Print(2^(64*1050), "\n");
948688599817481827330434634041214374947601458440559056883819756862448124537321\
8157000453169824579899621376


#
gap> STOP_TEST( "intarith.tst", 290000);

Expand Down
45 changes: 45 additions & 0 deletions tst/testinstall/numtheor.tst
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -10,5 +10,50 @@ gap> PValuation(100,3);
gap> PValuation(13/85,5);
-1

#
# Compare GAP and C implementations of Jacobi()
#
gap> JACOBI_INT_GAP := function ( n, m )
> local jac, t;
>
> # check the argument
> if m <= 0 then Error("<m> must be positive"); fi;
>
> # compute the Jacobi symbol similar to Euclid's algorithm
> jac := 1;
> while m <> 1 do
>
> # if the gcd of $n$ and $m$ is $>1$ Jacobi returns $0$
> if n = 0 or (n mod 2 = 0 and m mod 2 = 0) then
> jac := 0; m := 1;
>
> # $J(n,2*m) = J(n,m) * J(n,2) = J(n,m) * (-1)^{(n^2-1)/8}$
> elif m mod 2 = 0 then
> if n mod 8 = 3 or n mod 8 = 5 then jac := -jac; fi;
> m := m / 2;
>
> # $J(2*n,m) = J(n,m) * J(2,m) = J(n,m) * (-1)^{(m^2-1)/8}$
> elif n mod 2 = 0 then
> if m mod 8 = 3 or m mod 8 = 5 then jac := -jac; fi;
> n := n / 2;
>
> # $J(-n,m) = J(n,m) * J(-1,m) = J(n,m) * (-1)^{(m-1)/2}$
> elif n < 0 then
> if m mod 4 = 3 then jac := -jac; fi;
> n := -n;
>
> # $J(n,m) = J(m,n) * (-1)^{(n-1)*(m-1)/4}$ (quadratic reciprocity)
> else
> if n mod 4 = 3 and m mod 4 = 3 then jac := -jac; fi;
> t := n; n := m mod n; m := t;
>
> fi;
> od;
>
> return jac;
> end;;
gap> ForAll([-100 .. 100], a-> ForAll([1 .. 100], b -> JACOBI_INT_GAP(a,b)=JACOBI_INT(a,b)));
true

#
gap> STOP_TEST( "numtheor.tst", 290000);