File: gen_find_points_h.c

package info (click to toggle)
ratpoints 1:2.1.3-1
  • links: PTS, VCS
  • area: main
  • in suites: bullseye, buster, jessie, jessie-kfreebsd, sid, stretch
  • size: 572 kB
  • ctags: 230
  • sloc: ansic: 4,201; makefile: 155
file content (115 lines) | stat: -rw-r--r-- 4,145 bytes parent folder | download
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);
}