1 # Multiplikation ganzer Zahlen
2
3 # meldet Überlauf bei der Multiplikation:
mal_ueberlauf(void)4 local _Noreturn void mal_ueberlauf (void) {
5 pushSTACK(TheSubr(subr_self)->name); /* slot :OPERATION */
6 pushSTACK(NIL); /* slot :OPERANDS not available */
7 error(arithmetic_error,
8 GETTEXT("overflow during multiplication of large numbers"));
9 }
10
11 # karatsuba_threshold = Länge, ab der die Karatsuba-Multiplikation bevorzugt
12 # wird. Der Break-Even-Point bestimmt sich aus Zeitmessungen.
13 # Als Test dient (progn (time (! 5000)) nil), das viele kleine und einige
14 # ganz große Multiplikationen durchführt. Miss die Runtime.
15 # Unter Linux mit einem 80486: Auf einer Sparc 2:
16 # threshold time in 0.01 sec.
17 # 5 125
18 # 6 116
19 # 7 107
20 # 8 101
21 # 9 99
22 # 10 98
23 # 11 97
24 # 12 96
25 # 13 97
26 # 14 97
27 # 15 97
28 # 16 98
29 # 17 98
30 # 18 98
31 # 19 98
32 # 20 99
33 # 25 103
34 # 30 109
35 # 40 115
36 # 50 122
37 # 70 132
38 # 100 151
39 # 150 164
40 # 250 183
41 # 500 203
42 # 1000 203
43 # Das Optimum scheint bei karatsuba_threshold = 12 zu liegen.
44 # Da das Optimum aber vom Verhältnis
45 # Zeit für uintD-Multiplikation / Zeit für uintD-Addition
46 # abhängt und die gemessenen Zeiten auf eine Unterschreitung des Optimums
47 # empfindlicher reagieren als auf eine Überschreitung des Optimums,
48 # sind wir vorsichtig und wählen einen Wert etwas über dem Optimum:
49 #define karatsuba_threshold 16
50
51 # Quadriert eine Unsigned-Digit-sequence.
52 # UDS_square_UDS(len1,LSDptr1, MSDptr=,len=,LSDptr=);
53 # quadriert die UDS ../len1/LSDptr1.
54 # Dabei sollte len1>0 sein.
55 # Ergebnis ist die UDS MSDptr/len/LSDptr, mit len=2*len1, im Stack.
56 # Dabei wird num_stack erniedrigt.
57 #define UDS_square_UDS(len1,LSDptr1, MSDptr_assignment,len_assignment,LSDptr_assignment) \
58 { \
59 var uintL len_from_UDSmal = 2*(uintL)(len1); \
60 var uintD* LSDptr_from_UDSmal; \
61 if ((intWCsize < 32) && (len_from_UDSmal > (uintL)(bitc(intWCsize)-1))) { \
62 RESTORE_NUM_STACK; \
63 mal_ueberlauf(); \
64 } \
65 unused (len_assignment len_from_UDSmal); \
66 num_stack_need(len_from_UDSmal,MSDptr_assignment,LSDptr_assignment LSDptr_from_UDSmal =); \
67 square_2loop_down((LSDptr1),(len1),LSDptr_from_UDSmal); \
68 }
69
70 # Quadrier-Doppelschleife:
71 # Quadriert eine UDS und legt das Ergebnis in einer anderen UDS ab.
72 # square_2loop_down(sourceptr,len,destptr);
73 # quadriert die UDS sourceptr[-len..-1] (len1>0)
74 # und legt das Ergebnis in der UDS destptr[-2*len..-1] ab.
75 # Unterhalb von destptr werden 2*len Digits Platz benötigt.
76 local void square_2loop_down (const uintD* sourceptr, uintC len,
77 uintD* destptr);
78 local void square_2bigloop_down (const uintD* sourceptr, uintC len,
79 uintD* destptr);
square_2loop_down(const uintD * sourceptr,uintC len,uintD * destptr)80 local void square_2loop_down (const uintD* sourceptr, uintC len,
81 uintD* destptr)
82 {
83 if (len==1) {
84 var uintD digit = sourceptr[-1];
85 #if HAVE_DD
86 var uintDD prod = muluD(digit,digit);
87 destptr[-1] = lowD(prod); destptr[-2] = highD(prod);
88 #else
89 muluD(digit,digit, destptr[-2] =, destptr[-1] =);
90 #endif
91 } else if (len < karatsuba_threshold) {
92 # Multiplikation nach Schulmethode
93 # Gemischte Produkte:
94 # 2*( x[n-1..1] * x[0] * b^1
95 # + x[n-1..2] * x[1] * b^3
96 # + ...
97 # + x[n-1..n-1] * x[n-2] * b^(2*n-3))
98 {
99 var const uintD* sourceptr1 = sourceptr;
100 var uintD* destptr2 = destptr;
101 *--destptr2 = 0;
102 var uintC count = len-1;
103 {
104 var uintD digit = *--sourceptr1;
105 mulu_loop_down(digit,sourceptr1,destptr2,count);
106 }
107 var uintD* destptr1 = destptr - (len+1);
108 while (--count > 0) {
109 destptr2 -= 2;
110 var uintD digit = *--sourceptr1;
111 var uintD carry = muluadd_loop_down(digit,sourceptr1,destptr2,count);
112 *--destptr1 = carry;
113 }
114 var uintD carry = shift1left_loop_down(destptr-1,2*len-2);
115 destptr1[-1] = (carry==0 ? 0 : 1);
116 }
117 # Quadrate:
118 len = 2*len;
119 do {
120 len -= 2;
121 var uintD digit = *--sourceptr;
122 #if HAVE_DD
123 var uintDD prod = muluD(digit,digit);
124 var uintDD accu = highlowDD(destptr[-2],destptr[-1]);
125 accu += prod;
126 destptr[-1] = lowD(accu); destptr[-2] = highD(accu);
127 destptr -= 2;
128 if (accu < prod)
129 inc_loop_down(destptr,len);
130 #else
131 var uintD hi;
132 var uintD lo;
133 var uintD tmp;
134 muluD(digit,digit, hi=,lo=);
135 tmp = destptr[-1] + lo; destptr[-1] = tmp;
136 if (tmp < lo)
137 hi++;
138 tmp = destptr[-2] + hi; destptr[-2] = tmp;
139 destptr -= 2;
140 if (tmp < hi)
141 inc_loop_down(destptr,len);
142 #endif
143 } while (len > 0);
144 } else { # len groß
145 # Karatsuba-Quadrierung
146 # (ausgelagert, um die eigentliche Quadrierfunktion nicht
147 # durch zu viele Registervariablen zu belasten):
148 square_2bigloop_down(sourceptr,len,destptr);
149 }
150 }
square_2bigloop_down(const uintD * sourceptr,uintC len,uintD * destptr)151 local void square_2bigloop_down (const uintD* sourceptr, uintC len,
152 uintD* destptr)
153 # Karatsuba-Quadrierung
154 {
155 # Es ist 2 <= len.
156 SAVE_NUM_STACK
157 var uintC prod_len = 2*len;
158 var uintD* prod_LSDptr = destptr;
159 var uintC k_hi = floor(len,2); # Länge der High-Teile: floor(len/2) >0
160 var uintC k_lo = len - k_hi; # Länge der Low-Teile: ceiling(len/2) >0
161 # Es gilt k_hi <= k_lo <= len, k_lo + k_hi = len.
162 # Summe x1+x0 berechnen:
163 var uintD* sum_MSDptr;
164 var uintC sum_len = k_lo; # = max(k_lo,k_hi)
165 var uintD* sum_LSDptr;
166 num_stack_need(sum_len+1,sum_MSDptr=,sum_LSDptr=);
167 sum_MSDptr++; # 1 Digit vorne als Reserve
168 {
169 var uintD carry = # Hauptteile von x1 und x0 addieren:
170 add_loop_down(sourceptr-k_lo,sourceptr,sum_LSDptr,k_hi);
171 if (!(k_lo==k_hi)) {
172 # noch k_lo-k_hi = 1 Digits abzulegen
173 sum_MSDptr[0] = sourceptr[-(uintP)k_lo]; # = sourceptr[-(uintP)k_hi-1]
174 if (!(carry==0)) {
175 if (++(sum_MSDptr[0]) == 0)
176 carry=1;
177 else
178 carry=0;
179 }
180 }
181 if (carry) {
182 *--sum_MSDptr = 1;
183 sum_len++;
184 }
185 }
186 # Platz für Produkte x0*x0, x1*x1:
187 {
188 var uintC prodhi_len = 2*k_hi;
189 var uintD* prodhi_LSDptr = prod_LSDptr - 2*k_lo;
190 # prod_MSDptr/2*len/prod_LSDptr wird zuerst die beiden
191 # Produkte x1*x1 in prod_MSDptr/2*k_hi/prodhi_LSDptr
192 # und x0*x0 in prodhi_LSDptr/2*k_lo/prod_LSDptr,
193 # dann das Produkt (b^k*x1+x0)*(b^k*x1+x0) enthalten.
194 # Platz fürs Produkt (x1+x0)*(x1+x0) belegen:
195 var uintD* prodmid_MSDptr;
196 var uintC prodmid_len = 2*sum_len;
197 var uintD* prodmid_LSDptr;
198 num_stack_need(prodmid_len,prodmid_MSDptr=,prodmid_LSDptr=);
199 # Produkt (x1+x0)*(x1+x0) berechnen:
200 square_2loop_down(sum_LSDptr,sum_len,prodmid_LSDptr);
201 # Das Produkt beansprucht 2*k_lo + (0 oder 1) <= 2*sum_len = prodmid_len Digits.
202 # Produkt x0*x0 berechnen:
203 square_2loop_down(sourceptr,k_lo,prod_LSDptr);
204 # Produkt x1*x1 berechnen:
205 square_2loop_down(sourceptr-k_lo,k_hi,prodhi_LSDptr);
206 # Und x1*x1 abziehen:
207 {
208 var uintD carry =
209 subfrom_loop_down(prodhi_LSDptr,prodmid_LSDptr,prodhi_len);
210 # Carry um maximal prodmid_len-prodhi_len Digits weitertragen:
211 if (!(carry==0))
212 dec_loop_down(prodmid_LSDptr-prodhi_len,prodmid_len-prodhi_len);
213 }
214 # Und x0*x0 abziehen:
215 {
216 var uintD carry =
217 subfrom_loop_down(prod_LSDptr,prodmid_LSDptr,2*k_lo);
218 # Falls Carry: Produkt beansprucht 2*k_lo+1 Digits.
219 # Carry um maximal 1 Digit weitertragen:
220 if (!(carry==0))
221 prodmid_LSDptr[-2*(uintP)k_lo-1] -= 1;
222 }
223 # prodmid_LSDptr[-prodmid_len..-1] enthält nun 2*x0*x1.
224 # Dies ist < 2 * b^k_lo * b^k_hi = 2 * b^len,
225 # passt also in len+1 Digits.
226 # prodmid_len, wenn möglich, um maximal 2 verkleinern:
227 # (benutzt prodmid_len >= 2*k_lo >= len >= 2)
228 if (prodmid_MSDptr[0]==0) {
229 prodmid_len--;
230 if (prodmid_MSDptr[1]==0)
231 prodmid_len--;
232 }
233 # Nun ist k_lo+prodmid_len <= 2*len .
234 # (Denn es war prodmid_len = 2*sum_len <= 2*(k_lo+1)
235 # <= len+3, und nach 2-maliger Verkleinerung jedenfalls
236 # prodmid_len <= len+1. Wegen k_lo < len also
237 # k_lo + prodmid_len <= (len-1)+(len+1) = 2*len.)
238 # prodmid*b^k = 2*x0*x1*b^k zu prod = x1*x1*b^(2*k) + x0*x0 addieren:
239 {
240 var uintD carry =
241 addto_loop_down(prodmid_LSDptr,prod_LSDptr-k_lo,prodmid_len);
242 if (!(carry==0))
243 inc_loop_down(prod_LSDptr-(k_lo+prodmid_len),prod_len-(k_lo+prodmid_len));
244 }
245 }
246 RESTORE_NUM_STACK
247 }
248
249 # (* x x), wo x ein Integer ist. Ergebnis Integer.
250 # can trigger GC
251 # Methode:
252 # x=0 -> Ergebnis 0
253 # x Fixnum -> direkt quadrieren
254 # sonst: zu WS machen, quadrieren.
I_square_I(object x)255 local maygc object I_square_I (object x)
256 {
257 if (eq(x,Fixnum_0))
258 return Fixnum_0;
259 if (I_fixnump(x)) {
260 var sintV x_ = FN_to_V(x);
261 #if (oint_data_len+1 > 32)
262 # nur falls x ein Integer mit höchstens 32 Bit ist:
263 if ((uintV)((sintV)FN_sign(x) ^ x_) < vbit(31))
264 #endif
265 {
266 # Wert direkt quadrieren:
267 var uint32 hi;
268 var uint32 lo;
269 mulu32((uint32)x_,(uint32)x_,hi=,lo=); # erst unsigned multiplizieren
270 if (x_ < 0)
271 hi -= 2*(uint32)x_; # dann Korrektur für Vorzeichen
272 return L2_to_I(hi,lo);
273 }
274 }
275 var uintD* xMSDptr;
276 var uintC xlen;
277 var uintD* xLSDptr;
278 var uintD* ergMSDptr;
279 var uintC erglen;
280 var uintD* ergLSDptr;
281 I_to_NDS_nocopy(x, xMSDptr = , xlen = , xLSDptr = );
282 erglen = 2*xlen;
283 /* we use intWsize(16) instead of intWCsize(32) to avoid
284 __builtin_alloca running out of stack and causing a segfault */
285 if ((intWsize < 32) && (erglen > (uintL)(bitc(intWsize)-1)))
286 mal_ueberlauf();
287 {
288 SAVE_NUM_STACK # num_stack retten
289 num_stack_need(erglen, ergMSDptr = , ergLSDptr = );
290 begin_arith_call();
291 {
292 var uintC len = xlen;
293 var uintD MSD = xMSDptr[0];
294 if (MSD == 0) {
295 ergMSDptr[0] = 0; ergMSDptr[1] = 0; len--;
296 }
297 square_2loop_down(xLSDptr,len,ergLSDptr);
298 if ((sintD)MSD < 0) {
299 subfrom_loop_down(xLSDptr,ergLSDptr-xlen,xlen);
300 subfrom_loop_down(xLSDptr,ergLSDptr-xlen,xlen);
301 }
302 }
303 end_arith_call();
304 var object result = DS_to_I(ergMSDptr,erglen);
305 RESTORE_NUM_STACK # num_stack zurück
306 return result;
307 }
308 }
309
310 # Multipliziert zwei Unsigned-Digit-sequences.
311 # UDS_UDS_mult_UDS(len1,LSDptr1, len2,LSDptr2, MSDptr=,len=,LSDptr=);
312 # multipliziert die UDS ../len1/LSDptr1 und ../len2/LSDptr2.
313 # Dabei sollte len1>0 und len2>0 sein.
314 # Ergebnis ist die UDS MSDptr/len/LSDptr, mit len=len1+len2, im Stack.
315 # Dabei wird num_stack erniedrigt.
316 #define UDS_UDS_mult_UDS(len1,LSDptr1,len2,LSDptr2, MSDptr_assignment,len_assignment,LSDptr_assignment) \
317 { \
318 var uintL len_from_UDSmal = (uintL)(len1) + (uintL)(len2); \
319 var uintD* LSDptr_from_UDSmal; \
320 if ((intWCsize < 32) && (len_from_UDSmal > (uintL)(bitc(intWCsize)-1))) { \
321 RESTORE_NUM_STACK; \
322 mal_ueberlauf(); \
323 } \
324 unused (len_assignment len_from_UDSmal); \
325 num_stack_need(len_from_UDSmal,MSDptr_assignment,LSDptr_assignment LSDptr_from_UDSmal =); \
326 mulu_2loop_down((LSDptr1),(len1),(LSDptr2),(len2),LSDptr_from_UDSmal); \
327 }
328
329 # Multiplikations-Doppelschleife:
330 # Multipliziert zwei UDS und legt das Ergebnis in einer dritten UDS ab.
331 # mulu_2loop_down(sourceptr1,len1,sourceptr2,len2,destptr);
332 # multipliziert die UDS sourceptr1[-len1..-1] (len1>0)
333 # mit der UDS sourceptr2[-len1..-1] (len2>0)
334 # und legt das Ergebnis in der UDS destptr[-len..-1] (len=len1+len2) ab.
335 # Unterhalb von destptr werden len Digits Platz benötigt.
336 local void mulu_2loop_down (const uintD* sourceptr1, uintC len1,
337 const uintD* sourceptr2, uintC len2,
338 uintD* destptr);
339 local void mulu_2bigloop_down (const uintD* sourceptr1, uintC len1,
340 const uintD* sourceptr2, uintC len2,
341 uintD* destptr);
mulu_2loop_down(const uintD * sourceptr1,uintC len1,const uintD * sourceptr2,uintC len2,uintD * destptr)342 local void mulu_2loop_down (const uintD* sourceptr1, uintC len1,
343 const uintD* sourceptr2, uintC len2,
344 uintD* destptr)
345 {
346 # len1<=len2 erzwingen:
347 if (len1>len2) {
348 {
349 var const uintD* temp;
350 temp = sourceptr1; sourceptr1 = sourceptr2; sourceptr2 = temp;
351 }
352 {
353 var uintC temp;
354 temp = len1; len1 = len2; len2 = temp;
355 }
356 }
357 if (len1==1) {
358 # nur eine Einfachschleife
359 mulu_loop_down(*--sourceptr1,sourceptr2,destptr,len2);
360 } else if (len1 < karatsuba_threshold) {
361 # Multiplikation nach Schulmethode
362 # len2 Digits auf 0 setzen:
363 var uintD* destptr2 = clear_loop_down(destptr,len2);
364 # äußere Schleife läuft über source1 :
365 dotimespC(len1,len1, {
366 # innere Schleife läuft über source2 :
367 var uintD carry =
368 muluadd_loop_down(*--sourceptr1,sourceptr2,destptr,len2);
369 *--destptr2 = carry; # UDS um das Carry-Digit verlängern
370 destptr--;
371 });
372 } else { # len1 groß
373 # Karatsuba-Multiplikation
374 # (ausgelagert, um die eigentliche Multiplikationsfunktion nicht
375 # durch zu viele Registervariablen zu belasten):
376 mulu_2bigloop_down(sourceptr1,len1,sourceptr2,len2,destptr);
377 }
378 }
mulu_2bigloop_down(const uintD * sourceptr1,uintC len1,const uintD * sourceptr2,uintC len2,uintD * destptr)379 local void mulu_2bigloop_down (const uintD* sourceptr1, uintC len1,
380 const uintD* sourceptr2, uintC len2,
381 uintD* destptr)
382 {
383 # Karatsuba-Multiplikation
384 # Prinzip: (x1*b^k+x0) * (y1*b^k+y0)
385 # = x1*y1 * b^2k + ((x1+x0)*(y1+y0)-x1*y1-x0*y0) * b^k + x0*y0
386 # Methode 1 (Collins/Loos, Degel):
387 # source2 wird in floor(len2/len1) einzelne UDS mit je einer
388 # Länge len3 (len1 <= len3 < 2*len1) unterteilt,
389 # jeweils k=floor(len3/2).
390 # Methode 2 (Haible):
391 # source2 wird in ceiling(len2/len1) einzelne UDS mit je einer
392 # Länge len3 (0 < len3 <= len1) unterteilt, jeweils k=floor(len1/2).
393 # Aufwand für die hinteren Einzelteile:
394 # bei beiden Methoden jeweils 3*len1^2.
395 # Aufwand für das vorderste Teil (alles, falls len1 <= len2 < 2*len1)
396 # mit r = len1, s = (len2 mod len1) + len1 (>= len1, < 2*len1):
397 # bei Methode 1:
398 # | : | r
399 # | : | s
400 # (r-s/2)*s/2 + s/2*s/2 + s/2*s/2 = r*s/2 + s^2/4 .
401 # bei Methode 2:
402 # | : | r
403 # | | : | s
404 # (s-r)*r + r/2*r/2 + r/2*r/2 + r/2*r/2 = r*s - r^2/4 .
405 # Wegen (r*s/2 + s^2/4) - (r*s - r^2/4) = (r-s)^2/4 >= 0
406 # ist Methode 2 günstiger.
407 # Denkfehler! Dies gilt - wenn überhaupt - nur knapp oberhalb des
408 # Break-Even-Points.
409 # Im allgemeinen ist der Multiplikationsaufwand für zwei Zahlen der
410 # Längen u bzw. v nämlich gegeben durch min(u,v)^c * max(u,v),
411 # wobei c = log3/log2 - 1 = 0.585...
412 # Dadurch wird der Aufwand in Abhängigkeit des Parameters t = k,
413 # r/2 <= t <= s/2 (der einzig sinnvolle Bereich), zu
414 # (r-t)^c*(s-t) + t^c*(s-t) + t^(1+c).
415 # Dessen Optimum liegt (im Bereich r <= s <= 2*r)
416 # - im klassischen Fall c=1 tatsächlich stets bei t=r/2 [Methode 2],
417 # - im Karatsuba-Fall c=0.6 aber offenbar bei t=s/2 [Methode 1]
418 # oder ganz knapp darunter.
419 # Auch erweist sich Methode 1 im Experiment als effizienter.
420 # Daher implementieren wir Methode 1 :
421 #
422 # Es ist 2 <= len2 <= len1.
423 var bool first_part = true; # Flag, ob jetzt das erste Teilprodukt berechnet wird
424 if (len2 >= 2*len1) {
425 SAVE_NUM_STACK
426 # Teilprodukte von jeweils len1 mal len1 Digits bilden:
427 var uintC k_lo = floor(len1,2); # Länge der Low-Teile: floor(len1/2) >0
428 var uintC k_hi = len1 - k_lo; # Länge der High-Teile: ceiling(len1/2) >0
429 # Es gilt k_lo <= k_hi <= len1, k_lo + k_hi = len1.
430 # Summe x1+x0 berechnen:
431 var uintD* sum1_MSDptr;
432 var uintC sum1_len = k_hi; # = max(k_lo,k_hi)
433 var uintD* sum1_LSDptr;
434 num_stack_need(sum1_len+1,sum1_MSDptr=,sum1_LSDptr=);
435 sum1_MSDptr++; # 1 Digit vorne als Reserve
436 {
437 var uintD carry = # Hauptteile von x1 und x0 addieren:
438 add_loop_down(&sourceptr1[-(uintP)k_lo],sourceptr1,sum1_LSDptr,k_lo);
439 if (!(k_lo==k_hi)) {
440 # noch k_hi-k_lo = 1 Digits abzulegen
441 sum1_MSDptr[0] = sourceptr1[-(uintP)len1]; # = sourceptr1[-2*(uintP)k_lo-1]
442 if (!(carry==0)) {
443 if (++(sum1_MSDptr[0]) == 0)
444 carry=1;
445 else
446 carry=0;
447 }
448 }
449 if (carry) {
450 *--sum1_MSDptr = 1;
451 sum1_len++;
452 }
453 }
454 # Platz für Summe y1+y0 belegen:
455 var uintC sum2_maxlen = k_hi+1;
456 var uintD* sum2_LSDptr;
457 num_stack_need(sum2_maxlen,_EMA_,sum2_LSDptr=);
458 # Platz für Produkte x0*y0, x1*y1 belegen:
459 var uintD* prod_MSDptr;
460 var uintD* prod_LSDptr;
461 var uintD* prodhi_LSDptr;
462 num_stack_need(2*(uintL)len1,prod_MSDptr=,prod_LSDptr=);
463 prodhi_LSDptr = &prod_LSDptr[-2*(uintP)k_lo];
464 # prod_MSDptr/2*len1/prod_LSDptr wird zuerst die beiden
465 # Produkte x1*y1 in prod_MSDptr/2*k_hi/prodhi_LSDptr
466 # und x0*y0 in prodhi_LSDptr/2*k_lo/prod_LSDptr,
467 # dann das Produkt (b^k*x1+x0)*(b^k*y1+y0) enthalten.
468 # Platz fürs Produkt (x1+x0)*(y1+y0) belegen:
469 var uintD* prodmid_MSDptr;
470 var uintD* prodmid_LSDptr;
471 num_stack_need(sum1_len+sum2_maxlen,prodmid_MSDptr=,prodmid_LSDptr=);
472 # Schleife über die hinteren Einzelteile:
473 do {
474 # Produkt x0*y0 berechnen:
475 mulu_2loop_down(sourceptr1,k_lo,sourceptr2,k_lo,prod_LSDptr);
476 # Produkt x1*y1 berechnen:
477 mulu_2loop_down(&sourceptr1[-(uintP)k_lo],k_hi,&sourceptr2[-(uintP)k_lo],k_hi,prodhi_LSDptr);
478 # Summe y1+y0 berechnen:
479 {
480 var uintC sum2_len = k_hi; # = max(k_lo,k_hi)
481 var uintD* sum2_MSDptr = &sum2_LSDptr[-(uintP)sum2_len];
482 {
483 var uintD carry = # Hauptteile von y1 und y0 addieren:
484 add_loop_down(&sourceptr2[-(uintP)k_lo],sourceptr2,sum2_LSDptr,k_lo);
485 if (!(k_lo==k_hi)) {
486 # noch k_hi-k_lo = 1 Digits abzulegen
487 sum2_MSDptr[0] = sourceptr2[-(uintP)len1]; # = sourceptr2[-2*(uintP)k_lo-1]
488 if (!(carry==0)) {
489 if (++(sum2_MSDptr[0]) == 0)
490 carry=1;
491 else
492 carry=0;
493 }
494 }
495 if (carry) {
496 *--sum2_MSDptr = 1;
497 sum2_len++;
498 }
499 }
500 # Produkt (x1+x0)*(y1+y0) berechnen:
501 mulu_2loop_down(sum1_LSDptr,sum1_len,sum2_LSDptr,sum2_len,prodmid_LSDptr);
502 # Das Produkt beansprucht 2*k_hi + (0 oder 1) <= sum1_len + sum2_len Digits.
503 var uintC prodmid_len = sum1_len+sum2_len;
504 # Davon x1*y1 abziehen:
505 {
506 var uintD carry =
507 subfrom_loop_down(prodhi_LSDptr,prodmid_LSDptr,2*k_hi);
508 # Falls Carry: Produkt beansprucht 2*k_hi+1 Digits.
509 # Carry um maximal 1 Digit weitertragen:
510 if (!(carry==0))
511 prodmid_LSDptr[-2*(uintP)k_hi-1] -= 1;
512 }
513 # Und x0*y0 abziehen:
514 {
515 var uintD carry =
516 subfrom_loop_down(prod_LSDptr,prodmid_LSDptr,2*k_lo);
517 # Carry um maximal prodmid_len-2*k_lo Digits weitertragen:
518 if (!(carry==0))
519 dec_loop_down(&prodmid_LSDptr[-2*(uintP)k_lo],prodmid_len-2*k_lo);
520 }
521 # prodmid_LSDptr[-prodmid_len..-1] enthält nun x0*y1+x1*y0.
522 # Dies wird zu prod = x1*y1*b^(2*k) + x0*y0 addiert:
523 {
524 var uintD carry =
525 addto_loop_down(prodmid_LSDptr,&prod_LSDptr[-(uintP)k_lo],prodmid_len);
526 # (Benutze dabei k_lo+prodmid_len <= k_lo+2*(k_hi+1) = 2*len1-k_lo+2 <= 2*len1 .)
527 if (!(carry==0))
528 inc_loop_down(&prod_LSDptr[-(uintP)(k_lo+prodmid_len)],2*len1-(k_lo+prodmid_len));
529 }
530 }
531 # Das Teilprodukt zum Gesamtprodukt addieren:
532 if (first_part) {
533 copy_loop_down(prod_LSDptr,destptr,2*len1);
534 destptr -= len1;
535 first_part = false;
536 } else {
537 var uintD carry =
538 addto_loop_down(prod_LSDptr,destptr,len1);
539 destptr -= len1;
540 copy_loop_down(&prod_LSDptr[-(uintP)len1],destptr,len1);
541 if (!(carry==0))
542 inc_loop_down(destptr,len1);
543 }
544 sourceptr2 -= len1; len2 -= len1;
545 } while (len2 >= 2*len1);
546 RESTORE_NUM_STACK
547 }
548 # Nun ist len1 <= len2 < 2*len1.
549 # letztes Teilprodukt von len1 mal len2 Digits bilden:
550 {
551 SAVE_NUM_STACK
552 var uintD* prod_MSDptr;
553 var uintC prod_len = len1+len2;
554 var uintD* prod_LSDptr;
555 num_stack_need((uintL)prod_len,prod_MSDptr=,prod_LSDptr=);
556 {
557 var uintC k_hi = floor(len2,2); # Länge der High-Teile: floor(len2/2) >0
558 var uintC k_lo = len2 - k_hi; # Länge der Low-Teile: ceiling(len2/2) >0
559 # Es gilt k_hi <= k_lo <= len1 <= len2, k_lo + k_hi = len2.
560 var uintC x1_len = len1-k_lo; # <= len2-k_lo = k_hi <= k_lo
561 # Summe x1+x0 berechnen:
562 var uintD* sum1_MSDptr;
563 var uintC sum1_len = k_lo; # = max(k_lo,k_hi)
564 var uintD* sum1_LSDptr;
565 num_stack_need(sum1_len+1,sum1_MSDptr=,sum1_LSDptr=);
566 sum1_MSDptr++; # 1 Digit vorne als Reserve
567 {
568 var uintD carry = # x1 und unteren Teil von x0 addieren:
569 add_loop_down(&sourceptr1[-(uintP)k_lo],sourceptr1,sum1_LSDptr,x1_len);
570 # und den oberen Teil von x0 dazu:
571 copy_loop_down(&sourceptr1[-(uintP)x1_len],&sum1_LSDptr[-(uintP)x1_len],k_lo-x1_len);
572 if (!(carry==0)) {
573 carry = inc_loop_down(&sum1_LSDptr[-(uintP)x1_len],k_lo-x1_len);
574 if (carry) {
575 *--sum1_MSDptr = 1;
576 sum1_len++;
577 }
578 }
579 }
580 # Summe y1+y0 berechnen:
581 var uintD* sum2_MSDptr;
582 var uintC sum2_len = k_lo; # = max(k_lo,k_hi)
583 var uintD* sum2_LSDptr;
584 num_stack_need(sum2_len+1,sum2_MSDptr=,sum2_LSDptr=);
585 sum2_MSDptr++; # 1 Digit vorne als Reserve
586 {
587 var uintD carry = # Hauptteile von y1 und y0 addieren:
588 add_loop_down(&sourceptr2[-(uintP)k_lo],sourceptr2,sum2_LSDptr,k_hi);
589 if (!(k_lo==k_hi)) {
590 # noch k_lo-k_hi = 1 Digits abzulegen
591 sum2_MSDptr[0] = sourceptr2[-(uintP)k_lo]; # = sourceptr2[-(uintP)k_hi-1]
592 if (!(carry==0)) {
593 if (++(sum2_MSDptr[0]) == 0)
594 carry=1;
595 else
596 carry=0;
597 }
598 }
599 if (carry) {
600 *--sum2_MSDptr = 1;
601 sum2_len++;
602 }
603 }
604 # Platz für Produkte x0*y0, x1*y1:
605 var uintC prodhi_len = x1_len+k_hi;
606 var uintD* prodhi_LSDptr = &prod_LSDptr[-2*(uintP)k_lo];
607 # prod_MSDptr/len1+len2/prod_LSDptr wird zuerst die beiden
608 # Produkte x1*y1 in prod_MSDptr/x1_len+k_hi/prodhi_LSDptr
609 # und x0*y0 in prodhi_LSDptr/2*k_lo/prod_LSDptr,
610 # dann das Produkt (b^k*x1+x0)*(b^k*y1+y0) enthalten.
611 # Platz fürs Produkt (x1+x0)*(y1+y0) belegen:
612 var uintD* prodmid_MSDptr;
613 var uintC prodmid_len = sum1_len+sum2_len;
614 var uintD* prodmid_LSDptr;
615 num_stack_need(prodmid_len,prodmid_MSDptr=,prodmid_LSDptr=);
616 # Produkt (x1+x0)*(y1+y0) berechnen:
617 mulu_2loop_down(sum1_LSDptr,sum1_len,sum2_LSDptr,sum2_len,prodmid_LSDptr);
618 # Das Produkt beansprucht 2*k_lo + (0 oder 1) <= sum1_len + sum2_len = prodmid_len Digits.
619 # Produkt x0*y0 berechnen:
620 mulu_2loop_down(sourceptr1,k_lo,sourceptr2,k_lo,prod_LSDptr);
621 # Produkt x1*y1 berechnen:
622 if (!(x1_len==0)) {
623 mulu_2loop_down(&sourceptr1[-(uintP)k_lo],x1_len,&sourceptr2[-(uintP)k_lo],k_hi,prodhi_LSDptr);
624 # Und x1*y1 abziehen:
625 {
626 var uintD carry =
627 subfrom_loop_down(prodhi_LSDptr,prodmid_LSDptr,prodhi_len);
628 # Carry um maximal prodmid_len-prodhi_len Digits weitertragen:
629 if (!(carry==0))
630 dec_loop_down(&prodmid_LSDptr[-(uintP)prodhi_len],prodmid_len-prodhi_len);
631 }
632 } else {
633 # Produkt x1*y1=0, nichts abzuziehen
634 clear_loop_down(prodhi_LSDptr,prodhi_len);
635 }
636 # Und x0*y0 abziehen:
637 {
638 var uintD carry =
639 subfrom_loop_down(prod_LSDptr,prodmid_LSDptr,2*k_lo);
640 # Falls Carry: Produkt beansprucht 2*k_lo+1 Digits.
641 # Carry um maximal 1 Digit weitertragen:
642 if (!(carry==0))
643 prodmid_LSDptr[-2*(uintP)k_lo-1] -= 1;
644 }
645 # prodmid_LSDptr[-prodmid_len..-1] enthält nun x0*y1+x1*y0.
646 # Dies ist < b^k_lo * b^k_hi + b^x1_len * b^k_lo
647 # = b^len2 + b^len1 <= 2 * b^len2,
648 # passt also in len2+1 Digits.
649 # Im Fall x1_len=0 ist es sogar < b^k_lo * b^k_hi = b^len2,
650 # es passt also in len2 Digits.
651 # prodmid_len, wenn möglich, um maximal 2 verkleinern:
652 # (benutzt prodmid_len >= 2*k_lo >= len2 >= 2)
653 if (prodmid_MSDptr[0]==0) {
654 prodmid_len--;
655 if (prodmid_MSDptr[1]==0)
656 prodmid_len--;
657 }
658 # Nun ist k_lo+prodmid_len <= len1+len2 .
659 # (Denn es war prodmid_len = sum1_len+sum2_len <= 2*(k_lo+1)
660 # <= len2+3, und nach 2-maliger Verkleinerung jedenfalls
661 # prodmid_len <= len2+1. Im Falle k_lo < len1 also
662 # k_lo + prodmid_len <= (len1-1)+(len2+1) = len1+len2.
663 # Im Falle k_lo = len1 aber ist x1_len=0, sum1_len = k_lo, also
664 # war prodmid_len = sum1_len+sum2_len <= 2*k_lo+1 <= len2+2,
665 # nach 2-maliger Verkleinerung jedenfalls prodmid_len <= len2.)
666 # prodmid*b^k = (x0*y1+x1*y0)*b^k zu prod = x1*y1*b^(2*k) + x0*y0 addieren:
667 {
668 var uintD carry =
669 addto_loop_down(prodmid_LSDptr,&prod_LSDptr[-(uintP)k_lo],prodmid_len);
670 if (!(carry==0))
671 inc_loop_down(&prod_LSDptr[-(uintP)(k_lo+prodmid_len)],prod_len-(k_lo+prodmid_len));
672 }
673 }
674 # Das Teilprodukt zum Gesamtprodukt addieren:
675 if (first_part) {
676 copy_loop_down(prod_LSDptr,destptr,prod_len);
677 } else {
678 var uintD carry =
679 addto_loop_down(prod_LSDptr,destptr,len1);
680 destptr -= len1;
681 copy_loop_down(&prod_LSDptr[-(uintP)len1],destptr,len2);
682 if (!(carry==0))
683 inc_loop_down(destptr,len2);
684 }
685 RESTORE_NUM_STACK
686 }
687 }
688
689 # Multipliziert zwei Digit-sequences.
690 # DS_DS_mult_DS(MSDptr1,len1,LSDptr1, MSDptr2,len2,LSDptr2, MSDptr=,len=,LSDptr=);
691 # multipliziert die DS MSDptr1/len1/LSDptr1 und MSDptr2/len2/LSDptr2.
692 # Dabei sollte len1>0 und len2>0 sein. Alles sollten Variablen sein!
693 # Ergebnis ist die DS MSDptr/len/LSDptr, mit len=len1+len2, im Stack.
694 # Dabei wird num_stack erniedrigt.
695 # Methode:
696 # Erst unsigned multiplizieren. Dann bis zu zwei Subtraktionen.
697 # Sei b=2^intDsize, k=len1, l=len2, n=DS1, m=DS2.
698 # Gesucht ist n * m.
699 # Wir errechnen erst das unsigned-product p (mod b^(k+l)).
700 # n>0, m>0: p = n*m, n*m = p
701 # n<0, m>0: p = (n+b^k)*m, n*m + b^(k+l) = p - b^k * m (mod b^(k+l)).
702 # n>0, m<0: p = n*(m+b^l), n*m + b^(k+l) = p - b^l * n (mod b^(k+l)).
703 # n<0, m<0: p = (n+b^k)*(m+b^l),
704 # n*m = p - b^k * (m+b^l) - b^l * (n+b^k) (mod b^(k+l)).
705 #define DS_DS_mult_DS(MSDptr1,len1,LSDptr1,MSDptr2,len2,LSDptr2, MSDptr_assignment,len_assignment,LSDptr_assignment) \
706 { \
707 var uintD* LSDptr0; \
708 UDS_UDS_mult_UDS(len1,LSDptr1,len2,LSDptr2, MSDptr_assignment,len_assignment,LSDptr_assignment LSDptr0 = ); \
709 if ((sintD)(MSDptr1[0]) < 0) # n<0 ? \
710 # muss m bzw. m+b^l subtrahieren, um k Digits verschoben: \
711 subfrom_loop_down(LSDptr2,&LSDptr0[-(uintP)len1],len2); \
712 if ((sintD)(MSDptr2[0]) < 0) # m<0 ? \
713 # muss n bzw. n+b^k subtrahieren, um l Digits verschoben: \
714 subfrom_loop_down(LSDptr1,&LSDptr0[-(uintP)len2],len1); \
715 }
716
717 # (* x y), wo x und y Integers sind. Ergebnis Integer.
718 # can trigger GC
719 # Methode:
720 # x=0 oder y=0 -> Ergebnis 0
721 # x und y beide Fixnums -> direkt multiplizieren
722 # sonst: zu WS machen, multiplizieren.
I_I_mult_I(object x,object y)723 local maygc object I_I_mult_I (object x, object y)
724 {
725 if (eq(x,Fixnum_0) || eq(y,Fixnum_0))
726 return Fixnum_0;
727 if (I_fixnump(x) && I_fixnump(y)) {
728 var sintV x_ = FN_to_V(x);
729 var sintV y_ = FN_to_V(y);
730 #if (oint_data_len+1 > 32)
731 # nur falls x und y Integers mit höchstens 32 Bit sind:
732 if (((uintV)((sintV)FN_sign(x) ^ x_) < vbit(31))
733 && ((uintV)((sintV)FN_sign(y) ^ y_) < vbit(31)))
734 #endif
735 {
736 # Werte direkt multiplizieren:
737 var uint32 hi;
738 var uint32 lo;
739 mulu32((uint32)x_,(uint32)y_,hi=,lo=); # erst unsigned multiplizieren
740 if (x_ < 0)
741 hi -= (uint32)y_; # dann Korrektur für Vorzeichen
742 if (y_ < 0)
743 hi -= (uint32)x_; # (vgl. DS_DS_mult_DS)
744 return L2_to_I(hi,lo);
745 }
746 }
747 var uintD* xMSDptr;
748 var uintC xlen;
749 var uintD* xLSDptr;
750 I_to_NDS_nocopy(x, xMSDptr = , xlen = , xLSDptr = );
751 var uintD* yMSDptr;
752 var uintC ylen;
753 var uintD* yLSDptr;
754 I_to_NDS_nocopy(y, yMSDptr = , ylen = , yLSDptr = );
755 {
756 SAVE_NUM_STACK # num_stack retten
757 var uintD* ergMSDptr;
758 var uintC erglen;
759 begin_arith_call();
760 DS_DS_mult_DS(xMSDptr,xlen,xLSDptr,yMSDptr,ylen,yLSDptr, ergMSDptr=,erglen=,);
761 end_arith_call();
762 var object result = DS_to_I(ergMSDptr,erglen);
763 RESTORE_NUM_STACK # num_stack zurück
764 return result;
765 }
766 }
767
768 # (EXPT x y), wo x Integer, y Integer >0 ist.
769 # can trigger GC
770 local maygc object I_I_expt_I (object x, object y);
771 # Methode:
772 # a:=x, b:=y, c:=1. [a^b*c bleibt invariant, = x^y.]
773 # Solange b>1,
774 # falls b ungerade, setze c:=a*c,
775 # setze b:=floor(b/2),
776 # setze a:=a*a.
777 # Wenn b=1, setze c:=a*c.
778 # Liefere c.
779 # Oder optimiert:
780 # a:=x, b:=y.
781 # Solange b gerade, setze a:=a*a, b:=b/2. [a^b bleibt invariant, = x^y.]
782 # c:=a.
783 # Solange b:=floor(b/2) >0 ist,
784 # setze a:=a*a, und falls b ungerade, setze c:=a*c.
785 # Liefere c.
I_I_expt_I(object x,object y)786 local maygc object I_I_expt_I (object x, object y)
787 {
788 #if 0 # unoptimiert
789 pushSTACK(x); pushSTACK(Fixnum_1); pushSTACK(y);
790 # Stackaufbau: a, c, b.
791 while (!eq(STACK_0,Fixnum_1)) { # solange bis b=1
792 if (I_oddp(STACK_0)) # b ungerade?
793 STACK_1 = I_I_mult_I(STACK_2,STACK_1); # c:=a*c
794 STACK_0 = I_I_ash_I(STACK_0,Fixnum_minus1); # b := (ash b -1) = (floor b 2)
795 STACK_2 = I_square_I(STACK_2); # a:=a*a
796 }
797 skipSTACK(1);
798 var object c = popSTACK();
799 var object a = popSTACK();
800 return I_I_mult_I(a,c); # a*c als Ergebnis
801 #else # optimiert
802 pushSTACK(x); pushSTACK(y);
803 # Stackaufbau: a, b.
804 while (!I_oddp(y)) {
805 STACK_1 = I_square_I(STACK_1); # a:=a*a
806 STACK_0 = y = I_I_ash_I(STACK_0,Fixnum_minus1); # b := (ash b -1)
807 }
808 pushSTACK(STACK_1); # c:=a
809 # Stackaufbau: a, b, c.
810 while (!eq(y=STACK_1,Fixnum_1)) { # Solange b/=1
811 STACK_1 = I_I_ash_I(y,Fixnum_minus1); # b := (ash b -1)
812 var object a = STACK_2 = I_square_I(STACK_2); # a:=a*a
813 if (I_oddp(STACK_1))
814 STACK_0 = I_I_mult_I(a,STACK_0); # evtl. c:=a*c
815 }
816 var object erg = STACK_0;
817 skipSTACK(3);
818 return erg;
819 #endif
820 }
821
822 # Fakultät (! n), wo n Fixnum >=0 ist. Ergebnis Integer.
823 # can trigger GC
824 # Methode:
825 # n <= 10 -> Ergebnis (Fixnum) aus Tabelle
826 # Sonst:
827 # Zweierpotenzen extra am Schluss durch einen Shift um
828 # ord2(n!) = sum(k>=1, floor(n/2^k) ) = n - logcount(n) Bits.
829 # Für k>=1 wird jede ungerade Zahl m im Intervall n/2^k < m <= n/2^(k-1)
830 # genau k mal gebraucht (als ungerader Anteil von m*2^0,...,m*2^(k-1) ).
831 # Zur Bestimmung des Produkts aller ungeraden Zahlen in einem Intervall
832 # a < m <= b verwenden wir eine rekursive Funktion, die nach Divide-and-
833 # Conquer das Produkt über die Intervalle a < m <= c und c < m <= b
834 # (c := floor((a+b)/2)) bestimmt und beide zusammenmultipliziert. Dies
835 # vermeidet, dass oft große Zahlen mit ganz kleinen Zahlen multipliziert
836 # werden.
837 local maygc object FN_fak_I (object n);
838 # UP für Fakultät:
839 # Bilde das Produkt prod(a < i <= b, 2*i+1), wobei 0 <= a < b klein.
prod_ungerade(uintV a,uintV b)840 local object prod_ungerade (uintV a, uintV b)
841 {
842 var uintV diff = b-a; # Anzahl der Faktoren
843 if (diff <= 4) {
844 # Produkt iterativ bilden
845 var object faktor = fixnum(2*b+1); # 2*b+1 als letzter Faktor
846 var object produkt = faktor;
847 var uintC count;
848 dotimesC(count,diff-1, {
849 faktor = fixnum_inc(faktor,-2); # nächster Faktor
850 produkt = I_I_mult_I(faktor,produkt); # mit bisherigem Produkt multiplizieren
851 });
852 return produkt;
853 } else {
854 # Produkt rekursiv bilden
855 var uintV c = floor(a+b,2); # c:=floor((a+b)/2)
856 var object teil = prod_ungerade(a,c); # erstes Teilprodukt
857 pushSTACK(teil);
858 teil = prod_ungerade(c,b); # zweites Teilprodukt
859 return I_I_mult_I(popSTACK(),teil); # und beide multiplizieren
860 }
861 }
FN_fak_I(object n)862 local maygc object FN_fak_I (object n)
863 {
864 local var const uintV fakul_table [] = {
865 1UL,
866 1UL,
867 1UL*2,
868 #if (oint_data_len>=3)
869 1UL*2*3,
870 #if (oint_data_len>=5)
871 1UL*2*3*4,
872 #if (oint_data_len>=7)
873 1UL*2*3*4*5,
874 #if (oint_data_len>=10)
875 1UL*2*3*4*5*6,
876 #if (oint_data_len>=13)
877 1UL*2*3*4*5*6*7,
878 #if (oint_data_len>=16)
879 1UL*2*3*4*5*6*7*8,
880 #if (oint_data_len>=19)
881 1UL*2*3*4*5*6*7*8*9,
882 #if (oint_data_len>=22)
883 1UL*2*3*4*5*6*7*8*9*10,
884 #if (oint_data_len>=26)
885 1UL*2*3*4*5*6*7*8*9*10*11,
886 #if (oint_data_len>=29)
887 1UL*2*3*4*5*6*7*8*9*10*11*12,
888 #if (oint_data_len>=33)
889 ULL(1)*2*3*4*5*6*7*8*9*10*11*12*13,
890 #if (oint_data_len>=37)
891 ULL(1)*2*3*4*5*6*7*8*9*10*11*12*13*14,
892 #if (oint_data_len>=41)
893 ULL(1)*2*3*4*5*6*7*8*9*10*11*12*13*14*15,
894 #if (oint_data_len>=45)
895 ULL(1)*2*3*4*5*6*7*8*9*10*11*12*13*14*15*16,
896 #if (oint_data_len>=49)
897 ULL(1)*2*3*4*5*6*7*8*9*10*11*12*13*14*15*16*17,
898 #if (oint_data_len>=53)
899 ULL(1)*2*3*4*5*6*7*8*9*10*11*12*13*14*15*16*17*18,
900 #if (oint_data_len>=57)
901 ULL(1)*2*3*4*5*6*7*8*9*10*11*12*13*14*15*16*17*18*19,
902 #if (oint_data_len>=62)
903 ULL(1)*2*3*4*5*6*7*8*9*10*11*12*13*14*15*16*17*18*19*20,
904 #if (oint_data_len>=66)
905 ...
906 #endif
907 #endif
908 #endif
909 #endif
910 #endif
911 #endif
912 #endif
913 #endif
914 #endif
915 #endif
916 #endif
917 #endif
918 #endif
919 #endif
920 #endif
921 #endif
922 #endif
923 #endif
924 #endif
925 };
926 var uintV n_ = posfixnum_to_V(n);
927 if (n_ < sizeof(fakul_table)/sizeof(uintV)) {
928 return fixnum(fakul_table[n_]);
929 } else {
930 pushSTACK(Fixnum_1); # bisheriges Produkt := 1
931 pushSTACK(n); # n
932 pushSTACK(Fixnum_1); # k := 1
933 pushSTACK(n); # obere Intervallgrenze floor(n/2^(k-1))
934 while (1) {
935 # Stackaufbau: prod, n, k, floor(n/2^(k-1)).
936 # 'n' enthält floor(n/2^(k-1)).
937 n = I_I_ash_I(n,Fixnum_minus1); # untere Grenze floor(n/2^k)
938 # 'n' enthält floor(n/2^k).
939 # Bilde Teilprodukt prod(A < i <= B & oddp(i), i)
940 # = prod(floor((A-1)/2) < i <= floor((B-1)/2), 2*i+1)
941 # wobei B = floor(n/2^(k-1)), A = floor(n/2^k) = floor(B/2).
942 {
943 var uintV b = floor(posfixnum_to_V(STACK_0)-1,2);
944 if (b==0) # B=2 oder B=1 -> Produkt fertig
945 break;
946 var uintV a = floor(posfixnum_to_V(n)-1,2);
947 pushSTACK(n);
948 var object temp = prod_ungerade(a,b);
949 temp = I_I_expt_I(temp,STACK_2); # hoch k nehmen
950 STACK_4 = I_I_mult_I(temp,STACK_4); # und aufmultiplizieren
951 }
952 STACK_2 = fixnum_inc(STACK_2,1); # k:=k+1
953 n = popSTACK(); STACK_0 = n;
954 }
955 skipSTACK(2);
956 # Stackaufbau: prod, n.
957 var object temp = I_logcount_I(STACK_0); # (logcount n)
958 temp = I_I_minus_I(popSTACK(),temp); # (- n (logcount n))
959 return I_I_ash_I(popSTACK(),temp); # (ash prod (- n (logcount n)))
960 }
961 }
962
963