## Learn A New Language

### February 21, 2012

Here’s the beginning of the program with functions to store bitarrays as an array of characters and a basic list library; the idea is to keep things as simple as possible, so the list library is incomplete:

`/* prime.c -- programming with prime numbers`

* compile as: gcc -lgmp -o prime prime.c

* run as: ./prime -- takes no input */

```
```#include <stdio.h>

#include <malloc.h>

#include <string.h>

#include <gmp.h>

/* bit arrays */

#define ISBITSET(x, i) (( x[i>>3] & (1<<(i&7)) ) != 0)

#define SETBIT(x, i) x[i>>3] |= (1<<(i&7));

#define CLEARBIT(x, i) x[i>>3] &= (1<<(i&7)) ^ 0xFF;

/* basic list library */

typedef struct list {

void *data;

struct list *next;

} List;

List *insert(void *data, List *next)

{

List *new;

new = malloc(sizeof(List));

new->data = data;

new->next = next;

return new;

}

List *reverse(List *list) {

List *new = NULL;

List *next;

while (list != NULL)

{

next = list->next;

list->next = new;

new = list;

list = next;

}

return new;

}

`int length(List *xs)`

{

int len = 0;

while (xs != NULL)

{

len += 1;

xs = xs->next;

}

return len;

}

The Sieve of Eratosthenes we have seen before, but here I’ve used a somewhat different calculation of *j* to initialize the inner loop:

`/* sieve of eratosthenes */`

```
```List *primes(long n)

{

int m = (n-1) / 2;

char b[m/8+1];

int i = 0;

int p = 3;

List *ps = NULL;

int j;

ps = insert((void *) 2, ps);

memset(b, 255, sizeof(b));

while (p*p < n)

{

if (ISBITSET(b,i))

{

ps = insert((void *) p, ps);

j = (p*p - 3) / 2;

while (j < m)

{

CLEARBIT(b, j);

j += p;

}

}

i += 1; p += 2;

}

while (i < m)

{

if (ISBITSET(b,i))

{

ps = insert((void *) p, ps);

}

i += 1; p += 2;

}

` return reverse(ps);`

}

The `is_prime`

function doesn’t work; any help will be appreciated, either in the comments below or in private email. Fortunately, GMP provides this function natively under the name `mpz_probab_prime_p`

, so a non-working `is_prime`

doesn’t impact the integer factorization code:

`/* miller-rabin primality checker */`

```
```int is_even(mpz_t n)

{

mpz_t t;

mpz_init(t);

mpz_mod_ui(t, n, 2);

int r = mpz_cmp_ui(t, 0) ? 0 : 1;

mpz_clear(t);

return r;

}

void power_mod(mpz_t r, mpz_t b, mpz_t e, mpz_t m)

{

mpz_set_ui(r, 1);

while (mpz_cmp_ui(e, 0) > 0)

{

if (! is_even(e))

{

mpz_mul(r, r, b);

mpz_mod(r, r, m);

}

mpz_fdiv_q_ui(e, e, 2);

mpz_mul(b, b, b);

mpz_mod(b, b, m);

}

}

int is_spsp(mpz_t n, mpz_t a)

{

mpz_t d, n1, t;

mpz_inits(d, n1, t, NULL);

mpz_sub_ui(n1, n, 1);

mpz_set(d, n1);

int r = 0;

int s = 0;

while (is_even(d))

{

mpz_divexact_ui(d, d, 2);

s += 1;

}

power_mod(t, a, d, n);

if (mpz_cmp_ui(t, 1) == 0)

{

return 1;

}

while (r < s)

{

power_mod(t, a, d, n);

if (mpz_cmp(t, n1) == 0)

{

return 1;

}

r += 1;

mpz_mul_ui(d, d, 2);

}

return 0;

` mpz_clears(d, n1, t, NULL);`

}

Here is the function that provides integer factorization by trial division. It first extracts factors of 2, then divides by the odd numbers 3, 5, 7, 9, …:

`/* integer factorization by trial division */`

```
```void td_factors(List **fs, mpz_t n)

{

mpz_t d, d2, q, r;

mpz_inits(d, d2, q, r, NULL);

mpz_set_ui(d, 2);

while (is_even(n))

{

mpz_t *f = malloc(sizeof(*f));

mpz_init_set(*f, d);

*fs = insert(*f, *fs);

mpz_divexact(n, n, d);

}

mpz_set_ui(d, 3);

mpz_set_ui(d2, 9);

while (mpz_cmp(n, d2) >= 0)

{

mpz_fdiv_qr(q, r, n, d);

if (mpz_cmp_ui(r, 0) > 0)

{

mpz_add_ui(d, d, 2);

mpz_mul(d2, d, d);

}

else

{

mpz_t *f = malloc(sizeof(*f));

mpz_init_set(*f, d);

*fs = insert(*f, *fs);

mpz_divexact(n, n, d);

}

}

` if (mpz_cmp_ui(n, 1) != 0)`

{

*fs = reverse(insert(n, *fs));

}

mpz_clears(d, d2, q, r, NULL);

}

Integer factorization by Pollard’s rho method is somewhat more complicated. Although it appears here last, it was actually the first GMP function that I wrote, and I had all kinds of trouble getting it to work:

`/* integer factorization by pollard rho */`

```
```List *insert_in_order(void *x, List *xs)

{

if (xs == NULL || mpz_cmp(x, xs->data) < 0)

{

return insert(x, xs);

}

else

{

List *head = xs;

while (xs->next != NULL && mpz_cmp(x, xs->next->data) > 0)

{

xs = xs->next;

}

xs->next = insert(x, xs->next);

return head;

}

}

void rho_factor(mpz_t f, mpz_t n, long long unsigned c)

{

mpz_t t, h, d, r;

mpz_init_set_ui(t, 2);

mpz_init_set_ui(h, 2);

mpz_init_set_ui(d, 1);

mpz_init_set_ui(r, 0);

while (mpz_cmp_si(d, 1) == 0)

{

mpz_mul(t, t, t);

mpz_add_ui(t, t, c);

mpz_mod(t, t, n);

mpz_mul(h, h, h);

mpz_add_ui(h, h, c);

mpz_mod(h, h, n);

mpz_mul(h, h, h);

mpz_add_ui(h, h, c);

mpz_mod(h, h, n);

mpz_sub(r, t, h);

mpz_gcd(d, r, n);

}

if (mpz_cmp(d, n) == 0) /* cycle */

{

rho_factor(f, n, c+1);

}

else if (mpz_probab_prime_p(d, 25)) /* success */

{

mpz_set(f, d);

}

else /* found composite factor */

{

rho_factor(f, d, c+1);

}

mpz_clears(t, h, d, r, NULL);

}

void rho_factors(List **fs, mpz_t n)

{

while (! (mpz_probab_prime_p(n, 25)))

{

mpz_t *f = malloc(sizeof(*f));

mpz_init_set_ui(*f, 0);

rho_factor(*f, n, 1);

*fs = insert_in_order(*f, *fs);

mpz_divexact(n, n, *f);

}

` *fs = insert_in_order(n, *fs);`

}

Finally, here is a `main`

function that uses all of the functions given above:

`/* demonstration */`

```
```int main(int argc, char *argv[])

{

mpz_t n;

mpz_init(n);

List *ps = NULL;

List *fs = NULL;

printf("Primes less than a hundred:");

ps = primes(100);

while (ps != NULL)

{

printf(" %ld", (long) ps->data);

ps = ps->next;

}

printf("\n");

printf("Count of primes less than a million: %d\n",

length(primes(1000000)));

mpz_set_str(n, "600851475143", 10);

printf("Is_even(600851475143) is %d\n", is_even(n));

mpz_t b, e, m, t;

mpz_init_set_ui(b, 437);

mpz_init_set_ui(e, 13);

mpz_init_set_ui(m, 1741);

mpz_init(t);

power_mod(t, b, e, m);

printf("Power_mod of 437^13 (mod 1741) is %s\n", mpz_get_str(NULL, 10, t));

mpz_set_str(n, "600851475143", 10); /* composite */

printf("By 25 iterations of Miller-Random, 600851475143 is %s\n",

is_prime(n) ? "prime" : "composite");

mpz_set_str(n, "2305843009213693951", 10); /* prime */

printf("By 25 iterations of Miller-Random, 2305843009213693951 is %s\n",

is_prime(n) ? "prime" : "composite");

printf("Factors of 600851475143 by trial division:");

mpz_set_str(n, "600851475143", 10);

td_factors(&fs, n);

while (fs != NULL) {

printf(" %s", mpz_get_str(NULL, 10, fs->data));

fs = fs->next;

}

printf("\n");

` printf("Factors of 600851475143 by Pollard rho:");`

mpz_set_str(n, "600851475143", 10);

rho_factors(&fs, n);

while (fs != NULL) {

printf(" %s", mpz_get_str(NULL, 10, fs->data));

fs = fs->next;

}

printf("\n");

}

You can see the program at http://programmingpraxis.codepad.org/ZfVUDHVf.

I’m a little surprised (and humbled) to be the first to post an answer to this challenge. Back when I was in College, I learned about Touring Completeness. At the time, I didn’t give it much thought until about 5 or so years ago, I stumbled across a programming language called BrainF**k (BF for short). I didn’t care for the name, but the concept was truly fascinating. The language simulates a one tape Touring Machine using only 8 instructions which manipulate one pointer and tape data under that pointer. In theory, this is all that’s needed to write any program in the world. When I saw this 5 years ago, I wrote a Hello World to prove the concept to myself, and didn’t give it another thought … until now.

I decided to write a program which would output the HailStone sequence defined by the previous exercise, using BF. Well, almost. I actually used Procedural BF (pBrain for short) which slightly extends the instruction sets to add capability to define and call functions in 3 additional instructions. You can read more about pBrain at http://www.parkscomputing.com/code/pbrain/. To simplify my work, I made a couple of assumptions. The user input of the starting number must be between 01 and 99, always entered as two characters with a leading zero. The output will always be displayed as a 3 digit number with leading zeros. This touring machine has an 8-bit memory, so no single number within the sequence must exceed 255.

Without further ado, here is my Touring Machine (with procedural enhancements) implementation of a HailStone sequence. There are plenty of comments in the code. Frankly, I needed them to keep myself sane. Of course, the nice side effect of the comments is that everyone should be able to understand what is going on. As an example, when run with an input like 07, the output is 007 022 011 034 017 052 026 013 040 020 010 005 016 008 004 002 001

/**************************************/

/* HAILSTONE SEQUENCE PROGRAM */

/* in pBrain */

/**************************************/

// Memory map

// 0: ASCII_MSB

// 1: ASCII_LSB

// 2: User_Input

// 3: R1

// 4: R2

// 5: R3

// 6: R4

// 7: R5

// 8: R6

// 9: R7

// A: R8

// B: R9

// C: R10

// D: R11

// E: R12

// F: R13

// 10: R14

// 11: R15

// 12: R16

/**************************************/

/* Function 1 – Convert ASCII(n) to integer N by subtracting 6*8=ASCII(0) */

/* Input: @(-1) – the ASCII value of n */

/* Output: @(-1) – the integer value N */

/* Local: @(+1) and @(+2) */

+(

// loop 6*8 times, decrementing @(-1)

>++++++++[

>++++++[

<<<-

>>>-

]<-

]

// Exit with calling register set back to zero

<-

)

/*************************************/

/* Function 2 – Divide A/B */

/* Input: @(-2) – Divident A */

/* Input: @(-1) – Divisor B */

/* Local: @(+1) – Copy of dividend */

/* Local: @(+2) – Copy of dividend */

/* Local: @(+3) – Remainder */

/* Local: @(+4) – Copy of divisor */

/* Local: @(+5) – Quotient */

/* Local: @(+6) – Zero */

/* Local: @(+7) – Zero */

/* Output: @(+5) – Quotient */

/* Output: @(+3) – Remainder */

+(

// make sure the local variables are zero’d out

>[-]>[-]>[-]>[-]>[-]>[-]>[-]<<<<<<<

// copy dividend to @(+1) and @(+2) destroying @(-2)

<<[>>>+>+<<<<-]

// copy divisor to @(+4) destroying @(-1)

>[>>>>>+<<<<<-]

// point to first copy of dividend and start division

>>

[

>+>+>-[>>>]

<[[>+<-]>>+>]

<<<<<-

]

// Exit with calling register set back to zero

<–

)

/**************************************/

/* Function 3 – display integer 0-255 as ASCII */

/* Input: @(-1) – the ASCII value of n */

/* Local: @(+1), @(+2), @(+3) – ASCII result */

/* Local: @(+4) through @(+13) – used in division */

+(

// Make sure that we start with a clean slate for internal output registers

[-]>[-]>[-]>[-]

// Figure out the 100’s digit.

// Copy number to divident

<<<<[>+>>>>+<<<<<-]>[<+>-]

// Set divisor to 10*10=100

>>>>>>++++++++++[>++++++++++[<<+>>-]<-]<

// Call division function

>++:

// Process results of division

>>>>>

[

// Copy quotient into MSB

<<<<<<<<+

// subtract quotient*100 from number

>++++++++++ [>++++++++++ [<<<<<<->>>>>>-]<- ]

>>>>>>>-

]

// Convert to ascii by adding 6*8=48

<<<<<<[-]<[-]++++++ [ >++++++++ [<<+>>-]<- ]

// Figure out the 10’s digit

>>>[-]<[-]<[-]

// Copy number to divident

<<<<<<[>+>>>>+<<<<<-]>[<+>-]

// Set divisor to 2*5=10

>>>>>>++[>+++++[<<+>>-]<-]<

// Call division function

>++:

// Process results of division

>>>>>

[

// Copy quotient into Middle digit

<<<<<<<<<+

// subtract quotient*10 from number

>>++ [ >+++++ [<<<<<<->>>>>>-]<- ]

>>>>>>>-

]

// Convert to ascii by adding 6*8=48

<<<<<<[-]<[-]++++++ [ >++++++++ [<<<+>>>-]<- ]

// Copy the remainder (1’s digit) into LSB

>>>>>[<<<<<<<<+>>>>>>>>-]

// Convert to ascii by adding 6*8=48

<<<<[-]<[-]++++++ [ >++++++++ [<<<<+>>>>-]<- ]

// Display the number with a space after it

<.<.<.[-]>[-]>[-]< ++++ [ >++++++++ [<<+>>-]<- ] <.

// exit with calling register set to zero

<[-]

)

/**************************************/

/**************** MAIN ****************/

// Get user input for MSB and LSB stored in @ASCII_MSB and @ASCII_LSB

,>,

// Move @ASCII_MSB to R1

<[>>>+<<<-]

// Call Func1 from R2, using R1 as arg1, R3 and R4 as local vars

>>>>+:

// Move @ASCII_LSB to R2

<<<[>>>+<<<-]

// Call Func1 from R3, using R2 as arg1, R4 and R5 as local vars

>>>>+:

// @User_Input = R1*10 + R2, destroying R1, R2 contents

<<[<++++++++++>-]

>[<<+>>-]

// Copy @User_input to R1 using R2

// Display R1 by calling Func3 from R2, using R3-R16 as local vars

<< [>+>+<<-] >>[<<+>>-] +++:

// While @User_Input != 1

<<-

[

// Copy @User_Input to R1, using R2, set R2 to two

+>[-]<[>+>+<<-]>>[<<+>>-]++

// Call Func2 from R3, using R1 and R2 as arg1 and arg2, and R4-R10 as local vars

>[-]++:

// Test R6 (remainder) to see if it’s Zero, use R7 for starting else clause

>>>>[-]+<

[

// Remainder is not Zero, must have been an odd number

// @user_Input = 3*@User_input + 1

<<<<<[-]<[>+++<-]>+[<+>-]

// Set remainder and R7 to zero to skip else clause

>>>>>>-<-

]

>

[

// Remainder is Zero, must have been an even number

// Copy R8 (quotient) to @User_input

<<<<<<<[-]>>>>>>>>[<<<<<<<<+>>>>>>>>-]<[-]

]

// Copy @User_input to R1 using R2

// Display R1 by calling Func3 from R2, using R3-R16 as local vars

<<<<<<< [>+>+<<-] >>[<<+>>-] +++:

// Test if @User_Input == 1

<<-

]

Ran into SmallBasic yesterday: a world without scope is pretty interesting…

‘ Reference:

‘ Cormen, Thomas H.; Leiserson, Charles E.; and Rivest, Ronald L

‘ Introduction to Algorithms

‘ The MIT Press and McGraw-Hill Book Company, 1990

‘ ISBN 0-262-03141-8 (MIT Press), ISBN 0-07-013143-0 (McGraw-Hill)

arraySize = 40

maxValue = 99

TextWindow.WriteLine(“Begin execution”)

PopulateArray()

TextWindow.WriteLine(“Before sorting:”)

ShowArray()

BuildHeap()

ShowHeap()

HeapSort()

TextWindow.WriteLine(“After sorting:”)

ShowArray()

Validate()

‘ Ensure no dupes are in the array

Sub PopulateArray

For i = 1 To maxValue

selected[i] = 0

EndFor

For i = 1 To arraySize

random = Math.GetRandomNumber(maxValue)

While(selected[random] > 0)

random = Math.GetRandomNumber(maxValue)

EndWhile

selected[random] = 1

a[i] = random

EndFor

EndSub

‘ Output 4 lines of 10 each, with leading 0 if necessary

Sub ShowArray

For index = 1 To arraySize

If (a[index] 0 And heap[j] 1) Then

parent = 1

While(parent*2 <= heapsize)

left = parent * 2

right = left + 1

swapIndex = parent

If (heap[left] < heap[swapIndex]) Then

swapIndex = left

EndIf

If (heap[right] < heap[swapIndex] And right a[j]) then

TextWindow.BackgroundColor = “Red”

TextWindow.WriteLine(“OH NO!”)

Endif

Endfor

Endfor

EndSub

Still very, very new to Clojure. I thought I would do lazy primes using leftist heap, + the leftist heap exercise, which is therefore quite a lot for (laziness, data structures, etc.)

First the left heap:

Then the lazy primes with O’Neil’s algorithm:

Then to try it out: