home *** CD-ROM | disk | FTP | other *** search
/ Amiga ISO Collection / AmigaUtilCD2.iso / Programming / Misc / CLISP-1.LHA / CLISP960530-sr.lha / src / int2adic.d < prev    next >
Encoding:
Text File  |  1996-04-15  |  4.6 KB  |  141 lines
  1. # Operationen mit 2-adischen Integers
  2.  
  3. # Multipliziert zwei Zahlen mod 2^intDsize.
  4. # D_D_mal2adic_D(a,b)
  5. # > uintD a,b: Zahlen mod 2^intDsize
  6. # < ergebnis: Zahl c mod 2^intDsize mit c == a*b mod 2^intDsize
  7.   local uintD D_D_mal2adic_D (uintD a, uintD b);
  8. #if HAVE_DD
  9.   #define D_D_mal2adic_D(a,b)  lowD(muluD((uintD)(a),(uintD)(b)))
  10. #else
  11.   #ifdef GNU
  12.     #define D_D_mal2adic_D(a,b)  \
  13.       ({ var reg1 uintD __erg;    \
  14.          muluD((a),(b),_EMA_,__erg=); \
  15.          __erg;                   \
  16.        })
  17.   #else
  18.     #if (intDsize==32)
  19.       #define D_D_mal2adic_D(a,b)  mulu32_(a,b) # lo-Teil von mulu32(a,b)
  20.     #else
  21.       local uintD D_D_mal2adic_D(a,b)
  22.         var reg1 uintD a;
  23.         var reg2 uintD b;
  24.         { muluD(a,b,_EMA_,return); }
  25.     #endif
  26.   #endif
  27. #endif
  28.  
  29. # Potenziert eine Zahl mod 2^intDsize.
  30. # D_UL_expt_D(x,y)
  31. # > uintD x: Zahl mod 2^intDsize
  32. # > uintL y: Exponent >0
  33. # < uintD ergebnis: x^y mod 2^intDsize
  34.   local uintD D_UL_expt_D (uintD x, uintL y);
  35. # Methode:
  36. #   a:=x, b:=y, c:=1. [a^b*c bleibt invariant, = x^y.]
  37. #   Solange b>1,
  38. #     falls b ungerade, setze c:=a*c,
  39. #     setze b:=floor(b/2),
  40. #     setze a:=a*a.
  41. #   Wenn b=1, setze c:=a*c.
  42. #   Liefere c.
  43. # Oder optimiert:
  44. #   a:=x, b:=y.
  45. #   Solange b gerade, setze a:=a*a, b:=b/2. [a^b bleibt invariant, = x^y.]
  46. #   c:=a.
  47. #   Solange b:=floor(b/2) >0 ist,
  48. #     setze a:=a*a, und falls b ungerade, setze c:=a*c.
  49. #   Liefere c.
  50.   local uintD D_UL_expt_D(a,b)
  51.     var reg1 uintD a;
  52.     var reg2 uintL b;
  53.     { while ((b & bit(0)) ==0) { a = D_D_mal2adic_D(a,a); b = b>>1; }
  54.      {var reg3 uintD c = a;
  55.       until ((b = b>>1) == 0)
  56.         { a = D_D_mal2adic_D(a,a);
  57.           if (b & bit(0)) { c = D_D_mal2adic_D(a,c); }
  58.         }
  59.       return c;
  60.     }}
  61.  
  62. # Dividiert zwei Zahlen mod 2^intDsize.
  63. # D_D_durch2adic_D(a,b)
  64. # > uintD a: Zahl mod 2^intDsize
  65. # > uintD b: ungerade Zahl mod 2^intDsize
  66. # < ergebnis: Zahl c mod 2^intDsize mit b*c == a mod 2^intDsize
  67.   local uintD D_D_durch2adic_D (uintD a, uintD b);
  68. # Methode:
  69. # Konstruiere c Bit für Bit.
  70. # c := 0, d := a.
  71. # Für j=0,...,intDsize:
  72. #   [Hier b*c == a mod 2^j und d = (a-b*c)/2^j.] j=intDsize -> fertig.
  73. #   Falls d ungerade, setze c:=c+2^j und d:=(d-b)/2, sonst d:=d/2.
  74. # Ergebnis c.
  75. #if 1
  76.   local uintD D_D_durch2adic_D(a,b)
  77.     var reg1 uintD a;
  78.     var reg4 uintD b;
  79.     { ASSERT(!((b % 2) ==0))
  80.      {var reg3 uintD c = 0;
  81.       var reg2 uintD bit_j = 1; # 2^j
  82.       loop # Verwende a als Variable d
  83.         { if (a & bit(0)) { c = c+bit_j; a = a-b; }
  84.           a = a>>1;
  85.           bit_j = bit_j << 1;
  86.           if (bit_j == 0) break; # j=intDsize -> fertig
  87.         }
  88.       return c;
  89.     }}
  90. #else
  91.   local uintD D_D_durch2adic_D(a,b)
  92.     var reg1 uintD a;
  93.     var reg4 uintD b;
  94.     { ASSERT(!((b % 2) ==0))
  95.      {var reg2 uintD bit_j = 1; # 2^j
  96.       var reg3 uintD b_j = b-1; # (b-1)*2^j
  97.       loop # Verwende a als Variable d*2^j+c
  98.         { if (a & bit_j) { a = a - b_j; }
  99.           b_j = b_j << 1; bit_j = bit_j << 1;
  100.           if (bit_j == 0) break; # j=intDsize -> fertig
  101.         }
  102.       return a;
  103.     }}
  104. #endif
  105.  
  106. # UDS_UDS_durch2adic_UDS(len,a_LSDptr,b_LSDptr,dest_LSDptr);
  107. # dividiert die UDS a_LSDptr[-len..-1] mod 2^(intDsize*len)
  108. # durch die ungerade UDS b_LSDptr[-len..-1] mod 2^(intDsize*len) (len>0)
  109. # und liefert den Quotienten als UDS dest_LSDptr[-len..-1] mod 2^(intDsize*len).
  110.   local void UDS_UDS_durch2adic_UDS (uintC len, uintD* a_LSDptr, uintD* b_LSDptr, uintD* dest_LSDptr);
  111. # Methode: beta=2^intDsize. Schreibe jeweils x = x[0]*beta^0 + x[1]*beta^1 + ... .
  112. # Um b*c == a mod beta^m zu bestimmen:
  113. # Sei b' := (b mod beta)^(-1) = b[0]^(-1), d := a.
  114. # Für j=0,...,m:
  115. #   [Hier ist d = a - b*c == 0 mod beta^j, und c mod beta^j bekannt.]
  116. #   j=m -> fertig.
  117. #   Setze c[j] := b'*d[j] mod beta, also c := c + (b'*d[j] mod beta)*beta^j.
  118. #   Setze d := d - b * c[j] * beta^j.
  119. # Schließlich dest := c.
  120.   local void UDS_UDS_durch2adic_UDS(len,a_LSDptr,b_LSDptr,dest_LSDptr)
  121.     var reg3 uintC len;
  122.     var reg6 uintD* a_LSDptr;
  123.     var reg5 uintD* b_LSDptr;
  124.     var reg2 uintD* dest_LSDptr;
  125.     { var reg4 uintD b0inv = D_D_durch2adic_D(1,b_LSDptr[-1]); # b'
  126.       copy_loop_down(a_LSDptr,dest_LSDptr,len); # d := a
  127.       do { var reg1 uintD digit = dest_LSDptr[-1]; # nächstes d[j]
  128.            #if HAVE_DD
  129.              digit = lowD(muluD(b0inv,digit));
  130.            #else
  131.              muluD(b0inv,digit,_EMA_,digit=);
  132.            #endif
  133.            # digit = nächstes c[j]
  134.            mulusub_loop_down(digit,b_LSDptr,dest_LSDptr,len); # d := d - b * c[j] * beta^j
  135.            # Nun ist dest_LSDptr[-1] = 0.
  136.            *--dest_LSDptr = digit; len--; # c[j] ablegen, nächstes j
  137.          }
  138.          until (len==0);
  139.     }
  140.  
  141.