Reflections on Valentin's 2023 Pi day special

[EdS2]

Once again Valentin has posed an interesting programming challenge based on an interesting mathematical finding.
[VA] SRC #013 - Pi Day 2023 Special

The link to Pi this time is in the counting of square-free integers (I've a feeling we called them quadratfrei in my day - anyhow, numbers which have repeated prime factors, as compared to those which don't.)

I did a little investigating of this, in BBC Basic, one of my preferred languages, and of course this is outside Valentin's rules of engagement.

What's different in my case, then, is having none of the augmented 71B's augmentations: no factoring or prime number functions. It's certain to come out slower.

My related challenges, then, would be
- to respond to Valentin's challenge using an unaugmented 71B, or a similarly ordinary Basic
- to find values of N which incidentally deliver unfairly close values of Pi

In particular, I'd like to see (or write) a program which does the accurate counting by inclusion-exclusion, rather than trial division by prime squares. It should be much faster, but will be more complicated.

I'd also quite like to port V's small-memory solution to an unaugmented Basic.

[Valentin Albillo]

.
Hi, EdS2,

(03-22-2023 09:11 AM)EdS2 Wrote:  Once again Valentin has posed an interesting programming challenge based on an interesting mathematical finding.

You should work on diversifying your adjectives.

Quote:The link to Pi this time is in the counting of square-free integers (I've a feeling we called them quadratfrei in my day

No need to use a German term in an English-language forum when there's an exact English equivalent.

Quote:I did a little investigating of this, in BBC Basic, one of my preferred languages, and of course this is outside Valentin's rules of engagement.

Not so. Quoting my rules:

"Note: No hard rules so no need for a parallel thread, post here whatever you want as long as it's on topic and NO CODE PANELS, but I'd appreciate it if you'd use vintage HP calcs (physical/virtual), [...]"

and "I'd appreciate it if" is not the same as "I mandate that" so you could've posted your BBC BASIC code in the main thread. That is, if you actually wrote a solution.

Quote:What's different in my case, then, is having none of the augmented 71B's augmentations [...]

You should really work on diversifying your vocabulary.

Quote:My related challenges, then, would be
- to respond to Valentin's challenge using an unaugmented 71B, or a similarly ordinary Basic

Right now I'm going to have dinner but once it's over I'll concoct and post some code for you.

Quote:- to find values of N which incidentally deliver unfairly close values of Pi

Please define "unfairly close". Too vague as it stands.

Quote:In particular, I'd like to see (or write) a program which does the accurate counting by inclusion-exclusion, rather than trial division by prime squares.

It's fairly clear to me that you don't realize that this is exactly what the Möbius function accomplishes: it provides the required minus sign for the terms to be excluded, and the plus sign for the terms to be included. In other words, what both my solutions currently compute, so there you are.

Quote:I'd also quite like to port V's small-memory solution to an unaugmented Basic.

By "unaugmented Basic" you mean a non-expandable BASIC ? Please do. My code just about to be posted might be of help.

Thanks for your interest.
V.

[EdS2]

.

What might I mean by an unfairly close value of Pi? Well, perhaps the successive record-breakers which get closer to Pi than their predecessors, or those which get much closer.

Have I written some BBC Basic? Yes, but only exploratory, not presently something to share. I'll post it if and when it feels ready. I'd always be interested in anyone else's efforts, of course, including yours, in an unaugmented Basic of your choice. (I think we still hope for a response from Eric.)

Indeed, you're right, I've been thinking of square free numbers, more than thinking of the Möbius function, and had missed the usefulness of its varying sign. I'd even spotted the same sum which runs up to root N, but the penny hadn't dropped. You made good use of both!

(In passing, I noted the term quadratfrei precisely because it wasn't in a German context, but in an English-language context. No need to tell me that there's no need to choose the word I choose to choose. I make the choicest choices, excepting only the exceptions.)

[Valentin Albillo]

.
Hi,

(03-23-2023 10:43 PM)EdS2 Wrote:  What might I mean by an unfairly close value of Pi? Well, perhaps the successive record-breakers which get closer to Pi than their predecessors, or those which get much closer.

Still too vague, nothing operative here. Though perhaps you're referring to results such as

N = 33,000,000, Count = 20,061,593, Error = 0.00000011

and

N = 100,000,000, Count = 60,792,694, Error = 0.00000042

where we have 33 million ... An innocuous number in and of itself, though practically speaking, it could be considered an "unfairly close" approximation, i.e., ~ 1/3^rd the count but ~ 4x smaller error than those for the much bigger 100 million. Anyway, "unfair" surely isn't the proper word.

Quote:Have I written some BBC Basic? Yes, but only exploratory, not presently something to share. I'll post it if and when it feels ready. I'd always be interested in anyone else's efforts, of course, including yours, in an unaugmented Basic of your choice.

Thanks, very kind of you.

As for that mantra of yours, "augmented Basic, unaugmented Basic", it really grates. The "unaugmented" BASIC dialects of old are just the normal, plain-vanilla, run-of-the mill ones, which are the vast majority and so they need no exotic qualifier, they were the normal.

And the "augmented" ones, in the case of HP-71B BASIC and other technical HP BASIC dialects in general, they are just enhanced dialects, usually with subprogram and matrix capabilities, etc., among many others, and they catered for extensions (not "augmentations", as you called them) to the keyword set in the form of new keywords written in assembler, which add up seamlessly to the original instruction set.

Quote: (I think we still hope for a response from Eric.)

I'm not hoping for his reply. Eric's got too many more important matters in his always busy hands so he knows better than to allocate precious time for this, though he kindly suggested the possibility out of politeness.

Quote:Indeed, you're right, I've been thinking of square free numbers, more than thinking of the Möbius function, and had missed the usefulness of its varying sign. I'd even spotted the same sum which runs up to root N, but the penny hadn't dropped. You made good use of both!

Thanks. I take it that you now understand that the useful formula S(n) itself is what you come up with if you go on and apply the inclusion-exclusion principle, because it actually embodies it.

Quote:(In passing, I noted the term quadratfrei precisely because it wasn't in a German context, but in an English-language context. [...] I make the choicest choices, excepting only the exceptions.)

If only they were ...

Next stop, some code for normal, unexpandable BASIC dialects, such as perhaps the BBC's one (does it have the equivalent of LEX files, to seamlessly, transparently add new keywords to the language ?.)

V.

[EdS2]

.
Reflecting on how I might myself present a similar challenge, for example in a BBC Basic forum, it feels to me that the design problem is to bring to the attention of the reader a certain interesting property of integers, and then - surprise - to link it to pi. And then the challenge itself is to find ways of writing programs to explore this property.

In the case of the challenge on these forums, there's an additional aspect, which is to be able to make use of the marvellous augmentations to the 71B's Basic, the various ROMs. And as noted, they make a great deal of difference because of the relative slowness of interpreted language.

To keep with the topic of reflecting on Valentin's challenge, indeed in the HP Calculators subforum which is where we find this thread, I might wonder what can we do with the 15C, or the 42S, or perhaps the 41, on this problem. The 15C is of course not augmented - it knows nothing of primes or factorisations - whereas it's possible the 42S knows more, and surely there are excellent ROMs for the 41 which could help.

It is possible to augment BBC Basic with a ROM, but it wasn't very commonly done. I think it involves trapping a syntax error and re-parsing the failing statement. More common, perhaps, is to add some machine code (perhaps using the built-in assembler) and then CALL or USR that code, for some expensive inner function.

As to the unfairly accurate estimates of Pi which might come up from this counting challenge, it's clear to me that we won't find the best examples at round numbers. We might even find, like Buffon, that we can advantageously and unfairly batch up our results. (If indeed he did... ah, no, it was Lazzarini.)

[Valentin Albillo]

.
Hi, EdS2,

(03-24-2023 08:41 AM)EdS2 Wrote:  In the case of the challenge on these forums, there's an additional aspect, which is to be able to make use of the marvellous augmentations to the 71B's Basic, the various ROMs. And as noted, they make a great deal of difference because of the relative slowness of interpreted language.

First of all, if you continue to call them "augmentations" we'll have an additional issue. Frankly, it seems you do it on purpose to annoy me. FFS, you've used it and its friends 19 times already in this short thread. Did you learnt it recently and are you blabbering it all the time, baby-style ?.

Second, not that slow: getting rid of the FPRIME assembler keyword (not "predicate", gosh !) by substituting it by four lines of BASIC code incurs only a 24% runtime penalty, which isn't that much, actually.

Have a look at the promised code I've written for your first "related challenge", to run my solution 1 for the HP-71B but using an "unaugmented [sic] 71B, or a similarly ordinary Basic" (your words) instead.

Try it on your BBC BASIC micro and report how it goes, if you'd obligue.

V.

[EdS2]

(03-22-2023 09:11 AM)EdS2 Wrote:  My related challenges, then, would be
- to respond to Valentin's challenge using an unaugmented 71B, or a similarly ordinary Basic
- to find values of N which incidentally deliver unfairly close values of Pi

Having managed to produce a BBC Basic program - not by writing one, more by porting and augmenting solutions posted - here are some findings from my second challenge-to-self. I think it should be clear what I'm after, so here are some cherry-picked examples:

Code:

        28         17  3.14362099
       426        259  3.14145283
      1186        721  3.14159601
      5374       3267  3.14159277
     25684      15614  3.14159262
     89440      54373  3.14159265

(And in case it isn't clear, I keep a track of what the error in the estimate in Pi is, and look for a best-so-far lowest error, for each value of N. The above is a subset of the list of record-breakers. The final value listed above is rather close to pi: 3.1415926532...)

I used a fast BBC Basic engine for this - I simply ran the original program for successive values of N. Perhaps there's a cheaper way, as this does seem to be highly redundant.

[pier4r]

(04-01-2023 06:32 PM)EdS2 Wrote:  I used a fast BBC Basic engine for this - I simply ran the original program for successive values of N. Perhaps there's a cheaper way, as this does seem to be highly redundant.

Wouldn't be possible to keep track of the values as the computation goes?

Further, how much prevision does the BBC basic have? I am not sure but at a certain point the least significant digits may be correct due to computations errors if the precision is not that great.

[EdS2]

Yes, at some point we get so close to pi that our error calculation might be off. But the original BBC Basic I've been using, on a 6502 (emulator), has 5 byte floats with 32 bits of precision, and that's OK I think for the numbers I give above.

A modern BBC Basic on an Intel-based computer has 10 byte floats with 64 bits of precision, and can give us a few more record-holders:

Code:

    693570     421640  3.14159265
    793985     482685  3.14159265
    981973     596968  3.14159265

(I should have asked it to print more digits for the approximation, of course!)

And yes, I too suspect it should be possible to compute a running approximation somehow. But note that each term of the sum has a multiplier of N inside a floor function - it's not straightforward to factor N out.

[EdS2]

Nice to see Valentin has tackled this subsidiary challenge on his Pi Day thread.

Looks like a nice solution, probably very amenable to porting to BBC Basic, but I have a couple of ideas I'd like to follow before doing that.

(I need at least one large array, for the PRIM function, but I think I can re-use it for the Moebius values. On a plain original BBC Model B Micro with a floppy disk, I have about 25k bytes available. One can push that to nearly 30k. With a Master and a special (slow) version of BBC Basic (BAS128) one can get nearly 64k, whereas a Master Turbo runs faster than the original and offers 45k. Those are about the limits with the 6502 version of BBC Basic.)