1 /*
2  * DeepThought: This package provides functions for multiplication and other computations in finitely generated nilpotent groups based on the Deep Thought algorithm.
3  */
4 
5 #include "src/compiled.h"          /* GAP headers */
6 
7 #define IS_INTPOS(obj)          (TNUM_OBJ(obj) == T_INTPOS)
8 #define IS_INTNEG(obj)          (TNUM_OBJ(obj) == T_INTNEG)
9 #define IS_LARGEINT(obj)        (IS_INTPOS(obj) || IS_INTNEG(obj))
10 
11 #define IS_INT(obj)             (IS_INTOBJ(obj) || IS_LARGEINT(obj))
12 
13 #define IS_NEGATIVE(obj)        (IS_INTOBJ(obj) ? ((Int)obj < 0) : IS_INTNEG(obj))
14 #define IS_POSITIVE(obj)        (IS_INTOBJ(obj) ? ((Int)obj > 1) : IS_INTPOS(obj))
15 #define IS_ODD(obj)             (IS_INTOBJ(obj) ? ((Int)obj & 4) : (VAL_LIMB0(obj) & 1))
16 #define IS_EVEN(obj)            (!IS_ODD(obj))
17 
18 #ifdef IS_PREC_REP
19 // for compatibility with GAP <= 4.8
20 #define IS_PREC(x) IS_PREC_REP(x)
21 #endif
22 
23 static UInt RNleft, RNright, RNlength, RNnum, RNside;
24 
25 
DTP_Binomial(Obj self,Obj N,Obj K)26 Obj DTP_Binomial(Obj self, Obj N, Obj K)
27 {
28     // handle some special cases
29     if (K == INTOBJ_INT(1)) // K=1 is the most frequent case for us, so check it first
30         return N;
31 
32     // restrict input
33     if (!IS_INT(N) || !IS_INT(K))
34         return Fail;
35 
36     if (K == INTOBJ_INT(0))
37         return INTOBJ_INT(1);
38     if (IS_NEGATIVE(K))
39         return INTOBJ_INT(0);
40 
41     // from now on, k >= 2
42 
43     if (K == N)
44         return INTOBJ_INT(1);
45 
46     if (IS_NEGATIVE(N))
47         return Fail;    // TODO: implement this (does it happen?)
48 
49     if (LtInt(N, K))
50         return INTOBJ_INT(0);
51 
52     // We only support immediate integers for k. Anything else at this point
53     // would lead to output that's far too big for storage anyway.
54     if (!IS_INTOBJ(K))
55         return Fail;
56 
57     Int k = INT_INTOBJ(K);
58     Int i;
59     Obj res = N;
60     Int quot = 1;
61 
62     // "unroll" the first few operations, to avoid repeated divisions, while
63     // hopefully keeping res small enough to be represented as an immediate
64     // (if N wasn't too big)
65     for (i = 2; i <= k && i <= 6; ++i) {
66         N = DiffInt(N, INTOBJ_INT(1));
67         res = ProdInt(res, N);
68         quot *= i;
69     }
70     res = QuoInt(res, INTOBJ_INT(quot));
71 
72     // now the general case
73     for (; i <= k; ++i) {
74         N = DiffInt(N, INTOBJ_INT(1));
75         res = ProdInt(res, N);
76         res = QuoInt(res, INTOBJ_INT(i));
77     }
78 
79     return res;
80 }
81 
DTP_SequenceLetter(Obj self,Obj letter,Obj seq)82 Obj DTP_SequenceLetter(Obj self, Obj letter, Obj seq)
83 {
84     if (!IS_PREC(letter))
85         ErrorMayQuit("DTP_SequenceLetter: <letter> must be a plain record (not a %s)",
86                  (Int)TNAM_OBJ(letter), 0L);
87 
88     if (!IS_PLIST(seq))
89         ErrorMayQuit("DTP_SequenceLetter: <seq> must be a plain list (not a %s)",
90                  (Int)TNAM_OBJ(seq), 0L);
91 
92     if (IsbPRec(letter, RNleft))
93         DTP_SequenceLetter(self, ElmPRec(letter, RNleft), seq);
94 
95     if (IsbPRec(letter, RNright))
96         DTP_SequenceLetter(self, ElmPRec(letter, RNright), seq);
97 
98     UInt len = LEN_PLIST(seq);
99     AssPlist(seq, len+1, letter);
100 
101     return 0;
102 }
103 
DTP_Seq_i(Obj self,Obj letter,Obj i)104 Obj DTP_Seq_i(Obj self, Obj letter, Obj i)
105 {
106     if (!IS_PREC(letter))
107     ErrorMayQuit("DTP_Seq_I: <letter> must be a plain record (not a %s)",
108              (Int)TNAM_OBJ(letter), 0L);
109 
110     if (!IS_INT(i))
111     ErrorMayQuit("DTP_Seq_i: <i> must be an integer (not a %s)",
112          (Int)TNAM_OBJ(letter), 0L);
113 
114     // We only support immediate integers for i. Anything else at this point
115     // would lead to output that's far too big for storage anyway.
116     if (!IS_INTOBJ(i))
117         return Fail;
118 
119     while( ElmPRec(letter, RNlength) > i ){
120         Obj left = ElmPRec(letter, RNleft);
121         if( ElmPRec(left, RNlength) >= i )
122         {
123             letter = left;
124         } else {
125             i = DiffInt( i, ElmPRec(left, RNlength) );
126             letter = ElmPRec(letter, RNright);
127         }
128     }
129 
130     if(ElmPRec(letter, RNlength) != i){
131         ErrorMayQuit("DTP_Seq_i: Assertion failure, letter.l <> i", 0, 0);
132     }
133 
134     return letter;
135 }
136 
DTP_AreAlmostEqual(Obj self,Obj letter1,Obj letter2)137 Obj DTP_AreAlmostEqual(Obj self, Obj letter1, Obj letter2)
138 {
139     if (!IS_PREC(letter1))
140     ErrorMayQuit("DTP_AreAlmostEqual: <letter1> must be a plain record (not a %s)",
141              (Int)TNAM_OBJ(letter1), 0L);
142 
143     if (!IS_PREC(letter2))
144     ErrorMayQuit("DTP_AreAlmostEqual: <letter2> must be a plain record (not a %s)",
145              (Int)TNAM_OBJ(letter2), 0L);
146 
147     // We must have letter1.num = letter2.num
148     // and both must have the same length
149     if( !EQ( ElmPRec(letter1, RNnum), ElmPRec(letter2, RNnum) ) ||
150     !EQ( ElmPRec(letter1, RNlength), ElmPRec(letter2, RNlength) ) )
151     {
152         return False;
153     }
154 
155     // check whether both are atoms / non-atoms:
156     UInt is_atom = IsbPRec(letter1, RNside); // 1 <=> letter1 is atom
157     if( is_atom == IsbPRec(letter2, RNside) ){
158         if( is_atom == 1) // both are atoms
159         {
160             if( EQ( ElmPRec(letter1, RNside), ElmPRec(letter2, RNside) ) ){
161                 return True;
162             } else {
163                 return False;
164             }
165 
166         } else { // both are non-atoms
167             if( EQ( ElmPRec(letter1, RNleft), ElmPRec(letter2, RNleft) ) &&
168                 EQ( ElmPRec(letter1, RNright), ElmPRec(letter2, RNright) ) ){
169                 return True;
170             } else {
171                 return False;
172             }
173         }
174     } else { // one is an atom, the other one a non-atom
175         return False;
176     }
177 }
178 
179 
180 typedef Obj (* GVarFunc)(/*arguments*/);
181 
182 #define GVAR_FUNC_TABLE_ENTRY(srcfile, name, nparam, params) \
183   {#name, nparam, \
184    params, \
185    (GVarFunc)name, \
186    srcfile ":Func" #name }
187 
188 // Table of functions to export
189 static StructGVarFunc GVarFuncs [] = {
190     GVAR_FUNC_TABLE_ENTRY("DeepThought.c", DTP_Binomial, 2, "n, k"),
191     GVAR_FUNC_TABLE_ENTRY("DeepThought.c", DTP_SequenceLetter, 2, "letter, seq"),
192     GVAR_FUNC_TABLE_ENTRY("DeepThought.c", DTP_Seq_i, 2, "letter, i"),
193     GVAR_FUNC_TABLE_ENTRY("DeepThought.c", DTP_AreAlmostEqual, 2, "letter1, letter2"),
194     { 0 } /* Finish with an empty entry */
195 
196 };
197 
198 /******************************************************************************
199 *F  InitKernel( <module> )  . . . . . . . . initialise kernel data structures
200 */
InitKernel(StructInitInfo * module)201 static Int InitKernel( StructInitInfo *module )
202 {
203     /* init filters and functions                                          */
204     InitHdlrFuncsFromTable( GVarFuncs );
205 
206     /* return success                                                      */
207     return 0;
208 }
209 
PostRestore(StructInitInfo * module)210 static Int PostRestore( StructInitInfo *module )
211 {
212     RNleft = RNamName("left");
213     RNright = RNamName("right");
214     RNlength = RNamName("l");
215     RNnum = RNamName("num");
216     RNside = RNamName("side");
217 
218     return 0;
219 }
220 
221 /******************************************************************************
222 *F  InitLibrary( <module> ) . . . . . . .  initialise library data structures
223 */
InitLibrary(StructInitInfo * module)224 static Int InitLibrary( StructInitInfo *module )
225 {
226     /* init filters and functions */
227     InitGVarFuncsFromTable( GVarFuncs );
228 
229     PostRestore(module);
230 
231     /* return success                                                      */
232     return 0;
233 }
234 
235 
236 
237 /******************************************************************************
238 *F  InitInfopl()  . . . . . . . . . . . . . . . . . table of init functions
239 */
240 static StructInitInfo module = {
241  /* type        = */ MODULE_DYNAMIC,
242  /* name        = */ "DeepThought",
243  /* revision_c  = */ 0,
244  /* revision_h  = */ 0,
245  /* version     = */ 0,
246  /* crc         = */ 0,
247  /* initKernel  = */ InitKernel,
248  /* initLibrary = */ InitLibrary,
249  /* checkInit   = */ 0,
250  /* preSave     = */ 0,
251  /* postSave    = */ 0,
252  /* postRestore = */ PostRestore
253 };
254 
Init__Dynamic(void)255 StructInitInfo *Init__Dynamic( void )
256 {
257     return &module;
258 }
259