ejolson
Posts: 3406
Joined: Tue Mar 18, 2014 11:47 am

Re: A Final Fibonacci Challenge

Tue Jun 11, 2019 5:46 pm

jahboater wrote:
Tue Jun 11, 2019 5:32 pm
ejolson wrote:
Tue Jun 11, 2019 5:24 pm
I've never used valgrind. How does
valgrind wrote:definitely lost: 607,744,423 bytes in 41,100,299 blocks
indicate no memory leaks? Or was your declaration of "clean" meant to be humorous?
See the first line (the heap summary):

in use at exit: 607,744,744 bytes in 41,100,305 blocks

When a process exits all allocated memory is discarded and reclaimed by the OS; the address space dissapears. True for all operating systems except Novel's Netware.
Few programs bother to release memory before calling exit as may add considerably to the complexity and achieves nothing.
I'll take that as a yes. What needs to be changed if I want to compute a bunch of Fibonacci numbers, say 1000, in a loop? Or use GMP in an extension to a long running program written in Script Basic?

Heater
Posts: 13084
Joined: Tue Jul 17, 2012 3:02 pm

Re: A Final Fibonacci Challenge

Tue Jun 11, 2019 5:51 pm

jahboater,
Remove the letis() function and replace all the calls to it with strdup(), then it compiles cleanly with GCC 9.1
I'm keeping letis. It says what it does in the context of the program "let integer string". I have replaced the guts with strdup though. Thanks.

Code is in the repo: https://github.com/ZiCog/fibo_4784969/b ... _strings.c

jahboater
Posts: 4598
Joined: Wed Feb 04, 2015 6:38 pm

Re: A Final Fibonacci Challenge

Tue Jun 11, 2019 5:52 pm

ejolson wrote:
Tue Jun 11, 2019 5:46 pm
jahboater wrote:
Tue Jun 11, 2019 5:32 pm
ejolson wrote:
Tue Jun 11, 2019 5:24 pm
I've never used valgrind. How doesindicate no memory leaks? Or was your declaration of "clean" meant to be humorous?
See the first line (the heap summary):

in use at exit: 607,744,744 bytes in 41,100,305 blocks

When a process exits all allocated memory is discarded and reclaimed by the OS; the address space dissapears. True for all operating systems except Novel's Netware.
Few programs bother to release memory before calling exit as may add considerably to the complexity and achieves nothing.
I'll take that as a yes. What needs to be changed if I want to compute a bunch of Fibonacci numbers, say 1000, in a loop? Or use GMP in an extension to a long running program written in Script Basic?
There must be an mpz_free function, or your prior suggestion of making the local's static might be the answer.
mpz_t is tiny of course but it has pointers to the allocated memory.

Heater
Posts: 13084
Joined: Tue Jul 17, 2012 3:02 pm

Re: A Final Fibonacci Challenge

Tue Jun 11, 2019 5:53 pm

jahboater,
When a process exits all allocated memory is discarded and reclaimed by the OS; the address space dissapears. True for all operating systems except Novel's Netware.
Few programs bother to release memory before calling exit as may add considerably to the complexity and achieves nothing.
It's true that people often leave the OS to clean up their mess.

I have learned that it can lead to a world of pain in the future. And as such is very bad practice.

This is not acceptable.

Heater
Posts: 13084
Joined: Tue Jul 17, 2012 3:02 pm

Re: A Final Fibonacci Challenge

Tue Jun 11, 2019 5:56 pm

jahboater,
There must be an mpz_free function,...
According to the docs GMP uses normal malloc by default. So it should make no odds to use free.

One can change the allocator GMP uses and manage things ones self. I have not done that.

Heater
Posts: 13084
Joined: Tue Jul 17, 2012 3:02 pm

Re: A Final Fibonacci Challenge

Tue Jun 11, 2019 6:06 pm

How about we straighten out the fibo() to remove memory leaks. Like so:

Code: Select all

char* fibois (int n) {
    char* res;
    if (n <= 2) {
        return letis(memo[n]);
    }

    int k = (n / 2);
    char* fk = fibois(k);
    char* fk1 = fibois(k + 1);
    char* a;
    char* b;
    if ((n % 2) == 0) {
        a = mulis(fk1, "2");
        b = subis(a, fk);
        res = mulis(fk, b);
    } else {
        a = mulis(fk, fk);
        b = mulis(fk1, fk1);
        res = addis(a, b);
    }
    free(a);
    free(b);
    free(fk);
    free(fk1);
    return res;
}

jahboater
Posts: 4598
Joined: Wed Feb 04, 2015 6:38 pm

Re: A Final Fibonacci Challenge

Tue Jun 11, 2019 6:09 pm

Heater wrote:
Tue Jun 11, 2019 5:53 pm
I have learned that it can lead to a world of pain in the future. And as such is very bad practice.
I am only talking about just prior to a programs exit. On a proper operating system the entire address space of the process vanishes when you execute the exit system call.
How do any changes to the heap memory within that address space just prior to its disappearance help avoid problems in the future?
I would be genuinely interested to know.

(Of course it may be a problem with a weak OS like Netware, or possibly MSDOS I cant remember, and for debugging/verification, it might be useful to be able to prove all the allocated memory is accessible and accounted for at the end).

jahboater
Posts: 4598
Joined: Wed Feb 04, 2015 6:38 pm

Re: A Final Fibonacci Challenge

Tue Jun 11, 2019 6:14 pm

Heater wrote:
Tue Jun 11, 2019 6:06 pm
How about we straighten out the fibo() to remove memory leaks. Like so:
But what about the mulis() and addis() functions?

mpz_t is tiny but has pointers to the main array's. I thought that's where the leak was?
Declaring them static might help so they can be re-used.

I don't know why valgrind didn't complain - unless GMP has some global references to the allocated memory that kept it happy.
Last edited by jahboater on Tue Jun 11, 2019 6:15 pm, edited 1 time in total.

ejolson
Posts: 3406
Joined: Tue Mar 18, 2014 11:47 am

Re: A Final Fibonacci Challenge

Tue Jun 11, 2019 6:15 pm

jahboater wrote:
Tue Jun 11, 2019 6:09 pm
Heater wrote:
Tue Jun 11, 2019 5:53 pm
I have learned that it can lead to a world of pain in the future. And as such is very bad practice.
I am only talking about just prior to a programs exit. On a proper operating system the entire address space of the process vanishes when you execute the exit system call.
On the 8-bit personal computers of the golden age one can achieve the same effect by flipping the power switch off and on. Turning it off and on again is a popular technique even unto this day.
Last edited by ejolson on Tue Jun 11, 2019 6:42 pm, edited 1 time in total.

jalih
Posts: 75
Joined: Mon Apr 15, 2019 3:54 pm

Re: A Final Fibonacci Challenge

Tue Jun 11, 2019 6:29 pm

ejolson wrote:
Sat Jun 08, 2019 5:53 pm
From what I understand new PL/I and C codes are still under preparation. Is there anything else?
Here is PL/I conversion from your visual.bas source:

Code: Select all


*PROCESS MARGINS(1,120) LIBS(SINGLE,STATIC);
*PROCESS OPTIMIZE(2) DFT(REORDER);
*PROCESS LIMITS(FIXEDBIN(63));


/*
    Almost direct PL/I conversion from:

    visual.bas -- Compute the nth Fibonacci Number
    Written January 14, 2018 by Eric Olson
                                            */

fibonacci: package;

 dcl (cmult, rexp, nmax) fixed bin(31);
 dcl radix fixed bin(32) unsigned;
 dcl (cadd, cmax) fixed bin(64) unsigned;

 dcl (a(0:nmax), b(0:nmax), t1(0:nmax), t2(0:nmax)) fixed bin(32) unsigned ctl;
 dcl xy(0:nmax) fixed bin(64) unsigned ctl;


 bigprint: proc(x);
   dcl x(*) fixed bin(32) unsigned;

   dcl i fixed bin(31);
   if x(0) = 0 then
       put list('0') skip;
   else
     do;
       do i = x(0) to 1 by -1;
         dcl s0 pic'ZZZZZZZZZ';
         dcl s pic'999999999';
         dcl j fixed bin(31);
         j = x(i);
         if i = x(0) then
           do;
             s0 = j;
             put edit(trim(s0)) (A);
           end;
         else
           do;
             s = j;
             put edit(s) (A);
           end;
       end;
       put skip;
     end;
 end bigprint;


 biginc: proc(z);
   dcl z(*) fixed bin(32) unsigned;

   dcl i fixed bin(31);
   dcl (t, c) fixed bin(32) unsigned;

   i, c = 1;
   do while(c > 0);
     t = z(i) + 1;
     if t >= radix then
       z(i) = t - radix;
     else
       do;
         z(i) = t;
         c = 0;
       end;
     i += 1;
   end;
 end;


 bigadd: proc(z, x, y);
   dcl (z(*), x(*), y(*)) fixed bin(32) unsigned;

   dcl (i, max) fixed bin(31);

   if y(0) > x(0) then
     max = y(0);
   else
     max = x(0);

   do i = 1 to max + 1;
     z(i) = 0;
   end;

   do i = 1 to max;
     dcl (t, c) fixed bin(31);
     c = 0;
     if i <= x(0) then
       do;
         t = z(i) + x(i);
         if i <= y(0) then
           t += y(i);
       end;
     else
       t = z(i) + y(i);

     if t >= radix then
       do;
         c = 1;
         t -= radix;
       end;
     z(i) = t;
     z(i+1) += c;
   end;
   z(0)=max+1;
   do while(z(z(0)) = 0 & z(0) > 1);
     z(0) = z(0) - 1;
   end;
   return;
 end;


 bigdec: proc(z);
   dcl z(*) fixed bin(32) unsigned;

   dcl i fixed bin(31);
   dcl c fixed bin(32) unsigned;

   i, c = 1;
   do while(c > 0);
     if z(i) < 1 then
       z(i) = radix - 1;
     else
       do;
         z(i) -= 1;
         c = 0;
       end;
     i += 1;
   end;
 end;


 bigsub: proc(z, x, y);
   dcl (z(*), x(*), y(*)) fixed bin(32) unsigned;

   dcl i fixed bin(31);
   dcl (t, c) fixed bin(32) unsigned;

   c = 0;
   do i = 1 to y(0);
     t = y(i) + c;
     if x(i) < t then
       do;
         z(i) = x(i) + radix - t;
         c = 1;
       end;
     else
       do;
         z(i) = x(i) - t;
         c = 0;
       end;
   end;
   i = y(0) + 1;
   do while(c > 0);
     if x(i) < 1 then
       z(i) = x(i) + radix - 1;
     else
       do;
         z(i) = x(i) - 1;
         c = 0;
       end;
     i += 1;
   end;
   do while(i <= x(0));
     z(i) = x(i);
     i += 1;
   end;
   z(0) = x(0);
   do while(z(z(0)) = 0 & z(0) > 1);
     z(0) -= 1;
   end;
   return;
 end;



 bigmulw: proc(z, x, y);
   dcl (z(*), x(*), y(*)) fixed bin(32) unsigned;

   dcl (i, j, imax) fixed bin(31);

   imax = x(0) + y(0);
   do i = 1 to imax;
     xy(i) = 0;
   end;

   do i = 1 to x(0);
     do j = 1 to y(0);
       dcl (xt, yt) fixed bin(64) unsigned;
       xt = x(i);
       yt = y(j);
       xy(i+j-1) = xy(i+j-1) + xt * yt;
     end;
     if mod(i,cmult) = 0 then
       do;
         dcl (k, kmax) fixed bin(31);
         kmax = i + y(0) - 1;
         do k = 1 to kmax;
           if xy(k) >= cmax then
             do;
               xy(k) = xy(k) - cmax;
               xy(k+1) = xy(k+1) + cadd;
             end;
         end;
       end;
   end;
   do i = 1 to imax - 1;
     xy(i+1)=xy(i+1)+xy(i)/radix;
     z(i) = mod(xy(i), radix);
   end;
   z(imax) = xy(imax);
   do while(z(imax) = 0 & imax > 1);
     imax -= 1;
   end;
   z(0) = imax;
   return;
 end;


 bighigh: proc(z, x, n);
   dcl (z(*), x(*)) fixed bin(32) unsigned;
   dcl n fixed bin(31);

   if x(0) <= n then
     do;
       z(0) = 0;
       return;
     end;
   dcl i fixed bin(31);
   do i = n + 1 to x(0);
     z(i-n) = x(i);
   end;
   z(0) = x(0) - n;
   return;
 end bighigh;


 biglow: proc(z, x, n);
   dcl (z(*), x(*)) fixed bin(32) unsigned;
   dcl n fixed bin(31);

   if x(0) <= n then
     z(0) = x(0);
   else
     z(0) = n;

   dcl i fixed bin(31);

   do i = 1 to z(0);
     z(i) = x(i);
   end;
   return;
 end biglow;


 bigmul: proc(z, x, y) recursive;
   dcl (z(*), x(*), y(*)) fixed bin(32) unsigned;

   if x(0) < 44 | y(0) < 44 then
     do;
       call bigmulw(z,x,y);
       return;
     end;

   dcl (i, M, n) fixed bin(31);

   if x(0) < y(0) then
     M = y(0);
   else
     M = x(0);

   n = M / 2;

   dcl (xl(0:n+1), xh(0:n+1), yl(0:n+1), yh(0:n+1)) fixed bin(32) unsigned ctl;
   allocate xl, xh, yl, yh;
   call biglow(xl, x, n);
   call bighigh(xh, x, n);
   call biglow(yl, y, n);
   call bighigh(yh, y, n);

   dcl (z0(0:M+1), z1(0:M+3), z2(0:M+3)) fixed bin(32) unsigned ctl;
   allocate z0, z1, z2;
   call bigadd(z0,xl,xh);
   call bigadd(z2,yl,yh);
   call bigmul(z1,z0,z2);
   call bigmul(z0,xl,yl);
   call bigmul(z2,xh,yh);
   call bigsub(z1,z1,z0);
   call bigsub(z1,z1,z2);

   dcl imax fixed bin(31);
   imax = x(0) + y(0);
   dcl (t, c) fixed bin(32) unsigned;
   c = 0;
   do i = 1 to imax;
     t=c;
     if i <= z0(0) then t += z0(i);
     if i > n & i-n <= z1(0) then t += z1(i-n);
     if i > 2*n & i - 2 * n <= z2(0) then t += z2(i-2*n);
     c = 0;
     do while(t >= radix);
       t -= radix;
       c += 1;
     end;
     z(i) = t;
   end;
   do while(z(imax) = 0 & imax > 1);
     imax -= 1;
   end;
   z(0) = imax;
   free xl, xh, yl, yh, z0, z1, z2;
   return;
 end bigmul;


 fibowork: proc(n) recursive;
   dcl n fixed bin(31);

   if  n = 0 then
     do;
       a(0) = 0;
       b(0) = 1;
       b(1) = 1;
       return;
     end;
   call fibowork(n/2);
   if mod(n,2) = 0 then
     do;
       call bigadd(t1,b,b);
       call bigsub(t2,t1,a);
       call bigmul(a,a,t2);
       call bigmul(b,b,t2);
       if mod(n,4) = 0 then
         call bigdec(b);
       else
         call biginc(b);
     end;
   else
     do;
       call bigadd(t1,a,a);
       call bigadd(t2,t1,b);
       call bigmul(b,b,t2);
       call bigmul(a,a,t2);
       if mod(n,4) = 1 then
         call biginc(a);
       else
         call bigdec(a);
     end;
   return;
 end fibowork;


fibo: proc(n);
  dcl n fixed bin(31);

  if n < 2 then
    do;
      b(0) = 1;
      b(1) = n;
      return;
    end;
  call fibowork((n-1)/2);
  if mod(n,2) = 0 then
    do;
      call bigadd(t1,a,a);
      call bigadd(t2,t1,b);
      call bigmul(b,b,t2);
    end;
  else
    do;
      call bigadd(t1,b,b);
      call bigsub(t2,t1,a);
      call bigmul(b,b,t2);
      if mod(n,4) = 1 then
        call bigdec(b);
      else
        call biginc(b);
    end;
    return;
 end fibo;


 test: proc options(main);
   dcl n fixed bin(31) init(4784969);
   rexp = 9;
   cmult = 5;
   radix = 10**rexp;
   cadd = radix;
   cadd *= cmult;
   cmax = cadd * radix;
   nmax = 205100;
   allocate a, b, t1, t2, xy;

   call fibo(n);
   call bigprint(b);
   free a, b, t1, t2, xy;
 end test;

end fibonacci;
It takes about 6.5 seconds to run on my PC. Sadly, you cannot test this on Pi as IronSpring PL/I compiler is the only freely available PL/I compiler and not available for Pi. I used old IBM VisualAge PL/I compiler for testing verifying it works correctly.

jahboater
Posts: 4598
Joined: Wed Feb 04, 2015 6:38 pm

Re: A Final Fibonacci Challenge

Tue Jun 11, 2019 6:32 pm

ejolson wrote:
Tue Jun 11, 2019 5:46 pm
What needs to be changed if I want to compute a bunch of Fibonacci numbers, say 1000, in a loop? Or use GMP in an extension to a long running program written in Script Basic?
You should call mpz_clear() which frees the memory.

jahboater
Posts: 4598
Joined: Wed Feb 04, 2015 6:38 pm

Re: A Final Fibonacci Challenge

Tue Jun 11, 2019 6:35 pm

Heater wrote:
Tue Jun 11, 2019 5:56 pm
jahboater,
There must be an mpz_free function,...
According to the docs GMP uses normal malloc by default. So it should make no odds to use free.

One can change the allocator GMP uses and manage things ones self. I have not done that.
Its called mpz_clear(). If you were to use free() it would mean looking at the internals of mpz_t.

Code: Select all

Function: void mpz_clear (mpz_t x)

    Free the space occupied by x. Call this function for all mpz_t variables when you are done with them. 

Function: void mpz_clears (mpz_t x, ...)

    Free the space occupied by a NULL-terminated list of mpz_t variables. 


Heater
Posts: 13084
Joined: Tue Jul 17, 2012 3:02 pm

Re: A Final Fibonacci Challenge

Tue Jun 11, 2019 6:36 pm

jahboater,
You should call mpz_clear() which frees the memory.
It's done already whist you were typing :)

No leaks found by valgrind with this version:

Code: Select all

//
// An experiment in doing integer arithmetic using GMP with all numbers represented by strings.
//
// By heater.
// Modified June 11, 2019 to use base 32 strings for intermediate results.
//
#include <stdio.h>
#include <gmp.h>
#include <stdlib.h>
#include <memory.h>
#include <string.h>

// Number base used for internal calculations by GMP.
int BASE = 32;

// Functions letis, addis, subis and mulis do large integer arithmetic on integers represented by strings.

void writeis(const char *s) {
    mpz_t op1;
    mpz_init(op1);
    mpz_set_str (op1, s, BASE);
    char *buf=mpz_get_str (0, 10, op1);
    mpz_clear(op1);
    puts(buf);
    free(buf);
}

char* letis(const char* s) {
//    size_t size = strlen(s) + 1;
//    char* res = malloc(size);
//    strncpy (res, s, size);
    return strdup(s);
}

char* addis(const char* s1, const char* s2) {
    mpz_t op1;
    mpz_t op2;
    mpz_t res;
    mpz_init(op1);
    mpz_init(op2);
    mpz_init(res);

    mpz_set_str (op1, s1, BASE);
    mpz_set_str (op2, s2, BASE);
    mpz_add (res, op1, op2);  // result = x * y
    mpz_clear(op1);
    mpz_clear(op2);
    char* res_string = mpz_get_str (0, BASE, res); 
    mpz_clear(res);
    return res_string;
}

char* subis(const char* s1, const char* s2) {
    mpz_t op1;
    mpz_t op2;
    mpz_t res;
    mpz_init(op1);
    mpz_init(op2);
    mpz_init(res);

    mpz_set_str (op1, s1, BASE);
    mpz_set_str (op2, s2, BASE);
    mpz_sub (res, op1, op2);  // result = x * y
    mpz_clear(op1);
    mpz_clear(op2);
    char* res_string = mpz_get_str (0, BASE, res); 
    mpz_clear(res);
    return res_string;
}

char* mulis(const char* s1, const char* s2) {
    mpz_t op1;
    mpz_t op2;
    mpz_t res;
    mpz_init(op1);
    mpz_init(op2);
    mpz_init(res);

    mpz_set_str (op1, s1, BASE);
    mpz_set_str (op2, s2, BASE);
    mpz_mul (res, op1, op2);  // result = x * y
    mpz_clear(op1);
    mpz_clear(op2);
    char* res_string = mpz_get_str (0, BASE, res); 
    mpz_clear(res);
    return res_string;
}

char* memo[3];

void init_memo() {
    memo[0] = letis("0");
    memo[1] = letis("1");
    memo[2] = letis("1");
}

void clean_memo() {
    free(memo[0]);
    free(memo[1]);
    free(memo[2]);
}


// Return the n'th Fibonacci number as a decimal string for integer n
char* fibois (int n) {
    char* res;
    if (n <= 2) {
        return letis(memo[n]);
    }

    int k = (n / 2);
    char* fk = fibois(k);
    char* fk1 = fibois(k + 1);
    char* a;
    char* b;
    if ((n % 2) == 0) {
        a = addis(fk1, fk1);
        b = subis(a, fk);
        res = mulis(fk, b);
    } else {
        a = mulis(fk, fk);
        b = mulis(fk1, fk1);
        res = addis(a, b);
    }
    free(a);
    free(b);
    free(fk);
    free(fk1);
    return res;
}

int main(int argc, char* argv[]) {
    char* f;
    int n = 4784969;               // The first Fibonacci number with a million digits

    if (argc >= 2) {
        n = atoi(argv[1]);
    }

    init_memo();

    f = fibois(n);
    writeis(f);
    free(f);

    clean_memo();

    return (0);
}

Code: Select all

$ valgrind --leak-check=yes ./fibo_strings 100
==26619== Memcheck, a memory error detector
==26619== Copyright (C) 2002-2015, and GNU GPL'd, by Julian Seward et al.
==26619== Using Valgrind-3.12.0.SVN and LibVEX; rerun with -h for copyright info
==26619== Command: ./fibo_strings 100
==26619==
==26619== error calling PR_SET_PTRACER, vgdb might block
354224848179261915075
==26619==
==26619== HEAP SUMMARY:
==26619==     in use at exit: 0 bytes in 0 blocks
==26619==   total heap usage: 1,790 allocs, 1,790 frees, 18,818 bytes allocated
==26619==
==26619== All heap blocks were freed -- no leaks are possible
==26619==
==26619== For counts of detected and suppressed errors, rerun with: -v
==26619== ERROR SUMMARY: 0 errors from 0 contexts (suppressed: 0 from 0)

Heater
Posts: 13084
Joined: Tue Jul 17, 2012 3:02 pm

Re: A Final Fibonacci Challenge

Tue Jun 11, 2019 6:41 pm

As a bonus it runs twice as fast as it did this morning!

Code: Select all

$ time ./fibo_strings | head -c 32 ;  time ./fibo_strings | tail  -c 32
10727395641800477229364813596225
real    0m4.210s
user    0m4.141s
sys     0m0.047s
4856539211500699706378405156269

real    0m4.085s
user    0m4.000s
sys     0m0.047s
Code is in the repo.

Thanks everyone.

jahboater
Posts: 4598
Joined: Wed Feb 04, 2015 6:38 pm

Re: A Final Fibonacci Challenge

Tue Jun 11, 2019 6:52 pm

Heater wrote:
Tue Jun 11, 2019 6:41 pm
As a bonus it runs twice as fast as it did this morning!
Cool!

What about Ejolson's idea of making them static?

Code: Select all

static char* addis(const char* s1, const char* s2) {
    static mpz_t op1, op2, res;
    static bool first = true;
          
    if( first )
    {
      mpz_init(op1);
      mpz_init(op2);
      mpz_init(res);
      first = false;
    }

    mpz_set_str (op1, s1, BASE);
    mpz_set_str (op2, s2, BASE);
    mpz_add (res, op1, op2);  // result = x * y

    return mpz_get_str (0, BASE, res);
}
However, I bet the cost of converting to/from strings all time is the performance bottleneck.

Heater
Posts: 13084
Joined: Tue Jul 17, 2012 3:02 pm

Re: A Final Fibonacci Challenge

Tue Jun 11, 2019 7:11 pm

jahboater,
What about Ejolson's idea of making them static?
Oh my. Then it's 4 times faster than it was this morning!

Code: Select all

$ time ./fibo_strings | head -c 32 ;  time ./fibo_strings | tail  -c 32
10727395641800477229364813596225
real    0m2.009s
user    0m1.938s
sys     0m0.047s
4856539211500699706378405156269

real    0m2.022s
user    0m1.969s
sys     0m0.047s
This worries me as it's not thread safe anymore.

I just made op1, op2 and res global an initialized them in main. Like so:

Code: Select all

//
// An experiment in doing integer arithmetic using GMP with all numbers represented by strings.
//
// By heater.
// Modified June 11, 2019 to use base 32 strings for intermediate results.
//
#include <stdio.h>
#include <gmp.h>
#include <stdlib.h>
#include <memory.h>
#include <string.h>

// Number base used for internal calculations by GMP.
int BASE = 32;

mpz_t op1;
mpz_t op2;
mpz_t res;

// Functions letis, addis, subis and mulis do large integer arithmetic on integers represented by strings.

void writeis(const char *s) {
    mpz_set_str (op1, s, BASE);
    char *buf=mpz_get_str (0, 10, op1);
    puts(buf);
    free(buf);
}

char* letis(const char* s) {
    return strdup(s);
}

char* addis(const char* s1, const char* s2) {
    mpz_set_str (op1, s1, BASE);
    mpz_set_str (op2, s2, BASE);
    mpz_add (res, op1, op2);  // result = x * y
    char* res_string = mpz_get_str (0, BASE, res); 
    return res_string;
}

char* subis(const char* s1, const char* s2) {
    mpz_set_str (op1, s1, BASE);
    mpz_set_str (op2, s2, BASE);
    mpz_sub (res, op1, op2);  // result = x * y
    char* res_string = mpz_get_str (0, BASE, res); 
    return res_string;
}

char* mulis(const char* s1, const char* s2) {
    mpz_set_str (op1, s1, BASE);
    mpz_set_str (op2, s2, BASE);
    mpz_mul (res, op1, op2);  // result = x * y
    char* res_string = mpz_get_str (0, BASE, res); 
    return res_string;
}

char* memo[3];

void init_memo() {
    memo[0] = letis("0");
    memo[1] = letis("1");
    memo[2] = letis("1");
}

void clean_memo() {
    free(memo[0]);
    free(memo[1]);
    free(memo[2]);
}


// Return the n'th Fibonacci number as a decimal string for integer n
char* fibois (int n) {
    char* res;
    if (n <= 2) {
        return letis(memo[n]);
    }

    int k = (n / 2);
    char* fk = fibois(k);
    char* fk1 = fibois(k + 1);
    char* a;
    char* b;
    if ((n % 2) == 0) {
        a = addis(fk1, fk1);
        b = subis(a, fk);
        res = mulis(fk, b);
    } else {
        a = mulis(fk, fk);
        b = mulis(fk1, fk1);
        res = addis(a, b);
    }
    free(a);
    free(b);
    free(fk);
    free(fk1);
    return res;
}

int main(int argc, char* argv[]) {
    char* f;
    int n = 4784969;               // The first Fibonacci number with a million digits

    if (argc >= 2) {
        n = atoi(argv[1]);
    }

    mpz_init(op1);
    mpz_init(op2);
    mpz_init(res);

    init_memo();

    f = fibois(n);
    writeis(f);
    free(f);

    clean_memo();

    mpz_clear(op1);
    mpz_clear(op2);
    mpz_clear(res);

    return (0);
}
Last edited by Heater on Tue Jun 11, 2019 7:59 pm, edited 1 time in total.

User avatar
John_Spikowski
Posts: 1329
Joined: Wed Apr 03, 2019 5:53 pm
Location: Anacortes, WA USA
Contact: Website Twitter

Re: A Final Fibonacci Challenge

Tue Jun 11, 2019 7:20 pm

Is it safe to take another shot at getting this code converted to a ScriptBasic extension module?

I'm going to have to use the new GMP library to calculate the percent improvement from the last attempt at this. :D
Last edited by John_Spikowski on Tue Jun 11, 2019 7:32 pm, edited 2 times in total.

jahboater
Posts: 4598
Joined: Wed Feb 04, 2015 6:38 pm

Re: A Final Fibonacci Challenge

Tue Jun 11, 2019 7:25 pm

Heater wrote:
Tue Jun 11, 2019 7:11 pm
I just made op1, op2 and res global and initialized them in main. Like so:
That's a much better idea, I hate controlling logic with booleans.

I'd declare them static too (or use -fwhole-program which will declare everything static) which might shave a few insns.

ejolson
Posts: 3406
Joined: Tue Mar 18, 2014 11:47 am

Re: A Final Fibonacci Challenge

Tue Jun 11, 2019 7:28 pm

jalih wrote:
Tue Jun 11, 2019 6:29 pm
ejolson wrote:
Sat Jun 08, 2019 5:53 pm
From what I understand new PL/I and C codes are still under preparation. Is there anything else?
Here is PL/I conversion from your visual.bas source:

Code: Select all


*PROCESS MARGINS(1,120) LIBS(SINGLE,STATIC);
*PROCESS OPTIMIZE(2) DFT(REORDER);
*PROCESS LIMITS(FIXEDBIN(63));


/*
    Almost direct PL/I conversion from:

    visual.bas -- Compute the nth Fibonacci Number
    Written January 14, 2018 by Eric Olson
                                            */

fibonacci: package;

 dcl (cmult, rexp, nmax) fixed bin(31);
 dcl radix fixed bin(32) unsigned;
 dcl (cadd, cmax) fixed bin(64) unsigned;

 dcl (a(0:nmax), b(0:nmax), t1(0:nmax), t2(0:nmax)) fixed bin(32) unsigned ctl;
 dcl xy(0:nmax) fixed bin(64) unsigned ctl;


 bigprint: proc(x);
   dcl x(*) fixed bin(32) unsigned;

   dcl i fixed bin(31);
   if x(0) = 0 then
       put list('0') skip;
   else
     do;
       do i = x(0) to 1 by -1;
         dcl s0 pic'ZZZZZZZZZ';
         dcl s pic'999999999';
         dcl j fixed bin(31);
         j = x(i);
         if i = x(0) then
           do;
             s0 = j;
             put edit(trim(s0)) (A);
           end;
         else
           do;
             s = j;
             put edit(s) (A);
           end;
       end;
       put skip;
     end;
 end bigprint;


 biginc: proc(z);
   dcl z(*) fixed bin(32) unsigned;

   dcl i fixed bin(31);
   dcl (t, c) fixed bin(32) unsigned;

   i, c = 1;
   do while(c > 0);
     t = z(i) + 1;
     if t >= radix then
       z(i) = t - radix;
     else
       do;
         z(i) = t;
         c = 0;
       end;
     i += 1;
   end;
 end;


 bigadd: proc(z, x, y);
   dcl (z(*), x(*), y(*)) fixed bin(32) unsigned;

   dcl (i, max) fixed bin(31);

   if y(0) > x(0) then
     max = y(0);
   else
     max = x(0);

   do i = 1 to max + 1;
     z(i) = 0;
   end;

   do i = 1 to max;
     dcl (t, c) fixed bin(31);
     c = 0;
     if i <= x(0) then
       do;
         t = z(i) + x(i);
         if i <= y(0) then
           t += y(i);
       end;
     else
       t = z(i) + y(i);

     if t >= radix then
       do;
         c = 1;
         t -= radix;
       end;
     z(i) = t;
     z(i+1) += c;
   end;
   z(0)=max+1;
   do while(z(z(0)) = 0 & z(0) > 1);
     z(0) = z(0) - 1;
   end;
   return;
 end;


 bigdec: proc(z);
   dcl z(*) fixed bin(32) unsigned;

   dcl i fixed bin(31);
   dcl c fixed bin(32) unsigned;

   i, c = 1;
   do while(c > 0);
     if z(i) < 1 then
       z(i) = radix - 1;
     else
       do;
         z(i) -= 1;
         c = 0;
       end;
     i += 1;
   end;
 end;


 bigsub: proc(z, x, y);
   dcl (z(*), x(*), y(*)) fixed bin(32) unsigned;

   dcl i fixed bin(31);
   dcl (t, c) fixed bin(32) unsigned;

   c = 0;
   do i = 1 to y(0);
     t = y(i) + c;
     if x(i) < t then
       do;
         z(i) = x(i) + radix - t;
         c = 1;
       end;
     else
       do;
         z(i) = x(i) - t;
         c = 0;
       end;
   end;
   i = y(0) + 1;
   do while(c > 0);
     if x(i) < 1 then
       z(i) = x(i) + radix - 1;
     else
       do;
         z(i) = x(i) - 1;
         c = 0;
       end;
     i += 1;
   end;
   do while(i <= x(0));
     z(i) = x(i);
     i += 1;
   end;
   z(0) = x(0);
   do while(z(z(0)) = 0 & z(0) > 1);
     z(0) -= 1;
   end;
   return;
 end;



 bigmulw: proc(z, x, y);
   dcl (z(*), x(*), y(*)) fixed bin(32) unsigned;

   dcl (i, j, imax) fixed bin(31);

   imax = x(0) + y(0);
   do i = 1 to imax;
     xy(i) = 0;
   end;

   do i = 1 to x(0);
     do j = 1 to y(0);
       dcl (xt, yt) fixed bin(64) unsigned;
       xt = x(i);
       yt = y(j);
       xy(i+j-1) = xy(i+j-1) + xt * yt;
     end;
     if mod(i,cmult) = 0 then
       do;
         dcl (k, kmax) fixed bin(31);
         kmax = i + y(0) - 1;
         do k = 1 to kmax;
           if xy(k) >= cmax then
             do;
               xy(k) = xy(k) - cmax;
               xy(k+1) = xy(k+1) + cadd;
             end;
         end;
       end;
   end;
   do i = 1 to imax - 1;
     xy(i+1)=xy(i+1)+xy(i)/radix;
     z(i) = mod(xy(i), radix);
   end;
   z(imax) = xy(imax);
   do while(z(imax) = 0 & imax > 1);
     imax -= 1;
   end;
   z(0) = imax;
   return;
 end;


 bighigh: proc(z, x, n);
   dcl (z(*), x(*)) fixed bin(32) unsigned;
   dcl n fixed bin(31);

   if x(0) <= n then
     do;
       z(0) = 0;
       return;
     end;
   dcl i fixed bin(31);
   do i = n + 1 to x(0);
     z(i-n) = x(i);
   end;
   z(0) = x(0) - n;
   return;
 end bighigh;


 biglow: proc(z, x, n);
   dcl (z(*), x(*)) fixed bin(32) unsigned;
   dcl n fixed bin(31);

   if x(0) <= n then
     z(0) = x(0);
   else
     z(0) = n;

   dcl i fixed bin(31);

   do i = 1 to z(0);
     z(i) = x(i);
   end;
   return;
 end biglow;


 bigmul: proc(z, x, y) recursive;
   dcl (z(*), x(*), y(*)) fixed bin(32) unsigned;

   if x(0) < 44 | y(0) < 44 then
     do;
       call bigmulw(z,x,y);
       return;
     end;

   dcl (i, M, n) fixed bin(31);

   if x(0) < y(0) then
     M = y(0);
   else
     M = x(0);

   n = M / 2;

   dcl (xl(0:n+1), xh(0:n+1), yl(0:n+1), yh(0:n+1)) fixed bin(32) unsigned ctl;
   allocate xl, xh, yl, yh;
   call biglow(xl, x, n);
   call bighigh(xh, x, n);
   call biglow(yl, y, n);
   call bighigh(yh, y, n);

   dcl (z0(0:M+1), z1(0:M+3), z2(0:M+3)) fixed bin(32) unsigned ctl;
   allocate z0, z1, z2;
   call bigadd(z0,xl,xh);
   call bigadd(z2,yl,yh);
   call bigmul(z1,z0,z2);
   call bigmul(z0,xl,yl);
   call bigmul(z2,xh,yh);
   call bigsub(z1,z1,z0);
   call bigsub(z1,z1,z2);

   dcl imax fixed bin(31);
   imax = x(0) + y(0);
   dcl (t, c) fixed bin(32) unsigned;
   c = 0;
   do i = 1 to imax;
     t=c;
     if i <= z0(0) then t += z0(i);
     if i > n & i-n <= z1(0) then t += z1(i-n);
     if i > 2*n & i - 2 * n <= z2(0) then t += z2(i-2*n);
     c = 0;
     do while(t >= radix);
       t -= radix;
       c += 1;
     end;
     z(i) = t;
   end;
   do while(z(imax) = 0 & imax > 1);
     imax -= 1;
   end;
   z(0) = imax;
   free xl, xh, yl, yh, z0, z1, z2;
   return;
 end bigmul;


 fibowork: proc(n) recursive;
   dcl n fixed bin(31);

   if  n = 0 then
     do;
       a(0) = 0;
       b(0) = 1;
       b(1) = 1;
       return;
     end;
   call fibowork(n/2);
   if mod(n,2) = 0 then
     do;
       call bigadd(t1,b,b);
       call bigsub(t2,t1,a);
       call bigmul(a,a,t2);
       call bigmul(b,b,t2);
       if mod(n,4) = 0 then
         call bigdec(b);
       else
         call biginc(b);
     end;
   else
     do;
       call bigadd(t1,a,a);
       call bigadd(t2,t1,b);
       call bigmul(b,b,t2);
       call bigmul(a,a,t2);
       if mod(n,4) = 1 then
         call biginc(a);
       else
         call bigdec(a);
     end;
   return;
 end fibowork;


fibo: proc(n);
  dcl n fixed bin(31);

  if n < 2 then
    do;
      b(0) = 1;
      b(1) = n;
      return;
    end;
  call fibowork((n-1)/2);
  if mod(n,2) = 0 then
    do;
      call bigadd(t1,a,a);
      call bigadd(t2,t1,b);
      call bigmul(b,b,t2);
    end;
  else
    do;
      call bigadd(t1,b,b);
      call bigsub(t2,t1,a);
      call bigmul(b,b,t2);
      if mod(n,4) = 1 then
        call bigdec(b);
      else
        call biginc(b);
    end;
    return;
 end fibo;


 test: proc options(main);
   dcl n fixed bin(31) init(4784969);
   rexp = 9;
   cmult = 5;
   radix = 10**rexp;
   cadd = radix;
   cadd *= cmult;
   cmax = cadd * radix;
   nmax = 205100;
   allocate a, b, t1, t2, xy;

   call fibo(n);
   call bigprint(b);
   free a, b, t1, t2, xy;
 end test;

end fibonacci;
It takes about 6.5 seconds to run on my PC. Sadly, you cannot test this on Pi as IronSpring PL/I compiler is the only freely available PL/I compiler and not available for Pi. I used old IBM VisualAge PL/I compiler for testing verifying it works correctly.
Very nice! I wonder if it will run under QEMU on the Pi.

Unfortunately, it may be some time before I have an opportunity to try.

Heater
Posts: 13084
Joined: Tue Jul 17, 2012 3:02 pm

Re: A Final Fibonacci Challenge

Tue Jun 11, 2019 7:33 pm

ScriptBasic,
Is it safe to take another shot at getting this code converted to a ScriptBasic extension module?
I think so. It runs. It gives the right result. It shows no sign of leaking memory according to Valgrind.

There are still errors when trying to use the address and memory sanitizers. But that seems to be inside mpz_set_str and I have no idea what to do about that yet.

Be warned that the last version posted is not thread safe. It is using global op1, op2 and res.

User avatar
John_Spikowski
Posts: 1329
Joined: Wed Apr 03, 2019 5:53 pm
Location: Anacortes, WA USA
Contact: Website Twitter

Re: A Final Fibonacci Challenge

Tue Jun 11, 2019 10:47 pm

I'm still getting seg faults trying to run the code. Can the C pros on this thread have a peek and help me figure out what the problem is or the best way to debug it?

It gets to the first GMP::BI_ADD() function and seg faults. The original fibo function AIR helped out with works. My guess is SB doesn't like how Heater wants to return strings.

interface.c

Code: Select all

/* GMP Extension Module
UXLIBS: -lgmp
*/

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <unistd.h>
#include <memory.h>
#include <gmp.h>
#include "../../basext.h"


int BASE = 32;
mpz_t op1;
mpz_t op2;
mpz_t res;
char* memo[3];


static char* let(const char* s) {
    return strdup(s);
}


/**************************
 Extension Module Functions
**************************/

typedef struct _ModuleObject {
  void *HandleArray;
}ModuleObject,*pModuleObject;


besVERSION_NEGOTIATE
  return (int)INTERFACE_VERSION;
besEND


besSUB_START
  pModuleObject p;

  besMODULEPOINTER = besALLOC(sizeof(ModuleObject));
  if( besMODULEPOINTER == NULL )return 0;

  p = (pModuleObject)besMODULEPOINTER;
  return 0;
besEND


besSUB_FINISH
  pModuleObject p;

  p = (pModuleObject)besMODULEPOINTER;
  if( p == NULL )return 0;
  return 0;
besEND


/*************
 GMP Functions
*************/


besFUNCTION(fibo)
  int fval;

  besARGUMENTS("i")
    &fval
  besARGEND

  char buf[fval+1];
  memset(buf,0,1);
  mpz_t res;
  mpz_init(res);

  mpz_fib_ui(res, fval);

  gmp_snprintf( buf,sizeof(buf),"%Zd", res );

  besRETURN_STRING(buf);
besEND


besFUNCTION(writeis)
  const char *s;

  besARGUMENTS("z")
    &s
  besARGEND

  mpz_set_str (op1, s, BASE);
  char *buf=mpz_get_str (0, 10, op1);
  puts(buf);
  free(buf);
  besRETURNVALUE = NULL;

besEND


besFUNCTION(letis)
  const char* s;

  besARGUMENTS("z")
    &s
  besARGEND

  besRETURN_STRING(strdup(s));

besEND


besFUNCTION(addis)
  const char* s1;
  const char* s2;

  besARGUMENTS("zz")
    &s1, &s2
  besARGEND

  mpz_set_str (op1, s1, BASE);
  mpz_set_str (op2, s2, BASE);
  mpz_add (res, op1, op2);
  char* res_string = mpz_get_str (0, BASE, res);
  besRETURN_STRING(res_string);

besEND


besFUNCTION(subis)
  const char* s1;
  const char* s2;

  besARGUMENTS("zz")
    &s1, &s2
  besARGEND

  mpz_set_str (op1, s1, BASE);
  mpz_set_str (op2, s2, BASE);
  mpz_sub (res, op1, op2);
  char* res_string = mpz_get_str (0, BASE, res);
  besRETURN_STRING(res_string);

besEND


besFUNCTION(mulis)
  const char* s1;
  const char* s2;

  besARGUMENTS("zz")
    &s1, &s2
  besARGEND

  mpz_set_str (op1, s1, BASE);
  mpz_set_str (op2, s2, BASE);
  mpz_mul (res, op1, op2);
  char* res_string = mpz_get_str (0, BASE, res);
  besRETURN_STRING(res_string);

besEND


besFUNCTION(init_memo)
  memo[0] = let("0");
  memo[1] = let("1");
  memo[2] = let("1");
  besRETURNVALUE = NULL;

besEND

besFUNCTION(clean_memo)
  free(memo[0]);
  free(memo[1]);
  free(memo[2]);
  besRETURNVALUE = NULL;

besEND


besFUNCTION(init_globals)
  mpz_init(op1);
  mpz_init(op2);
  mpz_init(res);
  besRETURNVALUE = NULL;

besEND


besFUNCTION(clear_globals)
  mpz_clear(op1);
  mpz_clear(op2);
  mpz_clear(res);
  besRETURNVALUE = NULL;

besEND
gmp.bas

Code: Select all


MODULE GMP

DECLARE SUB ::FIBO              ALIAS "fibo"            LIB "gmp"
DECLARE SUB ::BI_WRITE          ALIAS "writeis"         LIB "gmp"
DECLARE SUB ::BI_LET            ALIAS "letis"           LIB "gmp"
DECLARE SUB ::BI_ADD            ALIAS "addis"           LIB "gmp"
DECLARE SUB ::BI_SUB            ALIAS "subis"           LIB "gmp"
DECLARE SUB ::BI_MUL            ALIAS "mulis"           LIB "gmp"
DECLARE SUB ::BI_INIT_MEMO      ALIAS "init_memo"       LIB "gmp"
DECLARE SUB ::BI_CLEAN_MEMO     ALIAS "clean_memo"      LIB "gmp"
DECLARE SUB ::BI_INIT_GLOBALS   ALIAS "init_globals"    LIB "gmp"
DECLARE SUB ::BI_CLEAR_GLOBALS  ALIAS "clear_globals"   LIB "gmp"

END MODULE
gmpfibo.sb

Code: Select all

' Return the n'th Fibonacci number as a decimal string for integer n

IMPORT gmp.bas

SPLITA STRING(4,"0") BY "" To memo

FUNCTION fibois(n)
  IF n <= 2 THEN
    fibis = GMP::BI_LET(memo[n])
  END IF

  k = (n / 2)
  fk = fibois(k)
  fk1 = fibois(k + 1)
  a = ""
  b = ""
  IF n % 2 = 0 THEN
    a = GMP::BI_ADD(fk1, fk1)
    b = GMP::BI_SUB(a, fk)
    res = GMP::BI_MUL(fk, b)
  ELSE
    a = GMP::BI_MUL(fk, fk)
    b = GMP::BI_MUL(fk1, fk1)
    res = GMP::BI_ADD(a, b)
  END IF
  fibois = res
END FUNCTION

' MAIN

n = 4784969
GMP::BI_INIT_GLOBALS()
GMP::BI_INIT_MEMO()
f = fibois(n)
GMP::BI_WRITE(f)
GMP::BI_CLEAN_MEMO()
GMP::BI_CLEAR_GLOBALS()

PRINT LEN(f),"\n"

Heater
Posts: 13084
Joined: Tue Jul 17, 2012 3:02 pm

Re: A Final Fibonacci Challenge

Tue Jun 11, 2019 11:39 pm

I have no idea what is going on there.

You seem to have rewritten the fibo_strings.c with a bunch of mysterious macros wrapping everything.

I was expecting that the original C code would remain intact and then the BASIC to C interface would only need to marshal strings between the BASIC and the C functions. That actual "glue" code should never need to have any mpz... anything in it.

Heater does not want to return strings as such. I return pointers to strings. Those strings are on the heap. Similarly for the string parameters.

All strings had better be on the heap. Not global constant strings or local strings on the stack. Else the free's will fail.

User avatar
Paeryn
Posts: 2634
Joined: Wed Nov 23, 2011 1:10 am
Location: Sheffield, England

Re: A Final Fibonacci Challenge

Tue Jun 11, 2019 11:40 pm

ScriptBasic wrote:
Tue Jun 11, 2019 10:47 pm
I'm still getting seg faults trying to run the code. Can the C pros on this thread have a peek and help me figure out what the problem is or the best way to debug it?

It gets to the first GMP::BI_ADD() function and seg faults. The original fibo function AIR helped out with works. My guess is SB doesn't like how Heater wants to return strings.
You haven't initialised the global variables op1, op2 and res. You need to do this to allocate the memory they will use, they will automatically reallocate more memory when needed. You should also clear them when you've done with them to free the memory, though if that only happens when the program quits you can get away with not clearing them.

I assume you'll need to do the init in besSUB_START (and the clear in besSUB_FINISH).

interface.c

Code: Select all

/* GMP Extension Module
UXLIBS: -lgmp
*/

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <unistd.h>
#include <memory.h>
#include <gmp.h>
#include "../../basext.h"


int BASE = 32;
mpz_t op1;
mpz_t op2;
mpz_t res;
char* memo[3];


static char* let(const char* s) {
    return strdup(s);
}


/**************************
 Extension Module Functions
**************************/

typedef struct _ModuleObject {
  void *HandleArray;
}ModuleObject,*pModuleObject;


besVERSION_NEGOTIATE
  return (int)INTERFACE_VERSION;
besEND


besSUB_START
  pModuleObject p;

  besMODULEPOINTER = besALLOC(sizeof(ModuleObject));
  if( besMODULEPOINTER == NULL )return 0;

  p = (pModuleObject)besMODULEPOINTER;
  
  // Initialise the three global mpz_t variables
  mpz_init(op1);
  mpz_init(op2);
  mpz_init(res);
  return 0;
besEND


besSUB_FINISH
  pModuleObject p;

  p = (pModuleObject)besMODULEPOINTER;
  if( p == NULL )return 0;
  
  // Release the memory used by the three global mpz_t variables
  mpz_clear(op1);
  mpz_clear(op2);
  mpz_clear(res);
  return 0;
besEND

//...
She who travels light — forgot something.

User avatar
John_Spikowski
Posts: 1329
Joined: Wed Apr 03, 2019 5:53 pm
Location: Anacortes, WA USA
Contact: Website Twitter

Re: A Final Fibonacci Challenge

Wed Jun 12, 2019 12:09 am

Thanks to the code improvements!

Can the memo init / clear also go in the START / FINISH?

I'm still getting seg faults calling GMP math functions and returning strings.
Should I be using the _ui versions of the functions like the fibo function that works?


interface.c

Code: Select all

/* GMP Extension Module
UXLIBS: -lgmp
*/

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <unistd.h>
#include <memory.h>
#include <gmp.h>
#include "../../basext.h"


int BASE = 32;
mpz_t op1;
mpz_t op2;
mpz_t res;
char* memo[3];


static char* let(const char* s) {
    return strdup(s);
}


/**************************
 Extension Module Functions
**************************/

typedef struct _ModuleObject {
  void *HandleArray;
}ModuleObject,*pModuleObject;


besVERSION_NEGOTIATE
  return (int)INTERFACE_VERSION;
besEND


besSUB_START
  pModuleObject p;

  besMODULEPOINTER = besALLOC(sizeof(ModuleObject));
  if( besMODULEPOINTER == NULL )return 0;

  p = (pModuleObject)besMODULEPOINTER;

  // Initialise the three global mpz_t variables
  mpz_init(op1);
  mpz_init(op2);
  mpz_init(res);
  memo[0] = let("0");
  memo[1] = let("1");
  memo[2] = let("1");

  return 0;
besEND


besSUB_FINISH
  pModuleObject p;

  p = (pModuleObject)besMODULEPOINTER;
  if( p == NULL )return 0;

  // Release the memory used by the three global mpz_t variables
  mpz_clear(op1);
  mpz_clear(op1);
  mpz_clear(op2);
  mpz_clear(res);
  free(memo[0]);
  free(memo[1]);
  free(memo[2]);

  return 0;
besEND


/*************
 GMP Functions
*************/


besFUNCTION(fibo)
  int fval;

  besARGUMENTS("i")
    &fval
  besARGEND

  char buf[fval+1];
  memset(buf,0,1);
  mpz_t res;
  mpz_init(res);

  mpz_fib_ui(res, fval);

  gmp_snprintf( buf,sizeof(buf),"%Zd", res );

  besRETURN_STRING(buf);
besEND


besFUNCTION(writeis)
  const char *s;

  besARGUMENTS("z")
    &s
  besARGEND

  mpz_set_str (op1, s, BASE);
  char *buf=mpz_get_str (0, 10, op1);
  puts(buf);
  free(buf);
  besRETURNVALUE = NULL;

besEND


besFUNCTION(letis)
  const char* s;

  besARGUMENTS("z")
    &s
  besARGEND

  besRETURN_STRING(strdup(s));

besEND


besFUNCTION(addis)
  const char* s1;
  const char* s2;

  besARGUMENTS("zz")
    &s1, &s2
  besARGEND

  mpz_set_str (op1, s1, BASE);
  mpz_set_str (op2, s2, BASE);
  mpz_add (res, op1, op2);
  char* res_string = mpz_get_str_ui (0, BASE, res);
  besRETURN_STRING(res_string);

besEND


besFUNCTION(subis)
  const char* s1;
  const char* s2;

  besARGUMENTS("zz")
    &s1, &s2
  besARGEND

  mpz_set_str (op1, s1, BASE);
  mpz_set_str (op2, s2, BASE);
  mpz_sub (res, op1, op2);
  char* res_string = mpz_get_str_ui (0, BASE, res);
  besRETURN_STRING(res_string);

besEND


besFUNCTION(mulis)
  const char* s1;
  const char* s2;

  besARGUMENTS("zz")
    &s1, &s2
  besARGEND

  mpz_set_str (op1, s1, BASE);
  mpz_set_str (op2, s2, BASE);
  mpz_mul (res, op1, op2);
  char* res_string = mpz_get_str_ui (0, BASE, res);
  besRETURN_STRING(res_string);

besEND

gmp.bas

Code: Select all


MODULE GMP

DECLARE SUB ::FIBO              ALIAS "fibo"            LIB "gmp"
DECLARE SUB ::BI_WRITE          ALIAS "writeis"         LIB "gmp"
DECLARE SUB ::BI_LET            ALIAS "letis"           LIB "gmp"
DECLARE SUB ::BI_ADD            ALIAS "addis"           LIB "gmp"
DECLARE SUB ::BI_SUB            ALIAS "subis"           LIB "gmp"
DECLARE SUB ::BI_MUL            ALIAS "mulis"           LIB "gmp"

END MODULE

gmpfibo.sb

Code: Select all

' Return the n'th Fibonacci number as a decimal string for integer n

IMPORT gmp.bas

SPLITA STRING(4,"0") BY "" To memo

FUNCTION fibois(n)
  IF n <= 2 THEN
    fibis = GMP::BI_LET(memo[n])
  END IF

  k = (n / 2)
  fk = fibois(k)
  fk1 = fibois(k + 1)
  a = ""
  b = ""
  IF n % 2 = 0 THEN
    a = GMP::BI_ADD(fk1, fk1)
PRINT "Got Here ADD\n"
    b = GMP::BI_SUB(a, fk)
    res = GMP::BI_MUL(fk, b)
  ELSE
    a = GMP::BI_MUL(fk, fk)
PRINT "Got Here MUL\n"
    b = GMP::BI_MUL(fk1, fk1)
    res = GMP::BI_ADD(a, b)
  END IF
  fibois = res
END FUNCTION

' MAIN

n = 4784969
f = fibois(n)
GMP::BI_WRITE(f)

PRINT LEN(f),"\n"

jalih
Posts: 75
Joined: Mon Apr 15, 2019 3:54 pm

Re: A Final Fibonacci Challenge

Wed Jun 12, 2019 6:02 am

ejolson wrote:
Tue Jun 11, 2019 7:28 pm
Very nice! I wonder if it will run under QEMU on the Pi.

Unfortunately, it may be some time before I have an opportunity to try.
Here is a cleaned up version of the PL/I code. I removed redudant returns that I copied directly and rewrote bigpprint procedure to eliminate tests inside loop.

Code: Select all


*PROCESS MARGINS(1,120) LIBS(SINGLE,STATIC);
*PROCESS OPTIMIZE(2) DFT(REORDER);
*PROCESS LIMITS(FIXEDBIN(63));


/*
    Almost direct PL/I conversion from:

    visual.bas -- Compute the nth Fibonacci Number
    Written January 14, 2018 by Eric Olson
                                            */

fibonacci: package;

 dcl (cmult, rexp, nmax) fixed bin(31);
 dcl radix fixed bin(32) unsigned;
 dcl (cadd, cmax) fixed bin(64) unsigned;

 dcl (a(0:nmax), b(0:nmax), t1(0:nmax), t2(0:nmax)) fixed bin(32) unsigned ctl;
 dcl xy(0:nmax) fixed bin(64) unsigned ctl;


 bigprint: proc(x);
   dcl x(*) fixed bin(32) unsigned;

   dcl i fixed bin(31);
   dcl s0 pic'ZZZZZZZZ9';
   dcl s pic'999999999';

   i = x(0);
   s0 = x(i);
   put edit(trim(s0)) (A);

   do i = x(0)-1 to 1 by -1;
     s = x(i);
     put edit(s) (A);
   end;
 end bigprint;


 biginc: proc(z);
   dcl z(*) fixed bin(32) unsigned;

   dcl i fixed bin(31);
   dcl (t, c) fixed bin(32) unsigned;

   i, c = 1;
   do while(c > 0);
     t = z(i) + 1;
     if t >= radix then
       z(i) = t - radix;
     else
       do;
         z(i) = t;
         c = 0;
       end;
     i += 1;
   end;
 end;


 bigadd: proc(z, x, y);
   dcl (z(*), x(*), y(*)) fixed bin(32) unsigned;

   dcl (i, max) fixed bin(31);

   if y(0) > x(0) then
     max = y(0);
   else
     max = x(0);

   do i = 1 to max + 1;
     z(i) = 0;
   end;

   do i = 1 to max;
     dcl (t, c) fixed bin(31);
     c = 0;
     if i <= x(0) then
       do;
         t = z(i) + x(i);
         if i <= y(0) then
           t += y(i);
       end;
     else
       t = z(i) + y(i);

     if t >= radix then
       do;
         c = 1;
         t -= radix;
       end;
     z(i) = t;
     z(i+1) += c;
   end;
   z(0)=max+1;
   do while(z(z(0)) = 0 & z(0) > 1);
     z(0) = z(0) - 1;
   end;
 end;


 bigdec: proc(z);
   dcl z(*) fixed bin(32) unsigned;

   dcl i fixed bin(31);
   dcl c fixed bin(32) unsigned;

   i, c = 1;
   do while(c > 0);
     if z(i) < 1 then
       z(i) = radix - 1;
     else
       do;
         z(i) -= 1;
         c = 0;
       end;
     i += 1;
   end;
 end;


 bigsub: proc(z, x, y);
   dcl (z(*), x(*), y(*)) fixed bin(32) unsigned;

   dcl i fixed bin(31);
   dcl (t, c) fixed bin(32) unsigned;

   c = 0;
   do i = 1 to y(0);
     t = y(i) + c;
     if x(i) < t then
       do;
         z(i) = x(i) + radix - t;
         c = 1;
       end;
     else
       do;
         z(i) = x(i) - t;
         c = 0;
       end;
   end;
   i = y(0) + 1;
   do while(c > 0);
     if x(i) < 1 then
       z(i) = x(i) + radix - 1;
     else
       do;
         z(i) = x(i) - 1;
         c = 0;
       end;
     i += 1;
   end;
   do while(i <= x(0));
     z(i) = x(i);
     i += 1;
   end;
   z(0) = x(0);
   do while(z(z(0)) = 0 & z(0) > 1);
     z(0) -= 1;
   end;
 end;



 bigmulw: proc(z, x, y);
   dcl (z(*), x(*), y(*)) fixed bin(32) unsigned;

   dcl (i, j, imax) fixed bin(31);

   imax = x(0) + y(0);
   do i = 1 to imax;
     xy(i) = 0;
   end;

   do i = 1 to x(0);
     do j = 1 to y(0);
       dcl (xt, yt) fixed bin(64) unsigned;
       xt = x(i);
       yt = y(j);
       xy(i+j-1) = xy(i+j-1) + xt * yt;
     end;
     if mod(i,cmult) = 0 then
       do;
         dcl (k, kmax) fixed bin(31);
         kmax = i + y(0) - 1;
         do k = 1 to kmax;
           if xy(k) >= cmax then
             do;
               xy(k) = xy(k) - cmax;
               xy(k+1) = xy(k+1) + cadd;
             end;
         end;
       end;
   end;
   do i = 1 to imax - 1;
     xy(i+1)=xy(i+1)+xy(i)/radix;
     z(i) = mod(xy(i), radix);
   end;
   z(imax) = xy(imax);
   do while(z(imax) = 0 & imax > 1);
     imax -= 1;
   end;
   z(0) = imax;
 end;


 bighigh: proc(z, x, n);
   dcl (z(*), x(*)) fixed bin(32) unsigned;
   dcl n fixed bin(31);

   if x(0) <= n then
     do;
       z(0) = 0;
       return;
     end;
   dcl i fixed bin(31);
   do i = n + 1 to x(0);
     z(i-n) = x(i);
   end;
   z(0) = x(0) - n;
 end bighigh;


 biglow: proc(z, x, n);
   dcl (z(*), x(*)) fixed bin(32) unsigned;
   dcl n fixed bin(31);

   if x(0) <= n then
     z(0) = x(0);
   else
     z(0) = n;

   dcl i fixed bin(31);

   do i = 1 to z(0);
     z(i) = x(i);
   end;
 end biglow;


 bigmul: proc(z, x, y) recursive;
   dcl (z(*), x(*), y(*)) fixed bin(32) unsigned;

   if x(0) < 44 | y(0) < 44 then
     do;
       call bigmulw(z,x,y);
       return;
     end;

   dcl (i, M, n) fixed bin(31);

   if x(0) < y(0) then
     M = y(0);
   else
     M = x(0);

   n = M / 2;

   dcl (xl(0:n+1), xh(0:n+1), yl(0:n+1), yh(0:n+1)) fixed bin(32) unsigned ctl;
   allocate xl, xh, yl, yh;
   call biglow(xl, x, n);
   call bighigh(xh, x, n);
   call biglow(yl, y, n);
   call bighigh(yh, y, n);

   dcl (z0(0:M+1), z1(0:M+3), z2(0:M+3)) fixed bin(32) unsigned ctl;
   allocate z0, z1, z2;
   call bigadd(z0,xl,xh);
   call bigadd(z2,yl,yh);
   call bigmul(z1,z0,z2);
   call bigmul(z0,xl,yl);
   call bigmul(z2,xh,yh);
   call bigsub(z1,z1,z0);
   call bigsub(z1,z1,z2);

   dcl imax fixed bin(31);
   imax = x(0) + y(0);
   dcl (t, c) fixed bin(32) unsigned;
   c = 0;
   do i = 1 to imax;
     t=c;
     if i <= z0(0) then t += z0(i);
     if i > n & i-n <= z1(0) then t += z1(i-n);
     if i > 2*n & i - 2 * n <= z2(0) then t += z2(i-2*n);
     c = 0;
     do while(t >= radix);
       t -= radix;
       c += 1;
     end;
     z(i) = t;
   end;
   do while(z(imax) = 0 & imax > 1);
     imax -= 1;
   end;
   z(0) = imax;
   free xl, xh, yl, yh, z0, z1, z2;
 end bigmul;


 fibowork: proc(n) recursive;
   dcl n fixed bin(31);

   if  n = 0 then
     do;
       a(0) = 0;
       b(0) = 1;
       b(1) = 1;
       return;
     end;
   call fibowork(n/2);
   if mod(n,2) = 0 then
     do;
       call bigadd(t1,b,b);
       call bigsub(t2,t1,a);
       call bigmul(a,a,t2);
       call bigmul(b,b,t2);
       if mod(n,4) = 0 then
         call bigdec(b);
       else
         call biginc(b);
     end;
   else
     do;
       call bigadd(t1,a,a);
       call bigadd(t2,t1,b);
       call bigmul(b,b,t2);
       call bigmul(a,a,t2);
       if mod(n,4) = 1 then
         call biginc(a);
       else
         call bigdec(a);
     end;
 end fibowork;


fibo: proc(n);
  dcl n fixed bin(31);

  if n < 2 then
    do;
      b(0) = 1;
      b(1) = n;
      return;
    end;
  call fibowork((n-1)/2);
  if mod(n,2) = 0 then
    do;
      call bigadd(t1,a,a);
      call bigadd(t2,t1,b);
      call bigmul(b,b,t2);
    end;
  else
    do;
      call bigadd(t1,b,b);
      call bigsub(t2,t1,a);
      call bigmul(b,b,t2);
      if mod(n,4) = 1 then
        call bigdec(b);
      else
        call biginc(b);
    end;
 end fibo;


 test: proc options(main);
   dcl n fixed bin(31) init(4784969);
   rexp = 9;
   cmult = 5;
   radix = 10**rexp;
   cadd = radix;
   cadd *= cmult;
   cmax = cadd * radix;
   nmax = 205100;
   allocate a, b, t1, t2, xy;

   call fibo(n);
   call bigprint(b);

   free a, b, t1, t2, xy;
 end test;

end fibonacci;
What I really like about PL/I is that there are no reserved words and most programming errors are catched on compile time. I have never needed debugger when writing PL/I code, usually dump and compiler diagnostic output is enough. Also, one can use "put data" to print all variable contents.

Return to “General programming discussion”