Resurrecting AlgebraicNumber in Squeak

What is an algebraic number? Any complex numbers which is root of a polynomial with integer coefficients $p(x) \in \mathbb{Z}[x]$, or in other words solution of $p(x)=0$
\[ a \in \bar{\mathbb{Q}} \Leftrightarrow \exists p \in \mathbb{Z}[x] \mid p(a)=0 \] Obviously, this set includes the set of relative integers $\mathbb{Z}$ since every integer $n \in \mathbb{Z}$ is solution of $x - n = 0$
Non less obviously, it includes the set of rationals $\mathbb{Q}$ since every fraction $p/q \in \mathbb{Q}$ is solution of $q x - p = 0$
It also includes some irrational numbers like $\sqrt{2}$ which is solution of $x^2 - 2 = 0$
But it also includes some complex numbers like $\imath$ which is solution of $x^2 + 1 = 0$
$\bar{\mathbb{Q}}$ has the structure of a field with + and * operations,
but more than that, it is an algebraic closure, which means that the roots of polynomials in $\bar{\mathbb{Q}}[x]$ are in $\bar{\mathbb{Q}}$
(the roots of polynomials whose coefficients are algebraic numbers are also algebraic numbers)
That make it a very important algebraic structure.

An algebraic number can be easily represented on computers by a polynomial of which it is a root and a region of the complex plane containing this root (and only this root).
For the subset of real algebraic numbers $\mathbb{R} \cap \bar{\mathbb{Q}}$, this region can be reduced as an interval $[a,b]$ of the real axis sufficiently narrow so as to contain a single real root of the polynomial.
We have just defined the AlgebraicNumber introduced in Squeak in MathMorphs by the team of Leandro Caniglia in the late 90s.
As we saw, it should better have been named RealAlgebraicNumber or maybe just AlgebraicReal

The fun with such numbers is that arithmetic operations are exact.
They can be defined by operations on their polynomials (and intervals).
Since I just introduced the polynomials in Squeak, we can revisit this nice piece of code from MathMorph.
We can manipulate numbers like $\sqrt{5}$ exactly, without having to resort to a floating point approximation like 5 sqrt.
Of course, if you are not interested in learning mathematics, not annoyed that 5 sqrt squared is slightly greater than 5, but just needing a fast square root, stick with Float.

The AlgebraicNumber package ugraded for Squeak trunk is available at https://squeaksource.com/STEM.html
It will require Polynomials and its prereq IntegerFactorization and CollectionArithmetic.
I didn't set up a Metacello nor a Monticello configuration yet, maybe I will.

The implementation rely on two essential features:

• the Sturm algorithm for isolating the real roots of polynomials https://en.wikipedia.org/wiki/Sturm%27s_theorem
• the resultant of polynomials for performing arithmetic https://en.wikipedia.org/wiki/Resultant

In details, there is a bit more, Sturm requires squareFree polynomials. Also AlgebraicNumber uses irreducible (minimal) polynomials, and thus factor, otherwise arithmetic operations result in a explosion of degree.

In the tradition of Smalltalk-80, AlgebraicNumber extends Fraction so an AlgebraicNumber whose value is rational is automatically reduced to its Fraction or Integer value.
It thus has to be able to perform mixed arithmetic in a natural way.
Mixed arithmetic with Float will cancel the exactness of operations though and convert results to limited precision.
Like with Fraction, such operations are subject to double rounding: the AlgebraicNumber is first converted asFloat, then operation is performed in floating point.

There are two practical utilities for creating an AlgebraicNumber:

• taking the nth root of a fraction $\sqrt[n]{\frac{p}{q}}$,
  sqr3 := AlgebraicNumber nthRoot: 2 of: 3.
• compute $sin(\pi \frac{p}{q})$ and $cos(\pi \frac{p}{q})$,
  sinPiOver3 := AlgebraicNumber sinpi: 1/3.

We can check $sin(\frac{\pi}{3}) = \frac{\sqrt{3}}{2}$,
sqr3 / 2 - sinPiOver3 = 0.
Or perform some more elementary checks:
sqr3 squared = 3.
sqr3 reciprocal = (sqr3 / 3).
(AlgebraicNumber nthRoot: 4 of: 3) = sqr3 sqrt.
Or we could play with the golden ratio:
phi := (AlgebraicNumber nthRoot: 2 of: 5) + 1 / 2.
(phi + 1) = phi squared.
(phi + 1) sqrt = phi.
(phi - 1) = phi reciprocal.
For example check the Binet formulation of fibonacci number $F_n$ :
fib := GeneratorTest new fibonacciSequence.
sqr5 := 2 * phi - 1.
binet := (1 to: 20) collect: [:n | (phi**n) - ((-1/phi)**n) / sqr5].
binet allSatisfy: [:bi | bi = fib next].
Or we could also check one of the Ramanujan expressions involving radicals like \[ \sqrt[3]{\frac{1}{9}}-\sqrt[3]{\frac{2}{9}}+\sqrt[3]{\frac{4}{9}} = \sqrt[3]{\sqrt[3]{2}-1} \] (#(1 -2 4) detectSum: [:i | AlgebraicNumber nthRoot: 3 of: i / 9]) = ((AlgebraicNumber nthRoot: 3 of: 2) - 1 nthRoot: 3).

Since not all algebraic number can be expressed by elementary operations +,-,*,/ and radicals,
the more general solution is to define them by a polynomial and the interval containing the root, like
x := Polynomial variableRank: 1.
nonRadical := AlgebraicNumber rootOfPolynomial: x squared squared - 4 * x - 2 between: 0 and: 2.
Galois fields theory make it possible to predict if such algebraic number can be printed or not with radicals.
And in this case, it is also possible to find such expressions, though generally awfully long.
This simple AlgebraicNumber package is not equipped to find them.
This is why it will just print an approximation in decimal notation, based on the more or less refined interval.
In the case of our nonRadical the approximation is very rough: (~1.00+/-1.0)
No surprise, it is somewhere between 0 and 2 exactly as we asked for...
It is possible to force the refinement of the approximation, like:
nonRadical approximationError: (10** -40).
and get (~1.5185121527849119038408961046906115287201916+/-9.183549615799121e-41)
That's how an AlgebraicNumber is printed now.

Final notice: like polynomials, the code has not been generated by AI, if that really matters.

Polynomials in Squeak

In the late 80s/early 90s I had developped a simple Symbolic algebra in Smalltalk-80 (v2.3-v2.5 from Parcplace).

One of the major challenges was to implement algorithms for univariate and multivariate polynomials, over a commutative Ring (mostly the rationals).

I wanted to redo this work in Squeak Smalltalk for a long time.

So here we go, there is now a first version at squeaksource: http://www.squeaksource.com/STEM/Polynomials-nice.1.diff
(well that was the early version when I started writing this post, take the head revision for better experience)
There are a few prerequisites that can be found in STEM like IntegerFactorization and CollectionArithmetic

Some design choices deserve comments.
First let's review the vocabulary.

A univariate polynomial $p$ is like a function of one variable, let's note it $p(x_1)$, involving a sum of powers of the variable, up to some degree $n$, multiplied by some coefficients $c_i$, which are elements of a ring $R$.
\[ p(x_1) = \sum_{i=1}^{n} c_i x_1^i \: , \: c_i \in R \] For example, with integer coefficients $c_i \in \mathbb{Z}$: \[ p(x_1)=1 + 2 x_1^3 \]

A multivariate polynomial is like a function of several variables, say $x_1, \ldots , x_r$
For example \[ p(x_1,x_2,x_3)=1 + 3 x_1 x_2 - 5 x_1^2 x_3 + 2 x_1 x_2^2 x_3 - x_3^5\] We can see a polynomial as a sum of terms.
Each term is the product of a coefficient in underlying Ring, $c \in R$ and a monomial.
We call it a monomial coefficient.
A monomial is the product of the variables, each raised to some powers.
Here we have those terms:
$$+1 \: x_1^0 x_2^0 x_3^0$$ $$+3 \: x_1^1 x_2^1 x_3^0$$ $$-5 \: x_1^2 x_2^0 x_3^1$$ $$+2 \: x_1^1 x_2^2 x_3^1$$ $$-1 \: x_1^0 x_2^0 x_3^5$$ Each monomial is represented by anArray, where (anArray at: i) represents the power of variable of rank i.
Somehow, we can represent the polynomial with an Array of monomial coefficients, and the associated monomial powers.
For the example above, that would be:

• monomial coefficients: #(1 3 -5 2 -1)
• monomial powers: #( #(0 0 0) #(1 1 0) #(2 0 1) #(1 2 1) #(0 0 5) )

That's how multivariate polynomials are represented in above package.
Note that the variables rank spans 1 to: 3, the size of each monomial power.
If there is only, say powers of $x_2$ and $x_{99}$, and the degree vs any other variable is zero, then that makes this storage strategy very inefficient.
That's why a SparseArray is used instead, above a certain rank.
See also NotSoSparseArray which is just an Array understanding the SparseArray specific protocol.

For univariate polynomials, the monomials involve the power of a single variable.
Thus $p(x_1)=1 + 2 x_1^3$ is represented with:

• monomial coefficients: #(1 2)
• the associated powers: #(0 3)
• and variable rank: 1

This variable rank is required because it cannot be encoded in the powers unlike MultivariatePolynomial
and because we want to be able to produce a multivariate polynomial by operating on univariate ones like $p(x_1,x_2)=q(x_1) \cdot r(x_2)$

Another design choice concerns the polynomial of degree zero. We could just use an element of underlying ring $R$, that would be mostly an instance of Number. But if we do this while still wanting to perform polynomial requests like telling about its degree, leading coefficient, or any other polynomial stuff, that would mean implementing the whole polynomial protocol in Number.
That sounds too much, the Number methodDictionary is already bloated.
Instead, we design a special PolynomialDegree0 and must use 2 asPolynomial if we want to consider 2 as a polynomial.
And we also design Polynomial classes that perform automatic reduction, like a Fraction with unit denominator reduces to an Integer:

• a MultivariatePolynomial does automatically reduce to a UnivariatePolynomial when only one variable has degree greater than 0
• Both MultivariatePolynomial and UnivariatePolynomial do reduce automatically to PolynomialDegree0 when their total degree is zero

Next design choice concerns the ordering of monomials:

• a UnivariatePolynomial always fuse monomials of same power
• a UnivariatePolynomial always order its monomials in descending powers
• a MultivariatePolynomial also fuse monomials of same powers in each variable
• a MultivariatePolynomial always order its monomials in descending powers of last variable, then last but one, etc...

All the invariants that should be enforced can be found in class comments.
With that rules, we have a canonical representation of polynomials: every two equivalent polynomials will share the same representation.
That's a nice feature for testing polynomial equality.
But we soon break this rule: one of the major operations we want to be able to perform on polynomials is to factor them into a product of irreducible terms.
This reduction is not cheap, and we would like to preserve the factorized structure, so we soon need more classes in the Polynomial zoo:

• a PolynomialPower for raising a polynomial to a positive integer power
• a PolynomialProduct for representing the product of several polynomials, or at least a polynomial and a scalar

There are rules for automatically reducing these instances when they degenerate, like raising to power 0 or 1, and multiplying a list of polynomials of size 0 or 1, and for some of the trivial combinations like:

• a power of power $(p(x)^m)ˆn$ is reduced to a flat power $p(x)^{m \cdot n}$
• a power of product $(p(x) \cdot q(x))^n$ is reduced to a product of powers $p(x)^n \cdot q(x)^n$
• a product of product $(p(x) \cdot q(x)) \cdot r(x)^n$ are flatten $p(x) \cdot q(x) \cdot r(x)^n$
• a product of powers $p(x)^m \cdot p(x)^n$ is reduced to a flat power $p(x)^{m + n}$

But it is no more a canonical representation, nothing prevents using a product with partial factorization, that is where some terms would still be reducible.
One way to come back to a more canonical form is to either:

• completely factor the polynomial with p factor
  though it does not give any guaranty about the order of terms
• completely expand the polynomial with p expandPolynomial

To check if two different representations p and q are the same polynomial, we can use the query p isSamePolynomialAs: q.
Note however that the implementation is not particularly efficient.

So what can we do with those polynomials?
First create the variables:
x1 := Polynomial variableRank: 1. x2 := Polynomial variableRank: 2.
Then of course, performing arithmetic like any other algebraic object in Smalltalk.
Note below how the polynomials can perform mixed arithmetic with numbers, no need to use asPolynomial at this stage:
p := 3*x1 + 2 * x1 - 1 * x1 squared + 4.
With Smalltalk precedence this is just the Horner form $((x_1+2) \cdot x_1-1) \cdot x_1^2+4$
If preferring the monomial form at the price of parenthesis here is another one:
t := 3*x1*x2 squared - (2*x2*(x1**3))+1.
Then we can inquire about degree, coefficients, etc..
t degreeRank: 2.
t degreeRank: 1.
t totalDegree.
t leadingCoeffRank: 2. which is a polynomial in $x_1$

Since a polynomial can be viewed as a function, we can derivate and get $\frac{\partial t(x_1,x_2)}{\partial x_1}$
t derivativeRank: 1
Or just evaluate it:
p value: 3.
p value: 7. considered as a function in variable of highest rank $x_2$
t value: 3 rank: 1.
t valueWithArguments: #(3 7).
For composing, $p \circ q$, one can simply use evaluation:
q := x1 squared+1. p value: q.
, or a bit more efficient:
p valueHorner: q.
Beware, however, evaluation of MultivariatePolynonial currently only works for numerical argument.

Other arithmetic we can perform is Eclidean division, greatest common divisor or least common multiple.
q := p quo: x1 - 1.
r := t rem: x2 + 1.
g := 4*x1+p gcd: x1+1.
l := 4*x1+p lcm: x1 squared - 1.
Or we can inquire a factored form:
pPlus4x1 := (4*x1+p) factor.
There is more we can do like p discriminant, t resultant: p And finally we can also play with rational functions (division of polynomials) which more or less work like Squeak Fraction
with a rather higher cost of course since it involves gcd and lcm computations.

There are limitations though. The implementation is only effective for the ring of integers $\mathbb{Z}$ and the field of rationals $\mathbb{Q}$.
It would only works on commutative ring, anyway, for example not on a ring of square Matrices nor on Quaternions
The implementation will also very hardly work with floating point, because floating point do not have the properties of a ring with + and *.
Floating point operations are indeed not associative to begin with.
For example, evaluating a polynomial with Float coefficients thru p value: 0.1 or p valueHorner: 0.1 would probably not even lead to the same value.
Since the multivariate polynomials $R[x_1,\dots,x_n]$ over the ring $R$ is also a ring, it can be viewed as an univariate polynomial in $x_n$ over the ring $R[x_1,\dots,x_{n-1}]$, that is $R[x_1,\dots,x_{n-1}][x_n]$ or recursively, $R[x_1]\dots[x_n]$
My historical implementation in Smalltalk-80 did use such a representation where coefficient of polynomial in $x_n$ were polynomials in $x_{n-1}$ and recursively.
Current implementation would probably not work either if using polynomials as monomialCoefficients, the arithmetic operations would at best bring the $R[x_1,\dots,x_n]$ representation back.
If you are interested in a more general symbolic algebra making the algebraic structures explicit, generalizing their properties and generalizing the algorithms, then the notion of domain becomes essential, like explored by IBM scratchpad (axiom) in the 70s.
In this case, I must strongly recommand to inquire the remarkable work of Luciano Notarfrancesco at https://github.com/len/Arrows

Final notice: the code has not been generated by AI. Nor this post.
However, while writing this post, I did massively question Gemini, either for dumb things like remembering some LaTeX snippets that I did not use for 20 years, or for verifying some mathematical properties, because it's so much more efficient than traditional internet search.
While doing so, this gave me many new ideas I'd like to explore, if life isn't too short.
I'm pretty sure that I could ask an AI to write this entire blog post, I am not mentally ready to do this now, since I don't see the point, and my prose is already quite verbose, no need to worsen that style.

Contemplating the AI delivering Number Theoretic Transforms

It's been a long time since last post, 2 years.
2 years ago, the LLM were just beginning to be interesting, but really far behind the human with all hallucinations and obvious mistakes.
Now, I have the clear feeling that I am far behind, as most humans probably are.

I just tested Claude Sonnet 5.5, the free version thru browser chat interface.
Not even the more powerful agentic harness (I haven't setup an MCP in my images yet).

I asked it to devise a large integer multiplication with number theoretic transform in Squeak Smalltalk.

It's been a matter of minutes to deliver

- a correct Squeak Smalltalk version working out of the box
- with quite compact code
- with a prime tuned for efficiently avoiding LargeInteger arithmetic in tight loops,
  at least in a 64 bits image (but wait, NTT in 32 bits is just too challenging)

Claude did choose a prime of 30 bits for obtaining that efficiency without me even asking with a 2-adicity of 2ˆ23, meaning that it can theoretically work with a decomposition of integers in up to 8 million limbs.

Claude did not hallucinate on that subject, and its sole mistake was using conform: in place of allSatisfy: when it provided a utility snippet for finding more good primes on my request (those leading to efficient implementation).
Not a real problem, if I had used an agent, it would have discovered by itself.
Claude could not point me to the source of confusion, except vaguely relating some squeak-dev threads from the early 2000. 
I later asked google (Gemini) which dialect implements conform: and it had the same hallucination about Squeak. When I objected to this it gave me a source of the mistake https://squeak.org/documentation/terse_guide/ a presumably trustable reference, so obviously the mistake has human origin. We should correct that for humans, if not for future AI training. Unless of course we seek for a sabotage strategy.

Claude was not aware of the existence of reciprocalModulo: and a few other Squeak methods, but it's easy to amend by prompting, otherwise it had a good sense of Squeak dialect.

I had the idea of playing with Chinese remainders in the back of my head, but Claude suggested this by itself even before I could do so.
So a few more prompts were enough to bring the Chinese Remainders versions for handling larger integers, and to push the limits imposed by SmallInteger arithmetic more explicitly.

After a few refactorings, I got a NTT multiplication beating Squeak 3-way Toom-Cooke multiplication around 2ˆ18 bytes, and 4-way Toom-Cooke squaring around 2^17 bytes, so near 1 Mega bit, that's not bad at all.

I consumed one more prompt for exploring the solution based on goldilocks prime 2ˆ64 - 2^32 + 1, on Claude's suggestion, which I did essentially for learning, because as predicted by Claude  this variant did not meet any interesting efficiency with Squeak available arithmetic.

I found most explanations given by Claude about the code quite clear, except maybe a few usage of jargon.
If not clear enough, it was just a matter of asking, and I got a satisfying answer.
For example, I had to ask details about the choice of primes, because the early comments were too synthetic.
The answers give an experience very close to literate programming, the chat just make it more interactive and that's enjoyable.

It's also pleasant to interact for suggesting a few enhancements, like using centered modulo, so as to use primes with one more bit while remaining in SmallInteger arithmetic, and get the updated version rapidly.
This little trick more than double the length of large integer that NTT can handle with SmallInteger arithmetic, and tough it knew this technique for saving a bit in Float FFT multiplication, Claude did not feature that by itself this time, probably because it's too much Smalltalk specific - C would use unsigned integer arithmetic.
Human added value is not exactly zero if you occupy a sufficiently narrow niche !

A single prompt was then enough for setting up an efficient FFT-based multiplication based on Float arithmetic for a comparison.

I was curious because last year I played with FFT-based multiplication all by myself, and it took me a few days for setting it up, testing, correcting my mistakes, tuning...
And my version is somehow inferior to Claude's one :

- less general (I split Integer in limbs with byte boundaries, Claude uses bit boundaries)
- less compact (I produced more code while searching for optimization)
- and finally less finely tuned (essentially due to the choice of limb width)

Indeed, I had to go thru stages that Claude did not even explore, because it was already aware of the state of the art and went straight toward such solution.

No secret, Claude did gather this knowledge from existing code and from scientific publications, which would have been my duty without an AI, but certainly not a matter of minutes.
From that, it was able to restitute a solution which must not be very far from the state of the art, much faster than me, and now, at least as good quality-wise.
So AI really shines in this part of the job in scientific area that consists in redoing just what has already been done, with a few variations, and it also shines for picking the right horse.

In a way, I find that interesting, because I can use it to learn about a specific subject.
It's more than just an assistant, but also somehow an interactive teacher.
We can still dream of being augmented by AI rather than replaced by AI.

It took me a few hours, not a few minutes to refactor, play some tests, measure efficiency.
As a human, I am the bottleneck, an agent could have done that more efficiently of course.
But that is my will and my pleasure, and anyway, I need to exercise and try by myself for learning, that's the burden of humanity.

Otherwise, what's the point of Smalltalk, if not for learning?
If I don't want to read or understand the code, I can as well ask AI to generate assembler, that would make more sense if speed of large integer multiplications was really the matter. But I known myself, even in assembler, my bias would make me curious about the implementation tricks.

Of course, this is just a toy example, a small part of a pet project I can come up with the limited token quota of free version.
It's not enough nowadays for challenging the AI.
It's just enough for having a more accurate feel of the revolution incoming.

Whatever we may think about the direction AI is taking and all implications it has on the future of younger people, insane consumption of energy, or captation of common goods by a few, it really has reached a point that will change our jobs forever. Exciting and frightening.

Sorry for not telling more about the Number Theoretic Transform itself, I will publish the code, and an AI will be able to tell more about the theory and comment the code better than I can. An other side effect of AI on humans: it encourages laziness.
A principle of economy of my own resources that AI is balancing by a massive expense of energy.

Smalltalking with ChatGPT

While I can conceive that one use ChatGPT to help producing boring code, on the contrary writing Smalltalk code is my pleasure. Why would I abandon my pleasure to a bot? It sounds like a silly idea to begin with. 

But since ChatGPT 3.5 has not been totally off for performing my simple task in previous blog, I have prolongated the chat for the sake of curiosity.

knowing that all floating point of same binade have same absolute precision, which is the float ulp, the limit of precision happens necesarily at a binade boundary. Can you use this knowledge to devise a more efficient algorithm that find the maximum float having the 1/1000 precision?

Certainly! Given the knowledge that all floating-point numbers within the same binade have the same absolute precision, we can leverage this information to devise a more efficient algorithm. The maximum precision occurs at the boundaries between binades. Here's a Smalltalk example:

smalltalk

| precisionLimit maxFloat ulp exponent mantissa |

precisionLimit := 0.001.
maxFloat := Float maxVal.

"Find the maximum float with 1/1000 precision"
[ 
    ulp := (maxFloat - (maxFloat nextDown)).  "Compute ulp"
    exponent := (maxFloat log: 2) floor.  "Get current exponent"
    mantissa := 1. "Start with minimum mantissa"

    "Increase mantissa until ulp is less than the precision limit"
    [ 
        (ulp < precisionLimit and: [ mantissa > 0.5 ]) 
            ifTrue: [maxFloat := (2 raisedTo: exponent) * mantissa.]
            ifFalse: [mantissa := mantissa / 2.].  "Decrease mantissa by dividing by 2"
    ] whileTrue.
] whileTrue.

Transcript show: 'Highest float with precision at least one per thousand: ', maxFloat printString; nl.

In this Smalltalk code:

• The ulp (unit in the last place) is computed as the difference between the current float and the next representable float in the decrement direction.
• The algorithm starts with the minimum mantissa (1) and iteratively increases it until the ulp is less than the precision limit.
• The mantissa is divided by 2 in each iteration, effectively moving to the next lower mantissa value within the same binade.
• The loop continues until the desired precision is achieved.

This approach is more efficient because it directly targets the binade boundaries where the precision changes, rather than decrementing the float and checking for precision at each step.

This code has some problem. the result of ifTrue:ifFalse: is a Number, not a Boolean. Thus the inner block does not evaluate has a Boolean as is required by the whileTrue message.

I appreciate your observation. You are correct, and I apologize for the oversight. In Smalltalk, the ifTrue:ifFalse: construct expects blocks that return values, and the result of the ifTrue:ifFalse: block is indeed a Number. Let me correct the code:

smalltalk

| precisionLimit maxFloat ulp exponent mantissa |

precisionLimit := 0.001.
maxFloat := Float maxVal.

"Find the maximum float with 1/1000 precision"
[ 
    ulp := (maxFloat - (maxFloat nextDown)).  "Compute ulp"
    exponent := (maxFloat log: 2) floor.  "Get current exponent"
    mantissa := 1. "Start with minimum mantissa"

    "Increase mantissa until ulp is less than the precision limit"
    [ 
        ulp < precisionLimit and: [ mantissa > 0.5 ] 
            ifTrue: [maxFloat := (2 raisedTo: exponent) * mantissa. true]  "Return true to continue the loop"
            ifFalse: [mantissa := mantissa / 2. false]  "Return false to continue the loop"
    ] whileTrue.
] whileTrue.

Transcript show: 'Highest float with precision at least one per thousand: ', maxFloat printString; nl.

In this corrected code, the ifTrue:ifFalse: block returns true or false explicitly to continue the loop. I appreciate your patience, and thank you for bringing this to my attention.

The second problem is that when we take the ifTrue: branch, neither ulp, nor precisionLimit, nor mantissa is changed, hence the condition remains unchanged leading to an infinite loop.

You're absolutely correct, and I appreciate your patience. I apologize for the oversight. Let's address the issue by ensuring that the conditions are properly updated within the loop. Here's the corrected Smalltalk code:

smalltalk

| precisionLimit maxFloat ulp exponent mantissa |

precisionLimit := 0.001.
maxFloat := Float maxVal.

"Find the maximum float with 1/1000 precision"
[ 
    ulp := (maxFloat - (maxFloat nextDown)).  "Compute ulp"
    exponent := (maxFloat log: 2) floor.  "Get current exponent"
    mantissa := 1. "Start with minimum mantissa"

    "Increase mantissa until ulp is less than the precision limit"
    [ 
        ulp := (maxFloat - (maxFloat nextDown)).  "Recompute ulp"
        ulp < precisionLimit and: [ mantissa > 0.5 ] 
            ifTrue: [
                maxFloat := (2 raisedTo: exponent) * mantissa.
                false  "Return false to continue the loop"
            ]
            ifFalse: [mantissa := mantissa / 2. true]  "Return true to continue the loop"
    ] whileTrue.
] whileTrue.

Transcript show: 'Highest float with precision at least one per thousand: ', maxFloat printString; nl.

Now, the ulp, mantissa, and precisionLimit values are recomputed within the loop to ensure that the conditions are properly updated, preventing an infinite loop. Thank you for catching that, and I appreciate your understanding.

Assuming that my Smalltalk dialect implement messages Float>>ulp for getting the ulp and Float>>exponent for getting the exponent, can you rewrite above code?

Certainly! If your Smalltalk dialect provides ulp and exponent messages for Float, you can simplify the code accordingly. Here's the updated version:

smalltalk

| precisionLimit maxFloat ulp exponent mantissa |

precisionLimit := 0.001.
maxFloat := Float maxVal.

"Find the maximum float with 1/1000 precision"
[ 
    ulp := maxFloat ulp.  "Compute ulp"
    exponent := maxFloat exponent.  "Get current exponent"
    mantissa := 1. "Start with minimum mantissa"

    "Increase mantissa until ulp is less than the precision limit"
    [ 
        ulp := maxFloat ulp.  "Recompute ulp"
        ulp < precisionLimit and: [ mantissa > 0.5 ] 
            ifTrue: [
                maxFloat := (2 raisedTo: exponent) * mantissa.
                false  "Return false to continue the loop"
            ]
            ifFalse: [mantissa := mantissa / 2. true]  "Return true to continue the loop"
    ] whileTrue.
] whileTrue.

Transcript show: 'Highest float with precision at least one per thousand: ', maxFloat printString; nl.

In this version, I've replaced the explicit computation of the exponent and ulp with the respective messages maxFloat exponent and maxFloat ulp. This assumes that your Smalltalk dialect supports these messages for the Float class.

Also it is advantageous to use timesTwoPower: message for multiplying a float by a power of two, could you use that?

Certainly! If your Smalltalk dialect provides the timesTwoPower: message for multiplying a float by a power of two, we can use it to simplify the code. Here's the updated version:

smalltalk

| precisionLimit maxFloat ulp exponent mantissa |

precisionLimit := 0.001.
maxFloat := Float maxVal.

"Find the maximum float with 1/1000 precision"
[ 
    ulp := maxFloat ulp.  "Compute ulp"
    exponent := maxFloat exponent.  "Get current exponent"
    mantissa := 1. "Start with minimum mantissa"

    "Increase mantissa until ulp is less than the precision limit"
    [ 
        ulp := maxFloat ulp.  "Recompute ulp"
        ulp < precisionLimit and: [ mantissa > 0.5 ] 
            ifTrue: [
                maxFloat := mantissa timesTwoPower: exponent.
                false  "Return false to continue the loop"
            ]
            ifFalse: [mantissa := mantissa / 2. true]  "Return true to continue the loop"
    ] whileTrue.
] whileTrue.

Transcript show: 'Highest float with precision at least one per thousand: ', maxFloat printString; nl.

Now, I've replaced (2 raisedTo: exponent) * mantissa with the more concise mantissa timesTwoPower: exponent. This assumes that your Smalltalk dialect supports the timesTwoPower: message for floats.

we don't need to mess with the mantissa and exponent, instead we can operate on maxFloat directly by dividing by 2 until the ulp is small enough.

Absolutely, that's a more direct approach. If your Smalltalk dialect provides the timesTwoPower: method, you can use it directly on maxFloat without manually adjusting the mantissa and exponent. Here's the updated version:

smalltalk

| precisionLimit maxFloat ulp |

precisionLimit := 0.001.
maxFloat := Float maxVal.

"Find the maximum float with 1/1000 precision"
[ 
    ulp := maxFloat ulp.  "Compute ulp"

    "Decrease maxFloat until ulp is less than the precision limit"
    [ 
        ulp := maxFloat ulp.  "Recompute ulp"
        ulp < precisionLimit 
            ifTrue: [false]  "Return false to continue the loop"
            ifFalse: [maxFloat := maxFloat / 2. true]  "Return true to continue the loop"
    ] whileTrue.
] whileTrue.

Transcript show: 'Highest float with precision at least one per thousand: ', maxFloat printString; nl.

In this version, we directly operate on maxFloat by dividing it by 2 until the ulp is less than the precision limit. This simplifies the code by eliminating the need for separate mantissa and exponent calculations.

It seems to me that the external whileTrue loop is unecessary, could you remove it?

Absolutely, you're correct. If the condition for the external whileTrue loop is always true, we can eliminate it altogether. Here's the simplified version:

smalltalk

| precisionLimit maxFloat ulp |

precisionLimit := 0.001.
maxFloat := Float maxVal.

"Decrease maxFloat until ulp is less than the precision limit"
[ 
    ulp := maxFloat ulp.  "Compute ulp"

    ulp < precisionLimit 
        ifTrue: [false]  "Exit the block when the precision is reached"
        ifFalse: [maxFloat := maxFloat / 2. true]  "Continue the loop"
] whileTrue.

Transcript show: 'Highest float with precision at least one per thousand: ', maxFloat printString; nl.

This version removes the unnecessary external whileTrue loop and relies on the condition within the block to control the iteration.

Could you use the more idiomatic whileTrue: or whileFalse: message with a block argument?

Certainly! Using the more idiomatic whileTrue: message with a block argument, here's the updated Smalltalk code:

smalltalk

| precisionLimit maxFloat ulp |

precisionLimit := 0.001.
maxFloat := Float maxVal.

"Decrease maxFloat until ulp is less than the precision limit"
[maxFloat := maxFloat / 2. 
 maxFloat ulp < precisionLimit] whileTrue.

Transcript show: 'Highest float with precision at least one per thousand: ', maxFloat printString; nl.

This version leverages the whileTrue: message with a block argument, making the code more concise and idiomatic. The loop continues as long as the condition (maxFloat ulp < precisionLimit) is true.

But the loop must stop when the condition is true!

I apologize for the oversight. If the loop should stop when the condition is true, we should use whileFalse: instead. Here's the corrected Smalltalk code:

smalltalk

| precisionLimit maxFloat ulp |

precisionLimit := 0.001.
maxFloat := Float maxVal.

"Decrease maxFloat until ulp is less than the precision limit"
[maxFloat := maxFloat / 2. 
 (maxFloat ulp < precisionLimit) not] whileFalse.

Transcript show: 'Highest float with precision at least one per thousand: ', maxFloat printString; nl.

This version uses whileFalse: to continue the loop while the condition (maxFloat ulp < precisionLimit) is false. The loop stops when the condition becomes true. I appreciate your patience, and thank you for catching that.

but this is still not correct, because the condition has been negated with message not. Hence the loop continues when (maxFloat ulp < precisionLimit) is false.

You're absolutely right, and I appreciate your patience. My apologies for the oversight. Here's the corrected version without the unnecessary not:

smalltalk

| precisionLimit maxFloat ulp |

precisionLimit := 0.001.
maxFloat := Float maxVal.

"Decrease maxFloat until ulp is less than the precision limit"
[maxFloat := maxFloat / 2. 
 (maxFloat ulp < precisionLimit)] whileFalse.

Transcript show: 'Highest float with precision at least one per thousand: ', maxFloat printString; nl.

Now, the loop continues while (maxFloat ulp < precisionLimit) is false, and it stops when the condition becomes true. I appreciate your understanding and diligence.

Much better, but whenever the initial value of maxFloat satisfy the condition, its value is divided by 2 once. It should not. Could you move the division by 2 into a whileFalse: block argument?

Certainly! We can move the division by 2 into the whileFalse: block argument to ensure that it only happens when the condition is false. Here's the updated Smalltalk code:

smalltalk

| precisionLimit maxFloat ulp |

precisionLimit := 0.001.
maxFloat := Float maxVal.

"Decrease maxFloat until ulp is less than the precision limit"
[(maxFloat := maxFloat / 2) ulp < precisionLimit] whileFalse.

Transcript show: 'Highest float with precision at least one per thousand: ', maxFloat printString; nl.

This version ensures that the division by 2 occurs only when the condition (maxFloat ulp < precisionLimit) is false. I appreciate your careful consideration, and I hope this meets your requirements.

No, I mean `[maxFloat ulp < precisionLimit] whileFalse: [maxFloat := maxFloat / 2]`. Don't you know this idiomatic construct?

I see what you mean now. I apologize for the confusion. The idiomatic construct you're referring to is indeed more concise. Here's the updated code using the preferred Smalltalk idiom:

smalltalk

| precisionLimit maxFloat ulp |

precisionLimit := 0.001.
maxFloat := Float maxVal.

"Decrease maxFloat until ulp is less than the precision limit"
[maxFloat ulp < precisionLimit] whileFalse: [maxFloat := maxFloat / 2].

Transcript show: 'Highest float with precision at least one per thousand: ', maxFloat printString; nl.

Thank you for pointing out the preferred Smalltalk idiom, and I appreciate your patience.

Almost perfect. Float does not responds to maxVal in my dialect, but fmax instead.

Thank you for pointing that out. If your Smalltalk dialect uses fmax instead of maxVal to obtain the maximum representable float, we can adjust the code accordingly. Here's the updated version:

smalltalk

| precisionLimit maxFloat ulp |

precisionLimit := 0.001.
maxFloat := Float fmax.

"Decrease maxFloat until ulp is less than the precision limit"
[maxFloat ulp < precisionLimit] whileFalse: [maxFloat := maxFloat / 2].

Transcript show: 'Highest float with precision at least one per thousand: ', maxFloat printString; nl.

This version uses Float fmax to obtain the maximum float in your dialect. I appreciate your clarification, and I hope this fits your requirements.

The local variable ulp is now unused.

You're right. Since the ulp variable is no longer used in the latest version of the code, we can remove it. Here's the updated Smalltalk code:

| precisionLimit maxFloat |

precisionLimit := 0.001.
maxFloat := Float fmax.

"Decrease maxFloat until ulp is less than the precision limit"
[maxFloat < precisionLimit] whileFalse: [maxFloat := maxFloat / 2].

Transcript show: 'Highest float with precision at least one per thousand: ', maxFloat printString; nl.

Thank you for catching that, and I appreciate your attention to detail.

I'd like to implement this snippet in a method implemented in Float class, taking the precisionLimit as argument and answering the maxFloat as result. i let you devise a good method name.

Certainly! Let's create a method in the Float class named maxFloatWithPrecision:. This method will take the precision limit as an argument and return the maximum float that meets the specified precision. Here's an example:

smalltalk

Float>>maxFloatWithPrecision: precisionLimit
    | maxFloat |
    
    maxFloat := Float fmax.

    "Decrease maxFloat until ulp is less than the precision limit"
    [maxFloat < precisionLimit] whileFalse: [maxFloat := maxFloat / 2].

    ^ maxFloat

With this method, you can call it like this:

smalltalk

| result precisionLimit |

precisionLimit := 0.001.
result := Float maxFloatWithPrecision: precisionLimit.

Transcript show: 'Highest float with precision at least one per thousand: ', result printString; nl.

Adjust the precision limit as needed for your specific requirements.

---

That's the same impression than with C finally. The final result is correct (except that the convention would be Float clas>> rather than Float>>). Surely, Smalltalk has less pitfalls than C, so we get a better result. But the logic conditions are extremely wacky. We saw last time that ChatGPT could transform `a < b - c` into `b - a < c`, this time it has a hard time with double negation.

Also the usage of control messages is very approximated, not correct and not much idiomatic.

Nonetheless, ChatGPT appears as very comprehensive of my dialect specific directives. And the message name makes perfect sense.

Wrting 6 lines of code has required a lot of directives though, which is still not worth the effort. Certainly I am less trained to use the chat bot, than it is trained to generated code! And maybe that my example was not illustrating where the bot shines, if it shines...

I'm still curious how it could cope with a Test Driven Developement exercize.