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.
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.
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
You mean is now!danjperron's Python solution using the Golden Ratio is not in the Gighub repository:
Yes it does. I believe this is a correct result:Does it work for smaller values of n under Python2?
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
No idea. Not sure I know what you are asking!Is there a confusion between bits of precision and decimal digits between versions of Python?
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
That what is not, now is!You mean is now!
No idea, I never use Python. But the "python" command here uses python2. After a while I was wondering why it was not finishing...python2!!! Why bother ??
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.No idea, I never use Python. But the "python" command here uses python2. After a while I was wondering why it was not finishing...
I would have thought, to be absolutely certain the rounding error doesn't contaminate the integer part of the result, thatdanjperron wrote: ↑Tue Jun 04, 2019 6:27 pmNow 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.No idea, I never use Python. But the "python" command here uses python2. After a while I was wondering why it was not finishing...
For fibo(478496) you need to calculate with a floating variable containing 1 Million digits + some digits to prevent truncation.
Me too.Damn right. I'm not a Numerics guy by any stretch but that really make my head spin.
You are doing that "Your holding it wrong" thing again.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".
Evidently. It happens. It's whart we invented words for though, right?
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...Heater wrote: ↑Tue Jun 04, 2019 6:41 pmI 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.
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.
Maybe so, but to use Python 3, your program should startdanjperron wrote: ↑Tue Jun 04, 2019 6:10 pmpython2!!! Why bother ?? You know that python2 will be end of life very soon
https://pythonclock.org/
Code: Select all
#!/usr/bin/env python3
Rather than numerical, here is an algebraic approach using the golden ratio along with exact big-number arithmetic. By simple algebra, we have
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.
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
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:danjperron wrote: ↑Tue Jun 04, 2019 6:27 pmNow 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.No idea, I never use Python. But the "python" command here uses python2. After a while I was wondering why it was not finishing...
For fibo(478496) you need to calculate with a floating variable containing 1 Million digits + some digits to prevent truncation.
Code: Select all
real 14m23,821s
user 14m21,188s
sys 0m1,229s
Code: Select all
real 0m10,878s
user 0m10,699s
sys 0m0,061s
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
$
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
$
There's a few major errors in your code,danjperron wrote: ↑Tue Jun 04, 2019 5:37 amok 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); }
Code: Select all
realloc(fibs, sizeof(Fibs_struct)*fibs_size);
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
}
Code: Select all
v = (*pta) * (*(ptb++));
ad 1) Perhaps one of the rare cases where things run much faster using a 64 Bit OS.Heater wrote: ↑Wed Jun 05, 2019 12:40 pmgkreidl,
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:So it would seem the decimal calculation using danjperron's method is only about 15 times slower. Not bad.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 $
But with print enabled:They are about the same speed!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 $
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....
You mean the fibo_phi.py results with and without printing?ad 2) Faster with printing than without? Looks strange.
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
Code: Select all
pi@pi64-aalto:~/fibo_4784969/python$ time python3 fibo_phi.py 4784969 | tail -c 100
real 1m15.625s
user 1m14.986s
sys 0m0.631s
$
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.
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
Then What do I do? if the fibo ask is only 1 and the other 10000000.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?
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.
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 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.MacBook-Air-de-Daniel:~ daniel$ ./t1
FFFFFFFF * 80000000 = 7FFFFFFF80000000
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);
}
Save that for the ponies, it's ghastly. Try the python.org docs instead.ejolson wrote: ↑Wed Jun 05, 2019 10:04 pmI 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