/****************************************************************************
	Kalkulon - A programmable calculator for programmers

	This is an example Kalkulon script.
	To load from within Kalkulon: Load("examples/scriptname.k")
	then you can call functions like: funcname() or funcname(arg0, ...)
	
	Author:		Juergen Holetzeck 2003 - 2008
	e-mail:		contact@kalkulon.de
	homepage:	www.kalkulon.de
****************************************************************************/

//////////////////////////////////////////////////////////////////////////////
// Simple implementation to support big numbers (positive integer numbers with
// arbitrary length)
// JH-2008-02-17
//
// Format of a big number
// 		a big number is a list of the following form: 
// 		{Nn, ..., N2, N1, N0} = Nn * BIG_NUM_BASE**n + ... + N1 * BIG_NUM_BASE**1 + N0 * BIG_NUM_BASE**0
// 		0 <= Nn < BIG_NUM_BASE
// 		e.g. BIG_NUM_BASE = 1e6:  {123,546} = 123000456
// 		big numbers must be positive, signs must be handled separately (see bigint.k)
// 
// Convention for parameter names:
// 		f(xBn, n) -> 	xBn must be of type big number: {Nn, ..., N0}
//				 		n   must be an ordinary number: 0 <= n < BIG_NUM_BASE
//////////////////////////////////////////////////////////////////////////////

//////////////////////////////////////////////////////////////////////////////
// CONFIGURATION
//
// BEWARE: BIG_NUM_BASE*BIG_NUM_BASE must fit into type double without loss!
//
   ::BIG_NUM_BASE = 1e6; // dec base
// ::BIG_NUM_BASE = 0x10000; // hex base 4 digits
// ::BIG_NUM_BASE = 1048576; // hex base 5 digits
// ::BIG_NUM_BASE = 10; // test base
//////////////////////////////////////////////////////////////////////////////

// return: string, convert bignum to string in base BIG_NUM_BASE
Bn2string(xBn) = (
    if( ::BIG_NUM_BASE == 1e6; f0="%i",f1="%06i";
    if( ::BIG_NUM_BASE == 0x10000; f0="%X",f1="%04X";
    if( ::BIG_NUM_BASE == 1048576; f0="%X",f1="%05X";
        Error("no format defined")
    ))),
    s=Size(xBn),
    out=Sprintf(f0,xBn[0]),
    i=1,
    while(i<s; out+=Sprintf(f1,xBn[i]), ++i),
    out
);

// prints bignum in BIG_NUM_BASE
// return: number of characters displayed
printBn(xBn)=(
    if( ::BIG_NUM_BASE == 1e6; prefix="0x";
    if( ::BIG_NUM_BASE == 0x10000; prefix="0x";
    if( ::BIG_NUM_BASE == 1048576; prefix="0x";
        Error("no format defined")
    ))),
	out=prefix+Bn2string(xBn),
	PrintLn(out),
    Size(out)
);


// return: string, convert bignum in base b; BEWARE: 2 <= b < 36 and b < BIG_NUM_BASE 
// e.g. Bn2string(x, 16) returns x as hex string
Bn2string(xBn, b) = (
	symbols="0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ",
	out="",
	b=Bn(b),
	while(xBn!={0}; 
		xBn=divBn(xBn, b),  mod=Bn2float(::divBn_remainder),
		out=Sub(symbols, mod, mod+1)+out
	),
	out
);

// prints bignum in base b with prefix 'prefix'; BEWARE: 2 <= b < 36 and b < BIG_NUM_BASE
// e.g. printBn(x, 16, "0x") prints x as hex number
// return: number of characters displayed
printBn(xBn, b, prefix)=(
	out=prefix+Bn2string(xBn, b),
    PrintLn(out),
    Size(out)
);


// return: bignum build from ordinary number x; x >= 0
Bn(x) =
(
    res = {},
    do( res = PushFront(res, fmod(x, ::BIG_NUM_BASE)), x = floor(x/::BIG_NUM_BASE);
        x
    ),
    res
);

// return: bignum build from number string in base b; 2 <= b < 36
// for b > 10 only capital letters are allowed
// BEWARE: not very efficient for long strings, choose BIG_NUM_BASE properly so 
// that big numbers can easily be built as list directly, e.g. {N100, N99, ..., N0}  
string2Bn(s, b) =
(
	symbols="0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ",
	sym_size=Size(symbols),
	res={0},
	bBn=Bn(b),
	i=0,
	while(Size(s);
		c=Back(s), r=0, while(Sub(symbols, r, r+1)!=c; ++r), s=PopBack(s),
		if(r>=b; Error("wrong format or base")),
		res=addBn(res, mulBn(Bn(r),powBn(bBn, i++)))
	),
	res
);


// return: float (ordinary number), convert bignum to float, 
// BEWARE: for big numbers some precision is lost
Bn2float(xBn) = (s=Size(xBn)-1, i=0, r=0, do(r+=xBn[s-i]*::BIG_NUM_BASE**i; ++i<=s));


// return: normalized bignum (deletes zeros at beginning)
normBn(xBn) =
(
    i=0, s=Size(xBn)-1,
    while(!xBn[i] && i<s; ++i),
    Sub(xBn, i)
);


// return: bignum chopped to n elements
chopBn(xBn, n) = if((s=Size(xBn))>n; Sub(xBn, s-n); xBn);
// return: bignum expanded to n elements (padding with zeros)
expandBn(xBn, n) = (while(Size(xBn)<n; xBn={0}+xBn), xBn);
// return: -1 for xBn < yBN; 0 for xBn==yBn, +1 for xBn > yBn, xBn and yBn must be normalized
compareBn(xBn, yBn) = 
(
	x_size=Size(xBn), y_size=Size(yBn),
	if(x_size<y_size; -1
		; //else
		if(x_size>y_size; 1
			; //else
			i=0, r=0,
			while(i<x_size && !r; if(xBn[i]<yBn[i]; r=-1;if(xBn[i]>yBn[i]; r=1)),++i),
			r
		)
	)
);


// add two bignums
// return: bignum
addBn(xBn, yBn) =
(
    x_size = Size(xBn), y_size = Size(yBn),
    i=x_size-1, j=y_size-1,
    
    res = {},
    carry = 0,
    while(  i>=0 || j>=0 || carry;
            if( i>=0; 
                if( j>=0; z = xBn[i]+yBn[j]+carry; z=xBn[i]+carry); 
                if( j>=0; z = yBn[j]+carry; z=carry)
            ),
            lo = fmod(z, ::BIG_NUM_BASE), carry = floor(z/::BIG_NUM_BASE),
            res = PushFront(res, lo),
            --i, --j
    ),
    res
);

// subtract two bignums
// return: 	(normalized/unnormalized) bignum, negative results are in BIG_NUM_BASE complement
//			 ::subBn_carry = 0, no underflow
//			 ::subBn_carry = 1, underflow (result is negative)
subBn(xBn, yBn) = normBn(subBn_u(xBn, yBn));
subBn_u(xBn, yBn) =
(
    x_size = Size(xBn), y_size = Size(yBn),
    if(x_size<y_size; xBn=Table(y_size-x_size,'0')+xBn, x_size=y_size),
    i=x_size-1, j=y_size-1,
    
    res = {},
    carry = 0,
    while(  i>=0 ;
            if( j>=0; z = xBn[i]-yBn[j]+carry; z=xBn[i]+carry),
            lo = fmod(z, ::BIG_NUM_BASE), carry = floor(z/::BIG_NUM_BASE),
            if(lo<0;lo+=::BIG_NUM_BASE),
            res = PushFront(res, lo),
            --i, --j
    ),
    ::subBn_carry = abs(carry),
    res
);



// multiply two bignums
// return: bignum
mulBn(xBn, yBn) =
(
    x_size = Size(xBn), y_size = Size(yBn),

    res=Table(x_size+y_size-1, '0'),
    j = y_size-1,
    while(  j >= 0;
            i=x_size-1, carry = 0,
            while(  i >= 0 || carry;
                    if( i>=0; 
                        z = yBn[j]*xBn[i]+res[i+j]+carry,
                        lo = fmod(z, ::BIG_NUM_BASE), carry = floor(z/::BIG_NUM_BASE);
                        // else
                        lo = carry, carry = 0
                    ),
                    if(i+j>=0; res[i+j] = lo; res=PushFront(res, lo)),
                    --i
            ),
            --j
    ),
    res
);

// divide two bignums
// return: 	bignum
//			::divBn_remainder = remainder of division (bignum)
divBn(xBn, yBn) =
(
	y_size = Size(yBn),
	// test for short division
 	if(y_size==1;
 		// short division
 		res = divBn_s(xBn,yBn[0]),
 		::divBn_remainder = {::divBn_s_remainder},
 		res
 		; // else long division
		res={},
		// normalize so that yBn[0] >= BIG_NUM_BASE/2
		d = floor(::BIG_NUM_BASE/(yBn[0]+1)),
		xBn = {0} + xBn,
		if(d > 1; xBn = mulBn(xBn, {d}), yBn = mulBn(yBn, {d})),
		x_size = Size(xBn),

		i=0,
		y = yBn[0], 							// left "digit" of yBn
		y1 = yBn[1],
		do(
			// calculate quotient q by making a guess first
			x = xBn[i] * ::BIG_NUM_BASE + xBn[i+1],	// left 2 "digits" of xBn
			q = floor(x/y), r = fmod(x, y),
///*debug*/	PrintLn(Sprintf("x = %g, y = %g, q1 = %g, r = %g", {x, y, q, r})),
			// make a better guess if necessary (max 2 steps are necessary)
/*P1a,P1b*/ if(q >= ::BIG_NUM_BASE || q*y1 > r*::BIG_NUM_BASE + xBn[i+2]; --q, r+=y,
/*P2*/			if(r<::BIG_NUM_BASE && q*y1 > r*::BIG_NUM_BASE + xBn[i+2]; --q)
			),
///*debug*/	PrintLn("\t q2 = "+(string) q),
			// calculate remainder
			rem = Sub(xBn, i, i+y_size+1),
			if(q;
///*debug*/		PrintLn("rem = \t"+(string)rem),
///*debug*/		PrintLn("mul = \t"+(string)mulBn(yBn, {q})),
				rem = subBn_u(rem, mulBn(yBn, {q})),
///*debug*/		PrintLn("res = \t"+(string)rem),
				// check for underflow,
				if(::subBn_carry;
					// guess was one to big
/*P3*/				--q,
					rem = Sub(addBn(rem, yBn), 1)
///*debug*/			,PrintLn("P3 = \t"+(string)rem)
				),
				xBn=Replace(xBn, i, rem)
			),
			res+={q};
			// while
			x_size - ++i > y_size
		),
		::divBn_remainder=divBn_s(normBn(rem), d),
		normBn(res)
	) // endif short/long division	
);


// short division, divides bignum by n
// return: 	bignum
//			::divBn_s_remainder = remainder of division
divBn_s(xBn, n)=
(
    x_size = Size(xBn),
    if( x_size == 1;
        res = {floor(xBn[0]/n)},
        rem = fmod(xBn[0],n)
        ; // else
        res={},
        i = 1,
        rem = xBn[0],	// remainder
        while(  i<x_size;
                t=rem * ::BIG_NUM_BASE+xBn[i],
                z=floor(t/n), rem=fmod(t,n),
                if( i==1;
                    lo = fmod(z, ::BIG_NUM_BASE), carry = floor(z/::BIG_NUM_BASE),
                    if(carry; res+={carry}), res+={lo}
                   	; //else
                    res+={z}
                ),
                ++i
        )
    ),
    ::divBn_s_remainder = rem,
    normBn(res)
);


// power, calculates xBn**n
// return: bignum
powBn(xBn, n) =
(
    result = {1},
    powarray={},
    while(  n>=1;
            if( p=floor(ld(n));
                 n-=2**p, 
                 if(p-1<Size(powarray);
                    tmp=powarray[p-1];
                    tmp=xBn, while(p; powarray+={tmp=mulBn(tmp,tmp)}, --p)
                 ),
                 result=mulBn(result,tmp);
                 // else
                 --n, result=mulBn(result,xBn)
            )
    ),
    result
);


// factorial (simple)
fakBn(n) =
(
    result={1},
    do(result=mulBn(result, Bn(n)); --n),
    result
);

///////////////////////////////
// factorial (swing double, see www.luschny.de/math)
fak1Bn(n) =
(
    ::fak1Bn_N = 1,
    ::fak1Bn_f = {1},
    fak1Bn_rec(n)
);

fak1Bn_rec(n) =
(
    if(n<2; res={1}; res=fak1Bn_rec(n/2), res=mulBn(res,res), res=mulBn(res, fak1Bn_swing(n))),
    res
);

fak1Bn_swing(n) =
(
    // make local copies
    N=::fak1Bn_N,
    f=::fak1Bn_f,
    
    s = N-1+((n-N+1)%4),
    oddN = !(N&1),
    while(  N<=s;
            if(oddN = !oddN; f=mulBn(f,Bn(N)); f=mulBn(f,Bn(4)), f=divBn_s(f,N)),       
            N++
    ),
    
    if( oddN;
        while(  N<=n;
                m = ((N+1)*(N+3)) << 1,
                d = (N*(N+2)) >> 3,
                
                f = mulBn(f, Bn(m)), f=divBn_s(f, d),
                
                N = N+4
        );
        while(  N<=n;
                m = (N*(N+2)) << 1,
                d = ((N+1)*(N+3)) >> 3,

                f = mulBn(f, Bn(m)), f=divBn_s(f, d),

                N = N+4
        )
    ),
    
    ::fak1Bn_N = N,
    ::fak1Bn_f = f
);