1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
|
/***********************************************************************
* ratpoints-2.1.2 *
* - A program to find rational points on hyperelliptic curves *
* Copyright (C) 2008, 2009 Michael Stoll *
* *
* This program is free software: you can redistribute it and/or *
* modify it under the terms of the GNU General Public License *
* as published by the Free Software Foundation, either version 2 of *
* the License, or (at your option) any later version. *
* *
* This program is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of version 2 of the GNU General *
* Public License along with this program. *
* If not, see <http://www.gnu.org/licenses/>. *
***********************************************************************/
/***********************************************************************
* gen_find_points_h.c *
* *
* This program writes the file find_points.h *
* *
* Michael Stoll, Mar 8, 2009 *
***********************************************************************/
#include "rp-private.h"
#include "primes.h"
long inv_mod_p(long p, long b)
{ /* invert b mod p */
/* doesn't have to be fast...
-- actually this will be faster in our range than the usual XGCD. */
long i = 1, n = b;
while(1)
{ if(n%p == 1) return(i);
i++;
n += b;
}
}
int main(int argc, char *argv[])
{
long n;
{ int work[RATPOINTS_MAX_PRIME];
printf("static const int "
"squares[RATPOINTS_NUM_PRIMES+1][RATPOINTS_MAX_PRIME] =\n{\n");
for(n = 0; n < RATPOINTS_NUM_PRIMES; n++)
{ long p = prime[n];
long i;
work[0] = 1;
for(i = 1; i < p; i++) work[i] = 0;
/* record non-zero squares mod p, p odd */
for(i = 1; i < p; i += 2) work[(i*i) % p] = 1;
printf("{");
for(i = 0; i < p; i++)
{ printf("%d", work[i]);
if(i < p-1) printf(",");
}
printf((n < RATPOINTS_NUM_PRIMES - 1) ? "},\n " : "}\n};\n");
}
}
printf("static const long offsets[RATPOINTS_NUM_PRIMES] =\n{");
for(n = 0; n < RATPOINTS_NUM_PRIMES; n++)
{ long p = prime[n];
{ printf("%ld", inv_mod_p(p, (2*RBA_LENGTH)%p)); }
printf((n < RATPOINTS_NUM_PRIMES - 1) ? "," : "};\n\n");
}
printf("static const long "
"inverses[RATPOINTS_NUM_PRIMES][RATPOINTS_MAX_PRIME] =\n{");
for(n = 0; n < RATPOINTS_NUM_PRIMES; n++)
{ long p = prime[n];
long i;
printf("{0");
for(i = 1; i < p; i++)
{ printf(",%ld", inv_mod_p(p, i)); }
printf((n < RATPOINTS_NUM_PRIMES - 1) ? "},\n " : "}\n};\n");
}
{ unsigned long work[RATPOINTS_MAX_PRIME];
printf("unsigned long "
"sieves0[RATPOINTS_NUM_PRIMES][2*RATPOINTS_MAX_PRIME_EVEN] =\n{\n");
for(n = 0; n < RATPOINTS_NUM_PRIMES; n++)
{ long p = prime[n];
long i;
for(i = 0; i < p; i++) work[i] = ~0UL;
for(i = 0; i < LONG_LENGTH; i++)
{ work[(p*i)>>LONG_SHIFT] &= ~(1UL<<((p*i) & LONG_MASK)); }
printf("{");
for(i = 0; i < p; i++)
{ printf("0x%*.*lx, ", 16, 16, work[i]); }
for(i = 0; i < p; i++)
{ printf("0x%*.*lx", 16, 16, work[i]);
if(i < p-1) printf(", ");
}
printf((n < RATPOINTS_NUM_PRIMES - 1) ? "},\n" : "}\n");
}
printf("};\n\n");
}
return(0);
}
|
