Extend your 50g with C (HPGCC 2.0SP2)

This lengthy article explains why you would and how you can extend the functionality of your 50g using C. Complete examples are provided to illustrate how to create high performance mathematical routines such as a complex LogGamma function, a sparse linear solver, and a 2D convex hull.

There are two reoccurring themes in this article:

Keep it simple and use the right tool for the right job. C is not a complete replacement for UserRPL. E.g. C it can require a lot of code to do what UserRPL does in one or two commands when processing stack arguments. One the other hand C is faster than Saturn assembly for intensive computation. C is also a lot easier to read, write, and debug with larger applications.
Do not reinvent the wheel. There is a lot of freely available mathematical C code available. Most of the examples will focus on using freely available C code.

This article is not a replacement for HPGCC, C, or 50g documentation. This article is structured more like a tutorial. In-depth knowledge of C is not required. It is possible for a person without any C experience to compile and run all the examples. However, you will get more out of this article if you know C. If you want to learn C then I suggest one or both of the follow texts:

The C Programming Language, 2nd ed., by Kernighan and Ritchie. This is the C book and is frequently referred to as K&R. This book is short and sweet (about 1cm thick). If you grasp the obvious then this is the book for you. I learned C from the first edition. The first edition is a collectors item now, I wish I had kept mine when I got the 2nd edition.
C Primer Plus, 5th ed., by Stephen Prata. Chatty and verbose, many find this bible sized text easier to understand than K&R. More demanding topics like pointers and structures are well documented in this text.

Why C?

C is fast, very fast.
There is an abundance of freely available quality mathematical C code and libraries. This will save you a lot of time.
Rapid development and prototyping. Every example here was written for my workstation first then ported to the 50g in seconds.
Skill reuse. C skills (unlike SysRPL and Saturn Assembly) have use almost anywhere. 50g C programming is an example of embedded programming. Embedded cross-compiling, object size optimization, and debugging skills have many applications elsewhere.
Portability. With little effort you can rebuild your 50g C applications for PCs, Macs, PDAs, and even TI calculators (TIGCC).
C is cool.

Why not C?

Serial and USB I/O is currently not supported.
Small UserRPL routines may be good enough. Why add a layer of complexity? Use the right tool for the right job.

HPGCC

HPGCC is a cross-compiler suite available for x86 Linux, Mac, and Windows capable of creating 50g ARM binaries. With some effort HPGCC can be compiled for almost any platform. HPGCC is based on GCC (GNU Compiler Collection). The GNU Compiler Collection includes front ends for C, C++, Objective-C, Fortran, Java, and Ada, as well as libraries for these languages (libstdc++, libgcj,...) [1].

Officially HPGCC 2.0SP2 only supports C, however Jean-Yves Avenard created Mac and Linux versions that also support C++ [2]. With little effort other languages like Fortran could also be added. Image compiling your existing Fortran code for your 50g. It would probably run faster than the computer it was originally coded for!

Unlike UserRPL, SysRPL, and Saturn Assembly, 50g ARM binaries run outside of the 50g UserRPL environment. This has a few advantages, e.g. speed and access to more memory. HPGCC creates a contiguous memory block of the unused port 0 and 1 memory as usable application memory, plus there is a bonus 90KB of RAM unused by the 50g OS. That's a total of 459KB (assuming the Ports 0 and 1 are empty). Although 50g ARM binaries run outside of UserRPL it is still possible to interact with the stack.

Disclaimer

Use at your own risk. Although nobody has reported bricking their 50g with HPGCC, I suppose anything is possible. I imagine that the worse that could happen is that you break the reset button because of the countless times you had to reset it. I've only had one scare. A crash put lines in the overscan area of the LCD that persisted after a hard reset, after 10 or so seconds it faded away. I wasn't too scared, because I did it over and over again to see what would happen.

I'm not an expert at C, the 50g, or HPGCC. I have tried to be as accurate and thorough as possible. Omissions and mistakes are most likely present. Again, use at your own risk.

Enjoy.

Development Environment

Hardware Requirements

A 50g. Believe it or not you can do a lot of development and testing without one, but eventually you are going to have to test it on the real thing. NOTE: EMU48 with 50g emulation only emulates the Saturn processor emulation of the 50g, i.e. EMU48 cannot execute ARM binaries.
A 50g USB or Serial cable. I find the 50g USB cable to be the better solution since it also provides power to the 50g. Some have had problems with USB file transfers and opt for the serial cable. The best serial cable for the 50g can be obtained here: http://commerce.hpcalc.org/serialcable.php.
An SD Card reader. I prefer the SanDisk MicroMate because I can easily remove the SD Card from the reader without having to snap it in to some type of shell. Any will do, I just happen to have this laying around. This one looks good: http://commerce.hpcalc.org/cardreader.php.
An SD Card. Any size will do with 10MB free space. I am using an old 128MB Lexar SD Card. 128MBs is like infinity.
A paper clip. You will never want to run a C program on your 50g without one handy. In a pinch you can just remove a single AAA battery for 1 second to get the same effect.
A workstation. Currently it is not possible to compile C code on the 50g. It must be cross compiled with HPGCC.

Software Requirements

A workstation OS. This article was written for Windows XP, but should be easy to adapt to OS/X or Linux. Although I am a UNIX and Linux user, Windows is still the dominate desktop OS, so I went for LCD.
Cygwin/X. This is not an HPGCC requirement, just a requirement for this article. I find developing with GNU tools easier with a Linux/UNIX shell. This also makes this article a bit more universal. Everything should work the same for OS/X and Linux too. Cygwin/X is also required when testing code on the PC. NOTE: All the HPGCC source code will compile and link just fine using the native HPGCC Windows tools.
The hpgcc.tgz tarball. This tarball contains:
- HPGCC 2.0SP2: The 50g C cross compiler for Windows.
- HPAPINE: A 50g simulator for testing 50g C code. This is not an emulator. More on this later.
- All the example code, makefiles, and binaries from this article.
- This document.

User Requirements

You should know how to:

Transfer files from your PC to your 50g.
Transfer files from your PC to your SD Card.
Operate an SD Card in a 50g. If you have 50g SD Card problems this read this: http://h20331.www2.hp.com/Hpsub/downloads/50gUsing_an_SD_Card.pdf.

Cygwin/X Installation

Skip if Linux or OS/X.

Cygwin/X is a GNU/Linux like environment for Windows. It does not install any services or drivers, it is self-contained user space software just like any other application. IOW, it is safe. (Standard disclaimer applies.)

Cygwin/X adds the X Windows System to your Windows desktop. X Windows is 10 (make that 11) times better than just MS Windows. For more information about X Windows read: http://en.wikipedia.org/wiki/X_Window_System. BTW, both Linux and OS/X use X.

Download http://cygwin.com/setup.exe, virus scan it, execute it, install all packages, then take a nap. Cygwin/X download and install takes hours. Example session:

Run setup.exe:

Click Next.

Click Next.

Click Next.

Click Next.

Click Next.

Select any mirror then click Next.

Click on the symbol between "All" and "Default" to change "Default" to "Install".

Before clicking Next, verify that your dialog matches the snapshot above. I.e. All packages to be installed.

Click Next.

Take a nap, this could take hours.

Click Finish. Congratulations you now have a GNU/Linux/UNIX like development environment for Windows.
Setup two icons anywhere you like.
1. StartX: C:\cygwin\usr\X11R6\bin\startxwin.bat
2. Xterm: C:\cygwin\bin\run.exe -p /usr/X11R6/bin xterm -display 127.0.0.1:0.0 -ls
Test. First run StartX. This only needs to be run once. You'll notice an X in your tray. Right-click, exit if you need to exit X. StartX also starts up an Xterm. If you want multiple Xterms, then run Xterm. If all went well you should see something like this:

Congratulations you have a very powerful GNU/Linux environment at your finger tips. NOTE: Cygwin/X paths are a different than Windows paths. First, Cygwin/X uses the forward slash (/) as the directory separator, and 2nd, C: is replaced with /cygdrive/c. E.g. my Windows home directory is C:\home\egan, but from a Cygwin/X perspective it is /cygdrive/c/home/egan. BTW, the ~ (tilde) is often used as a short cut for the home directory, and will be used frequently in this article.

HPGCC w/ Extras Installation

Open up an Xterm and type:

cd ~ (puts you in the root of your home directory)

wget http://sense.net/~egan/hpgcc/hpgcc.tgz

tar zxvf hpgcc.tgz

NOTE: Linux and OS/X users will need to remove the ~/hpgcc/2.0SP2 directory and replace it with the proper HPGCC for Linux (http://hpgcc.org) or OS/X (http://www.hydrix.com/Download/Hp/hpgcc/).

Everything you need is self-contained in the ~/hpgcc directory.

Sidebar: Play First!

In the ~/hpgcc directory are two files that contains all the wrappers, binaries, and libraries. To install:

Drag the HPGCC.hp folder to your 50g HOME using the Connectivity Kit.
Extract EXTEND.ZIP to the root of your SD card.

Do not want to install anything on your workstation? No problem:

http://sense.net/~egan/hpgcc/EXTEND.ZIP

http://sense.net/~egan/hpgcc/HPGCC.hp

ARM Toolbox Installation

To run ARM binaries on the 50g you must have the ARM Toolbox installed OR have had the binary converted to a standalone executable using the ARM Toolbox. As a developer you will need the ARM Toolbox, but if you convert the binaries to a standalone executable then you can distribute your binaries without requiring that the ARM Toolbox be installed.

Copy the file ~/hpgcc/2.0SP2/sources/ARMToolbox/SETUP.BIN to your 50g.

TIP: Windows users type in your Xterm window:

cd ~/hpgcc/2.0SP2/sources/ARMToolbox/
explorer .

This will open an explorer window in the ARMToolbox directory so that you can drag and drop SETUP.BIN to the Connectivity Kit.
On the 50g put the number 2 on the stack:

Next, put SETUP.BIN on the stack, not 'SETUP.BIN'. Your 50g should look like this:

Now press , then +.
Verify that the ARM Toolbox is installed by pressing: .

Setup Environmental Variables

Each time you open a new Xterm type:

cd ~/hpgcc

. hpgccenv (yes, that is . SPACE hpgccenv)

This will setup your HPGCC and PATH environment for that Window only.

Let the fun begin!

Getting Started with HPGCC

All the C source for this section is in ~/hpgcc/test.

Your 50g working directory should be home/hpgcc/test.

Hello, World

Lets start with the classic hello, world from K&R:

#include <hpgcc49.h>
       
int main(void)
{
    clear_screen();
    printf("hello, world\n");
    WAIT_CANCEL;
    return(0);
}

There are a few differences from the classic hello, world:

The header file hpgcc49.h replaces the ubiquitous stdio.h. hpgcc49.h is a collection of other header files specific to HPGCC. Just about every HPGCC specific call is defined in the hpgcc49.h collection.
clear_screen(), while not entirely necessary is generally a good idea if you plan to output text directly to the screen. E.g. hello, world with and without clear_screen():
WAIT_CANCEL, is not necessary either, but, if you want to see your output before returning to the 50g UserRPL environment you will need something that pauses for input. WAIT_CANCEL waits for the key to be pressed.

To build and run this example, type the following commands:

cd ~/hpgcc/test
make clean
make hello.hp

Next create a 50g home/hpgcc/test directory and transfer hello.hp to that directory. To run, press (you should see a lot of garbage), then . To exit press .

It is possible to make this a standalone executable that does not require the installation of the ARM Toolbox. Press , then . Then press . This will overwrite the that required with a that can run without it. To run just press . Easy! From this point forward all examples will be assumed to be standalone executables.

Hi Bob!

Next lets try something a bit more interesting:

1       #include <hpgcc49.h>
       
2       int main(void)
3       {
4           char *name;
5           char buf[40];
       
6           name = sat_stack_pop_string_alloc();
7           sprintf(buf,"Hi %s!",name);
8           sat_stack_push_string(buf);
9           return(0);
10      }

This program will take a string from the stack and put it back on the stack with "Hi " and "!" wrapped around it. In detail:

Line 6 will pop the string off the stack, allocate some memory to hold the string, then return the address of the string in memory to the character string pointer name allocated on line 4.
Line 7 stores in the character string buf the contents of the string that name points to with "Hi " and "!" wrapped around it.
Line 8 pushes the the character string buf to the stack.

There is no need to clear the screen or wait for input. Compile this like the previous example and create a standalone executable, then put a string on the stack and press , e.g.:

Crashes

If you have not already take your paper clip and reset your 50g using the reset hole on the back of the 50g. Expect to do this frequently. You should only lose the contents of the stack.

Crashes are usually caused yours or someone else's bug. Expect most of them to be your bugs. Crashes are frequent with C mostly because of the misuse of pointers. In the real world they are called segmentation faults. To create one simply try to access a region of memory not assigned to one of your variables. Most CLI workstation programs will error with some type of useful message, e.g. "Segmentation fault (core dumped)". HPGCC segfaults return nothing. HPGCC stuck in a loop also returns nothing. Debugging HPGCC can be a challenge, but a welcome challenge given the power C brings to the 50g.

To illustrate this take a look at the hi.c program above, can you spot the bug? What if name points to a string longer than 36 characters? buf is only allocated for 40 characters with 4 being used for "Hi !". Try it. Put 40 character string on the stack and press . Make sure your paper clip is handy. My workstation version of this program returns a segfault, I suspect HPGCC does too, but is more subtle. Hopefully future versions of HPGCC will create a core or stack dump on the SD card for analysis.

Fixing hi.c can be done a few different ways. E.g., you could take the first 36 characters and discard the rest (35 really since each string must end in '\0' a.k.a. NULL), or you could dynamically allocate the size of buf based on the length of name + 5, or you could reject long strings.

Hangs

Hangs are not crashes. Hangs are probably stuck in a loop or waiting for I/O. In the CLI world you send an interrupt (usually a CTRL-C), but that is not an option with HPGCC. But with a little effort you can put emergency exits in your code.

E.g., what not to do:

#include <hpgcc49.h>

int main(void)
{
    while(1)
        ;
    //never gets here
    return(0);
}

The above code will never exit the while loop and you will be forced to reset.

#include <hpgcc49.h>

int main(void)
{
    while(1)
        if(keyb_isAnyKeyPressed())
            break;
    return(0);
}

This code will also loop forever but will exit with any key press. If you are going to have large iterative loops consider giving yourself an out. If you are developing code for others consider some sort of indicator other than the hour glass icon that something is happening so that they do not prematurely exit loops or hard reset. The Convex Hull example has an example of a graphical progress bar.

Wrappers

As stated in the introduction; C is not a complete replacement for UserRPL. UserRPL provides some very powerful tools for input and stack manipulation that may be more difficult under C. Take forms for example. They are not going to run any faster under C, so why bother? There is a correlation between the number of bugs in your code and the length of the code. IOW, keep it short and simple whenever possible and use the right tool for the right job. Often a simple wrapper around your C code can make a big difference, e.g. try using using hi.c without a string:

hi.c cannot pop a string off the stack. NULL is returned and then pushed out as "Hi <NULL>!". A simple wrapper like: \<< \->STR hi \>> solves this problem for almost any case, e.g.:

Send nothing and get a "Too Few Arguments" error, send anything else and it gets converted to a string first. Coding object type identification and conversion to string in C is possible but more difficult with little benefit.

Wrappers are also convenient for returning C errors, e.g. hi2.c:

#include <hpgcc49.h>

int main(void)
{
    char *name;
    char buf[40];

    name = sat_stack_pop_string_alloc();
    if(name == NULL) {
        sat_stack_push_string("hi2 Error: String expected!");
        sat_push_real(1);
    }
    else {
        if(strlen(name) > 35) {
            sat_stack_push_string(name);
            sat_stack_push_string("hi2 Error: String too long!");
            sat_push_real(1);
        }
        else {
            sprintf(buf,"Hi %s!",name);
            sat_stack_push_string(buf);
            sat_push_real(0);
        }
    }
    return(0);
}

In this improved version, hi2.c checks for NULL and long strings and pushes an error message and an error number (1 in this case). If no error, then push the new string and 0. Note that if the name string is too long it is pushed back on the stack to be consistent with other types of 50g errors, i.e. if there is something wrong with the argument leave it on the stack. This is not necessary with the NULL check since nothing was popped off the stack.

The new wrapper and errors look like this:

The new wrapper code IF statement pops off the 0 or 1 return code, if nonzero the string returned is used by DOERR to put a nice error message on the screen. Otherwise, if zero, then END leaving the string "Hi ... !" on the stack. This is a nicer solution than blindly converting to a string before calling hi because you do not get unexpected results.

Another benefit of wrappers is that it is simple to create libraries for your newly created collection of C routines. This will be documented in detail at the end of this section.

Data Types

C is a statically typed language, i.e. variables are declared by type in advance at compile time. UserRPL is dynamically typed, i.e. variables are declared and typed at runtime. This static/dynamic dichotomy can be a challenge when passing data to and from the stack. This was clearly illustrated above with hi.c. To further compound the issue UserRPL has many more data types than C. However, C can use structures to create more types, e.g. a complex number could be represented as two real numbers.

In this article UserRPL wrappers will be used to precondition data passed to C programs and to interpret the results. This is a simple solution to the problem, but not the only solution, e.g. the freely available HPStack [3] HPGCC library provides C structures and routines to handle UserRPL data types such as lists, arrays, and complex numbers.

A detailed explanation of C data types, arrays, and structures is beyond the scope of this article, but it is an important topic that is discussed at length in any good C book. From a data management perspective it is good to know how to correlate C and basic UserRPL types and how much RAM each type takes so that you can budget your memory usage. The following program will return the number of bytes used for the common C types:

#include <hpgcc49.h>

int main(void) {
        clear_screen();

        printf("size_t: %d\n",(int)sizeof(size_t));
        printf("  char: %d\n",(int)sizeof(char));
        printf(" short: %d\n",(int)sizeof(short));
        printf("   int: %d\n",(int)sizeof(int));
        printf("  long: %d\n",(int)sizeof(long));
        printf(" float: %d\n",(int)sizeof(float));
        printf("double: %d\n",(int)sizeof(double));

        WAIT_CANCEL;
        return(0);
}

size_t is used for pointers and is 4 bytes (32 bits) long, this is expected since the ARM processor has a 32 bit address space. char is a single 8 bit character, it can also double as a very short signed or unsigned integer. Note: C has no string type, an array of characters is used instead. short, int, and long are signed integers. float and double are signed reals.

UserRPL integers can be popped and pushed as any C integer type as long as it will fit. Because the 50g has arbitrary precision (big num) integer support not all integers can be popped as C integers. If your C code understands large integers (see Unbounded Spigot, Next Prime example below), then conversions to and from strings would be necessary.

UserRPL reals can be popped and pushed as C doubles. Liberal use of the ->NUM function in wrappers can save a lot of frustration. E.g. 1/2 in exact mode will not be popped as a real. ->NUM will convert any number to a real (even integers, symbols, and variables).

C Binaries Are Big, Too Big

If you been following along and compiled/installed all the examples up to this point, you may have noticed that your amount of free RAM has been significantly reduced, e.g.:

76 KB of RAM has been used for a handful of small examples. A virgin 50g has less than 241KB available. The solution is simple, do what others have done since the digital dawn of time, dump it to mass storage.

Create the following UserRPL program and save it as INSTALL in home/hpgcc. Next create a Custom Menu CST variable with the contents {INSTALL} in the same directory.

%%HP: T(3)A(R)F(.);
\<< \-> p
  \<<
      p EVAL
      p PURGE
      S\->EXE
      "EXTEND/" p +
      3 \->TAG
      DUP
      PURGE
      STO
  \>>
\>>

This UserRPL program will push to the stack a compiled C program (but not a standalone), purge it from RAM, convert it to standalone, and save it to SD as EXTEND/PROGNAME. If PROGNAME exists, it will be replaced, and if the EXTEND directory does not exist, it will be created.

Purge everything from the home/hpgcc/test directory except HIWRAP, then transfer hi2.hp from your workstation to home/hpgcc/test. To install press .

Use the File Manager to examine the SD card. There should be an EXTEND directory. Within the EXTEND directory HI2 should have been installed.

NOTE: The SD FAT file system has a couple of restrictions:

A file and directory names must be in 8.3 format: SSSSSSSS.XXX. The extension can be omitted.
Names are not case sensitive. All the examples have been in lowercase (e.g. hi2.hp). The .hp gets dropped when transferred from workstation to 50g. Once saved to SD hi2 can be recalled as hi2 or HI2.

After moving your C binaries to SD all this is required to execute them is a minor update to the wrappers, e.g. for HIWRAP:

Old		New

Since there is no longer any naming conflict or confusion I recommend that your wrapper names closely match the C executable name.

Base Template

The above examples are a good foundation for your own code but are missing an important feature. Speed. Below is a better template to start with:

#include <hpgcc49.h>

int main(void)
{
    //your variables

    //optional clear screen for direct text output
    //clear_screen();

    //full speed ahead
    sys_slowOff();

    //your critical code here

    //back to slow
    sys_slowOn();

    //optional notify/pause for cancel
    //beep();
    //WAIT_CANCEL;

    return(0);
}

The 50g by default runs at 75 MHz. This is probably a requirement to get good Saturn emulation performance. HPGCC by default runs at 12 MHz, but can also run at 75 MHz. I did two different tests comparing native Saturn code to compiled C code and found that at 12 MHz, HPGCC was approximately 15 times faster, and at 75 MHz, approximately 90-100 times faster. There is a performance and battery life trade off you will need to consider. I speculate that most will want to run at full speed.

The calls sys_slowOff() and sys_slowOn() are used to change between 75 MHz and 12 MHz. After declaring all your variables call sys_slowOff(). This will clock the 50g to max speed and max battery burn. Write your critical code, then end with sys_slowOn().

If you decide to pause with WAIT_CANCEL or any other input request, I recommend that you call sys_slowOn() and beep() first. This is important because the 50g will not auto power off outside of UserRPL. Perhaps during a lengthy computation you get distracted, beep() will help you get refocused. And, if for some reason you left your 50g unattended at least its back to slow mode conserving battery power.

In general consider the uses of sys_slowOn(), sys_slowOff(), and beep() throughout your code if frequently pausing for input, e.g. a chess game program. Why burn batteries while you are thinking?

Extend your 50g with C

To truly extend the 50g you need to create a library. Once this library is created you can easily share your C programs with others.

In this last introductory example you will create an Arithmetic Geometric Mean (AGM) function for the 50g. This is an example of what not to do in C because this function converges very quickly (the UserRPL version is as fast as the C version), but, it is useful as an example.

agm.c

AGM UserRPL

#include <hpgcc49.h>

int main(void){
    double a,b,c;

    sys_slowOff();

    a = sat_pop_real();
    b = sat_pop_real();

    if(a < 0 || b < 0) {
        sat_push_real(a);
        sat_push_real(b);
        sat_push_real(1);
        sys_slowOn();
        return(0);
    }

    while (fabs(a-b) > 1e-10) {
        c = a;
        a = (a + b)/2;
        b = sqrt(b * c);
        if(keyb_isAnyKeyPressed())
            break;
    }

    sat_push_real((a+b)/2);
    sat_push_real(0);

    sys_slowOn();
    return(0);
}

\<< \-> a b
  \<<
    WHILE a b - ABS .0000000001 >
    REPEAT
      a b + 2. /
      a b * \v/
      'b' STO
      'a' STO
    END
    a b + 2. /
  \>>
\>>

The AGM function takes two reals, computes the AGM, and returns a single real. AGM approximately doubles in precision after each iteration. Only four iterations are needed to exhaust the precision of 64 bit (double) floats. However, do not make this assumption, the while test condition in both programs depends on the values of a and b. Certain values of a and b have the potential to loop forever. You have a few options to address this:

Give yourself an out. In the C program keyb_isAnyKeyPressed is used to check for input and break out of the loop and return intermediate results. Since UserRPL programs can be interrupted with CANCEL(), it is not necessary for the UserRPL version to explicitly check for input.
Prove mathematically that for any pair of 64 bit floating point numbers that there will be a small finite number of iterations required to converge.
Hard code the number of loops to four or hard code a maximum number of iterations. Taking the time to better understand the mathematical properties of the algorithm can go a long way in space and speed optimization. This has been critical throughout the ages in calculator algorithm design--keep with tradition.

Since the sqrt function cannot take the square root of a negative number it is necessary to check for negative values, this could have also been done in the wrapper.

Compile and install to SD agm.c using the previous instructions above. Next create the AGM wrapper:

%%HP: T(3)A(R)F(.); \<< "Usage: Input 2 positive reals." \-> u \<<	Define a local variable u as our universal error message for UserRPL and C.
DEPTH IF 2 < THEN u DOERR END	Check for a stack depth less than two, if true, then put u on the stack, call DOERR and end (DOERR implicitly ends the UserRPL program).
\->NUM SWAP \->NUM	Convert the first two stack elements to reals using ->NUM. ->NUM is very handy to precondition input for C programs. ->NUM never errors if it cannot convert to a real. And ->NUM will convert symbolic inputs to reals (e.g. ).
DUP2 TYPE SWAP TYPE +	Duplicates the two stack items, and sums their type values.
IF THEN u DOERR END	The type for reals is 0, if the sum is not 0, then put u on the stack and call DOERR.
:3: "EXTEND/AGM" EVAL	If you got this far you have 2 reals on the stack, call EXTEND/AGM.
IF THEN u DOERR END \>> \>>	Check for nonzero return code, if nonzero, then u DOERR. agm.c only returns 1 if a or b is negative. agm.c also pushes a and b back on stack to be consistent with other 50g errors, i.e. the values that causes the error are left on the stack.

Non-annotated version:

%%HP: T(3)A(R)F(.);
\<< "Usage: Input 2 positive reals." \-> u
  \<<
    DEPTH
    IF 2 < THEN u DOERR END
    \->NUM SWAP \->NUM
    DUP2 TYPE SWAP TYPE +
    IF THEN u DOERR END
    :3: "EXTEND/AGM" EVAL
    IF THEN u DOERR END
  \>>
\>>

Why not use C to do all of this?

Library needs a wrapper.
Currently the HPGCC code to validate reals on the stack is broken.
More work to convert symbolic numbers (e.g. ).
No performance benefit.
Could make C code longer and buggier.

I recommend that libraries for C programs only contain the wrappers and not the C executables. Libraries are stored in ports 1 or 2. Both have limited space.

To create a library first define and store the following variables (yes the $ ( in ALPHA mode) is part of the variable name):

Variable	Value	Description
$TITLE	"EXTEND"	This is the title of the library. The first five characters will be used for the name that is shown in the library menu.
$ROMID	1313	This is a number (real or integer) which must be unique to your library. This allows the calculator to identify the library, and should be in the range 769 to 1792.
$CONFIG	1	Default configuration.
$VISIBLE	{AGM}	This is a list of the commands in the current directory that you want to have in the library and visible in the library's menu.

Your home/hpgcc/test should look like this:

To create the library and store (not install) to SD card type into the 50g:

256 ATTACH

CRLIB

:3:EXTEND.LIB

Your stack should look like this:

Save to SD Card by pressing . Your SD card now has an EXTEND directory and an EXTEND.LIB file. If you wish to share your C code and wrappers then just zip up both EXTEND and EXTEND.LIB as EXTEND.ZIP and distribute. The installation instructions are very simple (you can skip 1 and 2 since you already have it on your 50g):

On your PC extract EXTEND.ZIP into the root of the SD Card.
Insert SD Card into 50g.
Type into the 50g:

:3:EXTEND.LIB
2

Now press , then press +.
Test. Try a few pairs of numbers, try negative numbers, try symbolics, try variables.

Congratulations you just added AGM to the 50g vocabulary. If you like, you can remove the HOME/HPGCC/TEST 50g directory, it is not needed to run AGM. If you enter 1313 MENU in your 50g you can see all new commands this library adds.

The Computational Geometry section has a larger example including how to add an entry to the APPS menu.

Debugging

Debugging C code is always a challenge with the proper tools and a real challenge with out them. HPGCC includes no traditional debuggers (e.g. GDB). Actually HPGCC does include the original early 80's debugger: printf. If you are not laughing then you are new to C.

Debugging Options:

Prototype on workstation first. That's the beauty of C, it runs on anything. Algorithm design, file I/O, etc... can be tested on your workstation first, then ported to the 50g limiting the debugging to HPGCC bugs.
printf. Print state information to the screen, this has been a tried and true method of debugging ever since computers could output to printer or screen.
fprintf. File printf. You can put a lot of information in a file for analysis on your workstation.
Use the HPAPINE simulator[4]. This tool can really speed up development and testing. It cannot debug all HPGCC bugs, but it can catch most of them, and, since it runs on your workstation, tools like GDB are available.

HPAPINE 101

HPAPINE is an implementation of the HPGCC's API. Usually, a HP49g+/50g program written in C is compiled by HPGCC. With HPAPINE, you can compile the same program so that it runs natively on your Unix system. [4]

HPAPINE prebuilt for Cygwin/X is included as part of the hpgcc.tgz tarball. Linux and OS/X users will need to obtain the source and build it themselves.

NOTE: XP users must have SP2 or later installed to build HPAPINE binaries.

To give HPAPINE a test drive type the following in a new Xterm session:

cd ~/hpgcc/test

make clean

make hello_local

make agm_local

Your session should look something like this (inputs are in bold):

$ cd ~/hpgcc/test
$ make clean
rm -f *.o
rm -f *.hp
rm -f *.exe
$ make hello_local
gcc -Wall -g -o hello_local hello.c -DHPGCC -I ~/hpgcc/2.0SP2-hpapine/include/ \
        -DHPAPINE -L ~/hpgcc/2.0SP2-hpapine/lib/ -lhpapine \
        -L/usr/X11R6/lib -lgdk-x11-2.0 -lgthread-2.0 -lgdk_pixbuf-2.0 -lpangoxft-1.0 \
        -lXft -lfreetype -lz -lXrender -lXext -lfontconfig -lpangox-1.0 -lX11 \
        -lpango-1.0 -lm -lgobject-2.0 -lgmodule-2.0 -lglib-2.0 -lintl -liconv  
$ make agm_local
gcc -Wall -g -o agm_local agm.c -DHPGCC -I ~/hpgcc/2.0SP2-hpapine/include/ \
        -DHPAPINE -L ~/hpgcc/2.0SP2-hpapine/lib/ -lhpapine \
        -L/usr/X11R6/lib -lgdk-x11-2.0 -lgthread-2.0 -lgdk_pixbuf-2.0 -lpangoxft-1.0 \
        -lXft -lfreetype -lz -lXrender -lXext -lfontconfig -lpangox-1.0 -lX11 \
        -lpango-1.0 -lm -lgobject-2.0 -lgmodule-2.0 -lglib-2.0 -lintl -liconv

Now run run hello_local. Session output (inputs are in bold):

$ ./hello_local
HPAPINE version 20071113_134757 (HPGCC 2.0)
  Compiled on Sat Dec  8 18:45:20 MST 2007
  CYGWIN_NT-5.1 1.5.24(0.156/4/2) i686 (oberon)
  STRICT_HPGCC: no
  GUI: GDK (X11)

A new window will pop up with the simulated 50g output:

This window is identical to the what was displayed on the 50g. If you still have hello.hp installed on your 50g see for yourself. Interacting with this window is identical to interacting with the 50g. The letters on your keyboard match up with the numeric, math operator, and alpha keys from the 50g. is simulated with <ESC>. Press <ESC> now to exit the application.

HPAPINE can also simulate the stack for input and output from 50g programs. To see how this works run agm_local. Session output (inputs are in bold):

$ ./agm_local
HPAPINE version 20071113_134757 (HPGCC 2.0)
  Compiled on Sat Dec  8 18:45:20 MST 2007
  CYGWIN_NT-5.1 1.5.24(0.156/4/2) i686 (oberon)
  STRICT_HPGCC: no
  GUI: GDK (X11)

--- Waiting for stack input on stdin...
1.456
7.4563
--- Stack reading: done.

--- RPL Stack:

 2: 1.456
 1: 7.4563


HPAPINE version 20071113_134757 (HPGCC 2.0)
  Compiled on Sat Dec  8 18:45:20 MST 2007
  CYGWIN_NT-5.1 1.5.24(0.156/4/2) i686 (oberon)
  STRICT_HPGCC: no
  GUI: GDK (X11)

--- RPL Stack:

 2: "          3.85362"
 1: "                0"

This is exactly what the 50g would have returned without the wrapper. The result in stack position 2 is our answer (verify on your 50g), and the 0 in position 1 was the return code. NOTE: When entering values into the stack, enter them as if you were using a 50g, to terminate entry and continue running the application use <Ctrl-D>.

Next try using GDB and hello_local. Session output (inputs are in bold):

$ gdb hello_local
GNU gdb 6.5.50.20060706-cvs (cygwin-special)
Copyright (C) 2006 Free Software Foundation, Inc.
GDB is free software, covered by the GNU General Public License, and you are
welcome to change it and/or distribute copies of it under certain conditions.
Type "show copying" to see the conditions.
There is absolutely no warranty for GDB.  Type "show warranty" for details.
This GDB was configured as "i686-pc-cygwin"...
(gdb) start
Breakpoint 1 at 0x401075: file hello.c, line 4.
Starting program: /cygdrive/c/home/egan/hpgcc/test/hello_local.exe 
Loaded symbols for /cygdrive/c/WINDOWS/system32/ntdll.dll
...
Loaded symbols for /usr/bin/cyggthread-2.0-0.dll
main ()
    at hello.c:4
4       {

At this point GDB displays the next line to be executed. Use 'n<Enter>' to step through GDB.

(gdb) n
5               clear_screen();
(gdb) n

After clear_screen(); executes a blank hpapine window pops up:

6               printf("hello, world\n");
(gdb) n
7               WAIT_CANCEL;
(gdb) n

You'll notice the GDB does not return a (gdb) prompt, it is waiting for program input. Press <ESC> in the hpapine window.

8               return(0);
(gdb) n
9       }
(gdb) n
0x61006198 in dll_crt0_1 () from /usr/bin/cygwin1.dll
(gdb) n
Single stepping until exit from function _Z10dll_crt0_1Pv, 
which has no line number information.
HPAPINE version 20071113_134757 (HPGCC 2.0)
  Compiled on Sat Dec  8 18:45:20 MST 2007
  CYGWIN_NT-5.1 1.5.24(0.156/4/2) i686 (oberon)
  STRICT_HPGCC: no
  GUI: GDK (X11)

Program exited normally.
(gdb) q

Great stuff! For more information on GDB read: http://sourceware.org/gdb/current/onlinedocs/gdb.pdf.gz.

Just a few warnings about HPAPINE:

It is alpha level code.
Sometimes code works with HPAPINE but does not work on the 50g.
Sometimes code works with the 50g but not HPAPINE.
The stack output can be a bit twitchy.
Some libraries may not have HPAPINE support (e.g. win, hpstack) and may need to be compiled or not used.

You mileage may vary. I will state that HPAPINE has saved me hours and hours of work when developing HPGCC applications. Use it.

Makefiles

Makefiles are plain text files with instructions on how to make something. All of the Makefiles in this article have been created for you. The easiest way to create a Makefile is to copy it from another source and modify it. The Makefile in ~/hpgcc/test is a Makefile I modified from the HPGCC distribution. This Makefile is all you need if you are developing single source file binaries using only HPGCC supplied libraries.

E.g. lets say you wrote the program foo.c. To compile it for the 50g you would copy foo.c to ~/hpgcc/test and then type:

$ make foo.hp

arm-elf-gcc -mtune=arm920t -mcpu=arm920t -mlittle-endian -fomit-frame-pointer -msoft-float -Wall -Os -I/home/egan/hpgcc/2.0SP2/include -L/home/egan/hpgcc/2.0SP2/lib -mthumb-interwork -mthumb -c hello.c
arm-elf-ld -L/home/egan/hpgcc/2.0SP2/lib -T VCld.script /home/egan/hpgcc/2.0SP2/lib/crt0.o hello.o -lhpg -lhplib -lgcc -o foo.exe
elf2hp -s3000 foo.exe foo.hp
rm foo.o foo.exe

Say you wanted to tested it first with HPAPINE, then type:

$ make foo_local

gcc -Wall -g -o foo_local foo.c -DHPGCC -I /cygdrive/c/home/egan/hpgcc/2.0SP2-hpapine/include/ \
-DHPAPINE -L /cygdrive/c/home/egan/hpgcc/2.0SP2-hpapine/lib/ -lhpapine \
-L/usr/X11R6/lib -lgdk-x11-2.0 -lgthread-2.0 -lgdk_pixbuf-2.0 -lpangoxft-1.0 -lXft -lfreetype -lz -lXrender -lXext -lfontconfig -lpangox-1.0 -lX11 -lpango-1.0 -lm -lgobject-2.0 -lgmodule-2.0 -lglib-2.0 -lintl -liconv

Throughout this article I will start with this Makefile and modify it as needed. I will comment on the changes. For detailed information about Makefiles consider reading Managing Projects with make by Oram and Talbott.

Support

You are not alone. There are a few avenues for support:

The official HPGCC home page: http://hpgcc.org.
The official HPGCC mailing list. To subscribe go to https://lists.sourceforge.net/lists/listinfo/hpgcc-devel. To read the archives go to http://sourceforge.net/mailarchive/forum.php?forum_name=hpgcc-devel.
Search first, then post to the Usenet group comp.sys.hp48. If you do not have an NNTP reader and account you can search and post freely here: http://groups.google.com/group/comp.sys.hp48/topics.
A bunch of helpful and friendly calculator nerds can be found here: http://www.hpmuseum.org/cgi-sys/cgiwrap/hpmuseum/forum.cgi. Go to http://www.cafepress.com and search for calculator if you want to see what these nerds wear. You may spot one in the wild. If you do, approach slowly, make no sudden movements, then politely ask for help.
Google.
Use the source. HPGCC is free as is free beer and free speech.

Summary

Hopefully you have enough information to be dangerous and possibility productive. Where do you go from here? I would suggest that you start by looking other people code. The ~/hpgcc/2.0SP2/examples directory on your system has a lot of great examples. The 2nd half of this article has more complex examples involving the use of existing open source C libraries and source code.

Good luck!

Extend your 50g with C - Part 2

Part one of this article assumes little to no knowledge of HPGCC and C. This section will focus on learning by example by illustrating how to create complete mathematical extensions to the 50g. The reader is assumed to know how to compile C code with HPGCC, install binaries to SD, and create wrappers. All of this was covered in Part 1.

How to Follow the Examples

All binaries should be installed on the SD Card in the EXTEND directory.
$ prompts indicate something you should type. The text to be entered will be in bold courier and the output in courier, e.g.:

$ hello
go away!
All the code is in the ~/hpgcc directory. Each example has it's own subdirectory. Files ending in .c.nl are line numbered version of the .c code. My comments will reference line numbers. Line numbers were generated with:

perl -pi -e 's/\t/ /g' < program.c | nl -nln -ba >program.c.nl
All C programs will be in lowercase courier, e.g. program. All RPL wrappers will be in uppercase courier, e.g. PROGRAM.
To maintain consistency most code was reformatted with:

indent -kr -ts4 *.c *.h

The Examples

Real and Complex LogGamma. This is an example of how to do complex operations in C. C, unlike C99, C++, and Fortran, has no native support for complex numbers. Furthermore, LogGamma is its own function and is not exactly the same as the 50g ln(gamma()).
Sparse Linear Solver. This is an example of a very fast sparse linear solver capable of solving sparse systems too large for the 50g internal dense linear solver. This example will also demonstrate how to work with lists and arrays.
π Shootout. This example adds no immediate benefit to the 50g, but its fun. This example illustrates how to use an arbitrary multiple precision math library and a large integer math library. This example will also demonstrate how to work with files and text screen output.
Computational Geometry Library. This is an example of how to create static and interactive CG applications from freely available CG C code (e.g., convex hull (10,000 points in seconds), Voronoi diagrams, etc...). This example will also demonstrate how to take screen shots (B&W and grayscale). This example ends by creating a redistributable CG library with an APPS menu.

Example: Real and Complex LogGamma

You might be asking yourself, "Why not use the built in 50g functions ln() and gamma() together to compute LogGamma()?"

There are two reasons:

ln(gamma(z)) is not exactly the same as LogGamma(z). E.g.:

ln(gamma(12.3+4.56i)) = 17.38 - 1.20i

whereas

LogGamma(12.3+4.56i) = 17.38 + 11.36i

A difference of 2πn where n = 2. Since ln(z) = ln|z| + iarg(z), and arg(z) is in the interval [-π,π), the imaginary part will always be mod 2π. IOW, ln(gamma(z)) and LogGamma(z) have different branch cut structures. Refer to [5] for details.
Because LogGamma() is its own function and not the combination of ln() and gamma() large values for z can be computed. E.g. 50g ln(gamma(500)) = 1151.29, whereas LogGamma(500) = 2605.12. Any value greater than 254 on the 50g will always return 1151.29.

The HPGCC working directory for this example is in ~/hpgcc/complex. The 50g working directory is assume to be HOME/HPGCC/COMPLEX.

Building c9x-complex

Unfortunately HPGCC C has no native support for complex numbers. Fortunately there are freely available complex C libraries. The c9x-complex library is very complete and includes a LogGamma() function. 99% of the work is already done!

Change to working directory:

$ cd ~/hpgcc/complex
Download and extract c9x-complex.tgz:

$ wget http://www.moshier.net/c9x-complex.tgz
$ tar zxvf c9x-complex.tgz
Create stubs. This is something that you will do frequently. You have two choices, you can edit every file and change stdio.h and possibly other files to hpgcc49.h or you can create your own stdio.h that includes hpgcc49.h. The second method is preferred. It is best to alter the code as little as possible so that you can benefit from updates and get more support from the authors and the community.

$ cd c9x-complex
$ mkdir stubs
$ echo '#include <hpgcc49.h>' >stubs/stdio.h

Create ~/hpgcc/complex/c9x-complex/Makefile.hpgcc. This is a hybrid of a standard HPGCC Makefile and the Makefile included with c9x-complex. By examining both Makefiles you can see that most of the differences are in the compiler and compiler flag settings. Also notice the -Istubs added to the standard HPGCC CFLAGS. When stdio.h is encounter by the C preprocessor it will search and find stubs/stdio.h and include that. That will then include hpgcc49.h.

INCLUDE_PATH= $(HPGCC)/include
CC= arm-elf-gcc
AR= arm-elf-ar
RANLIB= arm-elf-ranlib
CFLAGS= -mtune=arm920t -mcpu=arm920t \
    -mlittle-endian -fomit-frame-pointer -msoft-float -Wall \
    -Os -I$(INCLUDE_PATH) -I. -Istubs -mthumb-interwork -mthumb

DFILES = cmplx.o clog.o cgamma.o
LIBMCFILES = $(DFILES) stubs.o

all: libmc.a

libmc.a: $(LIBMCFILES)
        rm -f libmc.a
        $(AR) cr libmc.a $(LIBMCFILES)
        $(RANLIB) libmc.a

clean:
        rm -f *.o
        rm -f libmc.a

Make it:

$ make -f Makefile.hpgcc clean
rm -f *.o
rm -f libmc.a
$ make -f Makefile.hpgcc
arm-elf-gcc -mtune=arm920t -mcpu=arm920t -mlittle-endian -fomit-frame-pointer -msoft-float -Wall -Os -I/home/egan/hpgcc/2.0SP2/include -I. -Istubs -mthumb-interwork -mthumb   -c -o cmplx.o cmplx.c
arm-elf-gcc -mtune=arm920t -mcpu=arm920t -mlittle-endian -fomit-frame-pointer -msoft-float -Wall -Os -I/home/egan/hpgcc/2.0SP2/include -I. -Istubs -mthumb-interwork -mthumb   -c -o clog.o clog.c
arm-elf-gcc -mtune=arm920t -mcpu=arm920t -mlittle-endian -fomit-frame-pointer -msoft-float -Wall -Os -I/home/egan/hpgcc/2.0SP2/include -I. -Istubs -mthumb-interwork -mthumb   -c -o cgamma.o cgamma.c
arm-elf-gcc -mtune=arm920t -mcpu=arm920t -mlittle-endian -fomit-frame-pointer -msoft-float -Wall -Os -I/home/egan/hpgcc/2.0SP2/include -I. -Istubs -mthumb-interwork -mthumb   -c -o stubs.o stubs.c
In file included from mconf.h:177,
                 from stubs.c:5:
protos.h:85: error: expected ')' before '/' token
make: *** [stubs.o] Error 1

Doh! What's this hard to understand error? Line 85 in protos.h doesn't look out of the ordinary:

extern double log2 ( double x );

What's up? Well, log2 conflicts with the log2 macro in hpmath.h:

#define log2(x) (log(x) / M_LOG2_E)

One or the other has to yield. For now, just comment out line 85 in the protos.h with // and try to make again:

$ make -f Makefile.hpgcc
arm-elf-gcc -mtune=arm920t -mcpu=arm920t -mlittle-endian -fomit-frame-pointer -msoft-float -Wall -Os -I/home/egan/hpgcc/2.0SP2/include -I. -Istubs -mthumb-interwork -mthumb -c -o cmplx.o cmplx.c
arm-elf-gcc -mtune=arm920t -mcpu=arm920t -mlittle-endian -fomit-frame-pointer -msoft-float -Wall -Os -I/home/egan/hpgcc/2.0SP2/include -I. -Istubs -mthumb-interwork -mthumb -c -o clog.o clog.c
arm-elf-gcc -mtune=arm920t -mcpu=arm920t -mlittle-endian -fomit-frame-pointer -msoft-float -Wall -Os -I/home/egan/hpgcc/2.0SP2/include -I. -Istubs -mthumb-interwork -mthumb -c -o cgamma.o cgamma.c
arm-elf-gcc -mtune=arm920t -mcpu=arm920t -mlittle-endian -fomit-frame-pointer -msoft-float -Wall -Os -I/home/egan/hpgcc/2.0SP2/include -I. -Istubs -mthumb-interwork -mthumb -c -o stubs.o stubs.c
rm -f libmc.a
arm-elf-ar cr libmc.a cmplx.o clog.o cgamma.o stubs.o
arm-elf-ranlib libmc.a

c9x-complex is now built.

Get a Complete libm.a

The math library in HPGCC is missing a few things. Fortunately Jean-Yves Avenard has built us a proper libm.a. All that is needed is to extract it and link with it.

Change to working directory:

$ cd ~/hpgcc/complex
Obtain JYA toolchain:

$ wget http://www.hydrix.com/Download/Hp/hpgcc/arm-elf-toolchain.4.2.2.linux.tar.bz2
Extract libm.a

$ tar jxvf arm-elf-toolchain.4.2.2.linux.tar.bz2 /opt/arm-hp/arm-elf/lib/thumb/interwork/libm.a
$ cp opt/arm-hp/arm-elf/lib/thumb/interwork/libm.a .
Clean up:

$ rm -r ./opt arm-elf-toolchain.4.2.2.linux.tar.bz2

LogGamma

With libmc.a and libm.a in hand implementing a fast Complex LogGamma for the 50g is relatively easy in 29 lines of code:

1       #include <hpgcc49.h>
2       #include <complex.h>
3
4       int __errno;
5
6       int main(void) {
7           double r,i;
8           double complex w,z;
9
10          sys_slowOff();
11
12          i = sat_pop_real();
13          r = sat_pop_real();
14
15          if(i == 0 && (r == 1 || r == 2)) {
16              sat_push_real(0);
17              sat_push_real(0);
18          }
19          else {
20              z =  r + i * I;
21              w = clgam(z);
22
23              sat_push_real(creal(w));
24              sat_push_real(cimag(w));
25          }
26
27          sys_slowOn();
28          return (0);
29      }

Lines 2 and 4 add the declarations required by the c9x-complex and standard math libraries.
Line 8 defines 2 complex variables.
Lines 12 and 13 pop off 2 reals from the stack. HPGCC does not have native support for popping off a complex number. An even if it could, you still have the complexities of making sure both parts are reals and not symbolic. This is best left to the wrapper.
Lines 15-18 check for some known values that should return 0, but that clgam() returns as very small reals.
Lines 20 and 21 is where the magic happens.
Lines 23 and 24 push the real and imaginary values to the stack to be processed by the wrapper.

Makefile

The Makefile for LogGamma is the same used in Part 1 with the following changes in bold:

export C_FLAGS= -mtune=arm920t -mcpu=arm920t \
    -mlittle-endian -fomit-frame-pointer -msoft-float -Wall \
    -Os -I$(INCLUDE_PATH) -Ic9x-complex -L$(LIBS_PATH) -mthumb-interwork -mthumb
export LD= arm-elf-ld
export LD_FLAGS= -L$(LIBS_PATH) -Lc9x-complex -L. -T VCld.script $(LIBS_PATH)/crt0.o
export LIBS= -lmc -lm -lhpg -lhplib -lgcc

-Ic9x-complex tells the C processor to also look for include files in the directory c9x-complex. This is where complex.h exists.
-Lc9x-complex -L. tells the linker to look also look for libraries in the c9x-complex and current directory.
-lmc -lm tells the linker to link in libmc.a and libm.a. The order of the libraries is important. If two libraries share a common function, the first one listed gets used. If a library relies on functions in another library, then the library that provides the specific function must be placed after the library requesting the functions. IOW, libmc.a needs functions from libm.a and all libraries need base functions from libgcc.a.

To build it type:

make clean

make loggamma.hp

Wrapper

loggamma UserRPL wrapper:

%%HP: T(3)A(R)F(.);
\<< \->NUM DUP TYPE
  IF 1 ==
  THEN
    C\->R
  ELSE
    0.
  END
  :3: "EXTEND/LOGGAMMA" EVAL
  DUP
  IF 0 ==
  THEN
    DROP
  ELSE
    R\->C
  END
\>>

This wrapper starts out by converting the input to a number, duplicating it, then converting the duplicate to its type. If type 1 (complex number), then breakout the number into is real and imaginary components, otherwise put a 0 on the stack as the imaginary part. Up to this point all that has been done is to insure that there is a complex number Re and Im on the stack in that order.

Next, extend/loggamma is executed and returns to the stack a complex number as Re and Im as two stack elements. The Im is duplicated to preserve its value and checked for 0. If 0 the value is dropped and the wrapper ends with a real (Re) result, otherwise the two reals (Re and Im) are combined as a type 1 complex number on the stack.

LogGamma Test

z	ln Γ(z)	50g loggamma(z)	50g ln(gamma(z))
1	0	0	0
100	359.134205370	359.13420537	359.13420537
1000	5905.22042321	5905.22042321	1151.2925465
5555.4444	42343.2	42343.1698379	1151.2925465
12+34i	-11.7232 + 102.052i	(-11.7232,102.0521)	(-11.7232,1.5212)

Incorrect answers in bold.

Summary

You now have everything you need to create single source complex number C binaries for your 50g. Just copy them in ~/hpgcc/complex and type make.

Example: Sparse Linear Solver

Inspiration for this example came from the following challenge:

"Three panes of glass are arranged parallel to each other. Each pane allows 70% light through, reflects 20%, and absorbs 10%. What percentage of light will pass entirely through all panes of glass? Assume a beam of light is emitted to the left of the three panes and exits to the far right; how much emits out to the right?" -- http://www.hpmuseum.org/cgi-sys/cgiwrap/hpmuseum/archv017.cgi?read=118862

Constants:

T = 0.7 (light through)

R = 0.2 (light reflected)

The solution (C) to this problem can easily be obtained from the following system of equations:

A = T + R*Y	Initial light through + reflected Y light.
B = TA + RX	A light through + reflected X light.
C = T*B	B light through.
Y = TX + RA	X light through + reflected A light.
X = R*B	Reflected B light.

The above 5 equations can also be represented by the following sparse linear system:

Solving for C with 3 panes is simple, just put it in your 50g. But what about 4 panes or 100? That was the question I asked myself when I first read this challenge. The solution to the 4 panes problem is similar.

But, what about 100? Clearly there is a pattern that can be generated for any number of panes. The following RPL program will generate a coefficient matrix and constant term vector for any number of panes based on this pattern.

Dense Panes (DPANES):

%%HP: T(3)A(R)F(.);
\<< 0. \-> t r n m                        @ Get T, R, and N from stack
  \<<
    n 2 * 1 - 'm' STO                     @ Compute M (matrix size)
    t NEG 0 m 1 - NDUPN 1 + \->ARRY       @ Create constant term vector 
    -1 m NDUPN \->ARRY { m m } DIAG\->    @ Create -1 diagonal coefficient matrix
    n 1 + m FOR i                         @ Add R elements to diagonal coefficient matrix
      i n - i 2 \->LIST r PUT
      i i n - 2 \->LIST r PUT
    NEXT
    1 n 1 - FOR i                         @ Add T elements to diagonal coefficient matrix
      i 1 + i 2 \->LIST t PUT
    NEXT
    n 1 + m 1 - FOR i
      i i 1 + 2 \->LIST t PUT
    NEXT
  \>>
\>>

To solve using the 50g dense linear solver press / then n GET for the answer, e.g.:

Dense Panes Solver (DPS):

%%HP: T(3)A(R)F(.);
\<< \-> t r n
  \<<
     t r n DPANES
     / n GET
  \>>
\>>

DPS can be used to solve this problem for any (small) N.

The following table was constructed with output from DPANES/DPS:

N	DPS(.7,.2,N)	Coeff. Bytes 16(2N)-1)^2	Non-Zero Bytes 166(N-1)	Zeros Bytes 16(4N^2-10*N+7)	Time(s) Setup	Time(s) Solve
3	38.0266075388%	400	192	208	0.38	0.17
4	28.6686567164%	784	288	496	0.45	0.33
10	5.6721177756%	5776	864	4912	1.31	1.91
20	.3945314511%	24336	1824	22512	3.96	7.98
30	.0274666640%	55696	2784	52912	8.24	20.10
40	.0019121946%	99856	3744	96112	15.95	36.60
50	.0001331246%	156816	4704	152112	30.67	59.42
100		633616	9504	624112
200		2547216	19104<