Page 9 of 46

Re: A Final Fibonacci Challenge

Posted: Tue Jun 04, 2019 5:54 pm
by timrowledge
hippy wrote:
Tue Jun 04, 2019 2:13 pm
There's a place for "no talking" libraries. But I also believe there's a place for "scream your head off and bounce from the walls" libraries if that is what people want and it helps them advance.
Exactly. There's plenty of grump old gits like me & Heater; kids and other less jaundiced folk have a right to their stuff too. Scratch, LEGO, micro:bits, snap-circuits, all that stuff has its uses.

Re: A Final Fibonacci Challenge

Posted: Tue Jun 04, 2019 5:56 pm
by Heater
gkreidl,

danjperron's Python solution using the Golden Ratio is now in the Gighub repository:

https://github.com/ZiCog/fibo_4784969/t ... ter/python

Using python3 it produces the result in about 16 seconds on my PC. I have yet to run it under python2 long enough for it to finish!

Code: Select all

$ time python3 fibo_phi.py 4784969 | head -c 32 ; time python3 fibo_phi.py 4784969 | tail -c 32
1000013
107273956418004772293648Traceback (most recent call last):
  File "fibo_phi.py", line 45, in <module>
    print(Fn)
BrokenPipeError: [Errno 32] Broken pipe

real    0m15.840s
user    0m15.125s
sys     0m0.688s
4856539211500699706378405156269

real    0m15.628s
user    0m15.016s
sys     0m0.578s

Re: A Final Fibonacci Challenge

Posted: Tue Jun 04, 2019 6:02 pm
by ejolson
Heater wrote:
Tue Jun 04, 2019 5:56 pm
I have yet to run it under python2 long enough for it to finish!
Does it work for smaller values of n under Python2? Is there a confusion between bits of precision and decimal digits between versions of Python?

Re: A Final Fibonacci Challenge

Posted: Tue Jun 04, 2019 6:10 pm
by danjperron
danjperron's Python solution using the Golden Ratio is not in the Gighub repository:
You mean is now! ;-)

thanks!

python2!!! Why bother ?? You know that python2 will be end of life very soon
https://pythonclock.org/

Re: A Final Fibonacci Challenge

Posted: Tue Jun 04, 2019 6:13 pm
by Heater
ejolson,
Does it work for smaller values of n under Python2?
Yes it does. I believe this is a correct result:

Code: Select all

$ time python2 fibo_phi.py 20000 | head -c 32 ; time python2 fibo_phi.py 20000 | tail -c 32
4188
253116232373236124224015500close failed in file object destructor:
sys.excepthook is missing
lost sys.stderr

real    0m3.053s
user    0m2.453s
sys     0m0.047s
5174220460639857683971213093125

real    0m2.564s
user    0m2.484s
sys     0m0.063s
Is there a confusion between bits of precision and decimal digits between versions of Python?
No idea. Not sure I know what you are asking!

Re: A Final Fibonacci Challenge

Posted: Tue Jun 04, 2019 6:17 pm
by danjperron
this is on python2

Code: Select all

MacBook-Air-de-Daniel:~ daniel$ time python fibo2.py 10
55

real	0m0.056s
user	0m0.034s
sys	0m0.015s
MacBook-Air-de-Daniel:~ daniel$ time python fibo2.py 100
354224848179261915075

real	0m0.060s
user	0m0.036s
sys	0m0.015s
MacBook-Air-de-Daniel:~ daniel$ time python fibo2.py 1000
43466557686937456435688527675040625802564660517371780402481729089536555417949051890403879840079255169295922593080322634775209689623239873322471161642996440906533187938298969649928516003704476137795166849228875

real	0m0.054s
user	0m0.033s
sys	0m0.013s
MacBook-Air-de-Daniel:~ daniel$ time python fibo2.py 2000
4224696333392304878706725602341482782579852840250681098010280137314308584370130707224123599639141511088446087538909603607640194711643596029271983312598737326253555802606991585915229492453904998722256795316982874482472992263901833716778060607011615497886719879858311468870876264597369086722884023654422295243347964480139515349562972087652656069529806499841977448720155612802665404554171717881930324025204312082516817125

real	0m0.057s
user	0m0.034s
sys	0m0.016s

Re: A Final Fibonacci Challenge

Posted: Tue Jun 04, 2019 6:18 pm
by Heater
danjperron,
You mean is now!
That what is not, now is!
python2!!! Why bother ??
No idea, I never use Python. But the "python" command here uses python2. After a while I was wondering why it was not finishing...

Re: A Final Fibonacci Challenge

Posted: Tue Jun 04, 2019 6:27 pm
by danjperron
No idea, I never use Python. But the "python" command here uses python2. After a while I was wondering why it was not finishing...
Now you know why I put days or month on it. I was thinking that could be a good approach because is only one calculation but each variable contains a lot of significant numbers and it takes times. It is by far slower than the integer version.

For fibo(478496) you need to calculate with a floating variable containing 1 Million digits + some digits to prevent truncation.

Re: A Final Fibonacci Challenge

Posted: Tue Jun 04, 2019 6:39 pm
by ejolson
danjperron wrote:
Tue Jun 04, 2019 6:27 pm
No idea, I never use Python. But the "python" command here uses python2. After a while I was wondering why it was not finishing...
Now you know why I put days or month on it. I was thinking that could be a good approach because is only one calculation but each variable contains a lot of significant numbers and it takes times. It is by far slower than the integer version.

For fibo(478496) you need to calculate with a floating variable containing 1 Million digits + some digits to prevent truncation.
I would have thought, to be absolutely certain the rounding error doesn't contaminate the integer part of the result, that

extra digits >= 2 log(n)/log(2)

You appear to get away with fewer, which might mean my bound is overly pessimistic. It would be interesting to test more values to see.

Re: A Final Fibonacci Challenge

Posted: Tue Jun 04, 2019 6:41 pm
by Heater
timrowledge,
Damn right. I'm not a Numerics guy by any stretch but that really make my head spin.
Me too.

As far as I can tell the ratio between successive fibo's converges on the Golden Ratio, give or take a factor of root 5 somewhere.

That turns a series of summations into some kind of geometric series of multiplications. A power of n.

Not sure why the correct integer result pops out of this, rather than an "approximately about" real number kind of result.
Well yes, there's determination and there's "where's the padded armless jacket and the big strong folk that know how to dress you in it safely".
You are doing that "Your holding it wrong" thing again.

Evidently we were talking at cross purposes for a long time as I struggled to get a result out of Smalltalk.

I just wanted to do a simple thing. Run a program that took a single integer parameter in and put a single integer result out. In fact the input parameter was optional. Such a simple thing that we have twenty other languages/dialects/run times that do that with no fuss, even if one has to install them from source code.

In order to do that simple thing in the Smalltalk language you wanted me to buy into a whole graphical user interface, a whole different environment. Let's face it, a whole different operating system as that what Smalltalk was/is.

I would not say that is an example of using the correct "correct powertool". It's more like bringing in a huge excavator to scratch an itch!

By the way, I never did get an explaination as to why the Smalltalk I wrote that runs under GNU Smalltalk does not run under Squeak Smalltalk.

Re: A Final Fibonacci Challenge

Posted: Tue Jun 04, 2019 7:59 pm
by timrowledge
Heater wrote:
Tue Jun 04, 2019 6:41 pm
Evidently we were talking at cross purposes for a long time as I struggled to get a result out of Smalltalk.
Evidently. It happens. It's whart we invented words for though, right?
Heater wrote:
Tue Jun 04, 2019 6:41 pm
I just wanted to do a simple thing. Run a program that took a single integer parameter in and put a single integer result out. In fact the input parameter was optional. Such a simple thing that we have twenty other languages/dialects/run times that do that with no fuss, even if one has to install them from source code.
Now *I* see that latter part to be at least as much work as downloading Squeak, maybe more. I think we're (duh!) coming from very different daily lives, where we each see our normal daily experience as 'normal'. Which is hilarious if you compare it to most of the humans around...

I see a situation where someone asks "I want to draw a square and save a bitmap of it", which I answered by suggesting a Paint tool, downloading it, using a couple of menu options and done - wheras you wanted something more like "here's an OpenSCAD file, run it with `openscd -file ./thisFile.opn`. Which would mean downloading that, maybe conencting to a github thing, compiling etc. Which seems most sensible depends on a lot of personal baggage and needs. In this case the Paint tool is programmable enough that it can be changed to do what you wanted as well as its normal job.
Heater wrote:
Tue Jun 04, 2019 6:41 pm
By the way, I never did get an explaination as to why the Smalltalk I wrote that runs under GNU Smalltalk does not run under Squeak Smalltalk.
Really? Oh, well it's not really any different to the Python thing of 2.7 & 3.whatever. Time, preferences of the people developing the system, perceived needs, all that. Different libraries get bult/chosen and sometimes they don't work identically or include exactly the same things.

GNU Smalltalk was put together by some people that wanted something scripty and they emphasised using just textfiles etc. They chose some 'interesting' syntactical sugar to handle reading in code from those files; no idea why, or what their criteria might have been.
Most other Smalltalk systems have stuck with using the original 'chunk file' format that my old boss Glenn Krasner came up with, because it worked well enough for what it was meant to do. They also stick with the saveable image idea - I'm sure I mentioned that "Smalltalk is saved but not born again". Other code sharing tools have been developed to work different ways - the common Monticello tool uses chunk file stuff wrapped up in a zip file and an agreed convention for the structure of files within that. anothe tool that aims to make use of GitHub storage has been using other strategies in the hope of using git versioning for artefacts as well as plain code. One system uses a server Smalltalk image to serve up code in precompiled lumps so that the client system can be tiny and not need the compiler or many other tools; it does some very cool stuff with talking to html or javascript or Cairo to produce the UI when needed.

Re: A Final Fibonacci Challenge

Posted: Tue Jun 04, 2019 9:14 pm
by scruss
danjperron wrote:
Tue Jun 04, 2019 6:10 pm
python2!!! Why bother ?? You know that python2 will be end of life very soon
https://pythonclock.org/
Maybe so, but to use Python 3, your program should start

Code: Select all

#!/usr/bin/env python3
They're best thought of as two completely separate languages, united by a common fear of HTML text normalization.

Re: A Final Fibonacci Challenge

Posted: Tue Jun 04, 2019 9:34 pm
by ejolson
Heater wrote:
Tue Jun 04, 2019 6:41 pm
timrowledge,
I'm not a Numerics guy by any stretch but that really make my head spin.
As far as I can tell the ratio between successive fibo's converges on the Golden Ratio, give or take a factor of root 5 somewhere.
Rather than numerical, here is an algebraic approach using the golden ratio along with exact big-number arithmetic. By simple algebra, we have

(a+b√5)(c+d√5)=(ac+5bd)+(ad+bc)√5

so it is easy to compute powers of the golden ratio just by keeping track of the rational part and the coefficient of √5. A convenient simplification comes from the value of the golden ratio itself. To see this, take powers of ϕ to obtain

ϕ^1=(1+1√5)/2
ϕ^2=(3+1√5)/2
ϕ^3=(4+2√5)/2
ϕ^4=(7+3√5)/2
ϕ^5=(11+5√5)/2
ϕ^6=(18+8√5)/2
ϕ^7=(29+13√5)/2

Almost miraculously, the integer coefficients of √5 in the numerator are the Fibonacci numbers and the denominator is always two. In particular,

ϕ^n=(???+F(n)√5)/2

The following Pascal program employs the above observation to compute F(n) using only exact big-number arithmetic.

Code: Select all

program golden;

(*  golden.pas -- Compute the nth Fibonacci Number
    Written June 4, 2019 by Eric Olson

    This program demonstrates the expressiveness of FreePascal as
    measured by explicitly computing the Fibonacci numbers F(n) by
    using the formula

      F(n) = (phi^n - (-phi)^(-n))/sqrt(5)

    where phi = (1+sqrt(5))/2.  This program produces the powers
    of phi using conquer and divide but for simplicity does not
    use Karatsuba multiplication.  This results in an algorithm
    with asymptotic time complexity O(n^2). *)

type bignum = array of longint;
type ring5 = record a,b:bignum; end;
const rexp=4;
const radix=10000;
const radix2=radix div 2;

procedure ptrim(var a:bignum);
    var i,na:longint;
begin
    na:=high(a);
    for i:=na downto 0 do begin
        if a[i]<>0 then begin
            setlength(a,i+1);
            exit
        end
    end
end;

procedure pwrite(var a:bignum);
    var i,na:longint;
    var s:string[16];
begin
    na:=high(a);
    for i:=na downto 0 do begin
        writestr(s,a[i]);
        if length(s)>rexp then begin
            writeln('Error not normalized!');
            halt(1)
        end;
        if i=na then write(s)
        else write(copy('0000000000',0,rexp-length(s)),s)
    end;
    writeln
end;

function padd(var a,b:bignum):bignum;
    var na,nb,npadd:longint;
    var i,c,r:longint;
begin
    na:=high(a); nb:=high(b);
    if na>nb then setlength(padd,na+2)
    else setlength(padd,nb+2);
    npadd:=high(padd);
    c:=0;
    for i:=0 to npadd do begin
        r:=c;
        if i<=na then r:=r+a[i];
        if i<=nb then r:=r+b[i];
        c:=r div radix;
        padd[i]:=r mod radix
    end;
    ptrim(padd)
end;

function pmul(var a,b:bignum):bignum;
    var na,nb,npmul:longint;
    var i,j,c,r:longint;
begin
    na:=high(a); nb:=high(b); npmul:=na+nb+4;
    setlength(pmul,npmul+1);
    for i:=0 to npmul do pmul[i]:=0;
    c:=0;
    for i:=0 to na do begin
        for j:=0 to nb do begin
            r:=c+pmul[i+j]+a[i]*b[j];
            c:=r div radix;
            pmul[i+j]:=r mod radix
        end;
        j:=i+nb+1;
        while c>0 do begin
            r:=c+pmul[j];
            c:=r div radix;
            pmul[j]:=r mod radix;
            j:=j+1
        end
    end;
    ptrim(pmul)
end;

procedure pdiv2(var a:bignum);
    var na,i:longint;
begin
    na:=high(a);
    a[0]:=a[0] div 2;
    for i:=1 to na do begin
        if a[i] mod 2=1 then a[i-1]:=a[i-1]+radix2;
        a[i]:=a[i] div 2
    end;
    ptrim(a)
end;

procedure pmul5(var a:bignum);
    var na,i:longint;
    var c,r:longint;
begin
    na:=high(a);
    setlength(a,na+2);
    c:=0;
    for i:=0 to na do begin
        r:=c+5*a[i];
        c:=r div radix;
        a[i]:=r mod radix
    end;
    a[na+1]:=c;
    ptrim(a)
end;

function ring5mul(var x,y:ring5):ring5;
    var t1,t2:bignum;
begin
    t1:=pmul(x.a,y.a);
    t2:=pmul(x.b,y.b); pmul5(t2);
    ring5mul.a:=padd(t1,t2);
    t1:=pmul(x.a,y.b);
    t2:=pmul(x.b,y.a);
    ring5mul.b:=padd(t1,t2)
end;

function goldmul(var x,y:ring5):ring5;
begin
    goldmul:=ring5mul(x,y);
    pdiv2(goldmul.a);
    pdiv2(goldmul.b)
end;

function goldpow(var x:ring5;n:longint):ring5;
    var t1,t2:ring5;
begin
    if n <= 1 then goldpow:=x
    else begin
        t1:=goldpow(x,n div 2);
        t2:=goldmul(t1,t1);
        if n mod 2=0 then goldpow:=t2
        else goldpow:=goldmul(x,t2)
    end
end;

var phinum,z:ring5;
begin
    setlength(phinum.a,1); phinum.a[0]:=1;
    setlength(phinum.b,1); phinum.b[0]:=1;
    z:=goldpow(phinum,4784969);
    pwrite(z.b)
end.
Unfortunately, the code is not very fast because I left out the Karatsuba algorithm and didn't bother with any micro-optimizations similar to the ones in the FreeBasic fibo.bas program.

Since each divide-and-conquer step when calculating ϕ^n uses 4 or 8 big-number multiplications, even in the best case this algebraic method based on the golden ratio will be intrinsically 2 to 4 times slower compared to the doubling formula.

Running on a Ryzen 7 Pro 1700 results in the following output:

Code: Select all

$ fpc -O3 golden.pas
Free Pascal Compiler version 3.0.4 [2019/05/13] for x86_64
Copyright (c) 1993-2017 by Florian Klaempfl and others
$ time ./golden >golden.txt

real    7m3.387s
user    7m3.216s
sys     0m0.003s
$ head -c32 golden.txt ; echo
10727395641800477229364813596225
$ tail -c32 golden.txt
4856539211500699706378405156269
$ md5sum golden.txt 
c926caa99ed1e30fe0d3c966da901738  golden.txt

Re: A Final Fibonacci Challenge

Posted: Wed Jun 05, 2019 4:54 am
by Heater
Oh wait, is that a Pascal solution to the Fibo Challenge you are trying to sneak in there? Or is it not cooked yet?

That's a lot of P languages we are collecting now: pari-gp, perl, prolog, python.. pascal?

God forbid there will be a solution in PHP soon!

Re: A Final Fibonacci Challenge

Posted: Wed Jun 05, 2019 9:16 am
by gkreidl
danjperron wrote:
Tue Jun 04, 2019 6:27 pm
No idea, I never use Python. But the "python" command here uses python2. After a while I was wondering why it was not finishing...
Now you know why I put days or month on it. I was thinking that could be a good approach because is only one calculation but each variable contains a lot of significant numbers and it takes times. It is by far slower than the integer version.

For fibo(478496) you need to calculate with a floating variable containing 1 Million digits + some digits to prevent truncation.
I removed the print function to measure the speed of the algorithm. It's about 80 time slower for fibo(478496), although the decimal module is compiled C code. Using Python3:

Code: Select all

real	14m23,821s
user	14m21,188s
sys	0m1,229s
The bigint version (without timeit and printing the result) using Python3:

Code: Select all

real	0m10,878s
user	0m10,699s
sys	0m0,061s
Both on a RPi 3B+

Re: A Final Fibonacci Challenge

Posted: Wed Jun 05, 2019 12:40 pm
by Heater
gkreidl,

Interesting.

I have somewhat different results running on a Pi 3 (not +) under 64 bit Debian:

With printing of the output removed in both cases:

Code: Select all

$ time python3 fibo.py 4784969

real    0m5.096s
user    0m5.051s
sys     0m0.044s

$ time python3 fibo_phi.py 4784969

real    1m15.632s
user    1m15.019s
sys     0m0.578s
$
So it would seem the decimal calculation using danjperron's method is only about 15 times slower. Not bad.

But with print enabled:

Code: Select all

$ time python3 fibo.py | tail -c 100
180149736853685001275152076875379936330930391815964864885353407167474856539211500699706378405156269

real    1m39.553s
user    1m39.350s
sys     0m0.088s

$ time python3 fibo_phi.py 4784969 | tail -c 100
180149736853685001275152076875379936330930391815964864885353407167474856539211500699706378405156269

real    1m15.440s
user    1m14.838s
sys     0m0.595s
$ 
They are about the same speed!

Conclusion?

I'm not sure, they both end by printing a big integer. Seems the conversion from big decimal representation to integer is a lot faster than the conversion of whatever internal big integer representation to integer takes.

This calls for an adaption of original Python algorithm from using regular integers to big decimals....

Re: A Final Fibonacci Challenge

Posted: Wed Jun 05, 2019 4:44 pm
by Paeryn
danjperron wrote:
Tue Jun 04, 2019 5:37 am
ok I made one in C using tables of 32bits digits.

It is not finished since I need to convert the result from hex to decimal. But it is blazzing fast.

This version uses the same method than fibo.py. Since you don't want gmp, I made my own multiply and add fonction.

Code: Select all

#include <stdlib.h>
#include  <stdio.h>
#include <string.h>

typedef struct{
 int bytesize;
 unsigned int * pt;
}BigN_struct;

extern void print_BigN(BigN_struct *);

typedef struct{
  int index;
  BigN_struct number;
}Fibs_struct;


int fibs_size=100;
Fibs_struct * fibs;
int MaxFibs;

BigN_struct Big2;


Fibs_struct * fiboAdd(int n,BigN_struct * P)
{

  if((MaxFibs+1) >= fibs_size)
  {
     fibs_size+=100;
     realloc(fibs, sizeof(Fibs_struct)*fibs_size);
  }
 fibs[MaxFibs].index=n;
 fibs[MaxFibs].number= *P;
 MaxFibs++;
 return(&(fibs[MaxFibs-1]));
}


Fibs_struct * fiboFind(int n)
{
  int loop;
  for(loop = 0 ;loop < MaxFibs;loop++)
  {
     if(fibs[loop].index == n)
        return  &(fibs[loop]);
  }
  return NULL;
}



void deleteFibs()
 {
   int loop;
   for(loop=0;loop<MaxFibs;loop++)
   {
     if(fibs[loop].number.pt == NULL)
         continue;
     free(fibs[loop].number.pt);
     fibs[loop].number.pt=NULL;
   }
}


void Big_Reduce(BigN_struct *P)
{
   int loop;
   for(loop= P->bytesize-1 ; loop>0;loop--)
     if(P->pt[loop]) break;
   P->bytesize= loop+1;

}



BigN_struct *  Big_Mul(BigN_struct *p1, BigN_struct * p2)
{
   unsigned long long v;
   int loopa,loopb;

   int bytesize= p1->bytesize + p2->bytesize;
   BigN_struct * result =  malloc(sizeof(BigN_struct));

   result->bytesize=bytesize;
   result->pt= malloc(bytesize * sizeof(int));
   for(loopa=0;loopa<bytesize;loopa++)
       result->pt[loopa]=0;
//   memset(result->pt,bytesize,0); 


  int p1s= p1->bytesize;
  int p2s= p2->bytesize;

  unsigned int * pta= p1->pt;
  unsigned int * ptb;
  unsigned int * ptr;

   for(loopa=0;loopa<p1s;loopa++)
    {
      ptr = &(result->pt[loopa]);
      ptb = p2->pt;
     for(loopb=0;loopb<p2s;loopb++)
        {
           v = (*pta)  * (*(ptb++));
           *(ptr++)+=  v & 0xffffffff;
           v >>=32;
           if(v)
            (*ptr)+= v;
        }
      pta++;
    }


  Big_Reduce(result);


  return result;

}


BigN_struct *  Big_Add(BigN_struct *p1, BigN_struct * p2)
{

   int loop;
   int bytesize;


   BigN_struct * result =  malloc(sizeof(BigN_struct));

   if(p1->bytesize > p2->bytesize)
      bytesize = p1->bytesize;
   else
      bytesize = p2->bytesize;

   result->bytesize=bytesize;
   result->pt= malloc((bytesize+1) * sizeof(int));

   unsigned long long  v=0;

   for(loop=0;loop<bytesize;loop++)
    {
       if(p1->bytesize>= loop)
          v+= p1->pt[loop];
       if(p2->bytesize>= loop)
          v+= p2->pt[loop];
       result->pt[loop]= v & 0xffffffff;
        v >>=32;
    }
    if(v)
      {
       result->bytesize++;
       result->pt[result->bytesize-1]=v;
      }
    Big_Reduce(result);
    return(result);
}


void print_BigN(BigN_struct * P)
{
  int loop;
  printf("byte size = %d\n",P->bytesize);
  for(loop=(P->bytesize -1);loop>=0;loop--)
     printf("%08X ",P->pt[loop]);
  printf("\n");
}


Fibs_struct * fibo(int n)
{
  int k;
  BigN_struct * n1,* n2,*n3;
  Fibs_struct * fk = fiboFind(n);
  Fibs_struct * fk1;

  if(fk)
   {
     return fk;
   }


  k = (n + 1) / 2;

  fk = fibo(k);

  fk1 = fibo(k-1);

  if(n&1)
   {
//        result = fk ** 2 + fk1 ** 2
     n1 = Big_Mul(&(fk->number),&(fk->number));
     n2 = Big_Mul(&(fk1->number),&(fk1->number));
     n3 = Big_Add(n1,n2);
   }
   else
   {
//        result = (2 * fk1 + fk) * fk
     n1 = Big_Mul(&(fk1->number),&Big2);
     n2 = Big_Add(n1,&(fk->number));
     n3 = Big_Mul(n2,&(fk->number));
   }
   fk =	fiboAdd(n,n3);
   free(n1);
   free(n2);
   free(n3);
   return(fk);
}



int main(int argc, char * argv[])
{
  int loop;
  int * fibs_idx;
  fibs = malloc(fibs_size * sizeof(Fibs_struct));
  MaxFibs=2;
  // create First and second Number
  fibs[0].index=1;
  fibs[0].number.bytesize=1;
  fibs[0].number.pt = malloc(4);
  fibs[0].number.pt[0] = 1;

  fibs[1].index=2;
  fibs[1].number.bytesize=1;
  fibs[1].number.pt = malloc(4);
  fibs[1].number.pt[0] = 1;

  // create number 2 in BigN
  Big2.bytesize=1;
  Big2.pt = malloc(4);
  Big2.pt[0]=2;



  int value;
  Fibs_struct * V;

 if(argc == 2)
 {

   value =  atol(argv[1]);
   printf("F(%d)=\n",value);
   V = fibo(value);
 }
 else
  V = fibo(1);

  print_BigN(&(V->number));

  deleteFibs();
  free(Big2.pt);


}
There's a few major errors in your code,
Line 31:

Code: Select all

     realloc(fibs, sizeof(Fibs_struct)*fibs_size);
Do you really want to realloc fibs and throw away it's new location if it got moved?

Code: Select all

     Fibs_struct *new_fibs = realloc(fibs, sizeof(Fibs_struct)*fibs_size);
     if (new_fibs)
       fibs = new_fibs;
     else {
        printf("Error resizing fibs\n");
        // Cleanup
      }
Later on, line 105

Code: Select all

           v = (*pta)  * (*(ptb++));
You have v declared as unsigned long long, but you are multiplying two unsigned long values, that will result in an unsigned long value (just the lower half of the result) which is then promoted and stored in v, the upper half of the result is thrown away.

The return type of a multiply is the same type as its operands so to get an unsigned long long answer you need to promote both (or one, C will auto promote the lower type to match) of the operands to unsigned long long. With optimisations on gcc can detect that the upper halves of the operands are always zero and just emit a single 64b = 32b x 32b instruction.
Note, on ARM the smallest mul operates on 32 bits (not counting NEON) so any type smaller (like short) will be auto promoted to 32 bit for multiplication so you can get away with int = short × short but it isn't guaranteed to work on other architectures.

Re: A Final Fibonacci Challenge

Posted: Wed Jun 05, 2019 5:52 pm
by gkreidl
Heater wrote:
Wed Jun 05, 2019 12:40 pm
gkreidl,

Interesting.

I have somewhat different results running on a Pi 3 (not +) under 64 bit Debian:

With printing of the output removed in both cases:

Code: Select all

$ time python3 fibo.py 4784969

real    0m5.096s
user    0m5.051s
sys     0m0.044s

$ time python3 fibo_phi.py 4784969

real    1m15.632s
user    1m15.019s
sys     0m0.578s
$
So it would seem the decimal calculation using danjperron's method is only about 15 times slower. Not bad.

But with print enabled:

Code: Select all

$ time python3 fibo.py | tail -c 100
180149736853685001275152076875379936330930391815964864885353407167474856539211500699706378405156269

real    1m39.553s
user    1m39.350s
sys     0m0.088s

$ time python3 fibo_phi.py 4784969 | tail -c 100
180149736853685001275152076875379936330930391815964864885353407167474856539211500699706378405156269

real    1m15.440s
user    1m14.838s
sys     0m0.595s
$ 
They are about the same speed!

Conclusion?

I'm not sure, they both end by printing a big integer. Seems the conversion from big decimal representation to integer is a lot faster than the conversion of whatever internal big integer representation to integer takes.

This calls for an adaption of original Python algorithm from using regular integers to big decimals....
ad 1) Perhaps one of the rare cases where things run much faster using a 64 Bit OS.
ad 2) Faster with printing than without? Looks strange.

Re: A Final Fibonacci Challenge

Posted: Wed Jun 05, 2019 6:09 pm
by Heater
gkreidl,
ad 2) Faster with printing than without? Looks strange.
You mean the fibo_phi.py results with and without printing?

They are so close and it was only one run. I would instinctively put that down to random jitter caused by Linux scheduling and whatever other processes are going on at the time.

For example, I can get the reverse order just by running it again:

With printing:

Code: Select all

$ time python3 fibo_phi.py 4784969 | tail -c 100
180149736853685001275152076875379936330930391815964864885353407167474856539211500699706378405156269

real    1m15.645s
user    1m15.064s
sys     0m0.575s
$ vim fibo_phi.py 
Without printing:

Code: Select all

[email protected]:~/fibo_4784969/python$ time python3 fibo_phi.py 4784969 | tail -c 100

real    1m15.625s
user    1m14.986s
sys     0m0.631s
$ 

Re: A Final Fibonacci Challenge

Posted: Wed Jun 05, 2019 6:36 pm
by ejolson
Heater wrote:
Wed Jun 05, 2019 12:40 pm
I'm not sure, they both end by printing a big integer. Seems the conversion from big decimal representation to integer is a lot faster than the conversion of whatever internal big integer representation to integer takes.
From the name one might imagine the bigdecimal subroutine library is using a base-10^k internally. The underlying C code was almost surely optimised for 64-bit integers, which would explain why it runs so slow on Raspbian where 64-bit division is emulated in software. It seems likely the intended use was for financial data such as the national debt.

The difficulty with floating point is that all the multiplications needed to produce the quantity φ^n must be performed with more than million-digit precision for the final result to have at least that much. The divide-and-conquer way of computing an integer power is demonstrated by the goldpow function of the Pascal code. As evident, either one or two multiplications are performed per recursion. Assuming the worst case yields an upper bound of

2 log(n)/log(2)

on the total number of bigdecimal multiplications needed to compute the n-th Fibonacci number.

Since the lower bound also scales by log(n), then the runtime Tn needed to compute F(n) ≈ φ^n satisfies

Tn = O(n^α log(n))

assuming an O(n^α) multiplication algorithm is used.

There is no logarithm when using the doubling formulas. Therefore, provided α for the underlying multiplication routines is the same, it is theoretically guaranteed the floating-point method will run slower for n large enough.

When the computer has enough memory and is fast enough that computing with n large is a natural activity, then the use of efficient algorithms becomes more important. Stay tuned for the billion-digit Raspberry-Pi-4 Fibonacci challenge where the value α=log(3)/log(2) of Karatsuba multiplication is much too slow.

Re: A Final Fibonacci Challenge

Posted: Wed Jun 05, 2019 8:05 pm
by Heater
Kind of makes sense.

But, but...in the end both fibo.py and fibo_phi.pi are printing the final result from the same Python integer data type.

So I would have expected that they would both take the same time to do that print(). And fibo_phi.py would take twice as long to run as it does.

Or is my assumption wrong?

Re: A Final Fibonacci Challenge

Posted: Wed Jun 05, 2019 10:04 pm
by ejolson
Heater wrote:
Wed Jun 05, 2019 8:05 pm
But, but...in the end both fibo.py and fibo_phi.pi are printing the final result from the same Python integer data type.
I think fibo_phi is printing the final result as a Decimal type after rounding the fractional part to the nearest integer. Although the result quacks like a duck, the giveaway is that it prints much faster. If you don't mind the colour pink, see also

http://www.pythonlake.com/python/decima ... gral_value

Re: A Final Fibonacci Challenge

Posted: Wed Jun 05, 2019 11:07 pm
by danjperron
@Paeryn

Yes that version has some bugs .

But
Code: Select all

realloc(fibs, sizeof(Fibs_struct)*fibs_size);
Do you really want to realloc fibs and throw away it's new location if it got moved?
Then What do I do? if the fibo ask is only 1 and the other 10000000.
You have v declared as unsigned long long, but you are multiplying two unsigned long values, that will result in an unsigned long value (just the lower half of the result) which is then promoted and stored in v, the upper half of the result is thrown away.




I don't think so! The compiler will deal with it!

Code: Select all

#include <stdlib.h>
#include <stdio.h>

int main(void)
{

  unsigned long a,b;
  unsigned long long c;

  a = 0xffffffff;
  b = 0x80000000;

  c = a * b;

  printf("%8lX * %8lX = %16llX\n",a,b,c);

 return 0;
}
and the result
MacBook-Air-de-Daniel:~ daniel$ ./t1
FFFFFFFF * 80000000 = 7FFFFFFF80000000
And this is my new code with decimal conversion. Still need to work on the decimal conversion to print! it is the slowest part. I wonder if it is not because of memory cache.

I put Heater binary to decimal version from his code and I also add the dabble method.

Code: Select all

#include <stdlib.h>
#include  <stdio.h>
#include <string.h>
#include <stdint.h>



#ifndef nBIT
//  nBit set to 32 bit if not declare in compilation
  #define nBIT 32
#endif



#if nBIT == 32
  #pragma message "nBit set to 32 bit"
  #define bigInt  __uint32_t
  #define bigLong __uint64_t

#elif nBIT == 64
  #pragma message "nBit set to  64 bit"
  #define bigInt  __uint64_t
  #define bigLong __uint128_t

#else
  #pragma message "nBit set to wrong value force 64 bit"
  #define bigInt  __uint64_t
  #define bigLong __uint128_t
#endif


typedef __uint32_t digit;
typedef __uint64_t ddigit;

typedef struct{
 int isize;
 bigInt * pt;
}BigN_struct;

extern void print_BigN_Hex(BigN_struct *);
extern void print_BigN(BigN_struct *);

static digit halfbase=((bigInt)1<<(sizeof(digit)<<2));
static bigInt  base10,zero10;



typedef struct{
  int index;
  BigN_struct number;
}Fibs_struct;


int fibs_size=1000;
Fibs_struct * fibs;
int MaxFibs;

BigN_struct Big2;


Fibs_struct * fiboAdd(int n,BigN_struct * P)
{
 void * pointer;
  if((MaxFibs+1) >= fibs_size)
  {
     fibs_size+=1000;
     pointer=realloc(fibs, sizeof(Fibs_struct)*fibs_size);
  }
 fibs[MaxFibs].index=n;
 fibs[MaxFibs].number= *P;
 MaxFibs++;
 return(&(fibs[MaxFibs-1]));
}


Fibs_struct * fiboFind(int n)
{
  int loop;
  for(loop = 0 ;loop < MaxFibs;loop++)
  {
     if(fibs[loop].index == n)
        return  &(fibs[loop]);
  }
  return NULL;
}



void deleteFibs()
 {
   int loop;
   for(loop=0;loop<MaxFibs;loop++)
   {
     if(fibs[loop].number.pt == NULL)
         continue;
     free(fibs[loop].number.pt);
     fibs[loop].number.pt=NULL;
   }
}


void Big_Reduce(BigN_struct *P)
{
   int loop;
   for(loop= P->isize-1; loop>0;loop--)
     if(P->pt[loop]) break;

   if(loop<0)
    P->isize=1;
   else
    P->isize= loop+1;
  
}



BigN_struct *  Big_Mul(BigN_struct *p1, BigN_struct * p2)
{

   bigLong v;
   int loopa,loopb;

//   printf("====== mul =====\n");
//   printf("p1\n");
//   print_BigN_Hex(p1);
//   printf("p2\n");
//   print_BigN_Hex(p2);

   int isize= p1->isize + p2->isize;
   BigN_struct * result =  malloc(sizeof(BigN_struct));

   if(isize <1)
      isize=1;
   result->isize=isize;
   result->pt= malloc((isize+1) * sizeof(bigInt));

   memset((unsigned char *) result->pt,0, isize * sizeof(bigInt));

  int p1s= p1->isize;
  int p2s= p2->isize;

  bigInt * pta= p1->pt;
  bigInt * ptb;
  bigInt * ptr;
  bigInt * ptr_carry;

   for(loopa=0;loopa<p1s;loopa++)
    {
      ptr = &(result->pt[loopa]);
      ptb = p2->pt;
     for(loopb=0;loopb<p2s;loopb++)
        {
           v = (bigLong) (*pta)  * (bigLong) (*(ptb++));
           ptr_carry= ptr++;

           while(1)
           {
             v += (bigLong) *ptr_carry;
#if nBIT == 32
             *(ptr_carry++) = (bigInt) v & 0xffffffff;
             v>>=32;
#elif nBIT == 64
             *(ptr_carry++) = (bigInt) v & 0xffffffffffffffff;
             v>>=64;
#else
  printf("a - Nbit is not 32 bits\n");
#endif
             if(v==0) break;
           }

        }
      pta++;
    }


  Big_Reduce(result);

  return result;

}


BigN_struct *  Big_Add(BigN_struct *p1, BigN_struct * p2)
{

   int loop;
   int isize;

   BigN_struct * result =  malloc(sizeof(BigN_struct));

   if(p1->isize > p2->isize)
      isize = p1->isize;
   else
      isize = p2->isize;

   result->isize=isize;
   result->pt= malloc((isize+1) * sizeof(bigInt));

   bigLong v=0;

   for(loop=0;loop<isize;loop++)
    {
       if(p1->isize>= loop)
          v+= p1->pt[loop];
       if(p2->isize>= loop)
          v+= p2->pt[loop];
#if nBIT == 32
       result->pt[loop]= v & 0xffffffff;
        v >>=32;
#elif nBIT == 64
       result->pt[loop]= v & 0xffffffffffffffff;
        v >>=64;
#else
   printf("b - nBIT is note 32 bits\n");
#endif
    }
    if(v)
      {
       result->isize++;
       result->pt[result->isize-1]=v;
      }
    Big_Reduce(result);

    return(result);
}



void double_dabble(int n, unsigned short *arr, char **result)
{
    int nbits = 16*n;         /* length of arr in bits */
    int nscratch = nbits/3;   /* length of scratch in bytes */
    char *scratch = malloc((1 + nscratch));
    int i, j, k;
    int smin = nscratch-2;    /* speed optimization */
    unsigned short * pt;
    pt = &arr[n-1];

    for (i=0; i < n; ++i) {
        for (j=0; j < 16; ++j) {
            /* This bit will be shifted in on the right. */
            int shifted_in = ((*pt) & (1 << (15-j)))? 1: 0;

            /* Add 3 everywhere that scratch[k] >= 5. */
            for (k=smin; k < nscratch; ++k)
              scratch[k] += (scratch[k] >= 5)? 3: 0;

            /* Shift scratch to the left by one position. */
            if (scratch[smin] >= 8)
              smin -= 1;
            for (k=smin; k < nscratch-1; ++k) {
                scratch[k] <<= 1;
                scratch[k] &= 0xF;
                scratch[k] |= (scratch[k+1] >= 8);
            }

            /* Shift in the new bit from arr. */
            scratch[nscratch-1] <<= 1;
            scratch[nscratch-1] &= 0xF;
            scratch[nscratch-1] |= shifted_in;
        }
         pt--;
    }

    /* Remove leading zeros from the scratch space. */
    for (k=0; k < nscratch-1; ++k)
      if (scratch[k] != 0) break;
    nscratch -= k;
    memmove(scratch, scratch+k, nscratch+1);

    /* Convert the scratch space from BCD digits to ASCII. */
    for (k=0; k < nscratch; ++k)
      scratch[k] += '0';
    scratch[k]=0;


    /* Resize and return the resulting string. */
    *result = realloc(scratch, nscratch+1);
    return;
}

#if nBIT == 64
void double_dabble64(int n, __uint64_t *arr, char **result)
{
    int nbits = 64*n;         /* length of arr in bits */
    int nscratch = nbits/3;   /* length of scratch in bytes */
    char *scratch = malloc((1 + nscratch));

    int i, j, k;
    int smin = nscratch-2;    /* speed optimization */
    __uint64_t * pt;
    pt = &arr[n-1];

    memset(scratch,0,nscratch);
    for (i=0; i < n; ++i) {
        for (j=0; j < 64; ++j) {
            /* This bit will be shifted in on the right. */
            __uint64_t shifted_in = ((*pt) & ((__uint64_t)1 << (63-j)))? 1: 0;

            /* Add 3 everywhere that scratch[k] >= 5. */
            for (k=smin; k < nscratch; ++k)
              scratch[k] += (scratch[k] >= 5)? 3: 0;

            /* Shift scratch to the left by one position. */
            if (scratch[smin] >= 8)
              smin -= 1;
            for (k=smin; k < nscratch-1; ++k) {
                scratch[k] <<= 1;
                scratch[k] &= 0xF;
                scratch[k] |= (scratch[k+1] >= 8);
            }

            /* Shift in the new bit from arr. */
            scratch[nscratch-1] <<= 1;
            scratch[nscratch-1] &= (__uint64_t) 0xF;
            scratch[nscratch-1] |= shifted_in;
        }
         pt--;
    }

    /* Remove leading zeros from the scratch space. */
    for (k=0; k < nscratch-1; ++k)
      if (scratch[k] != 0) break;
    nscratch -= k;
    memmove(scratch, scratch+k, nscratch+1);

    /* Convert the scratch space from BCD digits to ASCII. */
    for (k=0; k < nscratch; ++k)
      scratch[k] += '0';
    scratch[k]=0;

    /* Resize and return the resulting string. */
    *result = realloc(scratch, nscratch+1);
    return;
}

#endif

// double dabble https://en.wikipedia.org/wiki/Double_dabble
// modified for 32 bits instead of 16 bits

void double_dabble32(int n, __uint32_t *arr, char **result)
{
    int nbits = 32*n;         /* length of arr in bits */
    int nscratch = nbits/3;   /* length of scratch in bytes */
    char *scratch = malloc(1 + nscratch);
    int i, j, k;
    int smin = nscratch-2;    /* speed optimization */

    __uint32_t *pt = &arr[n-1];

    if(scratch == NULL)
    {
       printf("double_dabble32. Unable to malloc\n");
       exit(-1);
    }

    memset(scratch,0,nscratch);

    for (i=0; i < n; ++i) {
        for (j=0; j < 32; ++j) {
            /* This bit will be shifted in on the right. */
            __uint32_t shifted_in = ((*pt) & (1 << (31-j)))? 1: 0;

            /* Add 3 everywhere that scratch[k] >= 5. */
            for (k=smin; k < nscratch; ++k)
              scratch[k] += (scratch[k] >= 5)? 3: 0;

            /* Shift scratch to the left by one position. */
            if (scratch[smin] >= 8)
              smin -= 1;
            for (k=smin; k < nscratch-1; ++k) {
                scratch[k] <<= 1;
                scratch[k] &= 0xF;
                scratch[k] |= (scratch[k+1] >= 8);
            }

            /* Shift in the new bit from arr. */
            scratch[nscratch-1] <<= 1;
            scratch[nscratch-1] &= 0xF;
            scratch[nscratch-1] |= shifted_in;
        }
       pt--;
    }
   
    /* Remove leading zeros from the scratch space. */
    for (k=0; k < nscratch-1; ++k)
      if (scratch[k] != 0) break;
    nscratch -= k;
    memmove(scratch, scratch+k, nscratch+1);

    /* Convert the scratch space from BCD digits to ASCII. */
    for (k=0; k < nscratch; ++k)
      scratch[k] += '0';
    scratch[k]='\0';
    /* Resize and return the resulting string. */
    *result = realloc(scratch, nscratch+1);
    return;
}



static void b10print(digit *d){
    int i;

#if nBIT == 64
    printf("%u",d[d[0]]);
#else
    printf("%u",d[d[0]]);
#endif
    for(i=d[0]-1;i>=1;i--){
#if nBIT == 64
        printf("%0*u",(int)zero10,d[i]);
#else
        printf("%0*u",(int)zero10,d[i]);
#endif
    }
    printf("\n");
}
static inline digit digitlow(digit x){
    return x&(((bigInt)1<<(sizeof(digit)<<2))-1);
}
static inline digit digithigh(digit x){
    return x>>(sizeof(digit)<<2);
}
static void b10add(digit *d,digit c){
    int i,imax;
    for(i=1,imax=d[0];i<=imax;i++){
        d[i]+=c;
        c=d[i]/base10;
        d[i]%=base10;
    }
    if(c){
        do {
            d[++imax]=c;
            c=d[imax]/base10;
            d[imax]%=base10;
        } while(c);
        d[0]=imax;
    }
}
static void b10mul(digit *d,digit x){
    int i,imax;
    for(i=1,imax=d[0];i<=imax;i++) d[i]=d[i]*x;
}


static void b2print(digit *a, int bytesize){
    int i = (bytesize+1)*3*4;
//    printf("b2print size=%d\n",i);
    digit * r=  malloc(i);
    
    if(r == NULL)
     {
      printf("b2print Unable to allocate");
      exit(-1);
     }

    r[0]=1; r[1]=0;
    digit  * pt = (digit  *)  &(a[bytesize -1]);
    for(i=0;i<bytesize;i++){
        b10mul(r,halfbase);
        b10add(r,digithigh(*pt));
        b10mul(r,halfbase);
        b10add(r,digitlow(*pt));
        pt--;
    }
    b10print(r);
   free(r);
}


void print_BigN(BigN_struct * P)
{

#define DABBLE

#ifdef DABBLE

// using double_dabble

    char *text = NULL;
#if nBIT == 64
    double_dabble64(P->isize,P->pt, &text);
#else
    double_dabble32(P->isize,P->pt, &text);
#endif
    printf("%s\n", text);
    free(text);
#else
// or using Heater method
   b2print((digit *) P->pt,P->isize*2);
    
#endif
}



void print_BigN_Hex(BigN_struct * P)
{
  int loop;
  printf("byte size = %d\n",P->isize);
  for(loop=(P->isize -1);loop>=0;loop--)
#if nBIT == 32
     printf("%08X ",P->pt[loop]);
#elif nBIT == 64
     printf("%16llX ",P->pt[loop]);

#endif
  printf("\n");
}


Fibs_struct * fibo(int n)
{
  int k;
  BigN_struct * n1,* n2,* n3;
  Fibs_struct * fk = fiboFind(n);
  Fibs_struct * fk1;

  if(fk)
   {
     return fk;
   }


  k = (n + 1) / 2;

  fk = fibo(k);

  fk1 = fibo(k-1);

  if(n&1)
   {
//        result = fk ** 2 + fk1 ** 2
     n1 = Big_Mul(&(fk->number),&(fk->number));
     n2 = Big_Mul(&(fk1->number),&(fk1->number));
     n3 = Big_Add(n1,n2);
   }
   else
   {
//        result = (2 * fk1 + fk) * fk
     n1 = Big_Mul(&(fk1->number),&Big2);
     n2 = Big_Add(n1,&(fk->number));
     n3 = Big_Mul(n2,&(fk->number));
   }
   fk =	fiboAdd(n,n3);
   if(n1->pt)   free(n1->pt);
   free(n1);
   if(n2->pt)   free(n2->pt);
   free(n2);
   free(n3);
   // n3->pt  is transfer to fibs table
   // do not freed
   return(fk);
}



int main(int argc, char * argv[])
{
  int loop;
  int * fibs_idx;

  for(zero10=0,base10=10;base10<halfbase;base10*=10) zero10++;
  base10/=10;

  fibs = malloc(fibs_size * sizeof(Fibs_struct));
  MaxFibs=2;
  // create First and second Number
  fibs[0].index=1;
  fibs[0].number.isize=1;
  fibs[0].number.pt = malloc(sizeof(bigInt));
  fibs[0].number.pt[0] = 1;

  fibs[1].index=2;
  fibs[1].number.isize=1;
  fibs[1].number.pt = malloc(sizeof(bigInt));
  fibs[1].number.pt[0] = 1;

  // create number 2 in BigN
  Big2.isize=1;
  Big2.pt = malloc(sizeof(bigInt));
  Big2.pt[0]=2;


  int value;
  Fibs_struct * V;

 if(argc == 2)
 {

   value =  atol(argv[1]);
//   printf("F(%d)=\n",value);
   V = fibo(value);
 }
 else
  V = fibo(1);

//  print_BigN_Hex(&(V->number));
  print_BigN(&(V->number));
  deleteFibs();
  free(Big2.pt);


}
This is work in progress.

if you want to compile it using 64 bits __uint64 you need to set nBIT to 64.

ex:
gcc -DnBIT=64 -o fiboC fiboC.c

Re: A Final Fibonacci Challenge

Posted: Wed Jun 05, 2019 11:21 pm
by davidcoton
ejolson wrote:
Wed Jun 05, 2019 10:04 pm
I think fibo_phi is printing the final result as a Decimal type after rounding the fractional part to the nearest integer. Although the result quacks like a duck, the giveaway is that it prints much faster. If you don't mind the colour pink, see also

http://www.pythonlake.com/python/decima ... gral_value
Save that for the ponies, it's ghastly. Try the python.org docs instead.
(After all you expect quacking ducks on a lake, don't you? :lol: )

Re: A Final Fibonacci Challenge

Posted: Wed Jun 05, 2019 11:26 pm
by danjperron
For PHP !

Another approach!

1T bytes SSD or a cheaper HHD.

Your current algorithm figure out fibo using pre-calculated fibonacci number with lower values. Then how long it will take to fill it with all possible number up to the maximum size of the drive.

This way the return value will be fast. and it should not be this long to fill it, The 1 Millions fibonacci long number should fit in the drive!

This is a good for a cluster system to fill the drive.

Once the drive is filled then the number will be output in less than 1 second.