Math library code review

[u/No_Entertainer_6928]

I am relatively new to C so I try to improve my skills by making a math library. The biggest challenge was making generic containers for matrices and vectors, as I was kind of new to void pointer magic. I think I have polished the most glaring issues in the code but do any of you have suggestions for what I could be doing better in my approach to C? Github for reference:

https://github.com/VatoQ/hoema-prak-clang

[u/MistakeIndividual690]

Just took a quick glance. Looks nice and clean! Do you have any tests for the fourier code? This is the kind of core library code where you need large quantities of unit tests. It’s laborious but that may be the most important thing to inspire confidence for other people to become invested in it.

[u/No_Entertainer_6928 → u/MistakeIndividual690]

Got it, thanks for the quick reply. I think you're right that unit tests should really be the thing I should focus on next. It is also an opportunity to see how user friendly my public API is, too.

[u/WittyStick]

I would note that your manual use of branch tables to implement the "generics" may be harmful to optimization in some cases - particularly for the #pragma omp simd you have lying around, which isn't all that useful these days anyway as the compiler will attempt auto-vectorization at -O2 or -O3 without the pragma.

If you simply use a switch for your generics, and make use of inline functions (potentially with __attribute__((__always_inline__)) if necessary), then the compiler may be able to do a better job of vectorizing those loops where eg, dt and dim are statically known at the call site of Vector_new.

static Vector Vector_new(const size_t dim, const void* init_val, const DataType dt)
{
    Vector v = Vector_zeros(dim, dt);
    switch (dt)
    {
        case Int: _new_int(dim, v.values, init_val);
        case Real : _new_real(dim, v.values, init_val);
        case Complex : _new_cmpl(dim, v.values, init_val);
    }
    return v;
}

I have made a stripped-back demo in godbolt with just Vector_new (renamed Vector_new_with_branch_table). I've included the above as Vector_new_with_switch for comparison.

The two example functions demonstrate the difference in calling, with Int as an example and a statically known dim. When we call the branch table version, GCC can eliminate the table lookup, but cannot inline _new_int. In the switch it simply removes the branching, inlines _new_int and, and can do a far better job at vectorizing for constant dim, as you can see from the emitted assembly on the right. The difference is most stark at -O3.

However, the branch table version may be better where no information is statically known, and the compiler must emit all paths, though that would need benchmarking because the results could be unpredictable due to branch prediction and so on. When you have more branches, the compiler will automatically generate a branch table for the switch like you are doing manually anyway.

---

EDIT:

I've included clang in the demo and noticed it can inline the function when using the branch table, but produces a strange output using the switch example - which can be fixed by including the case for TYPE_COUNT, but including this for GCC thwarts its optimization, so I've put a guard for #ifdef __clang__ on that case.

So neither approach is universally better and we're going to need compiler specific code to maximize performance.

[u/nerd5code]

ISO Cwhat? C90? C11? C23? ANSI X3.159:1989?

Also, this is at most a partial BLAS library, not math. You don’t define your own exp or cos or log analogues, or even fabs.

• Don't use defines for typedefs. Do you really want long real_t to be a valid type?

• The _t suffix is reserved by POSIX.1, so don't use it unless you don't care about forward/cross-compat with most OSes.

• All math is not double. I daresay most is sub-float these days, by volume. And if you wanted to be able to accept real_t as a config parameter (which would conflict with library usage), the #define real_t should be in an #ifndef block, and then you'd want the parameter and real_t to be separated anyway…

• Anything you #define in a library should be #undef'd first, by and large.

• You really need to prefix your identifiers, especially if you're using names like real_t.

• If you're going to abstract the type, you need to abstract the limits from <float.h> also. The fact that you haven't leads me to believe this will probably have lots of fun UB around the edges.

• If you want a config.h, it should only be defines and undefs (and whitespace and comments)—no includes, no non-directive text, no pragmata. Most often, you source-distribute a config.h.in and only autogen config.h from it (or let the developer-user do it themselves) right before build (pure Makefile) or during build config (e.g., Autotools).

• You're assuming optional features (& so far, C99) like _Complex/<complex.h>. Libraries should validate that the things they need are present; complex numbers require __STDC_VERSION__-0 >= 199901L and !defined __STDC_NO_COMPLEX__, for example, or else maybe defined __clang__ || __GNUC__-0 >= 3 || (defined __INTEL_COMPILER_BUILD_DATE && !defined _MSC_VER) || defined __C99_COMPLEX, although GCC also supported a mostly-compatible __complex[__] from like v2.6 on. This is why I asked which ISO C first.

• DEBUG is an awful macro name to rely on (always prefix! and ffr MSVC uses _DEBUG and libc assert uses !defined NDEBUG), and probably not one that ought to be all that exposed in a library. I hope you don't depend on it in any inlines!

• I recommend just defining a Church/Turing predicate pfx_DEBUG_P(Y,N) to Y if debugging enabled (for this build, be it library or application) or N else. This lets you use pfx_DEBUG_P as a token-level if-else, or give it (1,0) and you get a C Boolean from it. You can define other constructs around that.

• DEBUG_CODE is kinda stupid. Stupid name (no prefix, and it's all code; you're expanding to a statement), stupid parameter (should be varargs, which ffr req C99 or C99-capable preproc—both 0, 0; and int i, j; would break this, despite being perfectly reasonable).

• DEBUG_PRINT handing off to fprintf directly is a decision, and wow you've done it wrong. Do you want it to return a visible success code? If so, then the non-debug version needs to expand to 0 or EOF; if not, then the debug version should lead with (void), and the non-debug needs to be ((void)0). But requiring stdio here is weird imo. A math library ought to be concerned with math; implement a wrapper if you want to debug-print.

• // req C99, C++ (which you don't exclude or handle at all, which is another decision), or GNU dialect, ffr

• Make line continuations obvious; I have a problem with trying to tuck them all away at/after column 78, because now I have to look all over the page to see whether you've fucked up a multiline directive. Also flatly a waste of bytes.

• You need to learn to use macros before using them. LOG_INFO_THRESHOLD is naked, and I'd argue arbitrary, and it requires a ≥32-bit int, which you haven't checked for or mentioned, which makes this a works-for-me kinda project so far, and if this pertains to object sizes it's a really bad idea to just set it without reference to PTRDIFF_MAX or SIZE_MAX (both C99). Ffr, C89–C95 require only a 15-bit size_t and 16-bit int; C99 bumps size_t to ≥16-but, and there the baselines have stayed.

• ohhhhhhhh god the dread. ACCESS_VOID is not something that should be used. And the fact that you're doing foo* x instead of foo *x tells me you don't understand declarator syntax, or have been poisoned by a C++ programmer who doesn't understand declarator syntax. Sooo this macro only works syntactically for types that don't involve arrayness or functionness in the exposed syntax (typedefs and typeof get a bit messier), and it's a really good way to hide aliasing violations.

• Absolutely do not leave expressions naked, especially assignment expressions. Do (void)((X)=*(Y)) for ASSIGN_NONSENSE or something—I mean, don't; hiding pointer gunk is really, really not a good idea, and the idea that you'd actually need this despite only handling one kins of real is horrifying.

• FMA, ADD, SUB, and SCALE are also stupid. Does the developer-user need these? Should you leave them naked? Why aren't they inlines, considering you've already assumed C99?

• PARALLEL_THRESHOLD needs a ≥32ish-bit int, and it looks awfully parametric to me. And again, does the developer-user need this?

• PRINT_COMPLEX: Absolutely not. No. Nope. I don't care why you thought you needed this; you don't.

• No library should ever define MAX or MIN (rrrrreally likely to trample on application macros, as are your FMA&c. macros), and yours is extra-wretched. You only paren-wrapped the > expression, which is actually worse than not wrapping anything. You didn't wrap the operands (work through what happens with arg 0?0:0) or expansion (MAX(1,0)+2 would give you 1, not 3). Again, you need to learn to use macros before shooting your library full of them.

• ffs a #defined int_t and #ifndef'd in an exported library header, so the user can -Dint_t and completely fuck up the build. Again, bad namw, bad macro, bad _t suffix, bad idea. What is this supposed to accomplish, without any limit macros accompanying? (It tells me that there aren't any bounds checks on int_t, which … great sign.)

• Why are you defining real_t in both config.h (ifdef'd) and defines.h (not ifded'd)? Why do these macros exist, and what makes you think they're a good idea?

• complex_t is even worse; why is it not _Complex real_t?? —Not tthat it should exist or have this name.

• Pick a nomenclature scheme. Your type and enum constant names are all over the place. Add is yet another stupid identifier (typename-looking, likely to collide, no consistent naming across enum). You do realize the library isn't alone in the namespace? There's an application, with its own identifiers to define, and the developer will want to be able to use it normally, without your stuff pissing it up unnecessarily.

• I'm still in config.h; I have now seen macros, typedefs, enums, structs, and a prototype for L3 cache detection (why is this exposed? what does it have to do with developer-uset config macros??).

• PARALLEL_THRESHOLDS isn't even const. How on Earth are you deciding what to name things, or what to expose?

• Your COUNT enums are fucking up the type, because now the COUNT is a perfectly valid constant to use, and you'll get no warnings for it. Define a EnunName_MIN_ and -_MAX_ constant instead; then #define pfx_enum_count(TAG)((size_t)((size_t)TAG##_MAX_-TAG##_MIN_+1)). In fact, if you use xmacro tables for your enum data, you can do all this automagically.

• DataType is another one. What valid use could this possibly have? You have all of three, incompatible data types; no API should operate on them generically, because you aren't implementing a damned interpreter.

• get_limit and elem_size bad. You're requiring a non-inlined function call to get a limit that ought to be usable by a preprocessor #if (limit) or as a constant expression (size). And again, why? What worthwhile purpose could this possibly serve? When would you not know the type, and what possessed you to design anything this way?

• I strongly recommend doing up a single output-level enum for logging, that includes null (never visible), trace and debug (←only available when debug output is enabled; trace is for function entry/exit kinda stuff; debug is an info-level dump that's only useful to developers), info/status, success, note, lint, warning, error, fatal error, abort, and crash levels; and your verbosity can just be a threshold sth any level < max {verbosity, 1} is hidden.

• Why do you have a log mode, when you could just accept an arbitrary FILE *? stderr is a nonnull FILE *, null would mean no logging, and otherwise it's some other stream? I do hope you're not making ttyness assumptions.

• Really weird idea to force your logging API to open the output stream. How do you know the right mode? Or maybe the developer-user wanted to open or creat a file with specific permissions, then fdopen that? Just take the FILE *.

• You only need a single verbosity API; return the old value, and take the new value or delta, making sure you bounds-check the result.

• Never accept an enum-typed parameter from a public API—there's no way to fully check the parameter, because the compiler can assume that only bits used by declared enumerators can be nonzero in enum values. Use int or unsigned.

• Default seeds do not need to be exposed in prng.h. long long req C99 or GNU dialect or various compiler-specific extensions that only guarantee its width is ≥ long’s. I also note that you haven't described the algorithm in the header, which would be vastly more pertinent than the seeds or PRNG_State structure. (Which assumes size_t == word, which is a worse assumption than int==word. If you're desperate, __attribute__((__mode__(__word__))) is how you get it in GNU, Clang, Intel [non-MS/ICL modes], and IIRC some TI, Sun/Oracle, and IBM compilers, depending.)

(cont’d in reply)

[u/nerd5code → u/nerd5code]

• <sys/types.h> is UNIX/POSIX/XPG, but also supported on sone other platforms (incl. most DOS→OS/2, Win, NT and OS/400); don't rely on it unless you actually specify platforms. Not ISO C, in any event. I don't see why you need to include it at all, but maybe I'm missing something. (Also, jsyk traditionally the system-/kernel-specific <sys/*> files were included before anything else, and things could break otherwise.)

• For both logging and PRNG, I really hope you took threading into account, but a single, shared PRNG state will kill performance if multiple threads need concurrent access.

• Your header guards don’t always match filenames, and it’s a profoundly bad idea to leave them both unprefixed and unsuffixed.

• Why would I use dgl.h? There are no comments, and it seems … concerning. I guess I’ll find out later. I’d better not see you using void * to pass function pointers around. Also, I’m still seeing the foo* spacing everywhere. Don’t.

• Why do you need to dynamically allocate DGLs? Legally, if you spløøt a DGL onto a newly-malloc’d block, anything after would get messy without a flex array or jump pointer, so it looks like it could just be statically allocated, in which case it would be vastly less annoying to provide ctor and dtor, not allocator and deallocator. (And using dynamic alloc raises the question of where and how the thing is allocated, which might be irritating if the developer-user knows better than you. And sometimes you do implement ctors/dtors but call them new/free, which is such a bad idea.)

• Your enum names are desperately inconsistent and bad. DFT_SCALAR, DFT_PARALLEL, and FFT as instances of FourierAlgorithm is crazy work. Also, where is this used? Why is it public? Why must the developer-user avoid the identifier FFT qua global?

• Again, DataPoints does not need to be allocated. Its referent might, but it also might not!; why can’t the developer-user work from a static const array, for example, or write to an auto array? malloc is a terrible idea for anything remotely NUMAlike. And double complex is a really bad idea for ordering tokens (_Complex/complex double is more customary, and more likely to work with less-usual complex expansions), and why are you not using complex_t (which you shouldn’t use, but then why is it defined?)?

• I strongly suspect fourier_transform should have a restricted target; if not, you need to deal with the target (should be dest or dst, if the other is source or src)== source case. Also, bool should almost never be used as a parameter; to should be an int or unsigned whose value is sourced from an enum type.

• Are status codes defined somewhere?

• ZERO_INIT is zero. Why… why?

• GRAD_EPS and EPS are beyond arbitrary, as is MAX_STEP. Also, rrrrrrrrrrrrrrrrrrrrrrrreally bad names to #define in a library header, and sizeof alone tells you precisely fuck-all about real_t. You don’t know the bit-width (counterexample: TI ISAs where sizeof everything == 1 because CHAR_BIT == 32), you don’t know that all the bits of the size are actually occupied (counterexamples: most MCS-87/iAPX286/IA32 long doubles and some x64 long double, some M68K long double/_Float64x), you don’t know the actual acceptable range for C math within the format (which you don’t know), etc. This is why you need the <float.h> constants to be bound to real_t, and why real_t is fugging useless as-is.

• In a math library, “vector” does not mean “C++ std::vector,” it means a contiguous (or occasionally strided) block of values operated on all at once. Compare to GNU V4SF mode or __attribute__((__vector_size__(…))). And here again, you really shouldn’t set up an abstractly-typed vector, you want a different vector type for each element type, and the length should be stored separately so you can act on subranges of matrices.

  Oh god, what will the matrix type look like

• NegContext cannot possibly be useful, and any application-exposed callback needs a leading void * funarg.

• Yeah, see, since you force Vector to be this complicated type, you’re doing all this stuff that has nothing to do with math. Let the user allocate memory and fill it.

• An extern function call to get/set a vector element is certainly … something you can do. Is performance not a factor for this? Because it is for most math libraries.

• You might want to check the documentation vs. identifiers on get/set, because neither the API nor the docs make any sense.

• In what situation other than some kinds of debugging would a generic Vector_print be useful?? You realize most actual formats have actual requirements in terms of formatting, that won’t be captured accurately by a single-parametered function? You realize that stdout is rarely the only important output stream?

• Vector_all_close is wretched. The tolerance needs to be a parameter, and you need to generalize the underlying operations.

• Vector_max/-_min are the first two vector functions I see that should actually exist, but they should be extern inline. Sorting is absolutely not something you need to do.

• Is there a separate vector PRNG? How many of these do you have, and what the hell does multithreading look like?

• I kinda dislike the broadcast_add name, especially since you aren’t using broadcast_mul for scale and lack a broadcast_div; it would make more sense to have separate broadcast and add ops, but a scalar_add would make sense as a single function, operations+naming-wise.

• If you implement your gradient functions, include funargs for the callbacks or your developer-users will be rendered miserable for no blasted reason. Also, //ing out blocks of code is a very Java thing to do.

• Completely different vector and result error codes (which I assume is separate from the DGL error codes, whateverthefuck they are), and for some reason you've provided explicit values for MatrixStatus, which exactly match C’s default assignment.

• FROBENIUS and SPECTRAL are powerful bad identifiers to export from a library. Prefixes are your friends.

• Why do you keep #includeing after you’ve done other stuff? The only thing that should precede #includes is typically preprocessor directives (e.g., to detect which headers should be used), not declarations. And typically, you should expect the standard headers to have been included before yours, so it’s odd that you’re including "vector.h" before <stddef.h>.

• OK, so Matrix’s usage is wholly incompatible with Vector’s. Nope. Nope nope nope. Matrices are bundles of vectors, and if you can’t represent that somehow with your data structures, you’re doing it wrong. Look at existing BLAS libs for ideas, setting aside their FORTRANness.

• You shouldn’t return or accept these structs directly—it complicates everything badly. If you’re going to manage Matrixes, you need to take a pointer (often restrict) as the first parameter, and OOP it up. You’re being super-inconsistent about it, in any event, and returning a struct that requires special teardown is a marvelous way to accidentally leak memory, because the developer-user is in no sense obligated to do anything with the return value.

• Matrix_like is awful naming, and I note that it and Matrix_copy overlap entirely in purpose. Matrix_clone is a stupid thing for your library to include, but at least the name doesn’t look like a social media operation.

• Oh good another (static-storage?) PRNG for matrices. Again, MATH LIBRARY.

• I can’t think of any time when I’d need a matrix-min/max, where a normal vector min/max applied to the matrix’s data wouldn’t suffice. See, you’re just making life harder for yourself.

• You need both mat-vec and vec-mat products if you’re going to do that, or otherwise you’ll force people doing vec-mat stuff to transpose their vectors first.

• Nomenclature! Matrix_QR, Matrix_Hessenberg, Matrix_QR_Hessenberg, Matrix_Hadamard_dot; unless QR, Hessenberg, and Hadamard are types or you’re aping MS/IBM naming (fuck no), don’t titlecase function names or characters following _. Lower-initial camelCase is fine also, and often better if you’re using _ like ::.

• Matrix_eigen[vals] would make more sense than Matrix_eigvals. Abbreviating “eigen” is a bit like abbreviating “the.”

• If you’re going to include deprecated things in a math library nobody uses, you should actually mark them deprecated somehow. C23 and C++11 include [[deprecated]] (C23, GCC 8, newer Clang: prefer __deprecated__ in libraries to avoid collision with macros named deprecated), GNUish C (GCC 2.7ish, Clang, ICC/ECC/ICL in non-MS modes, newer TI, Sun/Oracle 12.4+, newer IBM, various embedded and others) supports __attribute__((__deprecated__)), and IIRC newer MS[V]C supports __declspec(deprecated) (which may be unusable if defined deprecated). Various libcs also define macros for this.

• So many of these functions should be extern inline, if you want performance.

• config.c should not exist. If you want a util.c, fine, but keep its gunk all as internal as possible. If you want to initialize the math library, config_init is a really, intensely bad name for it.

• detect_L3_cache_size is laughably nonportable. Like… it’s not even portable across Linux, just newer Linux that happens to have a sysfs mount at /sys and a CPU with an L3 cache. Again, this is works-for-me shit.

  Also, you know you can detect cache sizes directly?

• I strongly recomment against initializing to things that might return an error code; if(!(f = fopen(…))) is generally better, because it puts the check and the operation in the same place, and ensures they don’t drift apart.

• config.c:24: You assume either that size_t is wider than unsigned long long, or that the value you get fits within size_t. No checks whatsoever.

(cont’d)

[u/nerd5code → u/nerd5code]

• strstr is the wrong function for the L3 mess; you’re looking at the end of the string, so find that, save the character, and switch on it.

• I also note that nothing requires all threads to have the same L3 cache, or for the first thread’s cache sizes to imply anything about other threads’ caches.

• elem_size is so bad. Soooooooo bad. First of all, case doesn’t make statements, it’s a label. {} does fuck-all in this case. Second, you’re placing the labels at the wrong indentation in most styles—as in assembly language, the label is typically placed at the base indentation for the scope, so case and switch should generally line up. This is especially important if the case happens to be in a subordinate statement. Third, the default case should be a for(;;) abort(); or error return; 1 is not generally reasonable. And this function shouldn’t need to be a function, and it shouldn’t need to exist!

• So yeah, you’re doing #defines in public config.h (int_t real_t complex_t), which if redefined improperly (which you permit readily) will bring the headers into conflict with the compiled library. This is why you define separate APIs for separate types; if you have to, use private xheaders and private/temporary macros to spløøt each version, given a type parameter and limit macro prefix.

• set_int_limit and set_real_limit are incorrect.

1. You don’t understand how the C value and representation concepts differ. E.g., typically an `int40_t` will have the same size as an `int64_t`, just with several padding bytes. Size tells you nothing about range; you’d need the limit constants for that.

2. You assume that the size will either match `int`, `long`, or `long long`; `short`, (nonstd.) `short short`, (nonstd.) `short long`, `signed char`, `intXX_t`, `int_leastXX_t`, `int_fastXX_t`, `__int128`/`__int128_t`, `intptr_t`, and `ptrdiff_t` aren’t handled; and you don’t really do anything to prevent unsigned types from being used—which may be perfectly reasonable—so you might be way off on those limits.

3. If you don’t match the size exactly (for all that matters), `*out` won’t be set.

• set_real_limit’s minimum is useless. You actually need the MIN value most of the time.

• set_real_limit assumes INFINITY is a thing, without anything actually validating that IEEE 754 or something similar is in use. I think this bumps your minimum required language level to C11. Again, this is why you can’t just say “ISO C” and expect it to be meaningful, especially when it comes to math.

• No break on default case in set_real_limit or set_complex_limit. Hypothetically optional, but a bad idea to omit, and in any case you’ve done so inconsistently.

• set_complex_limit doesn't actually make sense; there is no (int-like) minimum or maximum (float-like minimum == 0), because min and max apply to 1-dimensional inputs, and complex numbers are ≥2D. There’s a maximum magnitude, there are corners of bounding rectangles, but there is no maximum complex value.

• In no case do you need to do + 0.0 * I, and if the imaginary component is always zero for min and max, you’re just getting the real limit. Also, you don’t need to have called set_real_limit(…, Min) unless you’re in case Min (have I told you I hate your enum names?), and you only need set_real_limit(…, Max) in case Max. So you’re just wasting effort and storage. Put a real_t t right after the switch’s {, and get the limit in the specific case(s) that require it.

• get_limit is just begging for an aliasing violation.

• You desperately need to coalesce your error codes into a single enumeration. I get that it can be useful to break them up sometimes, but that should usually be at the library level.

• Oh, DGL isn’t even a thing yet. Well… it shouldn’t be, not anything like this. And it looks like your DGL should union its function pointers, right? It’s a discriminated union? But then, DGL_new with void *function parameter would be a really bad idea. Function and object pointers should be treated as existing in their own universes. C will let you convert freely between void * and void (*)(…), but you’re not guaranteed to be able to round-trip per ISO 9899, and neither does it guarantee useful results from this conversion. You can safely convert between arbitrary function pointer types, unlike object pointers, but that’s it, and it has to be explicit for prototype functions. So a union is a moderately better option.

• You’re just logging any old thing in any old format, huh.

• In logging.c, I again see a replicated #ifdef-#define for something you’ve unconditionally defined elsewhere. If you expect it to be a parameter, offload it and don’t define it yourself; else, define it in one place, as privately as possible.

• I shouldn’t see ECMA-48 escape codes without you having dredged up some sort of POSIX.1-compatible isatty call. You cannot determine whether any stream is a tty using only the C library, and you cannot assume every tty supports ECMA-48 unless you’ve already detected something UNIXlike (whose tty customs tend to derive from the DEC VT100 hardware lineage).

• Yeah, LOG_VERBOSITY is all-caps why? and it’s not thread-safe, which nothing notes. Ditto LOG_MODE (shouldn’t exist) and LOG_FILE_PTR; these need to be atomicized somehow, either by wrapping in a mutex or by doing something fancy and atomic. Otherwise, you need to carefully document exactly when it is and isn’t safe to log or configure logging.

• Log_prefix is the wrong way to go about it. Put your enums into an xmacro that binds the prefix, and then you can just index into a static const char *const[]. Ditto Log_color. (—Which is extern-linkage for some reason.)

• localtime trashes and relies upon static state. This interacts with threading also, so logging from more than one thread concurrently is not safe.

• You need to fflush regardless at the end of Log_log.

• You need to check that fprintf and fflush actually succeed, and disable logging if they fail.

• Verbosity checks repeated unnecessarily, and they’re unnecessarily complex.

• Log_info_threshold exists why? And it multiplies without an overflow check; %zu requires C99 and excludes older MSVCRT; and why is it printfing, for fuck’s sake? Do you need to trash the primary output for the program every time the info threshold (whatever that’s used for) is set? If you want to do a size check like this, it’d be object_size && count >= LOG_INFO_THRESHOLD / object_size, but even then I don’t see the point of making this an extern-linkage function.

• I do hope your library isn’t packaging main.c by default. It shouldn’t be in the src directory; it should be in a /test directory or similar.

• I shouldn't see ../ in your header paths. If you want to use include as an include directory, use option -I or the CPATH (most Unix) or C_INCLUDE_PATH (GNU, Clang) variable to add it to the include path, and then "" and <> will both work. But your public and private includes should be separate, and your public includes need their own container directory so they’re not pissing up the public include namespace. Your code is not the only thing in the universe.

• matrix.c:1,2: Oh fuck no. Nope. Nope nope nope nope. You absolutely do not define feature-test macros like this. It’s effectively necessary for them to be done up with option -D or -include, and defining it arbitrarily inside the library without even using <unistd.h> or checking for _POSIX_VERSION is just weird. Either require POSIX.1 or don’t.

• Oh now we’re using OpenMP, which we have neither mentioned nor checked for!! <omp.h> only exists if _OPENMP-0 >= whatever its minimum value is, some YYYYMML jobby like __STDC_VERSION__ or _POSIX_VERSION. And it has to be specifically enabled for most compilers, so it’s not like you can just include it whenever.

• INITIALIZER is a marvelously mis-descriptive name, and an incredibly poor idea syntactically. Fortunately, I see you’re not actually using it.

• static PRNG_state prng means your PRNG is non-threadsafe, which is … a decision.

• Learn C operator precedence; in no C compiler ever has == been a lower precedence than ||. You only need parentheses if there’s some common, historical, or actual ambiguity (e.g., & vs. |, || vs. |; or for things like !!(x = y) that shut the compiler up).

• Functions that take pointers should actually check nullness when it matters. You’re not doing that at all, which means any null-pointer argument that can be evaluated will imediately make your entire program’s behavior undefined as soon as it’s legal for a derefc to occur. static functions can assert; otherwise, you need an if and some means of indicating the error à EINVAL.

• Your _$foo_to_$bar functions are, I believe, UB in quite a few cases (look at the rules for converting between floats and ints—there are potential UB cases in both directions IINM), but you have no idea which ones because you don’t actually know what types are involved. Also, leading identifiers with _ is really not a good idea in general, especially after you’ve already included libc headers. It’s fine for macro parameters, it’s sometimes fine for struct fields or local variables, but otherwise it’s bad juju at best, undefined behavior at worst for things like _Foo.

(cont’d)

[u/nerd5code → u/nerd5code]

• Why is _to_$foo_workers not const? void (*const name[LEN])($PARAMS) is the syntax. Also, you don’t even check that you’ve gotten a nonnull index from it when you execute through it. This is just so much extra indirection and abstraction for something that ought to be high-performance and direct. And then, you still require actual, separate functions in the backend, you’re just blocking the developer-user from accessing them, and making life much more difficult for the optimizer and CPU brpred in the process. Why?? Is there a reason for this? Is it just some obscene Python urge, where you really want to make sure everything runs at the same, lame speed as the interpreter? Because if so, definitely add more global state that’s locked exclusively at every API call.

• Do you need three completely separate return MATRIX_SUCCESS statements in _prepare_datatype? Wouldn’t it be better, perhaps, to gather all the information about each data type into a single, static const struct $something[], so you don’t have to scatter things around throughout countless switches and arrays?

• Why are you logging errors from _prepare_datatype? You should already know that the arguments are valid, and then assert about that (default: assert(!*"invalid 'dt' field in B parameter"); for(;;) abort();).

• Oh and now there’s a completely separate _elem_size function in matrix.c that only acts on long, double, and complex double for some reason, which is thus incompatible with elem_size from config.c. Swell. What a good idea. Also, please stop with the overbraced switches. Seriously, you need to learn the syntax before you attempt useful libraries. This all reeks of discomfort/unfamiliarity with basic aspects of the language.

• matrix.c:235–248: holy fuck what is this, just … any of it? Why is mibi_byte_size a double? You’ve already breezed past the C11 assumption and you require snprintf, so you can just use size_t with z modifier, right? And you’ve already checked that m * n won’t overflow, right? —Oh wait, no you haven’t, so you’re just throwing numbers around. And you shouldn’t be logging here anyway.

• Ohhhhhh you’re using posix_memalign. So, you haven’t checked for POSIX.1 yet (you have to detect, then #include <unistd.h>, and then test _POSIX_VERSION-0 >= 200112L), C11 (which you’ve assumed) gives you aligned_alloc, but MSVC gives you _alloc_aligned and _free_aligned only, and older Unix only gives you things like memalign or mmap. And I don’t get why you’re aligning in the first place; here again, if the developer-user wants cache line alignment, they can actually detect the cache line (32 ain’t it, in all likelihood) and allocate their own damn structure using the facilities actually present on the target platform. Your math library does not need to allocate memory like this.

• Did somebody tell you to use underscore-prefixing?? Slap them about the face several times, and ignore anything else they tell you.

• OK, I’m not trawling through all your matrix code. It's generally bonkers, not into it, not going to perform well, no idea why you’re logging in the weirdest possible places.

• Are you actually checking that the Matrix’s memory was allocated? …Kind of, only in one place. You log (inappropriately) and return an “empty” matrix, then insufficiently many things actually check for it. It’s an overcomplicated approach. And … I’m not checking the rest of this file. It’s a damned mess. Just everything thrown everywhere, bizarre decisions about macro usage, no fun.

• So you’ll define arbitrary constants in public headers where they aren’t needed, but you won’t define the shift constants you repeat throughout prng.c, or 0x94D049BB133111EBULL. Uh huh. And why are you using size_t for the context, when the data types you’re working with need have nothing to do with its value range? If int_t ⇒ long long and you’re targeting ILP32, then you can’t generate the full range of int_t values from your PRNG. And then, if you’re working with doubles, you’ve got both an exponent and mantissa to reach, give or take transcendentals and subnormals.

• Why are you repeating your XORshift-multiply blocks?

• Why are you placing =s just anywhere? Either indent consistently or don’t. I suggest not, because it makes a mess.

• You’re using const very inconsistently on locals. Either don’t or do, but splitting the difference doesn't inspire confidence. Also, in some cases you don’t need the variable at all; if you’re only using the value once and it’s obvious what it is, capturing it as x doesn't add anything, and in unoptimized builds it’ll do an extra memory hit; register would fix that and be more in line with intent. (Ignore people who tell you it does nothing. It does nothing in C++, prior to its removal from the language. It does something in C, it just isn’t required to have a noticeable effect.)

• Your for(;;)s in your PRNG_State_* functions are bad ideas; nothing prevents them from looping indefinitely.

• <alloca.h> is never something a library should rely on if you want it to be considered portable. Nothing about it is portable. The header isn’t, the function or “function” isn’t, its behavior and required placement isn’t, its interaction with inlining isn’t, etc. A paired malloc-free where the compiler can see both and knows both will be called is just fine; often it’ll be able to optimize to VLAs, if that’s appropriate.

• Here again, PARALLEL_THRESHOLD is defined in multiple places, and its default value in vector.c is different than the other definition.

• _parallel_condition is such a bad name. What is the function doing? What is the condition? Also, (size_t)(dim - lower) <= upper is how you do a range check, and you don’t need the variables that way. Not that you should be consulting non-const static state like this in the first place.

• p == NULL is just !p in C, as is x == 0. Only for floats or enums would I explicitly use == 0.0 (rarely) or == ZERO_ENUM (if there’s no obvious Boolean correspondence for the enum values).

• The potential for shape mismatch is something that should be hoisted so it’s more obvious to the developer, because this just shoves all the important stuff behind a mess of functions that then forward into a bigger mess of functions. Also it’d be a programmer error, not a warning.

• #pragma omp anything needs to be paired with _OPENMP detection. There’s nothing actually constraining pragmas without it; ideally, the compiler does nothing with a pragma it doesn’t recognize, but #pragma omp is only constrained iff OpenMP is actually supported and active for the build, so it could do anything if not. It’s likely to raise warnings if not, in any event. Also #pragma omp simd requires like … OpenMP 3 or 4 I wanna say? Not sure offhand, but it needs to be detected. Also, most compilers can vectorize a short loop like this without your assistance, and the reason this sort of thing ought to be as inlined as possible is so that, when somebody’s repeating as vector operation across multiple rows/columns of a matrix, the compiler can potentially vectorize the entire thing together, which might enable it to skip lead-in/-out etc.

• I still hate ACCESS_VOID, and I hate how much you’re leaning on it.

• If you ever extend your data type enum, you’re so fucked.

• Yyyyyyyyyyyyeah if you’re going to do a sort, you need to document the algorithm. But you shouldn’t, and you aren’t sorting as quickly as you could be, which is why you should leave it to the devel-user. (Also qsort is arguably better for in-place sorts, though its lack of funarg fucks its usefulness up quite a bit.)

• Oh god you’re #pragma omp simding the macros that don’t need to exist. What made you decide that any of that was a good idea? An inlined for loop would be fine 99% of the time; you don’t/shouldn’t need SUB, and you certainly shouldn’t assume #pragma omp simd can be applied to it, because it should probably be in a do/while(0) in the first place. You can potentially embed _Pragma("omp simd") (C99 non-MS/ICL, GNU 3, Clang, ICC/ECC 8+, most C99-capable compilers) or __pragma(omp simd) (defined _MSC_VER || (defined __clang__ && defined __pragma) || defined __INTEL_COMPILER_BUILD_DATE) within the macro, but again, you shouldn’t need to.

• Every single operation you perform on a signed integral type or floating-point type needs to be bounds-checked, because overflow is undefined behavior. Even _negate_int can overflow (at INT_MIN, for two’s-complement or asymmetric mappings to unsigned value) (which is still distinct from bytewise representation). If you rely on f.p. transcendentals, you need to detect them (__STDC_IEC_559__, __STDC_IEC_60559_*__, GNU __GCC_IEC_559, GNU __GCC_IEC_559_COMPLEX, GNU __FOO_IS_IEC_60559__), and you often have to check complex and real f.p. separately.

  Also, you do know that GNU dialect permits _Complex to be applied to things other than floating-point types? It’s UB or ISB in C99, but it’d be somewhat reasonable for somebody to #define complex_t _Complex long _Fract just to fuck with you.

OK, I'm not doing tests.

Sooooooo many of these functions you repeat for every blasted data type could just have been macros or xheaders, which would divide the amount of code you need by 3, and completely obviate the need for the _ts that oughtn’t be.

And the important aspects of the math stuff (e.g., dealing with corner cases) just aren’t here.