Smallissimo

Which fraction has these digits ?

2012-01-27T10:32:00.001+01:00

Sven Van Caekenberghe posted these nice example on Pharo mailing list:

1/9801 asScaledDecimal: 200.
1/998001 asScaledDecimal: 3000.

This also works with:

1/9801 printShowingDecimalPlaces: 200.

We could conjecture that 1/99980001 asScaledDecimal: 40000 should give a nice serie of digits too.

But if we want to form our own fraction, how can we construct it ?

Let us consider a fraction in interval ]0,1[

In "decimal" notation (I say decimal, but it is valid in any other base b), this fraction print in base b with a head composed of head size digits, and a repeated tail composed of tail size digits 0.head_tail_tail_tail...

For example, b := 10, head := '123' and tail := '40', will lead to (fraction asScaledDecimal: tail size * 3 + head size) = '0.123404040s9'

I once thought we could extend the syntax with such notation to denote infinite repetition:

fraction := 0.123(45).

But we didn't, so what is the value of the fraction?

Let us work with numbers rather than strings:

h := head ifEmpty: [0] ifNotEmpty: [Number readFrom: head base: b].

t := tail ifEmpty: [0] ifNotEmpty: [Number readFrom: tail base: b].

nh := head size.

nt := tail size.

If we had infinite series (lazy one would be preferred), we could define:

positiveIntegers := 1 to: Integer infinity.

infiniteTail := (positiveIntegers collect: [:i | t / (b raisedTo: nt * i)]).

fraction := (infiniteTail + h) / (b raisedTo: nh).

We have no such ready made infinity, nor lazy collections...

But even with such objects, we cannot evaluate this infiniteTail that easily.

Really? If we shift it by tail size digits, we will obtain this number with same infiniteTail:

infiniteTail * (b raisedTo: nt) = (t + infiniteTail).

From this property, it is easy to get:

infiniteTail := t / ((b raisedTo: nt) - 1).

If we want a leading head, we simply have to shift above fraction:

fraction := (t / ((b raisedTo: nt) - 1) + h) / (b raisedTo: nh).

So, let us reconstruct Sven example:

| base tail |
base := 10.
tail := ((0 to: 97) , #(99) inject: '' writeStream into: [:s :n | n printOn: s base: base length: 2 padded: true. s]) contents.
^(Number readFrom: tail base: base) / ((base raisedTo: tail size) - 1)

-> 1/9801

Good. Now, this expression won't lead to nice fraction for any tail...

If we want a fraction with numerator = 1, we need to find a divisor of the denominator 999....99. That's the nice thing with Sven example. I don't know how or where he found it, but I guess it'll be hard to be as elegant. Oh, maybe thanks to the extension to any base you can play with it too (you will need to implement printShowingDecimalPlaces:base:).

About Complex in Squeak

2011-10-12T01:34:00.000+02:00

I noticed that Complex is present in Squeak trunk and not in Pharo, though not used. So I suggested to create a shared external package common to the two distributions in this thread with a few other technical questions. This can benefit to Pharo users, and since the base of interested people did not seem that large to me, it can't be bad to group forces.

I got quite a few answers, thanks to all participants. Scattering my answers in many posts would lead to a lot of repetitions. Grouping the threads into a big string would be misbehaved. Though all subjects are tied together: this is a pretext for blogging. So here are my comments.

Why Complex was removed from Pharo?
I don't have a mailing list ref available, but I think it was on a path to a reduced image (not to say minimal). Complex is not used in a base image so it could be loaded on demand. It's the never ending discussions of whether constructing an image by assembling modules (packages?) or by stripping... Full vs Core...
Another argument was that Complex implementation was not that good (understand not universal - see below). So we should better let it leave the image and give a chance to alternate solutions (driven by different trade off/paradigms - see below).

Why Quaternion?

Quaternion are to SO3 what Complex are to SO2. They are often used in engineering, specially in robot arm control or spatial vehicle control because not suffering from Euler/Cardan angle singularities. They also replace trigonometric operations by algebraic ones (an advantage if there is no hard-wired trigonometry in FPU they just require sqrt). You can find Quaternion in squeaksource.
Octonions are generalization in spaces of dimension 4, so can be used in generalized relativity... I'm not aware of any Smalltalk implementation yet, maybe this niche is too narrow ;)

Why I've suggested to rename Complex as ComplexNumber?

Because more expressive. I think it was also to match isComplexNumber. But, sure, we can dissociate the two...
And also for stupid digression about an eventual distinction required by a symbolic Complex - apologize for the follow ups in the thread ;) .

Why I suggested to implement isComplexNumber rather than isComplex?

Because the collisions and are not just imaginary ;) but happening for real as Juan reported.
There is an issue opened in mantis.

Why -1 sqrt raise an Error instead of answering one of the complex roots?
This raises the question of the paradigm we want to use...
It's very much like Fraction and Integer.
In Smalltalk, every Fraction with a denominator = 1 is automatically reduced to an Integer. Thus we can consider that every Integer is a special kind of Fraction (with an implicit 1 denominator). It's thus logical to have Integer behave like a Fraction and respond to numerator, denominator, fractionPart, integerPart and even answer true to isFraction and to selfasFraction (if it quacks like a duck... I made the changes recently in trunk).
So we could consider the same with ComplexNumber and Number. Every Number is a Complex with a null imaginary part. Thus we could let Number behave like a Complex and extend its mathematical functions over complex domain (-1 sqrt, 2 arcCos, -1 argCosh etc...), and also let it answer to some Complex specific protocol like conjugated , imaginaryPart, realPart... In such paradigm, every Complex with a null imaginary part would be automatically reduced to a Number.

This discussion already took place before...

Above generalization is quite nice. But unlike the case of Integer and Fraction we are changing the behaviour of sqrt and potentially breaking compatibility.
Some applications do expect and require an Exception. If the Exception does not happen, then it will lead to delayed failure, or worse, incorrect results.
One possible incorrect result would be to fill your bank account we imaginary €.
Especially, knowing that 1 i isNumber now answers true, defensive "type" checks that worked once are now by-passed, which worsen the situation. See digression below.

The opposite point of view is to preserve compatibility and implement both behaviours.
The solution adopted in Squeak is to let the programmer explicitly force a Number to behave as a Complex with -1 asComplex sqrt. In such case, the Complex must not automatically be reduced to a Number when its imaginary part is null. this way we can also chain complex behaviors, and that's what we currently get with 2 asComplex sqrt arcSin we don't need to force a second time 2 asComplex sqrt asComplex arcSin.
Another solution would be to distinguish the selectors (like -1 sqrt versus -1 complexSqrt) and let Number and Complex respond to both protocols... But that's many mathematical functions to duplicate.

I don't think Squeak solution is that bad, especially if we consider that Complex are rarely used. It is somehow very much like statically typed languages FORTRAN/C++ so some of you might not like it ;)
The biggest grief is that we don't have static typing (not every receiver is literal), so we are forced to always send asComplex when in doubt. This can lead to scattering more such messages than manageable... That might happen in case of large usage of complex, but until now, nobody never complained.
So I would not change implementation without very careful thoughts, and probably not change it at all in trunk. Because such change is a complete shift of paradigm with potentially nasty side effects for Complex-unaware-apps.
If we want such change, I vote for removing Complex from image (the arguments already raised in Pharo).

Why 1 i isNumber answers true?
This was just to let this assertion work (0 i = 0) = (0 = 0 i), originated here... OK, they don't behave the same, but still can be equal (or shall they not?).
The second argument was that a Complex is a Number as the name doesn't tell - well ComplexNumber is more expressive, I warned you ;)

Why I want to revert this change?
1) This change did pierce now obsolete isNumber defence... And expectations: belongs to R suddenly became belongs to C.
If you analyze senders of isNumber, just tell me how many are equipped to handle a Complex? I tried it once.
It's the same compatibility problem as -1 sqrt...
2) Since both Number and Complex answer true to isNumber, then we logically shifted to the more general point of view that every Number is a ComplexNumber.
Unfortunately, this is not true, it does not match current paradigm. A Number does not behave like a Complex, so isNumber distinction is now quite useless.

The side effects are not imaginary ;), again this broke Quaternion implementation for example. I had to commit a quick workaround (x isNumber and: [x isComplex not]) but it's overkill, and just a smell that isNumber does not discriminate enough.
We now badly need isRealNumber - where real must be understood as not imaginary - see below.

Why reverting Complex>>isNumber rather than creating Number>>isRealNumber?
1) Please don't confuse isRealNumber and isFloat, that happens each time I open this can of worms. So I prefer to repeat myself, we have no representation of reals in Squeak, only of some special kinds of reals, Integer, Fraction and the degenerated inexact rounded fractions (I named Float).
We could represent many others (like AlgebraicNumber in MathMorph).
But real numbers are uncountable... See http://en.wikipedia.org/wiki/Cardinality
And we can only represent a countable subset of the real numbers... See http://en.wikipedia.org/wiki/Computability.
That's my first reason: this selector is confusing.
2) With current implementation of Complex, a Complex currently means two things:
- if imaginary part is not null, then it is a real ComplexNumber, well I mean unreal...
- if imaginary part is null, then it is a real number that behaves like a Complex
Shall 0 i isRealNumber answer true? The question will arise inevitably...
That's my second reason: this selector is confusing.
What we really seek is behavesLikeARealNumber. Ouch! I wish I never see such specialization in trunk ;)

Ah, the good old time when

isComplex did mean isKindOf: Complex
isNumber did mean isKindOf: Number (what Complex are not).

The alternative is to shift the paradigm and really let Number behave like Complex as told above, and assume all implications or reject the problem out the image... I have no easy solution handy.

Reviewing fraction asFloat

2011-09-25T12:31:00.000+02:00

After changing Integer>>asFloat, let's inspect Fraction>>asFloat. The spirit is quite simple: perform an Euclidean division after shifting the numerator or the denominator enough to obtain a certain precision in the quotient. The quotient will be the mantissa of the Float and the remainder are the truncated bits. Then come a dissertation on the conditions when to round upper or not. Let's see the code:

asFloat
    "Answer a Float that closely approximates the value of the receiver.
    This implementation will answer the closest floating point number to
    the receiver.
    It uses the IEEE 754 round to nearest even mode
    (can happen in case denominator is a power of two)"

    | a b q r exponent floatExponent n ha hb hq q1 |
    a := numerator abs.
    b := denominator abs.
    ha := a highBit.
    hb := b highBit.

    "If both numerator and denominator are represented exactly in floating point number,
    then fastest thing to do is to use hardwired float division"
    (ha < 54 and: [hb < 54]) ifTrue: [^numerator asFloat / denominator asFloat].

    "Try and obtain a mantissa with 54 bits.
    First guess is rough, we might get one more bit or one less"
    exponent := ha - hb - 54.
    exponent > 0
        ifTrue: [b := b bitShift: exponent]
        ifFalse: [a := a bitShift: exponent negated].
    q := a quo: b.
    r := a - (q * b).
    hq := q highBit.

    "check for gradual underflow, in which case we should use less bits"
    floatExponent := exponent + hq - 1.
    n := floatExponent > -1023
        ifTrue: [54]
        ifFalse: [54 + floatExponent + 1022].

    hq > n
        ifTrue: [exponent := exponent + hq - n.
            r := (q bitAnd: (1 bitShift: hq - n) - 1) * b + r.
            q := q bitShift: n - hq].
    hq < n
        ifTrue: [exponent := exponent + hq - n.
            q1 := (r bitShift: n - hq) quo: b.
            q := (q bitShift: n - hq) bitAnd: q1.
            r := (r bitShift: n - hq) - (q1 * b)].

    "check if we should round upward.
    The case of exact half (q bitAnd: 1) isZero not & (r isZero)
    will be handled by Integer>>asFloat"
    ((q bitAnd: 1) isZero or: [r isZero])
        ifFalse: [q := q + 1].

    ^ (self positive
        ifTrue: [q asFloat]
        ifFalse: [q = 0
            ifTrue: [Float negativeZero]
            ifFalse: [q asFloat negated]])
        timesTwoPower: exponent

Let's review the three parts in yellow : first the comment (can happen in case denominator is a power of two) is obscure. It is supposed to explain that in case of tie, the algorithm will round to nearest even. The case of tie can only happen if the denominator is a power of two indeed, any other factor in the denominator will lead to an infinite series of digits and thus cannot produce a tie. But this information is not essential enough to be placed in the header. After all it is not essential at all.

The second part in yellow tells that the quotient could be one bit less. Let's see how false it is.

Let's write that b have hb bits, a ha bits and q hq bits, and the conditions of Euclidean division:

(1<< (ha-1)) <= a, a < (1 << ha)
b < (1 << hb)
q < (1 << hq)
a=(b*q+r)
r < b

From the last two relations we have a < (b*q+b), that is a < (b * (q+1)).

We still have, q + 1 <= (1 << hq), thus b * (q+1) < ((1 < hb) * (1 << hq)), that is a < (1 < (hb+hq)).

But we also have (1 << (ha - 1)) <= a, thus ha - 1 < (hb + hq).

This can be written hq > (ha - hb - 1), or hq >= (ha - hb). So the quotient will have at least ha - hb bits, maybe more, but never less.

The third part in yellow is trying to perform an addition with bitAnd: - ouch ! It could eventually be bitOr:, but why not simply writing q := (q bitShift: n - q) + q1...

Fortunately, this is the case of one bit less which can never happen.

There are other things to be enhanced. The numbers 54 and 1022 are kind of magical and use implicit knowledge that 54 = (Float precision + 1) and -1022 = Float emin, the minimum possible exponent of a Float before gradual underflow occurs.

So we should better re-write this code and:

change the heading comment;
remove the unreachable branch;
use the new floating point inquiry messages precision and emin;
name the mantissa mantissa rather than q;
avoid to compute the remainder r, but just check if it is zero or hasTruncatedBits;
let the floating point hardware handle the case of negative zero by itself
(we just have to keep at least one bit in the mantissa).

asFloat
    "Answer a Float that closely approximates the value of the receiver.
    This implementation will answer the closest floating point number to the receiver.
    In case of a tie, it will use the IEEE 754 round to nearest even mode.
    In case of overflow, it will answer +/- Float infinity."

    | a b mantissa exponent hasTruncatedBits lostBit n ha hb hm |
    a := numerator abs.
    b := denominator.    "denominator is always positive"
    ha := a highBit.
    hb := b highBit.

    "Number of bits to keep in mantissa plus one to handle rounding."
    n := 1 + Float precision.

    "If both numerator and denominator are represented exactly in floating point number,
    then fastest thing to do is to use hardwired float division."
    (ha < n and: [hb < n]) ifTrue: [^numerator asFloat / denominator asFloat].

    "Shift the fraction by a power of two exponent so as to obtain a mantissa with n bits.
    First guess is rough, the mantissa might have n+1 bits."
    exponent := ha - hb - n.
    exponent >= 0
        ifTrue: [b := b bitShift: exponent]
        ifFalse: [a := a bitShift: exponent negated].
    mantissa := a quo: b.
    hasTruncatedBits := a > (mantissa * b).
    hm := mantissa highBit.

    "Check for gradual underflow, in which case the mantissa will loose bits.

    Keep at least 1 bit to let the underflow preserve the sign of zero."
    lostBit := Float emin - (exponent + hm - 1).
    lostBit > 0 ifTrue: [n := n - lostBit max: 1].

    "Remove excess bits in the mantissa."
    hm > n
        ifTrue:
            [exponent := exponent + hm - n.
            hasTruncatedBits := hasTruncatedBits or: [mantissa anyBitOfMagnitudeFrom: 1 to: hm - n].
            mantissa := mantissa bitShift: n - hm].

    "Check if mantissa must be rounded upward.
    The case of tie (mantissa odd & hasTruncatedBits not)
    will be handled by Integer>>asFloat."
    (hasTruncatedBits and: [mantissa odd])
        ifTrue: [mantissa := mantissa + 1].

    ^ (self positive
            ifTrue: [mantissa asFloat]
            ifFalse: [mantissa asFloat negated])
        timesTwoPower: exponent

That's it. No revolution, but I hope it's a bit clearer, and for the same price, a tiny bit more efficient.

I could also explain in this blog why it is better to use floating point division if both numerator and denominator can be converted asFloat exactly... Or why the mantissa loose bits in case of underflow... But, oh, I'm so lazy, ain't that boring enough?

Clarifying and optimizing Integer>>asFloat

2011-09-22T00:09:00.000+02:00

I wanted to explain some notions about floating point and thought using Integer>>asFloat in Squeak as example. Unfortunately, I found the code was both not that clear and not that efficient and some of the comments were inaccurate (I can tell so without any diplomacy since it's my own production).

So I decided to change it once again.

First, let's see the comment. Three lines should be enough. Note that I didn't mention nearest Float, because in case of overflow we will answer Float infinity, while the nearest Float would be Float fmax.

Integer>>asFloat
   "Answer a Float that best approximates the value of the receiver.
   This algorithm is optimized to process only the significant digits of a LargeInteger.
   And it does honour IEEE 754 round to nearest even mode in case of excess precision (see details below)."

Then, it might be useful to detail what round to nearest even means exactly. It's quite long and I'm not sure this belongs to the method or if it should just be replaced by a link... But it does also help understanding the code below.

   "How numbers are rounded in IEEE 754 default rounding mode:
   A shift is applied so that the highest 53 bits are placed before the floating point to form a mantissa.
   The trailing bits form the fraction part placed after the floating point.
   This fractional number must be rounded to the nearest integer.
   If fraction part is 2r0.1, exactly between two consecutive integers, there is a tie.
   The nearest even integer is chosen in this case.
   Examples (First 52bits of mantissa are omitted for brevity):
   2r0.00001 is rounded downward to 2r0
   2r1.00001 is rounded downward to 2r1
   2r0.1 is a tie and rounded to 2r0 (nearest even)
   2r1.1 is a tie and rounded to 2r10 (nearest even)
   2r0.10001 is rounded upward to 2r1
   2r1.10001 is rounded upward to 2r10
   Thus, if the next bit after floating point is 0, the mantissa is left unchanged.
   If next bit after floating point is 1, an odd mantissa is always rounded upper.
   An even mantissa is rounded upper only if the fraction part is not a tie."

Then we can come with a cleaner implementation by just naming the variables more explicitly to fit the comments:

| mantissa shift sum excess mask trailingBits nextBitIsSet tie |

mantissa := self abs.

"Check how many bits excess the maximum precision of a Float mantissa."

excess := mantissa highBit - Float precision.

excess > 0

ifTrue:

["Remove the excess bits"

mask := (1 bitShift: excess) - 1.

trailingBits := mantissa bitAnd: mask.

"trailingBits isZero ifFalse: [Inexact signal]."

mantissa := mantissa bitShift: excess negated.

shift := excess.

"But care to honour IEEE 754 round to nearest even mode."

nextBitIsSet := trailingBits highBit = excess.

(nextBitIsSet and: [mantissa odd or:

[tie := trailingBits isPowerOfTwo.

tie not]])

ifTrue: [mantissa := mantissa + 1]]

ifFalse: [shift := 0].

"Now that mantissa has no more excess precision, the following floating point operations will be exact."

sum := 0.0.

1 to: mantissa digitLength do:

[:byteIndex |

sum := sum + ((mantissa digitAt: byteIndex) asFloat timesTwoPower: shift)

shift := shift + 8].

^ self positive

ifTrue: [sum]

ifFalse: [sum negated]

Note that in this case, the rounding is performed entirely in Smalltalk and the algorithm will round to nearest even whatever the hardware rounding mode is, because all successive floating point operations will be exact.

The hardware inexact flag thus won't be set. We could replace this by a soft Inexact signal, but this is commented out because we don't care of and don't have such exception in Smalltalk until now. It's just a Smalltalkish comment.

But this code still performs many bit operations on LargeIntegers in case of excess precision. We could try alternatives that won't produce intermediate LargeIntegers.

excess > 0

ifTrue:

["Take care to honour IEEE 754 round to nearest even mode."

mustRoundUpper := (mantissa bitAt: excess) > 0 "The bit after floating point is set"

and: [(mantissa bitAt: 1 + excess) > 0 "The mantissa is odd"

or: [excess > 1 and: [self anyBitOfMagnitudeFrom: 1 to: excess - 1 "This is not a tie"]]].

"Remove the excess bits"

mantissa := mantissa bitShift: excess negated.

shift := excess.

mustRoundUpper ifTrue: [mantissa := mantissa + 1]]

But see how we must assist the code with many comments... Not that good. And not that faster...

We can still write more efficient code: we just have to let the hardware perform the rounding by itself. In effect, if we just add one excess bit, we are sure that there will be at most one inexact operation, and thus that we won't cumulate round off errors. But in this case, it is necessary to explain in preamble:

"Algorihm details:

Floating point hardware will correctly handle the rounding by itself with a single inexact operation if mantissa has one excess bit of precision.

Except in the last case when extra bits are present after an even mantissa, we must round upper by ourselves.

Note 1: the inexact flag in floating point hardware must not be trusted because it won't take into account the bits we truncated by ourselves.

Note 2: the floating point hardware is presumed configured in default rounding mode."

    | mantissa shift sum excess |
   mantissa := self abs.

   "Check how many bits excess the maximum precision of a Float mantissa."
   excess := mantissa highBit - Float precision.
   excess > 1
       ifTrue:
           ["Remove the excess bits but one"
           mantissa := mantissa bitShift: 1 - excess.
           shift := excess - 1.
           "Handle the case of extra bits truncated after an even mantissa."
           ((mantissa bitAnd: 2r11) = 2r01 and: [self anyBitOfMagnitudeFrom: 1 to: shift])
               ifTrue: [mantissa := mantissa + 1]]
       ifFalse: [shift := 0].

   "Now that mantissa has at most 1 excess bit of precision, let floating point operations perform the final rounding."
   sum := 0.0.
   1 to: mantissa digitLength do:
       [:byteIndex |
       sum := sum + ((mantissa digitAt: byteIndex) asFloat timesTwoPower: shift).
       shift := shift + 8].
   ^ self positive
       ifTrue: [sum]
       ifFalse: [sum negated]

This is shorter and faster, and a bit more fragile because now depending on hardware settings... But we never change these settings in Smalltalk, so we don't mind. Eventually, we could also come with a version that will always honour hardware rounding mode in the future - only if we can control it from within Smalltalk. The most annoying thing is that we trade some didactic properties now hidden into hardware black box. Now I must decide wether to commit immediately or wait for further comments...

ArbitraryPrecisionFloat

2011-08-30T02:35:00.000+02:00

After Smallapack, I decided to open another google code project for another already existing cross-dialect package - ArbitraryPrecisionFloat. This time, no FFI complications, the package is 100% "pure" Smalltalk (it does not mean 100% beautiful, I can't judge myself but severly). I will try to update the Wikipedia page because my favourite language deserves a bit more visibility.

Smallapack on google code

2011-08-26T00:26:00.001+02:00

Finally I decided to gather notes about Smallapack on google code https://code.google.com/p/smallapack
The site is very young and does not contain much material yet, but it will be completed with installation notes, code snippets and a few insights on design.

Smallapack progress on Squeak

2011-08-23T22:14:00.002+02:00

This post is about Smallapack, a Smalltalk interface to LAPACK.

Smallapack deals with basic linear algebra (the BLAS) or more ellaborated linear algebra, like solving least squares problem, eigenvalues problems, singular values decompositions etc... (LAPACK).

It enables Smalltalk to deal efficiently with number crunching by delegating these operations to specially optimized external libraries, at the express condition that operations can be arrayed (a lot like Matlab interpreter for example).

Smallapack implementation was working for a long time in Cincom VW and Dolphin dialects, but I never was able to have it running reliably on Squeak. Until yesterday night, when I learned the stupid little detail that caused that mess:

Modern C ABI expect an 8 byte alignment for double * data (especially for some library optimized thru SSE2 instructions or such).

Squeak has a super feature that enables passing a ByteArray allocated in Object memory rather than an ExternalAddress pointing to some memory allocated on external heap. I abused this feature, because heap management requires malloc and free, and why bother with so much basic tasks, when I can just let the garbage collector take care of Object memory... Unfortunately, this excess of trust was my mistake:

Squeak objects are 4 bytes aligned, at least on all 32 bit VMs available to date.

Hence, some function calls will randomly crash the VM with a protection fault (OSX gdb says EXC_BAD_ACCESS). The garbage collector has a license to allocate on 4 byte boundary or to rellocate on 4 bytes boundary at any time, and thus there is not much we can do about it. Eliot told me that a new garbage collector with proper 8 byte alignment was in his plans, but not soon, Once more, I will depend on the tremendous amount of work this guy put in Smalltalk, but who doesn't? Viva Miranda!

So my bug was happening by example with cblas_dgemv on Mac OSX veclib if the matrix is transposed, and by extension with any LAPACK function calling a transposing DGEMV...

GOOD NEWS! Once discovered last night, the bug was quickly fixed this evening, it was just a matter of redefining two methods (arrayPointer) so that they transfer the data in external heap (storeInCSpace) before answering the pointer object...

IN THE INTERIM, FOR THOSE WHO USE Squeak FFI: BEWARE.

It is very dangerous to use a (ByteArray new: byteSize) as a memory for holding a table of double and passing it as argument to an external procedure... Some external procedures will tolerate a 4 bytes alignment, some won't... You should always prefer an (ExternalAddress gcallocate: byteSize) for passing a double *. By virtue of libc malloc(), the ExternalAddress shall always contain a correctly aligned pointer.

Funny design decision

2011-08-16T23:58:00.000+02:00

I write this for sharing but also for reminding me later: if you ever have to develop in FORTRAN there is a gfortran option that SHOULD be your companion

-finit-real=snan

With this option, every real variable will automatically be initialized with a signalling nan, and your program should stop when it will incorrectly use the variable. The SUN FORTRAN compiler did have such option for long, but it takes time to reach this original quality.
And indeed, g77 did not provide such option but only a poor man's

-finit-local-zero

The former helps you finding your bugs while the later is more to help you hiding your bugs. What a funny design decision...

I had to code array/matrix crunching code for tight loops and wanted to use Lapack/Blas for that. And I wanted the latest version with bug corrections. However, I discovered that release 3.3 of Lapack now contains FORTRAN 90 specific instructions. Of course, there are plenty sources of errors in FORTRAN 77 that have been addressed by FORTRAN 90... The ones that would have been most useful to me today are:

verification of subroutine interface signature (number of arguments, types and input/output intent)
passing an object (with number of dimensions, bounds and stride) rather than just the address of the raw array storage.

The former was traditionally verified by 3rd party analyzers (forcheck, ftnchek).
But the later is a real improvement. There was a compiler option that did perform runtime checks

-fbounds-check

But it was more a joke than something useful because it uses the declared bounds, not the really allocated ones.
Unfortunately, if you want to profit by these features, your code has to adhere to FORTRAN 90 interface definitions. Lapack does not. It was written in FORTRAN 77 for historical reasons and I suspect it remains like this both for economical reasons and efficiency reasons (like passing an address is just faster than passing an array and going thru runtime checks). Since the FORTRAN 90 interface declarations wrappers to Lapack are not even kept up to date, there is not much to gain by using FORTRAN 90... That just prevents me to use g77, f2c, or ftnchek for no good reason (except a recursive routine), so I thought, what a funny design decision...

Finally, I went thru edit/compile/print (gdb will not debug FORTRAN that easily, especially if you want to inspect multi-dimensional arrays, you first have to debug the debugger). As usual, some error was on my test code, while I was stupidly focusing on the more complex tested code. At the end of the day I had to admit that most of the time lost was in:

translating a few f90-ism back to f77 (to please g77);
and not testing uninitialized variables (no such option in g77).

All this, just because I decided to use g77... What a funny design decision!

And next time, I SHALL remember that ftnchek SHOULD always be the natural companion of such funny decisions.

Well, we are very far from Smalltalk. Our debugger rarely has to be debugged and our instance and temp variables are initialized to nil, so design quality was again ahead and the subject is closed. But would Smalltalk perform your number crunching? Not mine. Not even, C, C++, or FORTRAN will. That's why the Blas were invented. I eventually have Smallapack, a functioning interface to Blas/Lapack in Visualworks, but it's not up to date either (Lapack 3.0...) and debugging bare bone allocation problems through an additional DLLCC/VM layer would not help. What I would need (and maybe program one day in Smalltalk) is a FORTRAN interpreter with all the necessary bound checks...

Lazy initialization of Shared Variable Bindings

2011-08-04T14:53:00.000+02:00

One problem we encounter with package management in Squeak/Pharo, but also in other Smalltalk dialects is the initialisation of shared variables.

Generally a class initialization message is defined like this:

MyClass class>>initialize

MyShared1 := SomeOtherClass new.

self snipLotOfMoreCode

And by convention, this initialize method is executed at package load.

But what if the initialization of a package require another package which is not yet loaded (for example containing SomeOtherClass or less trivially, containing a message sent during the initialization process).

One work around often encountered is to use a message indirection coupled with a lazy initialization to access the shared variable value, like this:

MyClass class>>someSharedInfo

^SomeSharedInfo ifNil: [self initializeSomeSharedInfo]

But why putting pressure on coder when the system could do it by itself ?

The idea is to create a new class for holding uninitialized variable bindings or value and associate an initialization method to be executed when first access will be attempted.

Associating an initialization method could be performed through use of pragmas, for example:

MyClass class>>initializeSomeSharedVariables

#<initialize>

MyShared1 := 5.

MyShared2 := #yourself.

In Smalltalk, a shared variable is a binding (an Association) with the variable name as key (generally a Symbol) and a value (the value of the variable). The Smalltalk Compiler arrange so that binding is unique across all method accessing the variable, thus the value can effectively be shared.

Now, what happens if we access MyShared1 ? In Smalltalk, shared variable are generally not accessed by sending value/value: messages to the binding. Instead, the compiler produces a byte code that directly fetch the second inst.var. for reading and another similar byte code for writing the value.

But this can be controlled in Squeak/Pharo by a compiler hook: just define your own binding class and two methods:

YourOwnSharedVariableBinding>>isSpecialReadBinding
"Return true if this variable binding is read protected, e.g., should not be accessed primitively but rather by sending #value messages"
^true

YourOwnSharedVariableBinding>>isSpecialWriteBinding
"Return true if this variable binding is write protected, e.g., should not be accessed primitively but rather by sending #value: messages"
^true

So we now have many options to implement a lazy initialization...

Solution 1) instead of initializing value with nil, create an UnitializedSharedValue class pointing to the binding and an initialization message to be performed

And implement doesNotUnderstand: like this

ProtoObject subclass: #UnitializedSharedValue
    instanceVariableNames: 'binding initializer'
    classVariableNames: ''
    poolDictionaries: ''
    category: 'Utility'

setInitializer: anEvaluator
"Set the evaluator responsible of initializing variable value.
anEvaluator should understand the message value."

initializer := anEvaluator

setBinding: aVariableBinding
binding := aVariableBinding

doesNotUnderstand: aMessage

initializer

ifNil:

[binding error: 'This variable binding does not have any known initializer'.

^nil].

initializer value.

binding value == self

ifTrue:

[binding error: 'The initializer failed to initialize this binding'.

^nil].

^binding value perform: aMessage selector withArguments: aMessage arguments

The problem with this kind of trick is that it is hard to catch all message sends (in Squeak/Pharo, class is not a true message send for example). And if you catch them, you then can't debug easily, because the debugger uses message sends too.

Solution 2) create AutoInitializedSharedVariableBinding that is Special Read Binding

and use an indirection to access the value:

LookupKey subclass: #AutoInitializedSharedVariableBinding
    instanceVariableNames: 'value undefined initializer'
    classVariableNames: ''
    poolDictionaries: ''
    category: 'Collections-Support'

initialize

value := undefined := Object new.

isVariableBinding
^true

isSpecialReadBinding
^true

setInitializer: anEvaluator

initializer := anEvaluator

value

^value == undefined

ifTrue: [value]

ifFalse: [self initializeValue]

initializeValue

initializer

ifNil:

[self error: 'This variable binding does not have any known initializer'.

^nil].

initializer value.

value == undefined

ifTrue:

[self error: 'The initializer failed to initialize this binding'.

^nil].

^value

This will cost an indirection (the #value send) and one test at each access. Not worse than manual lazy initialization we saw first.

Solution 3) avoid the value == undefined test by using a become: trick. There will still be a value/value: indirection cost though, so it might not be worth the complications...

LookupKey subclass: #UninitializedVariableBinding
    instanceVariableNames: 'value'
    classVariableNames: ''
    poolDictionaries: ''
    category: 'Collections-Support'

key: aKey value: anEvaluator
key := aKey.
value := anEvaluator

setInitializer: anEvaluator

value := anEvaluator

isVariableBinding
^true

isSpecialReadBinding
^true

isSpecialWriteBinding
^true

value
    ^value
        ifNil: [UnitializedVariableException signal]
        ifNotNil: [self value: value value]

value: anObject
self become: (InitializedVariableBinding key: key value: anObject).
^anObject

LookupKey subclass: #InitializedVariableBinding
   instanceVariableNames: 'value'
   classVariableNames: ''
   poolDictionaries: ''
   category: 'Collections-Support'

key: aKey value: anObject
key := aKey.
value := anObject

setInitializer: anEvaluator

self become: (UninitializedVariableBinding key: key value: anEvaluator).
^self

isVariableBinding
^true

isSpecialReadBinding
^true

isSpecialWriteBinding
^true

value
^value

value: anObject
^value := anObject

Note the return in last method so that we can chain MyShared1 := MyShared2 := someValue...

Solution 4) like in Cincom VisualWorks, let all SharedVariableBinding carry their own initialization code (which is like a clean block). But this does not solve the initialization order problems and this requires tools support (and also source code management support for storing/retrieving/executing the code associated to the shared variable binding).

The pragma handling common to all solutions:

Just create this class

InstructionClient subclass: #SharedInitializationFinder
    instanceVariableNames: 'usedBeforeDefined initializedBindingsByFetchingValue initializedBindingsBySendingValue stack stackPerPC pc'
    classVariableNames: ''
    poolDictionaries: ''
    category: 'Utility'`

With these class methods:

initialize
"SharedInitializationFinder initialize."
self registerForEvents

registerForEvents
SystemChangeNotifier uniqueInstance noMoreNotificationsFor: self.
SystemChangeNotifier uniqueInstance notify: self ofAllSystemChangesUsing: #initializeEvent:.

scanMethod: aCompiledMethod
"Answer with a list of variable bindings initialized by aCompiledMethod via sending #value: message."
^self new scanMethod: aCompiledMethod

initializeEvent: anEvent
    "Check if this system event defines or removes an automatic initialization."
    | aClass aSelector method |
    (anEvent itemKind = SystemChangeNotifier classKind and: [anEvent isRemoved])
        ifTrue: ["We should track initializations bound to this class..."].
    anEvent itemKind = SystemChangeNotifier methodKind ifTrue: [
        aClass := anEvent itemClass.
        aClass isMeta ifFalse: [^self]. "ignore instance methods"
        aClass := aClass theNonMetaClass.
        aSelector := anEvent itemSelector.
        aSelector numArgs > 0 ifTrue: [^self]. "ignore methods with arguments"
        (anEvent isRemoved) ifTrue: [
            "We should track initializations bound to this class and method..."].
        (anEvent isAdded or: [anEvent isModified]) ifTrue: [
            method := anEvent item.
            method pragmas do: [:pragma |
                | shared |
                pragma keyword == #initialize ifTrue: [
                    shared := self scanMethod: method.
                    shared do: [:binding | binding setInitializer: (MessageSend receiver: aClass selector: aSelector)]]]]].

On instance side, let the class process the bytecodes decoding (most are omitted), the most important being:

initializeBindings
    usedBeforeDefined := Set new.
    initializedBindingsByFetchingValue := Set new.
    initializedBindingsBySendingValue := Set new.

initializedBindingsByFetchingValue
    "Answer with the list of variable bindings being initialized by target CompiledMethod by directly storing into value."
    ^initializedBindingsByFetchingValue
initializedBindingsBySendingValue
    "Answer with the list of variable bindings being initialized by target CompiledMethod by directly sending value: message."
    ^initializedBindingsBySendingValue

pushLiteralVariable: anAssociation
"Push Contents Of anAssociation On Top Of Stack bytecode."
stack addLast: anAssociation

popIntoLiteralVariable: anAssociation
    "Remove Top Of Stack And Store Into Literal Variable bytecode."
    initializedBindingsByFetchingValue add: anAssociation.
    stack removeLast

send: aSelector super: supered numArgs: numberArguments
    | theReceiver |
    1 to: numberArguments do: [:i | stack removeLast].
    theReceiver := stack removeLast.
    theReceiver ifNotNil: [
        aSelector == #value ifTrue: [usedBeforeDefined add: theReceiver].
        aSelector == #value: ifTrue: ["Protect against trivial infinite loops..."
            (usedBeforeDefined includes: theReceiver)
                ifFalse: [initializedBindingsBySendingValue add: theReceiver]]].
    stack addLast: nil

scanMethod: method
    | scanner end |
    self initializePC.
    scanner := InstructionStream on: method.
    end := method endPC.
    [(pc := scanner pc) <= end] whileTrue:
        [stackPerPC at: pc ifPresent: [:restoredStack | stack := restoredStack].
        scanner interpretNextInstructionFor: self].
    ^self initializedBindingsBySendingValue

Last, we need to define a hook, such that any binding could become an uninitialized or auto initialized variable binding, for example, solution 3) would be:

LookUpKey>>setInitializer: anEvaluator
    "Set the evaluator responsible of initializing variable value.
    anEvaluator should understand the message value."
    self become: (UninitializedVariableBinding key: key value: anEvaluator).
    ^self

The nice thing with above hook is that each time you change the initialization method, the variables will be reset by the setInitializer: mechanism. No more manual doIt : self initializeThisOrThat.

Of course, there are more things to handle before getting a fully functioning framework:

1) what if you remove an initializer method/class ?

2) what if you define two or more initializer methods ?

Especially, if you define a second initialization method then, oops, remove it, the variables should switch back to the first initialization...

3) detecting unchanged initialization method and avoiding unecessary re-initialization would be smart.

But these are boring little details... (where the devil is hidden as you may know).

Xtreams version of the optimized prime factors factorial

2011-07-12T14:15:00.000+02:00

I now want to use Xtreams on the optimized version of primeFactorFactorial:

optimizedPrimeFactorFactorialViaXtreams

| primes uniqPowers highestPower prevPower primeGroup primesGroupedByPower bitIterator terms numGroupsEstimate |

Forming a generator of primes:

primes := [:out | self orLessAllPrimesDo: [:prime | out put: prime ] ] reading.

Removing first prime, because we will handle power of 2 by a final bitShift:

primes get.

Forming two collections, the group of primes having the same power, and their unique powers.

   numGroupsEstimate := 1 << (self highBit // 2).
   uniqPowers := (Array new: numGroupsEstimate) writing.
   primesGroupedByPower := (Array new: numGroupsEstimate) writing.
   prevPower := nil.
    primeGroup := Array new writing.
   primes do: [:p |
       | primePower |
       primePower := self - (self sumDigitsInBase: p) // (p - 1).
       prevPower
           ifNil:
               [uniqPowers put: primePower.
               prevPower := highestPower := primePower]
           ifNotNil: [prevPower ~= (prevPower := primePower)
               ifTrue:
                   [primesGroupedByPower put: primeGroup conclusion reading.
                   uniqPowers put: primePower.
                    primeGroup := Array new writing]].
        primeGroup put: p].
   primesGroupedByPower put: primeGroup conclusion reading.
   primesGroupedByPower := primesGroupedByPower conclusion reading.
   uniqPowers := uniqPowers conclusion reading.

Note that this part is not lazy, because I later want to reuse these two collections several times. This part is not the most elegant I came up with. It is much like using a (ending:) slicing sequence, but ending: is not splitting the groups where I want to...

Another option would be to count the number of equal powers, and use that in a (limiting:) slicing sequence, but limiting: is not taking a block argument...

Forming an iterator on the powers bits:

    bitIterator := [:out |
       | bitProbe |
       bitProbe := 1.
       [ bitProbe <= highestPower ]
           whileTrue:
           [out put: bitProbe.
           bitProbe := bitProbe bitShift: 1 ] ] reading.

For each bit, grouping primes having the bit set in their associated power, and raise the product of this group to the power of corresponding bit:

    terms := bitIterator collect: [:bit |
        primesGroupedByPower += 0.
        uniqPowers += 0.
        (primesGroupedByPower selecting: [:pp |
            pp += 0.
            (uniqPowers get bitAnd: bit) ~= 0])
                stitching rest karatsubaProduct raisedToInteger: bit].

Note that stitching is concatenating the substreams into a single stream, while += 0 is resetting the streams at the beginning.

Performing the product of products and shifting to account the power of 2:

^terms karatsubaProduct bitShift: self - self bitCount

There is still an operation I didn't Xtreamize, the karatsubaProduct. It consists in a recursive split. Not sure how to do it...

Now let us see if timing is decent:

Smalltalk garbageCollect.
[50000 optimizedPrimeFactorFactorialViaXtreams] timeToRun

-> 1321

Approximately the same as the non Xtreams version (+4%). That's pretty good.

But this version has a tricky bug: stitching is raising Incomplete if there is no substream, and this happens for some numbers. For example,

9 optimizedPrimeFactorFactorialViaXtreams

produces the groups ((3) (5 7)) with the powers (4 1). When the bitIterator generator is evaluating with bit = 2, there is no group having a power with this bit set, Incomplete is raised by stitching, and caught by the generator which closes the stream without the final term (3 raisedTo: 4). Thus it returns a wrong result... Ouch! That's one reason for I dislike Exceptions.

My single test on a huge number did not fail, and I only discovered the bug with the generator free version:

optimizedPrimeFactorFactorialViaXtreams2
    | uniqPowers highestPower prevPower primeGroup primesGroupedByPower terms numGroupsEstimate |
    numGroupsEstimate := 1 << (self highBit >> 1).
    uniqPowers := (Array new: numGroupsEstimate) writing.
    primesGroupedByPower := (Array new: numGroupsEstimate) writing.
    prevPower := nil.
    highestPower := 0.
    primeGroup := Array new writing.
    Integer primesUpTo: self + 1 do: [:p | p = 2 ifFalse: [
        | primePower |
        primePower := self - (self sumDigitsInBase: p) // (p - 1).
        prevPower
            ifNil:
                [uniqPowers put: primePower.
                prevPower := highestPower := primePower]
            ifNotNil: [prevPower ~= (prevPower := primePower)
                ifTrue:
                    [primesGroupedByPower put: primeGroup conclusion reading.
                    uniqPowers put: primePower.
                    primeGroup := Array new writing]].
        primeGroup put: p]].
    primesGroupedByPower put: primeGroup conclusion reading.
    primesGroupedByPower := primesGroupedByPower conclusion reading.
    uniqPowers := uniqPowers conclusion reading.
    terms := (1 to: highestPower highBit) collect: [:bit |
        primesGroupedByPower += 0.
        uniqPowers += 0.
        [| product |
        product := (primesGroupedByPower selecting: [:pp |
            pp += 0.
            (uniqPowers get bitAt: bit) ~= 0])
                stitching rest karatsubaProduct.
        1 to: bit - 1 do: [:i | product := product squared].
        product] on: Incomplete do: [1]].
    ^terms karatsubaProduct bitShift: self - self bitCount

Overall, I find Xtreams powerfull, far more than Smalltalk Stream. But there are still things lacking. For example, the grouping was not obvious. The stitching of empty streams was surprising. I might need a readingStitching that works with an empty stream.

Xtreams is interesting both for its functional like aspects and for its potential use in parallelism, a lot like this. From this later point of view, I find it easy to pipe serial operations with Xtreams in case of a single assembly line. It also seems possible to fork processing for each element when they are independent (this is data parallelism - and could be applied for performing karatsubaProduct of each term in above example, and thus get a very efficient parallel version of factorial).
On the other hand, having two parallel assembly lines does not seem that easy. There is a duplicating: operation that should allow such parallel streams, but it necessarily uses a writeStream. I'd like the assembly chain to be as easy as:

fork := [:in :out1 :out2 | | item | item := in get. out1 put: item. out2 put: item] forking.
line1 := [:in :out | out put: in get op1] filtering.

line2 := [:in :out | out put: in get op2] filtering.

join := [:in1 :in2 :out | out put: (in1 get op3: in2 get)] joining.
factory := generator | fork | (line1 & line2) | join.

Funily, I used unix pipe | to denote serial and ampersand & to denote parallel, but I find the contrary more natural: || makes me think of parallelism, and & of sequencing. Just a point of view...

From above point of view, the differences between read and write streams seems a bit arbitrary. Filters are inherently both. The Generator pattern could be used to turn any WriteStream (push) into a ReadStream (pull), but it's current implementation has potential dead lock conditions. But wait, it's possible I missed some features, or miss-used some others, Xtreams is evolving rapidly. There is a new post about using Xtreams with parallel VM, I have to inquire a bit more!

Testing new Xtreams generators

2011-07-09T02:14:00.000+02:00

So, I upgraded Squeak port of Xtreams, and just tested it on the primeFactorFactorial example:

primeFactorFactorialViaXtreams
    | primes powers factorial |
    primes := [:out | Integer primesUpTo: self + 1 do: [:prime | out put: prime ] ] reading.
    powers := primes collecting: [:p | p raisedToInteger: self - (self sumDigitsInBase: p) // (p - 1)].
    factorial := powers injecting: 1 into: [:product :power | product * power].
    ^(factorial positioning buffer: (XTRingBuffer new: 1 class: Array))
        -= 1;
        get

This looks good, except the final sentence which is not simple to read to my taste, remember it just extract the last element...

Now let's see performance:

Smalltalk garbageCollect.
Time millisecondsToRun: [10000 primeFactorFactorial].
-> 297

Smalltalk garbageCollect.
Time millisecondsToRun: [10000 primeFactorFactorialViaGeneratorAndXtream].
-> 390

Smalltalk garbageCollect.
Time millisecondsToRun: [10000 primeFactorFactorialViaXtreams].
-> 945

Huh?

MessageTally spyOn: [10000 primeFactorFactorialViaXtreams].

- 820 tallies, 822 msec.

**Tree**
--------------------------------
Process: other processes
--------------------------------
47.8% {393ms} WeakArray class>>finalizationProcess
39.4% {324ms} WeakFinalizationList class>>initTestPair
    |39.4% {324ms} Object class(Behavior)>>new
8.4% {69ms} primitives
17.0% {140ms} EventSensor>>eventTickler
17.0% {140ms} Delay>>wait
--------------------------------
Process: (40s) 86874: nil
--------------------------------
34.0% {279ms} XTPositionReadStream(XTReadStream)>>-=
34.0% {279ms} XTPositionReadStream>>length
    34.0% {279ms} XTPositionReadStream>>++
      34.0% {279ms} XTPositionReadStream(XTReadStream)>>++
        34.0% {279ms} XTPositionReadStream>>read:into:at:
          34.0% {279ms} XTCollectReadStream>>read:into:at:
            34.0% {279ms} XTCollectReadStream>>directRead:into:at:

**Leaves**
39.4% {324ms} Object class(Behavior)>>new
34.0% {279ms} XTCollectReadStream>>directRead:into:at:
17.0% {140ms} Delay>>wait
8.4% {69ms} WeakArray class>>finalizationProcess

**Memory**
    old            -1,959,252 bytes
    young        +1,059,496 bytes
    used        -899,756 bytes
    free        -1,059,496 bytes

**GCs**
    full            3 totalling 455ms (55.0% uptime), avg 152.0ms
    incr        19 totalling 65ms (8.0% uptime), avg 3.0ms
    tenures        0
    root table    0 overflows

MessageTally is not reliable with COG VM, but clearly, only one third of execution time is spent in the primeFactorFactorialViaXtreams. It seems like something is wrong with this image, and implementing the generator via interprocess synchronisation did favour those CPU hijackers!

Optimizing the product of prime powers

2011-07-06T06:00:00.106+02:00

Remember how we generated the factorial from the product of its prime factors?

Let's try and optimize the product.

First trick is to use a Karatsuba multiplication, because multiplication is the major contributor in execution time.

The Squeak variant takes care to evaluate Karatsuba only on well balanced operands, and replace the recursion by an iteration if not well balanced:

LargePositiveInteger>>karatsubaTimes: anInteger
    "eventually use Karatsuba algorithm to perform the multiplication"

    | half xHigh xLow yHigh yLow low high mid xLen yLen |
    (anInteger isLargeEnoughForKaratsuba
        and: [self isLargeEnoughForKaratsuba]) ifFalse: [^self timesInteger: anInteger].

    "Check if length ratio is more than 2, and engage a loop
    to operate on integers with well balanced lengths.
    Note that we only add overhead at this level,
    but we hope to gain in lower level recursion"
    (xLen := self digitLength) >= (yLen := anInteger digitLength)
        ifTrue: [(half := xLen bitShift: -1) >= yLen
            ifTrue: [^(0 to: xLen by: yLen) detectSum: [:yShift |
                ((self copyDigitsFrom: yShift + 1 to: yShift + yLen)
                    karatsubaTimes: anInteger)
                        bitShift: 8 * yShift]]]
        ifFalse: [(half := yLen bitShift: -1) >= xLen
            ifTrue: [^(0 to: yLen by: xLen) detectSum: [:xShift |
                (self karatsubaTimes:
                    (anInteger copyDigitsFrom: xShift + 1 to: xShift + xLen))
                        bitShift: 8 * xShift]]].

    "At this point, lengths are well balanced, divide each integer in two halves"
    xHigh := self bitShift: -8 * half.
    xLow := self lowestNDigits: half.
    yHigh := anInteger bitShift: -8 * half.
    yLow := anInteger lowestNDigits: half.

    "Karatsuba trick: perform with 3 multiplications instead of 4"
    low := xLow karatsubaTimes: yLow.
    high := xHigh karatsubaTimes: yHigh.
    mid := (xHigh + xLow karatsubaTimes: yHigh + yLow) - (low + high).

    "Sum the parts of decomposition"
    ^low + (mid bitShift: 8*half) + (high bitShift: 16*half)

With a bit of tuning for the COG VM:

LargePositiveInteger>>isLargeEnoughForKaratsuba
^self digitLength >= 160

SmallInteger is not large enough for Karatsuba and falls back to regular multiplication *, we omit the code here.
The copyDigitsFrom:to: and lowestNDigits: are just using a primitive for splitting a LargePositiveInteger, they are omitted too.

We'd also better optimize squared which is a simple case:

LargePositiveInteger>>squared
    "Eventually use a divide and conquer algorithm to perform the multiplication"

    | half xHigh xLow low high mid |
    self isLargeEnoughForKaratsuba ifFalse: [^self * self].

    "Divide digits in two halves"
    half := self digitLength bitShift: -1.
    xHigh := self bitShift: -8 * half.
    xLow := self lowestNDigits: half.

    "Use Karatsuba"
    low := xLow squared.
    high := xHigh squared.
    mid := xLow karatsubaTimes: xHigh.

    "Sum the parts of decomposition"
    ^low + (mid bitShift: 8*half+1) + (high bitShift: 16*half)

Then use squared in raisedToInteger:

Number>>raisedToInteger: anInteger
    | bitProbe result |

    anInteger negative ifTrue: [^(self raisedToInteger: anInteger negated) reciprocal].
    bitProbe := 1 bitShift: anInteger highBit - 1.
    result := self class one.
   [
        (anInteger bitAnd: bitProbe) = 0 ifFalse: [result := result * self].
        bitProbe := bitProbe bitShift: -1.
        bitProbe > 0 ]
    whileTrue: [result := result squared].

    ^result

We have a reasonable implementation for Karatsuba, but we cannot use a naïve:

Integer>>karatsubaPrimeFactorFactorial

    "Recompose the factorial from the prime factors, knowing the power of each prime."
    ^((Integer primesUpTo: self + 1) collect: [:p |
        p raisedToInteger: self - (self sumDigitsInBase: p) // (p - 1)]) karatsubaProduct

Collection>>karatsubaProduct

^self inject: 1 into: [:product :element | product karatsubaTimes: element]

Karatsuba is just a drag in this case because we always multiply a large integer with a small one. The only advantage comes from using squared in raisedToInteger:, and instead of 19136ms, we get:

Smalltalk garbageCollect.
Time millisecondsToRun: [1 to: 3000 do: [:x | x karatsubaPrimeFactorFactorial]].
-> 18408 18503

So the second trick is essential: evaluate the product of terms in a divide and conquer fashion, dividing the numbers to multiply in two groups and recursively:

SequenceableCollection>>karatsubaProduct

    "Compute the product of self elements, using Karatsuba multiplication"
    self isEmpty ifTrue: [^1].
    ^self karatsubaProductFrom: 1 to: self size by: 1

SequenceableCollection>>karatsubaProductFrom: startIndex to: stopIndex by: inc
    | nextInc nextIndex |
    (nextIndex := startIndex + inc) > stopIndex ifTrue: [^self at: startIndex].
    nextInc := inc * 2.
    ^(self karatsubaProductFrom: startIndex to: stopIndex by: nextInc) karatsubaTimes:
      (self karatsubaProductFrom: nextIndex to: stopIndex by: nextInc)

Smalltalk garbageCollect.
Time millisecondsToRun: [1 to: 3000 do: [:x | x karatsubaPrimeFactorFactorial]].
-> 12627 12539

Better, but there is room toward the prime swing variant which performs the job in less than 7000ms. Let's adapt our strategy and use a third trick: try and square as many LargeInteger as possible. One idea is to group the terms having same powers, multiply them together, then raise to the prescribed power. But if term p is raised to the power of 7, we will have to evaluate (p squared squared) * p squared * p. So we will store p into three collections,

the collection of terms raisedTo: 1;
the collection of terms raisedTo: 2;
the collection of terms raisedTo: 4.

The last point, is that we won't bother with powers of 2, because the fastest variants also don't. Here we come with this algorithm:

Integer>>optimizedPrimeFactorFactorial
    "This is the optimized version of primeFactorFactorial"

    | powers primes |
    self <= 1 ifTrue: [^1].
    primes := (Integer primesUpTo: self + 1) allButFirst.
    powers := primes collect: [:e | self - (self sumDigitsInBase: e) // (e - 1)].
    ^(primes karatsubaProductPowers: powers) bitShift: self - self bitCount

SequenceableCollection>>karatsubaProductPowers: powers

    "Compute the product of self elements, each raised to the corresponding powers"
    | lastPower rank terms |
    (lastPower := powers size) = 0 ifTrue: [^1].
    rank := 1.
    terms := OrderedCollection new: (powers at: 1) highBit.
    [rank <= (powers at: 1)]
        whileTrue:
            [ [ rank > (powers at: lastPower) ] whileTrue: [lastPower := lastPower - 1].
            terms add: ((((1 to: lastPower) select: [:e | ((powers at: e) bitAnd: rank) ~= 0]) collect: [:i | self at: i]) karatsubaProduct raisedToInteger: rank).
            rank := rank bitShift: 1].
    ^terms karatsubaProduct

Smalltalk garbageCollect.
Time millisecondsToRun: [1 to: 3000 do: [:x | x optimizedPrimeFactorFactorial]].
-> 8436 8433

The optimized prime factor algorithm is about 2.25 times faster than the naive prime factors on small integers up to 3000.

Let's instrument the algorithm to check if the terms are well balanced. We add the line
(terms collect: [:e | e highBit]) inspect.
and check 50000 optimizedPrimeFactorFactorial

-> an OrderedCollection(49775 49532 49490 49130 49165 48285 49407 48651 41929 40562 36299 26336 38838 55004 25969)

It's remarkable, that we decomposed the factorial of 50000 into terms having about 50000 bits!

The score on a larger integer :

[50000 factorial] timeToRun. -> 10288
[50000 primeFactorFactorial] timeToRun. -> 8853

[50000 optimizedPrimeFactorFactorial] timeToRun. -> 1263
[50000 primeSwingFactorial] timeToRun. -> 1004

Now a challenge: write optimizedPrimeFactorFactorial with Xtreams instead of SequenceableCollection!

Xtreams and Generators

2011-07-05T22:11:00.000+02:00

In my last post I reported the lack of Generator in Xtreams framework.
Well, I just received this clarification from Michael and Martin, the authors of Xtreams.
Definitely, Xtreams has a lot of potential!

Using generators on the primeFactorsFactorial example

2011-07-04T23:59:00.000+02:00

Last example algorithm, primeFactorsFactorial did involve enumerating of several collections:

The collection of primes up to n
The collection of powers of these primes

An alternative is to work with a stream that generates each item on the fly, rather than collecting the whole sequence of values, and then pipe to other kind of streams that transform the values. This technique has two advantages:

avoiding allocation of large auxiliary collections that will later be thrown out, can help some algorithms to scale up.
piping could theoretically enable parallel execution if we had a multi-threaded VM (though in fact we keep some collections, the optimal performance being often reached by tuning some buffers size, since the cost of inter thread synchronisation is not null)

This does not really apply to our prime factors decomposition since the collection of primes are not huge, but we have this example available, let's use it and see how it would be possible in Squeak.

First, Squeak provides an primesUpTo:do: taking a block as argument. That's all we need.

We can transform this block into a stream using a Generator. Let see in class comment how to use it:

    | generator |
    generator := Generator on: [:g| Integer primesUpTo: 100 do:[:prime| g yield: prime]].
    [generator atEnd] whileFalse:[Transcript show: generator next].

OK, so a new Generator is associated to an outer block, [:g| Integer primesUpTo: 100 do:[:prime| g yield: prime]], which tells which method to execute to produce the sequence of values.

The inner block [:prime| g yield: prime] just tells the generator that a new prime is available.

Each time we ask for generator next, the generator machinery arranges to advance the execution of outer block until it reaches execution of the inner block, and then switch back and return the prime provided to this block.

This is known as the coroutine technique, and implemented by using a clever a simple hack manipulating the call stack. Look at the implementation by yourself.

Note that a Generator now understand value: and can be used directly in place of the inner block, try (Generator on: [:g| Integer primesUpTo: 100 do: g]) contents.

Can we apply this to our primeFactorFactorial?

primeFactorFactorialViaGenerator
   | generator |
   generator := Generator on: [:g | Integer primesUpTo: self + 1 do: g].
   ^ (generator collect: [:p | p raisedToInteger: self - (self sumDigitsInBase: p) // (p - 1)]) product

Hmm, it will not work, because a Generator understands only one iteration operator, do:, but not collect: Squeak streams are a pity.

More over, the collect: operation if implemented by means of do: would not be lazy but rather produce the whole collection of powers. We would have only eliminated the first collection of primes, but not the collection of powers of primes.

One good and recent framework - Xtream - It is an interesting alternative and provides the lazy operations we're looking for; follow the squeaksource link to get download instructions. Let's see how we can lazily pipe the operations through streams:

primeFactorFactorialViaGeneratorAndXtream
    | generator stream factorial |
    generator := Generator on: [:g | Integer primesUpTo: self + 1 do: g].
    stream := ((XTBlockClosureReadStream
                on: [generator atEnd
                        ifTrue: [Incomplete raise]
                        ifFalse: [generator next]])
                collecting: [:p | p raisedToInteger: self - (self sumDigitsInBase: p) // (p - 1)])
                injecting: 1 into: [:product :power | product * power].
    ^ [[factorial := stream get] repeat]
        on: Incomplete
        do: [factorial]

Unfortunately Xtream has no Generator, so we have to wrap it up in a BlockClosureReadStream.

Another problem is that injecting:into: is not reducing the stream to a single element, instead it answers a stream of all the intermediate results. Since we only need the last one, we have to repeat the operation by ourself, capture and answer the last element when the stream starves (we catch this Incomplete exception).

Once I proposed an alternative selector ^[stream get] repeatUntil: Incomplete, with a slight modification, this would have been useful; or maybe we could just extend XTReadStream to respond to last.

Let see how the generator variant performs:

Smalltalk garbageCollect.
Time millisecondsToRun: [1 to: 3000 do: [:x | x primeFactorFactorialViaGeneratorAndXtream] ]
-> 22835 22529

This is to be compared to the 19136ms of naive form. Well, that was expected, we eliminated the intermediate collections, but they did not cost that much in this case. Instead we pay a price for indirection and maybe also the Generator trick forces COG VM to un-optimize some portions. But that doesn't matter, we saw some smart streaming techniques.

A more difficult exercise would be to transform a divide and conquer recursive product evaluation variant into a streaming equivalence... Yet I do not know how, but if I find it, this will be another post.

Evaluating a factorial from its prime factors

2011-07-03T06:00:00.032+02:00

Again, I will evaluate a suboptimal algorithm because it is nice and fun in Smalltalk.

I guess in other languages too, especially functional languages, but you know my favourite path.

The exploration started with this little factorial contest submitted by Paul Baumann.

We already know the fastest algorithm are the prime swing variant (read this wiki to get an explanation of swing by Peter Luschny), and simple and clever split recursive algorithm submitted by John Brant.

In both cases, only the odd products are performed, in a divide and conquer fashion which produces well balanced large integers (this is important for accelerated products like Karatsuba Toom or FFT). The final power of two being replaced by a bitShift: which saves many costly LargeInteger multiplications.

Our suboptimal variant of interest consists in evaluating the factorial from its decomposition in prime factors. We know that each prime <= n will be compose n factorial, but we also know the power of each of these primes thanks to Legendre. This algorithm can be naively transcribed in Smalltalk:

Integer>>primeFactorFactorial
    "Recompose the factorial from the prime factors, knowing the power of each prime."
    ^((Integer primesUpTo: self + 1) collect: [:p |
        p raisedToInteger: self - (self sumDigitsInBase: p) // (p - 1)]) product

primesUpTo: is a method found in Squeak which use a sieve of Eratosthenes, up to means not inclusive, so we need the disgracious + 1. Then Squeak lacks two methods, sumDigitsInBase: and product.

A naive product, could be simply:

Collection>>product

^self inject: 1 into: [:product :element | product * element]

With this definition, Collection new product is answering 1, and we don't need any protection for evaluating 0 primeFactorFactorial.

For sure, it won't be efficient and should better be replaced by a divide and conquer variant, but we'll see this in another post.

A naive sumDigitsInBase: could be:

Integer>>sumDigitsInBase: base

^(self printStringBase: base) inject: 0 into: [:sum :digit | sum + digit digitValue]

OK, we're done...

1000 primeFactorFactorial. -> Error: subscript is out of bounds: 303

Oops! Unfortunately, it's not only suboptimal, it won't work for primes above 36 because Squeak has no idea which character to use to represent digit 37. Anyway, even with unicode we would only have a limited set of characters. Huge, but limited.

It does not matter anyway, we can just reimplement printStringBase: and instead of streaming out digits converted to characters, just sum the digits.

A naive algorithm would be to just remove the last digit (the remainder by a division by base):

Integer>>sumDigitsInBase: base
"Answer the sum

| sum |

sum := 0.

q := self abs.

[q = 0]

whileFalse:

[sum := sum + (q - ((q := q // base) * base))].

^sum

Since we probably won't compute the factorial of a LargeInteger, this shall be enough.

For LargeInteger, we would use a divide and conquer method in order to avoid repeated LargeInteger divisions, and replace them by as many SmallInteger divisions as possible. This would be the exact image of printStringBase:, so we won't develop it here.

Now we are done with this very naive implementation. Let's see how it performs on COG VM:

Smalltalk garbageCollect.
Time millisecondsToRun: [1 to: 3000 do: [:x | x factorial]].
-> 29753 29780

Smalltalk garbageCollect.
Time millisecondsToRun: [1 to: 3000 do: [:x | x primeFactorFactorial]].
-> 19204 19136

Incredible, with this very naive code we already gain 33%. In fact, any implementation will probably beat factorial performance, but not its simplicity.

Revisiting the sieve of Eratosthenes

2011-07-02T00:58:00.002+02:00

Welcome, this is my first blogger blogging about Smalltalk, and its Squeak flavour.

Today, I will revisit the sieve of Eratosthenes for fun.

I don't think I need to introduce this famous algorithm.

The first idea is to use integers as arrays of bits, because it is fun.

If p is a composite integer, then the bit of rank p will be set to 1.

The second idea is to use an integer division to produce the repeated bit pattern using this property:

(2r11111111111 // 2r11) printStringBase: 2. -> '1010101010'

(2r11111111111111 // 2r111) printStringBase: 2. -> '100100100100'

(2r1111111111111111111 // 2r1111) printStringBase: 2. -> '1000100010001000'

OK, I think you get the idea of division. Now we need two things:

get the right number of 1 left of //
remove the trailing 1 in the quotient

We can do both with bitShift: or using the binary selector <<

If p is prime, and we want all primes up to n, then the adequate bit pattern is produced by:

1 << (n // p - 1 * p) - 1 // (1 << p - 1) << (p << 1 - 1).

Mind the parenthesis, above expression is a bit tricky...

Remember, Smalltalk binary operators read from left to right, all with same precedence.

First operation produces the 1s left of //:

1 << (n // p - 1 * p) - 1

Indeed, we need a length which is multiple of p in order to get a trailing 1 in the quotient

And we need enough bits to obtain a quotient with a leading 1 of rank <= n - (2*p-1)

Second operation is the division by appropriate number of 1s (p):

// (1 << p - 1)

Third operation shift by 2*p - 1 to obtain the trailing 0s:

<< (p << 1 - 1)

Indeed, we need p - 1 to align the 1s on rank multiples of p, + p to clear the trailing 1.

Of course, I can propose a few variants, for example a single 1 is enough left of //:

1 << (n // p - 1 * p) // (1 << p - 1) << (p << 1 - 1).

We can just complete the missing trailing 0s:

1 << (n - p) // (1 << p - 1) << (p << 1 - 1 - (n \\ p)).

Or arrange to have first p - 1 trailing 0s produced by the division:

1 << (n // p * p - 1) // (1 << p - 1) << p.

The last one looks good, we are ready to produce the composite mask up to n:

p := 1.

composite := 1.

[ p <= n ]

whileTrue:

[(composite bitAt: p) = 0

ifTrue:

["p is prime, mark all multiple of p up to n as composite"

composite := composite bitOr: 1 << (n // p * p - 1) // (1 << p - 1) << p ].

p := p + 1].

Of course we could stop the loop at n sqrtFloor, increment by 2 from second loop etc...

We could also enumerate primes by feeding a block with p inside the ifTrue: [ ].

But let's have fun by enumerating primes from the composite mask with bit tricks.

For example, we want to collect primes in an Array.

To know the number of primes up to n, we just need to count the number of bits set to 0.

Or subtract n with number of bit set to 1.

numberOfPrimes := n - composite bitCount.

Then, we just enumerate by finding and removing next bit set to zero, starting at least significant bit.

(Array new: numberOfPrimes) collect: [:void |

| isolatedBit |

"find next zero in composite"

isolatedBit := composite bitXor: composite + 1.

"remove next zero in composite"

composite := composite bitOr: isolatedBit.

isolatedBit highBit]

Indeed, 2rxxxxx01111 + 1 will generate a carry up to trailing 0, and result in 2rxxxxx10000, the leading xxxxx being unchanged.

The bitXor: then will remove the leading xxxxx, and result in 2r11111.

We then just need to get the rank of the highest bit, which is highBit's job...

An alternative is to first reverse the bits of composite mask.

Here we will just print the primes in the Transcript:

composite := 1 << n - 1 - composite.

[composite > 0] whileTrue:

[Transcript cr; show: composite lowBit.

composite := composite bitAnd: composite - 1]

Same here, 2rxxxxx10000 - 1 will generate a carry up to trailing 1 and result in 2rxxxxx01111, the leading xxxxx being unchanged. Then bitAnd: will just clear all but the leading xxxxx.

Of course, this method is not really efficient because each bitOr: and bitXor: operation will allocate a new LargePositiveInteger, not counting all the bitShift:, additions and subtractions, but do we care? Squeak Integer class>>primesUpTo: is already fast, we don't need another one...

Nastily, I removed the comments, obfuscated the names, and factored out repeated code with additional temps. It's not a good idea because some operations leaked out the primality condition block, but just for the pleasure to quizz colleagues:

quizz: aBlock
"Could you tell which serie will feed aBlock when this message is sent to an integer n?

    What the return value means?"
    | s t x y z |
    s := (x := y := 1) << self - 2.
    [(x := x + 1) <= self]
        whileTrue:
            [y := (z := y) << 1 + 1.
            (s bitAnd: (t := z + 1)) = 0
                ifFalse:
                    [aBlock value: x.
                    s := s bitAnd: z << (self // x * x + x) // y + t]].
    ^s

Did you notice the variant?