1 # Division ganzer Zahlen
2
3 # Dividiert zwei Unsigned Digit sequences durcheinander.
4 # UDS_divide(a_MSDptr,a_len,a_LSDptr, b_MSDptr,b_len,b_LSDptr, &q,&r);
5 # Die UDS a = a_MSDptr/a_len/a_LSDptr (a>=0) wird durch
6 # die UDS b = b_MSDptr/b_len/b_LSDptr (b>=0) dividiert:
7 # a = q * b + r mit 0 <= r < b. Bei b=0 Error.
8 # q der Quotient, r der Rest.
9 # q = q_MSDptr/q_len/q_LSDptr, r = r_MSDptr/r_len/r_LSDptr beides
10 # Normalized Unsigned Digit sequences.
11 # Vorsicht: q_LSDptr <= r_MSDptr,
12 # Vorzeichenerweiterung von r kann q zerstören!
13 # Vorzeichenerweiterung von q ist erlaubt.
14 # a und b werden nicht modifiziert.
15 # num_stack wird erniedrigt.
16 #define UDS_divide(a_MSDptr,a_len,a_LSDptr,b_MSDptr,b_len,b_LSDptr,q_,r_) \
17 { # Platz fürs Ergebnis machen. Brauche maximal a_len+1 Digits. \
18 var uintC _a_len = (a_len); \
19 var uintD* roomptr; num_stack_need_1(_a_len+1,roomptr=,); \
20 UDS_divide_(a_MSDptr,_a_len,a_LSDptr,b_MSDptr,b_len,b_LSDptr,roomptr,q_,r_); \
21 }
22 local void UDS_divide_ (uintD* a_MSDptr, uintC a_len, uintD* a_LSDptr,
23 uintD* b_MSDptr, uintC b_len, uintD* b_LSDptr,
24 uintD* roomptr, DS* q_, DS* r_);
25 # Methode:
26 # erst a und b normalisieren: a=[a[m-1],...,a[0]], b=[b[n-1],...,b[0]]
27 # mit m>=0 und n>0 (Stellensystem der Basis beta=2^intDsize).
28 # Falls m<n, ist q:=0 und r:=a.
29 # Falls m>=n=1, Single-Precision-Division:
30 # r:=0, j:=m,
31 # while j>0 do
32 # {Hier (q[m-1]*beta^(m-1)+...+q[j]*beta^j) * b[0] + r*beta^j =
33 # = a[m-1]*beta^(m-1)+...+a[j]*beta^j und 0<=r<b[0]<beta}
34 # j:=j-1, r:=r*beta+a[j], q[j]:=floor(r/b[0]), r:=r-b[0]*q[j].
35 # Normalisiere [q[m-1],...,q[0]], liefert q.
36 # Falls m>=n>1, Multiple-Precision-Division:
37 # Es gilt a/b < beta^(m-n+1).
38 # s:=intDsize-1-(Nummer des höchsten Bits in b[n-1]), 0<=s<intDsize.
39 # Schiebe a und b um s Bits nach links und kopiere sie dabei. r:=a.
40 # r=[r[m],...,r[0]], b=[b[n-1],...,b[0]] mit b[n-1]>=beta/2.
41 # Für j=m-n,...,0: {Hier 0 <= r < b*beta^(j+1).}
42 # Berechne q* :
43 # q* := floor((r[j+n]*beta+r[j+n-1])/b[n-1]).
44 # Bei Überlauf (q* >= beta) setze q* := beta-1.
45 # Berechne c2 := ((r[j+n]*beta+r[j+n-1]) - q* * b[n-1])*beta + r[j+n-2]
46 # und c3 := b[n-2] * q*.
47 # {Es ist 0 <= c2 < 2*beta^2, sogar 0 <= c2 < beta^2 falls kein
48 # Überlauf aufgetreten war. Ferner 0 <= c3 < beta^2.
49 # Bei Überlauf und r[j+n]*beta+r[j+n-1] - q* * b[n-1] >= beta,
50 # das heißt c2 >= beta^2, kann man die nächste Abfrage überspringen.}
51 # Solange c3 > c2, {hier 0 <= c2 < c3 < beta^2} setze
52 # q* := q* - 1, c2 := c2 + b[n-1]*beta, c3 := c3 - b[n-2].
53 # Falls q* > 0:
54 # Setze r := r - b * q* * beta^j, im einzelnen:
55 # [r[n+j],...,r[j]] := [r[n+j],...,r[j]] - q* * [b[n-1],...,b[0]].
56 # also: u:=0, for i:=0 to n-1 do
57 # u := u + q* * b[i],
58 # r[j+i]:=r[j+i]-(u mod beta) (+ beta, falls Carry),
59 # u:=u div beta (+ 1, falls bei der Subtraktion Carry)
60 # r[n+j]:=r[n+j]-u.
61 # {Da stets u = (q* * [b[i-1],...,b[0]] div beta^i) + 1
62 # < q* + 1 <= beta, läuft der Übertrag u nicht über.}
63 # Tritt dabei ein negativer Übertrag auf, so setze q* := q* - 1
64 # und [r[n+j],...,r[j]] := [r[n+j],...,r[j]] + [0,b[n-1],...,b[0]].
65 # Setze q[j] := q*.
66 # Normalisiere [q[m-n],..,q[0]] und erhalte den Quotienten q,
67 # Schiebe [r[n-1],...,r[0]] um s Bits nach rechts, normalisiere und
68 # erhalte den Rest r.
69 # Dabei kann q[j] auf dem Platz von r[n+j] liegen.
70 local void UDS_divide_ (uintD* a_MSDptr, uintC a_len, uintD* a_LSDptr,
71 uintD* b_MSDptr, uintC b_len, uintD* b_LSDptr,
72 uintD* roomptr, # ab roomptr kommen a_len+1 freie Digits
73 DS* q_, DS* r_)
74 {
75 # a normalisieren (a_MSDptr erhöhen, a_len erniedrigen):
76 while ((a_len>0) && (a_MSDptr[0]==0)) {
77 a_MSDptr++; a_len--;
78 }
79 # b normalisieren (b_MSDptr erhöhen, b_len erniedrigen):
80 while (1) {
81 if (b_len==0)
82 divide_0();
83 if (b_MSDptr[0]==0) {
84 b_MSDptr++; b_len--;
85 } else
86 break;
87 }
88 # jetzt m=a_len >=0 und n=b_len >0.
89 if (a_len < b_len) {
90 # m<n: Trivialfall, q=0, r:= Kopie von a.
91 var uintD* r_MSDptr = roomptr;
92 var uintD* r_LSDptr = &roomptr[a_len];
93 # Speicheraufbau: r_MSDptr/0/r_MSDptr/a_len/r_LSDptr
94 # | q | r |
95 copy_loop_down(a_LSDptr,r_LSDptr,a_len);
96 q_->MSDptr = r_MSDptr; q_->len = 0; q_->LSDptr = r_MSDptr; # q = 0, eine NUDS
97 r_->MSDptr = r_MSDptr; r_->len = a_len; r_->LSDptr = r_LSDptr; # r = Kopie von a, eine NUDS
98 return;
99 } else if (b_len==1) {
100 # n=1: Single-Precision-Division
101 # beta^(m-1) <= a < beta^m ==> beta^(m-2) <= a/b < beta^m
102 var uintD* q_MSDptr = roomptr;
103 var uintD* q_LSDptr = &roomptr[a_len];
104 var uintD* r_MSDptr = q_LSDptr;
105 var uintD* r_LSDptr = &r_MSDptr[1];
106 # Speicheraufbau: q_MSDptr/a_len/q_LSDptr r_MSDptr/1/r_LSDptr
107 # | q | | r |
108 var uintD rest = divucopy_loop_up(b_MSDptr[0],a_MSDptr,q_MSDptr,a_len); # Division durch b[0]
109 var uintC r_len;
110 if (!(rest==0)) {
111 r_MSDptr[0] = rest; r_len=1; # Rest als r ablegen
112 } else {
113 r_MSDptr = r_LSDptr; r_len=0; # Rest auf 0 normalisieren
114 }
115 if (q_MSDptr[0]==0) {
116 q_MSDptr++; a_len--; # q normalisieren
117 }
118 q_->MSDptr = q_MSDptr; q_->len = a_len; q_->LSDptr = q_LSDptr; # q ablegen
119 r_->MSDptr = r_MSDptr; r_->len = r_len; r_->LSDptr = r_LSDptr; # r ablegen
120 return;
121 } else {
122 # n>1: Multiple-Precision-Division
123 # beta^(m-1) <= a < beta^m, beta^(n-1) <= b < beta^n ==>
124 # beta^(m-n-1) <= a/b < beta^(m-n+1).
125 var uintL s;
126 SAVE_NUM_STACK # num_stack retten
127 # s bestimmen:
128 {
129 var uintD msd = b_MSDptr[0]; # b[n-1]
130 #if 0
131 s = 0;
132 while ((sintD)msd >= 0) {
133 msd = msd<<1; s++;
134 }
135 #else # ein wenig effizienter, Abfrage auf s=0 vorwegnehmen
136 if ((sintD)msd < 0) {
137 s = 0; goto shift_ok;
138 } else {
139 integerlengthD(msd, s = intDsize - ); goto shift;
140 }
141 #endif
142 }
143 # 0 <= s < intDsize.
144 # Kopiere b und schiebe es dabei um s Bits nach links:
145 if (!(s==0))
146 shift: {
147 var uintD* old_b_LSDptr = b_LSDptr;
148 num_stack_need(b_len,b_MSDptr=,b_LSDptr=);
149 shiftleftcopy_loop_down(old_b_LSDptr,b_LSDptr,b_len,s);
150 }
151 shift_ok: {
152 # Wieder b = b_MSDptr/b_len/b_LSDptr.
153 # Kopiere a und schiebe es dabei um s Bits nach links, erhalte r:
154 var uintD* r_MSDptr = roomptr;
155 var uintD* r_LSDptr = &roomptr[a_len+1];
156 # Speicheraufbau: r_MSDptr/ a_len+1 /r_LSDptr
157 # | r |
158 # später: q_MSDptr/a_len-b_len+1/r_MSDptr/b_len/r_LSDptr
159 # | q | r |
160 if (s==0) {
161 copy_loop_down(a_LSDptr,r_LSDptr,a_len); r_MSDptr[0] = 0;
162 } else {
163 r_MSDptr[0] = shiftleftcopy_loop_down(a_LSDptr,r_LSDptr,a_len,s);
164 }
165 # Nun r = r_MSDptr/a_len+1/r_LSDptr.
166 var uintC j = a_len-b_len; # m-n
167 var uintD* r_ptr = &r_LSDptr[-(uintP)j]; # Pointer oberhalb von r[j]
168 var uintD* q_MSDptr = r_MSDptr;
169 var uintC q_len = j = j+1; # q wird m-n+1 Digits haben
170 var uintD b_msd = b_MSDptr[0]; # b[n-1]
171 var uintD b_2msd = b_MSDptr[1]; # b[n-2]
172 #if HAVE_DD
173 var uintDD b_msdd = highlowDD(b_msd,b_2msd); # b[n-1]*beta+b[n-2]
174 #endif
175 # Divisions-Schleife: (wird m-n+1 mal durchlaufen)
176 # j = Herabzähler, b_MSDptr/b_len/b_LSDptr = [b[n-1],...,b[0]], b_len=n,
177 # r_MSDptr = Pointer auf r[n+j] = Pointer auf q[j],
178 # r_ptr = Pointer oberhalb von r[j].
179 do {
180 var uintD q_star;
181 var uintD c1;
182 if (r_MSDptr[0] < b_msd) { # r[j+n] < b[n-1] ?
183 # Dividiere r[j+n]*beta+r[j+n-1] durch b[n-1], ohne Überlauf:
184 #if HAVE_DD
185 divuD(highlowDD(r_MSDptr[0],r_MSDptr[1]),b_msd, q_star=,c1=);
186 #else
187 divuD(r_MSDptr[0],r_MSDptr[1],b_msd, q_star=,c1=);
188 #endif
189 } else {
190 # Überlauf, also r[j+n]*beta+r[j+n-1] >= beta*b[n-1]
191 q_star = bitm(intDsize)-1; # q* = beta-1
192 # Teste ob r[j+n]*beta+r[j+n-1] - (beta-1)*b[n-1] >= beta
193 # <==> r[j+n]*beta+r[j+n-1] + b[n-1] >= beta*b[n-1]+beta
194 # <==> b[n-1] < floor((r[j+n]*beta+r[j+n-1]+b[n-1])/beta) {<= beta !} ist.
195 # Wenn ja, direkt zur Subtraktionschleife.
196 # (Andernfalls ist r[j+n]*beta+r[j+n-1] - (beta-1)*b[n-1] < beta
197 # <==> floor((r[j+n]*beta+r[j+n-1]+b[n-1])/beta) = b[n-1] ).
198 if ((r_MSDptr[0] > b_msd) || ((c1 = r_MSDptr[1]+b_msd) < b_msd)) {
199 # r[j+n] >= b[n-1]+1 oder
200 # r[j+n] = b[n-1] und Addition r[j+n-1]+b[n-1] gibt Carry ?
201 goto subtract; # ja -> direkt in die Subtraktion
202 }
203 }
204 # q_star = q*,
205 # c1 = (r[j+n]*beta+r[j+n-1]) - q* * b[n-1] (>=0, <beta).
206 #if HAVE_DD
207 {
208 var uintDD c2 = highlowDD(c1,r_MSDptr[2]); # c1*beta+r[j+n-2]
209 var uintDD c3 = muluD(b_2msd,q_star); # b[n-2] * q*
210 # Solange c2 < c3, c2 erhöhen, c3 erniedrigen:
211 # Rechne dabei mit c3-c2:
212 # solange >0, um b[n-1]*beta+b[n-2] erniedrigen.
213 # Dies kann wegen b[n-1]*beta+b[n-2] >= beta^2/2
214 # höchstens zwei mal auftreten.
215 if (c3 > c2) {
216 q_star = q_star-1; # q* := q* - 1
217 if (c3-c2 > b_msdd)
218 q_star = q_star-1; # q* := q* - 1
219 }
220 }
221 #else
222 # Wie oben, nur mit zweigeteilten c2=[c2hi|c2lo] und c3=[c3hi|c3lo]:
223 #define c2hi c1
224 {
225 var uintD c2lo = r_MSDptr[2]; # c2hi*beta+c2lo = c1*beta+r[j+n-2]
226 var uintD c3hi;
227 var uintD c3lo;
228 muluD(b_2msd,q_star, c3hi=,c3lo=); # c3hi*beta+c3lo = b[n-2] * q*
229 if ((c3hi > c2hi) || ((c3hi == c2hi) && (c3lo > c2lo))) {
230 q_star = q_star-1; # q* := q* - 1
231 c3hi -= c2hi; if (c3lo < c2lo) c3hi--; c3lo -= c2lo; # c3 := c3-c2
232 if ((c3hi > b_msd) || ((c3hi == b_msd) && (c3lo > b_2msd)))
233 q_star = q_star-1; # q* := q* - 1
234 }
235 }
236 #undef c2hi
237 #endif
238 if (!(q_star==0))
239 subtract:
240 {
241 # Subtraktionsschleife: r := r - b * q* * beta^j
242 var uintD carry = mulusub_loop_down(q_star,b_LSDptr,r_ptr,b_len);
243 # Noch r_ptr[-b_len-1] -= carry, d.h. r_MSDptr[0] -= carry
244 # durchführen und danach r_MSDptr[0] vergessen:
245 if (carry > r_MSDptr[0]) {
246 # Subtraktion ergab Übertrag
247 q_star = q_star-1; # q* := q* - 1
248 addto_loop_down(b_LSDptr,r_ptr,b_len); # Additionsschleife
249 # r[n+j] samt Carry kann vergessen werden...
250 }
251 }
252 # Berechnung von q* ist fertig.
253 *r_MSDptr++ = q_star; # als q[j] ablegen
254 r_ptr++;
255 } while (--j != 0);
256 # Nun ist q = [q[m-n],..,q[0]] = q_MSDptr/q_len/r_MSDptr
257 # und r = [r[n-1],...,r[0]] = r_MSDptr/b_len/r_LSDptr.
258 # q normalisieren und ablegen:
259 if (q_MSDptr[0]==0) {
260 q_MSDptr++; q_len--;
261 }
262 q_->MSDptr = q_MSDptr; q_->len = q_len; q_->LSDptr = r_MSDptr;
263 # Schiebe [r[n-1],...,r[0]] um s Bits nach rechts:
264 if (!(s==0)) {
265 shiftright_loop_up(r_MSDptr,b_len,s);
266 }
267 # r normalisieren und ablegen:
268 while ((b_len>0) && (r_MSDptr[0]==0)) {
269 r_MSDptr++; b_len--;
270 }
271 r_->MSDptr = r_MSDptr; r_->len = b_len; r_->LSDptr = r_LSDptr;
272 RESTORE_NUM_STACK # num_stack zurück
273 return;
274 }}
275 }
276
277 # Dividiert zwei Integers x,y >=0 und liefert Quotient und Rest
278 # der Division x/y. Bei y=0 Error.
279 # I_I_divide_I_I(x,y);
280 # > x,y: Integers >=0
281 # < STACK_1: Quotient q
282 # < STACK_0: Rest r
283 # Erniedrigt STACK um 2
284 # can trigger GC
I_I_divide_I_I(object x,object y)285 local maygc void I_I_divide_I_I (object x, object y)
286 {
287 if (I_fixnump(x)) {
288 # x Fixnum >=0
289 if (I_fixnump(y)) {
290 # auch y Fixnum >=0
291 var uintV x_ = posfixnum_to_V(x);
292 var uintV y_ = posfixnum_to_V(y);
293 if (y_==0) {
294 divide_0();
295 } else if (x_ < y_) {
296 # Trivialfall: q=0, r=x
297 goto trivial;
298 } else {
299 #if (intVsize>32)
300 if (x_ >= vbit(32)) {
301 if (y_ < vbit(32)) {
302 # 64-durch-32-Bit-Division
303 var uint64 q;
304 var uint32 r;
305 divu_6432_6432(x_,y_,q=,r=);
306 pushSTACK(UQ_to_I(q));
307 pushSTACK(UL_to_I(r));
308 } else {
309 # volle 64-durch-64-Bit-Division
310 var uint64 q;
311 var uint64 r;
312 divu_6464_6464(x_,y_,q=,r=);
313 pushSTACK(UQ_to_I(q));
314 pushSTACK(UQ_to_I(r));
315 }
316 } else
317 #endif
318 {
319 if (y_ < bit(16)) {
320 # 32-durch-16-Bit-Division
321 var uint32 q;
322 var uint16 r;
323 divu_3216_3216(x_,y_,q=,r=);
324 pushSTACK(UL_to_I(q));
325 pushSTACK(fixnum((uintL)r));
326 } else {
327 # volle 32-durch-32-Bit-Division
328 var uint32 q;
329 var uint32 r;
330 divu_3232_3232(x_,y_,q=,r=);
331 pushSTACK(UL_to_I(q));
332 pushSTACK(UL_to_I(r));
333 }
334 }
335 }
336 } else {
337 # y Bignum >0
338 trivial:
339 # Trivialfall: q=0, r=x
340 pushSTACK(Fixnum_0);
341 pushSTACK(x);
342 }
343 } else {
344 # x Bignum -> allgemeine Division:
345 var uintD* x_MSDptr;
346 var uintC x_len;
347 var uintD* x_LSDptr;
348 var uintD* y_MSDptr;
349 var uintC y_len;
350 var uintD* y_LSDptr;
351 # x in NDS umwandeln, als UDS auffassen:
352 BN_to_NDS_nocopy(x, x_MSDptr=,x_len=,x_LSDptr=);
353 {
354 # y in NDS umwandeln, als UDS auffassen:
355 I_to_NDS_nocopy(y, y_MSDptr=,y_len=,y_LSDptr=);
356 # dividieren:
357 {
358 SAVE_NUM_STACK # num_stack retten
359 var DS q;
360 var DS r;
361 begin_arith_call();
362 UDS_divide(x_MSDptr,x_len,x_LSDptr,y_MSDptr,y_len,y_LSDptr, &q,&r);
363 end_arith_call();
364 # q in Integer umwandeln:
365 pushSTACK(NUDS_to_I(q.MSDptr,q.len));
366 # r in Integer umwandeln (jetzt erst, nachdem q verwertet ist!):
367 pushSTACK(NUDS_to_I(r.MSDptr,r.len));
368 RESTORE_NUM_STACK # num_stack zurück
369 }
370 }
371 }
372 }
373
374 /* Error, when the quotient is not an integer
375 > STACK_1: numerator x
376 > STACK_0: denominator y */
error_exquo(void)377 local _Noreturn void error_exquo (void) {
378 pushSTACK(S(exquo)); /* ARITHMETIC-ERROR slot OPERATION */
379 pushSTACK(STACK_(1+1)); pushSTACK(STACK_(0+2));
380 { var object tmp = listof(2); pushSTACK(tmp); } /* ARITHMETIC-ERROR slot OPERANDS */
381 pushSTACK(STACK_(1+2)); /* x */
382 pushSTACK(STACK_(0+3)); /* y */
383 error(arithmetic_error,GETTEXT("quotient ~S / ~S is not an integer"));
384 }
385
386 # Dividiert zwei Integers x,y >=0 und liefert den Quotienten x/y >=0.
387 # Bei y=0 Error. Die Division muss aufgehen, sonst Error.
388 # I_I_exquopos_I(x,y)
389 # > x,y: Integers >=0
390 # < result: Quotient x/y, ein Integer >=0
391 # can trigger GC
392 local object I_I_exquopos_I (object x, object y);
393 # Methode:
394 # (exquopos x y) :==
395 # (DIVIDE x y) -> q,r
396 # Falls r<>0, Error.
397 # Liefere q.
I_I_exquopos_I(object x,object y)398 local maygc object I_I_exquopos_I (object x, object y)
399 {
400 pushSTACK(y);
401 pushSTACK(x);
402 # Stackaufbau: y, x.
403 I_I_divide_I_I(x,y); # q,r auf den Stack
404 # Stackaufbau: y, x, q, r.
405 if (!eq(STACK_0,Fixnum_0)) {
406 skipSTACK(2); error_exquo();
407 }
408 var object q = STACK_1;
409 skipSTACK(4); return q;
410 }
411
412 # Dividiert zwei Integers x,y und liefert den Quotienten x/y.
413 # Bei y=0 Error. Die Division muss aufgehen, sonst Error.
414 # I_I_exquo_I(x,y)
415 # > x,y: Integers
416 # < result: Quotient x/y, ein Integer
417 # can trigger GC
418 local object I_I_exquo_I (object x, object y);
419 # Methode:
420 # (exquo x y) :==
421 # (DIVIDE (abs x) (abs y)) -> q,r
422 # Falls r<>0, Error.
423 # Falls x,y verschiedene Vorzeichen haben, liefere -q, sonst q.
I_I_exquo_I(object x,object y)424 local maygc object I_I_exquo_I (object x, object y)
425 {
426 pushSTACK(y);
427 pushSTACK(x);
428 pushSTACK(I_abs_I(y));
429 # Stackaufbau: y, x, (abs y).
430 x = I_abs_I(STACK_1); # (abs x)
431 I_I_divide_I_I(x,STACK_0); # q,r auf den Stack
432 # Stackaufbau: y, x, (abs y), q, r.
433 if (!eq(STACK_0,Fixnum_0)) {
434 skipSTACK(3); error_exquo();
435 }
436 var object q = STACK_1;
437 if (!same_sign_p(STACK_3,STACK_4)) {
438 # x,y haben verschiedene Vorzeichen
439 skipSTACK(5); return I_minus_I(q);
440 } else {
441 # x,y haben gleiche Vorzeichen
442 skipSTACK(5); return q;
443 }
444 }
445
446 # I_I_mod_I(x,y) = (mod x y), wo x,y Integers sind.
447 # can trigger GC
448 local object I_I_mod_I (object x, object y);
449 # Methode:
450 # (mod x y) :==
451 # (DIVIDE (abs x) (abs y)) -> q,r
452 # Falls r=0, liefere 0.
453 # Falls x,y verschiedene Vorzeichen haben, setze r:=r-abs(y).
454 # Falls x<0, setze r:=-r.
455 # Liefere r.
I_I_mod_I(object x,object y)456 local maygc object I_I_mod_I (object x, object y)
457 {
458 pushSTACK(y);
459 pushSTACK(x);
460 pushSTACK(I_abs_I(y));
461 # Stackaufbau: y, x, (abs y).
462 x = I_abs_I(STACK_1); # (abs x)
463 I_I_divide_I_I(x,STACK_0); # q,r auf den Stack
464 # Stackaufbau: y, x, (abs y), q, r.
465 var object r = STACK_0;
466 if (!eq(r,Fixnum_0)) {
467 if (!same_sign_p(STACK_3,STACK_4)) {
468 # x,y haben verschiedene Vorzeichen
469 r = I_I_minus_I(r,STACK_2); # r := (- r (abs y))
470 }
471 if (R_minusp(STACK_3))
472 r = I_minus_I(r); # x<0 -> r := (- r)
473 }
474 skipSTACK(5); return r;
475 }
476
477 # I_I_rem_I(x,y) = (rem x y), wo x,y Integers sind.
478 # can trigger GC
479 local object I_I_rem_I (object x, object y);
480 # Methode:
481 # (rem x y) :==
482 # (DIVIDE (abs x) (abs y)) -> q,r
483 # Falls x<0, setze r:=-r.
484 # Liefere r.
I_I_rem_I(object x,object y)485 local maygc object I_I_rem_I (object x, object y)
486 {
487 pushSTACK(y);
488 pushSTACK(x);
489 pushSTACK(I_abs_I(y));
490 # Stackaufbau: y, x, (abs y).
491 x = I_abs_I(STACK_1); # (abs x)
492 I_I_divide_I_I(x,STACK_0); # q,r auf den Stack
493 # Stackaufbau: y, x, (abs y), q, r.
494 var object r = STACK_0;
495 if (!eq(r,Fixnum_0)) {
496 if (R_minusp(STACK_3))
497 r = I_minus_I(r); # x<0 -> r := (- r)
498 }
499 skipSTACK(5); return r;
500 }
501
502 # Dividiert zwei Integers x,y und liefert Quotient und Rest
503 # (q,r) := (floor x y)
504 # I_I_floor_I_I(x,y);
505 # > x,y: Integers
506 # < STACK_1: Quotient q
507 # < STACK_0: Rest r
508 # Erniedrigt STACK um 2
509 # can trigger GC
510 local void I_I_floor_I_I (object x, object y);
511 # Methode:
512 # (floor x y) :==
513 # (DIVIDE (abs x) (abs y)) -> q,r
514 # Falls x,y verschiedene Vorzeichen haben und r<>0,
515 # setze q:=q+1 und r:=r-abs(y).
516 # Falls x<0, setze r:=-r.
517 # Falls x,y verschiedene Vorzeichen haben, setze q:=-q.
518 # Liefere q,r.
I_I_floor_I_I(object x,object y)519 local maygc void I_I_floor_I_I (object x, object y)
520 {
521 pushSTACK(y);
522 pushSTACK(x);
523 pushSTACK(I_abs_I(y));
524 # Stackaufbau: y, x, (abs y).
525 x = I_abs_I(STACK_1); # (abs x)
526 I_I_divide_I_I(x,STACK_0); # q,r auf den Stack
527 # Stackaufbau: y, x, (abs y), q, r.
528 if (!same_sign_p(STACK_3,STACK_4))
529 # x,y haben verschiedene Vorzeichen
530 if (!eq(STACK_0,Fixnum_0)) {
531 # r/=0, also r>0.
532 STACK_1 = I_1_plus_I(STACK_1); # q := (1+ q)
533 STACK_0 = I_I_minus_I(STACK_0,STACK_2); # r := (- r (abs y))
534 }
535 if (R_minusp(STACK_3)) {
536 # x<0
537 STACK_0 = I_minus_I(STACK_0); # r := (- r)
538 if (!R_minusp(STACK_4)) # y>=0 ?
539 goto negate_q; # q := (- q)
540 } else {
541 # x>=0
542 if (R_minusp(STACK_4)) # y<0 ?
543 negate_q: { STACK_1 = I_minus_I(STACK_1); } # q := (- q)
544 }
545 STACK_4 = STACK_1; STACK_3 = STACK_0; skipSTACK(3); # Stack aufräumen
546 }
547
548 # Dividiert zwei Integers x,y und liefert Quotient und Rest
549 # (q,r) := (ceiling x y)
550 # I_I_ceiling_I_I(x,y);
551 # > x,y: Integers
552 # < STACK_1: Quotient q
553 # < STACK_0: Rest r
554 # Erniedrigt STACK um 2
555 # can trigger GC
556 local void I_I_ceiling_I_I (object x, object y);
557 # Methode:
558 # (ceiling x y) :==
559 # (DIVIDE (abs x) (abs y)) -> q,r
560 # Falls x,y selbes Vorzeichen haben und r<>0,
561 # setze q:=q+1 und r:=r-abs(y).
562 # Falls x<0, setze r:=-r.
563 # Falls x,y verschiedene Vorzeichen haben, setze q:=-q.
564 # Liefere q,r.
I_I_ceiling_I_I(object x,object y)565 local maygc void I_I_ceiling_I_I (object x, object y)
566 {
567 pushSTACK(y);
568 pushSTACK(x);
569 pushSTACK(I_abs_I(y));
570 # Stackaufbau: y, x, (abs y).
571 x = I_abs_I(STACK_1); # (abs x)
572 I_I_divide_I_I(x,STACK_0); # q,r auf den Stack
573 # Stackaufbau: y, x, (abs y), q, r.
574 if (same_sign_p(STACK_3,STACK_4))
575 # x,y haben selbes Vorzeichen
576 if (!eq(STACK_0,Fixnum_0)) {
577 # r/=0, also r>0.
578 STACK_1 = I_1_plus_I(STACK_1); # q := (1+ q)
579 STACK_0 = I_I_minus_I(STACK_0,STACK_2); # r := (- r (abs y))
580 }
581 if (R_minusp(STACK_3)) {
582 # x<0
583 STACK_0 = I_minus_I(STACK_0); # r := (- r)
584 if (!R_minusp(STACK_4)) # y>=0 ?
585 goto negate_q; # q := (- q)
586 } else {
587 # x>=0
588 if (R_minusp(STACK_4)) # y<0 ?
589 negate_q: { STACK_1 = I_minus_I(STACK_1); } # q := (- q)
590 }
591 STACK_4 = STACK_1; STACK_3 = STACK_0; skipSTACK(3); # Stack aufräumen
592 }
593
594 # Dividiert zwei Integers x,y und liefert Quotient und Rest
595 # (q,r) := (truncate x y)
596 # I_I_truncate_I_I(x,y);
597 # > x,y: Integers
598 # < STACK_1: Quotient q
599 # < STACK_0: Rest r
600 # Erniedrigt STACK um 2
601 # can trigger GC
602 local void I_I_truncate_I_I (object x, object y);
603 # Methode:
604 # (truncate x y) :==
605 # (DIVIDE (abs x) (abs y)) -> q,r
606 # Falls x<0, setze r:=-r.
607 # Falls x,y verschiedene Vorzeichen haben, setze q:=-q.
608 # Liefere q,r.
I_I_truncate_I_I(object x,object y)609 local maygc void I_I_truncate_I_I (object x, object y)
610 {
611 pushSTACK(y);
612 pushSTACK(x);
613 pushSTACK(I_abs_I(y));
614 # Stackaufbau: y, x, (abs y).
615 x = I_abs_I(STACK_1); # (abs x)
616 I_I_divide_I_I(x,popSTACK()); # q,r auf den Stack
617 # Stackaufbau: y, x, q, r.
618 if (R_minusp(STACK_2)) {
619 # x<0
620 STACK_0 = I_minus_I(STACK_0); # r := (- r)
621 if (!R_minusp(STACK_3)) # y>=0 ?
622 goto negate_q; # q := (- q)
623 } else {
624 # x>=0
625 if (R_minusp(STACK_3)) # y<0 ?
626 negate_q: { STACK_1 = I_minus_I(STACK_1); } # q := (- q)
627 }
628 STACK_3 = STACK_1; STACK_2 = STACK_0; skipSTACK(2); # Stack aufräumen
629 }
630
631 # Dividiert zwei Integers x,y und liefert Quotient und Rest
632 # (q,r) := (round x y)
633 # I_I_round_I_I(x,y);
634 # > x,y: Integers
635 # < STACK_1: Quotient q
636 # < STACK_0: Rest r
637 # Erniedrigt STACK um 2
638 # can trigger GC
639 local void I_I_round_I_I (object x, object y);
640 # Methode:
641 # (round x y) :==
642 # (DIVIDE (abs x) (abs y)) -> q,r
643 # Setze s:=abs(y)-r.
644 # Falls (r>s) oder (r=s und q ungerade),
645 # (d.h. falls r>abs(y)/2 oder r=abs(y)/2 und q ungerade),
646 # setze q:=q+1 und r:=-s (d.h. r:=r-abs(y)).
647 # {Nun ist abs(r) <= abs(y)/2, bei abs(r)=abs(y)/2 ist q gerade.}
648 # Falls x<0, setze r:=-r.
649 # Falls x,y verschiedene Vorzeichen haben, setze q:=-q.
650 # Liefere q,r.
I_I_round_I_I(object x,object y)651 local maygc void I_I_round_I_I (object x, object y)
652 {
653 pushSTACK(y);
654 pushSTACK(x);
655 pushSTACK(I_abs_I(y));
656 # Stackaufbau: y, x, (abs y).
657 x = I_abs_I(STACK_1); # (abs x)
658 I_I_divide_I_I(x,STACK_0); # q,r auf den Stack
659 # Stackaufbau: y, x, (abs y), q, r.
660 var object s = I_I_minus_I(STACK_2,STACK_0); # (- (abs y) r)
661 var signean comp_r_s = I_I_comp(STACK_0,s); # vergleiche r und s
662 if ((comp_r_s>0) || ((comp_r_s==0) && (I_oddp(STACK_1)))) { # (r>s) oder (r=s und q ungerade) ?
663 STACK_0 = I_minus_I(s); # r := (- s) = (- r (abs y))
664 STACK_1 = I_1_plus_I(STACK_1); # q := (1+ q)
665 }
666 if (R_minusp(STACK_3)) {
667 # x<0
668 STACK_0 = I_minus_I(STACK_0); # r := (- r)
669 if (!R_minusp(STACK_4)) # y>=0 ?
670 goto negate_q; # q := (- q)
671 } else {
672 # x>=0
673 if (R_minusp(STACK_4)) # y<0 ?
674 negate_q: { STACK_1 = I_minus_I(STACK_1); } # q := (- q)
675 }
676 STACK_4 = STACK_1; STACK_3 = STACK_0; skipSTACK(3); # Stack aufräumen
677 }
678
679 /* (EXT:MOD-EXPT base exponent modulo) = (MOD (EXPT base exponent) modulus)
680 where exponent >= 0. */
I_I_I_mod_expt_I(object base,object exponent,object modulus)681 local maygc object I_I_I_mod_expt_I (object base, object exponent, object modulus)
682 {
683 /* See I_I_expt_I() in intmal.d for algorithm description.
684 A much better algorithm is [Cohen, Algorithm 1.2.1.], implemented in
685 cln-1.1.5/src/modinteger/cl_MI_std.h. */
686 pushSTACK(base); pushSTACK(exponent); pushSTACK(modulus);
687 if (eq(exponent,Fixnum_0)) {
688 pushSTACK(Fixnum_1); /* result := 1 */
689 } else {
690 while (!I_oddp(exponent)) {
691 STACK_2 = I_square_I(STACK_2); /* base := base * base */
692 STACK_2 = I_I_mod_I(STACK_2,STACK_0); /* base := base mod modulus */
693 STACK_1 = exponent = I_I_ash_I(STACK_1,Fixnum_minus1); /* exponent /= 2 */
694 }
695 pushSTACK(STACK_2); /* result := base */
696 /* STACK layout: base exponent modulus result */
697 while (!eq(exponent=STACK_2,Fixnum_1)) {
698 STACK_2 = I_I_ash_I(exponent,Fixnum_minus1); /* exponent /= 2 */
699 var object a = STACK_3 = I_square_I(STACK_3); /* base := base * base */
700 a = STACK_3 = I_I_mod_I(a,STACK_1); /* base := base mod modulus */
701 if (I_oddp(STACK_2))
702 /* we do not take MOD here because this will not grow too much anyway */
703 STACK_0 = I_I_mult_I(a,STACK_0); /* result *= base */
704 }
705 }
706 var object result = I_I_mod_I(STACK_0,STACK_1);
707 skipSTACK(4);
708 return result;
709 }
710