src/rational/transcendental/cl_I_logp.cc

   1 // logp().
   2
   3 // General includes.
   4 #include "cl_sysdep.h"
   5
   6 // Specification.
   7 #include "cln/rational.h"
   8
   9
  10 // Implementation.
  11
  12 #include "cl_I.h"
  13 #include "cl_RA.h"
  14 #include "cl_xmacros.h"
  15
  16 namespace cln {
  17
  18 bool logp (const cl_I& a, const cl_I& b, cl_RA* l)
  19 {
  20 // Methode:
  21 //   log(a,b) soll Bruch c/d mit teilerfremdem c>=0,d>0 ergeben.
  22 //   a=1 -> c=0, d=1.
  23 //   a>=b -> Dividiere a durch b. Rest da -> geht nicht.
  24 //           Sonst log(a,b) = 1+log(a/b,b).
  25 //           Berechne also c/d := log(a/b,b) und setze c:=c+d.
  26 //   1<a<b -> log(a,b) = 1/log(b,a).
  27 //           Berechne also c/d := log(b,a) und vertausche c und d.
  28 // Man konstruiert hierbei eigentlich die Kettenbruchentwicklung von c/d.
  29 // Wegen a>=2^c, b>=2^d sind c,d < (integer-length a,b) < intDsize*2^intCsize.
  30 // In Matrizenschreibweise:
  31 //   Wenn eine Folge von Divisionsschritten D und Vertauschungsschritten V
  32 //   ausgeführt werden muß, z.B. (a,b) V D D = (1,*), so ist
  33 //     ( c )           ( 0 )             ( 1 1 )           ( 0 1 )
  34 //     ( d )  =  V D D ( 1 )  wobei  D = ( 0 1 )  und  V = ( 1 0 ).
  35 //   Man baut diese Matrizen nun von links nach rechts auf, zum Schluß von
  36 //              ( 0 )
  37 //   rechts mit ( 1 ) multiplizieren.
  38 // Entrekursiviert:
  39 //   Wir werden (a,b) und damit auch c/d = log(a/b) verändern.
  40 //   Invariante: Statt (c,d) wollen wir (uc*c+ud*d,vc*c+vd*d) zurückliefern.
  41 //                                           ( uc ud )
  42 //   D.h. die bisherige Matrix von links ist ( vc vd ).
  43 //   uc:=1, ud:=0, vc:=0, vd:=1.
  44 //   Solange a>1,
  45 //     a>=b -> Dividiere a durch b. Rest da -> geht nicht.
  46 //             Sonst a:=a/b, und (für später c:=c+d) ud:=uc+ud, vd:=vc+vd.
  47 //     1<a<b -> vertausche a und b, uc und ud, vc und vd.
  48 //   Liefere (ud,vd), der Bruch ud/vd ist gekürzt.
  49  {      Mutable(cl_I,a);
  50         Mutable(cl_I,b);
  51         var uintL uc = 1;
  52         var uintL ud = 0;
  53         var uintL vc = 0;
  54         var uintL vd = 1;
  55         loop {
  56                 if (eq(a,1)) // a=1 -> Rekursion zu Ende
  57                         break;
  58                 if (a >= b) {
  59                         var cl_I_div_t div = cl_divide(a,b); // a durch b dividieren
  60                         if (!eq(div.remainder,0)) // Rest /=0 ?
  61                                 return false; // -> fertig
  62                         a = div.quotient;  // a := a/b
  63                         ud = uc + ud; vd = vc + vd;
  64                 } else {
  65                         // 1<a<b -> a und b vertauschen
  66                         swap(cl_I, a, b);
  67                         swap(uintL, uc, ud); swap(uintL, vc, vd);
  68                 }
  69         }
  70         // a=1 -> c=0,d=1 -> Ergebnis ud/vd
  71         *l = I_I_to_RA(UL_to_I(ud),UL_to_I(vd)); return true;
  72 }}
  73
  74 }  // namespace cln