1 #include "cado.h" // IWYU pragma: keep
2 #include <stdio.h>
3 #include <string.h>
4 #include "mod_ul.h" // modulusul_t
5 #include "galois_utils.h"
6 #include "macros.h" // ASSERT_ALWAYS
7
automorphism_init(int * ord,int mat[4],const char * galois_autom)8 void automorphism_init(int *ord, int mat[4], const char *galois_autom){
9 int A, B, C, D; // x -> (A*x+B)/(C*x+D)
10
11 if(strcmp(galois_autom, "autom2.1") == 0
12 || strcmp(galois_autom, "autom2.1g") == 0
13 || strcmp(galois_autom, "1/y") == 0){
14 *ord = 2; A = 0; B = 1; C = 1; D = 0;
15 }
16 else if(strcmp(galois_autom, "autom2.2") == 0
17 || strcmp(galois_autom, "autom2.2g") == 0
18 || strcmp(galois_autom, "_y") == 0){
19 *ord = 2; A = -1; B = 0; C = 0; D = 1;
20 }
21 else if(strcmp(galois_autom, "autom3.1") == 0
22 || strcmp(galois_autom, "autom3.1g") == 0){
23 // 1-1/x = (x-1)/x
24 *ord = 3; A = 1; B = -1; C = 1; D = 0;
25 }
26 else if(strcmp(galois_autom, "autom3.2") == 0){
27 // -1-1/x = (-x-1)/x
28 *ord = 3; A = -1; B = -1; C = 1; D = 0;
29 }
30 else if(strcmp(galois_autom, "autom4.1") == 0){
31 // -(x+1)/(x-1)
32 *ord = 4; A = -1; B = -1; C = 1; D = -1;
33 }
34 else if(strcmp(galois_autom, "autom6.1") == 0){
35 // -(2*x+1)/(x-1)
36 *ord = 6; A = -2; B = -1; C = 1; D = -1;
37 }
38 else{
39 fprintf(stderr, "Unknown automorphism: %s\n", galois_autom);
40 ASSERT_ALWAYS(0);
41 }
42 mat[0] = A; mat[1] = B; mat[2] = C; mat[3] = D;
43 }
44
45 /* OUTPUT: sigma(r).
46 Since we can have r == qq, we need unsigned long's.
47 */
automorphism_apply(residueul_t mat[4],unsigned long r,const modulusul_t mm,const unsigned long qq)48 unsigned long automorphism_apply(residueul_t mat[4], unsigned long r,
49 const modulusul_t mm, const unsigned long qq)
50 {
51 residueul_t xx;
52 unsigned long sigma_r;
53
54 modul_init(xx, mm);
55 if(r == qq){
56 // sigma(oo) = A/C
57 // TODO: test C w.r.t. +/- 1
58 modul_inv(xx, mat[2], mm);
59 modul_mul(xx, xx, mat[0], mm);
60 sigma_r = modul_get_ul(xx, mm);
61 }
62 else{
63 residueul_t rr;
64
65 modul_init(rr, mm);
66 modul_set_ul(rr, r, mm);
67 // denominator: C*r+D
68 modul_mul(xx, mat[2], rr, mm);
69 modul_add(xx, xx, mat[3], mm);
70 if(modul_is0(xx, mm))
71 // 1/0 = oo
72 sigma_r = qq;
73 else{
74 residueul_t yy;
75
76 modul_init(yy, mm);
77 modul_inv(yy, xx, mm);
78 // numerator: A*r+B
79 modul_mul(xx, mat[0], rr, mm);
80 modul_add(xx, xx, mat[1], mm);
81 modul_mul(yy, yy, xx, mm);
82 sigma_r = modul_get_ul(yy, mm);
83 modul_clear(yy, mm);
84 }
85 modul_clear(rr, mm);
86 }
87 modul_clear(xx, mm);
88 return sigma_r;
89 }
90