Accelerating Large Integer arithmetic in Squeak Smalltalk

I recently introduced some changes in Kernel package of Squeak development trunk http://source.squeak.org/trunk.html so as to accelerate Large Integer arithmetic involving huge numbers.

The operations boosted by this set of changes are:

• multiplications *
• divisions / // \\ quo: and rem:
• squaring squared
• square root sqrt sqrtFloor

All these changes can easily be made optionals. There is only one mandatory change in the Kernel package, the introduction of an indirection multiplyByInteger: and divideByInteger: in Integer so that subclass can define their own hook. Default behaviour of the hook is to call the schoolbook primitives primDigitMultiplyNegative via digitMultiply:neg: and primDigitDivNegative via digitDiv:neg:. I believe that Pharo and Cuis Smalltalk should adopt these basic changes, as the penalty is really neglectable, and make the accelerated arithmetic a package on its own. I've also inquired about Visualworks, but the changes required to the base classes are a bit more invasive.

For accelerating those operations, I have used these classical divide and conquer algorithms:

• Karatsuba multiplication in digitMul22:
• Toom-Cook multiplication in digitMul23: and digitMul33:
• the Fast Recursive Division - MPI-I-98-1-22 from Christoph Burnikel & Joachim Ziegler in digitDiv21: digitDiv32: and driver digitDivSplit:
•  and beautiful Karatsuba Square Root from Paul-Zimmerman - Rapport de recherche n° 3805
  
  lk*jFiGhFg6=f42/eHeW>4@/eH+4YW􏰓4dc`IRWGCba2aFJ`
  

Here, the 22, 23, 33, 21, and 32 in selector names tells in how many parts the operands (receiver and argument) are divided, 22 means two parts each (Karatsuba).

The multiplyByInteger: method is refined in LargePositiveInteger and serves as a dispatcher to the available variants, based on heuristics on the operands digitLengh. Remind that those digits are 1-byte long (8 bits) in traditional Smalltalk, though the LargeIntegersPlugin of opensmalltalk-vm uses 32-bits digits under the hood - that is why the split is made at four bytes boundary (notice the bitClear: 2r11 used in most digit-splitting methods).

When the operands sizes are unbalanced (the larger is twice longer than the smaller or more), a digitMulSplit: is used with a strategy of splitting the larger in chunks of about 1.5 times the smaller. I have tried other combinations, 1-1 and 1-2, and measured that 1-1.5 was faster, thus the addition of a digitMul23: variant of Toom-Cook. It's interesting to notice how I reconstruct the result from collected parts in O(N) by in-place modification. Interlacing is used so as to be sure to get non-overlapping digits (see inplaceAddNonOverlapping:digitShiftBy:). My first iteration did use a O(N log N) of the form (a3 * x + a2)*x^2 + (a1 * x + a0), where multiplication by x correspond to bitShift: of the appropriate digit length, but each operation means allocating memory, moving bytes, and ends up in non neglectable cost noticeable in profiling...

The 3-way Toom-Cook implemented in digitMul33: is a variant from Marco Bodrato & Alberto Zanoni, and it is interesting to read What About Toom-Cook Matrices Optimality? on this subject.

Squaring also use Karatsuba and Toom-Cook, and these are asymetrical versions which outperform the symetrical ones. Specifically, squaredByFourth use the 4-way Toom-Cook variant of  Jaewook Chung and M. Anwar Hasan - Asymmetric Squaring Formulae   https://www.lirmm.fr/arith18/papers/Chung-Squaring.pdf.

The new square root algorithms introduces a new message sqrtRem which answers an Array with truncated square root and remainder for almost the same price
self sqrtRem first squared + self sqrtRem last = self
 
This is just a starting activity, there is more to do: I have played with exponentation (raisedToInteger:), cubing (cubed) and cube root (cbrtRem as a generalization of sqrtRem) and maybe I'll blog about them, or publish if worth. There is also a need to auto-adapt the thresholds hardcoded here and there, but this is quite involved to get it right, benchmarking is an art.

Why bother with those classical algorithms without much innovating, when the state of the art GNU Multi Precision library (GMP) is far more advanced, optimized, fast and ready to use? Because I have an interest in supporting ArbitraryPrecisionFloat, and more than all, because it's very easy and fun to experiment in Squeak. It's fascinating how straight forward it is to transcribe those algorithms in Smalltalk and I believe that researchers should prototype new algorithms in such a language before paying the ultimate tribute to optimization. If you find that extending GMP is equally fun, which means programming in C, dealing with memory optimization, edit-compile-run-segfault cycles, then OK, you're even more a weirdo than I am!

 

Swimming against the tide

This is an old article of 2018 that I forgot to publish, about how argh-woops-dang the programming experience can be sometimes, maybe the title should have been swimming in bitumen.

HDF5 is an essential format for exchanging large chunks of data. I don't much like it. It's overly complex. But not having an interface to it is a severe handicap for a language aiming to be generalist, and crunching data is rather a generalist task these days.

Of course, I could be using Python like anyone else on earth, but I would prefer some more productive and lively environment. So a few months ago, I've started to write such an interface for exchanging data between HDF5 and Smalltalk (Visualworks for a beginning). I'll publish on public store when having minimal usable core features, and it will be licensed MIT. But it's not yet ready. It's an activity taken on my free time.

I've started development in windows because having a not-supported-by-apple library working on Macosx is rarely an "out of the box" experience. This week end, I decided to give it a kick and to switch back to my preferred platform.

Building the library with the cmake build configuration provided by the hdfgroup was fairly easy. Just follow the instructions from http://support.hdfgroup.org/HDF5/release/cmakebuild.html.... and produce what you don't need: static libraries.

Building the dynamic libraries has been more involved, because neither setting the STATIC_ONLY=NO nor NO_MAC_FORTRAN=YES in build-unix.sh did the trick, whatever the recommendations found in http://support.hdfgroup.org/HDF5/release/chgcmkbuild.html.
The worse is that I don't know why. Did I forget to remove the previous build?
It's certainly not the right way, but brute force hardcoding of set(ADD_BUILD_OPTIONS "${ADD_BUILD_OPTIONS} -DBUILD_SHARED_LIBS:BOOL=ON") in the various cmake files produced the libhdf5.dylib that I wanted...

Or not... These were x86_64 libraries.  With my dated Smalltalk environment, I wanted an i386, or better a fat 96bits universal version. Not a problem. As usual, Stackoverflow has THE solution: http://stackoverflow.com/questions/5334095/cmake-multiarchitecture-compilation. I stupidely passed option "-DCMAKE_OSX_ARCHITECTURES=x86_64;i386" to ctest, which was not the right incantation, and without any flashing illumination, but rather after laborious attempts, finally added it as build option inside CTestScript.cmake

set (BUILD_OPTIONS "${BUILD_OPTIONS} \"-DCMAKE_OSX_ARCHITECTURES=x86_64;i386\" "). 

This was the right syntax, but not the right way...
My move was anticipated for a long time by another chess player, see those lines in
 hdf5-1.10.1/config/cmake_ext_mod/ConfigureChecks.cmake:

if (APPLE)
  list (LENGTH CMAKE_OSX_ARCHITECTURES ARCH_LENGTH)
  if (ARCH_LENGTH GREATER 1)
    set (CMAKE_OSX_ARCHITECTURES "" CACHE STRING "" FORCE)
    message(FATAL_ERROR "Building Universal Binaries on OS X is NOT supported by the HDF5 project. This is"
    "due to technical reasons. The best approach would be build each architecture in separate directories"
    "and use the 'lipo' tool to combine them into a single executable or library. The 'CMAKE_OSX_ARCHITECTURES'"
    "variable has been set to a blank value which will build the default architecture for this system.")
  endif ()
  set (${HDF_PREFIX}_AC_APPLE_UNIVERSAL_BUILD 0)
endif ()

Ouch! I don't care much of the best approach. I just want something that works! 
And I don't even agree. The best approach is to follow the guidelines of the target platform. 
Everything else is swimming against the tide...
For the time being, it's me who is swimming against the tide, too far away from the mainstream.
Even if I don't build for universal but just i386 architecture, the option will be blanked.
So what are these technical reasons exactly? Where can I read about that?

Without a clue, I started to hack those lines if (ARCH_LENGTH GREATER 1000), and also  set (${HDF_PREFIX}_AC_APPLE_UNIVERSAL_BUILD 1) but of course, as advertised this did not work. OK, with -VV verbose verbose option, we know that the link failed because ZLib and SZip were compiled for x86_64 only too. Untar-ing the provided ZLib.tar.gz and SZip.tar.gz, modifying the ConfigureChecks.cmake inside which had the very same guard as above, re-taring-re-gzipping ... equally failed! Same problem, the external libraries were NOT compiled with universal support. Adding the DCMAKE_OSX_ARCHITECTURES=x86_64;i386 to various CMAKE_ARGS definitions that I found did not help either...

Hacking randomly like this is vain and probably doomed to fail, but how are you supposed to understand what happens and where it comes from exactly among those thousand lines scattered all around the .cmake files?

You think I'm exagerating? 

find . -name '*cmake*' | wc -l
-> 101

wc -l `find . -name '*cmake*'`
-> 23729

More than twenty thousand lines in a hundred files! With all that stuff, HDF5 should certainly compile on your android tablet and your connected fridge! But not on my Mac...

If you count the generated build directory, this roughly doubles.
But if you remove the duplicates (and blanks) it's only half:

cat `find . -name '*cmake*'` | sort |  uniq | wc -l
-> 12799

Believe me now?

Back to the technical reasons they spoke about, I now have a clue: could it be the choice of the right tools? Ugly scripts, ugly environment variables, ugly macros, duplicated code and no debugger for getting a chance to dive inside the machinery. Welcome in the 21st century.

That's what sucks with configure/cmake and all those meta levels that generates configuration files for yet another tool (make): you gradually loose control on the lower levels where things has to be solved. The problem is fairly simple: compile/link ZLib and SZip with the universal build flags. The solution of going across all these layers is really involved, even if I don't even want to make it work for others, but just for me!

I could have talked about Smalltalk and the implementation, but half the week-end is gone in smoke now without the slightest progress.

Slow HDF5 progress

I've been struggling with DLLCC, and found another bizarre behavior that I reported on S.O., http://stackoverflow.com/questions/49564024/isnt-pointer-type-checking-disabled-in-dll-c-connect-and-is-that-ok. I then took some time to reread the F... manual.

I think that I have deeper understanding of DLLCC now and published a new patch using double dispatching on the public store. This enables passing a Smallapack.CArrayAccessor as argument to the external function call. What's this thing (the name is not so nice)? It's a proxy to some C data (on heap), somehow like a CPointer (a CType plus a pointer on C data), but with two differences: 

1. it carries a length, and thus enables safe access of contents from within the Smalltalk world. If you think of it two minutes, every pointer should be bounded, that's clearly the way to go. Alas, in 2018 we still are at assembler level with those C API, so the pointer coming from external world still are unbounded (a pity for safety...).
2. it is using 1-based index and enables handling from within Smalltalk world like any other collection without too much mindstorm.

With the time spent on DLLCC, the progress was slow on HDF5, I just corrected the H5TString transfer. The documentation was not clear about fixed size null terminated strings. I first thought (or read?) that the terminating null was mandatory. But then the type H5T_C_S1 - a C type string of fixed size 1 - makes no sense! So I speculated that an extra byte was allocated for the terminator... This was wrong! If not enough place is allocated, then the null terminator is omitted. The difference of NullTerminated and NullPadded HDF5 Strings is then germane: the former ends before first null, the later at last non null, but in neither case null is mandatory...

Along with the CArrayAccessor change, I can now explore the contents of a HDF5 file created by Matlab save -v7.3. a Small step in the right direction.

HDF5 interface requires a patch in DLLCC

Apologies to those who would have tried the HDF5 interface (if any): I completely forgot that I had to patch argument coercion in DLL/C-Connect to let it work.

I suspect that this is a bug of DLLCC, see http://stackoverflow.com/questions/49544642/why-cant-i-pass-an-uninterpretedbytes-to-a-void-thru-dll-c-connect

I have published a new bundle with the override. Alternate solution would be to either pretend that the function can take an _oopref* buffer argument, or to go thru copyToHeap copyFromHeap: complications that I absolutely do not want (HDF5 data can be big, thus I do not wish any extra/unecessary copy).

Interfacing Smalltalk with HDF5

I am going to speak about connecting Smalltalk with external (scientific) data again.

This time, it's not with Matlab, but HDF5, the hierarchical data format supported by http://www.hdfgroup.org 

Every language targetting science/engineering niche must have an interface to HDF5. That is the case of Python and Matlab to only cite two. And you know that I'd like to promote the usage of Smalltalk in this area too. So lets do it: here comes the HDF5 bundle in Cincom public store.

If interfacing with Matlab mat-file format was a piece of cake, HDF5 is much more involved. First, because HDF5 is like a file system. There are recursive named Group of data, like directories. And groups are not necessarily arranged as a tree, but can form arbitrary graphs (circular) thanks to links. Second, because HDF5 comes with a type system. It can hold arbitrary types, whether atomic (integer, floating point,bit fields) or composite (structure, arrays). The types can even be references to other objects, which means that it is sufficiently general to describe heterogeneous collections of dynamically typed data.

There are two kind of way to store data: the first is Dataset which are named entries in Groups (a bit like files in directories). a Dataset is a rectilinear multiple dimension array (like the MultipleDimensionArray that I recently promoted in Visualworks public store and Squeak STEM). The second is Attributes. All named entries (Group, Dataset and named types) can have attributes, which are kind of string-key arbitrary-value pairs. Attributes generally hold meta data. They lack the ability to perform the read/write of sub-regions of data. We somehow can compare attributes to the property list of Morphic.

For additional complexity from the interfacing point of view, arrays can be of variable length. Thus the buffer for holding complex data cannot be preallocated from within Smalltalk before the data transfer, but has to be allocated on the fly by the HDF5 library. Thats potentially means either memory leaks or dangling references to freed memory, or the two if the programmer worked too late at night.

HDF5 comes with lot of documentation including reference, user guide, tutorial and examples if you want to learn more.

Among the many features of HDF5, let's focus of some of the most essential:

• the ability to perform type transformations during transfer operations. This enables language interoperability, since the type system is flexible enough to accomodate many languages (if not all thru C interface).
• the ability to read/write sub-regions of dataset. This enables handling of huge data, the whole dataset does not have to reside in memory.
• it scales well (or at least it can if used adequately).

That being said, this comes with a price: HDF5 is complex, and like the manipulated data, the API is bigger than big too.

It's clear that the target languages that HDF5 creators had in mind were statically typed. So the mapping to dynamically types objects is not straightforward, nor optimized, especially for the compound type (structure). The transfer necessarily involves two steps: the transfer or raw data, followed by a transformation to Smalltalk data for read, et vice et versa for write. For handling all cases, including references, a visitor pattern will be necessary (references can be cyclic too).

For huge data, the way to go is to use proxy to HDF5 objects. That means giving minimal behavior to the proxys, at least for storing/retrieving whole data or subregions. Application specific behavior should better lie in application specific objects, and those objects would use a generic HDF5 proxy. Since data structure is hierarchical (recursive), an application specific visitor should construct the application specific object graph. The early implementation that I just published is not yet there. It does not even use a visitor, but hardcode an arbitrary visit in the HDF5 proxies. It's more a minimal proof of concept at this stage, but a usable proof of concept though.

So, how do we use it? For creating a HDF5 file, first do this:

out:= H5File create: 'foo.h5'.
H5Dataset createPath: 'float1' parent: out rootGroup value: 1.3e0.
H5Dataset createPath: 'double2' parent: out rootGroup value: Double pi.
H5Dataset createPath: 'int3' parent: out rootGroup value: -357.
out close.

The close operation is not strictly necessary, because I implemented a registry of opened hdf5 objects with auto-close facility when the entries are reclaimed. But forcing a close fushes the file, otherwise the close will be delayed until all opened HDF5 entities have been reclaimed by garbage collector.

For reading, it's like this:

in:= H5File readOnly: 'foo.h5'.
float1 := (in / 'float1') value.
double2 := (in / 'double2') value.
int3 := (in / 'int3') value.
in close.

This example is a bit poor. We can also create/query groups, attributes and handle more complex data, thanks to MultipleDimensionArray, RawArray for atomic data types, the CArrayAccessor of Smallapack for arbitrary data. 

For writing Smalltalk objects on HDF5 files, one must define these 3 essential selectors:

• h5mem returns a buffer containing the data to transfer and usable by HDF5 API (must be compatible with DLL-C-Connect interface)
• h5type gives a HDF5 type description of buffer contents
• h5space gives  a HDF5 description of buffer layout (dimensions of the dataset)

Last note: I've used HDF5 version 1.8.10 for the interface. Unfortunately, HDF5 use macros for enabling backward compatibility, but the way DLL/C-connect works, we will have to subclass the H5Interface in order to support various versions.

That's all for this post, I may add more implementation details later on.

MatFileReader is ported to Squeak

I've ported the MatFileReader package initially written in Visualworks to Squeak. It should work unchanged in Pharo but I did not test. Usage in Squeak is simply:

    (MatFileReader on: (StandardFileStream readOnlyFileNamed: 'test.mat') binary) decode.

This will answer a Dictionary with workspace variable names as keys and MxArray as values which are essentially stubs by now. With a visitor pattern or by subclassing MatFileReader, it is possible to map the values to other Smalltalk classes.

Code repository is on squeaksource, http://www.squeaksource.com/STEM.html

    MCHttpRepository
        location: 'http://www.squeaksource.com/STEM'
        user: ''
        password: ''.

STEM stands for Smalltalk Tools for Engineering and Mathematics (or Squeak Tools if you're tainted). The name is short, close to the main acronym for Science Technology Engineering Mathematics, and otherwise have quite many meanings (http://en.wikipedia.org/wiki/Stem) . I think I will use it by now.

I've not assembled any MonticelloConfigurationMap nor prepared any Metacello ConfigurationOfSTEM, so by now you'll have to install all the packages in this repository by yourself.

Note that STEM is not a concurrent of Polymath (http://github.com/PolyMathOrg/PolyMath). Polymath is for Pharo, and I don't yet use Pharo, I'm far more fluent in Squeak. License is MIT, and volunteers for porting/integrating into Polymath are welcome if the library sound interesting enough. I just ask to have a little respect for authorship, if ever the environment still permits it.

Importing Mat-file in Smalltalk

A very long time since I did not write anything here...

I've been busy, but finally found time to release a small utility in Cincom Smalltalk public store: it's a package named MatFileReader.

The package is capable of importing data from the .mat files produced by Matlab, into Visualworks. Nothing spectacular, but it may be useful for scientific applications. Usage is fairly simple, just do-it:

     (MatFileReader on: 'foobar.mat' asFilename readStream) decode.

I had developped the package around 2005, but it was messy and full of personnal dependencies. The dependencis are now reduced to the BinaryStream,  RawArray and MultipleDimensionArray packages.

• RawArray are collections holding bits of data efficiently in space - bits are translated into/from Smalltalk objects with at: and at:put: operations like a ByteArray or a WordArray.
• BinaryStream is just a Stream wrapper with convenient messages for reading chunks of binary data, especially into RawArray.
• MultipleDimensionArray is just a SequenceableCollection wrapper which let the single dimension raw storage be viewed with multiple dimensions. The layout is that of FORTRAN and Matlab, column-wise (first index varies first).

Also two other dependencies are the unzip from PNG readers and EphemeralValueDictionary from Glorp which have been used for reading compressed data and automatically reclaiming classes respectively (classes are created on the fly for holding Matlab struct data). When Ephemerons will be ready in Squeak, this package can provide a good test.

The package would be more useful if data were imported into some application specific classes (like Smallapack for example), but for that one just have to either subclass MatFileReader, or post-process with a tree visitor.

License is M.I.T.

When I'll have time, I'll also handle the export to mat-file, and try to port into Squeak/Pharo. Such package could be useful for Polymath (https://github.com/PolyMathOrg/PolyMath).