1 // FILE SYMB.CC: Implementations for symbols
2 //////////////////////////////////////////////////////////////////////////
3 //
4 // Copyright 1990-2012 John Cremona
5 //
6 // This file is part of the eclib package.
7 //
8 // eclib is free software; you can redistribute it and/or modify it
9 // under the terms of the GNU General Public License as published by the
10 // Free Software Foundation; either version 2 of the License, or (at your
11 // option) any later version.
12 //
13 // eclib is distributed in the hope that it will be useful, but WITHOUT
14 // ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
15 // FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
16 // for more details.
17 //
18 // You should have received a copy of the GNU General Public License
19 // along with eclib; if not, write to the Free Software Foundation,
20 // Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA
21 //
22 //////////////////////////////////////////////////////////////////////////
23
24 #include <eclib/symb.h>
25
26 // Friend of class symb:
operator <<(ostream & s,const symb & sy)27 ostream& operator<< (ostream& s, const symb& sy)
28 {
29 s << "(" << sy.c << ":" << sy.d << ")";
30 return s;
31 }
32
33 //#define DEBUG_NORMALIZE
normalize() const34 symb symb::normalize() const
35 {
36 #ifdef DEBUG_NORMALIZE
37 cout<<"Normalizing symbol "<<(*this)<<endl;
38 #endif
39 long n=N->modulus;
40 long u=N->unitdiv(c);
41 #ifdef DEBUG_NORMALIZE
42 cout<<"scaling by u = "<<u<<endl;
43 #endif
44 long cc=N->reduce(xmodmul(c,u,n));
45 long dd=N->reduce(xmodmul(d,u,n))%(n/cc);
46 #ifdef DEBUG_NORMALIZE
47 cout<<"new c = "<<cc<<endl;
48 cout<<"new d = "<<dd<<endl;
49 #endif
50 symb ans(cc,dd,N);
51 #ifdef DEBUG_NORMALIZE
52 cout<<"Returning normalized symbol "<<ans;
53 int ok = (ans==(*this)) && ::div(ans.cee(),n);
54 if(ok) cout<<" ok"; else cout<<" wrong!";
55 cout<<endl;
56 #endif
57 return ans;
58 }
59
60 // Constructor for modsym, converting from symb:
modsym(const symb & s)61 modsym::modsym(const symb& s)
62 {
63 long c,d,h,x,y;
64 c = s.cee(); d = s.dee();
65 h = bezout(c , d, x, y);
66 a=rational(-x , d/h);
67 b=rational( y , c/h);
68 }
69
70 // Friend of class modsym:
operator <<(ostream & s,const modsym & m)71 ostream& operator<< (ostream& s, const modsym& m)
72 {
73 s << "{" << (m.a) << "," << (m.b) << "}";
74 return s;
75 }
76
77 //Members of class symblist:
78
symblist(long n)79 symblist::symblist(long n)
80 {
81 maxnum=n;
82 num=0;
83 list=new symb[n];
84 }
85
~symblist()86 symblist::~symblist()
87 {
88 delete[] list;
89 }
90
91
add(const symb & s,long start)92 void symblist::add(const symb& s, long start)
93 {
94 if (index(s,start)==-1)
95 {
96 if (num<maxnum)
97 {
98 list[num]=s;
99 long c = s.cee(), d=posmod(s.dee(),s.modulus()/c);
100 hashtable[pair<long,long>(c,d)]=num;
101 num++;
102 // cout<<"Adding symbol "<<s<<" as special number "<<num<<endl;
103 }
104 else
105 {
106 cerr << "Error in symblist::add: attempt to add too many symbols to list!"<<endl;
107 }
108 }
109 }
110
111 /*
112 long symblist::index(const symb& s, long start) const
113 {
114 long i,ans;
115 for (i=start,ans=-1; ((i<num)&&(ans==-1)); i++) if (list[i]==s) ans=i;
116 return ans;
117 }
118 */
119
index(const symb & s,long start) const120 long symblist::index(const symb& s, long start) const
121 {
122 // cout<<"index of "<<s;
123 symb ss = s.normalize();
124 long c = ss.cee(), d=ss.dee();
125 map<pair<long,long>,long>::const_iterator
126 j = hashtable.find(pair<long,long>(c,d));
127 if(j==hashtable.end())
128 return -1;
129 // cout<<" is "<<j->second<<endl;
130 return j->second;
131 }
132
133
item(long n) const134 symb symblist::item(long n) const
135 {
136 if ((n>num)||(n<0))
137 {
138 cerr<<"Error in symblist::item: index out of range!"<<endl;
139 return symb();
140 }
141 else return list[n];
142 }
143
144 //Member functions for class symbdata:
symbdata(long n)145 symbdata::symbdata(long n) :moddata(n),specials(nsymb2)
146 {
147 // cout << "In constructor symbdata::symbdata.\n";
148 // cout << "nsymb2 = " << nsymb2 << "\n";
149 if (nsymb2>0)
150 { long ic,id,c,d,start; symb s;
151 //N.B. dlist include d=1 at 0 and d=mod at end, which we don't want here
152 for (ic=1; (ic<ndivs-1)&&(specials.count()<nsymb2); ic++)
153 { c=dlist[ic];
154 dstarts[ic]=start=specials.count();
155 for (id=1; (id<modulus-phi)&&(specials.count()<nsymb2); id++)
156 { d = noninvlist[id];
157 if (::gcd(d,c)==1)
158 { s = symb(c,d,this);
159 specials.add(s,start); //only adds it if not there already!
160 }
161 } // end of d loop
162 } // end of c loop
163 if (specials.count()<nsymb2)
164 {
165 cout << "Problem: makesymbols found only " << specials.count() << " symbols ";
166 cout << "out of " << nsymb2 << endl;
167 }
168 // cout << "Special symbols: "; specials.display();
169 }
170 }
171
index2(long c,long d) const172 long symbdata::index2(long c, long d) const
173 { long kd = code(d);
174 // cout<<"index2("<<c<<":"<<d<<"):"<<endl;
175 if (kd>0) // d invertible, with inverse kd
176 return reduce(xmodmul(c,kd,modulus)); // (c:d) = (c*kd:1)
177 else
178 { long kc = code(c);
179 if (kc>0) // (c:d) = (1:kc*d) if c invertible
180 return modulus-code(xmodmul(kc,d,modulus));
181 else
182 {
183 long start = dstarts[noninvdlist[-kc]];
184 symb s(c,d,this);
185 long ind = specials.index(s,start);
186 if(ind<0)
187 {
188 cout<<"error in index(): symbol "<<s<<" not in list!"<<endl;
189 }
190 return nsymb1+ind;
191 }
192 }
193 }
194
symbol(long i) const195 symb symbdata::symbol(long i) const
196 { if (i<modulus) return symb(i,1,this);
197 else if (i<nsymb1) return symb(1,noninvlist[i-modulus],this);
198 else return specials[i-nsymb1]; // specials.item[i-nsymb1];
199 }
200
display() const201 void symbdata::display() const
202 { moddata::display();
203 cout << "Number of special symbols = " << nsymb2 << "\n";
204 specials.display();
205 }
206
check(void) const207 void symbdata::check(void) const
208 {long i,j; int ok=1; symb s;
209 for (i=0; i<nsymb; i++)
210 {j = index(s=symbol(i));
211 if (i!=j)
212 {
213 cout << i << "-->" << s << "-->" << j << "\n";
214 ok=0;
215 }
216 }
217 if (ok) cout << "symbols check OK!\n";
218 else cout << "symbols check found errors!\n";
219 }
220
jumpsymb(symb s1,symb s2)221 modsym jumpsymb(symb s1, symb s2)
222 {
223 //Assuming s1==s2, returns closed modular symbol {g1(0),g2(0)} where gi<->si
224 long c1=s1.cee(), c2=s2.cee(), d1=s1.dee(), d2=s2.dee();
225 return modsym(rational(-invmod(c1,d1),d1),rational(-invmod(c2,d2),d2));
226 }
227
228