Skip to content

Commit 992b43f

Browse files
mikestillmanmikestillman
committed
Add 3 functions that Frank Schreyer wrote: nextPrime, getPrimeWithRootOfUnity, randomKRationalPoint
1 parent 325f40d commit 992b43f

File tree

8 files changed

+266
-78
lines changed

8 files changed

+266
-78
lines changed

‎M2/Macaulay2/m2/exports.m2‎

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -309,6 +309,7 @@ export {
309309
"Quotient",
310310
"QuotientRing",
311311
"RR",
312+
"Range",
312313
"RealField",
313314
"Reduce",
314315
"Reload",
@@ -665,6 +666,7 @@ export {
665666
"getNetFile",
666667
"getNonUnit",
667668
"getPackage",
669+
"getPrimeWithRootOfUnity",
668670
"getSymbol",
669671
"getWWW",
670672
"getc",
@@ -868,6 +870,7 @@ export {
868870
"newRing",
869871
"newline",
870872
"nextkey",
873+
"nextPrime",
871874
"norm",
872875
"not",
873876
"notImplemented",
@@ -950,6 +953,7 @@ export {
950953
"quotientRemainder'",
951954
"radical",
952955
"random",
956+
"randomKRationalPoint",
953957
"randomMutableMatrix",
954958
"rank",
955959
"read",

‎M2/Macaulay2/m2/quotring.m2‎

Lines changed: 50 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -323,6 +323,56 @@ getNonUnit = R -> if R.?Engine and R.Engine then (
323323
r := rawGetNonUnit raw R;
324324
if r != 0 then new R from r)
325325

326+
-- nextPrime and getPrimeWithRootOfUnity, written by Frank Schreyer.
327+
nextPrime=method(TypicalValue=>ZZ)
328+
nextPrime Number := n -> (
329+
n0 := ceiling n;
330+
if n0 <= 2 then return 2;
331+
if even n0 then n0=n0+1;
332+
while not isPrime n0 do n0=n0+2;
333+
n0
334+
)
335+
336+
337+
getPrimeWithRootOfUnity = method(Options=>{Range=>(10^4,3*10^4)})
338+
getPrimeWithRootOfUnity(ZZ,ZZ) := opt-> (n,r1) -> (
339+
a := opt.Range_0;
340+
b := opt.Range_1;
341+
--(a,b)=(100,200),n=12,r=12
342+
if b<a+2*log a or b<2 or not ( class a === ZZ and class b === ZZ) then
343+
error "the range (a,b) should be an integral interval of diameter b-a > 2*log a";
344+
if n==2 then return (nextPrime(a+random(b-a)),-1);
345+
r2 := r1-1;
346+
primeFactors := apply(toList factor n, f-> first f); -- the prime factors of n
347+
ds := apply(primeFactors, d->lift(n/d,ZZ)); -- the largest factors of n
348+
L := toList factor(r2^n-1);
349+
q := last L;
350+
p := first q;
351+
while (
352+
while p>b or p<a do (
353+
if #L>1 and p>=a then (
354+
L = delete(q,L);
355+
q = last L;
356+
p = first q
357+
)
358+
else (
359+
r2 = r2+1;
360+
L = toList factor(r2^n-1);
361+
q = last L;
362+
p = first q
363+
);
364+
);
365+
-- r2 is now a root of unity in ZZ/p with a <= p <= b
366+
#select(ds, d-> (r2^d)%p == 1) !=0) --check for primitive
367+
do (
368+
r2 = r2+1;
369+
L = toList factor(r2^n-1);
370+
q = last L;
371+
p = first q
372+
);
373+
(p,r2)
374+
);
375+
326376
-- Local Variables:
327377
-- compile-command: "make -C $M2BUILDDIR/Macaulay2/m2 "
328378
-- End:

‎M2/Macaulay2/m2/varieties.m2‎

Lines changed: 39 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -443,6 +443,45 @@ euler ProjectiveVariety := X -> (
443443
d := dim X;
444444
sum(0 .. d, j -> hh^(j,j) X + 2 * sum(0 .. j-1, i -> (-1)^(i+j) * hh^(i,j) X)))
445445

446+
------------------------------------
447+
-- Code donated by Frank Schreyer --
448+
------------------------------------
449+
450+
randomKRationalPoint = method()
451+
randomKRationalPoint Ideal := I -> (
452+
R:=ring I;
453+
if char R == 0 then error "expected a finite ground field";
454+
if not class R === PolynomialRing then error "expected an ideal in a polynomial ring";
455+
if not isHomogeneous I then error "expected a homogenous ideal";
456+
n:=dim I;
457+
if n<=1 then error "expected a positive dimensional scheme";
458+
c:=codim I;
459+
Rs:=R;
460+
Re:=R;
461+
f:=I;
462+
if not c==1 then (
463+
-- projection onto a hypersurface
464+
parametersystem:=ideal apply(n,i->R_(i+c));
465+
if not dim(I+parametersystem)== 0 then return print "make coordinate change";
466+
kk:=coefficientRing R;
467+
Re=kk(monoid[apply(dim R,i->R_i),MonomialOrder => Eliminate (c-1)]);
468+
rs:=(entries selectInSubring(1,vars Re))_0;
469+
Rs=kk(monoid[rs]);
470+
f=ideal substitute(selectInSubring(1, generators gb substitute(I,Re)),Rs);
471+
if not degree I == degree f then return print "make coordinate change"
472+
);
473+
H:=0;pts:=0;pts1:=0;trial:=1;pt:=0;ok:=false;
474+
while (
475+
H=ideal random(Rs^1,Rs^{dim Rs-2:-1});
476+
pts=decompose (f+H);
477+
pts1=select(pts,pt-> degree pt==1 and dim pt ==1);
478+
ok=( #pts1>0);
479+
if ok then (pt=saturate(substitute(pts1_0,R)+I);ok==(degree pt==1 and dim pt==0));
480+
not ok) do (trial=trial+1);
481+
pt
482+
)
483+
484+
446485
-- Local Variables:
447486
-- compile-command: "make -C $M2BUILDDIR/Macaulay2/m2 "
448487
-- End:

‎M2/Macaulay2/packages/Macaulay2Doc/changes.m2‎

Lines changed: 9 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -58,9 +58,15 @@ document {
5858
-- UL {
5959
-- }
6060
-- },
61-
}
62-
}
63-
61+
LI { "useful functions involving prime numbers, submitted by Frank Schreyer:",
62+
UL {
63+
LI { TO "nextPrime", ", a simple function to find the first prime number at least as large as a given number"},
64+
LI { TO "getPrimeWithRootOfUnity", ", used to find a prime number p s.t. ZZ/p contains a n-th root of unity"},
65+
LI { TO "randomKRationalPoint", ", a function to find a random rational point on a variety over a finite field"}
66+
}
67+
}
68+
}
69+
}
6470

6571
document {
6672
Key => "list of obsolete functions",

‎M2/Macaulay2/packages/Macaulay2Doc/functions.m2‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -131,6 +131,7 @@ load "./functions/koszul-doc.m2"
131131
load "./functions/lift-doc.m2"
132132
load "./functions/liftable-doc.m2"
133133
load "./functions/map-doc.m2"
134+
load "./functions/nextPrime-doc.m2"
134135
load "./functions/replace-doc.m2"
135136
load "./functions/match-doc.m2"
136137
load "./functions/regex-doc.m2"

‎M2/Macaulay2/packages/Macaulay2Doc/functions/nextPrime-doc.m2‎

Lines changed: 157 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,157 @@
1+
doc ///
2+
Key
3+
nextPrime
4+
(nextPrime,Number)
5+
Headline
6+
compute the smallest prime greater than or equal to a given number
7+
Usage
8+
nextPrime n
9+
Inputs
10+
n: Number
11+
Outputs
12+
: ZZ
13+
the smallest prime $\ge n$
14+
Description
15+
Example
16+
nextPrime 10000
17+
nextPrime 3.5678
18+
nextPrime (3/7)
19+
SeeAlso
20+
isPrime
21+
getPrimeWithRootOfUnity
22+
///
23+
24+
doc ///
25+
Key
26+
getPrimeWithRootOfUnity
27+
(getPrimeWithRootOfUnity, ZZ,ZZ)
28+
[getPrimeWithRootOfUnity,Range]
29+
Headline
30+
find a prime p with a primitive n-th root of unity r in ZZ/p
31+
Usage
32+
(p,r)=getPrimeWithRootOfUnity(n,r1)
33+
Inputs
34+
n: ZZ
35+
the exponent of the root of unity
36+
r1: ZZ
37+
tentative root of unity
38+
Range => Sequence
39+
of integers (a,b)
40+
Outputs
41+
p: ZZ
42+
a prime in the specified range. The default range is (10^4,3*10^4)
43+
r: ZZ
44+
a primitive n-th root of unity mod p
45+
Description
46+
Text
47+
We compute the prime p as a larger prime factor of r1^n-1.
48+
If the largest p in the desired range does not work, we pass to r1+1 and repeat.
49+
Example
50+
n = 12
51+
(p,r) = getPrimeWithRootOfUnity(n,5)
52+
factor(r^n-1)
53+
r^12%p==1, r^6%p==1, r^4%p==1
54+
(p,r) = getPrimeWithRootOfUnity(12,11,Range=>(100,200))
55+
factor(r^n-1)
56+
r^12%p==1, r^6%p==1, r^4%p==1
57+
SeeAlso
58+
nextPrime
59+
Range
60+
///
61+
62+
doc ///
63+
Key
64+
Range
65+
Headline
66+
can be assigned a integral interval
67+
Usage
68+
getPrimeWithRootOfUnity(p,r,Range=>(a,b))
69+
Inputs
70+
a: ZZ
71+
b: ZZ
72+
Description
73+
Text
74+
Specifies an integral interval in which we search for a prime number with desired properties.
75+
If $b< a+2*log a$ an error message will be returned. The default value is (10^4,3*10^4)
76+
SeeAlso
77+
nextPrime
78+
getPrimeWithRootOfUnity
79+
///
80+
81+
doc ///
82+
Key
83+
randomKRationalPoint
84+
(randomKRationalPoint,Ideal)
85+
Headline
86+
Pick a random K rational point on the scheme X defined by I
87+
Usage
88+
randomKRationalPoint I
89+
Inputs
90+
I: Ideal
91+
in a polynomial ring over a finite ground field K
92+
Outputs
93+
: Ideal
94+
of a K-rational point on V(I)
95+
Description
96+
Text
97+
If X has codimension 1, then we intersect X with a randomly choosen
98+
line, and hope that the decomposition of the the intersection contains a
99+
K-rational point. If n=degree X then the probability P that this happens, is the
100+
proportion of permutations in $S_n$ with a fix point on $\{1,\ldots,n \}$,
101+
i.e. $$P=\sum_{j=1}^n (-1)^{j-1} binomial(n,j)(n-j)!/n! = 1-1/2+1/3! + \ldots $$
102+
which approachs $1-exp(-1) = 0.63\ldots$. Thus a probabilistic approach works.
103+
104+
For higher codimension we first project X birationally onto a
105+
hypersurface Y, and find a point on Y. Then we take the preimage of this point.
106+
Example
107+
p=nextPrime(random(2*10^4))
108+
kk=ZZ/p;R=kk[x_0..x_3];
109+
I=minors(4,random(R^5,R^{4:-1}));
110+
codim I, degree I
111+
time randomKRationalPoint(I)
112+
R=kk[x_0..x_5];
113+
I=minors(3,random(R^5,R^{3:-1}));
114+
codim I, degree I
115+
time randomKRationalPoint(I)
116+
Text
117+
The claim that 63 % of the intersections contain a K-rational point can be experimentally tested:
118+
Example
119+
p=10007,kk=ZZ/p,R=kk[x_0..x_2]
120+
n=5; sum(1..n,j->(-1)^(j-1)*binomial(n,j)*(n-j)!/n!)+0.0
121+
I=ideal random(n,R);
122+
time (#select(apply(100,i->(degs=apply(decompose(I+ideal random(1,R)),c->degree c);
123+
#select(degs,d->d==1))),f->f>0))
124+
///
125+
126+
127+
-- tests for nextprime
128+
TEST ///
129+
assert( nextPrime(-10) == 2)
130+
assert( nextPrime 100 == 101)
131+
assert( nextPrime 1000 == 1009)
132+
///
133+
134+
TEST ///
135+
setRandomSeed("getPrimeOfUnity")
136+
(p,r)=getPrimeWithRootOfUnity(2,3,Range=>(10^3,10^4))
137+
assert( (p,r)==(3511,-1)) -- works if the random number generator is not unchanged
138+
(p,r)=getPrimeWithRootOfUnity(15,20)
139+
assert((p,r)==(18181,21))
140+
(p,r)=getPrimeWithRootOfUnity(12,2,Range=>(100,200))
141+
assert(r^12%p==1 and r^6%p !=1 and r^4%p != 1)
142+
///
143+
144+
TEST ///
145+
setRandomSeed 0
146+
p=10007,kk=ZZ/p
147+
R=kk[x_0..x_2]
148+
I=ideal(random(4,R));
149+
time randomKRationalPoint(I)
150+
R=kk[x_0..x_4]
151+
I=minors(3,random(R^5,R^{3:-1}));
152+
codim I
153+
time pt = randomKRationalPoint(I)
154+
-- The following is the answer given the above random seed.
155+
ans = ideal(x_3+74*x_4,x_2+2336*x_4,x_1-4536*x_4,x_0-4976*x_4)
156+
assert(pt == ans)
157+
///

‎M2/Macaulay2/packages/RandomPlaneCurves.m2‎

Lines changed: 3 additions & 38 deletions
Original file line numberDiff line numberDiff line change
@@ -20,11 +20,10 @@ newPackage(
2020
DebuggingMode => false
2121
)
2222

23-
if not version#"VERSION" >= "1.4" then
24-
error "this package requires Macaulay2 version 1.4 or newer"
23+
if not version#"VERSION" >= "1.8" then
24+
error "this package requires Macaulay2 version 1.8 or newer"
2525

26-
export{"nextPrime",
27-
"distinctPlanePoints",
26+
export{"distinctPlanePoints",
2827
"constructDistinctPlanePoints",
2928
"certifyDistinctPlanePoints",
3029
"nodalPlaneCurve",
@@ -48,14 +47,6 @@ undocumented {
4847
-- (for complex numbers c this is next
4948
-- prime number of ceiling(Re(c))
5049

51-
nextPrime=method(TypicalValue=>ZZ)
52-
nextPrime Number:=n->(
53-
n0:=ceiling n;
54-
if n0 <= 2 then return 2;
55-
if even n0 then n0=n0+1;
56-
while not isPrime n0 do n0=n0+2;
57-
n0)
58-
5950

6051
-- construction of general points in the plane
6152
-- via their Hilbert-Burch matrix that occurs
@@ -154,26 +145,6 @@ imageUnderRationalMap(Ideal,Matrix):=(J,L)->(
154145

155146
beginDocumentation()
156147

157-
doc ///
158-
Key
159-
nextPrime
160-
(nextPrime,Number)
161-
Headline
162-
compute the smallest prime greater than or equal to a given number
163-
Usage
164-
nextPrime n
165-
Inputs
166-
n: Number
167-
Outputs
168-
: ZZ
169-
the smallest prime $\ge n$
170-
Description
171-
Example
172-
nextPrime 10000
173-
nextPrime 3.5678
174-
nextPrime (3/7)
175-
///
176-
177148
doc ///
178149
Key
179150
distinctPlanePoints
@@ -322,12 +293,6 @@ doc ///
322293

323294

324295
------------- TESTS --------------
325-
-- tests for nextprime
326-
TEST ///
327-
assert( nextPrime(-10) == 2)
328-
assert( nextPrime 100 == 101)
329-
assert( nextPrime 1000 == 1009)
330-
///
331296

332297
-- tests for distinct plane curves
333298
TEST ///

0 commit comments

Comments
 (0)