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
| // This is an old function that was used before the port to GMP
// It is here to show how the GCD function below used to look.
// Unfortunately I lost the original version of it after the port.
// The GMP really makes the code alot harder to read
long long powmod(long long base, long long power, long long mod) {
//2^7 = 2^2 * 2^2 * 2^2 * 2^1
long long result = 1;
long long i;
for(i=0; i < power/2; i++) {
result *= base*base%mod;
}
//odd power?
if (power%2 != 0)
result *= base;
return result%mod;
}
//GMP Port:
//data structures
typedef struct gcdNode {
mpz_t j;
mpz_t k;
mpz_t q;
mpz_t r;
mpz_t x;
mpz_t y;
struct gcdNode* next;
} gcdNode;
typedef struct gcdReturn {
mpz_t mult_inv;
mpz_t gcd;
} gcdReturn;
//function
gcdReturn* gcdCalc(mpz_t j, mpz_t k) {
gcdNode* help_ptr = (gcdNode*)malloc(sizeof(gcdNode));
mpz_init_set(help_ptr->j, j);
mpz_init_set(help_ptr->k, k);
mpz_inits(help_ptr->q, help_ptr->r, help_ptr->x, help_ptr->y, NULL);
help_ptr->next = NULL;
//calculate gcd by using euclids division algorithm
while (1) {
//multiplier and remainder
mpz_divexact(help_ptr->q, help_ptr->k, help_ptr->j);
mpz_mod(help_ptr->r, help_ptr->k, help_ptr->j);
//remainder 0?
if (mpz_cmp_ui(help_ptr->r, 0) == 0) break;
//create next row
gcdNode* new_ptr = (gcdNode*)malloc(sizeof(gcdNode));
mpz_inits(new_ptr->j, new_ptr->k, new_ptr->q, new_ptr->r, new_ptr->x, new_ptr->y, NULL);
mpz_set(new_ptr->j, help_ptr->r);
mpz_set(new_ptr->k, help_ptr->j);
//link and traverse
new_ptr->next = help_ptr;
help_ptr = new_ptr;
}//gcd calc loop
//save gcd
gcdReturn* returnVal = (gcdReturn*) malloc(sizeof(gcdReturn));
mpz_init_set(gcdReturn->gcd, help_ptr->j);
//set intial x & y values for last row in table
mpz_set_ui(help_ptr->x, 1);
mpz_set_ui(help_ptr->y, 0);
//column headings
dbg("j\tk\tq\tr\tx\ty");
//fill in the x & y columns of the table
while (help_ptr->next != NULL) {
gmp_printf("%Zd\t%Zd\t%Zd\t%Zd\t%Zd\t%Zd\n", help_ptr->j, help_ptr->k, help_ptr->q, help_ptr->r, help_ptr->x, help_ptr->y);
gcdNode* next_ptr = help_ptr->next;
mpz_t numerator, kx, one;
mpz_inits(numerator, kx, one, NULL);
mpz_set_ui(one, 1);
mpz_mul(kx, next_ptr->k, help_ptr->x);
mpz_sub(numerator, one, kx);
mpz_divexact(next_ptr->x, numerator, next_ptr->j);
mpz_set(next_ptr->y, help_ptr->x);
//cleanup
mpz_clears(numerator, kx, one, help_ptr->j, help_ptr->k, help_ptr->q, help_ptr->r, help_ptr->x, help_ptr->y, NULL);
free(help_ptr);
help_ptr = next_ptr;
}
gmp_printf("%Zd\t%Zd\t%Zd\t%Zd\t%Zd\t%Zd\n", help_ptr->j, help_ptr->k, help_ptr->q, help_ptr->r, help_ptr->x, help_ptr->y);
//set gcd
//take care of negative multiplicative inverse
mpz_init(returnVal->mult_inv);
if (mpz_sgn (help_ptr->x) < 0)
mpz_add(returnVal->mult_inv, help_ptr->x, help_ptr->k);
else
mpz_set(returnVal->mult_inv, help_ptr->x);
return returnVal;
}
|