Branch data Line data Source code
1 : : /*
2 : : * SPDX-License-Identifier: MIT
3 : : * Copyright (c) 2024 IUPAC and InChI Trust
4 : : */
5 : :
6 : : #include <string.h>
7 : :
8 : : #include "mode.h"
9 : : #include "ichicomn.h"
10 : : #include "ichicant.h"
11 : :
12 : : #include "bcf_s.h"
13 : :
14 : : #if 0
15 : : #define RET_MAX 32767
16 : : #define CUTOFF 8 /* testing shows that this is good value */
17 : : #endif
18 : :
19 : : /* Note: the theoretical number of stack entries required is
20 : : no more than 1 + log2(num). But we switch to insertion
21 : : sort for CUTOFF elements or less, so we really only need
22 : : 1 + log2(num) - log2(CUTOFF) stack entries. For a CUTOFF
23 : : of 8, that means we need no more than 30 stack entries for
24 : : 32 bit platforms, and 62 for 64-bit platforms. */
25 : : #define STKSIZ (8*sizeof(void*) - 2)
26 : :
27 : :
28 : : /****************************************************************************
29 : : inchi's qsort
30 : : ****************************************************************************/
31 : 749 : void inchi_qsort( void *pParam,
32 : : void *base,
33 : : size_t num,
34 : : size_t width,
35 : : int( *comp )( const void *, const void *, void * ) )
36 : : {
37 : : char *lo, *hi; /* ends of sub-array currently sorting */
38 : : char *mid; /* points to middle of subarray */
39 : : char *loguy, *higuy; /* traveling pointers for partition step */
40 : : size_t size; /* size of the sub-array */
41 : : char *lostk[STKSIZ], *histk[STKSIZ];
42 : : int stkptr; /* stack for saving sub-array to be processed */
43 : :
44 [ + + ]: 749 : if (num < 2)
45 : 354 : return; /* nothing to do */
46 : :
47 : 395 : stkptr = 0; /* initialize stack */
48 : :
49 : 395 : lo = (char *) base;
50 : 395 : hi = (char *) base + width * ( num - 1 ); /* initialize limits */
51 : :
52 : : /* This entry point is for pseudo-recursion calling: setting
53 : : lo and hi and jumping to here is like recursion, but stkptr is
54 : : preserved, locals aren't, so we preserve stuff on the stack */
55 : 1623 : recurse:
56 : :
57 : 2018 : size = ( hi - lo ) / width + 1; /* number of el's to sort */
58 : :
59 : : /* First we pick a partitioning element. The efficiency of the
60 : : algorithm demands that we find one that is approximately the median
61 : : of the values, but also that we select one fast. We choose the
62 : : median of the first, middle, and last elements, to avoid bad
63 : : performance in the face of already sorted data, or data that is made
64 : : up of multiple sorted runs appended together. Testing shows that a
65 : : median-of-three algorithm provides better performance than simply
66 : : picking the middle element for the latter case. */
67 : :
68 : 2018 : mid = lo + ( size / 2 ) * width; /* find middle element */
69 : :
70 : : /* Sort the first, middle, last elements into order */
71 [ + + ]: 2018 : if (comp( lo, mid, pParam ) > 0)
72 : : {
73 : 191 : inchi_swap( lo, mid, width );
74 : : }
75 [ + + ]: 2018 : if (comp( lo, hi, pParam ) > 0)
76 : : {
77 : 97 : inchi_swap( lo, hi, width );
78 : : }
79 [ + + ]: 2018 : if (comp( mid, hi, pParam ) > 0)
80 : : {
81 : 248 : inchi_swap( mid, hi, width );
82 : : }
83 : :
84 : : /* We now wish to partition the array into three pieces, one consisting
85 : : of elements <= partition element, one of elements equal to the
86 : : partition element, and one of elements > than it. This is done
87 : : below; comments indicate conditions established at every step. */
88 : :
89 : 2018 : loguy = lo;
90 : 2018 : higuy = hi;
91 : :
92 : : /* Note that higuy decreases and loguy increases on every iteration,
93 : : so loop must terminate. */
94 : : for (;;)
95 : : {
96 : : /* lo <= loguy < hi, lo < higuy <= hi,
97 : : A[i] <= A[mid] for lo <= i <= loguy,
98 : : A[i] > A[mid] for higuy <= i < hi,
99 : : A[hi] >= A[mid] */
100 : :
101 : : /* The doubled loop is to avoid calling comp(mid,mid), since some
102 : : existing comparison funcs don't work when passed the same
103 : : value for both pointers. */
104 : :
105 [ + + ]: 2265 : if (mid > loguy)
106 : : {
107 : : do
108 : : {
109 : 4247 : loguy += width;
110 : : }
111 [ + + + + ]: 4247 : while (loguy < mid && comp( loguy, mid, pParam ) <= 0);
112 : : }
113 [ + + ]: 2265 : if (mid <= loguy)
114 : : {
115 : : do
116 : : {
117 : 2216 : loguy += width;
118 : : }
119 [ + + + + ]: 2216 : while (loguy <= hi && comp( loguy, mid, pParam ) <= 0);
120 : : }
121 : :
122 : : /* lo < loguy <= hi+1, A[i] <= A[mid] for lo <= i < loguy,
123 : : either loguy > hi or A[loguy] > A[mid] */
124 : :
125 : : do
126 : : {
127 : 3906 : higuy -= width;
128 : : }
129 [ + + + + ]: 3906 : while (higuy > mid && comp( higuy, mid, pParam ) > 0);
130 : :
131 : : /* lo <= higuy < hi, A[i] > A[mid] for higuy < i < hi,
132 : : either higuy == lo or A[higuy] <= A[mid] */
133 : :
134 [ + + ]: 2265 : if (higuy < loguy)
135 : : {
136 : 2018 : break;
137 : : }
138 : :
139 : : /* if loguy > hi or higuy == lo, then we would have exited, so
140 : : A[loguy] > A[mid], A[higuy] <= A[mid],
141 : : loguy <= hi, higuy > lo */
142 : :
143 : 247 : inchi_swap( loguy, higuy, width );
144 : :
145 : : /* If the partition element was moved, follow it. Only need
146 : : to check for mid == higuy, since before the swap,
147 : : A[loguy] > A[mid] implies loguy != mid. */
148 : :
149 [ + + ]: 247 : if (mid == higuy)
150 : : {
151 : 119 : mid = loguy;
152 : : }
153 : :
154 : : /* A[loguy] <= A[mid], A[higuy] > A[mid]; so condition at top
155 : : of loop is re-established */
156 : : }
157 : :
158 : : /* A[i] <= A[mid] for lo <= i < loguy,
159 : : A[i] > A[mid] for higuy < i < hi,
160 : : A[hi] >= A[mid]
161 : : higuy < loguy
162 : : implying:
163 : : higuy == loguy-1
164 : : or higuy == hi - 1, loguy == hi + 1, A[hi] == A[mid] */
165 : :
166 : : /* Find adjacent elements equal to the partition element. The
167 : : doubled loop is to avoid calling comp(mid,mid), since some
168 : : existing comparison funcs don't work when passed the same value
169 : : for both pointers. */
170 : :
171 : 2018 : higuy += width;
172 [ + + ]: 2018 : if (mid < higuy)
173 : : {
174 : : do
175 : : {
176 : 1339 : higuy -= width;
177 : : }
178 [ + + - + ]: 1339 : while (higuy > mid && comp( higuy, mid, pParam ) == 0);
179 : : }
180 [ + + ]: 2018 : if (mid >= higuy)
181 : : {
182 : : do
183 : : {
184 : 1908 : higuy -= width;
185 : : }
186 [ + + - + ]: 1908 : while (higuy > lo && comp( higuy, mid, pParam ) == 0);
187 : : }
188 : :
189 : : /* OK, now we have the following:
190 : : higuy < loguy
191 : : lo <= higuy <= hi
192 : : A[i] <= A[mid] for lo <= i <= higuy
193 : : A[i] == A[mid] for higuy < i < loguy
194 : : A[i] > A[mid] for loguy <= i < hi
195 : : A[hi] >= A[mid] */
196 : :
197 : : /* We've finished the partition, now we want to sort the subarrays
198 : : [lo, higuy] and [loguy, hi].
199 : : We do the smaller one first to minimize stack usage.
200 : : We only sort arrays of length 2 or more.*/
201 : :
202 [ + + ]: 2018 : if (higuy - lo >= hi - loguy)
203 : : {
204 [ + + ]: 1906 : if (lo < higuy)
205 : : {
206 : 950 : lostk[stkptr] = lo;
207 : 950 : histk[stkptr] = higuy;
208 : 950 : ++stkptr;
209 : : } /* save big recursion for later */
210 : :
211 [ + + ]: 1906 : if (loguy < hi)
212 : : {
213 : 521 : lo = loguy;
214 : 521 : goto recurse; /* do small recursion */
215 : : }
216 : : }
217 : : else
218 : : {
219 [ + - ]: 112 : if (loguy < hi)
220 : : {
221 : 112 : lostk[stkptr] = loguy;
222 : 112 : histk[stkptr] = hi;
223 : 112 : ++stkptr; /* save big recursion for later */
224 : : }
225 : :
226 [ + + ]: 112 : if (lo < higuy)
227 : : {
228 : 40 : hi = higuy;
229 : 40 : goto recurse; /* do small recursion */
230 : : }
231 : : }
232 : :
233 : : /* We have sorted the array, except for any pending sorts on the stack.
234 : : Check if there are any, and do them. */
235 : :
236 : 1457 : --stkptr;
237 [ + + ]: 1457 : if (stkptr >= 0)
238 : : {
239 : 1062 : lo = lostk[stkptr];
240 : 1062 : hi = histk[stkptr];
241 : 1062 : goto recurse; /* pop subarray from stack */
242 : : }
243 : : else
244 : : {
245 : 395 : return; /* all subarrays done */
246 : : }
247 : : }
248 : :
249 : :
250 : : /****************************************************************************/
251 : 4322 : void inchi_swap( char *a, char *b, size_t width )
252 : : {
253 : : char tmp;
254 [ + - ]: 4322 : if (a != b)
255 : : {
256 [ + + ]: 15248 : while (width--)
257 : : {
258 : 10926 : tmp = *a;
259 : 10926 : *a++ = *b;
260 : 10926 : *b++ = tmp;
261 : : }
262 : : }
263 : 4322 : }
264 : :
265 : :
266 : : /****************************************************************************
267 : : Sort by insertions
268 : : ****************************************************************************/
269 : 4089 : int insertions_sort( void *pCG,
270 : : void *base,
271 : : size_t num, size_t width,
272 : : int( *compare )( const void *, const void *, void * ) ) /* djb-rwth: types of variables are sufficient */
273 : : {
274 : 4089 : char *i, *j, *pk = (char*) base;
275 : 4089 : int num_trans = 0;
276 : : size_t k;
277 [ + + ]: 12795 : for (k = 1; k < num; k++, pk += width)
278 : : {
279 : : /*for( i = pk, j = pk + width; j > (char*)base && (*compare)(i,j) > 0; j=i, i -= width )*/
280 : 8706 : for (i = j = pk + width;
281 [ + + + + ]: 11731 : j > ( char* )base && ( i -= width, ( *compare )( i, j, pCG ) ) > 0;
282 : 3025 : j = i) /* changed to keep BoundsChecker happy 2007-09-24 DT */
283 : : {
284 : 3025 : inchi_swap( i, j, width );
285 : 3025 : num_trans++;
286 : : }
287 : : }
288 : :
289 : 4089 : return num_trans;
290 : : }
291 : :
292 : :
293 : : /****************************************************************************
294 : : Sort by insertions
295 : : ****************************************************************************/
296 : 593 : int insertions_sort_AT_NUMBERS( void *pCG,
297 : : AT_NUMB *base,
298 : : int num,
299 : : int( *compare )( const void *e1, const void *e2, void * ) )
300 : : {
301 : : AT_NUMB *i, *j, *pk, tmp;
302 : 593 : int k, num_trans = 0;
303 [ + + ]: 1354 : for (k = 1, pk = base; k < num; k++, pk++)
304 : : {
305 [ + + + + ]: 865 : for (j = ( i = pk ) + 1, tmp = *j; j > base && ( *compare )( i, &tmp, pCG ) > 0; j = i, i--)
306 : : {
307 : 104 : *j = *i;
308 : 104 : num_trans++;
309 : : }
310 : 761 : *j = tmp;
311 : : }
312 : :
313 : 593 : return num_trans;
314 : : }
315 : :
316 : :
317 : : /****************************************************************************
318 : : Sort neighbors according to ranks in ascending order
319 : : ****************************************************************************/
320 : 2176 : void insertions_sort_NeighList_AT_NUMBERS( NEIGH_LIST base, AT_RANK *nRank )
321 : : {
322 : : AT_NUMB *i, *j, *pk, tmp;
323 : : AT_RANK rj; /* optimization */
324 : 2176 : int k, num = (int) *base++;
325 [ + + ]: 5723 : for (k = 1, pk = base; k < num; k++, pk++)
326 : : {
327 [ + + + + ]: 3984 : for (j = ( i = pk ) + 1, rj = nRank[(int) *j]; j > base && nRank[(int) *i] > rj; j = i, i--)
328 : : {
329 : 437 : tmp = *i;
330 : 437 : *i = *j;
331 : 437 : *j = tmp;
332 : : }
333 : : }
334 : 2176 : }
335 : :
336 : :
337 : : /****************************************************************************
338 : : Sort neighbors according to ranks in ascending order
339 : : ****************************************************************************/
340 : 0 : int insertions_sort_AT_RANK( AT_RANK *base, int num )
341 : : {
342 : : AT_RANK *i, *j, *pk, tmp;
343 : 0 : int k, num_trans = 0;
344 [ # # ]: 0 : for (k = 1, pk = base; k < num; k++, pk++)
345 : : {
346 [ # # # # ]: 0 : for (j = ( i = pk ) + 1, tmp = *j; j > base && *i > tmp; j = i, i--)
347 : : {
348 : 0 : *j = *i;
349 : 0 : num_trans++;
350 : : }
351 : 0 : *j = tmp;
352 : : }
353 : :
354 : 0 : return num_trans;
355 : : }
356 : :
357 : :
358 : : /****************************************************************************
359 : : Sort neighbors according to ranks in ascending order
360 : : ****************************************************************************/
361 : 924 : int insertions_sort_NeighList_AT_NUMBERS3( NEIGH_LIST base, AT_RANK *nRank )
362 : : {
363 : : AT_NUMB *i, *j, *pk, tmp;
364 : : AT_RANK rj;
365 : 924 : int k, n, num = (int) *base++;
366 [ + + ]: 2128 : for (k = 1, pk = base, n = 0; k < num; k++, pk++)
367 : : {
368 : 1204 : for (j = ( i = pk ) + 1, rj = nRank[(int) ( tmp = *j )];
369 [ + + + + ]: 1219 : j > base && nRank[(int) *i] > rj;
370 : 15 : j = i, i--)
371 : : {
372 : 15 : *j = *i;
373 : 15 : n++;
374 : : }
375 : 1204 : *j = tmp;
376 : : }
377 : :
378 : 924 : return n;
379 : : }
380 : :
381 : :
382 : : /****************************************************************************
383 : : Sort neighbors according to symm. ranks (primary key) and canon.
384 : : ranks (secondary key), in descending order
385 : : ****************************************************************************/
386 : 3318 : void insertions_sort_NeighListBySymmAndCanonRank( NEIGH_LIST base,
387 : : const AT_RANK *nSymmRank,
388 : : const AT_RANK *nCanonRank )
389 : : {
390 : : AT_NUMB *i, *j, *pk, tmp;
391 : : int diff;
392 : 3318 : int k, num = (int) *base++;
393 [ + + ]: 6120 : for (k = 1, pk = base; k < num; k++, pk++)
394 : : {
395 : 2802 : for (j = ( i = pk ) + 1;
396 [ + + ]: 4959 : j > base && /* always j > i */
397 [ + + + + ]: 3702 : ( 0 > ( diff = (int) nSymmRank[(int) *i] - (int) nSymmRank[(int) *j] ) ||
398 [ + - ]: 96 : (!diff && nCanonRank[(int) *i] < nCanonRank[(int) *j]) ); /* djb-rwth: addressing LLVM warning */
399 : 2157 : j = i, i--)
400 : : {
401 : 2157 : tmp = *i;
402 : 2157 : *i = *j;
403 : 2157 : *j = tmp;
404 : : }
405 : : }
406 : :
407 : 3318 : }
408 : :
409 : :
410 : : /****************************************************************************
411 : : *
412 : : * Comparison functions
413 : : *
414 : : ****************************************************************************/
415 : :
416 : :
417 : : /****************************************************************************/
418 : 2091 : int CompNeighborsAT_NUMBER( const void* a1, const void* a2, void *p )
419 : : {
420 : 2091 : CANON_GLOBALS *pCG = (CANON_GLOBALS *) p;
421 : : #ifdef CT_NEIGH_INCREASE
422 : 4182 : return (int) pCG->m_pn_RankForSort[pCG->m_pNeighborsForSort[( int )*(const AT_NUMB*) a1]] -
423 : 2091 : (int) pCG->m_pn_RankForSort[pCG->m_pNeighborsForSort[( int )*(const AT_NUMB*) a2]];
424 : : #else
425 : : return (int) ( (CANON_GLOBALS *) pCG )->m_pn_RankForSort[pNeighborsForSort[( int )*(const AT_NUMB*) a2]] -
426 : : (int) ( (CANON_GLOBALS *) pCG )->m_pn_RankForSort[pNeighborsForSort[( int )*(const AT_NUMB*) a1]];
427 : : #endif
428 : : }
429 : :
430 : :
431 : : /**********************************************************************************/
432 : 741 : int comp_AT_RANK( const void* a1, const void* a2, void *p )
433 : : {
434 : 741 : return ( int )*(const AT_RANK*) a1 - ( int )*(const AT_RANK*) a2;
435 : : }
436 : :
437 : :
438 : : /**********************************************************************************/
439 : : /* Compare for sorting Ranks only */
440 : 1964 : int CompRank( const void* a1, const void* a2, void *p )
441 : : {
442 : 1964 : CANON_GLOBALS *pCG = (CANON_GLOBALS *) p;
443 : 1964 : int ret = (int) pCG->m_pn_RankForSort[( int )*(const AT_RANK*) a1] -
444 : 1964 : (int) pCG->m_pn_RankForSort[( int )*(const AT_RANK*) a2];
445 : :
446 : 1964 : return ret;
447 : : }
448 : :
449 : :
450 : : /**********************************************************************************/
451 : 9426 : int CompRanksOrd( const void* a1, const void* a2, void *p )
452 : : {
453 : : int ret;
454 : 9426 : CANON_GLOBALS *pCG = (CANON_GLOBALS *) p;
455 : 9426 : ret = (int) pCG->m_pn_RankForSort[( int )*(const AT_RANK*) a1] -
456 : 9426 : (int) pCG->m_pn_RankForSort[( int )*(const AT_RANK*) a2];
457 [ + + ]: 9426 : if (!ret)
458 : : {
459 : 2292 : ret = ( int )*(const AT_RANK*) a1 - ( int )*(const AT_RANK*) a2;
460 : : }
461 : :
462 : 9426 : return ret;
463 : : }
464 : :
465 : :
466 : : /**********************************************************************************/
467 : 4145 : int CompAtomInvariants2Only( const void* a1, const void* a2, void *p )
468 : : {
469 : 4145 : CANON_GLOBALS *pCG = (CANON_GLOBALS *) p;
470 : 4145 : const ATOM_INVARIANT2 *pAI1 = pCG->m_pAtomInvariant2ForSort + ( int )*(const AT_RANK*) a1;
471 : 4145 : const ATOM_INVARIANT2 *pAI2 = pCG->m_pAtomInvariant2ForSort + ( int )*(const AT_RANK*) a2;
472 : : int i;
473 [ + + ]: 20633 : for (i = 0; i < AT_INV_BREAK1; i++)
474 : : {
475 [ + + ]: 18374 : if (pAI1->val[i] == pAI2->val[i])
476 : 16488 : continue;
477 : 1886 : return (int) pAI1->val[i] - (int) pAI2->val[i];
478 : : }
479 [ - + ]: 2259 : if (pAI1->iso_sort_key != pAI2->iso_sort_key)
480 : : {
481 [ # # ]: 0 : return ( pAI1->iso_sort_key > pAI2->iso_sort_key ) ? 1 : -1;
482 : : }
483 [ + - ]: 2259 : for (; i < AT_INV_LENGTH; i++)
484 : : {
485 [ - + ]: 2259 : if (pAI1->val[i] != pAI2->val[i])
486 : : {
487 : 0 : continue;
488 : : }
489 : 2259 : return (int) pAI1->val[i] - (int) pAI2->val[i];
490 : : }
491 [ # # ]: 0 : if (pAI1->iso_aux_key != pAI2->iso_aux_key)
492 : : {
493 [ # # ]: 0 : return ( pAI1->iso_aux_key > pAI2->iso_aux_key ) ? 1 : -1;
494 : : }
495 : :
496 : 0 : return 0;
497 : : }
498 : :
499 : :
500 : : /**********************************************************************************/
501 : 3398 : int CompAtomInvariants2( const void* a1, const void* a2, void *p )
502 : : {
503 : : /* Warning: the following line may be compiler implementation dependent */
504 : 3398 : int ret = CompAtomInvariants2Only( a1, a2, p );
505 [ + + ]: 3398 : if (!ret)
506 : : {
507 : 1755 : ret = ( int )*(const AT_RANK*) a1 - ( int )*(const AT_RANK*) a2;
508 : : }
509 : :
510 : 3398 : return ret;
511 : : }
512 : :
513 : :
514 : : /**********************************************************************************/
515 : : /* Compare two elements lexicographically */
516 : 50 : int CompChemElemLex( const void *a1, const void *a2 )
517 : : {
518 : 50 : return memcmp( a1, a2, 2 );
519 : : }
520 : :
521 : :
522 : : /****************************************************************************
523 : : Lexicographic compare
524 : : ****************************************************************************/
525 : 3974 : int CompareNeighListLex( NEIGH_LIST pp1, NEIGH_LIST pp2, const AT_RANK *nRank )
526 : : {
527 : 3974 : int len1 = (int) *pp1++;
528 : 3974 : int len2 = (int) *pp2++;
529 : 3974 : int len = inchi_min( len1, len2 );
530 : 3974 : int diff = 0;
531 : : int ret;
532 : : /* djb-rwth: fixing oss-fuzz issue #25642 */
533 [ + + + + ]: 10000 : while ((len > 0) && !diff)
534 : : {
535 : 6026 : len--;
536 : 6026 : diff = (int)nRank[*pp1++] - (int)nRank[*pp2++];
537 : : };
538 : :
539 [ + + ]: 3974 : ret = diff ? diff : (len1 - len2);
540 : 3974 : return ret;
541 : : }
542 : :
543 : :
544 : : /****************************************************************************
545 : : Lexicographic compare
546 : : ****************************************************************************/
547 : 578 : int CompareNeighListLexUpToMaxRank( NEIGH_LIST pp1, NEIGH_LIST pp2, const AT_RANK *nRank, AT_RANK nMaxAtNeighRank )
548 : : {
549 : 578 : int len1 = (int) *pp1++;
550 : 578 : int len2 = (int) *pp2++;
551 : 578 : int diff = 0;
552 : : int len;
553 [ + - + + ]: 1019 : while (0 < len1 && nRank[pp1[len1 - 1]] > nMaxAtNeighRank)
554 : : {
555 : 441 : len1--;
556 : : }
557 [ + - + + ]: 1019 : while (0 < len2 && nRank[pp2[len2 - 1]] > nMaxAtNeighRank)
558 : : {
559 : 441 : len2--;
560 : : }
561 : 578 : len = inchi_min( len1, len2 );
562 [ + + + + ]: 1244 : while (len-- > 0 && !( diff = (int) nRank[*pp1++] - (int) nRank[*pp2++] ))
563 : : {
564 : : ;
565 : : }
566 : :
567 [ + + ]: 578 : return diff ? diff : ( len1 - len2 );
568 : : }
569 : :
570 : :
571 : : /****************************************************************************/
572 : 3492 : int compare_NeighLists( const NEIGH_LIST *op1, const NEIGH_LIST *op2, void *p )
573 : : {
574 : 3492 : CANON_GLOBALS *pCG = (CANON_GLOBALS *) p;
575 : 3492 : return CompareNeighListLex( *op1, *op2, pCG->m_pn_RankForSort );
576 : : }
577 : :
578 : :
579 : : /****************************************************************************/
580 : 7698 : int CompNeighListRanks( const void* a1, const void* a2, void *p )
581 : : {
582 : 7698 : CANON_GLOBALS *pCG = (CANON_GLOBALS *) p;
583 : : int ret;
584 : 7698 : ret = (int) pCG->m_pn_RankForSort[*( (const AT_RANK*) a1 )] -
585 : 7698 : (int) pCG->m_pn_RankForSort[*( (const AT_RANK*) a2 )];
586 [ + + ]: 7698 : if (!ret)
587 : : {
588 : 2985 : ret = compare_NeighLists( pCG->m_pNeighList_RankForSort + *( (const AT_RANK*) a1 ),
589 : 2985 : pCG->m_pNeighList_RankForSort + *( (const AT_RANK*) a2 ), p );
590 : : }
591 : :
592 : 7698 : return ret;
593 : : }
594 : :
595 : :
596 : : /****************************************************************************/
597 : 507 : int CompNeighLists( const void* a1, const void* a2, void *p )
598 : : {
599 : 507 : CANON_GLOBALS *pCG = (CANON_GLOBALS *) p;
600 : : int ret;
601 : 507 : ret = compare_NeighLists( pCG->m_pNeighList_RankForSort + *( (const AT_RANK*) a1 ),
602 : 507 : pCG->m_pNeighList_RankForSort + *( (const AT_RANK*) a2 ), p );
603 : :
604 : 507 : return ret;
605 : : }
606 : :
607 : :
608 : : /****************************************************************************/
609 : 299 : int CompNeighListsUpToMaxRank( const void* a1, const void* a2, void *p )
610 : : {
611 : 299 : CANON_GLOBALS *pCG = (CANON_GLOBALS *) p;
612 : : int ret;
613 : 299 : ret = CompareNeighListLexUpToMaxRank( pCG->m_pNeighList_RankForSort[*( (const AT_RANK*) a1 )],
614 : 299 : pCG->m_pNeighList_RankForSort[*( (const AT_RANK*) a2 )],
615 : 299 : pCG->m_pn_RankForSort, pCG->m_nMaxAtNeighRankForSort );
616 : :
617 : 299 : return ret;
618 : : }
619 : :
620 : :
621 : : /****************************************************************************/
622 : 4219 : int CompNeighListRanksOrd( const void* a1, const void* a2, void *p )
623 : : {
624 : 4219 : int ret = CompNeighListRanks( a1, a2, p );
625 [ + + ]: 4219 : if (!ret)
626 : : {
627 : 811 : ret = ( int )*( (const AT_RANK*) a1 ) - ( int )*( (const AT_RANK*) a2 ); /* keep original order if identical */
628 : : }
629 : :
630 : 4219 : return ret;
631 : : }
632 : :
633 : :
634 : : /****************************************************************************/
635 : 0 : int CompRanksInvOrd( const void* a1, const void* a2, void *p )
636 : : {
637 : 0 : return ( int )*(const AT_RANK*) a2 - ( int )*(const AT_RANK*) a1;
638 : : }
639 : :
640 : :
641 : : /****************************************************************************/
642 : 0 : int CompNeighborsRanksCountEql( const void* a1, const void* a2, void *p )
643 : : {
644 : 0 : CANON_GLOBALS *pCG = (CANON_GLOBALS *) p;
645 : : #ifdef CT_NEIGH_INCREASE
646 : 0 : int ret = (int) pCG->m_pn_RankForSort[( int )*(const AT_RANK*) a1] -
647 : 0 : (int) pCG->m_pn_RankForSort[( int )*(const AT_RANK*) a2];
648 : : #else
649 : : int ret = (int) pCG->m_pn_RankForSort[( int )*(const AT_RANK*) a2] -
650 : : (int) pCG->m_pn_RankForSort[( int )*(const AT_RANK*) a1];
651 : : #endif
652 : 0 : pCG->m_nNumCompNeighborsRanksCountEql += !ret;
653 : :
654 : 0 : return ret;
655 : : }
656 : :
657 : :
658 : : /****************************************************************************
659 : : *
660 : : * In this neighbor list the (vertex number) = (canonical number) - 1
661 : : * Since LinearCT is sorted so that parents are in ascending order
662 : : * and all neighbors of a parent are smaller than the parent and are
663 : : * in ascending order, the neighbors in the NEIGH_LIST are automatically
664 : : * sorted in ascending order
665 : : ****************************************************************************/
666 : 70 : NEIGH_LIST *CreateNeighListFromLinearCT( AT_NUMB *LinearCT, int nLenCT, int num_atoms )
667 : : {
668 : : /* Atom numbers in LinearCT are canonical numbers
669 : : * order: parent[i] > neigh[i][0] < neigh[i][1]...<neigh[i][n] < parent[i+1] > neigh[i+1][0] < ...
670 : : * parent[i] < parent[i+1]
671 : : */
672 : : int i, j;
673 : 70 : S_CHAR *valence = NULL;
674 : 70 : NEIGH_LIST *pp = NULL;
675 : 70 : AT_NUMB *pAtList = NULL;
676 : : AT_RANK n_vertex, n_neigh;
677 : 70 : int err = 1, num_bonds;
678 : : int length, start;
679 : :
680 [ - + ]: 70 : if ((int) LinearCT[0] > num_atoms)
681 : : {
682 : 0 : goto exit_function;
683 : : }
684 : 70 : valence = (S_CHAR*)inchi_calloc((long long)num_atoms + 1, sizeof(valence[0]));
685 [ - + ]: 70 : if (!valence) /* djb-rwth: cast operator added */
686 : : {
687 : 0 : goto exit_function;
688 : : }
689 : :
690 [ + + ]: 1219 : for (i = 1, num_bonds = 0, n_vertex = LinearCT[0]; i < nLenCT; i++)
691 : : {
692 [ + + ]: 1149 : if (( n_neigh = LinearCT[i] ) < n_vertex)
693 : : {
694 : 587 : valence[n_neigh] ++;
695 : 587 : valence[n_vertex] ++;
696 : 587 : num_bonds += 2;
697 : : }
698 : : else
699 : : {
700 [ - + ]: 562 : if ((int) ( n_vertex = n_neigh ) > num_atoms)
701 : : {
702 : 0 : goto exit_function;
703 : : }
704 : : }
705 : : }
706 [ - + ]: 70 : if ((int) n_vertex != num_atoms)
707 : : {
708 : 0 : goto exit_function;
709 : : }
710 : 70 : length = num_bonds + num_atoms + 1;
711 : 70 : pp = (NEIGH_LIST*)inchi_calloc(((long long)num_atoms + 1), sizeof(NEIGH_LIST));
712 : 70 : pAtList = (AT_NUMB*)inchi_malloc(length * sizeof(AT_NUMB));
713 [ - + - + ]: 70 : if (pp && pAtList) /* djb-rwth: cast operator added; addressing LLVM warning */
714 : : {
715 : : /* Create empty connection table */
716 [ + + ]: 702 : for (i = 1, length = 0; i <= num_atoms; i++)
717 : : {
718 : 632 : start = length;
719 : 632 : length += ( valence[i] + 1 );
720 : 632 : pp[i - 1] = pAtList + start;
721 : 632 : pp[i - 1][0] = 0;
722 : : }
723 : : /* Fill out the CT */
724 [ + + ]: 1219 : for (i = 1, n_vertex = LinearCT[0] - 1; i < nLenCT; i++)
725 : : {
726 [ + + ]: 1149 : if (( n_neigh = LinearCT[i] - 1 ) < n_vertex)
727 : : {
728 : : /* Vertex - neighbor connection */
729 : 587 : j = (int) ( ++pp[(int) n_vertex][0] );
730 : 587 : pp[(int) n_vertex][j] = n_neigh;
731 : : /* neighbor - vertex connection */
732 : 587 : j = (int) ( ++pp[(int) n_neigh][0] );
733 : 587 : pp[(int) n_neigh][j] = n_vertex;
734 : : }
735 : : else
736 : : {
737 [ - + ]: 562 : if ((int) ( n_vertex = n_neigh ) >= num_atoms)
738 : : {
739 : 0 : goto exit_function;
740 : : }
741 : : }
742 : : }
743 : 70 : err = 0;
744 : : }
745 : :
746 : 0 : exit_function:
747 [ + - ]: 70 : if (valence)
748 : : {
749 [ + - ]: 70 : inchi_free( valence );
750 : : }
751 [ - + ]: 70 : if (err) /* djb-rwth: ignoring LLVM warning */
752 : : {
753 [ # # ]: 0 : if (pAtList)
754 : : {
755 [ # # ]: 0 : inchi_free( pAtList );
756 : : }
757 [ # # ]: 0 : if (pp)
758 : : {
759 [ # # ]: 0 : inchi_free( pp );
760 : 0 : pp = NULL;
761 : : }
762 : : }
763 : :
764 : 70 : return pp; /* djb-rwth: ignoring LLVM warning: since a pointer is returned, memory should be freed in a function which calls *CreateNeighListFromLinearCT */
765 : : }
766 : :
767 : :
768 : : /****************************************************************************
769 : : NEIGH_LIST pp[] is an array of pointers to the lists of neighboring
770 : : atoms numbers. The first number in each list is a number of neighbors.
771 : : In case of bDoubleBondSquare != 0 neighbors connected by the double bond
772 : : appear 2 times. The first element pp[0] is a pointer to be deallocated
773 : : to free all the lists.
774 : : ****************************************************************************/
775 : 502 : NEIGH_LIST *CreateNeighList( int num_atoms,
776 : : int num_at_tg,
777 : : sp_ATOM* at,
778 : : int bDoubleBondSquare,
779 : : T_GROUP_INFO *t_group_info )
780 : : {
781 : : /* +1 to add NULL termination */
782 : 502 : NEIGH_LIST *pp = (NEIGH_LIST *) inchi_calloc( ( (long long)num_at_tg + 1 ), sizeof( NEIGH_LIST ) ); /* djb-rwth: cast operator added */
783 : 502 : T_GROUP *t_group = NULL;
784 : 502 : AT_NUMB *nEndpointAtomNumber = NULL;
785 : 502 : int num_t_groups = 0;
786 : : int nFirstEndpointAtNoPos;
787 : :
788 : 502 : AT_NUMB *pAtList = NULL;
789 : : int length, start, val, i, j;
790 [ + - ]: 502 : if (pp)
791 : : {
792 [ + + ]: 502 : if (num_at_tg > num_atoms)
793 : : {
794 : 7 : t_group = t_group_info->t_group;
795 : 7 : num_t_groups = t_group_info->num_t_groups;
796 : 7 : nEndpointAtomNumber = t_group_info->nEndpointAtomNumber;
797 : : }
798 : :
799 [ + - ]: 502 : if (!bDoubleBondSquare)
800 : : {
801 [ + + ]: 4681 : for (i = 0, length = 0; i < num_atoms; i++)
802 : : {
803 [ + + + + ]: 4179 : length += (int) at[i].valence + ( num_t_groups && at[i].endpoint );
804 : : }
805 : 502 : length += num_atoms;
806 [ + + ]: 514 : for (i = 0; i < num_t_groups; i++)
807 : : {
808 : 12 : length += (int) t_group[i].nNumEndpoints;
809 : : }
810 : 502 : length += num_t_groups;
811 : : }
812 : : else
813 : : {
814 [ # # ]: 0 : for (i = 0, length = 0; i < num_atoms; i++)
815 : : {
816 : 0 : val = (int) at[i].valence;
817 [ # # ]: 0 : for (j = 0; j < val; j++)
818 : : {
819 [ # # # # ]: 0 : length += 1 + ( bDoubleBondSquare && BOND_DOUBLE == at[i].bond_type[j] );
820 : : }
821 [ # # # # ]: 0 : length += ( num_t_groups && at[i].endpoint );
822 : : }
823 : 0 : length += num_atoms;
824 [ # # ]: 0 : for (i = 0; i < num_t_groups; i++)
825 : : {
826 : 0 : length += (int) t_group[i].nNumEndpoints;
827 : : }
828 : 0 : length += num_t_groups;
829 : : }
830 : 502 : length++; /* +1 to save number of neighbors */
831 : 502 : pAtList = (AT_NUMB*)inchi_malloc(length * sizeof(AT_NUMB));
832 [ + - ]: 502 : if (pAtList) /* djb-rwth: addressing LLVM warning */
833 : : {
834 [ + - ]: 502 : if (!bDoubleBondSquare)
835 : : {
836 [ + + ]: 4681 : for (i = 0, length = 0; i < num_atoms; i++)
837 : : {
838 : 4179 : val = at[i].valence;
839 : 4179 : start = length++;
840 [ + + ]: 11723 : for (j = 0; j < val; j++)
841 : : {
842 : 7544 : pAtList[length++] = at[i].neighbor[j];
843 : : }
844 : : /* add endpoint */
845 [ + + + + ]: 4179 : if (num_t_groups && at[i].endpoint)
846 : : {
847 : 28 : pAtList[length++] = num_atoms + (int) at[i].endpoint - 1;
848 : : }
849 : 4179 : pAtList[start] = length - start - 1; /* number of neighbors before the list of neighbors */
850 : 4179 : pp[i] = pAtList + start; /* pointer to the <num.neigh.><list of neigh> */
851 : : }
852 : : }
853 : : else
854 : : {
855 [ # # ]: 0 : for (i = 0, length = 0; i < num_atoms; i++)
856 : : {
857 : 0 : val = at[i].valence;
858 : 0 : start = length++;
859 [ # # ]: 0 : for (j = 0; j < val; j++)
860 : : {
861 : 0 : pAtList[length++] = at[i].neighbor[j]; /* djb-rwth: buffer overrun avoided implicitly */
862 [ # # # # ]: 0 : if (bDoubleBondSquare && BOND_DOUBLE == at[i].bond_type[j])
863 : : {
864 : 0 : pAtList[length++] = at[i].neighbor[j]; /* a list of neighbor orig. numbers */
865 : : }
866 : : }
867 : : /* Add endpoint */
868 [ # # # # ]: 0 : if (num_t_groups && at[i].endpoint)
869 : : {
870 : 0 : pAtList[length++] = num_atoms + (int) at[i].endpoint - 1;
871 : : }
872 : 0 : pAtList[start] = length - start - 1; /* number of neighbors before the list of neighbors */
873 : 0 : pp[i] = pAtList + start; /* pointer to the <num.neigh.><list of neigh> */
874 : : }
875 : : }
876 : :
877 : : /* Add t-groups */
878 [ + + ]: 514 : for (i = 0; i < num_t_groups; i++)
879 : : {
880 : 12 : val = (int) t_group[i].nNumEndpoints;
881 : 12 : start = length++;
882 : 12 : nFirstEndpointAtNoPos = (int) t_group[i].nFirstEndpointAtNoPos;
883 [ + + ]: 40 : for (j = 0; j < val; j++)
884 : : {
885 : 28 : pAtList[length++] = nEndpointAtomNumber[nFirstEndpointAtNoPos + j];
886 : : }
887 : 12 : pAtList[start] = length - start - 1; /* number of neighbors before the list of neighbors */
888 : 12 : pp[num_atoms + i] = pAtList + start; /* pointer to the <num.neigh.><list of neigh> */
889 : : }
890 : : }
891 : : else
892 : : {
893 [ # # ]: 0 : inchi_free(pAtList); /* djb-rwth: fixing coverity ID #499598 */
894 [ # # ]: 0 : inchi_free( pp );
895 : 0 : return NULL;
896 : : }
897 : : } /* djb-rwth: ignoring LLVM warning */
898 : :
899 : : /* djb-rwth: fixing coverity ID #499598 -- pp uses pAtList values */
900 : :
901 : 502 : return pp; /* djb-rwth: ignoring LLVM warning: since a pointer is returned, memory should be freed in a function which calls *CreateNeighList */
902 : : }
903 : :
904 : :
905 : : /****************************************************************************/
906 : 1363 : void FreeNeighList( NEIGH_LIST *pp )
907 : : {
908 [ + + ]: 1363 : if (pp)
909 : : {
910 [ + - ]: 572 : if (pp[0])
911 : : {
912 [ + - ]: 572 : inchi_free( pp[0] );
913 : : }
914 [ + - ]: 572 : inchi_free( pp );
915 : : }
916 : :
917 : 1363 : }
918 : :
919 : :
920 : : /****************************************************************************/
921 : 218 : int BreakAllTies( CANON_GLOBALS *pCG,
922 : : int num_atoms,
923 : : int num_max,
924 : : AT_RANK **pRankStack,
925 : : NEIGH_LIST *NeighList,
926 : : AT_RANK *nTempRank,
927 : : CANON_STAT *pCS )
928 : : {
929 : 218 : int i, nRet = -1, nNumRanks = 1 /* value does not matter*/;
930 : :
931 : 218 : AT_RANK *nPrevRank = *pRankStack++;
932 : 218 : AT_RANK *nPrevAtomNumber = *pRankStack++;
933 : :
934 : 218 : AT_RANK *nNewRank = NULL;
935 : 218 : AT_RANK *nNewAtomNumber = NULL;
936 : :
937 [ + + ]: 218 : if (!pRankStack[0])
938 : : {
939 : 96 : pRankStack[0] = (AT_RANK *) inchi_malloc( num_max * sizeof( *nNewRank ) );
940 : : }
941 [ + + ]: 218 : if (!pRankStack[1])
942 : : {
943 : 96 : pRankStack[1] = (AT_RANK *) inchi_malloc( num_max * sizeof( *nNewAtomNumber ) );
944 : : }
945 [ + - - + ]: 218 : if (!pRankStack[0] || !pRankStack[1])
946 : : {
947 : 0 : return CT_OUT_OF_RAM; /* <BRKPT> */
948 : : }
949 : 218 : nNewRank = pRankStack[0];
950 : 218 : nNewAtomNumber = pRankStack[1];
951 : :
952 [ + - + - ]: 218 : if (nNewRank && nNewAtomNumber)
953 : : {
954 : 218 : memcpy(nNewAtomNumber, nPrevAtomNumber, num_atoms * sizeof(nNewAtomNumber[0]));
955 : 218 : memcpy(nNewRank, nPrevRank, num_atoms * sizeof(nNewRank[0]));
956 [ + + ]: 2224 : for (i = 1, nRet = 0; i < num_atoms; i++)
957 : : {
958 : : /* 12-12-2001: replaced Prev... with New... */
959 [ + + ]: 2006 : if (nNewRank[(int) nNewAtomNumber[i - 1]] == nNewRank[(int) nNewAtomNumber[i]])
960 : : {
961 : 58 : nNewRank[nNewAtomNumber[i - 1]] = (AT_RANK) i;
962 : 58 : nNumRanks = DifferentiateRanks2( pCG, num_atoms, NeighList,
963 : : nNumRanks, nNewRank, nTempRank,
964 : : nNewAtomNumber,
965 : : &pCS->lNumNeighListIter, 1 );
966 : 58 : pCS->lNumBreakTies++;
967 : 58 : nRet++;
968 : : }
969 : : }
970 : : }
971 : :
972 : 218 : return nRet;
973 : : }
974 : :
975 : :
976 : : /****************************************************************************
977 : : Int insertions sort
978 : : ****************************************************************************/
979 : 2 : int * iisort( int *list, int num )
980 : : {
981 : : int i;
982 [ + + ]: 6 : for (i = 1; i < num; i++)
983 : : {
984 : 4 : int tmp = list[i];
985 : 4 : int j = i - 1;
986 [ + + + + ]: 6 : while (j >= 0 && list[j] > tmp)
987 : : {
988 : 2 : list[j + 1] = list[j];
989 : 2 : j--;
990 : : }
991 : 4 : list[j + 1] = tmp;
992 : : }
993 : :
994 : 2 : return list;
995 : : }
|
