4 #include "base/cl_sysdep.h"
7 #include "cln/rational.h"
12 #include "integer/cl_I.h"
13 #include "rational/cl_RA.h"
14 #include "base/cl_xmacros.h"
18 bool logp (const cl_I& a, const cl_I& b, cl_RA* l)
21 // log(a,b) soll Bruch c/d mit teilerfremdem c>=0,d>0 ergeben.
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
37 // rechts mit ( 1 ) multiplizieren.
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.
42 // D.h. die bisherige Matrix von links ist ( vc vd ).
43 // uc:=1, ud:=0, vc:=0, vd:=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.
56 if (eq(a,1)) // a=1 -> Rekursion zu Ende
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;
65 // 1<a<b -> a und b vertauschen
67 swap(uintL, uc, ud); swap(uintL, vc, vd);
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;