1 /*
2 This is an interface to masquerade ordinary long integers in C
3 language behind an interface mimicing arbitrary precision arithmetic
4 interger (mpz) interface in The GNU MP Bignum library.
5
6 Only the functions needed by Kbipack are available.
7
8 THIS INTERFACE IS NOT PERFORMING ANY ARBITRARY PRECISION ARITHMETIC.
9 You should always use the GMP library (http://gmplib.org/) when possible.
10
11 Some basic techniques to detect integer overflows are implemented.
12 However, they don't interfere the program flow IN ANY WAY, they just write
13 an error message to stdout.
14
15 */
16
17 #include "mpz.h"
18 #include "GmshConfig.h"
19
20 #if ! defined(HAVE_GMP)
21
22 #include "GmshMessage.h"
23 #include "limits.h"
24
overflow()25 void overflow()
26 {
27 printf("ERROR: Integer overflow detected! Compile with GMP library to fix this.\n");
28 Msg::Error("Integer overflow detected! Compile with GMP library to fix this.");
29 }
30
addcheck(long int a,long int b)31 long int addcheck(long int a, long int b){
32 long int c = a + b;
33 if (b >= 0){
34 if (c < a) overflow();
35 }
36 else if (b < 0){
37 if (c > a) overflow();
38 }
39 return c;
40 }
41
multcheck(long int a,long int b)42 long int multcheck(long int a, long int b){
43 long long int c = (long long int)a*b;
44 if(c > LONG_MAX || c < LONG_MIN) overflow();
45 return c;
46 }
47
mpz_init(mpz_ptr x)48 void mpz_init(mpz_ptr x)
49 {
50 x->z = 0;
51 }
52
mpz_init_set(mpz_ptr rop,mpz_ptr op)53 void mpz_init_set(mpz_ptr rop, mpz_ptr op)
54 {
55 rop->z = op->z;
56 }
57
mpz_init_set_si(mpz_ptr rop,signed long int op)58 void mpz_init_set_si(mpz_ptr rop, signed long int op)
59 {
60 rop->z = op;
61 }
62
mpz_set_ui(mpz_ptr rop,unsigned long int op)63 void mpz_set_ui(mpz_ptr rop, unsigned long int op)
64 {
65 rop->z = op;
66 }
67
mpz_set(mpz_ptr rop,mpz_ptr op)68 void mpz_set(mpz_ptr rop, mpz_ptr op)
69 {
70 rop->z = op->z;
71 }
mpz_set_si(mpz_ptr rop,signed long int op)72 void mpz_set_si(mpz_ptr rop, signed long int op)
73 {
74 rop->z = op;
75 }
76
mpz_clear(mpz_ptr x)77 void mpz_clear(mpz_ptr x)
78 {
79
80 }
81
82 // arithmethic
mpz_add(mpz_ptr rop,mpz_ptr op1,mpz_ptr op2)83 void mpz_add(mpz_ptr rop, mpz_ptr op1, mpz_ptr op2)
84 {
85 rop->z = addcheck(op1->z, op2->z);
86 }
87
mpz_mul(mpz_ptr rop,mpz_ptr op1,mpz_ptr op2)88 void mpz_mul(mpz_ptr rop, mpz_ptr op1, mpz_ptr op2)
89 {
90 rop->z = multcheck(op1->z, op2->z);
91 }
92
mpz_addmul(mpz_ptr rop,mpz_ptr op1,mpz_ptr op2)93 void mpz_addmul(mpz_ptr rop, mpz_ptr op1, mpz_ptr op2)
94 {
95 rop->z = addcheck(rop->z,multcheck(op1->z,op2->z));
96 }
97
mpz_neg(mpz_ptr rop,mpz_ptr op)98 void mpz_neg(mpz_ptr rop, mpz_ptr op)
99 {
100 rop->z = -op->z;
101 }
102
103
104 // division
mpz_divexact(mpz_ptr q,mpz_ptr n,mpz_ptr d)105 void mpz_divexact(mpz_ptr q, mpz_ptr n, mpz_ptr d)
106 {
107 ldiv_t temp;
108 temp = ldiv(n->z, d->z);
109 q->z = temp.quot;
110 }
111
mpz_cdiv_q(mpz_ptr q,mpz_ptr n,mpz_ptr d)112 void mpz_cdiv_q(mpz_ptr q, mpz_ptr n, mpz_ptr d)
113 {
114 ldiv_t temp;
115 temp = ldiv(n->z, d->z);
116 q->z = temp.quot;
117 }
mpz_cdiv_qr(mpz_ptr q,mpz_ptr r,mpz_ptr n,mpz_ptr d)118 void mpz_cdiv_qr(mpz_ptr q, mpz_ptr r, mpz_ptr n, mpz_ptr d)
119 {
120 ldiv_t temp;
121 temp = ldiv(n->z, d->z);
122 q->z = temp.quot;
123 r->z = temp.rem;
124 }
125
mpz_tdiv_r(mpz_ptr r,mpz_ptr n,mpz_ptr d)126 void mpz_tdiv_r(mpz_ptr r, mpz_ptr n, mpz_ptr d)
127 {
128 ldiv_t temp;
129 temp = ldiv(n->z, d->z);
130 if(n->z < 0) r->z = -temp.rem;
131 else r->z = temp.rem;
132 }
133
134 // compare
mpz_cmp_si(mpz_ptr op1,signed long int op2)135 int mpz_cmp_si(mpz_ptr op1, signed long int op2)
136 {
137 if(op1->z > op2) return 1;
138 else if(op1->z < op2) return -1;
139 else return 0;
140 return 0;
141 }
142
mpz_cmpabs(mpz_ptr op1,mpz_ptr op2)143 int mpz_cmpabs(mpz_ptr op1, mpz_ptr op2)
144 {
145 if(abs(op1->z) > abs(op2->z)) return 1;
146 else if(abs(op1->z) < abs(op2->z)) return -1;
147 else return 0;
148 return 0;
149 }
150
mpz_sgn(mpz_ptr op)151 int mpz_sgn(mpz_ptr op)
152 {
153 if(op->z > 0) return 1;
154 else if(op->z < 0) return -1;
155 else return 0;
156 }
157
mpz_get_si(mpz_ptr op)158 signed long int mpz_get_si(mpz_ptr op)
159 {
160 signed long int temp = op->z;
161 return temp;
162 }
163
gcd(long int a,long int b)164 long int gcd(long int a, long int b)
165 {
166 long int temp = 0;
167 long int remainder = 0;
168 if (a < b) {
169 temp = a;
170 a = b;
171 b = temp;
172 }
173 remainder = a % b;
174 if(remainder == 0) return b;
175 else return gcd(b, remainder);
176 }
177
extended_gcd(long int * g,long int * s,long int * t,long int a,long int b)178 void extended_gcd(long int* g, long int* s, long int* t,
179 long int a, long int b)
180 {
181 long int x = 0;
182 long int lastx = 1;
183 long int y = 1;
184 long int lasty = 0;
185 while (b != 0){
186 ldiv_t divt = ldiv(a,b);
187 long int quotient = divt.quot;
188
189 long int temp = b;
190 b = a % b;
191 a = temp;
192
193 temp = x;
194 x = addcheck(lastx, multcheck(-quotient,x));
195 lastx = temp;
196
197 temp = y;
198 y = addcheck(lasty, multcheck(-quotient,y));
199 lasty = temp;
200 }
201 *s = lastx;
202 *t = lasty;
203 *g = a;
204 }
205
mpz_gcdext(mpz_ptr g,mpz_ptr s,mpz_ptr t,mpz_ptr a,mpz_ptr b)206 void mpz_gcdext(mpz_ptr g, mpz_ptr s, mpz_ptr t, mpz_ptr a, mpz_ptr b)
207 {
208
209 extended_gcd(&g->z, &s->z, &t->z, a->z, b->z);
210 /*printf("g: %ld, s: %ld, t: %ld, a: %ld, b: %ld. \n",
211 g->z, s->z, t->z, a->z, b->z); */
212 }
213
mpz_out_str(FILE * stream,int base,mpz_ptr op)214 size_t mpz_out_str(FILE *stream, int base, mpz_ptr op)
215 {
216 return 0;
217 }
218
219 /*
220 main()
221 {
222
223 mpz_t a;
224 mpz_t b;
225 mpz_t c;
226 mpz_init(a);
227 mpz_init(b);
228 mpz_init(c);
229
230 mpz_set_si(a, 2147483647);
231 mpz_set_si(b, 2);
232 mpz_mul(c, a, b);
233
234 printf("a: %ld, b: %ld, c: %ld. \n", a->z, b->z, c->z);
235
236 return EXIT_SUCCESS;
237
238 }*/
239 #endif
240