]]> https://programmingpraxis.com/2019/03/15/the-henry-used-car-dealership/#comment-179227 Mon, 19 Aug 2024 13:09:24 +0000 http://programmingpraxis.com/?p=12640#comment-179227 In response to @matthew, the book contains such gems as “In Chapter 2, you learned about many features of program modules,
also called methods.” It should come as no surprise that it makes HEAVY use of flowcharts, rather than say Nassi-Shneiderman charts. By the 1970s, we had learned a lot about program design and logic. They’re in the title of the book. But the book does not mention ‘assertion’, ‘precondition’, ‘postcondition’, (data structure) ‘invariant’ or (loop) ‘invariant’. ‘Termination’ is absent from the index, so good luck searching for advice on how to make sure your loops and recursions stop.

As a free speech advocate, I am pleased that such a book can find a publisher in the 21st century. But it makes me very sad that anyone would try to learn programming from such a book.

RosettaCode reasonably recently added this same task (https://rosettacode.org/wiki/Legendre_prime_counting_function) to which I added “Comments to the Task”. It turns out that there is a very simple optimizaton of stopping ” “splitting” the φ(x, a) “tree” when ‘x’ is zero” that stops the exponential growth of the recursion tree and will give much better performance even without memoization. For instance, the non-memoizing Haskell version I contributed there takes only about 273 milliseconds to count the primes to 1e9: https://rosettacode.org/wiki/Legendre_prime_counting_function#Non_Memoized_Haskell_Versions.

There are other optimizations covered there, such as reducing recursion by using a “bottom-up” rather than “top-down” recursion and a look-up table for small values of “phi”, but the major advance is to use a scheme that doesn’t use recursion at all but partial sieving with linear arrays as implemented in several languages there, including Haskell and Nim. Using this technique, the primes to 1e9 can be calculated in just a millisecond or two, to 1e11 in about 10 milliseconds, and to the quite large range of 1e16 in minutes. Note that not using recursion is a form of “LMO” as applied to the Meissel (rather than the Legendre) technique but converted back to apply to the Legendre mathematics, and that this technique is no longer linear time complexity with range (or exponential without memoization or the short-circuiting optimization). However, using the Legendre partial-sieving technique still has its limits as it uses considerably more memory (proportional to the square root of the counting range rather than the cube root) and therefore is limited to the range of about 1e16 on many modern desktop computers.

It is also more difficult to adapt even the Legendre partial sieving technique to multi-threading than the “LMO” and later techniques as applied to the Meissel mathematics, which have been used in the Kim Walich “primecount” (https://github.com/kimwalisch/primecount) implementations to count the number of primes to very large ranges of 1e29 and likely advancing…

It is not too hard to adapt the non-threading partial-sieving “LMO” technique to even use with JavaScript (https://stackoverflow.com/questions/43432760/looking-for-a-fast-prime-counting-function/70049343#70049343) which allows one to count the number of primes up to almost 1e16 in a minute or two…

float invlog(float x) { return 1.0 / log(x); }

float cube(float x) { return pow(x, 3); }

float trapezoid(float (*f)(float), float a, float b, int steps) {
if (steps <= 0) { printf(“Invalid number of intervals\n”); return 0; }

float delta = fabs(b - a) / (float)steps;
float integral = delta * (f(a) + f(b)) / 2;

for (int i = 1; i < steps; ++i) {
    float x_n = a + i * delta;
    integral += delta * f(x_n);
}

return integral;

}

float simpson(float (*f)(float), float a, float b, int steps) {
if (steps <= 0 || steps % 2 != 0) { printf(“Invalid steps. Must be even and positive\n”); }

float delta = fabs(b - a) / (float)steps;
float integral = (delta / 3) * (f(a) + f(b));

for (int i = 1; i < steps; ++i) {
    float x_n = a + i * delta;
    integral += (delta / 3) * ((i % 2 != 0) * 4 * f(x_n) + (i % 2 == 0) * 2 * f(x_n));
}

return integral;

}

float quadrature(float (*f)(float), float a, float b) {
float eps = 1e-7;
float integral = 0;
float midpoint = fabs(b – a) / 2.0;
float error_margin = fabs(simpson(f, a, midpoint, 2) + simpson(f, midpoint, b, 2) – simpson(f, a, b, 2));

if (error_margin > eps) {
    integral += quadrature(f, a, midpoint);
    integral += quadrature(f, midpoint, b);
} else {
    integral += simpson(f, a, b, 2);
}
return integral;

}

int main() {
float (invlog_ptr)(float) = &invlog;
float (cube_ptr)(float) = &cube;
float limit = pow(10, 21);

float integral = quadrature(invlog_ptr, 2, limit);
float simp = simpson(cube_ptr, 0, 1, 10000);
float traps = trapezoid(cube_ptr, 0, 1, 10000);

printf("%f\n", integral);
printf("%f\n", simp);
printf("%f\n", traps);

invlog_ptr = NULL;
cube_ptr = NULL;

return 0;

}

https://www.dreamsongs.com/10ideas.html contains John McCarthy’s (himself!) Common Lisp solution using the “gopher” function, which surgically rotates each tree in turn to the right until its left child is a leaf.

Similarly, the Haskell solution above is problematic. The implied constraint in the problem is that the leaves enumeration be done in O(n) i.e. without any unnecessary overhead, but that flatten implementation will be quadratic for trees which are left-degenerate. Instead, and similar to the “gopher” idea, it should be defined as

    flatten t = go t []  where
        go (Leaf x) xs = [x] ++ xs
        go (Node lt rt) xs = go lt (go rt xs)

which can be seen as doing the said rotations implicitly, if we squint hard enough, since just as a gopher it burrows down to the leftmost leaf first, thanks to Haskell’s lazy evaluation.

What’s not obvious from that is what the standard deviations are.
A simulation can give that quite simply.
Monday => (count: 1000000 min: 3 max: 112 mean: 9.996582000000235 sd: 7.614711570447746 skew: 2.029387841503218)
Tuesday => (count: 1000000 min: 3 max: 71 mean: 7.994588000000178 sd: 4.892717307278136 skew: 1.833047061117329)

I tackled problem 2 by brute force.
(1000 to: 3000) select: [:year | (year // 100) + (year \ 100) = (year // 10 \ 100)]