## 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.

Pages: 1 2 3

### 3 Responses to “Learn A New Language”

1. Paul G. said

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
<<-
]

2. Joe Taylor said

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

3. David said

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:

```(defn first
"Return the first value in a priority queue (heap)  If the
queue is empty, throw an exception."
[t] (t :value))

(defn merge
"Merge two heaps together, assuming both heaps may be compared
with the same ordering predicate.  It is used internally to implement
the other heap operations, but may be called by user code if desired."
[ord t1 t2]
(let [rank (fn [t]
(if (nil? t)
0
(t :rank))),

node (fn [val left right]
(if (> (rank right) (rank left))
{:value val, :rank (inc (rank left)),  :left right, :right left}
{:value val, :rank (inc (rank right)), :left left,  :right right}))
]
(cond
(nil? t1) t2
(nil? t2) t1
(ord (t1 :value) (t2 :value))
(node (t1 :value) (t1 :left) (merge ord (t1 :right) t2))
:else
(node (t2 :value) (t2 :left) (merge ord t1 (t2 :right))))))

(defn remove
"Return a new priority queue with the highest priority (based
on ordering predicate) removed.  Note: the predicate is used to
reorder the collection, but is NOT used to determine the top.
If the collection is empty, throws an exception.  Returns nil
when the last element is removed."
[ord t]  (merge ord (t :left) (t :right)))

(defn insert
"Insert an item in a priority queue, the ordering predicate is
used to determine the proper order for items in the queue."
[ord t val] (merge ord {:value val, :rank 1, :left nil, :right nil} t))
```

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

```; lazy sieve using priority queue

(defn ord [x y]  (< (first x) (first y)))

(def heap-insert    (partial leftist-heap/insert ord))
(def heap-remove    (partial leftist-heap/remove ord))
(def heap-create    (partial heap-insert nil))
(defn heap-top [h]  (first (leftist-heap/first h)))

(loop [h heap]
(if (= (heap-top h) n)
(let [[n step] (leftist-heap/first h)]
(recur (heap-insert (heap-remove h) [(+ n step) step])))
h)))

(defn first-prime []
[3 (heap-create [9 6])])

(defn next-prime [[n heap]]
(loop [n (+ 2 n), h heap]
(if (not= (heap-top h) n)
[n (heap-insert h [(* n n) (+ n n)])]
(recur (+ n 2) (adjust-heap h n)))))

(def lazy-primes
(cons 2 (map first (iterate next-prime (first-prime)))))
```

Then to try it out:

```user=> (take 168 lazy-primes)
(2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97 101 103 107 109 113 127 131 137 139 149 151 157
163 167 173 179 181 191 193 197 199 211 223 227 229 233 239 241 251 257 263 269 271 277 281 283 293 307 311 313 317 331
337 347 349 353 359 367 373 379 383 389 397 401 409 419 421 431 433 439 443 449 457 461 463 467 479 487 491 499 503 509
521 523 541 547 557 563 569 571 577 587 593 599 601 607 613 617 619 631 641 643 647 653 659 661 673 677 683 691 701 709
719 727 733 739 743 751 757 761 769 773 787 797 809 811 821 823 827 829 839 853 857 859 863 877 881 883 887 907 911 919
929 937 941 947 953 967 971 977 983 991 997)
```