6
\$\begingroup\$

THE TASK

We are dealing with a string of symbols and need quick responses to queries of the following types:

  1. What is the position of the k-th occurrence of symbol X in the string?
  2. Reading from position k to the right, which symbol do we encounter first, X or Y?
  3. How many X symbols are there between positions i and j?

Input format:

In the first line: space-separated positive integers N and Q, corresponding to the length of the symbol string and the number of queries, N ≤ 1000000, Q ≤ 100000.

In the second line: a string of N printable ASCII characters (with codes from 33 to 126), without white spaces (except for the last newline, which is not part of the string and will not be queried). The string positions are indexed from 0.

In the following Q lines: first, a number r, equal to 1, 2, or 3, indicating the type of query.

When r=1, it is followed by the number k, a space, and a symbol s.

When r=2, it is followed by the number k and two symbols s1 and s2, separated by spaces.

When r=3, it is followed by symbol s and numbers k1 and k2k1, separated by spaces.

Numbers k, k1, k2 will always be less than N and not less than 0, and the symbols will be non-white printable ASCII characters.

Output:

You should print Q lines, each of which should contain the answer to query number i.

For type 1 queries, print a single number which is the position of the k-th occurrence of symbol s, or BRAK if such a position does not exist.

For type 2 queries, print the symbol s1 or s2, depending on which occurs first to the right of position k (inclusive), or BRAK if neither occurs. For clarity: if exactly one of them occurs, that is the one to print; if s1=s2, print this symbol if it occurs to the right of position k, and BRAK otherwise.

For type 3 queries, print the number of occurrences of symbol s between positions k1 and k2 (inclusive).

Examples

Input:

6 4
abcacb
1 1 c
1 2 c
1 3 c
2 2 a b

Output:

2
4
BRAK
a

The first 'c' symbol is at position 2, the second at position 4, and the third does not occur. Counting from position 2 (the first 'c' symbol), the 'a' symbol occurs before 'b'.


Input:

8 4
abccab1c
1 1 1
2 3 c b
2 3 4 d
3 c 0 7

Output:

6
c
BRAK
3

The first '1' symbol is at position 6. Counting from position 3, 'c' occurs before 'b' (because it's at position 3 itself). The string does not have symbols '4' or 'd'. There are three 'c' symbols in the string.


Input:

20 6
K1fpH;hjPeGg;ObjFvQh
2 6 f b
3 a 8 18
3 ; 6 16
1 1 b
3 P 4 10
1 1 O

Output:

16
b
0
1
14
1
13

Note

In this task, solutions are compiled with an additional -O2 flag, which results in better performance of the generated machine code but may also lead to additional warnings. For example, even though the input data in the tests have the correct format, it may be necessary to check the return value from scanf() in this (or a similar) way:

if (scanf("%f %d %c", &x, &k, &b) != 3) return 1;

MY CODE:

#include <stdio.h>
#include <string.h>

int findKthOccurrence(char *str, char s, int k) {
    int count = 0;
    for (long int i = 0; str[i] != '\0'; ++i) {
        if (str[i] == s) {
            ++count;
            if (count == k) {
                return i;
            }
        }
    }
    return -1;
}

char findFirstOccurrence(char *str, int k, char s1, char s2) {
    for (int i = k; str[i] != '\0'; ++i) {
        if (str[i] == s1 || str[i] == s2) {
            return str[i];
        }
    }
    return '0';
}

int countOccurrencesInRange(char *str, char s, int k1, int k2) {
    int count = 0;
    for (int i = k1; i <= k2 && str[i] != '\0'; ++i) {
        if (str[i] == s) {
            ++count;
        }
    }
    return count;
}

int main() {
    long int N, Q;
    int r, k, k1, k2;

    if (scanf("%ld %ld", &N, &Q) != 2)
        {
            printf("Invalid Input");
            return 0;
        }
    char s, s1, s2, str[N];
    if (scanf(" %s", str) != 1)
    {
        printf("Invalid Input");
        return 0;
    }
    for (long int i = 0; i < Q; ++i) {
        if (scanf(" %d", &r) != 1)
            {
                printf("Invalid Input");
                return 0;
            }
        switch (r) {
            case 1:
                if (scanf(" %d %c", &k, &s) != 2)
                {
                    printf("Invalid Input");
                    return 0;
                }
                int pos = findKthOccurrence(str, s, k);
                if (pos == -1) {
                    printf("BRAK\n");
                } else {
                    printf("%d\n", pos);
                }
                break;

            case 2:
                if (scanf(" %d %c %c", &k, &s1, &s2) != 3)
                {
                    printf("Invalid Input");
                    return 0;
                }
                char result = findFirstOccurrence(str, k, s1, s2);
                if (result == '0') {
                    printf("BRAK\n");
                } else {
                    printf("%c\n", result);
                }
                break;

            case 3:
                if (scanf(" %c %d %d", &s, &k1, &k2) != 3)
                {
                    printf("Invalid Input");
                    return 0;
                }
                printf("%d\n", countOccurrencesInRange(str, s, k1, k2));
                break;
        }
    }

    return 0;
}

I'm supposed to write this as an assignment for some C programming classes, and my code is checked by an automatic tool that I have no access to. It seems that it scored rather poorly, but I have no idea what might be wrong with it. Maybe someone here can help?

\$\endgroup\$
3
  • \$\begingroup\$ Is this assignment about finding an algorithm for string queries? It is written like a competition programming question. \$\endgroup\$
    – qwr
    Commented Jan 22 at 4:20
  • 1
    \$\begingroup\$ Are you allowed to sort the queries, as in Mo's algorithm? \$\endgroup\$
    – qwr
    Commented Jan 22 at 4:25
  • \$\begingroup\$ Is the problem-statement part of the question copy/pasted from a public web site? If so, please link it. Especially if there's an input-generator for testing, or some big sample test-cases for performance. (For uniformly-distributed random strings, How to generate a random string?, but interesting worst cases for brute force involve non-uniform inputs with some very rare characters.) But even if not, you're copying someone else's writing so that should be attributed for at least copyright purposes. \$\endgroup\$ Commented Jan 22 at 21:37

4 Answers 4

9
\$\begingroup\$

index type

You seem to choose (signed) long or int apparently at random when writing a loop. Prefer (unsigned) size_t i for array indexes.

names

I'm not super fond of uppercase identifiers like Q, but I suppose you were trying to match the problem statement, so OK, fine.

In code like this it is conventional for s to point at a string, and for c or ch to be a character within a string. So char s was slightly jarring, but again I suppose you were going for consistency with the problem statement.

dynamic allocation

    char ..., str[N];

The problem gave us a bound of a million on N. Prefer a static allocation here, to give the optimizer a little more to work with.

use libc

You wrote some simple loops, and they're very nice. An -O3 optimizing compiler likely manipulates them, perhaps by unrolling.

But memchr() has definitely been benchmarked to within an inch of its life. It certainly won't run slower than your loop. It might exploit SIMD vector processing if that's available on the target system which runs the queries. So use such standard building blocks where feasible.

BTW, I do hope you specified -march, or otherwise told the compiler it is allowed to use non-portable features of a particular chipset.


it scored rather poorly

I assume the score mostly depends on elapsed running time.

You mentioned there are some compiler warnings to be cleaned up, which might possibly cost demerits.

algorithm

set membership

We repeatedly test whether str[i] is within a set of length two:

        if (str[i] == s1 || str[i] == s2) {

IDK, maybe allocate a vector of length 256, clear it to zeros, and then set the s1 and s2 elements to 1. So we do a single dereference rather than a pair of comparisons. Bench it. Maybe it's faster?

pre-processing

We issue Q queries (many queries) against a megabyte string. So taking a moment to pre-process the string could be helpful, as Boyer & Moore teach us, when they tackle a different kind of problem. The idea here is to anticipate possible queries and memoize the results, just in case the result becomes relevant when we see the queries.

Also, we assume chief cost of these extremely simple queries will be dragging that giant megabyte string through the memory hierarchy, so the more we can do with each retrieved character, the better.

type 1

print a single number which is the position of the k-th occurrence of symbol s

(Note that 20 bits is enough to represent an index of a million.)

  • Define a tuple struct, using 8 bits for s and 24 bits for position.
  • Scan the N characters of the input string, outputting N tuples.
  • qsort() the tuples, costing O(N log N) time.

At query time create a (s, 0) tuple, and use binary search to identify where the s entries start. Advance the index by k to immediately read off the answer. Total query cost is O(log N) time, a significant savings over the O(N) time of your proposed solution.

type 2

print the symbol s1 or s2, depending on which occurs first to the right of position k

Exploit that same sorted datastructure. Do a pair of binary searches to find the beginning of the s1 entries and of the s2 entries.

Simple approach: linearly scan through the entries, looking for an index that exceeds k. Do this again. Compare the indexes you found, and return the appropriate character.

Fancy approach: Let succ(s1) be the successor of that character, so e.g. succ('X') == 'Y'. Do a pair of binary searches for the successors of s1 and s2. Now you know their exact ranges, and hence are able to do binary searches for indexes exceeding k.

Both approaches have time complexity that improves on the O(N) time of your proposed solution.

type 3

print the number of occurrences of symbol s between positions k1 and k2

Again we can adopt a simple or fancy way of exploiting our sorted datastructure. Assume we choose fancy.

Do a pair of binary searches to identify where s entries and succ(s) entries begin. Do a pair of binary searches to identify where the k1 and k2 indexes appear within those entries. Use subtraction to report the number of occurrences.

Again the time complexity significantly improves on the O(N) time of your proposed solution.

\$\endgroup\$
9
  • 1
    \$\begingroup\$ Note that current compilers won't auto-vectorize loops whose iteration-count couldn't be calculated ahead of the first iteration. So that rules out search loops that have an early-out. Library implementations of memchr and strchr using hand-written asm can go 16 or 32x faster than current compilers with byte-at-a-time C, on typical x86 with AVX2. (glibc detects CPU features at run time, so can also take advantage of AVX2 in a program compiled without -march=native). \$\endgroup\$ Commented Jan 21 at 23:28
  • \$\begingroup\$ I ended up posting a whole answer about SIMD instead of cluttering up comments. memchr as a building block for "magical" fast checking of many bytes could maybe give a 10x speedup for query type 2 vs. the OP's byte-at-a-time loop if we're lucky. Or maybe 2x to 4x with a portable scalar bithack to check 8 bytes at a time for v ^ x or v ^ y having any zero bytes (i.e. a match) for ISAs with SIMD, or without efficient branching on SIMD. \$\endgroup\$ Commented Jan 22 at 1:42
  • 1
    \$\begingroup\$ BTW, Boyer Moore and similar strstr algorithms preprocess the "needle", not the "haystack". A better analogy might be database indexes, but it gets the point across. \$\endgroup\$ Commented Jan 22 at 1:45
  • 2
    \$\begingroup\$ @Peter: given the input constraints, we don't even need to care about CHAR_BIT - we can legitimately abort if we get anything outside the range 33-126. \$\endgroup\$ Commented Jan 22 at 7:52
  • 2
    \$\begingroup\$ For fun, I updated my answer with a Godbolt link with a simple branchy binary search. For random k type-2 queries on a uniform random string (so hits are typically within 100 bytes or so), it's 3x slower for 960kB strings than the OP's byte-at-a-time search, on Skylake, with the same 2 chars so L2 hits on their pos lists, but L2 miss for linear. And 10.5x slower than a manually-vectorized AVX2 linear search. (For throughput; I didn't time latency. The AVX2 version usually avoids mispredicts so can pipeline L2 miss / L3 hits across queries, with out of order exec doing its job.) @MatthieuM. \$\endgroup\$ Commented Jan 26 at 2:53
7
\$\begingroup\$

The first thing that strikes me is that you haven't shown your tests of the code. Although you don't have access to the scoring system, you can easily prepare some basic tests of functionality that can be selectively compiled:

#ifdef TESTING

#include <stdio.h>
#include <stdlib.h>

int main(void)
{
    size_t errors = 0;
    errors += findKthOccurrence("aaa", 'a', 4) != -1;
    /* etc */
    return errors ? EXIT_FAILURE : EXIT_SUCCESS;
}

#else

/* your main here */

#endif

Then tell your compiler to define TESTING (e.g. gcc -DTESTING). And change the functions to accept the string as a const char* (which we should be doing anyway, since we don't intend to modify the content).

We can make our tests more sophisticated to let us know exactly which test failed, but this is enough to get started.

Once we have some robust tests, we can start making improvements.


The next thing I notice is that we're not making use of standard library functions such as strchr()¹ and strcspn(), but instead writing our own loops. While a good optimising compiler might be able to recognise the patterns involved, it's likely that we'll get lower performance than the highly-tuned functions in the library.

¹ or memchr() - see Peter Cordes's answer

In any case, we're using poor (linear) algorithms here. Since we will be performing several operations on each string, it would be better to transform the string into a more suitable representation for the queries.

For these queries, the appropriate structure would be an array of vectors (one vector for each character in our input character set), with each list containing the indexes (in order) where that character appears. Then all three types of query can be more quickly serviced (O(1) for the first query, O(log V) for the binary searches in the other queries, where V is the length of the particular character's vector).

We'll need an expandable array type, perhaps something like

struct array
{
    size_t capacity;
    size_t size;
    size_t content[];
};

static struct array *array_append(struct array *arr, size_t element)
{
    if (arr->size == arr->capacity) {
        /* allocate another 50% */
        size_t new_capacity = arr->capacity * 3 / 2;
        struct array *new_arr = realloc(arr, sizeof *arr + sizeof *arr->content * new_capacity);
        if (!new_arr) { return new_arr; }
        arr = new_arr;
    }
    arr->content[arr->size++] = element;
    return arr;
}

For reasonable implementations of realloc() which round up to the next block size, we can probably omit the capacity member and always reallocate (in most cases, we'll get the same pointer back).

Note that we don't need to store the whole line of input when making this conversion - we can do it character by character as we read using getchar().


Some review of main() (which should be written int main(void) to be specific that it takes no arguments:

        printf("Invalid Input");
        return 0;

It's good that we are validating the input. But we shouldn't be returning a success status when reading fails - we should use EXIT_FAILURE macro from <stdlib.h>.

Error messages should go to the error stream, and comprise complete lines (ending in a newline).

I would replace this with:

        fprintf(stderr, "Invalid Input\n");
        return EXIT_FAILURE;

This line is very dangerous:

if (scanf(" %s", str) != 1)

%s conversion with no field width is unbounded and can do great harm. However, we can simply consume newline and then use fgets(str, N, stdin) which does have a length limit.

main() is inconsistent with where it declares its local variables. Prefer to declare as close as possible to where they are initialised, and ideally in the smallest scope:

int main(void)
{
    int N, Q;
    if (scanf("%d %d", &N, &Q) != 2) {
        ⋮
    }
    char str[N];
    if (!fgets(str, N, stdin)) {
        ⋮
    }
    for (long int i = 0; i < Q; ++i) {
        int r;
        if (scanf(" %d", &r) != 1) {
            ⋮
        }
        switch (r) {
        case 1:
            {
                int k;
                char s;
                if (scanf(" %d %c", &k, &s) != 2) {
                    ⋮
                }
                ⋮
                break;
            }
        case 2:
            {
                int k;
                char s1, s2;
                if (scanf(" %d %c %c", &k, &s1, &s2) != 3) {
                    ⋮
                }
                ⋮
                break;
            }
        case 3:
            {
                int k1, k2;
                char s;
                if (scanf(" %c %d %d", &s, &k1, &k2) != 3) {
                    ⋮
                }
                ⋮
                break;
            }
        }
    }
}
\$\endgroup\$
8
  • \$\begingroup\$ I forgot to mention that on systems like Linux that normally allow unlimited overcommit of memory allocations, there's not much downside to just allocating enough space for the worst case with malloc(sizeof(uint_least32_t) * N) where N is the string length. Any virtual pages we don't touch won't ever get physical pages allocated to back them, and virtual address-space isn't a limiting factor even on a 32-bit system for a string with max length 1MiB so 93 * sizeof(T) = 372 MiB. (93 = 126-33, the ASCII range length). One downside is L1dTLB miss (L2 hit) from random access to different pages. \$\endgroup\$ Commented Jan 22 at 14:17
  • 1
    \$\begingroup\$ You can realloc each array to its size when we're done, but that won't compact them. So answering queries will also suffer TLB misses which we maybe could have avoided if we compacted the allocations, especially on x86-64 Linux where the hugepage size is 2MiB and Linux will opportunistically use transparent hugepages when most of a 2M region is populated. (Modern x86 L2 TLBs are big enough that the L1dTLB misses wouldn't cause page walks, but hit in the second-level TLB. Low-end ARM, IDK.) \$\endgroup\$ Commented Jan 22 at 14:24
  • 1
    \$\begingroup\$ Also, since glibc malloc uses mmap for whole pages for big allocs, all the arrays would have the same alignment relative to a 4K boundary. So their first bytes would all alias the same set in L1d cache. Keeping the size (and capacity) outside the dynamic allocation would avoid that and avoid a level of indirection to access them. With different sizes the binary search startpoints would be different. I like my answer's version with an array of structs of size/pointer (but a 3rd member makes the size not a power of 2. A fixed growth increment of 8192B can avoid separate capacity?) \$\endgroup\$ Commented Jan 22 at 14:30
  • \$\begingroup\$ I think the array of vectors is definitely on the right path, but honestly I'd recommend just a giant allocation: all max-sized vectors catenated one after the other, for a total of 372 MiB. It may not be faster but it would be simpler than implementing a resizable array. Once the sizes are known, it could be compacted... but I would not bother unless explicitly told to. Also, at 372 MiB, this massive array can be statically allocated, no need to even mess with malloc (it can't be stack allocated though). \$\endgroup\$ Commented Jan 22 at 16:12
  • 1
    \$\begingroup\$ @MatthieuM.: I updated my answer with some experimental results from playing around with a build_idx function. Indeed, static storage is much faster (like 3x) than 94 separate malloc/free, since for a repeat loop I'm freeing every time but the static storage version gets to reuse the space. That spends most of its time in system calls. But even worse is one big malloc/free, since that tempts Linux into zeroing some hugepages for the pos arrays for a uniformly distributed 240000 byte string (so each pos list is on average about 2.5 pages long, very sparse.) \$\endgroup\$ Commented Jan 23 at 19:56
6
\$\begingroup\$

General Observations

For the most part the code is portable which is good, however, variable length arrays that were introduced in C99 were made optional in C11 which means they may not compile in all versions of C and definitely won't compile in any version of C older than C99. In my Visual Studio 2022 version of C the following line does not compile:

    char s, s1, s2, str[N];

To make the code even more portable include stdlib.h and use the system defined constants EXIT_SUCCESS and EXIT_FAILURE to indicate the applications status on exit. While 0 for EXIT_SUCCESS and 1 for EXIT_FAILURE are fairly common, to ensure system compatability you should really use EXIT_SUCCESS and EXIT_FAILURE.

When a program fails it is better to return EXIT_FAILURE rather than EXIT_SUCCESS, such as the input code is currently doing.

Checking the return value of system functions such as scanf() is definitely a best practice.

Variable Names

Use meaningful variable names, single letter variable names such as N, Q, r, k, k1, k2, s, s1, and s2 are not meaningful variable names. The function names are descriptive, the same method used to generate the function names should be used to generaet the variable names.

DRY Code

There is a programming principle called the Don't Repeat Yourself Principle sometimes referred to as DRY code. If you find yourself repeating the same code mutiple times it is better to encapsulate it in a function. If it is possible to loop through the code that can reduce repetition as well.

The input code is repetitive.

Complexity

The function main() is too complex (does too much). As programs grow in size the use of main() should be limited to calling functions that parse the command line, calling functions that set up for processing, calling functions that execute the desired function of the program, and calling functions to clean up after the main portion of the program.

There is also a programming principle called the Single Responsibility Principle that applies here. The Single Responsibility Principle states:

that every module, class, or function should have responsibility over a single part of the functionality provided by the software, and that responsibility should be entirely encapsulated by that module, class or function.

The code is missing at least 1 function which is an input function that does error checking. Quite possibly each case in the switch statement should be a separate function call.

\$\endgroup\$
1
  • 1
    \$\begingroup\$ The return 0; from main() is in C99, and OP uses other C99 features as well. \$\endgroup\$
    – Davislor
    Commented Jan 21 at 3:08
4
\$\begingroup\$

(This isn't so much a code-review as what I was going to comment on other answers that mentioned optimized asm libc implementation of memchr and so on. And/or how I'd go about solving it, especially in non-portable C that could use <immintrin.h> or <arm_neon.h>. It's disappointing that compilers are almost completely unable to make truly efficient code from simple source like yours for simple problems like these.)

Performance of brute-force loops

There's no portable way with C to take full advantage of the SIMD brute force modern CPUs can apply to problems like this, comparing 16 or 32 bytes at once and doing one test/branch to check for any matches. For non-huge problems, carefully applied brute-force should be very good, like for strings that stay hot in L2 cache. (often 256K on older CPUs, 1M or larger on some recent x86). You have the option of building some index data structures, but depending on typical queries (how far we have to scan to get a match), it might not be worth it if we could write a program that was anywhere near as good as what we could do with x86 AVX2 or AArch64 ASIMD intrinsics (or hand-written asm).

(SIMD = Single Instruction Multiple Data, like SSE2 pcmpeqb to do 16 byte-compares with one instruction that modern CPUs can run at 3 per clock cycle, not much worse than integer add. Processing the results so we can branch on them is non-trivial and should be amortized over multiple SIMD vectors if we don't expect a match very soon. See Why does glibc's strlen need to be so complicated to run quickly? for discussion of both the the SIMD hand-written asm it uses on x86, and the "portable" bithack long-at-a-time C version (but not really; strict-aliasing UB) used as a fallback on a few targets like MIPS.)

Current compilers unfortunately won't auto-vectorize loops whose iteration-count couldn't be calculated ahead of the first iteration. So that rules out search loops that have an early-out. Library implementations of memchr and strchr using hand-written asm can go 16 or 32x faster than current compilers with byte-at-a-time C, on typical x86 with AVX2. (glibc detects CPU features at run time, so can also take advantage of AVX2 in a program compiled without -march=native).

Take advantage of hand-written asm library functions

C's library exposes a few of the things one can do fast in non-portable asm, like finding the position of a given byte (memchr). Unfortunately missing some very useful things, like position of the first mismatch between two strings / arrays, but fortunately we don't need that.

You know N, the length of the string, so there's no need to be searching for the terminating 0 all the time. Use memchr not strchr: it has lower startup overhead since it can just check the size and know it's allowed to read a whole vector, vs. strchr having to worry that a 32 byte load could go off the end of a page into an unmapped page and segfault. So it has to check pointer alignment and branch on a bunch of special cases.

To amortize call overhead for query type 2 (X or Y next), you might use memchr(pos, x, 128) and memchr(pos, y, 128) as building-blocks in a loop doing pos+=128 (and checking if we're getting near the end of the string). If there's actually a Y 20 bytes from pos, we'll check up to 128 bytes looking for an X before we call the second memchr to find the Y, but there is some startup overhead in each memchr.

The block size is a tunable parameter you can play with. 512 bytes is only 8x 64-byte cache lines, or 16x 32-byte vectors. 128 bytes (2 lines, 4 YMM vectors) or even 64 or 96 seems good with glibc on Skylake (AVX2), with uniformly distributed symbols in the string (so it's very rare for the first occurrence to be farther than 128.

For uniform symbols, this is about 1.4x faster than the naive pure C loop (which also includes checks for \0 instead of using the size, but if we're checking every character separately that's fine) on Skylake with glibc, compiled with clang 16 -fno-plt -O3 -march=skylake -mbranches-within-32B-boundaries. (Average about 136 cycles per random query vs. 100 +- 3 for calling memchr with a chunk size of 64, vs. about 103 cycles for a chunk size of 128. No I/O, just a fast random number to generate a k, timing 100M queries, with a 240000 byte string that mostly stays hot in L2 cache.)

If the next-patch position is a few thousand bytes away, I'd guess this might run 10x faster with a chunk size of 512 than the simple version on a typical modern x86 with AVX2 like Skylake or Zen 2, for data hot in L2 cache. (i.e. in cases designed to be slow with a naive algorithm, where you have to scan far.)

char findFirstOccurrence_memchr(const char *str, size_t len, size_t k, char x, char y)
{
    const int chunk_size = 128;

    size_t len_remain = len-k;  // bytes from pos to end of string
    for (const char *pos = str+k ; len_remain >= chunk_size ; pos += chunk_size, len_remain -= chunk_size) {
        const char *xpos = memchr(pos, x, chunk_size);
        const char *ypos = memchr(pos, y, chunk_size);
        if (xpos && ypos) {
            return xpos < ypos ? x : y;
        } else if (xpos) {
            return x;
        } else if (ypos) {
            return y;
        }
         // else neither found a match
    }

    if (len_remain){
        // same as loop body but with len_remain instead of a full chunk_size
        // we could have done min(len_remain, chunk_size) inside the loop instead of peeling this
        const char *xpos = memchr(pos, x, len_remain);
        const char *ypos = memchr(pos, y, len_remain);
      // For production use I'd probably factor this into a helper function or macro.
        if (xpos && ypos) {
            return xpos < ypos ? x : y;
        } else if (xpos) {
            return x;
        } else if (ypos) {
            return y;
        }
        // else neither found a match
    }
    // return '0';   // why ASCII decimal digit '0'?
    return 0;        // ASCII nul '\0' makes more sense to me.
}

Toby suggested a good way to avoid the if (len_remain){ cleanup part: round up the allocation for the string buffer to a multiple of the chunk size, padding the end with something like 0 (ASCII '\0') that can't be a false-positive match. This works because we're only handling ASCII string data, not binary where every char value is possible.

This is significantly worse than what you could do with SSE2 intrinsics, like this sketch that omits some of the pointer math. It's like memchr but with bitwise OR of two compare results. If desired, we can unroll the way glibc does. An AVX2 version would be identical but with __m256i and _mm256_... intrinsics. An AVX-512 / AVX10 version using compare-into-mask could use kortest to branch on either compare result having a set bit.

// See also the Godbolt link in the last section for a working AVX2 version
#include <immintrin.h>

int findFirstOccurrence_SSE2(const char *str, size_t len, size_t k, char x, char y)
{
  __m128i vx = _mm_set1_epi8(x);  // vector with all 16 elements = x
  __m128i vy = _mm_set1_epi8(y);

  const char *pos = str + k;
  const char *endpos = str + len;  // With a padded input buffer, ok to read up to chunk_size past string[N-1]

  for(; pos < endpos ; pos += 16){
      __m128i v = _mm_loadu_si256((__m256i*)pos);  // unaligned load
      __m128i cmpx = _mm_cmpeq_epi8(v, vx);  // vector of 0 or -1 bytes
      __m128i cmpy = _mm_cmpeq_epi8(v, vy);
      __m128i match = _mm_or_si128(cmpx, cmpy);  // true in elements where either had a match
      // and then the rest is identical to a memchr loop
      // optionally do a couple vectors unconditionally to amortize test+branch, and then sort out which vector had the match later.
      unsigned mask = _mm_movemask_epi8(match);   // extract the high bits of  each element to an integer bitmask
      if (mask) {
          // test eax,eax / jcc  then  tzcnt or BSF
          return pos[ __builtin_ctz(mask) ];  // find the match position and get the char there
      }
  }

  // if we didn't pad the input buffer, the full-vector loop might exit early
  if (the total search region was at least 16 bytes) {
     // check the last 16 bytes of the string, 
     // partially overlapping with the last vector from the loop
     __m128i v = _mm_loadu_si128((__m128i*)str+len-16);
     compare and check mask as in the loop.
  } else {
     scalar fallback, or _mm_loadu_si64 then scalar
  }
  return 0;
}

ARM / AArch64 don't have a pmovmskb instruction that pulls the high bit from each vector. For AArch64 memchr and similar, you normally reduce 128-bit to 64-bit with a shift right and insert instruction, then move that from SIMD to general-purpose integer and check it for non-zero. (If so, then bit-scan it for the match position with rbit (bit-reverse) / clz (count leading zero), or a bithack to isolate the lowest set bit + clz.)

AArch64 CPUs usually have efficient SIMD → integer transfers, but many 32-bit ARMv7 CPUs did not, so it wasn't as easy or profitable to use NEON for memchr type problems. (The core had to stall and wait for SIMD instructions to finish because they normally weren't tightly coupled enough to track dependencies.) Counting matches over a range is fine, but branching every 32 or 64 bytes on compare results was a problem. But not for AArch64 usually, and not for x86. On some x86 CPUs, SIMD vector to integer data transfers have highish latency, but it's always pipelined so not big throughput problem.


Counting matching chars in a range

Efficient manual vectorization looks like this with AVX2: https://stackoverflow.com/questions/54541129/how-to-count-character-occurrences-using-simd - pcmpeqb / psubb in the inner loop, using byte compare results as 0/-1 integers to increment a vector of uint8_t counts. And psadbw against zero in an outer loop to horizontally sum unsigned bytes to uint64_t in 8-byte chunks, before the uint8_t counters can overflow. This can run about 1.5x 32 bytes per cycle on modern Intel, or 2x 32 byte per cycle on modern AMD (4 SIMD ALU ports), if data is hot in L1d cache.

You can get compilers to vectorize things like counting matching bytes. (Between two arrays, or with a fixed character against arr[i] is basically the same problem, accumulating the compare results is the hard part. a[i] == b[i] vs. a[i] == c is just the difference between another vector load vs. broadcasting a scalar ahead of the loop). Compilers can sometimes do an ok job, like what @harold managed to get clang to do in Count the number of mismatches between two arrays with some careful programming to hand-hold the compiler toward the asm he knew he wanted it to make. Otherwise expect quite a bit worse, but still hopefully better than scalar.

But you said programs are only compiled with -O2. If that's GCC, you won't get auto-vectorization except of the easiest loops, like a constant iteration count that's a multiple of 32 or something.

It might help to write the source like count += (str[i] == c); to encourage the compiler to make branchless asm. This pattern of conditionally adding based on a compare result is basically the same as Stack Overlow's branch-prediction canonical Q&A, https://stackoverflow.com/questions/11227809/why-is-processing-a-sorted-array-faster-than-processing-an-unsorted-array . At least for scalar code.

(For vectorization, that code is adding a 32-bit array element to a 64-bit counter, but here we're just incrementing by 1 or not so you want to work with 16x 1-byte elements as long as possible, not unpacking to four vectors of 4x 32-bit elements right away.)

Query type 1: kth occurrence

You could just call memchr k times. But if the density of occurrences is high-ish, like more than one per 32 bytes, this is much worse for large k than we can do with SIMD brute force. You'd maybe want to build this out of the counting code, perhaps counting over a range like 512 or 1024 bytes and then continuing forward if we haven't reached k, or backtracking if we've gone past.

It's non-trivial to do it even with manual vectorization with SIMD intrinsics. To count efficiently, you normally want to avoid reducing to scalar very often. Probably work in large batches, and when getting close use a reduce-every-16-bytes loop.

For all of these, building arrays of symbol-occurrence positions can be a big win, but type 1 maybe moreso, especially if you're trying to hack around the lack of portable SIMD search in languages like C.


Tuning an array of symbol_occurrences positions

C23 unsigned _BitInt(24) might be useful to save some cache footprint on arrays of occurrence position, if compilers pack that into 3 bytes. Or not, since spatial locality is poor when binary searching and instruction latency from load to compare is critical (for branchy to catch mispredicts earlier, for branchless it's a loop-carried dependency.) But if it lets more of your "index" stay hot in cache more of the time across queries it could be worth making each access cost more CPU instructions.

Speaking of searchable sorted data structures, an array isn't necessarily optimal, especially if you can use SIMD efficiently for searching (like modern x86 can, and also AArch64). https://stackoverflow.com/questions/56521133/what-is-the-most-efficient-way-to-implement-a-bst-in-such-a-way-the-findvalue - an implicit tree has good locality without wasting any space on pointers. And it doesn't have to be a binary tree; each node can be a whole cache line of an n-ary tree, so you could have a 17-ary tree with 32-bit integers and 64-byte cache lines. (Use aligned_alloc to get alignas(64) for dynamic allocation.) It takes a bit of computation to answer type-1 queries (kth element) from a tree instead of array, though. Perhaps it helps to make a tree where every element is present as a leaf, so internal nodes aren't the only place you can find those elements. (For a balanced tree, does that make the final depth into a sorted array?)

Turning the arrays of symbol positions into implicit trees would add some preprocesing work, so might only be worth it if Q is large. But we know Q before we even read the string so we can decide ahead of time how much preprocessing it's worth doing on the string!

For low Q we might even just brute-force search the original string, especially if N is large. And especially if stdin is a regular file we can mmap, so we might end up never touching some of its pages if we don't pre-process and none of the queries make us read a byte in that range. (That means writing multiple versions of the rest of the program with if(strategy) or switch(strategy), or using function pointers to plug in different implementations to dispatch after parsing a query.)

Presumably for this assignment, a fixed strategy is totally fine, and struct { uint32_t *pos; size_t len; } symbol_occurrences[UCHAR_MAX+1]; would be plenty fast, sorted arrays of indices. (No actual sorting work necessary, one pass over the string appending positions to each symbol array will make sorted output.)

Actually we're guaranteed that the symbols are non-whitespace printable ASCII, so the range of symbols is 33..126 ('!' = ' '+1 to '~' but in this case that's not a more meaningful way to write the constants.) Instead of UCHAR_MAX+1ULL = 1ULL<<CHAR_BIT (typically 256), our array size can be 126-33 = 93, and we index it with arr[c - 33] (where we'd want to put the magic numbers into a named constant like #define SYMBOL_MIN 33 // one past ASCII space.)

#include <stdint.h>
#include <stdlib.h>

#define SYMBOL_MIN   33  // one past ASCII space
#define SYMBOL_MAX   126 // one before ASCII  DEL
#define SYMBOL_RANGE (SYMBOL_MAX + 1 - SYMBOL_MIN)  // array size

// TODO: use as narrow a type as possible for pos, perhaps unsigned on a machine where that's 20 or more bits (enough for max N) but less than 32, like some DSPs
// or with C23, perhaps unsigned _BitInt(20) or unsigned _BitInt(24)
//  .len might as well also be uint32_t, but size_t will usually be pointer width and structs will be padded for alignment anyway if we use a narrower type.
static struct { uint32_t *pos; size_t len; } 
     symbol_occurrences[SYMBOL_RANGE] = {};

// we use size_t for local vars, only the narrow position type in the array.
// Using uint32_t everywhere for the string size would also be fine, but not unsigned; it could be too short.
void build_idx(const char *str, size_t str_N)
{
    // allocate max size, resize later.
    //  Lazy allocation by the kernel avoids wasting actual RAM for pages we never touch
    for (int i = 0 ; i < SYMBOL_RANGE ; i++){
        symbol_occurrences[i].pos = malloc(str_N * sizeof(symbol_occurrences[0].pos[0]));
        // .len already zeroed
        if (! symbol_occurrences[i].pos) {
            abort();
        }
    }
    // Alternative: could start with one big allocation and carve it up with pointers into it.

    for (size_t pos = 0 ; pos < str_N ; pos++){
        unsigned c = (unsigned char)str[pos];  // array indices in 0..255 not -128..127 range!
        if (c < SYMBOL_MIN || c > SYMBOL_MAX){
            return;
            // error, abort() or make this non-void and return 0
        }
        size_t list_idx = symbol_occurrences[c - SYMBOL_MIN].len++;
        symbol_occurrences[c - SYMBOL_MIN].pos[list_idx] = pos;
    }

    // Optionally free up space in the tails of our arrays.
    // This is unlikely to actually compact them closer together for dTLB locality or transparent hugepages
    // If we wanted that, we'd have to manually memcpy, or memmove if we started with one huge allocation.
    for (int i = 0 ; i < SYMBOL_RANGE ; i++){
        void *new_alloc = realloc(symbol_occurrences[i].pos, 
                    sizeof(symbol_occurrences[0].pos[0]) * symbol_occurrences[i].len);
        if (new_alloc) {
            symbol_occurrences[i].pos = new_alloc;
        } // else realloc failed, leaving the old allocation still there.
          // For most implementations, impossible when shrinking
    }
}

Notice the use of unsigned char for indexing arrays. In cases where you want a lookup table of all UCHAR_MAX+1 possible byte values instead of a range-check, make sure you use unsigned char so your indexes are non-negative even with unexpected input.

OSes like Linux normally allow unlimited overcommit of memory allocations, so there's not much downside to just allocating enough space for the worst case with malloc(sizeof(uint_least32_t) * N) where N is the string length. Any virtual pages we don't touch won't ever get physical pages allocated to back them, and virtual address-space isn't a limiting factor even on a 32-bit system for a string with max length 1MiB so 94 * sizeof(T) = 376 MiB of virtual address-space. (But even rare symbols have their arrays in a page by itself.)

Instead of overcommit, another way to avoid separate size and capacity struct members (like std::vector) would be a fixed growth increment of 8192B. Check .len for crossing a multiple of 8192 and realloc. i.e. capacity is inferred from .len as (.len + 8191) & -8192, i.e. .len rounded up to the next multiple of 8192.

One downside to big initial allocations is L1dTLB misses (L2 TLB hits) from random access to different pages. Modern x86 L2 TLBs are big enough that the L1dTLB misses would hit in the second-level TLB instead of causing a page-walk so it's not a total disaster. On low-end ARM, IDK, could be worse. (Skylake stats - L1dTLB size of 64 entries for 4K or 2M pages, L2 TLB 1536 entries.) And if some symbols are rare, those arrays can just stay cold.

You can realloc each array to its size when we're done, but that won't compact them. So answering queries will also suffer TLB misses which we maybe could have avoided if we compacted the allocations, especially on x86-64 Linux where the hugepage size is 2MiB and Linux will opportunistically use transparent hugepages when most of a 2M region is populated. Compaction into a single allocation with memmove is possible. (Or instead of just copying, compact into an implicit tree for cache locality of searching. As discussed above in the SIMD section, an implicit 17-ary tree is good if you're SIMD searching it, otherwise maybe a binary or 3-ary or 5-ary tree for scalar code that doesn't have SIMD's fixed overhead for finding which element of a group matched.)

Glibc malloc uses mmap for whole pages for big allocs, returning a pointer that's 16 bytes past the start of a page. So all the arrays would have the same alignment relative to a 4K boundary. So their first bytes would all alias the same set in L1d cache.

But since (unlike Toby's answer), I'm keeping the size outside the dynamic allocation, we don't need to access that position for most queries. Presumably the array of pointer/len structs would stay hot in cache, and with that the first access to the dynamic allocations would be at a position that depends on the query, or in the middle (for a binary search) which depends on the length of that symbol-position entry. With different counts for different symbols, we hopefully don't get too many conflict-misses in cache.


TODO: make a single large allocation (or even a static uint32_t pos_storage[sizeof(uint32_t) * MAX_N];). As discussed in comments with @MatthieuM., you don't even need an array of pointers, you can just calculate the start of each symbol's portion with idx = (c - RANGE_MIN) * N. (Using current N, not MAX_N, to get good locality for small N, and to avoid having all arrays start a multiple of 4K from each other for cache aliasing and other effects.)

This is also very good for very small N like the test-cases. Zero or one calls to a library allocation function and minimal malloc bookkeeping overhead, and not too sparse use of the space.

Playing around with this on my Arch GNU/Linux desktop (gcc 13.2, i7-6700k), one big allocations can be slower (for with str_N = 240000) because it tempts Linux into zeroing 2M hugepages. But we don't have this problem using BSS space (static or global variables), unlike with dynamic allocation where my test code was was freeing and allocating between calls of a repeat loop. Static storage runs about 250 ms for 1024 repeats, separate allocations for each symbol-position array about 750 ms, one big allocation that we free between each build_idx runs about 6800 ms.

Just for fun, I was playing around with optimizing that part. Godbolt link (messy, lots of commented-out and ifdefed versions).

I was assuming compilers would be smart enough to do symbol_counts[c - SYMBOL_MIN] with a subtraction on the pointer outside the loop. But no, even if you write (symbol_counts-SYMBOL_MIN)[c], or even uint32_t *offset_counts = symbol_counts - SYMBOL_MIN; outside the loop, GCC really wants to put the subtract inside the loop. (It's UB in ISO C to form a pointer outside the bounds of the array, even if you don't deref. But it's totally fine in x86 asm, so the compiler could do it. Trying to hand-hold the compiler into doing it for us treads on thin ice. I ended up having to use asm("" : "+r"(offset_counts)); (like Benchmark::DoNotOptimize, a black box for the compiler) to get GCC to make something like the asm I wanted to benchmark.

That did save an instruction in the loop, but barely any speedup for N=240000, the front-end seems not to be the bottleneck on my i7-6700k Skylake. Rather L1dTLB misses (which indeed hit in L2TLB, significant counts for dtlb_store_misses.stlb_hit) and memory access. Using uint16_t instead of uint32_t for SymPos_type i.e. for the position arrays, does speed things up from about 245 ms to about 180 ms, for 1024 repeats of building the position data at 4.2GHz.

(Fun fact: POSIX reserves all symbols of the form foo_t, so that's not a safe portable choice for your own custom type names, unfortunately. I ended up making a lot of local vars size_t while trying to get the compiler to optimize the subtract into the original array address, so they didn't have to potentially do wrapping to some narrower unsigned type, but I still wanted a typedef for the position size to be able to see how uint16_t performs for cases that don't overflow, or unsigned _BitInt(24) which is slower for building the indices on my system (290 ms vs. 240 ms for unsigned _BitInt(32) or uint32_t, vs. 181 ms for uint16_t, for 1024 iters of build_idx with a 240000 byte string, with arrays in static storage, with clang 16 -O3 -march=skylake -fPIE. With the character range check commented out so there's no early-out in the loop, letting clang unroll). But I haven't tried answering queries from it which would be random access that benefits more from staying hot in cache.)


Benchmark results from optimized versions

Playing around with this some more (Godbolt), optimizing the type-2 query function using memchr or AVX2 _mm256_cmpeq_epi8, it seems that unconditionally doing 2 vectors (64 bytes) is enough to have a high chance of a match in the first iteration for at least one of the characters. (2x 64 = 128 chances to match, from a range of 94 different symbols). With only one vector (32 bytes) we get a 20% overall branch mispredict rate for a test loop that generates random queries, vs. 10% for this version. (1350 ms overall for 100 M queries, vs. 860 to 900 ms overall. Or about 3.56 Gclocks in user-space, so about 36 clocks per query vs. 56, using core clock cycle performance counters from perf stat --all-user -etask-clock,context-switches,cpu-migrations,page-faults,cycles,instructions,uops_issued.any,uops_retired.retire_slots,idq.mite_uops,branches,branch-misses when compiled with clang.)

Also, 64 seems to be a good chunk size for memchr; apparently its startup overhead is not bad. Slightly better than 96 or 128 byte sizes, getting about 100 user-space cycles per query in the same test as the AVX2 version. vs. about 137 to to 143 cycles for the naive loop on a uniformly distributed random string with two different symbols, for uniform random k that hits in L2 cache.

Hard to imagine that binary search over a sorted array would be much faster, either branchy or branchless, especially if that array has 4x the footprint so we get L2 misses. Searching an index has a worst case that's much less bad, though, and there's no guarantee the input is uniformly distributed like I was testing with. An index scales to strings that don't fit in L2 or even L3 cache, and handles queries about rare symbols efficiently, without having to scan many KiB of data to find them.

Update: findFirstOccurrence_index (Godbolt) with two branchy binary searches (with no fallback to linear SIMD scan for a small range) is about 11x slower (45.4 Gcycles for 100M iters) than findFirstOccurrence_AVX2, and 3x slower than naive (OP's code). On Skylake, for 100M type-2 queries on a string of 960000 bytes, with random k but fixed c1='A' and c2='b' so the two index buckets (of about 48KiB each of uint32_t) can stay hot in L2 cache (256K) but the raw string is almost 4x too big. (So it's intentionally biased in favour of the binary search, which has asymptotic O(log(N)) runtime, but it still loses badly with all those branch mispredicts. It might perhaps do better with cmov data dependencies, but then it wouldn't have speculative execution to prefetch and you might want SW prefetch. And then L2 load-use latency would be on the critical path. Or L3 for a more realistic case with random choices of characters so it's not the same 2 position-lists every time.)
I didn't implement any clever initial-guess strategy yet, based on the assumption of uniform random data (like size_t mid = (unsigned)(bucketlen * k) / (unsigned)str_len to use 32-bit division so it's not disastrously slow on Intel before Ice Lake, or better with libdivide for repeated division by the same integer.)
Interestingly, even with this near-max string size that misses to L3 for each query with linear searching, the AVX2 linear search version doesn't slow down much like 1.05 seconds vs. 0.86 sec for 240kB. (Skylake-client has pretty good L3 latency and bandwidth.) I am testing query throughput not latency so there's very good overlap of L3 misses between queries, out-of-order exec doing exactly what it was designed for. Uncommenting k1 |= loop_carried; doesn't change binary search speed (because it has mispredicts within each iteration), but does slow down the AVX2 linear search by about a factor of 2 (less than 2x for L2 hit sizes (62c / query), a bit more than 2x for L3 hit sizes (91c / query).)

An index is good for count-in-range over large ranges, though. It still needs two (binary) searches of the index, but looking at the string directly doesn't allow any early-out for large ranges. (Still, AVX2 can count matching bytes fast enough to keep up with L2 or L3 bandwidth). And obviously perfect for type-1 queries, just one load from a trivially-calculated index.


Benchmarking loop: scanf text queries was way slower than the queries.

Reading queries from a file with fscanf, rewinding the same 10k-query file with ftell / fseek, led to only about 16% to 10% of the program's total time spent in main, the rest in libc I/O and kernel code. And this was after replacing printf output with assignment to volatile int sink, which makes the compiler do the work but doesn't waste huge amounts of time formatting numbers into decimal strings and calling C stdio library code, which is costly even with output redirected to /dev/null (so it's full-buffered).

So I replaced that with a low-overhead very simplified 32-bit xorshift PRNG that's sufficient to distribute k randomly and get L1d cache misses (but mostly L2 hits because I was testing with a 240000 byte array on Skylake.) Read the file once, then randomize k for 100M queries, and store a result to a volatile int sink. Then we can slap a repeat loop around just the work we want to measure, so running the whole thing under perf stat or perf record gives us performance counters for just what we want to measure. Over 99% of the time spent in main (after callees inline into it, for the version without calls to memchr.)

\$\endgroup\$
3
  • \$\begingroup\$ Wow! So many things... Thank you. \$\endgroup\$
    – J_H
    Commented Jan 22 at 2:21
  • 2
    \$\begingroup\$ For this particular case, we could avoid "checking if we're getting near the end of the string" by over-allocating its storage (rounding up to the next whole block). We could avoid even having to check the return value is in range if we NUL-fill the excess. \$\endgroup\$ Commented Jan 22 at 8:05
  • \$\begingroup\$ @TobySpeight: Thanks, great idea. Updated my answer to mention that, and a section about allocating pos arrays like I was commenting about under your answer. \$\endgroup\$ Commented Jan 22 at 16:16

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.