X-Git-Url: https://ginac.de/CLN/cln.git//cln.git?a=blobdiff_plain;f=src%2Fbase%2Fdigitseq%2Fcl_DS_mul_kara.h;h=3682a6d7c868ea7e24123cf5e2303a622ce5a231;hb=22549ef70fee95faab1e9f2adaf710ba9e0bdabf;hp=8531b5b1256cba0a43d89d377659105e36df8062;hpb=5370ad8054201cf23d4f94a52f4d3f7f9f3cd511;p=cln.git diff --git a/src/base/digitseq/cl_DS_mul_kara.h b/src/base/digitseq/cl_DS_mul_kara.h index 8531b5b..3682a6d 100644 --- a/src/base/digitseq/cl_DS_mul_kara.h +++ b/src/base/digitseq/cl_DS_mul_kara.h @@ -6,14 +6,14 @@ // = x1*y1 * b^2k + ((x1+x0)*(y1+y0)-x1*y1-x0*y0) * b^k + x0*y0 // Methode 1 (Collins/Loos, Degel): // source2 wird in floor(len2/len1) einzelne UDS mit je einer - // Länge len3 (len1 <= len3 < 2*len1) unterteilt, + // Länge len3 (len1 <= len3 < 2*len1) unterteilt, // jeweils k=floor(len3/2). // Methode 2 (Haible): // source2 wird in ceiling(len2/len1) einzelne UDS mit je einer - // Länge len3 (0 < len3 <= len1) unterteilt, jeweils k=floor(len1/2). - // Aufwand für die hinteren Einzelteile: + // Länge len3 (0 < len3 <= len1) unterteilt, jeweils k=floor(len1/2). + // Aufwand für die hinteren Einzelteile: // bei beiden Methoden jeweils 3*len1^2. - // Aufwand für das vorderste Teil (alles, falls len1 <= len2 < 2*len1) + // Aufwand für das vorderste Teil (alles, falls len1 <= len2 < 2*len1) // mit r = len1, s = (len2 mod len1) + len1 (>= len1, < 2*len1): // bei Methode 1: // | : | r @@ -24,23 +24,23 @@ // | | : | s // (s-r)*r + r/2*r/2 + r/2*r/2 + r/2*r/2 = r*s - r^2/4 . // Wegen (r*s/2 + s^2/4) - (r*s - r^2/4) = (r-s)^2/4 >= 0 - // ist Methode 2 günstiger. - // Denkfehler! Dies gilt - wenn überhaupt - nur knapp oberhalb des + // ist Methode 2 günstiger. + // Denkfehler! Dies gilt - wenn überhaupt - nur knapp oberhalb des // Break-Even-Points. - // Im allgemeinen ist der Multiplikationsaufwand für zwei Zahlen der - // Längen u bzw. v nämlich gegeben durch min(u,v)^c * max(u,v), + // Im allgemeinen ist der Multiplikationsaufwand für zwei Zahlen der + // Längen u bzw. v nämlich gegeben durch min(u,v)^c * max(u,v), // wobei c = log3/log2 - 1 = 0.585... - // Dadurch wird der Aufwand in Abhängigkeit des Parameters t = k, + // Dadurch wird der Aufwand in Abhängigkeit des Parameters t = k, // r/2 <= t <= s/2 (der einzig sinnvolle Bereich), zu // (r-t)^c*(s-t) + t^c*(s-t) + t^(1+c). // Dessen Optimum liegt (im Bereich r <= s <= 2*r) - // - im klassischen Fall c=1 tatsächlich stets bei t=r/2 [Methode 2], + // - im klassischen Fall c=1 tatsächlich stets bei t=r/2 [Methode 2], // - im Karatsuba-Fall c=0.6 aber offenbar bei t=s/2 [Methode 1] // oder ganz knapp darunter. // Auch erweist sich Methode 1 im Experiment als effizienter. // Daher implementieren wir Methode 1 : { // Es ist 2 <= len1 <= len2. - // Spezialfall Quadrieren abfangen (häufig genug, daß sich das lohnt): + // Spezialfall Quadrieren abfangen (häufig genug, daß sich das lohnt): if (sourceptr1 == sourceptr2) if (len1 == len2) { mulu_karatsuba_square(sourceptr1,len1,destptr); return; } @@ -48,8 +48,8 @@ if (len2 >= 2*len1) { CL_SMALL_ALLOCA_STACK; // Teilprodukte von jeweils len1 mal len1 Digits bilden: - var uintC k_lo = floor(len1,2); // Länge der Low-Teile: floor(len1/2) >0 - var uintC k_hi = len1 - k_lo; // Länge der High-Teile: ceiling(len1/2) >0 + var uintC k_lo = floor(len1,2); // Länge der Low-Teile: floor(len1/2) >0 + var uintC k_hi = len1 - k_lo; // Länge der High-Teile: ceiling(len1/2) >0 // Es gilt k_lo <= k_hi <= len1, k_lo + k_hi = len1. // Summe x1+x0 berechnen: var uintD* sum1_MSDptr; @@ -65,11 +65,11 @@ } if (carry) { lsprefnext(sum1_MSDptr) = 1; sum1_len++; } } - { // Platz für Summe y1+y0 belegen: + { // Platz für Summe y1+y0 belegen: var uintC sum2_maxlen = k_hi+1; var uintD* sum2_LSDptr; num_stack_small_alloc(sum2_maxlen,,sum2_LSDptr=); - // Platz für Produkte x0*y0, x1*y1 belegen: + // Platz für Produkte x0*y0, x1*y1 belegen: { var uintD* prod_MSDptr; var uintD* prod_LSDptr; var uintD* prodhi_LSDptr; @@ -79,11 +79,11 @@ // Produkte x1*y1 in prod_MSDptr/2*k_hi/prodhi_LSDptr // und x0*y0 in prodhi_LSDptr/2*k_lo/prod_LSDptr, // dann das Produkt (b^k*x1+x0)*(b^k*y1+y0) enthalten. - // Platz fürs Produkt (x1+x0)*(y1+y0) belegen: + // Platz fürs Produkt (x1+x0)*(y1+y0) belegen: {var uintD* prodmid_MSDptr; var uintD* prodmid_LSDptr; num_stack_small_alloc(sum1_len+sum2_maxlen,prodmid_MSDptr=,prodmid_LSDptr=); - // Schleife über die hinteren Einzelteile: + // Schleife über die hinteren Einzelteile: do { // Produkt x0*y0 berechnen: cl_UDS_mul(sourceptr1,k_lo,sourceptr2,k_lo,prod_LSDptr); // Produkt x1*y1 berechnen: @@ -118,7 +118,7 @@ if (!(carry==0)) { dec_loop_lsp(prodmid_LSDptr lspop 2*k_lo,prodmid_len-2*k_lo); } } - // prodmid_LSDptr[-prodmid_len..-1] enthält nun x0*y1+x1*y0. + // prodmid_LSDptr[-prodmid_len..-1] enthält nun x0*y1+x1*y0. // Dies wird zu prod = x1*y1*b^(2*k) + x0*y0 addiert: {var uintD carry = addto_loop_lsp(prodmid_LSDptr,prod_LSDptr lspop k_lo,prodmid_len); @@ -151,8 +151,8 @@ var uintC prod_len = len1+len2; var uintD* prod_LSDptr; num_stack_small_alloc(prod_len,prod_MSDptr=,prod_LSDptr=); - { var uintC k_hi = floor(len2,2); // Länge der High-Teile: floor(len2/2) >0 - var uintC k_lo = len2 - k_hi; // Länge der Low-Teile: ceiling(len2/2) >0 + { var uintC k_hi = floor(len2,2); // Länge der High-Teile: floor(len2/2) >0 + var uintC k_lo = len2 - k_hi; // Länge der Low-Teile: ceiling(len2/2) >0 // Es gilt k_hi <= k_lo <= len1 <= len2, k_lo + k_hi = len2. var uintC x1_len = len1-k_lo; // <= len2-k_lo = k_hi <= k_lo // Summe x1+x0 berechnen: @@ -183,14 +183,14 @@ } if (carry) { lsprefnext(sum2_MSDptr) = 1; sum2_len++; } } - // Platz für Produkte x0*y0, x1*y1: + // Platz für Produkte x0*y0, x1*y1: { var uintC prodhi_len = x1_len+k_hi; var uintD* prodhi_LSDptr = prod_LSDptr lspop 2*k_lo; // prod_MSDptr/len1+len2/prod_LSDptr wird zuerst die beiden // Produkte x1*y1 in prod_MSDptr/x1_len+k_hi/prodhi_LSDptr // und x0*y0 in prodhi_LSDptr/2*k_lo/prod_LSDptr, // dann das Produkt (b^k*x1+x0)*(b^k*y1+y0) enthalten. - // Platz fürs Produkt (x1+x0)*(y1+y0) belegen: + // Platz fürs Produkt (x1+x0)*(y1+y0) belegen: {var uintD* prodmid_MSDptr; var uintC prodmid_len = sum1_len+sum2_len; var uintD* prodmid_LSDptr; @@ -220,13 +220,13 @@ // Carry um maximal 1 Digit weitertragen: if (!(carry==0)) { lspref(prodmid_LSDptr,2*k_lo) -= 1; } } - // prodmid_LSDptr[-prodmid_len..-1] enthält nun x0*y1+x1*y0. + // prodmid_LSDptr[-prodmid_len..-1] enthält nun x0*y1+x1*y0. // Dies ist < b^k_lo * b^k_hi + b^x1_len * b^k_lo // = b^len2 + b^len1 <= 2 * b^len2, - // paßt also in len2+1 Digits. + // paßt also in len2+1 Digits. // Im Fall x1_len=0 ist es sogar < b^k_lo * b^k_hi = b^len2, - // es paßt also in len2 Digits. - // prodmid_len, wenn möglich, um maximal 2 verkleinern: + // es paßt also in len2 Digits. + // prodmid_len, wenn möglich, um maximal 2 verkleinern: // (benutzt prodmid_len >= 2*k_lo >= len2 >= 2) if (mspref(prodmid_MSDptr,0)==0) { prodmid_len--;