5
\$\begingroup\$

While working on Project Euler, I have come across a few questions involving usage of the little omega function, big omega function, and Möbius function.

My prior code wrote factors into an array and stepped over each factor to check values for the little and big omega functions, then counted the array's length and calculated the Liouville function. This seemed very inefficient to me as a number could have a lot of very small prime factors which would cause ω(x) to not equal Ω(x), then a large prime factor (e.g. 2 * 2 * 8191), but every prime factor would have to be calculated first.

My code uses an "iterator"; since C does not have iterators, I have a state which stores all information necessary for a trial division algorithm. The function which implements the algorithm will return a truthy value if a factor exists, but will return a falsy value if no factor exists. This system is implemented in a way where the function can be called in a loop to generate every factor present. I see several positives of this, particularly that the factorization can be stopped early (and picked up, if desired) as well as that no heap- or stack-allocated memory is required. If the function is being called many, many times, speed may be a concern, but no speed concerns would be present for a factorization algorithm.

Currently, my code implements trial division factorization. I know that trial division is very inefficient, but a sophisticated wheel factorization or Pollard's rho algorithm would not focus on the crux of this question.

Is this solution common practice? Is it being implemented in a commonly-implemented manner? Are there any other improvements that should be done to my code?

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

// these integer types would be included by a local header
#include <stdint.h>
typedef uint32_t u32;
typedef uint64_t u64;
typedef int_fast8_t fs8;
typedef uint_fast8_t fu8;

typedef struct {
  u64 on, original;
  u32 currentFactor;
} FactorState;

FactorState Factor_new(
  u64 x
) {
  if (x == 0) {
    fprintf(stderr, "cannot factor zero!");
    abort();
  }

  return (FactorState){ x, x, 2 };
}

fu8 Factor_factorAll(
  FactorState *const state,
  u64 *const factorOut
) {
  for (
    u64 cf;
    cf = state->currentFactor, state->on > 1 && cf * cf <= state->original;
    ++state->currentFactor
  )
    if (state->on % cf == 0) {
      state->on /= cf;
      *factorOut = cf;
      return 1;
    }
    
  if (state->on > 1) {
    *factorOut = state->on;
    state->on = 1;
    return 1;
  }

  return 0;
}

fu8 Factor_factorDistinct(
  FactorState *const state,
  u64 *const factorOut
) {
  for (
    u64 cf;
    cf = state->currentFactor, state->on > 1 && cf * cf <= state->original;
    ++state->currentFactor
  )
    if (state->on % cf == 0) {
      do { state->on /= cf; } while (state->on % cf == 0);
      *factorOut = cf;
      return 1;
    }
    
  if (state->on > 1) {
    *factorOut = state->on;
    state->on = 1;
    return 1;
  }

  return 0;
}

fs8 mobius(
  u64 x
) {
  if (x == 1)
    return 1;
  
  FactorState all = Factor_new(x), distinct = Factor_new(x);
  u64 allFactor, distinctFactor;
  size_t allCount = 0;
  
  do {
    fu8 allRes = Factor_factorAll(&all, &allFactor),
        distinctRes = Factor_factorDistinct(&distinct, &distinctFactor);
    
    if ((allRes || distinctRes) && !(allRes && distinctRes))
      return 0;
    if (!allRes && !distinctRes)
      return (allCount & 1) ? -1 : +1;
    
    if (allFactor != distinctFactor)
      return 0;
    
    ++allCount;
  } while (1);
}

int main(
  int argc,
  char **argv
) {
  for (u64 v = 1; v < 20; ++v)
    printf("%lu = %d\n", v, mobius(v));
}
\$\endgroup\$
3
  • 1
    \$\begingroup\$ Why do you feel the need to make your own aliases for the standard C types? It makes things harder to read. Furthermore, could you imagine if you had an application using five libraries, and every library wanted to make its own u64 typedef? \$\endgroup\$
    – JohnScott
    Commented Nov 22, 2022 at 23:40
  • \$\begingroup\$ @JohnScott I used those integer typedefs because it's easier for me to type and visually process. However, I would imagine that the five libraries would be able to agree on u64 meaning uint64_t. \$\endgroup\$ Commented Nov 23, 2022 at 1:10
  • 1
    \$\begingroup\$ "because it's easier for me to type and visually process." --> better to code with clarity for others as the goal. \$\endgroup\$
    – chux
    Commented Nov 26, 2022 at 23:03

1 Answer 1

1
\$\begingroup\$

Is this solution common practice?

Scant Comments

As code lack code comments, what is truly happening is harder to discern than needed. Hard to tell what is really happing let alone determine if it is common practice. Code golf (minimal test) should not be the goal - clarity is.


Is it being implemented in a commonly-implemented manner?

Shorthand for fixed width types

Just stick with uint64_t, etc., rather than make up your own shorter versions.

bool

Rather than invent a Boolean-like type fu8, use bool.


Are there any other improvements that should be done to my code?

Overflow risk

cf * cf in cf * cf <= state->original may overflow when state->original is a large prime near UINT64_MAX.

cf <= state->original/cf accomplishes the same test without overflow concerns.

Matching specifier

Code uses "%lu" with u64 (aka uint64_t). C does not specify that match. Use "%" PRIu64 from <inttypes.h> to portability match uint64_t.

Similar issue for printing a fs8 (aka int_fast8_t) that should use "%" PRIdFAST8.

Why size_t

allCount is used simple as a bool. Consider that type instead.

\$\endgroup\$
3
  • \$\begingroup\$ On many platforms, division is slower than multiplication, so there might be a case for a smarter function to determine whether cf * cf <= state->original using division when necessary and multiplication otherwise. \$\endgroup\$ Commented Dec 8, 2022 at 17:54
  • 1
    \$\begingroup\$ @TobySpeight I suspect code could be reworked to only need a a/b and nearby a%b. As good compilers perform these 2 efficiently as one, a b*b <= a would slow things down. IAC, correct functionality first over the range of possible inputs, performance 2nd. \$\endgroup\$
    – chux
    Commented Dec 8, 2022 at 18:01
  • 1
    \$\begingroup\$ Ah yes, I didn't look at the surrounding context! \$\endgroup\$ Commented Dec 8, 2022 at 18:09

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.