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));
}
u64
typedef? \$\endgroup\$u64
meaninguint64_t
. \$\endgroup\$