Sunday, 1 January 2017

BigIntegers and BigDecimals now on GitHub

To make it easier for me, as well as for downloaders, I put my BigIntegers and BigDecimals code on GitHub.

You can find these and related files (tests, debug visualizers, test data generators, etc.) on my BigNumbers project on GitHub.

Monday, 18 July 2016

New: Velthuis.AutoConsole unit

I often (really, very often) write console programs. Most of them are simple test programs to test algorithms, calculations or language features, or stuff copied from StackOverflow.

But console programs, if run from the IDE, or from the Windows Explorer, have one annoying feature: they close immediately if they are finished. So if my program has many Writelns in it, they simply rush by and the console window closes before I can inspect the output. To avoid that, I often put Readln; in the last line of the program. If the program is finished, it pauses and waits for me to press Enter.

But such a Readln; is annoying if I really want to run the program from an existing console window (aka a Windows command line). Then I don't want to have to press Enter, because the output is shown in the current console window already, and that won't close until I tell it to.

To get the best of both worlds, I wrote my simple Velthuis.AutoConsole unit. Put it in the uses clause of your console program and, if the program is a console program and it was started from the command line, or if it is not a console program, it will do nothing. But if it is a console program and it was started any other way (from the Windows Explorer, from the Delphi IDE with or without debugger, from another non-console program), then, before the console window can close, it will display:

Press any key...

Then I really only have to press a key (it does not have to be Enter, it can be any key) and the program will end, closing the console Window.

I placed the new unit in the zip for my Console unit. You can download it from the page I linked to.

Here is the code of this simple unit. I think I may have to add a ConsoleCtrlHandler, so the window is not closed when the user presses Ctrl+C or closes the console window with the  x  close button on the window frame. But I am not sure if that is necessary. If I press those, I usually want to stop the program anyway. Comments on this are welcome.

unit Velthuis.AutoConsole;



  Velthuis.Console,    // as described in the link above 

function GetConsoleWindow: HWnd; stdcall; 
  external 'kernel32.dll' name 'GetConsoleWindow';

function StartedFromConsole: Boolean;
  ConsoleHWnd: THandle;
  ProcessId: DWORD;
  ConsoleHwnd := GetConsoleWindow;
  if ConsoleHWnd <> 0 then
    GetWindowThreadProcessId(ConsoleHWnd, ProcessId);
    Result := GetCurrentProcessId <> ProcessId;
    Result := False;

procedure Pause;
  if IsConsole and not StartedFromConsole then
    Write('Press any key... ');





I removed the reliance on the Velthuis.Console unit. I added this piece, a simplified version of ReadKey, which does not try to translate the key, making things a lot simpler:

procedure WaitForInput;
  InputRec: TInputRecord;
  NumRead: Cardinal;
  OldKeyMode: DWORD;
  StdIn: THandle;
  StdIn := GetStdHandle(STD_INPUT_HANDLE);
  GetConsoleMode(StdIn, OldKeyMode);
  SetConsoleMode(StdIn, 0);
    ReadConsoleInput(StdIn, InputRec, 1, NumRead);
  until (InputRec.EventType and KEY_EVENT <> 0) and InputRec.Event.KeyEvent.bKeyDown;
  SetConsoleMode(StdIn, OldKeyMode);

and changed the Pause procedure a little:

procedure Pause;
  if IsConsole and not StartedFromConsole then
    Write('Press any key... ');

Friday, 20 May 2016

New blogging site

The blog that I had on the TeamB site is not functional anymore. I guess that is because it was still hosted by Borland, as far as I know. My blog pages were moved to the new Embarcadero Community site, but all comments were lost (sorry folks, but I had nothing to do with this), my user name is listed as "Rudy Velthuis Velthuis (TeamB)" and there seems to be nothing I can do to change this. And a search for "Rudy" or "Velthuis" on the Bloggers page doesn't find it either, or well, it says it found one instance, but it doesn't show it.

But what is worse, I can't write, edit or administer the posts there. That is why I decided to reactivate this blog site, which I had reserved ages ago. I have complete control over what I post here. I managed to copy most posts from the Community site over, and updated the formatting where that was necessary.

So please, if you had any links to my old blog, please update them to point to this site. The old URL will redirect to the Community site, but I will probably not be using that anymore. I hope that it will be picked up by DelphiFeeds and similar sites soon, because I intend to write regular blog posts again.


I was told that the glitches and problems with the moved blogs on the Community server are going to be solved, and that the site will be made faster as well. They'll even try to re-install the comments. If that is the case, I will move back to the Community site again. Hang on.

Sunday, 10 April 2016

New: BigDecimals

After BigIntegers, the next logical step was the implementation of BigDecimals. I implemented them using BigIntegers, because a BigDecimal is not much more than a BigInteger divided (or multiplied) by - scaled by - a power of 10.

Unlike BigIntegers, my implementation of BigDecimals does not use any assembler, so it should be usable anywhere BigIntegers can be used.

The BigDecimals are mostly modelled after the interface of the BigDecimal class in Java, but, like BigIntegers, made more Delphi-like. Java BigDecimals do not support operator overloading, while my BigDecimals do. This makes the handling of rounding and precision a little different too.

Anyway, I think that every Delphi should not only have BigIntegers, but BigDecimals too. Have fun with them. They can be found in the same .zip file as BigIntegers. Just read more about them on my website.

The next step is BigDecimal math, like trigonometric, exponential functions etc. One other step will be to host this on GitHub. I am also thinking of writing prime, factorization and other mathematically useful functions for BigIntegers.

Monday, 1 February 2016

In Delphi, is addition commutative?

Answer: no, not necessarily. Despite the existence of operator precedence, i.e. the fact that the following

  X := 3 + 4 * 5;

results in 23 and not in 35, the order of operands can still have an effect.

In my BigInteger code, I discovered an odd error, that only happened in some very rare cases, and only in PUREPASCAL code, i.e. code that did not use assembler, and only in 64 bit.

It took me several hours to find out that this was the problematic expression:

  Value := NormDividend[I + J] + NormDivisor[I] + Value shr CDivLimbBits;

CDivLimbBits is a constant with value 32, Value is an Int64, NormDividend[] and NormDivisor[] are arrays of UInt32 (Cardinal). Only in some very special circumstances, this caused an error.

What happened?

In this unit, which does lots of odd and ugly things to UInt32s and to speed things up, I turned off range and overflow checks, so it went unnoticed that NormDividend[I + J] + NormDivisor[I] caused an overflow. Since overflow and range checks were off, the 33rd bit simply got cut off.

But you might say: "Hey, the third operand is an Int64, so why were these two operands not promoted to 64 bit?" It turns out that this only happens once it is required, so what the compiler actually compiles is:

  UInt32(Intermediate) := UInt32(NormDividend[I + J]) + UInt32(NormDivisor[I]);
  Value := Int64(Intermediate) + Value shr 32;

while I expected:

  Value := Int64(NormDividend[I + J]) + Int64(NormDivisor[I]) + Value shr 32;

Now, if you rearrange the Int64 to come first, like:

  Value := Value shr 32 + NormDividend[I + J] + NormDivisor[I];

then all is well. The first operand is an Int64, so all following operands are promoted too and you really get:

  Value := Value shr 32 + Int64(NormDividend[I + J]) + Int64(NormDivisor[I]);

Note that this error did not happen in 32 bit code. There, NormDividend[] and NormDivisor[] are arrays of UInt16, and Value is an Int32. In other words, in 32 bit code (and even in 64 bit code on Windows), everything seems to be promoted to Int32 (signed 32 bit integer) anyway, probably because that type is somehow the preferred type for integer expressions (most integer code uses 32 bit registers, in 32 bit as well as in 64 bit).

So take care to either cast to the required type, or to put the largest operand first, otherwise you might be in for a surprise. It certainly was a surprise to me, especially because the data I had used for unit testing had not caught this.

Only the fact I wanted to improve the speed of ToString (converting the largest known prime to a string of 22 million decimal digits still takes approx. 2'30", while converting it to a string of — much more easily convertible — hex digits only takes 135 ms), and the coincidence that in one test I had to divide by exactly 10^57, made me catch this error. Note that the assembler code did not suffer from this. There I can control exactly what gets promoted and when.

This also made me aware again of the fact that testing can only show the presence of errors, and never the absence, and that it is extremely hard to find test cases that cover everything. The fact I had to divide by a number that caused the error was sheer coincidence.

Thursday, 19 November 2015

New: BigIntegers

I wrote a new implementation of BigIntegers (almost unlimited size integers) for Delphi. It uses assembler for the low level stuff, although there is a PUREPASCAL version for each of these routines too. Instead of using the GPL-licensed GNU Multi-Precision (GMP) library, I wrote everything from scratch, without using any code from others.

I think it is pretty fast, especially because I optimized many of the routines, either by techniques like loop unrolling or using 64 bits at once, or by using bitwise tricks, or by high-level algorithms like Burnikel-Ziegler, Karatsuba, Toom-Cook 3-way, etc. Results are checked with results generated by .NET. The C# test program is included.

In the course of writing this, I read an enormous numbers of publications on these subjects and I learned a lot, not only about the maths involved, but also about many other mathematical notions. Not so easy for a dentist. I also had some good help on StackOverflow.

Go take a look at my website. I think every Delphi should have BigIntegers.

Sunday, 4 October 2015

Speed problems caused by code that never ran

Let me start out by telling you that I am currently implementing a (simple) BigInteger implementation for Delphi. This implementation tries to be compatible with the BigIntegers as used in .NET, so I can test my results to be the same as the results of my tests in C#.

For this, RightShift must use two's complement semantics, and shifting a negative integer must shift in 1-bits from the top, so -100 ($FFFFFF9C) shifts to -50 (hex $FFFFFFCE). For normal Integer values, the result in Delphi would be 2147483598 (hex $7FFFFFCE), because Delphi does a simple shift and does not preserve the sign bit. But my BigIntegers should, because most other Biginteger implementations (including the .NET one, my reference) do that too.

To achieve this, in my routine, I did something like this:

class operator BigInteger.RightShift(const Value: BigInteger; 
  Shift: Integer): BigInteger;
  LSize: Integer;
  ShiftOffset: Integer;
  RSize: Integer;
  P: PLimb;
  if Value.IsZero then

  if Value.IsPositive then
    // Call internal methods that allocate the result, shift the value right, etc.
    // Works fine, no problem.
    // Recursive call
    Result := MinusOne - ((MinusOne - Value) shr Shift);

That gave me the right results, so I was happy until I started benchmarking the new routine. I noticed that the right shift was more than twice as slow as the corresponding left shift. I did not understand this, because a right shift actually has to move fewer bits and the result is smaller. Even if I only passed in positive values, it would still be slow. The code was almost the same as the code for left shift, though. Well, except for the else clause.

People cleverer than I am probably already see the problem: the "negative" branch (the else clause) of the if clause. When I removed it, code was indeed faster than that for left shifts. But when I put it back in, even though it would never run, it slowed down everything. Only after some hard thinking and some debugging I noticed, in the CPU view, that there were calls to InitializeArray and FinalizeArray in the compiled routine, and that everything was surrounded by a try-finally block. The expression in the "negative" branch had some intermediate results, and records for these were allocated and managed by the runtime. That made the entire code slow.

First solution

My first solution was to put the code for the "negative" branch into a nested procedure of its own. The hidden local BigIntegers were now confined in that procedure and would not cause the entire routine to be slow. The "positive" part was indeed up to speed now.

I finally deconstructed Result := MinusOne - ((MinusOne - Value) shr Shift); into internal calls entirely, so no hidden BigIntegers were allocated anymore. Now I could put the "negative" code back from the nested procedure to the negative branch, and it was quite a bit faster as well.


I learned two things from this:

  • Beware of hidden code. Here, it was caused by the expression with intermediate results.
  • Inline with caution. If I had inlined the nested routine, the hidden BigIntegers would be put back in the outer scope, and things would have been slow again.