diff --git a/modules/bigz/bigz.wx b/modules/bigz/bigz.wx new file mode 100644 index 00000000..9cc4c251 --- /dev/null +++ b/modules/bigz/bigz.wx @@ -0,0 +1,175 @@ + +Namespace bigz + +'#Import "native/*.h" 'ERROR in final exe ( problem in order of importation??? ) +#Import "native/bign.h" +#Import "native/bigq.h" +#Import "native/bigz.h" +#Import "native/bzf.h" + +#Import "makefile" + +#Import "structs" + +Extern + +'bign.h + +Enum BigNumBool +End +Const BN_FALSE:BigNumBool +Const BN_TRUE:BigNumBool + +Enum BigNumCmp +End +Const BN_LT:BigNumCmp +Const BN_EQ:BigNumCmp +Const BN_GT:BigNumCmp + +Enum BigNumCarry +End +Const BN_NOCARRY:BigNumCarry +Const BN_CARRY:BigNumCarry + +Alias BigNumDigit:ULong +Alias BigNum:BigNumDigit Ptr +Alias BigNumLength:UInt + +Function BnnAdd:BigNumCarry( mm:BigNum,ml:BigNumLength,nn:BigNum,nl:BigNumLength,carryin:BigNumCarry ) +Function BnnAddCarry:BigNumCarry( nn:BigNum,nl:BigNumLength,carryin:BigNumCarry ) +Function BnnAndDigits( n:BigNum,d:BigNumDigit ) +Function BnnAssign( mm:BigNum,nn:BigNum,nl:BigNumLength ) +Function BnnCompare:BigNumCmp( mm:BigNum,ml:BigNumLength,nn:BigNum,nl:BigNumLength ) +Function BnnCompareDigits:BigNumCmp( d1:BigNumDigit,d2:BigNumDigit ) +Function BnnComplement( nn:BigNum,nl:BigNumLength ) +Function BnnComplement2( nn:BigNum,nl:BigNumLength ) +Function BnnDivide( nn:BigNum,nl:BigNumLength,dd:BigNum,dl:BigNumLength ) +Function BnnDivideDigit:BigNumDigit( qq:BigNum,nn:BigNum,nl:BigNumLength,d:BigNumDigit ) +Function BnnGetDigit:BigNumDigit( nn:BigNum ) +Function BnnIsPower2:BigNumBool( nn:BigNum,nl:BigNumLength ) +Function BnnIsDigitEven:BigNumBool( d:BigNumDigit ) +Function BnnIsDigitOdd:BigNumBool( d:BigNumDigit ) +Function BnnIsDigitNormalized:BigNumBool( d:BigNumDigit ) +Function BnnIsDigitZero:BigNumBool( d:BigNumDigit ) +Function BnnIsZero:BigNumBool( nn:BigNum,nl:BigNumLength ) +Function BnnMultiply:BigNumCarry( pp:BigNum,pl:BigNumLength,mm:BigNum,ml:BigNumLength,nn:BigNum,nl:BigNumLength ) +Function BnnMultiplyDigit:BigNumCarry( pp:BigNum,pl:BigNumLength,mm:BigNum,ml:BigNumLength,d:BigNumDigit ) +Function BnnNumDigits:BigNumLength( nn:BigNum,nl:BigNumLength ) +Function BnnNumLength:BigNumLength( nn:BigNum,nl:BigNumLength ) +Function BnnNumCount:BigNumLength( nn:BigNum,nl:BigNumLength ) +Function BnnNumLeadingZeroBitsInDigit:BigNumLength( d:BigNumDigit ) +Function BnnOrDigits( n:BigNum,d:BigNumDigit ) +Function BnnSetDigit( nn:BigNum,d:BigNumDigit ) +Function BnnSetToZero( nn:BigNum,nl:BigNumLength ) +Function BnnShiftLeft:BigNumDigit( mm:BigNum,ml:BigNumLength,nbits:BigNumLength ) +Function BnnShiftRight:BigNumDigit( mm:BigNum,ml:BigNumLength,nbits:BigNumLength ) +Function BnnSubtract:BigNumCarry( mm:BigNum,ml:BigNumLength,nn:BigNum,nl:BigNumLength,carryin:BigNumCarry ) +Function BnnSubtractBorrow:BigNumCarry( nn:BigNum,nl:BigNumLength,carryin:BigNumCarry ) +Function BnnXorDigits( n:BigNum,d:BigNumDigit ) + +'bigq.h +'TODO + +'bigz.h + +Enum BzSign +End +Const BZ_MINUS:BzSign +Const BZ_ZERO:BzSign +Const BZ_PLUS:BzSign + +Enum BzCmp +End +Const BZ_LT:BzCmp +Const BZ_EQ:BzCmp +Const BZ_GT:BzCmp + +Enum BzStrFlag +End +Const BZ_UNTIL_END:BzStrFlag +Const BZ_UNTIL_INVALID:BzStrFlag +Const BZ_UNTIL_SLASH:BzStrFlag +Const BZ_UNTIL_SPACE:BzStrFlag + +Struct BigZHeader + Field Size:BigNumLength + Field Sign:BzSign +End + +Struct BigZStruct + Field Header:BigZHeader + Field Digits:BigNumDigit[] +End + +Alias BigZ:BigZStruct Ptr +Alias BzChar:Byte +Alias BzUInt:UInt +Alias BzInt:Int +Alias BzLDouble:Double +Alias BzSeed:UInt + +Function BzVersion:CString() +Function BzCreate:BigZ( Size:BigNumLength ) +Function BzNumDigits:BigNumLength( z:BigZ ) +Function BzLength:BigNumLength( z:BigZ ) +Function BzCopy:BigZ( z:BigZ ) +Function BzNegate:BigZ( z:BigZ ) +Function BzAbs:BigZ( z:BigZ ) +Function BzCompare:BzCmp( y:BigZ,z:BigZ ) +Function BzAdd:BigZ( y:BigZ,z:BigZ ) +Function BzSubtract:BigZ( y:BigZ,z:BigZ ) +Function BzMultiply:BigZ( y:BigZ,z:BigZ ) +Function BzDivide:BigZ( y:BigZ,z:BigZ,r:BigZ Ptr ) +Function BzDiv:BigZ( y:BigZ,z:BigZ ) +Function BzTruncate:BigZ( y:BigZ,z:BigZ ) +Function BzFloor:BigZ( y:BigZ,z:BigZ ) +Function BzCeiling:BigZ( y:BigZ,z:BigZ ) +Function BzRound:BigZ( y:BigZ,z:BigZ ) +Function BzMod:BigZ( y:BigZ,z:BigZ ) +Function BzRem:BigZ( y:BigZ,z:BigZ ) +Function BzPow:BigZ( base:BigZ,exponent:BzUInt ) +Function BzIsEven:BigNumBool( y:BigZ ) +Function BzIsOdd:BigNumBool( y:BigZ ) +Function BzToString:CString( z:BigZ,base:BigNumDigit,sign:Int ) +Function BzStrLen:Int( s:CString ) 'size_t +Function BzToStringBuffer:CString( z:BigZ,base:BigNumDigit,sign:Int,buf:CString,len:Int Ptr ) +Function BzToStringBufferExt:CString( z:BigZ,base:BigNumDigit,sign:Int,buf:CString,len:Int Ptr,slen:Int Ptr ) +Function BzFromStringLen:BigZ( s:CString,len:Int,base:BigNumDigit,flag:BzStrFlag ) +Function BzFromString:BigZ( s:CString,base:BigNumDigit,flag:BzStrFlag ) +Function BzFromInteger:BigZ( i:BzInt ) +Function BzFromUnsignedInteger:BigZ( i:BzUInt ) +Function BzToInteger:BzInt( z:BigZ ) +Function BzToDouble:Double( z:BigZ ) +Function BzToLongDouble:BzLDouble( z:BigZ ) +Function BzToIntegerPointer:Int( z:BigZ,p:BzInt Ptr ) +Function BzToUnsignedInteger:BzUInt( z:BigZ ) +Function BzToUnsignedIntegerPointer:Int( z:BigZ,p:BzUInt Ptr ) +Function BzFromBigNum:BigZ( n:BigNum,nl:BigNumLength ) +Function BzToBigNum:BigNum( z:BigZ,nl:BigNumLength Ptr ) +Function BzTestBit:BigNumBool( bit:BigNumLength,z:BigZ ) +Function BzBitCount:BigNumLength( z:BigZ ) +Function BzNot:BigZ( z:BigZ ) +Function BzAnd:BigZ( y:BigZ,z:BigZ ) +Function BzOr:BigZ( y:BigZ,z:BigZ ) +Function BzXor:BigZ( y:BigZ,z:BigZ ) +Function BzNand:BigZ( x:BigZ,y:BigZ ) +Function BzNor:BigZ( x:BigZ,y:BigZ ) +Function BzEqv:BigZ( x:BigZ,y:BigZ ) +Function BzAndC1:BigZ( x:BigZ,y:BigZ ) +Function BzAndC2:BigZ( x:BigZ,y:BigZ ) +Function BzOrC1:BigZ( x:BigZ,y:BigZ ) +Function BzOrC2:BigZ( x:BigZ,y:BigZ ) +Function BzAsh:BigZ( y:BigZ,n:Int ) +Function BzSqrt:BigZ( z:BigZ ) +Function BzLcm:BigZ( y:BigZ,z:BigZ ) +Function BzGcd:BigZ( y:BigZ,z:BigZ ) +Function BzRandom:BigZ( n:BigZ,seed:BzSeed Ptr ) +Function BzModExp:BigZ( base:BigZ,exponent:BigZ,modulus:BigZ ) +Function BzHash:BzUInt( z:BigZ ) + +'Function BnDebug( m:CString,bzstr:CString,n:BigNum,nl:BigNumLength,sign:BzSign ) +'Function BzDebug( m:CString,y:BigZ ) + +'bzf.h + +Function BzFactorial:BigZ( z:BigZ ) diff --git a/modules/bigz/makefile.wx b/modules/bigz/makefile.wx new file mode 100644 index 00000000..613459cb --- /dev/null +++ b/modules/bigz/makefile.wx @@ -0,0 +1,8 @@ + +Namespace bigz + +#If __TARGET__="linux" + + #Import "makefile_linux" + +#Endif diff --git a/modules/bigz/makefile_linux.wx b/modules/bigz/makefile_linux.wx new file mode 100644 index 00000000..912a52f5 --- /dev/null +++ b/modules/bigz/makefile_linux.wx @@ -0,0 +1,7 @@ + +Namespace bigz + +#Import "native/bign.c" +#Import "native/bigq.c" +#Import "native/bigz.c" +#Import "native/bzf.c" diff --git a/modules/bigz/module.json b/modules/bigz/module.json new file mode 100644 index 00000000..3e9fc888 --- /dev/null +++ b/modules/bigz/module.json @@ -0,0 +1,8 @@ +{ + "module":"bigz", + "about":"BIGZ library wrapper", + "author":"Denise Amiga", + "version":"1.0.0", + "support":"http://github.com/denise-amiga/bigz-module", + "depends":[] +} diff --git a/modules/bigz/native/bign.c b/modules/bigz/native/bign.c new file mode 100644 index 00000000..19436310 --- /dev/null +++ b/modules/bigz/native/bign.c @@ -0,0 +1,1233 @@ +/* + * Simplified BSD License + * + * Copyright (c) 1988-1989, Digital Equipment Corporation & INRIA. + * Copyright (c) 1992-2023, Eligis + * All rights reserved. + * + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions are + * met: + * + * o Redistributions of source code must retain the above copyright + * notice, this list of conditions and the following disclaimer. + * o Redistributions in binary form must reproduce the above copyright + * notice, this list of conditions and the following disclaimer in the + * documentation and/or other materials provided with the distribution. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS + * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT + * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR + * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT + * HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, + * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT + * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, + * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY + * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT + * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE + * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + */ + +/** + * @file bign.c + * @brief the kernel written in pure C (it uses no C library) + * Description of types and constants. + * + * Several conventions are used in the commentary: + * - A "BigNum" is the name for an infinite-precision number. + * - Capital letters (e.g., "N") are used to refer to the value of BigNums. + * - The word "digit" refers to a single BigNum digit. + * - The notation "Size(N)" refers to the number of digits in N, + * which is typically passed to the subroutine as "nl". + * - The notation "Length(N)" refers to the number of digits in N, + * not including any leading zeros. + * - The word "Base" is used for the number 2 ** BN_DIGIT_SIZE, where + * BN_DIGIT_SIZE is the number of bits in a single BigNum digit. + * - The expression "BBase(N)" is used for Base ** NumDigits(N). + * - The term "leading zeros" refers to any zeros before the most + * significant digit of a number. + * + * In the code, we have: + * - "nn" is a pointer to a big number, + * - "nl" is the number of digits from nn, + * - "d" is a digit. + * @version 2.1.0 + * @copyright Digital Equipment Corporation & INRIA, 1988-1989. + * @copyright Eligis, 1992-2023 + * @authors B. Serpette + * @authors J. Vuillemin + * @authors J.C. Hervé + * @authors C. Jullien + * $Revision: 392 $ + * $Date: 2023-02-04 06:21:06 +0000 (Sat, 04 Feb 2023) $ + */ + +#if !defined(__BIGN_H) +#include "./bign.h" +#endif + +static void +BnnDivideHelper(BigNum nn, BigNumLength nl, BigNum dd, BigNumLength dl); + +/** + * BnnSetToZero. + * Sets all the specified digits of the BigNum to BN_ZERO (0). + * @param [in, out] nn BigNum + * @param [in] nl BigNumLength + */ +void +BnnSetToZero(BigNum nn, BigNumLength nl) { + BigNumLength d; + + for (d = 0; d < nl; ++d) { + nn[d] = BN_ZERO; + } +} + +/** + * BnnAssign. + * Copies N => M + * @param [in, out] mm BigNum + * @param [in] nn const BigNum + * @param [in] nl BigNumLength + */ +void +BnnAssign(BigNum mm, const BigNum nn, BigNumLength nl) { + int d; + + if ((mm < nn) || (mm > (nn + nl))) { + /* + * no memory overlap using classic loop + */ + for (d = 0; d < (int)nl; ++d) { + mm[d] = nn[d]; + } + } else if (mm > nn) { + /* + * memory overlap, loop starting from most significant digit + */ + for (d = (int)(nl - 1); d >= 0; --d) { + mm[d] = nn[d]; + } + } +} + +/** + * BnnSetDigit. + * Sets a single digit of N to the passed value + * @param [in, out] nn BigNum + * @param [in] d BigNumDigit + */ +void +BnnSetDigit(BigNum nn, BigNumDigit d) { + *nn = d; +} + +/** + * BnnGetDigit. + * Returns the single digit pointed by N + * @param [in] nn BigNum + * @return BigNumDigit + */ +BigNumDigit +BnnGetDigit(const BigNum nn) { + return *nn; +} + +/** + * BnnNumDigits. + * Returns the number of digits of N, not counting leading zeros. + * @param [in] nn BigNum + * @param [in] nl BigNumLength + * @return BigNumLength + */ +BigNumLength +BnnNumDigits(const BigNum nn, BigNumLength nl) { + int d; + + /* + * loop starting from most significant digit + */ + + for (d = (int)(nl - 1); d >= 0; --d) { + if (nn[d] != BN_ZERO) { + /* + * length = d+1 + */ + return (BigNumLength)(d + 1); + } + } + + return (BigNumLength)1; +} + +/** + * BnnNumLength. + * Returns the number of bits of N, not counting leading zeros. + * @param [in] nn BigNum + * @param [in] nl BigNumLength + * @return BigNumLength + */ +BigNumLength +BnnNumLength(const BigNum nn, BigNumLength nl) { + const BigNumDigit d = nn[nl - 1]; + int i; + + for (i = (int)(BN_DIGIT_SIZE - 1); i >= 0; --i) { + if ((d & (BN_ONE << (BigNumLength)i)) != 0) { + BigNumLength bits = (BigNumLength)i + 1; + return (BigNumLength)(((nl-1) * BN_DIGIT_SIZE) + bits); + } + } + + return 0; +} + +/** + * BnnNumCount. + * Returns the count of bits set of N. + * @param [in] nn BigNum + * @param [in] nl BigNumLength + * @return BigNumLength + */ +BigNumLength +BnnNumCount(const BigNum nn, BigNumLength nl) { + BigNumLength count = 0; + int j; + + for (j = 0; j < (int)nl; ++j) { + const BigNumDigit d = nn[j]; + int i; + + for (i = (int)(BN_DIGIT_SIZE - 1); i >= 0; --i) { + if ((d & (BN_ONE << (BigNumLength)i)) != 0) { + ++count; + } + } + } + + return count; +} + +/** + * BnnNumLeadingZeroBitsInDigit. + * Returns the number of leading zero bits in a digit + * @param [in] d BigNumDigit + * @return BigNumLength + */ +BigNumLength +BnnNumLeadingZeroBitsInDigit(BigNumDigit d) { + BigNumDigit mask = (BigNumDigit)(BN_ONE << (BN_DIGIT_SIZE - 1)); + BigNumLength p; + + if (d == BN_ZERO) { + return (BigNumLength)BN_DIGIT_SIZE; + } + + for (p = 0; (d & mask) == 0; ++p) { + mask >>= 1; + } + + return p; +} + +/** + * BnnIsPower2. + * Returns BN_TRUE iff nn is a power of 2. + * @param [in] nn BigNum + * @param [in] nl BigNumLength + * @return BigNumBool + */ +BigNumBool +BnnIsPower2(const BigNum nn, BigNumLength nl) { + BigNumLength i; + BigNumLength nbits; + BigNumDigit d; + + /* + * The n-1 digits must be 0 + */ + + for (i = 0; i < (nl - 1); ++i) { + if (nn[i] != BN_ZERO) { + return BN_FALSE; + } + } + + /* + * There must be only 1 bit set on the last Digit. + */ + + d = nn[i]; + nbits = 0; + + for (i = 0; i < (BigNumLength)BN_DIGIT_SIZE; ++i) { + if ((d & (BN_ONE << i)) != 0) { + if (nbits++ > 0) { + /* + * More than two digits. + */ + return BN_FALSE; + } + } + } + + return BN_TRUE; +} + +/** + * BnnIsDigitZero. + * Returns BN_TRUE iff digit = 0 + * @param [in] d igNumDigit + * @return BigNumBool + */ +BigNumBool +BnnIsDigitZero(BigNumDigit d) { + return (BigNumBool)(d == 0); +} + +/** + * BnnIsDigitNormalized. + * Returns BN_TRUE iff Base/2 <= digit < Base i.e. if digit's leading + * bit is 1 + * @param [in] d BigNumDigit + * @return BigNumBool + */ +BigNumBool +BnnIsDigitNormalized(BigNumDigit d) { + if ((d & (BN_ONE << (BN_DIGIT_SIZE - 1))) != 0) { + return BN_TRUE; + } else { + return BN_FALSE; + } +} + +/** + * BnnIsDigitOdd. + * Returns BN_TRUE iff digit is odd. + * @param [in] d BigNumDigit + * @return BigNumBool + */ +BigNumBool +BnnIsDigitOdd(BigNumDigit d) { + if ((d & 1) != 0) { + return BN_TRUE; + } else { + return BN_FALSE; + } +} + +/** + * BnnIsDigitEven. + * Returns BN_TRUE iff digit is even + * @param [in] d BigNumDigit + * @return BigNumBool + */ +BigNumBool +BnnIsDigitEven(BigNumDigit d) { + if ((d & 1) == 0) { + return BN_TRUE; + } else { + return BN_FALSE; + } +} + +/** + * BnnCompareDigits. + * Returns: + * - BN_GT if d1 > d2 + * - BN_EQ if d1 = d2 + * - BN_LT if d1 < d2 + * @param [in] d1 BigNumDigit + * @param [in] d2 BigNumDigit + * @return BigNumCmp + */ +BigNumCmp +BnnCompareDigits(BigNumDigit d1, BigNumDigit d2) { + return (BigNumCmp)((d1 > d2) ? BN_GT : (d1 == d2 ? BN_EQ : BN_LT)); +} + +/** + * BnnComplement. + * Performs the computation BBase(N) - N - 1 => N + * @param [in, out] nn BigNum + * @param [in] nl BigNumLength + */ +void +BnnComplement(BigNum nn, BigNumLength nl) { + BigNumLength d; + + for (d = 0; d < nl; ++d) { + nn[d] ^= BN_COMPLEMENT; + } +} + +/** + * BnnComplement2. + * Performs the computation neg(N) => N + * @param [in, out] nn BigNum + * @param [in] nl BigNumLength + */ +void +BnnComplement2(BigNum nn, BigNumLength nl) { + BigNumDigit one = BN_ONE; + + /* + * Initialize constants + */ + + BnnComplement(nn, nl); + (void)BnnAdd(nn, nl, &one, (BigNumLength)1, BN_NOCARRY); +} + +/** + * BnnAndDigits. + * Returns the logical computation n[0] AND d in n[0] + * @param [in, out] n BigNum + * @param [in] d BigNumDigit + */ +void +BnnAndDigits(BigNum n, BigNumDigit d) { + *n &= d; +} + +/** + * BnnOrDigits. + * Returns the logical computation n[0] OR d2 in n[0]. + * @param [in, out] n BigNum + * @param [in] d BigNumDigit + */ +void +BnnOrDigits(BigNum n, BigNumDigit d) { + *n |= d; +} + +/** + * BnnXorDigits. + * Returns the logical computation n[0] XOR d in n[0]. + * @param [in, out] n BigNum + * @param [in] d BigNumDigit + */ +void +BnnXorDigits(BigNum n, BigNumDigit d) { + *n ^= d; +} + +/* + * Shift operations + */ + +/** + * BnnShiftLeft. + * Shifts M left by "nbits", filling with 0s. + * Returns the leftmost "nbits" of M in a digit. + * @pre 0 <= nbits < BN_DIGIT_SIZE. + * @param [in, out] mm BigNum + * @param [in] ml BigNumLength + * @param [in] nbits BigNumLength + * @return BigNumDigit + */ +BigNumDigit +BnnShiftLeft(BigNum mm, BigNumLength ml, BigNumLength nbits) { + BigNumDigit res = BN_ZERO; + + if (nbits != 0) { + BigNumLength rnbits = (BigNumLength)(BN_DIGIT_SIZE - nbits); + BigNumLength evenlen = (ml & ~(BigNumLength)1); + BigNumLength d; + + /* + * Loop is now unrooled two BigNumDigit at a time. + */ + + for (d = 0; d < evenlen; ++d) { + BigNumDigit save0; + BigNumDigit save1; + save0 = mm[d]; + mm[d] = (save0 << nbits) | res; + save1 = mm[++d]; + mm[d] = (save1 << nbits) | (save0 >> rnbits); + res = save1 >> rnbits; + } + + if (ml != evenlen) { + BigNumDigit save = mm[d]; + mm[d] = (save << nbits) | res; + res = save >> rnbits; + } + } + + return res; +} + +/** + * BnnShiftRight. + * Shifts M right by "nbits", filling with 0s. + * Returns the rightmost "nbits" of M in a digit. + * @pre 0 <= nbits < BN_DIGIT_SIZE. + * @param [in, out] mm BigNum + * @param [in] ml BigNumLength + * @param [in] nbits BigNumLength + * @return BigNumDigit + */ +BigNumDigit +BnnShiftRight(BigNum mm, BigNumLength ml, BigNumLength nbits) { + BigNumDigit res = BN_ZERO; + + if (nbits != 0) { + const BigNumLength lnbits = (BigNumLength)BN_DIGIT_SIZE - nbits; + + /* + * loop starting from most significant digit + */ + + if ((ml & (BigNumLength)1) != (BigNumLength)0) { + /* + * Odd number of digits, start with most significant + * digit, then loop on even number of digits. + */ + BigNumDigit save = mm[--ml]; + mm[ml] = (save >> nbits); /* res==0, no need to | res */ + res = save << lnbits; + } + + /* + * Loop is now unrooled two digits at a time. + */ + + while (ml != (BigNumLength)0) { + BigNumDigit save0; + BigNumDigit save1; + save0 = mm[--ml]; + mm[ml] = (save0 >> nbits) | res; + save1 = mm[--ml]; + mm[ml] = (save1 >> nbits) | (save0 << lnbits); + res = save1 << lnbits; + } + } + + return res; +} + +/* + * Additions + */ + +/** + * BnnAddCarry. + * Performs the sum N + CarryIn => N. + * Returns the CarryOut. + * @param [in, out] nn BigNum + * @param [in] nl BigNumLength + * @param [in] carryin BigNumCarry + * @return BigNumCarry + */ +BigNumCarry +BnnAddCarry(BigNum nn, BigNumLength nl, BigNumCarry carryin) { + if (carryin == BN_NOCARRY) { + return BN_NOCARRY; + } else if (nl == 0) { + return BN_CARRY; + } else { + BigNumLength d; + + for (d = 0; d < nl; ++d) { + if (++nn[d] != 0) { + return BN_NOCARRY; + } + } + + return BN_CARRY; + } +} + +/** + * BnnAdd. + * Performs the sum M + N + CarryIn => M. + * Returns the CarryOut. + * @pre Size(M) >= Size(N). + * @param [in, out] mm BigNum + * @param [in] ml BigNumLength + * @param [in] nn BigNum + * @param [in] nl BigNumLength + * @param [in] carryin BigNumCarry + * @return BigNumCarry + */ +BigNumCarry +BnnAdd(BigNum mm, + BigNumLength ml, + const BigNum nn, + BigNumLength nl, + BigNumCarry carryin) { + BigNumProduct c = (BigNumProduct)carryin; + BigNumLength i; + + for (i = 0; i < nl; ++i) { + BigNumProduct save = (BigNumProduct)*mm; + + c += save; + if (c < save) { + *(mm++) = nn[i]; + c = (BigNumProduct)1; + } else { + save = (BigNumProduct)nn[i]; + c += save; + *(mm++) = (BigNumDigit)c; + c = (BigNumProduct)((c < save) ? 1 : 0); + } + } + + return BnnAddCarry(mm, ml-nl, ((c == 0) ? BN_NOCARRY : BN_CARRY)); +} + +/* + * Subtraction + */ + +/** + * BnnSubtractBorrow. + * Performs the difference N + CarryIn - 1 => N. + * Returns the CarryOut. + * @param [in, out] nn BigNum + * @param [in] nl BigNumLength + * @param [in] carryin BigNumCarry + * @return BigNumCarry + */ +BigNumCarry +BnnSubtractBorrow(BigNum nn, BigNumLength nl, BigNumCarry carryin) { + if (carryin == BN_CARRY) { + return BN_CARRY; + } else if (nl == 0) { + return BN_NOCARRY; + } else { + BigNumLength d; + + for (d = 0; d < nl; ++d) { + if (nn[d]-- != 0) { + return BN_CARRY; + } + } + + return BN_NOCARRY; + } +} + +/** + * BnnSubtract. + * Performs the difference M - N + CarryIn - 1 => M. + * Returns the CarryOut. + * @pre Size(M) >= Size(N). + * @param [in, out] mm BigNum + * @param [in] ml BigNumLength + * @param [in] nn BigNum + * @param [in] nl BigNumLength + * @param [in] carryin BigNumCarry + * @return BigNumCarry + */ +BigNumCarry +BnnSubtract(BigNum mm, + BigNumLength ml, + const BigNum nn, + BigNumLength nl, + BigNumCarry carryin) { + BigNumProduct c = (BigNumProduct)((carryin == BN_CARRY) ? 1 : 0); + BigNumDigit invn; + BigNumProduct save; + BigNumLength i; + + for (i = 0; i < nl; ++i) { + save = (BigNumProduct)*mm; + invn = nn[i] ^ BN_COMPLEMENT; + c += save; + + if (c < save) { + *(mm++) = invn; + c = (BigNumProduct)1; + } else { + c += invn; + *(mm++) = (BigNumDigit)c; + c = (BigNumProduct)((c < invn) ? 1 : 0); + } + } + + if (c == 0) { + return BnnSubtractBorrow(mm, ml-nl, BN_NOCARRY); + } else { + return BnnSubtractBorrow(mm, ml-nl, BN_CARRY); + } +} + +/* + * Multiplication + */ + +/** @cond */ +#define LOW(x) (BigNumDigit)(x & ((BN_ONE << (BN_DIGIT_SIZE / 2)) - 1)) +#define HIGH(x) (BigNumDigit)(x >> (BN_DIGIT_SIZE / 2)) +#define L2H(x) (BigNumDigit)(x << (BN_DIGIT_SIZE / 2)) + +#define UPDATE_S(c, V, X3) c += V; X3 += (BigNumDigit)(c < V); +/** @endcond */ + +/** + * BnnMultiplyDigit. + * Performs the product: + * - Q = P + M * d + * - BB = BBase(P) + * - Q mod BB => P + * - Q div BB => CarryOut + * Returns the CarryOut. + * @pre Size(P) >= Size(M) + 1. + * @param [in, out] pp BigNum + * @param [in] pl BigNumLength + * @param [in] mm BigNum + * @param [in] ml BigNumLength + * @param [in] d BigNumDigit + * @return BigNumCarry + */ +BigNumCarry +BnnMultiplyDigit(BigNum pp, + BigNumLength pl, + const BigNum mm, + BigNumLength ml, + BigNumDigit d) { + BigNumLength i; + BigNumProduct c = 0; + BigNumDigit save; + + if (d == BN_ZERO) { + return BN_NOCARRY; + } + + if (d == BN_ONE) { + return BnnAdd(pp, pl, mm, ml, BN_NOCARRY); + } + + for (i = 0; i < ml; ++i) { + BigNumDigit Lm; + BigNumDigit Hm; + BigNumDigit Ld; + BigNumDigit Hd; + BigNumDigit X0; + BigNumDigit X1; + BigNumDigit X2; + BigNumDigit X3; + + Ld = LOW(d); + Hd = HIGH(d); + Lm = LOW(mm[i]); + Hm = HIGH(mm[i]); + X0 = Ld * Lm; + X1 = Ld * Hm; + X2 = Hd * Lm; + X3 = Hd * Hm; + + UPDATE_S(c, X0, X3); + UPDATE_S(c, L2H(X1), X3); + UPDATE_S(c, L2H(X2), X3); + UPDATE_S(c, *pp, X3); + + --pl; + *(pp++) = (BigNumDigit)c; + c = X3 + HIGH(X1) + HIGH(X2); + } + + if (pl == 0) { + return BN_NOCARRY; + } + + save = *pp; + c += save; + *pp = (BigNumDigit)c; + + if (c >= save) { + return BN_NOCARRY; + } + + ++pp; + --pl; + + while (pl != 0 && (++(*pp++)) == 0) { + pl--; + } + + return (pl != 0) ? BN_NOCARRY : BN_CARRY; +} + +/* + * Division + */ + +/** @cond */ +/* xh:xl -= yh:yl */ + +#define SUB(xh, xl, yh, yl) \ + xh -= (BigNumDigit)(yh + (BigNumDigit)(yl > xl)); \ + xl -= yl; +/** @endcond */ + +/** + * BnnDivideDigit. + * Performs the quotient: N div d => Q. + * Returns R = N mod d. + * @pre Assumes leading digit of N < d, and d > 0. + * @param [in, out] qq BigNum + * @param [in, out] nn BigNum + * @param [in] nl BigNumLength + * @param [in] d BigNumDigit + * @return BigNumDigit + */ +BigNumDigit +BnnDivideDigit(BigNum qq, BigNum nn, BigNumLength nl, BigNumDigit d) { + BigNumLength k; + BigNumLength orig_nl; + BigNumDigit rh; /* Two halves of current remainder */ + BigNumDigit rl; /* Correspond to quad above */ + BigNumDigit ph; + BigNumDigit pl; /* product of c and qa */ + BigNumDigit ch; + BigNumDigit cl; + BigNumDigit prev_qq; + + /* + * Normalize divisor + */ + + k = BnnNumLeadingZeroBitsInDigit(d); + + if (k != 0) { +#if defined(__GNUC__) && (__GNUC__ >= 9) +#pragma GCC diagnostic push +#pragma GCC diagnostic ignored "-Warray-bounds" +#endif /* __GNUC__ >= 9 */ + prev_qq = qq[-1]; +#if defined(__GNUC__) && (__GNUC__ >= 9) +#pragma GCC diagnostic pop +#endif /* __GNUC__ >= 9 */ + orig_nl = nl; + if (k < (BigNumLength)BN_DIGIT_SIZE) { + /* + * k should always be < BN_DIGIT_SIZE + * when calling this low level function. + */ + d <<= k; + } + (void)BnnShiftLeft(nn, nl, k); + } else { + prev_qq = 0; + orig_nl = 0; + } + + nn += nl; + nl--; + qq += nl; + + ch = HIGH(d); + cl = LOW(d); + + if (ch == 0) { + /* + * At this point ch can't be == 0; d has been shifted by k + * (the number of leading 0). The test is there to make happy + * static code analysers but is never executed. + */ + return (BigNumDigit)0; + } + + rl = *(--nn); + + while (nl-- != 0) { + BigNumDigit qa; /* Current appr. to quotient */ + + rh = rl; + rl = *(--nn); + qa = rh / ch; /* appr. quotient */ + + /* + * Compute ph, pl + */ + + pl = cl * qa; + ph = ch * qa; + ph += HIGH(pl); + pl = L2H(pl); + + /* + * While ph:pl > rh:rl, decrement qa, adjust qh:ql + */ + + while ((ph > rh) || ((ph == rh) && (pl > rl))) { + qa--; + SUB(ph, pl, ch, L2H(cl)); + } + + SUB(rh, rl, ph, pl); + + /* + * Top half of quotient is correct; save it + */ + + *(--qq) = L2H(qa); + qa = (L2H(rh) | HIGH(rl)) / ch; + + /* + * Approx low half of q. Compute ph, pl, again + */ + + pl = cl * qa; + ph = ch * qa; + ph += HIGH(pl); + pl = LOW(pl) | L2H(LOW(ph)); + ph = HIGH(ph); + + /* + * While ph:pl > rh:rl, decrement qa, adjust qh:ql + */ + + while ((ph > rh) || ((ph == rh) && (pl > rl))) { + qa--; + SUB(ph, pl, 0, d); + } + + /* + * Subtract ph:pl from rh:rl; we know rh will be 0 + */ + + rl -= pl; + *qq |= qa; + } + + /* + * Denormalize dividend + */ + + if (k != 0) { + if ((qq > nn) && (qq < &nn[orig_nl])) { + /* + * Overlap between qq and nn. Care of *qq! + */ + orig_nl = (BigNumLength)(qq - nn); + (void)BnnShiftRight(nn, orig_nl, k); + nn[orig_nl - 1] = prev_qq; + } else if (qq == nn) { + (void)BnnShiftRight(&nn[orig_nl-1], (BigNumLength)1, k); + } else { + (void)BnnShiftRight(nn, orig_nl, k); + } + } + return rl >> k; +} + +/** + * BnnIsZero. + * Returns BN_TRUE iff N = 0 + * @param [in] nn BigNum + * @param [in] nl BigNumLength + * @return BigNumBool + */ +BigNumBool +BnnIsZero(const BigNum nn, BigNumLength nl) { + if ((BnnNumDigits(nn, nl) == (BigNumLength)1) + && (nl == 0 || BnnIsDigitZero(*nn) != BN_FALSE)) { + return BN_TRUE; + } else { + return BN_FALSE; + } +} + +/** + * BnnMultiplyDigit. + * Performs the product: + * - Q = P + M * N + * - BB = BBase(P) + * - Q mod BB => P + * - Q div BB => CarryOut + * + * Returns the CarryOut. + * + * @pre + * - Size(P) >= Size(M) + Size(N), + * - Size(M) >= Size(N). + * @param [in, out] pp BigNum + * @param [in] pl BigNumLength + * @param [in] mm BigNum + * @param [in] ml BigNumLength + * @param [in] nn BigNum + * @param [in] nl BigNumLength + * @return BigNumCarry + */ +BigNumCarry +BnnMultiply(BigNum pp, + BigNumLength pl, + const BigNum mm, + BigNumLength ml, + const BigNum nn, + BigNumLength nl) { + BigNumLength i; + BigNumCarry c = BN_NOCARRY; + + /* + * Multiply one digit at a time + */ + + for (i = 0; i < nl; ++i) { + if (BnnMultiplyDigit(&pp[i], pl--, mm, ml, nn[i]) == BN_CARRY) { + c = BN_CARRY; + } + } + + return c; +} + +/** @cond */ +#define BNN_COMPARE_DIGITS(d1, d2) (d1 == d2) +/** @endcond */ + +/** + * BnnDivideHelper + * + * In-place division. + * + * Input (N has been EXTENDED by 1 PLACE; D is normalized): + * ~~~{.unparsed} + * +-----------------------------------------------+----+ + * | N EXT| + * +-----------------------------------------------+----+ + * + * +-------------------------------+ + * | D 1| + * +-------------------------------+ + * ~~~ + * Output (in place of N): + * ~~~{.unparsed} + * +-------------------------------+---------------+----+ + * | R | Q | + * +-------------------------------+---------------+----+ + * ~~~ + * @pre + * - N > D + * - Size(N) > Size(D) + * - last digit of N < last digit of D + * - D is normalized (Base/2 <= last digit of D < Base) + * @param [in, out] nn BigNum + * @param [in] nl BigNumLength + * @param [in, out] dd BigNum + * @param [in] dl BigNumLength + */ +static void +BnnDivideHelper(BigNum nn, BigNumLength nl, BigNum dd, BigNumLength dl) { + BigNumDigit DDigit; + BigNumDigit BaseMinus1; + BigNumDigit QApp; + BigNumLength ni; + + /* + * Initialize constants + */ + + /* + * BaseMinus1 = BN_COMPLEMENT; + * -> + * BnnSetDigit(&BaseMinus1, BN_ZERO); + * BnnComplement(&BaseMinus1, (BigNumLength)1); + */ + + BaseMinus1 = BN_COMPLEMENT; + + /* + * Save the most significant digit of D + */ + + DDigit = BN_ZERO; + BnnAssign(&DDigit, dd + dl - 1, (BigNumLength)1); + + /* + * Replace D by Base - D + */ + + BnnComplement(dd, dl); + (void)BnnAddCarry(dd, dl, BN_CARRY); + + /* + * For each digit of the divisor, from most significant to least: + */ + + QApp = BN_ZERO; + nl += 1; + ni = nl - dl; + + while (ni != 0) { + /* + * Compute the approximate quotient + */ + + ni--; + nl--; + + /* + * If first digits of numerator and denominator are the same, + */ + + if (BNN_COMPARE_DIGITS(*(nn + nl), DDigit)) { + /* + * Use "Base - 1" for the approximate quotient + */ + BnnAssign(&QApp, &BaseMinus1, (BigNumLength)1); + } else { + /* + * Divide the first 2 digits of N by the + * first digit of D + */ + (void)BnnDivideDigit(&QApp, + nn + nl - 1, + (BigNumLength)2, + DDigit); + } + + /* + * Compute the remainder + */ + + (void)BnnMultiplyDigit(nn + ni, dl + 1, dd, dl, QApp); + + /* + * Correct the approximate quotient, in case it was too large + */ + + while (!BNN_COMPARE_DIGITS(*(nn + nl), QApp)) { + /* + * Subtract D from N + */ + + (void)BnnSubtract(nn+ni, dl + 1, dd, dl, BN_CARRY); + + /* + * Q -= 1 + */ + + (void)BnnSubtractBorrow(&QApp, + (BigNumLength)1, + BN_NOCARRY); + } + } + + /* + * Restore original D + */ + + BnnComplement(dd, dl); + (void)BnnAddCarry(dd, dl, BN_CARRY); +} + +/** + * BnnDivide. + * Performs the quotient: + * - N div D => high-order bits of N, starting at N[dl] + * - N mod D => low-order dl bits of N + * @pre + * - Size(N) > Size(D), + * - last digit of N < last digit of D (if N > D). + * @param [in, out] nn BigNum + * @param [in] nl BigNumLength + * @param [in, out] dd BigNum + * @param [in] dl BigNumLength + */ +void +BnnDivide(BigNum nn, BigNumLength nl, BigNum dd, BigNumLength dl) { + BigNumLength nshift; + + /* + * Take care of easy cases first + */ + + switch (BnnCompare(nn, nl, dd, dl)) { + case BN_LT: /* n < d */ + /* nop */ /* N => R */ + BnnSetToZero(nn + dl, nl - dl); /* 0 => Q */ + return; + case BN_EQ: /* n == d */ + BnnSetToZero(nn, nl); /* 0 => R */ + BnnSetDigit(nn + nl - 1, BN_ONE); /* 1 => Q */ + return; + case BN_GT: /* n > d */ + /* + * If divisor is just 1 digit, use a special divide + */ + + if (dl == (BigNumLength)1) { + /* + * note: nn+1 = nn+dl + */ + + *nn = BnnDivideDigit(nn + 1, nn, nl, *dd); + + /* + * Otherwise, divide one digit at a time + */ + } else { + /* + * Normalize + */ + + nshift = BnnNumLeadingZeroBitsInDigit(*(dd + dl - 1)); + (void)BnnShiftLeft(dd, dl, nshift); + (void)BnnShiftLeft(nn, nl, nshift); + + /* + * Divide + */ + + BnnDivideHelper(nn, nl - 1, dd, dl); + + /* + * Unnormalize + */ + + (void)BnnShiftRight(dd, dl, nshift); + (void)BnnShiftRight(nn, dl, nshift); + + /* + * note: unnormalize N <=> unnormalize R (with R < D) + */ + } + } +} + +/** + * BnnCompare. + * Compare two bignums M and N and returns: + * - BN_GT iff M > N + * - BN_EQ iff N = N + * - BN_LT iff N < N + * @param [in] mm BigNum + * @param [in] ml BigNumLength + * @param [in] nn BigNum + * @param [in] nl BigNumLength + * @return BigNumCmp + */ +BigNumCmp +BnnCompare(const BigNum mm, + BigNumLength ml, + const BigNum nn, + BigNumLength nl) { + ml = BnnNumDigits(mm, ml); + nl = BnnNumDigits(nn, nl); + + if (ml != nl) { + return (ml > nl) ? BN_GT : BN_LT; + } else { + int d; + + /* + * loop starting from most significant digit + */ + + for (d = (int)(nl - 1); d >= 0; --d) { + if (mm[d] != nn[d]) { + return (mm[d] > nn[d]) ? BN_GT : BN_LT; + } + } + + return BN_EQ; + } +} diff --git a/modules/bigz/native/bign.h b/modules/bigz/native/bign.h new file mode 100644 index 00000000..6e907237 --- /dev/null +++ b/modules/bigz/native/bign.h @@ -0,0 +1,243 @@ +/* + * Simplified BSD License + * + * Copyright (c) 1988-1989, Digital Equipment Corporation & INRIA. + * Copyright (c) 1992-2023, Eligis + * All rights reserved. + * + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions are + * met: + * + * o Redistributions of source code must retain the above copyright + * notice, this list of conditions and the following disclaimer. + * o Redistributions in binary form must reproduce the above copyright + * notice, this list of conditions and the following disclaimer in the + * documentation and/or other materials provided with the distribution. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS + * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT + * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR + * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT + * HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, + * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT + * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, + * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY + * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT + * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE + * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + */ + +/** + * @file bign.h + * @brief Types and structures for clients of BigN. + * Bignum internal representation + * ~~~{.unparsed} + * <--------------------------- nl ----------------------------> + * | Least Most | + * |Significant| | | |Significant| + * |BigNumDigit| | | |BigNumDigit| + * |___________|___________|___________|___________|___________| + * ^ (sometimes + * | is zero) + * nn + * ~~~ + * @version 2.1.0 + * @copyright Digital Equipment Corporation & INRIA, 1988-1989. + * @copyright Eligis, 1992-2023 + * @authors B. Serpette + * @authors J. Vuillemin + * @authors J.C. Hervé + * @authors C. Jullien + * $Revision: 392 $ + * $Date: 2023-02-04 06:21:06 +0000 (Sat, 04 Feb 2023) $ + */ + +#if !defined(__BIGN_H) +#define __BIGN_H + +#if defined(HAVE_CONFIG_H) && !defined(__CONFIG_H) +#include "config.h" +#define __CONFIG_H +#endif + +#if defined(_WIN32) && !defined(HAVE_STDINT_H) +#define HAVE_STDINT_H +#endif + +#if defined(_WIN64) +#if !defined(SIZEOF_VOID_P) +#define SIZEOF_VOID_P 8 +#endif +#if !defined(SIZEOF_LONG) +#define SIZEOF_LONG 8 +#endif +#endif /* _WIN64 */ + +#if defined(HAVE_STDINT_H) +#include +#endif + +#if defined(__cplusplus) && !defined(CPP_MODULE) +extern "C" { +#endif + +/** @cond */ +#if defined(BN_EXPERIMENTAL_128BIT) +#define BN_NUM_DIGIT_TYPE +typedef __uint128_t BigNumDigit; +#endif + +/* + * Internal digit type. + */ + +#if !defined(BN_NUM_DIGIT_TYPE) +#define BN_NUM_DIGIT_TYPE +#if (defined(HAVE_STDINT_H) && defined(SIZEOF_VOID_P) && (SIZEOF_VOID_P >= 8)) +typedef uint64_t BigNumDigit; +#else +typedef unsigned long BigNumDigit; +#endif +#endif + +#if defined(__GNUC__) && (__GNUC__ >= 3) +#if !defined(BN_CONST_FUNCTION) +#define BN_CONST_FUNCTION __attribute__((const)) +#endif +#if !defined(BN_PURE_FUNCTION) +#define BN_PURE_FUNCTION __attribute__((pure)) +#endif +#endif /* __GNUC__ >= 3 */ + +#if !defined(BN_PURE_FUNCTION) +#define BN_PURE_FUNCTION +#endif + +#if !defined(BN_CONST_FUNCTION) +#define BN_CONST_FUNCTION +#endif + +/* + * bignum types: digits, big numbers, carries ... + */ + +typedef BigNumDigit * BigNum; /* A big number is a digit pointer */ +typedef BigNumDigit BigNumProduct; /* The product of two digits */ +#if defined(BN_EXPERIMENTAL_128BITX) +typedef __uint128_t BigNumLength; /* The length of a bignum */ +#else +typedef unsigned int BigNumLength; /* The length of a bignum */ +#endif +/** @endcond */ + +/** + * BigNum boolean enum. + */ +typedef enum { + /** false value (0). */ + BN_FALSE = 0, + /** true value (1). */ + BN_TRUE = 1 +} BigNumBool; + +/** + * Results of compare functions. + */ +typedef enum { + /** Less than value (-1). */ + BN_LT = -1, + /** Equal than == value (0). */ + BN_EQ = 0, + /** Greater than comparison value (1). */ + BN_GT = 1 +} BigNumCmp; + +/** + * Carry pseudo-boolean enum type. + */ +typedef enum { + /** No carry (0). */ + BN_NOCARRY = 0, + /** Carry (1). */ + BN_CARRY = 1 +} BigNumCarry; + +/* + * sizes + * + * BN_BYTE_SIZE: number of bits in a byte + * BN_DIGIT_SIZE: number of bits in a digit of a BigNum + */ + +#if !defined(BN_BYTE_SIZE) +/** + * A byte size is generally 8 on most modern machines + * but may be 9! on old 36bit computers. + */ +#define BN_BYTE_SIZE ((BigNumLength)8) +#endif + +/** + * Digit size. + */ +#define BN_DIGIT_SIZE (sizeof(BigNumDigit) * BN_BYTE_SIZE) + +/* + * some constants + */ + +/** + * Digit ZERO constant (0) + */ +#define BN_ZERO ((BigNumDigit)0) +/** + * Digit ONE constant (1). + */ +#define BN_ONE ((BigNumDigit)1) +/** + * Digit COMPLEMENT constant (~0) + */ +#define BN_COMPLEMENT (~(BigNumDigit)0) + +/* + * functions of bign.c + */ + +extern BigNumCarry BnnAdd(BigNum mm, BigNumLength ml, const BigNum nn, BigNumLength nl, BigNumCarry carryin); +extern BigNumCarry BnnAddCarry(BigNum nn, BigNumLength nl, BigNumCarry carryin); +extern void BnnAndDigits(BigNum n, BigNumDigit d); +extern void BnnAssign(BigNum mm, const BigNum nn, BigNumLength nl); +extern BigNumCmp BnnCompare(const BigNum mm, BigNumLength ml, const BigNum nn, BigNumLength nl) BN_PURE_FUNCTION; +extern BigNumCmp BnnCompareDigits(BigNumDigit d1, BigNumDigit d2) BN_CONST_FUNCTION; +extern void BnnComplement(BigNum nn, BigNumLength nl); +extern void BnnComplement2(BigNum nn, BigNumLength nl); +extern void BnnDivide(BigNum nn, BigNumLength nl, BigNum dd, BigNumLength dl); +extern BigNumDigit BnnDivideDigit(BigNum qq, BigNum nn, BigNumLength nl, BigNumDigit d); +extern BigNumDigit BnnGetDigit(const BigNum nn) BN_PURE_FUNCTION; +extern BigNumBool BnnIsPower2(const BigNum nn, BigNumLength nl) BN_PURE_FUNCTION; +extern BigNumBool BnnIsDigitEven(BigNumDigit d) BN_CONST_FUNCTION; +extern BigNumBool BnnIsDigitOdd(BigNumDigit d) BN_CONST_FUNCTION; +extern BigNumBool BnnIsDigitNormalized(BigNumDigit d) BN_CONST_FUNCTION; +extern BigNumBool BnnIsDigitZero(BigNumDigit d) BN_CONST_FUNCTION; +extern BigNumBool BnnIsZero(const BigNum nn, BigNumLength nl) BN_PURE_FUNCTION; +extern BigNumCarry BnnMultiply(BigNum pp, BigNumLength pl, const BigNum mm, BigNumLength ml, const BigNum nn, BigNumLength nl); +extern BigNumCarry BnnMultiplyDigit(BigNum pp, BigNumLength pl, const BigNum mm, BigNumLength ml, BigNumDigit d); +extern BigNumLength BnnNumDigits(const BigNum nn, BigNumLength nl) BN_PURE_FUNCTION; +extern BigNumLength BnnNumLength(const BigNum nn, BigNumLength nl) BN_PURE_FUNCTION; +extern BigNumLength BnnNumCount(const BigNum nn, BigNumLength nl) BN_PURE_FUNCTION; +extern BigNumLength BnnNumLeadingZeroBitsInDigit(BigNumDigit d) BN_CONST_FUNCTION; +extern void BnnOrDigits(BigNum n, BigNumDigit d); +extern void BnnSetDigit(BigNum nn, BigNumDigit d); +extern void BnnSetToZero(BigNum nn, BigNumLength nl); +extern BigNumDigit BnnShiftLeft(BigNum mm, BigNumLength ml, BigNumLength nbits); +extern BigNumDigit BnnShiftRight(BigNum mm, BigNumLength ml, BigNumLength nbits); +extern BigNumCarry BnnSubtract(BigNum mm, BigNumLength ml, const BigNum nn, BigNumLength nl, BigNumCarry carryin); +extern BigNumCarry BnnSubtractBorrow(BigNum nn, BigNumLength nl, BigNumCarry carryin); +extern void BnnXorDigits(BigNum n, BigNumDigit d); + +#if defined(__cplusplus) && !defined(CPP_MODULE) +} +#endif + +#endif /* __BIGN_H */ diff --git a/modules/bigz/native/bigq.c b/modules/bigz/native/bigq.c new file mode 100644 index 00000000..dd5a4cfe --- /dev/null +++ b/modules/bigz/native/bigq.c @@ -0,0 +1,1242 @@ +/* + * Simplified BSD License + * + * Copyright (c) 1992-2023, Eligis + * All rights reserved. + * + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions are + * met: + * + * o Redistributions of source code must retain the above copyright + * notice, this list of conditions and the following disclaimer. + * o Redistributions in binary form must reproduce the above copyright + * notice, this list of conditions and the following disclaimer in the + * documentation and/or other materials provided with the distribution. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS + * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT + * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR + * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT + * HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, + * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT + * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, + * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY + * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT + * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE + * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + */ + +/** + * @file bigq.c + * @brief provides an implementation of "unlimited-precision" + * arithmetic for signed rational. + * + * Several conventions are used in the commentary: + * - A "BigQ" is the name for an arbitrary-precision signed rational. + * - Capital letters (e.g., "Q") are used to refer to the value of BigQs. + * + * Rational numbers are stored in their canonical form a/b where a and + * b are coprime integers (i.e. the only positive integer that + * divides both of them is 1), and b > 0. + * + * @note If any BigQ parameter is passed as BZNULL, function returns BZNULL + * which should be considered as an error. + * @version 2.1.0 + * @copyright Eligis, 1992-2023 + * @author C. Jullien + * $Revision: 394 $ + * $Date: 2023-02-04 15:29:09 +0000 (Sat, 04 Feb 2023) $ + */ + +#include +#include +#include + +#if !defined(__BIGQ_H) +#include "./bigq.h" +#endif + +/** @cond */ +typedef enum { + BQ_COPY, + BQ_SET +} BqCreateMode; +/** @endcond */ + +static BigQ BqCanonicalize(BigQ q); +static BigQ BqCreateInternal(const BigZ n, const BigZ d, BqCreateMode mode); + +/** + * BqCreateInternal. Internally create new BigQ using two BigZ for + * numerator and denominator. + * @param [in] n numerator. + * @param [in] d denominator. + * @param [in] mode creation mode. BQ_SET uses n and d to build the new BigZ. + * BZ_COPY makes a copy of n and d. + * @return a new BigQ. + * @pre + * - n in Z + * - d in (N \ {0}) + */ +static BigQ +BqCreateInternal(const BigZ n, const BigZ d, BqCreateMode mode) { + BigQ q; + BigZ cn; + BigZ cd; + + if ((n == BZNULL) || (d == BZNULL)) { + return BQNULL; + } + +#if !defined(BQ_NEGATIVE_DENOMINATOR) + /* + * By default, when BQ_NEGATIVE_DENOMINATOR is *not* defined, + * bigq module only accepts strictly positive denominator. + */ + if (BzGetSign(d) != BZ_PLUS) { + if (mode == BQ_SET) { + BzFree(d); + BzFree(n); + } + return BQNULL; + } +#else + /* + * Allow denominator to be negative. + */ + if (BzGetSign(d) == BZ_ZERO) { + if (mode == BQ_SET) { + BzFree(d); + BzFree(n); + } + return BQNULL; + } +#endif + + if ((q = (BigQ)BqAlloc()) == 0) { + return BQNULL; + } + + BqSetNumerator(q, BZNULL); + BqSetDenominator(q, BZNULL); + + if (BzGetSign(n) == BZ_ZERO) { + BigZ zero; + BigZ one; + + if (mode == BQ_SET) { + if (BzToInteger(d) == (BzInt)1) { + /* + * n=0, d=1. Already normalised as 0/1. + */ + BqSetNumerator(q, n); + BqSetDenominator(q, d); + return q; + } + /* + * Only free denominator as numerator is already 0. + */ + BzFree(d); + zero = n; + } else if ((zero = BzFromInteger((BzInt)0)) == BZNULL) { + BqFree(q); + return BQNULL; + } + + if ((one = BzFromInteger((BzInt)1)) == BZNULL) { + BzFree(zero); + BqFree(q); + return BQNULL; + } + + BqSetNumerator(q, zero); + BqSetDenominator(q, one); + return q; + } else if (mode == BQ_COPY) { + if ((cn = BzCopy(n)) == BZNULL) { + BqFree(q); + return BQNULL; + } + + if ((cd = BzAbs(d)) == BZNULL) { + BzFree(cn); + BqFree(q); + return BQNULL; + } + + if (BzGetSign(n) != BzGetSign(d)) { + BzSetSign(cn, BZ_MINUS); + } + } else { + cn = n; + cd = d; + + if (BzGetSign(n) != BzGetSign(d)) { + BzSetSign(cn, BZ_MINUS); + } + + BzSetSign(cd, BZ_PLUS); + } + + if (BzGetSign(n) != BzGetSign(d)) { + BzSetSign(cn, BZ_MINUS); + } else { + BzSetSign(cn, BZ_PLUS); + } + + BqSetNumerator(q, cn); + BqSetDenominator(q, cd); + + if (BzLength(cd) != (BigNumLength)1) { + q = BqCanonicalize(q); + } + + return q; +} + +/** + * BqCanonicalize q which is physically modified. + * This function is only called by BqCreateInternal. + * @param [in,out] q BigQ + * @return normalized rational. + * @pre q != BZNULL. + */ +static BigQ +BqCanonicalize(BigQ q) { + const BigZ n = BqGetNumerator(q); + const BigZ d = BqGetDenominator(q); + BigZ gcd; + + if ((gcd = BzGcd(n, d)) == BZNULL) { + BqDelete(q); + return BQNULL; + } + + if (BzToInteger(gcd) != (BzInt)1) { + BigZ nn; + BigZ nd; + + if ((nn = BzDiv(n, gcd)) == BZNULL) { + BzFree(gcd); + BqDelete(q); + return BQNULL; + } + + if ((nd = BzDiv(d, gcd)) == BZNULL) { + BzFree(nn); + BzFree(gcd); + BqDelete(q); + return BQNULL; + } + + BzFree(d); + BzFree(n); + + BqSetNumerator(q, nn); + BqSetDenominator(q, nd); + } + + BzFree(gcd); + + return q; +} + +/* + * Public interface + */ + +/** + * BqCreate a new rational number. + * @param [in] n const BigZ numerator. + * @param [in] d const BigZ + * @return BigQ + */ +BigQ +BqCreate(const BigZ n, const BigZ d) { + return BqCreateInternal(n, d, BQ_COPY); +} + +/** + * BqDelete q. It does nothing if q == BQNULL. + * @param [in] q BigQ + */ +void +BqDelete(BigQ q) { + if (q == BQNULL) { + return; + } else { + const BigZ n = BqGetNumerator(q); + const BigZ d = BqGetDenominator(q); + + if (n != BZNULL) { + BzFree(n); + } + if (d != BZNULL) { + BzFree(d); + } + BqFree(q); + } +} + +/** + * BqAdd. + * Create a new canonicalized BigQ: a + b. + * @param [in] a left hand side BigQ + * @param [in] b right hand side BigQ + * @return BigQ. If either a or b is BQNULL, returns BQNULL. + */ +BigQ +BqAdd(BigQ a, BigQ b) { + if (a == BQNULL || b == BQNULL) { + return BQNULL; + } else { + const BigZ an = BqGetNumerator(a); + const BigZ ad = BqGetDenominator(a); + const BigZ bn = BqGetNumerator(b); + const BigZ bd = BqGetDenominator(b); + BigZ n; + BigZ d; + + if (BzCompare(ad, bd) == BZ_EQ) { + /* + * Easy, same denominator. Only add numerators. + */ + if ((n = BzAdd(an, bn)) == BZNULL) { + return BQNULL; + } + + if ((d = BzCopy(ad)) == BZNULL) { + BzFree(n); + return BQNULL; + } + + return BqCreateInternal(n, d, BQ_SET); + } else { + BigZ tmp1; + BigZ tmp2; + + if ((tmp1 = BzMultiply(an, bd)) == BZNULL) { + return BQNULL; + } + + if ((tmp2 = BzMultiply(ad, bn)) == BZNULL) { + BzFree(tmp1); + return BQNULL; + } + + n = BzAdd(tmp1, tmp2); + + BzFree(tmp1); + BzFree(tmp2); + + if (n == BZNULL) { + return BQNULL; + } + + if ((d = BzMultiply(ad, bd)) == BZNULL) { + BzFree(n); + return BQNULL; + } + + return BqCreateInternal(n, d, BQ_SET); + } + } +} + +/** + * BqSubtract. + * Create a new canonicalized BigQ: a - b. + * @param [in] a left hand side BigQ + * @param [in] b right hand side BigQ + * @return BigQ. If either a or b is BQNULL, returns BQNULL. + */ +BigQ +BqSubtract(BigQ a, BigQ b) { + if (a == BQNULL || b == BQNULL) { + return BQNULL; + } else { + const BigZ an = BqGetNumerator(a); + const BigZ ad = BqGetDenominator(a); + const BigZ bn = BqGetNumerator(b); + const BigZ bd = BqGetDenominator(b); + BigZ n; + BigZ d; + + if (BzCompare(ad, bd) == BZ_EQ) { + /* + * Easy, same denominator. Only subtract numerators. + */ + if ((n = BzSubtract(an, bn)) == BZNULL) { + return BQNULL; + } + + if ((d = BzCopy(ad)) == BZNULL) { + BzFree(n); + return BQNULL; + } + + return BqCreateInternal(n, d, BQ_SET); + } else { + BigZ tmp1; + BigZ tmp2; + + if ((tmp1 = BzMultiply(an, bd)) == BZNULL) { + return BQNULL; + } + + if ((tmp2 = BzMultiply(ad, bn)) == BZNULL) { + BzFree(tmp1); + return BQNULL; + } + + if ((n = BzSubtract(tmp1, tmp2)) == BZNULL) { + BzFree(tmp2); + BzFree(tmp1); + return BQNULL; + } + + d = BzMultiply(ad, bd); + + BzFree(tmp2); + BzFree(tmp1); + + if (d == BZNULL) { + BzFree(n); + return BQNULL; + } + + return BqCreateInternal(n, d, BQ_SET); + } + } +} + +/** + * BqMultiply. + * Create a new canonicalized BigQ: a * b. + * @param [in] a left hand side BigQ + * @param [in] b right hand side BigQ + * @return BigQ. If either a or b is BQNULL, returns BQNULL. + */ +BigQ +BqMultiply(BigQ a, BigQ b) { + if (a == BQNULL || b == BQNULL) { + return BQNULL; + } else { + const BigZ an = BqGetNumerator(a); + const BigZ ad = BqGetDenominator(a); + const BigZ bn = BqGetNumerator(b); + const BigZ bd = BqGetDenominator(b); + BigZ n; + BigZ d; + + if ((n = BzMultiply(an, bn)) == BZNULL) { + return BQNULL; + } + + if ((d = BzMultiply(ad, bd)) == BZNULL) { + BzFree(n); + return BQNULL; + } + + return BqCreateInternal(n, d, BQ_SET); + } +} + +/** + * BqDiv. + * Create a new canonicalized BigQ: a / b. + * @param [in] a left hand side BigQ + * @param [in] b right hand side BigQ + * @return BigQ. If either a or b is BQNULL, returns BQNULL. + */ +BigQ +BqDiv(BigQ a, BigQ b) { + if (a == BQNULL || b == BQNULL) { + return BQNULL; + } else { + const BigZ an = BqGetNumerator(a); + const BigZ ad = BqGetDenominator(a); + const BigZ bn = BqGetNumerator(b); + const BigZ bd = BqGetDenominator(b); + BigZ n; + BigZ d; + + if ((n = BzMultiply(an, bd)) == BZNULL) { + return BQNULL; + } + + if ((d = BzMultiply(ad, bn)) == BZNULL) { + BzFree(n); + return BQNULL; + } + + if (BzGetSign(n) != BzGetSign(d)) { + BzSetSign(n, BZ_MINUS); + } else { + BzSetSign(n, BZ_PLUS); + } + + BzSetSign(d, BZ_PLUS); + + return BqCreateInternal(n, d, BQ_SET); + } +} + +/** + * BqCompare. + * Compare a and b. It returns BQ_EQ (0) if a == b, BQ_LZ (-1) if a < b and + * BQ_GT (1) if a > b. When a and b have a different denominator, two BigZ + * are internally created which may fail. In this case BQ_ERR (100) is returned. + * This may fool functions like qsort waiting for negative, 0 or positive value. + * @param [in] a left hand side BigQ + * @param [in] b right hand side BigQ + * @return BqCmp. If either a or b is BQNULL, returns BQ_ERR. + */ +BqCmp +BqCompare(BigQ a, BigQ b) { + if (a == BQNULL || b == BQNULL) { + return BQ_ERR; + } else { + const BigZ an = BqGetNumerator(a); + const BigZ ad = BqGetDenominator(a); + const BigZ bn = BqGetNumerator(b); + const BigZ bd = BqGetDenominator(b); + BigZ tmp1; + BigZ tmp2; + BzCmp cmp; + + if (BzGetSign(an) != BzGetSign(bn)) { + /* + * Sign differs, easy case! + */ + if (BzGetSign(an) == BZ_MINUS) { + return BQ_LT; + } else if (BzGetSign(an) == BZ_ZERO + && BzGetSign(bn) == BZ_PLUS) { + return BQ_LT; + } else { + return BQ_GT; + } + } else if (BzCompare(ad, bd) == BZ_EQ) { + /* + * Easy, same denominator. Only compare numerators. + * This works also if a == b == 0 because a and b + * are canonicalized and their denominator is 1. + */ + cmp = BzCompare(an, bn); + } else { + if ((tmp1 = BzMultiply(an, bd)) == BZNULL) { + return BQ_ERR; + } + + if ((tmp2 = BzMultiply(ad, bn)) == BZNULL) { + BzFree(tmp1); + return BQ_ERR; + } + + cmp = BzCompare(tmp1, tmp2); + + BzFree(tmp2); + BzFree(tmp1); + } + + switch (cmp) { + case BZ_LT: + return BQ_LT; + case BZ_GT: + return BQ_GT; + case BZ_EQ: + default: + return BQ_EQ; + } + } +} + +/** + * BqNegate. + * Create a new canonicalized BigQ: -a. + * @param [in] a BigQ + * @return BigQ. If a is BQNULL, returns BQNULL. + */ +BigQ +BqNegate(BigQ a) { + if (a == BQNULL) { + return BQNULL; + } else { + const BigZ an = BqGetNumerator(a); + const BigZ ad = BqGetDenominator(a); + BigQ res = BqCreateInternal(an, ad, BQ_COPY); + + if (res == BQNULL) { + return res; + } + + switch (BzGetSign(an)) { + case BZ_MINUS: + BzSetSign(BqGetNumerator(res), BZ_PLUS); + return res; + case BZ_PLUS: + BzSetSign(BqGetNumerator(res), BZ_MINUS); + return res; + case BZ_ZERO: + default: + return res; + } + } +} + +/** + * BqAbs. + * Create a new canonicalized BigQ: |a|. + * @param [in] a BigQ + * @return BigQ + */ +BigQ +BqAbs(BigQ a) { + if (a == BQNULL) { + return BQNULL; + } else { + const BigZ an = BqGetNumerator(a); + const BigZ ad = BqGetDenominator(a); + BigQ res = BqCreateInternal(an, ad, BQ_COPY); + + if (res == BQNULL) { + return res; + } + + if (BzGetSign(BqGetNumerator(res)) == BZ_MINUS) { + /* + * Force positive sign only when a res is negative. + */ + BzSetSign(BqGetNumerator(res), BZ_PLUS); + } + + return res; + } +} + +/** + * BqInverse. + * Create a new canonicalized BigQ: an/ad => ad/an. The result sign is the sign + * of an. + * @param [in] a BigQ + * @return BigQ. If a is BQNULL, returns BQNULL. + */ +BigQ +BqInverse(BigQ a) { + if (a == BQNULL) { + return BQNULL; + } else { + const BigZ an = BqGetNumerator(a); + const BigZ ad = BqGetDenominator(a); + + return BqCreateInternal(ad, an, BQ_COPY); + } +} + +/* + * Define QNaN as an array (not a pointer!!!) to let sizeof returns + * the null terminating string length (including '\000'). + */ + +static const char BqNaN[] = "#.QNaN"; /* it includes NUL terminated char */ + +/** + * BqToStringBufferExt. + * Returns a pointer to a string that represents Q number in the specified + * base. Assumes BZ_MIN_BASE <= base <= BZ_MAX_BASE. If optional + * buffer is supplied, len is a pointer of this buffer size. If there + * is enough room to print the number buf is used otherwise function + * returns NULL and len contains the required size. If buf is passed + * as NULL, this string is allocated on the heap, so it must be + * desallocated by the user. + * @param [in] q BigQ + * @param [in] base BigNumDigit + * @param [in] sign flag. When set to BQ_FORCE_SIGN and q != 0, + * a sign +/- is always present. + * @param [out] buf output buffer. When passed NULL, buffer is allocated. + * @param [in, out] buflen on input it contains the maximal buffer length. When + * buffer is too small to represent this BigQ number, function returns + * NULL and buflen is set the the required buffer size. + * @param [out] slen when not NULL, slen is a pointer of size_t which contains + * the actual string length of q. + * @return NUL terminated string + */ +BzChar * +BqToStringBufferExt(BigQ q, + BigNumDigit base, + int sign, + BzChar *buf, + size_t *buflen, + size_t *slen) { + BzChar *n; + BzChar *d; + BzChar *res; + size_t len; + int i; + + if ((buf != (BzChar *)NULL) && (buflen == (size_t *)NULL)) { + /* + * A buffer is passed without its associated length. + */ + return (BzChar *)NULL; + } else if (q == BQNULL) { + len = sizeof(BqNaN); /* works because BqNaN is an array */ + + if (buf != (BzChar *)NULL) { + if (*buflen >= len) { + res = buf; + } else { + res = (BzChar *)NULL; + } + } else { + /* + * contract is to allocate a new string, even + * for #.QNaN error. + */ + res = (BzChar *)BzStringAlloc(len); + } + + if (slen != (size_t *)NULL) { + *slen = len; + } + + if (res != (BzChar *)NULL) { + for (i = 0; BqNaN[i] != '\000'; ++i) { + res[i] = (BzChar)BqNaN[i]; + } + res[i] = (BzChar)'\000'; + } + + return res; + } else if (BzLength(BqGetDenominator(q)) == (BigNumLength)1) { + /* + * n/1 is printed just n + */ + return BzToStringBufferExt(BqGetNumerator(q), + base, + sign, + buf, + buflen, + slen); + } else if (buf != (BzChar *)NULL) { + /* + * Get numerator string. + */ + n = BzToString(BqGetNumerator(q), base, sign); + + if (n == (BzChar *)NULL) { + return n; + } + + /* + * Get denominator string. + */ + d = BzToString(BqGetDenominator(q), base, BZ_DEFAULT_SIGN); + + if (d == (BzChar *)NULL) { + BzFreeString(n); + return d; + } + + /* + * Compute total length + */ + len = BzStrLen(n) + BzStrLen(d) + 1; /* +1 for '/' */ + + if (slen != (size_t *)NULL) { + *slen = len; + } + + res = &buf[0]; + + /* + * catenate n/d in return buffer. + */ + + for (i = 0; n[i] != (BzChar)'\000'; ++i) { + *res++ = n[i]; + } + + *res++ = (BzChar)'/'; + + for (i = 0; d[i] != (BzChar)'\000'; ++i) { + *res++ = d[i]; + } + + *res++ = (BzChar)'\000'; + + BzFreeString(d); + BzFreeString(n); + + return &buf[0]; + } else { + /* + * Get numerator string. + */ + n = BzToString(BqGetNumerator(q), base, sign); + + if (n == (BzChar *)NULL) { + return n; + } + + /* + * Get denominator string. + */ + d = BzToString(BqGetDenominator(q), base, BQ_DEFAULT_SIGN); + + if (d == (BzChar *)NULL) { + BzFreeString(n); + return d; + } + + /* + * Compute total length + */ + len = BzStrLen(n) + BzStrLen(d) + 1; /* +1 for '/' */ + + /* + * Alloc result string. + */ + res = (BzChar *)BzStringAlloc(len + 1); /* +1 for nul */ + + if (slen != (size_t *)NULL) { + *slen = len; + } + + /* + * catenate n/d in return buffer. + */ + + if (res != (BzChar *)NULL) { + int pos = 0; + + for (i = 0; n[i] != (BzChar)'\000'; ++i) { + res[pos++] = n[i]; + } + + res[pos++] = (BzChar)'/'; + + for (i = 0; d[i] != (BzChar)'\000'; ++i) { + res[pos++] = d[i]; + } + + res[pos] = (BzChar)'\000'; + } + + BzFreeString(d); + BzFreeString(n); + + return res; + } +} + +/** + * BqToStringBuffer. + * Returns a pointer to a string that represents Q number in the specified + * base. Assumes BZ_MIN_BASE <= base <= BZ_MAX_BASE. If optional + * buffer is supplied, len is a pointer of this buffer size. If there + * is enough room to print the number buf is used otherwise function + * returns NULL and len contains the required size. If buf is passed + * as NULL, this string is allocated on the heap, so it must be + * desallocated by the user. + * @param [in] q BigQ + * @param [in] base BigNumDigit + * @param [in] sign flag. When set to BQ_FORCE_SIGN and q != 0, + * a sign +/- is always present. + * @param [out] buf output buffer + * @param [in, out] buflen buffer length + * @return NUL terminated string + */ +BzChar * +BqToStringBuffer(BigQ q, + BigNumDigit base, + int sign, + BzChar *buf, + size_t *buflen) { + return BqToStringBufferExt(q, + base, + sign, + buf, + buflen, + (size_t *)0); +} + +/** + * BqToString. + * Convert q to string. The result string is allocated by BzStringAlloc + * and must be free by caller using BzFreeString. + * When q is BQNULL, it returns "#.QNaN". + * @param [in] q BigQ. + * @param [in] sign flag. When set to BQ_FORCE_SIGN and q != 0, + * a sign +/- is always present. + * @return BzChar string. + * @see BqToStringBufferExt for more options. + */ +BzChar * +BqToString(BigQ q, int sign) { + return BqToStringBufferExt(q, + (BigNumDigit)10, + sign, + (BzChar *)NULL, + (size_t *)NULL, + (size_t *)NULL); +} + +/** + * BqFromString. + * Create a new BigQ from string. + * It returns BQNULL is s does not represent a valid rational number or if + * denominator is not a strictly positive number. The new rational is + * returned in its canonicalized form (e.g. "3/6" => 1/2). + * @param [in] s const BzChar + * @param [in] base int in [1 .. 36]. + * @return BigQ. If s is NULL or the empty string "", returns BQNULL. + */ +BigQ +BqFromString(const BzChar *s, int base) { + BigZ n; + BigZ d; + BigQ q; + const BzChar *p; + + if (s == (const BzChar *)NULL || *s == (BzChar)0) { + return BQNULL; + } + + /* + * Throw away any initial space + */ + + while ((*s == (BzChar)' ') + || (*s == (BzChar)'\t') + || (*s == (BzChar)'\n') + || (*s == (BzChar)'\r')) { + s++; + } + + /* + * search for '/' + */ + + p = s; + + if (*p == (BzChar)'+' || *p == (BzChar)'-') { + ++p; + if ((*p == (BzChar)'+') || (*p == (BzChar)'-')) { + /* + * optional + or - must be unique and already skipped. + */ + return BQNULL; + } + } + + if (*p == (BzChar)0) { + return BQNULL; + } + + while (*p != (BzChar)'\000') { + if (*p == (BzChar)'/') { + break; + } else { + ++p; + } + } + + if (*p == (BzChar)'\000') { + /* + * simply an integer in Z (no denominator). + */ + n = BzFromString(s, (BigNumDigit)base, BZ_UNTIL_SPACE); + if (n == BZNULL) { + return BQNULL; + } + + d = BzFromInteger((BzInt)1); + if (d == BZNULL) { + BzFree(n); + return BQNULL; + } + + q = BqCreateInternal(n, d, BQ_SET); + return q; + } else { + ++p; /* skip slash */ + + if ((*p == (BzChar)'+') + || (*p == (BzChar)'-') + || (*p == (BzChar)' ')) { + /* + * Denominator is always positive with no sign and + * must follow without spaces. + */ + return BQNULL; + } + + n = BzFromString(s, (BigNumDigit)base, BZ_UNTIL_SLASH); + if (n == BZNULL) { + return BQNULL; + } + + d = BzFromString(p, (BigNumDigit)base, BZ_UNTIL_SPACE); + if (d == BZNULL) { + BzFree(n); + return BQNULL; + } + + q = BqCreateInternal(n, d, BQ_SET); + return q; + } +} + +/** + * BqFromLongDouble. + * Find rational approximation to given real number (Farey's method). + * @param [in] num long double. + * @param [in] maxd maximal denominator value. + * @return BigQ approximation. + */ +BigQ +BqFromLongDouble(BzLDouble num, BzInt maxd) { + BzInt ln = (BzInt)0; /* lower value = 0/1 */ + BzInt ld = (BzInt)1; + BzInt un = (BzInt)1; /* upper value = 1/0 = oo */ + BzInt ud = (BzInt)0; + BzInt rn = (BzInt)1; + BzInt rd = (BzInt)0; + BigZ n; + BigZ d; + BigQ q; + int sign; + + /* + * See: http://en.wikipedia.org/wiki/Farey_series + * http://wiki.cs.princeton.edu/index.php/Rational.ck + */ + + if (num < (BzLDouble)0.0) { + sign = -1; + num *= (BzLDouble)-1.0; + } else { + sign = 1; + } + + for (;;) { + const BzInt mn = ln + un; + const BzInt md = ld + ud; + + if ((num * (BzLDouble)md) > (BzLDouble)mn) { + if (maxd < md) { + /* + * return upper. + */ + rn = un; + rd = ud; + break; + } else { + /* + * set lower to median and continue + */ + ln = mn; + ld = md; + continue; + } + } else if ((num * (BzLDouble)md) == (BzLDouble)mn) { + if (maxd >= md) { + /* + * return median. + */ + rn = mn; + rd = md; + break; + } else if (ld < ud) { + /* + * return lower. + */ + rn = ln; + rd = ld; + break; + } else { + /* + * return upper. + */ + rn = un; + rd = ud; + break; + } + } else { + if (maxd < md) { + /* + * return lower. + */ + rn = ln; + rd = ld; + break; + } else { + /* + * set lower to median and continue + */ + un = mn; + ud = md; + continue; + } + } + } + + if ((n = BzFromInteger(sign * rn)) == BZNULL) { + return BQNULL; + } + + if ((d = BzFromInteger(rd)) == BZNULL) { + BzFree(n); + return BQNULL; + } + + q = BqCreateInternal(n, d, BQ_SET); + + return q; +} + +/** + * BqFromDouble. + * Find rational approximation to given real number (Farey's method). + * @param [in] num double. + * @param [in] maxd maximal denominator value. + * @return BigQ approximation. + */ +BigQ +BqFromDouble(double num, BzInt maxd) { + return BqFromLongDouble((BzLDouble)num, maxd); +} + +#if defined(isnan) +/* + * isnan exists since C99, otherwise (x != x) should work the same. + * It is said that "NaN values never compare equal to themselves or to + * other NaN values." + */ +#define BqIsNan(x) (x != x) /* isnan(x) */ +#else +#define BqIsNan(x) (x != x) +#endif + +/** + * BqToLongDouble. + * Convert q to a long double approximation. + * @param [in] q BigQ + * @return BigQ. If q == BQNULL and if the implementation supports quiet NaNs + * it retruns NAN (as defined in math.h). Otherwise it returns 0.0. + */ +BzLDouble +BqToLongDouble(BigQ q) { + if (q == BQNULL) { +#if defined(NAN) + return (BzLDouble)NAN; +#else + return ((BzLDouble)0/(BzLDouble)0); +#endif + } else { + BzLDouble res; + + res = BzToLongDouble(BqGetNumerator(q)) + / BzToLongDouble(BqGetDenominator(q)); + + if (BqIsNan(res)) { + /* + * Naive code above does not work because numerator + * and/or denominator is too big to be converted + * to (long?) double, try harder as the expected + * result may be well below max (long?) double. + * + * (defun rational-to-float (x) + * (let ((zdiv (ash 1 52)) + * (zval nil) + * (qtmp nil) + * (ztmp nil) + * (dres nil)) + * (setf zval (floor (numerator x) + * (denominator x))) + * (setf qtmp (* (- x zval) zdiv)) + * (setf ztmp (floor (numerator qtmp) + * (denominator qtmp))) + * (setf qtmp (/ ztmp zdiv)) + * (setf dres (+ (float zval) (float qtmp))) + * dres)) + */ + BigZ zone; + BigZ zdiv; + BigQ qdiv; + BigQ qtmp1; + + zone = BzFromInteger((BzInt)1); + zdiv = BzAsh(zone, 52); /* mantissa is on 53 bits */ + qdiv = BqCreate(zdiv, zone); + + { + BigZ zfloor; + BigQ qfloor; + + zfloor = BzFloor(BqGetNumerator(q), + BqGetDenominator(q)); + qfloor = BqCreate(zfloor, zone); + res = BzToLongDouble(zfloor); + qtmp1 = BqSubtract(q, qfloor); + + BqFree(qfloor); + BzFree(zfloor); + } + + { + BigQ qtmp2; + + qtmp2 = BqMultiply(qtmp1, qdiv); + BqFree(qtmp1); + { + BigZ ztmp; + + ztmp = BzFloor(BqGetNumerator(qtmp2), + BqGetDenominator(qtmp2)); + qtmp1 = BqCreate(ztmp, zdiv); + res += BqToLongDouble(qtmp1); + + BqFree(qtmp1); + BzFree(ztmp); + } + BqFree(qtmp2); + } + + BqFree(qdiv); + BzFree(zdiv); + BzFree(zone); + } + + return res; + } +} + +/** + * BqToDouble. + * Convert q to a double approximation. + * @param [in] q BigQ + * @return BigQ. If q == BQNULL and if the implementation supports quiet NaNs + * it retruns NAN (as defined in math.h). Otherwise it returns 0.0. + */ +double +BqToDouble(BigQ q) { + return (double)BqToLongDouble(q); +} + +/** + * BqHash + * The hash value of n/q is computed as hash(n) + ~hash(q). + * @param [in] q BigQ + * @return BzUint hash value + */ +BzUInt +BqHash(BigQ q) { + const BigZ n = BqGetNumerator(q); + const BigZ d = BqGetDenominator(q); + + return BzHash(n) + ~BzHash(d); +} diff --git a/modules/bigz/native/bigq.h b/modules/bigz/native/bigq.h new file mode 100644 index 00000000..6a647d3b --- /dev/null +++ b/modules/bigz/native/bigq.h @@ -0,0 +1,174 @@ +/* + * Simplified BSD License + * + * Copyright (c) 1992-2023, Eligis + * All rights reserved. + * + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions are + * met: + * + * o Redistributions of source code must retain the above copyright + * notice, this list of conditions and the following disclaimer. + * o Redistributions in binary form must reproduce the above copyright + * notice, this list of conditions and the following disclaimer in the + * documentation and/or other materials provided with the distribution. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS + * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT + * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR + * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT + * HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, + * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT + * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, + * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY + * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT + * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE + * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + */ + +/** + * @file bigq.h + * @brief Types and structures for clients of BigQ. + * @version 2.1.0 + * @copyright Eligis, 1992-2023 + * @author C. Jullien + * $Revision: 394 $ + * $Date: 2023-02-04 15:29:09 +0000 (Sat, 04 Feb 2023) $ + */ + +#if !defined(__BIGQ_H) +#define __BIGQ_H + +#if !defined(__BIGZ_H) +#include "./bigz.h" +#endif + +#if defined(__cplusplus) +extern "C" { +#endif + +/** @cond */ +#define BQ_PURE_FUNCTION BN_PURE_FUNCTION +#define BQ_CONST_FUNCTION BN_CONST_FUNCTION +/** @endcond */ + +/* + * BigQ.h: Types and structures for clients of BigQ + */ + +/** + * BqToString and related functions add sign only when number is negative ('-'). + */ +#define BQ_DEFAULT_SIGN BZ_DEFAULT_SIGN + +/** + * Force BqToString to always add a sign [+/-] when number is not 0. + */ +#define BQ_FORCE_SIGN BZ_FORCE_SIGN + +/** + * BigQ compare result + */ +typedef enum { + /** Less than value (-1). */ + BQ_LT = BN_LT, + /** Equal than == value (0). */ + BQ_EQ = BN_EQ, + /** Greater than comparison value (1). */ + BQ_GT = BN_GT, + /** Rerror */ + BQ_ERR = 100 +} BqCmp; + +/** @cond */ +/** + * BigQ number is a pair of two bignums (numerator / denominator). + * The BigQ sign is the sign of numerator, denominator is always greater than + * 0. + */ +typedef struct { + /** numerator, a signed BigZ. */ + BigZ N; + /** denumerator, a strictly positive unsigned BigZ. */ + BigZ D; +} BigQStruct; + +typedef BigQStruct * __BigQ; + +#if !defined(BQ_RATIONAL_TYPE) +#define BQ_RATIONAL_TYPE +typedef const BigQStruct * BigQ; +#endif +/** @endcond */ + +#if !defined(__EXTERNAL_BIGQ_MEMORY) +/** + * User overloadable macro that gets native BigQ implementation from + * a high level object (for example managed in C++ or Lisp). + */ +#define __toBqObj(q) ((__BigQ)q) +/** + * NULL BigQ. + */ +#define BQNULL ((BigQ)0) +/** + * User overloadable macro called to allocate size bytes to store a BigQ. + */ +#define BqAlloc() malloc(sizeof(BigQStruct)) +/** + * User overloadable macro called to free a BigQ allocated by BqAlloc. + */ +#define BqFree(q) free((void *)q) /* free(__toBqObj(q)) */ +/** + * Get BigQ numerator. + */ +#define BqGetNumerator(q) (__toBqObj(q)->N) +/** + * Get BigQ denominator. + */ +#define BqGetDenominator(q) (__toBqObj(q)->D) +/** + * Set BigQ numerator. + */ +#define BqSetNumerator(q, n) (__toBqObj(q)->N = (n)) +/** + * Set BigQ denominator. + */ +#define BqSetDenominator(q, d) (__toBqObj(q)->D = (d)) +#endif + +/* + * functions of bigq.c + */ + +extern BqCmp BqCompare(BigQ a, BigQ b) BQ_PURE_FUNCTION; +extern BzChar * BqToString(BigQ q, int sign); +extern BzChar * BqToStringBuffer(BigQ q, BigNumDigit base, int sign, /*@null@*/ BzChar *buf, /*@null@*/ size_t *buflen); +extern BzChar * BqToStringBufferExt(BigQ q, BigNumDigit base, int sign, /*@null@*/ BzChar *buf, /*@null@*/ size_t *buflen, /*@null@*/ size_t *slen); +extern BigQ BqCreate(const BigZ n, const BigZ d); +extern BigQ BqFromDouble(double num, BzInt maxd); +extern BigQ BqFromLongDouble(BzLDouble num, BzInt maxd); +extern BigQ BqFromString(const BzChar *s, int base); +extern double BqToDouble(BigQ a); +extern BzLDouble BqToLongDouble(BigQ a); +extern void BqDelete(BigQ a); + +extern BigQ BqAbs(BigQ a); +extern BigQ BqAdd(BigQ a, BigQ b); +extern BigQ BqDiv(BigQ a, BigQ b); +extern BigQ BqInverse(BigQ a); +extern BigQ BqMultiply(BigQ a, BigQ b); +extern BigQ BqNegate(BigQ a); +extern BigQ BqSubtract(BigQ a, BigQ b); +extern BzUInt BqHash(BigQ z); + +#if 0 +extern BzChar * BqToStringBuffer(BigQ q, int sign, BzChar *buf, size_t *len); +#endif + +#if defined(__cplusplus) +} +#endif + +#endif /* __BIGQ_H */ diff --git a/modules/bigz/native/bigz.c b/modules/bigz/native/bigz.c new file mode 100644 index 00000000..b3bf616a --- /dev/null +++ b/modules/bigz/native/bigz.c @@ -0,0 +1,3379 @@ +/* + * Simplified BSD License + * + * Copyright (c) 1988-1989, Digital Equipment Corporation & INRIA. + * Copyright (c) 1992-2023, Eligis + * All rights reserved. + * + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions are + * met: + * + * o Redistributions of source code must retain the above copyright + * notice, this list of conditions and the following disclaimer. + * o Redistributions in binary form must reproduce the above copyright + * notice, this list of conditions and the following disclaimer in the + * documentation and/or other materials provided with the distribution. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS + * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT + * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR + * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT + * HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, + * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT + * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, + * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY + * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT + * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE + * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + */ + +/** + * @file bigz.c + * @brief provides an implementation of "unlimited-precision" + * arithmetic for signed integers. + * + * Several conventions are used in the commentary: + * - A "BigZ" is the name for an arbitrary-precision signed integer. + * - Capital letters (e.g., "Z") are used to refer to the value of BigZs. + * + * The logical operations provide a convenient way to represent an + * infinite vector of bits. Let such a conceptual vector be indexed + * by the non-negative integers. Then bit j is assigned a ``weight'' + * 2**j. Assume that only a finite number of bits are 1's or only a + * finite number of bits are 0's. A vector with only a finite number + * of one-bits is represented as the sum of the weights of the + * one-bits, a positive integer. A vector with only a finite number + * of zero-bits is represented as -1 minus the sum of the weights of + * the zero-bits, a negative integer. + * + * This method of using integers to represent bit-vectors can in turn + * be used to represent sets. Suppose that some (possibly countably + * infinite) universe of discourse for sets is mapped into the + * non-negative integers. Then a set can be represented as a bit + * vector; an element is in the set if the bit whose index corresponds + * to that element is a one-bit. In this way all finite sets can be + * represented (by positive integers), as well as all sets whose + * complements are finite (by negative integers). + * @version 2.1.0 + * @copyright Digital Equipment Corporation & INRIA, 1988-1989. + * @copyright Eligis, 1992-2023 + * @authors B. Serpette + * @authors J. Vuillemin + * @authors J.C. Hervé + * @authors C. Jullien + * $Revision: 394 $ + * $Date: 2023-02-04 15:29:09 +0000 (Sat, 04 Feb 2023) $ + */ + +/** @cond */ +#if !defined(_CRT_SECURE_NO_DEPRECATE) +#define _CRT_SECURE_NO_DEPRECATE 1 +#endif +#if !defined(_CRT_NONSTDC_NO_DEPRECATE) +#define _CRT_NONSTDC_NO_DEPRECATE 1 +#endif +/** @endcond */ + +#include + +#include +#include +#include +#include + +#if !defined(__BIGZ_H) +#include "./bigz.h" +#endif + +/** @cond */ +#define BzMaxInt(a, b) (((a) < (b)) ? (b) : (a)) +#define BzAbsInt(x) (((x) >= 0) ? (x) : -(x)) + +#define BZMAXINT ((BzInt)((~(BzUInt)0) >> 1)) +#define BZMAXUINT (~(BzUInt)0) + +#define BzFreeIf(cond, ptr) if (cond) BzFree(ptr) +/** @endcond */ + +/* + * See ./etc/hextable.c if you need to change BigHexToDigit tables. + */ + +/** @cond */ +#if !defined(OLEBCDIC) +/** + * @brief An array for converting hexadecimal digits to integer values + * using ASCII encoding characters. + * The value at index `i` is the integer value of the hexadecimal + * digit represented by the character with ASCII code `i`. If the + * character is not a valid hexadecimal digit, the value is -1. + * + * This array is used for converting hexadecimal strings to BigNum values. + */ +static const int BigHexToDigit[] = { + -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, + -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, + -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, + 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, -1, -1, -1, -1, -1, -1, + -1, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, + 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, -1, -1, -1, -1, -1, + -1, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, + 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, -1, -1, -1, -1, -1 +}; + +#define CTOI(c) ((((unsigned int)c) < (unsigned int)127) \ + ? BigHexToDigit[(unsigned int)c] \ + : -1) +#else +/** + * @brief An array for converting hexadecimal digits to integer values + * using EBCDIC encoding characters (as most IBM manframes like z/OS, OS/360...) + * The value at index `i` is the integer value of the hexadecimal + * digit represented by the character with ASCII code `i`. If the + * character is not a valid hexadecimal digit, the value is -1. + * + * This array is used for converting hexadecimal strings to BigNum values. + */ +static const int BigHexToDigit[] = { + -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, + -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, + -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, + -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, + -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, + -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, + -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, + -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, + -1, 10, 11, 12, 13, 14, 15, 16, 17, 18, -1, -1, -1, -1, -1, -1, + -1, 19, 20, 21, 22, 23, 24, 25, 26, 27, -1, -1, -1, -1, -1, -1, + -1, -1, 28, 29, 30, 31, 32, 33, 34, 35, -1, -1, -1, -1, -1, -1, + -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, + -1, 10, 11, 12, 13, 14, 15, 16, 17, 18, -1, -1, -1, -1, -1, -1, + -1, 19, 20, 21, 22, 23, 24, 25, 26, 27, -1, -1, -1, -1, -1, -1, + -1, -1, 28, 29, 30, 31, 32, 33, 34, 35, -1, -1, -1, -1, -1, -1, + 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, -1, -1, -1, -1, -1, -1 +}; + +#define CTOI(c) ((((unsigned int)c) < (unsigned int)255) \ + ? BigHexToDigit[(unsigned int)c] \ + : -1) +#endif + +/** + * @brief Check if a character is a white space. + * @param [in] c the character to check. + * @return true if the character is a white space, false otherwise. + */ +#define BZ_ISSPACE(c) \ + ((c == (BzChar)' ') \ + || (c == (BzChar)'\t') \ + || (c == (BzChar)'\n') \ + || (c == (BzChar)'\r')) + +/** + * @brief Check if a character is a sign (+ or -). + * @param [in] c the character to check. + * @return true if the character is a sign, false otherwise. + */ +#define BZ_ISSIGN(c) \ + ((c == (BzChar)'+') \ + || (c == (BzChar)'-')) + +/** + * @brief Check if a character is a slash (/). + * @param [in] c the character to check. + * @return true if the character is a slash, false otherwise. + */ +#define BZ_ISSLASH(c) \ + (c == (BzChar)'/') + +/** @endcond */ + +static BzSign BzGetOppositeSign(const BigZ z); + +#if defined(BZ_DEBUG) +static void BzShowBits(BigNumDigit n); +static void BzShowUnsingnedInt(unsigned int n); + +/** + * @brief Prints the binary representation of a BigNum digit to standard output. + * @param [in] n the BigNum digit to print. + * This function is used for debugging purposes. It prints the binary + * representation of the input digit, with the most significant bit on + * the left and the least significant bit on the right. + */ + +static void +BzShowBits(BigNumDigit n) { + int i; + + for (i = (int)(BN_DIGIT_SIZE - 1); i >= 0; i--) { + if ((n & (BN_ONE << (unsigned int)i)) != 0) { + (void)printf("1"); + } else { + (void)printf("0"); + } + } +} + +/** + * @brief Prints the hexadecimal representation of an unsigned integer + * to standard output. + * + * This function is used for debugging purposes. It prints the + * hexadecimal representation of the input integer, with the most + * significant nibble on the left and the least significant nibble on + * the right. The number of digits printed depends on the size of the + * input integer in bytes. + * + * @param [in] n The unsigned integer to print. + * @note This function handles integers of different sizes, up to 16 bytes. + */ +static void + +static void +BzShowUnsingnedInt(unsigned int n) { + /* + * Deal with different integer sizes. This avoid to build format + * string dynamically which results to potential buffer overflow. + * This also makes splint happy. + */ + + size_t size = sizeof(n); + + switch (size) { + case 2: (void)printf("%04x", n); break; + case 4: (void)printf("%08x", n); break; + case 8: (void)printf("%016x", n); break; + case 16: (void)printf("%032x", n); break; + default: (void)printf("%016x", n); break; + } +} + +/** @cond */ +#define BZ_INT_CHUNKS ((unsigned int)8) /* Assume BigNumDigit < (8 * int) */ +/** @endcond */ + +/** + * @brief Prints debugging information for a BigNum. + * This function prints the following information: + * - A message `m` (if provided). + * - The BigNum string representation `bzstr` and its address. + * - The size of a BigNum digit and a pointer in bytes. + * - The high `nl` digits of the BigNum `n`. + * - The sign of the BigNum `sign`. + * @param [in] m the message to print (optional, can be NULL). + * @param [in] bzstr the BigNum string representation. + * @param [in] n the BigNum. + * @param [in] nl the number of high digits to print. + * @param [in] sign the sign of the BigNum. + */ +void +BnDebug(const char *m, + const BzChar *bzstr, + const BigNum n, + BigNumLength nl, + BzSign sign) { + BigNumLength l; + BzChar c; + + if (m != NULL) { + (void)printf("%-20s\n", m); + } + + (void)printf("\t: BigZ = %s at 0x%p, ", (char *)bzstr, n); + (void)printf("digit = %d, word = %d bytes\n", + (int)BN_DIGIT_SIZE, + (int)sizeof(void *)); + (void)printf("\t: <-- high %02d digit(s) low -->\n", (int)nl); + + if (sign == BZ_ZERO) { + c = (BzChar)'0'; + } else if (sign == BZ_MINUS) { + c = (BzChar)'-'; + } else { + c = (BzChar)'+'; + } + + (void)printf(" %c : ", c); + + for (l = nl; l-- != 0;) { + BigNumDigit d = n[l]; + unsigned int dsize = (unsigned int)sizeof(d); + unsigned int isize = (unsigned int)sizeof(dsize); /* int */ + + (void)printf("|"); + + if ((dsize > isize) && ((dsize / isize) <= BZ_INT_CHUNKS)) { + /* + * sizeof(BigNumDigit) > sizeof(int). + * (can be 64bit BigNumDigit and 32bit int). + * NOTE: this code also works if isize == dsize. + */ + unsigned int chunk[BZ_INT_CHUNKS]; + unsigned int mask; + unsigned int shift; + unsigned int i; + + mask = ~((unsigned int)0); /* 0xff..ff */ + shift = (unsigned int)(isize * BN_BYTE_SIZE); + + /* + * Split BigNumDigit in int chunks. + */ + + for (i = 0; i < (dsize / isize); ++i) { + chunk[i] = (unsigned int)(d & mask); + d >>= shift; + } + + while (i-- != 0) { + BzShowUnsingnedInt(chunk[i]); + } + } else { + /* + * sizeof(BigNumDigit) <= sizeof(int). + */ + BzShowUnsingnedInt((unsigned int)d); + } + } + + (void)printf("|\n"); + (void)printf(" %c : ", c); + + for (l = nl; l-- != 0;) { + (void)printf("|"); + BzShowBits(n[l]); + } + + (void)printf("|\n"); +} + +/** + * BzDebug + * @param [in] m C string + * @param [in] y BigZ + */ +void +BzDebug(const char *m, const BigZ y) { + BzChar *s = BzToString(y, (BigNumDigit)10, 0); + + BnDebug(m, s, BzToBn(y), BzNumDigits(y), BzGetSign(y)); + BzFreeString(s); +} +#endif + +/** + * BzVersion + * @return C string + */ +const char * +BzVersion(void) { + return BZ_VERSION; +} + +/* + * constants used by BzToString() and BzFromString() + */ + +/** @cond */ +#define BZ_MIN_BASE 2 +#define BZ_MAX_BASE 36 +/** @endcond */ + +/* + * following table is computed using: + * + * for (i = 1; i <= BZ_MAX_BASE; ++i) { + * printf("\t%16.16f, // log(%2d)\n", log((double)i), i); + * } + */ + +/** + * @brief An array of precomputed logarithms for integers 1 to 36. + * This array is used for optimization purposes. + */ +static const double BzLog[] = { + 0.0000000000000000, + 0.0000000000000000, /* log(1) */ + 0.6931471805599453, /* log(2) */ + 1.0986122886681098, /* log(3) */ + 1.3862943611198906, /* log(4) */ + 1.6094379124341003, /* log(5) */ + 1.7917594692280550, /* log(6) */ + 1.9459101490553132, /* log(7) */ + 2.0794415416798357, /* log(8) */ + 2.1972245773362196, /* log(9) */ + 2.3025850929940459, /* log(10) */ + 2.3978952727983707, /* log(11) */ + 2.4849066497880004, /* log(12) */ + 2.5649493574615367, /* log(13) */ + 2.6390573296152584, /* log(14) */ + 2.7080502011022101, /* log(15) */ + 2.7725887222397811, /* log(16) */ + 2.8332133440562162, /* log(17) */ + 2.8903717578961645, /* log(18) */ + 2.9444389791664403, /* log(19) */ + 2.9957322735539909, /* log(20) */ + 3.0445224377234230, /* log(21) */ + 3.0910424533583161, /* log(22) */ + 3.1354942159291497, /* log(23) */ + 3.1780538303479458, /* log(24) */ + 3.2188758248682006, /* log(25) */ + 3.2580965380214821, /* log(26) */ + 3.2958368660043291, /* log(27) */ + 3.3322045101752038, /* log(28) */ + 3.3672958299864741, /* log(29) */ + 3.4011973816621555, /* log(30) */ + 3.4339872044851463, /* log(31) */ + 3.4657359027997265, /* log(32) */ + 3.4965075614664802, /* log(33) */ + 3.5263605246161616, /* log(34) */ + 3.5553480614894135, /* log(35) */ + 3.5835189384561099, /* log(36) */ +}; + +/** + * BzCreate + * Allocates a zeroed BigZ of the desired size. + * @param [in] size BigNumLength + * @return BigZ + */ +BigZ +BzCreate(BigNumLength size) { + BigZ z; + size_t chunk; + + /* + * Compute BigZ allocation size taking care of aligment. + */ + chunk = sizeof(BigZStruct) + - (BZ_DUMMY_SIZE * sizeof(BigNumDigit)) + + (size * sizeof(BigNumDigit)); + + if ((z = (BigZ)(BzAlloc(chunk))) != BZNULL) { + /* + * reset digits + */ + + BnnSetToZero(BzToBn(z), size); + + /* + * init header + */ + + BzSetSize(z, size); + BzSetSign(z, BZ_ZERO); + } + + return z; +} + +/** + * BzNumDigits + * Returns the number of digits used by z. + * @param [in] z BigZ + * @return BigNumLength + */ +BigNumLength +BzNumDigits(const BigZ z) { + return BnnNumDigits(BzToBn(z), BzGetSize(z)); +} + +/** + * BzLength + * Returns the number of bits used by z. + * @param [in] z BigZ + * @return BigNumLength + */ +BigNumLength +BzLength(const BigZ z) { + BigNumLength nl; + + switch (BzGetSign(z)) { + case BZ_MINUS : + nl = BnnNumLength(BzToBn(z), BzNumDigits(z)); + if (BnnIsPower2(BzToBn(z), BzNumDigits(z)) == BN_TRUE) { + return nl - 1; + } else { + return nl; + } + case BZ_PLUS : + return BnnNumLength(BzToBn(z), BzNumDigits(z)); + case BZ_ZERO : + default : + return (BigNumLength)0; + } +} + +/** + * BzCopy + * Creates a copy of the passed BigZ. + * @param [in] z BigZ + * @return BigZ + */ +BigZ +BzCopy(const BigZ z) { + BigZ y; + BigNumLength zl; + + zl = BzNumDigits(z); + + if ((y = BzCreate(zl)) != BZNULL) { + /* + * copy the digits + */ + BnnAssign(BzToBn(y), BzToBn(z), zl); + + /* + * copy the header WITHOUT the size !! + */ + BzSetSign(y, BzGetSign(z)); + } + + return y; +} + +/** + * BzGetOppositeSign + * @param [in] z BigZ + * @return BzSign + */ +static BzSign +BzGetOppositeSign(const BigZ z) { + switch (BzGetSign(z)) { + case BZ_MINUS: + return BZ_PLUS; + case BZ_ZERO: + return BZ_ZERO; + case BZ_PLUS: + default: + return BZ_MINUS; + } +} + +/** + * BzNegate + * Negates the passed BigZ. + * @param [in] z BigZ + * @return BigZ + * @pre z != BZNULL. + */ +BigZ +BzNegate(const BigZ z) { + BigZ y; + + if ((y = BzCopy(z)) != BZNULL) { + switch (BzGetSign(z)) { + case BZ_MINUS: + BzSetSign(y, BZ_PLUS); + break; + case BZ_ZERO: + BzSetSign(y, BZ_ZERO); + break; + case BZ_PLUS: + BzSetSign(y, BZ_MINUS); + break; + } + } + + return y; +} + +/** + * BzAbs + * Takes the absolute value of the passed BigZ. + * @param [in] z BigZ + * @return BigZ + */ +BigZ +BzAbs(const BigZ z) { + BigZ y; + + if ((y = BzCopy(z)) != BZNULL) { + if (BzGetSign(z) == BZ_MINUS) { + BzSetSign(y, BZ_PLUS); + } + } + + return y; +} + +/** + * BzCompare. + * Returns: + * - BZ_GT if y > z, + * - BZ_LT if y < z, + * - BZ_EQ otherwise. + * @param [in] y BigZ + * @param [in] z BigZ + * @return BzCmp + * @pre y != BZNULL. + * @pre z != BZNULL. + */ +BzCmp +BzCompare(const BigZ y, const BigZ z) { + if (BzGetSign(y) > BzGetSign(z)) { + return BZ_GT; + } else if (BzGetSign(y) < BzGetSign(z)) { + return BZ_LT; + } else if (BzGetSign(y) == BZ_PLUS) { + return (BzCmp)BnnCompare(BzToBn(y), BzGetSize(y), + BzToBn(z), BzGetSize(z)); + } else if (BzGetSign(y) == BZ_MINUS) { + return (BzCmp)BnnCompare(BzToBn(z), BzGetSize(z), + BzToBn(y), BzGetSize(y)); + } else { + return BZ_EQ; + } +} + +/** + * BzAdd + * Returns y + z. + * @param [in] y BigZ + * @param [in] z BigZ + * @return BigZ + * @pre y != BZNULL. + * @pre z != BZNULL. + */ +BigZ +BzAdd(const BigZ y, const BigZ z) { + BigZ n; + BigNumLength yl; + BigNumLength zl; + + BzHash(z); + + + + yl = BzNumDigits(y); + zl = BzNumDigits(z); + + if (BzGetSign(y) == BzGetSign(z)) { + /* + * Add magnitudes if signs are the same + */ + switch (BnnCompare(BzToBn(y), yl, BzToBn(z), zl)) { + case BN_EQ: + case BN_GT: /* |Y| >= |Z| */ + if ((n = BzCreate(yl + 1)) != BZNULL) { + BnnAssign(BzToBn(n), BzToBn(y), yl); + (void)BnnAdd(BzToBn(n), + yl + 1, + BzToBn(z), + zl, + BN_NOCARRY); + BzSetSign(n, BzGetSign(y)); + } + break; + case BN_LT: + default: /* |Y| < |Z| */ + if ((n = BzCreate(zl+1)) != BZNULL) { + BnnAssign(BzToBn(n), BzToBn(z), zl); + (void)BnnAdd(BzToBn(n), + zl + 1, + BzToBn(y), + yl, + BN_NOCARRY); + BzSetSign(n, BzGetSign(z)); + } + break; + } + } else { + /* + * Subtract magnitudes if signs are different + */ + switch (BnnCompare(BzToBn(y), yl, BzToBn(z), zl)) { + case BN_EQ: /* Y = -Z */ + n = BzCreate((BigNumLength)1); + break; + case BN_GT: /* |Y| > |Z| */ + if ((n = BzCreate(yl)) != BZNULL) { + BnnAssign(BzToBn(n), BzToBn(y), yl); + (void)BnnSubtract(BzToBn(n), + yl, + BzToBn(z), + zl, + BN_CARRY); + BzSetSign(n, BzGetSign(y)); + } + break; + case BN_LT: + default: /* |Y| < |Z| */ + if ((n = BzCreate(zl)) != BZNULL) { + BnnAssign(BzToBn(n), BzToBn(z), zl); + (void)BnnSubtract(BzToBn(n), + zl, + BzToBn(y), + yl, + BN_CARRY); + BzSetSign(n, BzGetSign(z)); + } + break; + } + } + + return n; +} + +/** + * BzSubtract + * Returns y - z. + * @param [in] y BigZ + * @param [in] z BigZ + * @return BigZ + * @pre y != BZNULL. + * @pre z != BZNULL. + */ +BigZ +BzSubtract(const BigZ y, const BigZ z) { + BigZ diff; + + if (y != z) { + BzSetSign(z, BzGetOppositeSign(z)); + diff = BzAdd(y, z); + BzSetSign(z, BzGetOppositeSign(z)); + } else { + diff = BzFromInteger((BzInt)0); + } + + return diff; +} + +/** + * BzMultiply + * Returns y * z. + * @param [in] y BigZ + * @param [in] z BigZ + * @return BigZ + * @pre y != BZNULL. + * @pre z != BZNULL. + */ +BigZ +BzMultiply(const BigZ y, const BigZ z) { + BigZ n; + BigNumLength yl; + BigNumLength zl; + + if (BzGetSign(y) == BZ_ZERO || BzGetSign(z) == BZ_ZERO) { + return BzFromInteger((BzInt)0); + } + + yl = BzNumDigits(y); + zl = BzNumDigits(z); + + if ((n = BzCreate(yl + zl)) == BZNULL) { + return BZNULL; + } + + (void)BnnMultiply(BzToBn(n), yl + zl, BzToBn(y), yl, BzToBn(z), zl); + + if (BzGetSign(y) == BzGetSign(z)) { + BzSetSign(n, BZ_PLUS); + } else { + BzSetSign(n, BZ_MINUS); + } + + return n; +} + +/** + * BzDivide. + * - y div z => q + * - y mod z => r, + * + * such that y = zq + r + * @pre z != NULL, *z == BZNULL + * @post *z != BZNULL + * @param [in] y BigZ + * @param [in] z BigZ + * @param [out] r BigZ + * @return BigZ + * @pre y != BZNULL. + * @pre z != BZNULL. + * @pre r != NULL. + */ +BigZ +BzDivide(const BigZ y, const BigZ z, BigZ *r) { + BigZ q; + BigNumLength yl; + BigNumLength zl; + BigNumLength ql; + BigNumLength rl; + + if (BzGetSign(z) == BZ_ZERO) { + return BZNULL; + } + + yl = BzNumDigits(y); + zl = BzNumDigits(z); + ql = (BigNumLength)BzMaxInt((int)yl - (int)zl + 1, 1) + 1; + rl = (BigNumLength)BzMaxInt(zl, yl) + 1; + + /* + * Set up quotient, remainder + */ + + if ((q = BzCreate(ql)) == BZNULL) { + return BZNULL; + } + + if ((*r = BzCreate(rl)) == BZNULL) { + BzFree(q); + return BZNULL; + } + + BnnAssign(BzToBn(*r), BzToBn(y), yl); + + /* + * Do the division + */ + + BnnDivide(BzToBn(*r), rl, BzToBn(z), zl); + BnnAssign(BzToBn(q), BzToBn(*r) + zl, rl - zl); + BnnSetToZero(BzToBn(*r) + zl, rl - zl); + rl = zl; + + /* + * Correct the signs, adjusting the quotient and remainder + */ + + if (BnnIsZero(BzToBn(*r), rl) == BN_FALSE) { + /* + * R<>0 + */ + if (BzGetSign(y) != BzGetSign(z)) { + /* + * (Z-R) => R + */ + BnnComplement(BzToBn(*r), rl); + (void)BnnAdd(BzToBn(*r), rl, BzToBn(z), zl, BN_CARRY); + /* + * (Q+1) => Q + */ + (void)BnnAddCarry(BzToBn(q), ql, BN_CARRY); + } + /* + * The sign of the result (mod) is the sign of z. + */ + BzSetSign(*r, BzGetSign(z)); + } else { + /* + * Set sign to BZ_ZERO + * (already made by BzCreate but makes it clear) + */ + BzSetSign(*r, BZ_ZERO); + } + + /* + * Set the sign of the quotient. + */ + + if (BnnIsZero(BzToBn(q), ql) == BN_TRUE) { + BzSetSign(q, BZ_ZERO); + } else if (BzGetSign(y) == BzGetSign(z)) { + BzSetSign(q, BZ_PLUS); + } else { + BzSetSign(q, BZ_MINUS); + } + + return q; +} + +/** + * BzDiv. + * Returns div(y, z). + * @param [in] y BigZ + * @param [in] z BigZ + * @return BigZ + */ +BigZ +BzDiv(const BigZ y, const BigZ z) { + BigZ q; + BigZ r = BZNULL; + + q = BzDivide(y, z, &r); + + if (r != BZNULL) { + BzFree(r); + } + + return q; +} + +/** + * BzTruncate. + * Returns truncate(y, z). + * @param [in] y BigZ + * @param [in] z BigZ + * @return BigZ + * @pre y != BZNULL. + * @pre z != BZNULL. + */ +BigZ +BzTruncate(const BigZ y, const BigZ z) { + BigZ q; + BigZ r = BZNULL; + BigNumLength ql; + + q = BzDivide(y, z, &r); + + if (r == BZNULL) { + /* + * This should never happend. + */ + BzFree(q); + return BZNULL; + } + + ql = BzNumDigits(q); + + if ((BzGetSign(q) == BZ_MINUS) && (BzGetSign(r) != BZ_ZERO)) { + /* + * Q < 0, R <> 0, 2*R>= Z : Q-1 => Q + */ + (void)BnnSubtractBorrow(BzToBn(q), ql, BN_NOCARRY); + + if (BnnIsZero(BzToBn(q), ql) == BN_TRUE) { + BzSetSign(q, BZ_ZERO); + } + } else if ((BnnIsZero(BzToBn(q), ql) == BN_TRUE) + && (BzGetSign(y) == BzGetSign(z))) { + /* + * Q == 0, sign(Y) == sign(Z) : 0 => Q + */ + BzFree(q); + q = BzFromInteger((BzInt)0); + } + + BzFree(r); + + return q; +} + +/** + * BzFloor. + * Returns floor(y, z). + * @param [in] y BigZ + * @param [in] z BigZ + * @pre y != BZNULL. + * @pre z != BZNULL. + * @return BigZ + */ +BigZ +BzFloor(const BigZ y, const BigZ z) { + BigZ q; + BigZ r = BZNULL; + + q = BzDivide(y, z, &r); + + if (r != BZNULL) { + BzFree(r); + } + + return q; +} + +/** + * BzCeiling. + * Returns ceiling(y, z). + * @param [in] y BigZ + * @param [in] z BigZ + * @pre y != BZNULL. + * @pre z != BZNULL. + * @return BigZ + */ +BigZ +BzCeiling(const BigZ y, const BigZ z) { + BigZ q; + BigZ r = BZNULL; + BigNumLength ql; + + if ((q = BzDivide(y, z, &r)) == BZNULL) { + return BZNULL; + } + + if (r == BZNULL) { + /* + * This should never happend. + */ + BzFree(q); + return BZNULL; + } + + ql = BzNumDigits(q); + + if (BzGetSign(q) == BZ_PLUS && BzGetSign(r) != BZ_ZERO) { + /* + * Q > 0, R <> 0 : Q+1 => Q + */ + BigNumDigit one = BN_ONE; + (void)BnnAdd(BzToBn(q), ql, &one, (BigNumLength)1, BN_NOCARRY); + } else if (BzGetSign(q) == BZ_MINUS && BzGetSign(r) != BZ_ZERO) { + /* + * Q < 0, R <> 0 : Q-1 => Q + */ + (void)BnnSubtractBorrow(BzToBn(q), ql, BN_NOCARRY); + if (BnnIsZero(BzToBn(q), ql) == BN_TRUE) { + BzSetSign(q, BZ_ZERO); + } + } else if (BzGetSign(q) == BZ_ZERO + && BzGetSign(y) == BzGetSign(z)) { + /* + * Q == 0, sign(Y) == sign(Z) : 1 => Q + */ + BzFree(q); + q = BzFromInteger((BzInt)1); + } + + BzFree(r); + + return q; +} + +/** + * BzRound. + * Returns round(y, z). + * @param [in] y BigZ + * @param [in] z BigZ + * @return BigZ + * @pre y != BZNULL. + * @pre z != BZNULL. + */ +BigZ +BzRound(const BigZ y, const BigZ z) { + BigZ q; + BigZ r = BZNULL; + BigNumLength ql; + + if ((q = BzDivide(y, z, &r)) == BZNULL) { + return BZNULL; + } + + if (r == BZNULL) { + /* + * This should never happend. + */ + BzFree(q); + return BZNULL; + } + + ql = BzNumDigits(q); + + if (BzGetSign(q) == BZ_PLUS && BzGetSign(r) != BZ_ZERO) { + BigNumDigit one = BN_ONE; + BigZ roundz; + BzSign sign = BzGetSign(z); + + BzSetSign(r, BZ_PLUS); + BzSetSign(z, BZ_PLUS); + + roundz = BzAsh(r, 1); + + switch (BzCompare(roundz, z)) { + case BZ_LT : + break; + case BZ_EQ : + /* + * Q > 0, R <> 0, 2*R= Z : + * + * roundz is exactly halfway between two integers, + * choose even number. + */ + if (BzIsOdd(q) == BN_TRUE) { + (void)BnnAdd(BzToBn(q), + ql, + &one, + (BigNumLength)1, + BN_NOCARRY); + } + break; + case BZ_GT : + /* + * Q > 0, R <> 0, 2*R>= Z : Q+1 => Q + */ + (void)BnnAdd(BzToBn(q), + ql, + &one, + (BigNumLength)1, + BN_NOCARRY); + break; + } + + BzSetSign(z, sign); + BzFree(roundz); + } else if (BzGetSign(q) == BZ_MINUS && BzGetSign(r) != BZ_ZERO) { + /* + * Q < 0, R <> 0, 2*R>= Z : Q-1 => Q + */ + BigZ roundz; + BzSign sign = BzGetSign(z); + + BzSetSign(r, BZ_PLUS); + BzSetSign(z, BZ_PLUS); + + roundz = BzAsh(r, 1); + + switch (BzCompare(roundz, z)) { + case BZ_LT : + break; + case BZ_EQ : + /* + * roundz is exactly halfway between two integers, + * choose even number. + */ + if (BzIsOdd(q) == BN_TRUE) { + (void)BnnSubtractBorrow(BzToBn(q), + ql, + BN_NOCARRY); + } + break; + case BZ_GT : + (void)BnnSubtractBorrow(BzToBn(q), ql, BN_NOCARRY); + break; + } + + BzSetSign(z, sign); + BzFree(roundz); + + if (BnnIsZero(BzToBn(q), ql) == BN_TRUE) { + BzSetSign(q, BZ_ZERO); + } + } else if ((BzGetSign(q) == BZ_ZERO) + && (BzGetSign(y) == BzGetSign(z))) { + /* + * Q == 0, sign(Y) == sign(Z): + */ + BigZ roundz; + BzSign sign = BzGetSign(z); + + BzFree(q); + + BzSetSign(r, BZ_PLUS); + BzSetSign(z, BZ_PLUS); + + roundz = BzAsh(r, 1); + + if (BzCompare(roundz, z) == BZ_LT) { + /* + * 2*R< Z : 0 => Q + */ + q = BzFromInteger((BzInt)0); + } else { + /* + * 2*R>= Z : 1 => Q + */ + q = BzFromInteger((BzInt)1); + } + + BzSetSign(z, sign); + BzFree(roundz); + } + + BzFree(r); + + return q; +} + +/** + * BzMod. + * Returns y mod z. + * @param [in] y BigZ + * @param [in] z BigZ + * @return BigZ + * @pre y != BZNULL. + * @pre z != BZNULL. + */ +BigZ +BzMod(const BigZ y, const BigZ z) { + BigZ r = BZNULL; + + BzFree(BzDivide(y, z, &r)); + + return r; +} + +/** + * BzRem. + * Returns y rem z. + * @param [in] y BigZ + * @param [in] z BigZ + * @return BigZ + * @pre y != BZNULL. + * @pre z != BZNULL. + */ +BigZ +BzRem(const BigZ y, const BigZ z) { + BigZ q; + BigZ r = BZNULL; + BigZ rem; + + if ((q = BzDivide(y, z, &r)) != BZNULL) { + BzFree(q); + + if (r == BZNULL) { + /* + * This should never happend. + */ + return BZNULL; + } + + if (BzGetSign(r) == BZ_ZERO) { + return r; + } else if (BzGetSign(y) == BzGetSign(z)) { + return r; + } else if (BzGetSign(y) == BZ_MINUS) { + if ((rem = BzSubtract(z, r)) != BZNULL) { + BzSetSign(rem, BZ_MINUS); + } + BzFree(r); + return rem; + } else { + if ((rem = BzSubtract(r, z)) != BZNULL) { + BzSetSign(rem, BZ_PLUS); + } + BzFree(r); + return rem; + } + } else { + return BZNULL; + } +} + +/** + * BzIsEven. + * Returns BN_TRUE iff y is even + * @param [in] y BigZ + * @return BigNumBool + * @pre y != BZNULL. + */ +BigNumBool +BzIsEven(const BigZ y) { + return BnnIsDigitEven(BzGetDigit(y, 0)); +} + +/** + * BzIsOdd. + * Returns BN_TRUE iff y is odd + * @param [in] y BigZ + * @return BigNumBool + * @pre y != BZNULL. + */ +BigNumBool +BzIsOdd(const BigZ y) { + return BnnIsDigitOdd(BzGetDigit(y, 0)); +} + +/** @cond */ +#if defined(BZ_OPTIMIZE_PRINT) +typedef struct { + int MaxDigits; + BigNumDigit MaxValue; +} BzPrintTable; + +#if (BZ_BUCKET_SIZE == 32) +static const BzPrintTable BzPrintBase[] = { + { 0, (BigNumDigit)0U }, /* 0 */ + { 0, (BigNumDigit)0U }, /* 1 */ + { 31, (BigNumDigit)2147483648U }, /* 2 */ + { 20, (BigNumDigit)3486784401U }, /* 3 */ + { 15, (BigNumDigit)1073741824U }, /* 4 */ + { 13, (BigNumDigit)1220703125U }, /* 5 */ + { 12, (BigNumDigit)2176782336U }, /* 6 */ + { 11, (BigNumDigit)1977326743U }, /* 7 */ + { 10, (BigNumDigit)1073741824U }, /* 8 */ + { 10, (BigNumDigit)3486784401U }, /* 9 */ + { 9, (BigNumDigit)1000000000U }, /* 10 */ + { 9, (BigNumDigit)2357947691U }, /* 11 */ + { 8, (BigNumDigit)429981696U }, /* 12 */ + { 8, (BigNumDigit)815730721U }, /* 13 */ + { 8, (BigNumDigit)1475789056U }, /* 14 */ + { 8, (BigNumDigit)2562890625U }, /* 15 */ + { 7, (BigNumDigit)268435456U }, /* 16 */ + { 7, (BigNumDigit)410338673U }, /* 17 */ + { 7, (BigNumDigit)612220032U }, /* 18 */ + { 7, (BigNumDigit)893871739U }, /* 19 */ + { 7, (BigNumDigit)1280000000U }, /* 20 */ + { 7, (BigNumDigit)1801088541U }, /* 21 */ + { 7, (BigNumDigit)2494357888U }, /* 22 */ + { 7, (BigNumDigit)3404825447U }, /* 23 */ + { 6, (BigNumDigit)191102976U }, /* 24 */ + { 6, (BigNumDigit)244140625U }, /* 25 */ + { 6, (BigNumDigit)308915776U }, /* 26 */ + { 6, (BigNumDigit)387420489U }, /* 27 */ + { 6, (BigNumDigit)481890304U }, /* 28 */ + { 6, (BigNumDigit)594823321U }, /* 29 */ + { 6, (BigNumDigit)729000000U }, /* 30 */ + { 6, (BigNumDigit)887503681U }, /* 31 */ + { 6, (BigNumDigit)1073741824U }, /* 32 */ + { 6, (BigNumDigit)1291467969U }, /* 33 */ + { 6, (BigNumDigit)1544804416U }, /* 34 */ + { 6, (BigNumDigit)1838265625U }, /* 35 */ + { 6, (BigNumDigit)2176782336U } /* 36 */ +}; +#endif /* BZ_BUCKET_SIZE == 32 */ + +#if (BZ_BUCKET_SIZE == 64) +static const BzPrintTable BzPrintBase[] = { + { 0, (BigNumDigit)0UL }, /* 0 */ + { 0, (BigNumDigit)0UL }, /* 1 */ + { 63, (BigNumDigit)9223372036854775808UL }, /* 2 */ + { 40, (BigNumDigit)12157665459056928801UL }, /* 3 */ + { 31, (BigNumDigit)4611686018427387904UL }, /* 4 */ + { 27, (BigNumDigit)7450580596923828125UL }, /* 5 */ + { 24, (BigNumDigit)4738381338321616896UL }, /* 6 */ + { 22, (BigNumDigit)3909821048582988049UL }, /* 7 */ + { 21, (BigNumDigit)9223372036854775808UL }, /* 8 */ + { 20, (BigNumDigit)12157665459056928801UL }, /* 9 */ + { 19, (BigNumDigit)10000000000000000000UL }, /* 10 */ + { 18, (BigNumDigit)5559917313492231481UL }, /* 11 */ + { 17, (BigNumDigit)2218611106740436992UL }, /* 12 */ + { 17, (BigNumDigit)8650415919381337933UL }, /* 13 */ + { 16, (BigNumDigit)2177953337809371136UL }, /* 14 */ + { 16, (BigNumDigit)6568408355712890625UL }, /* 15 */ + { 15, (BigNumDigit)1152921504606846976UL }, /* 16 */ + { 15, (BigNumDigit)2862423051509815793UL }, /* 17 */ + { 15, (BigNumDigit)6746640616477458432UL }, /* 18 */ + { 15, (BigNumDigit)15181127029874798299UL }, /* 19 */ + { 14, (BigNumDigit)1638400000000000000UL }, /* 20 */ + { 14, (BigNumDigit)3243919932521508681UL }, /* 21 */ + { 14, (BigNumDigit)6221821273427820544UL }, /* 22 */ + { 14, (BigNumDigit)11592836324538749809UL }, /* 23 */ + { 13, (BigNumDigit)876488338465357824UL }, /* 24 */ + { 13, (BigNumDigit)1490116119384765625UL }, /* 25 */ + { 13, (BigNumDigit)2481152873203736576UL }, /* 26 */ + { 13, (BigNumDigit)4052555153018976267UL }, /* 27 */ + { 13, (BigNumDigit)6502111422497947648UL }, /* 28 */ + { 13, (BigNumDigit)10260628712958602189UL }, /* 29 */ + { 13, (BigNumDigit)15943230000000000000UL }, /* 30 */ + { 12, (BigNumDigit)787662783788549761UL }, /* 31 */ + { 12, (BigNumDigit)1152921504606846976UL }, /* 32 */ + { 12, (BigNumDigit)1667889514952984961UL }, /* 33 */ + { 12, (BigNumDigit)2386420683693101056UL }, /* 34 */ + { 12, (BigNumDigit)3379220508056640625UL }, /* 35 */ + { 12, (BigNumDigit)4738381338321616896UL } /* 36 */ +}; +#endif /* BZ_BUCKET_SIZE == 64 */ +#endif /* BZ_OPTIMIZE_PRINT */ +/** @endcond */ + +/** + * BzToString + * wrapper to BzToStringBuffer that always allocate buffer. + * The returned string must be deallocated using BzFreeString. + * @param [in] z BigZ + * @param [in] base BigNumDigit + * @param [in] sign flag. When set to BZ_FORCE_SIGN and z != 0, + * a sign +/- is always present. + * @return NUL terminated string + * @pre z != BZNULL. + */ +BzChar * +BzToString(const BigZ z, BigNumDigit base, int sign) { + return BzToStringBuffer(z, base, sign, (BzChar *)0, (size_t *)0); +} + +/** + * BzToStringBuffer. + * If buf is passed as NULL, a new string is allocated using BzStringAlloc and + * must be deallocated by BzFreeString. + * @param [in] z BigZ + * @param [in] base BigNumDigit + * @param [in] sign flag. When set to BZ_FORCE_SIGN and z != 0, + * a sign +/- is always present. + * @param [out] buf output buffer + * @param [in] len buffer length + * @return NUL terminated string + * @pre z != BZNULL. + */ +BzChar * +BzToStringBuffer(const BigZ z, + BigNumDigit base, + int sign, + BzChar * const buf, + size_t *len) { + return BzToStringBufferExt(z, + base, + sign, + buf, + len, + (size_t *)0); +} + +/** + * BzToStringBufferExt. + * Returns a pointer to a string that represents Z in the specified + * base. Assumes BZ_MIN_BASE <= base <= BZ_MAX_BASE. If optional + * buffer is supplied, len is a pointer of this buffer size. If there + * is enough room to print the number buf is used otherwise function + * returns NULL and len contains the required size. If buf is passed + * as NULL, this string is allocated on the heap, so it must be + * desallocated by the user. + * @param [in] z BigZ + * @param [in] base BigNumDigit + * @param [in] sign flag. When set to BZ_FORCE_SIGN and z != 0, + * a sign +/- is always present. + * @param [out] buf output buffer + * @param [in, out] len on input it contains the maximal buffer length. When + * buffer is too small to represent this BigZ number, this function returns + * NULL and len is set the the required buffer size. + * @param [out] slen when not NULL, slen is a pointer of size_t which contains + * the actual string length of z. + * @return NUL terminated string + * @pre z != BZNULL. + */ +BzChar * +BzToStringBufferExt(const BigZ z, + BigNumDigit base, + int sign, + BzChar * const buf, + size_t *len, + size_t *slen) { + static const BzChar Digit[] = { + (BzChar)'0', (BzChar)'1', (BzChar)'2', (BzChar)'3', + (BzChar)'4', (BzChar)'5', (BzChar)'6', (BzChar)'7', + (BzChar)'8', (BzChar)'9', (BzChar)'a', (BzChar)'b', + (BzChar)'c', (BzChar)'d', (BzChar)'e', (BzChar)'f', + (BzChar)'g', (BzChar)'h', (BzChar)'i', (BzChar)'j', + (BzChar)'k', (BzChar)'l', (BzChar)'m', (BzChar)'n', + (BzChar)'o', (BzChar)'p', (BzChar)'q', (BzChar)'r', + (BzChar)'s', (BzChar)'t', (BzChar)'u', (BzChar)'v', + (BzChar)'w', (BzChar)'x', (BzChar)'y', (BzChar)'z' + }; + + BigZ y; + BigZ q; + BigNumLength zl; + BigNumLength sl; + BzChar * s; + BzChar * slast; + BzChar * strg; + + if (base < (BigNumDigit)BZ_MIN_BASE + || base > (BigNumDigit)BZ_MAX_BASE) { + if (len != 0) { + *len = 0; + } + return (BzChar *)NULL; + } + + /* + * Allocate BigNums and set up string + */ + + zl = BzNumDigits(z) + 1; + sl = (BigNumLength)((BzLog[2] * BN_DIGIT_SIZE * zl) / BzLog[base] + 3); + + if (buf != (BzChar *)NULL + && len != (size_t *)NULL + && (sl > (BigNumLength)*len)) { + /* + * a buffer is passed but there is not enough room, + * return NULL and set required size in len. + */ + *len = (size_t)sl; + return (BzChar *)NULL; + } + + if (len != (size_t *)NULL) { + /* + * set len to 0 in case BzToStringBuffer returns NULL + * because of allocation failure. It should be checked by + * caller that may otherwise allocate a bigger buffer. + */ + *len = 0; + } + + if ((y = BzCreate(zl)) == BZNULL) { + return (BzChar *)NULL; + } + + if ((q = BzCreate(zl)) == BZNULL) { + BzFree(y); + return (BzChar *)NULL; + } + + if (buf != (BzChar *)NULL) { + strg = buf; + } else { + strg = (BzChar *)BzStringAlloc((size_t)sl); + if (strg == (BzChar *)NULL) { + BzFree(y); + BzFree(q); + return (BzChar *)NULL; + } + if (len != (size_t *)NULL) { + /* + * a buffer is allocated and caller wants to know + * allocated size. + */ + *len = (size_t)sl; + } + } + + BnnAssign(BzToBn(y), BzToBn(z), zl - 1); + s = strg + sl; + slast = s; /* remember last position in order to compute size */ + + /* + * Divide Z by base repeatedly; successive digits given by remainders + */ + + *--s = (BzChar)'\0'; + + if (BzGetSign(z) == BZ_ZERO) { + *--s = (BzChar)'0'; +#if defined(BZ_OPTIMIZE_PRINT) + } else { + /* + * Compute maxval and digits that can be used with + * this base. + */ + + BigNumDigit maxval = (BigNumDigit)BzPrintBase[base].MaxValue; + BigNumLength digits = (BigNumLength)BzPrintBase[base].MaxDigits; + + /* + * This optimization makes BigZ output 10 to 20x faster. + */ + do { + BigZ v; + BigNumDigit r; + /* + * compute: y div maxval => q, + * returns r = y mod maxval + * + * maxval is the greatest integer in base 'base' + * that fits in a BigNumDigit. + */ + + r = BnnDivideDigit(BzToBn(q), + BzToBn(y), + zl, + maxval); + + if (BnnIsZero(BzToBn(q), zl) == BN_FALSE) { + /* + * More digits to come on left, add exactly + * the number of digits with possible + * leading 0 (when r becomes 0). + */ + int i; + for (i = 0; i < (int)digits; ++i) { + if (r == 0) { + /* + * No need to divide, fill + * the rest with '0'. + */ + *--s = (BzChar)'0'; + } else { + *--s = Digit[r % base]; + r = r / base; + } + } + } else { + /* + * Last serie (top left). Print only available + * digits (stop when r becomes 0). + */ + while (r != 0) { + *--s = Digit[r % base]; + r = r / base; + } + } + + /* + * exchange y and q (to avoid BzMove(y, q)) + */ + + v = q; + q = y; + y = v; + } while (BnnIsZero(BzToBn(y), zl) == BN_FALSE); + } +#else /* BZ_OPTIMIZE_PRINT */ + } else { + do { + BigZ v; + BigNumDigit r; + /* compute: y div base => q, returns r = y mod base */ + + r = BnnDivideDigit(BzToBn(q), BzToBn(y), zl, base); + *--s = Digit[r]; + + /* + * exchange y and q (to avoid BzMove(y, q)) + */ + + v = q; + q = y; + y = v; + } while (BnnIsZero(BzToBn(y), zl) == BN_FALSE); + } +#endif /* BZ_OPTIMIZE_PRINT */ + + /* + * Add sign if needed. + */ + + switch (BzGetSign(z)) { + case BZ_MINUS: + /* + * z < 0, always add '-' sign. + */ + *--s = (BzChar)'-'; + break; + case BZ_ZERO: + /* + * z = 0 -> no sign even if BZ_FORCE_SIGN is passed. + */ + break; + case BZ_PLUS: + /* + * z > 0, add '+' only if sign is required. + */ + if (sign == BZ_FORCE_SIGN) { + *--s = (BzChar)'+'; + } + } + + if (buf != (BzChar *)NULL) { + strg = &s[0]; + } else if ((s - strg) > 0) { + /* + * and move string into position as no buffer + * has been supplied. + */ + BigNumLength i; + + for (i = 0; s[i] != (BzChar)'\000'; ++i) { + strg[i] = s[i]; + } + +#if defined(CPP_MODULE) && defined(_MSC_VER) +/* + * Microsoft VC++ incorrectly finds a Buffer overrun. + */ +#pragma warning(push) +#pragma warning(disable:6386) + strg[i] = (BzChar)'\000'; +#pragma warning(pop) +#else + strg[i] = (BzChar)'\000'; +#endif + } + + + /* + * Free temporary BigNums and return the string + */ + + BzFree(y); + BzFree(q); + + if (slen != (size_t *)NULL) { + /* + * A non null pointer was passed to get string length + * (which is not the same a buffer length used to build string). + */ + *slen = (size_t)(slast - s - 1); + } + + return strg; +} + +/** + * BzStrLen. + * @param [in] s const BzChar + * @return size_t + * @pre s != NULL. + */ +size_t +BzStrLen(const BzChar *s) { + size_t len; + + for (len = (size_t)0; *s++ != (BzChar)'\000'; ++len) { + continue; + } + + return len; +} + +/** + * BzFromStringLen. + * Creates a BigZ whose value is represented by "string" in the + * specified base. The "string" may contain leading spaces, followed + * by an optional sign, followed by a series of digits. Assumes + * BZ_MIN_BASE<=base<=BZ_MAX_BASE. When called from C, only the + * first 2 arguments are passed. The BzStrFlag flag controls the characters + * that can be ignored. + * Depending on base valid characters are + * "0123456789abcdefghijklmnopqrstuvwxy". Alphabetic characters can be either + * upper or lower case. + * @param [in] s const BzChar + * @param [in] len string length. + * @param [in] base BigNumDigit + * @param [in] flag BzStrFlag to tell how we stop reading the string. By + * default, BZ_UNTIL_END reads until the end of string. + * The BZ_UNTIL_SLASH value is used by bigq to read rational numbers + * and stop on '/' character. + * @return BigZ + * @pre base in range [BZ_MIN_BASE .. BZ_MAX_BASE]. + */ +BigZ +BzFromStringLen(const BzChar *s, size_t len, BigNumDigit base, BzStrFlag flag) { + BigZ z; + BigZ p; + BzSign sign; + BigNumLength zl; + size_t i; + + if (s == (const BzChar*)NULL) { + return BZNULL; + } + + if ((base < 2) || (base > 36)) { + return BZNULL; + } + + /* + * Throwaway any initial space + */ + + while (BZ_ISSPACE(*s)) { + ++s; + --len; + } + + /* + * Set up sign, base, initialize result + */ + + switch (*s) { + case ((BzChar)'-'): + sign = BZ_MINUS; + ++s; + --len; + break; + case ((BzChar)'+'): + sign = BZ_PLUS; + ++s; + --len; + break; + default: + sign = BZ_PLUS; + } + + if (BZ_ISSIGN(*s) || (*s == (BzChar)0)) { + /* + * There must be at most one sign. + */ + return BZNULL; + } + + /* + * Allocate BigNums + */ + + zl = (BigNumLength)(((double)len * BzLog[base]) + / (BzLog[2] * BN_DIGIT_SIZE) + 1); + + if ((z = BzCreate(zl)) == BZNULL) { + return BZNULL; + } + + if ((p = BzCreate(zl)) == BZNULL) { + BzFree(z); + return BZNULL; + } + + /* + * Multiply in the digits of the string, one at a time + */ + + for (i = 0; i < len; ++i) { + BzChar c = s[i]; + int val = CTOI(c); + BigNumDigit next = (BigNumDigit)val; + BigZ v; + + if ((val == -1) || (next >= base)) { + if (i != 0) { + if ((flag == BZ_UNTIL_SPACE) && BZ_ISSPACE(c)) { + /* + * At least one digit has been used + * and we stop on first space. The + * rest of string must be with spaces. + */ + size_t j; + for (j = i + 1; j < len; ++j) { + c = s[j]; + if (!BZ_ISSPACE(c)) { + /* + * non-space if found. + */ + BzFree(p); + BzFree(z); + return BZNULL; + } + } + break; + } else if ((flag == BZ_UNTIL_SLASH) + && (BZ_ISSLASH(c) + || BZ_ISSPACE(c))) { + /* + * At least one digit has been used + * and we stop on first '/' which is + * + */ + break; + } else if (flag == BZ_UNTIL_INVALID) { + /* + * At least one digit has been used + * and we stop on first invalid digit. + */ + break; + } + } + /* + * Invalid syntax for base. + */ + BzFree(p); + BzFree(z); + return BZNULL; + } + + BnnSetToZero(BzToBn(p), zl); + BnnSetDigit(BzToBn(p), next); + (void)BnnMultiplyDigit(BzToBn(p), zl, BzToBn(z), zl, base); + + /* + * exchange z and p (to avoid BzMove (z, p) + */ + + v = p; + p = z; + z = v; + } + + /* + * Set sign of result + */ + + BzSetSign(z, (BnnIsZero(BzToBn(z), zl) == BN_TRUE) ? BZ_ZERO : sign); + + /* + * Free temporary BigNums + */ + + BzFree(p); + + return z; +} + +/** + * BzFromString. + * Build a new BigZ form a string. Depending on base valid characters are + * "0123456789abcdefghijklmnopqrstuvwxy". Alphabetic characters can be either + * upper or lower case. + * @param [in] s const BzChar + * @param [in] base BigNumDigit between 2 and 36. + * @param [in] flag BzStrFlag + * @return BigZ + * @pre base in range [BZ_MIN_BASE .. BZ_MAX_BASE]. + */ +BigZ +BzFromString(const BzChar *s, unsigned long base, unsigned int flag) { + if (s == (const BzChar*)NULL) { + return BZNULL; + } else { + return BzFromStringLen(s, BzStrLen(s), base, flag); + } +} + +/** + * BzFromInteger. + * @param [in] i zInt + * @return BigZ + */ +BigZ +BzFromInteger(BzInt i) { + BigZ z = BzCreate((BigNumLength)1); + + if (z != BZNULL) { + BzSetDigit(z, 0, (BigNumDigit)BzAbsInt(i)); + + if (i > 0) { + BzSetSign(z, BZ_PLUS); + } else if (i < 0) { + BzSetSign(z, BZ_MINUS); + } else { + BzSetSign(z, BZ_ZERO); + } + } + + return z; +} + +/** + * BzFromUnsignedInteger + * @param [in] i zUInt + * @return BigZ + */ +BigZ +BzFromUnsignedInteger(BzUInt i) { + BigZ z = BzCreate((BigNumLength)1); + + if (z != BZNULL) { + BzSetDigit(z, 0, (BigNumDigit)i); + BzSetSign(z, ((i > (BzInt)0) ? BZ_PLUS : BZ_ZERO)); + } + + return z; +} + +/** + * BzToInteger + * Return the value of z if it can fit in an BzInt. In case of overflow, + * it returns BZMAXUINT. + * @param [in] z BigZ + * @return BzInt + * @pre z != BZNULL. + */ +BzInt +BzToInteger(const BigZ z) { + if (BzNumDigits(z) > (BigNumLength)1) { + return BZMAXINT; + } else if (BzGetSign(z) == BZ_MINUS) { + return -(BzInt)BzGetDigit(z, 0); + } else { + return (BzInt)BzGetDigit(z, 0); + } +} + +/** + * BzToLongDouble. + * It converts BigZ z to long double approximation. Conversion is + * inaccurate and looses precision if number most significant digits + * exceed long double mantissa size. + * It returns +inf/-inf if z exceeds max double values. + * @param [in] z BigZ + * @return long double + * @pre z != BZNULL. + */ +BzLDouble +BzToLongDouble(const BigZ z) { + const size_t nl = (size_t)BzNumDigits(z); + const BigNum nn = BzToBn(z); + BzLDouble res = (BzLDouble)nn[nl - 1]; /* highest bucket */ + + if (nl > 1) { + /* + * n spans on more than one bucket. + * + * Caution: ideally, we would like to compute + * (1 << BZ_BUCKET_SIZE) as the factor to apply to combine + * buckets. As it overflows, we use two steps: + * (1 << (BZ_BUCKET_SIZE - 1)) => factor + * factor * 2 => factor. + */ + BzLDouble factor = (BzLDouble)((BzUInt)1 + << (BzUInt)(BZ_BUCKET_SIZE - 1)) + * 2.0; + size_t d; + + /* + * Improve accuracy by adding the second highest bucket. + */ + res = (res * factor) + (BzLDouble)nn[nl - 2]; + + for (d = nl - 2; d > 0; --d) { + /* + * Don't waste our time summing buckets after the two + * first ones as they are below mantissa. + * It's sufficient to ajust res exponent by the number + * of remaining buckets. + */ + res = (res * factor); + + if (res > BZ_DBL_MAX) { + /* + * reached inifinty. + */ + break; + } + } + } + + if (BzGetSign(z) == BZ_MINUS) { + return -res; + } else { + return res; + } +} + +/** + * BzToDouble. + * It converts BigZ z to double approximation. Conversion is + * inaccurate and looses precision if number most significant digits + * exceed double mantissa size. + * It returns +inf/-inf if z exceeds max double values. + * @param [in] z BigZ + * @return double + * @pre z != BZNULL. + */ +double +BzToDouble(const BigZ z) { + return (double)BzToLongDouble(z); +} + +/** + * BzToIntegerPointer + * @param [in] z BigZ + * @param [in] p BzInt + * @return int + * @pre z != BZNULL. + */ +int +BzToIntegerPointer(const BigZ z, BzInt *p) { + BzInt value; + + if (BzNumDigits(z) > (BigNumLength)1) { + *p = BZMAXINT; + return 0; + } + + value = (BzInt)BzGetDigit(z, 0); + + if (value < 0) { + *p = BZMAXINT; + return 0; + } else if (BzGetSign(z) == BZ_MINUS) { + *p = -value; + } else { + *p = value; + } + + return 1; +} + +/** + * BzToUnsignedInteger. + * Return the value of z if it can fit in an BzUInt. In case of overflow, + * it returns BZMAXUINT. + * @param [in] z BigZ + * @return BzUInt + * @see BzToUnsignedIntegerPointer which tells if operation overflows. + */ +BzUInt +BzToUnsignedInteger(const BigZ z) { + if (BzNumDigits(z) > (BigNumLength)1) { + return BZMAXUINT; + } else { + return (BzUInt)BzGetDigit(z, 0); + } +} + +/** + * BzToUnsignedIntegerPointer. + * Store in p the value of z if it can fit in an BzUInt. In case of overflow, + * p is set to BZMAXUINT and BzToUnsignedIntegerPointer returns 0. + * @param [in] z BigZ + * @param [in] p BzUInt + * @return 1 if BigZ does not overlow BzUInt max value. 0 otherwise. + * @pre z != BZNULL. + * @pre p != NULL. + */ +int +BzToUnsignedIntegerPointer(const BigZ z, BzUInt *p) { + if (BzNumDigits(z) > (BigNumLength)1) { + *p = BZMAXUINT; + return 0; + } else { + *p = (BzUInt)BzGetDigit(z, 0); + + return 1; + } +} + +/** + * BzFromBigNum + * @param [in] n BigNum + * @param [in] nl BigNumLength + * @return BigZ + */ +BigZ +BzFromBigNum(const BigNum n, BigNumLength nl) { + BigZ z; + BigNumLength i; + + z = BzCreate(nl); + + if (z != BZNULL) { + /* + * set the sign of z such that the pointer n is unchanged yet + */ + + if (BnnIsZero(n, nl) == BN_TRUE) { + BzSetSign(z, BZ_ZERO); + } else { + BzSetSign(z, BZ_PLUS); + } + + for (i = 0; i < nl; ++i) { + BzSetDigit(z, i, n[i]); + } + } + + return z; +} + +/** + * BzToBigNum + * @param [in] z BigZ + * @param [in] nl BigNumLength + * @return BigNum + * @pre z != BZNULL. + * @pre nl != NULL. + */ +BigNum +BzToBigNum(const BigZ z, BigNumLength *nl) { + BigNum n; + BigNum m; + BigNumLength i; + size_t size; + + if (BzGetSign(z) == BZ_MINUS) { + return (BigNum)NULL; + } + + *nl = BzNumDigits(z); + size = ((size_t)(*nl) + 1) * sizeof(BigNumDigit); + + if ((n = (BigNum)(BzAlloc(size))) != NULL) { + *n = (BigNumDigit)*nl; /* set size */ + + for (i = 0, m = ++n; i < *nl; i++, m++) { + *m = BzGetDigit(z, i); + } + } + + return n; +} + +/* + * Logical oprations. + */ + +/** @cond */ +#define BZ_SIGN1 0x1 /* 01 */ +#define BZ_SIGN2 0x2 /* 10 */ +/** @endcond */ + +/** + * BzNot + * @param [in] z BigZ + * @return BigZ + * @pre z != BZNULL. + */ +BigZ +BzNot(const BigZ z) { + /* + * Negates the passed BigZ. + */ + + BigZ y; + + switch (BzGetSign(z)) { + case BZ_MINUS: + if ((y = BzCopy(z)) != BZNULL) { + BnnComplement2(BzToBn(y), BzNumDigits(y)); + BnnComplement(BzToBn(y), BzNumDigits(y)); + if (BnnIsZero(BzToBn(y), BzNumDigits(y)) == BN_TRUE) { + /* + * ~(-1) -> 0 + */ + BzSetSign(y, BZ_ZERO); + } else { + BzSetSign(y, BZ_PLUS); + } + } + break; + case BZ_ZERO: + y = BzFromInteger((BzInt)-1); + break; + case BZ_PLUS: + default: /* case BZ_PLUS: */ + if ((y = BzCopy(z)) != BZNULL) { + BnnComplement(BzToBn(y), BzNumDigits(y)); + BnnComplement2(BzToBn(y), BzNumDigits(y)); + BzSetSign(y, BZ_MINUS); + } + break; + } + + return y; +} + +/** + * BzAnd. + * Returns y & z. + * @param [in] y BigZ + * @param [in] z BigZ + * @return BigZ + * @pre y != BZNULL. + * @pre z != BZNULL. + */ +BigZ +BzAnd(const BigZ y, const BigZ z) { + /* + * Returns Y & Z. + */ + + BigZ n; + BigZ yy; + BigZ zz; + BigNumLength yl; + BigNumLength zl; + BigNumLength l; + unsigned int sign = 0; + + yl = BzNumDigits(y); + zl = BzNumDigits(z); + + if (BzGetSign(y) == BZ_MINUS) { + if ((yy = BzCopy(y)) == BZNULL) { + return BZNULL; + } + BnnComplement2(BzToBn(yy), yl); + sign |= BZ_SIGN1; + } else { + yy = y; + } + + if (BzGetSign(z) == BZ_MINUS) { + if ((zz = BzCopy(z)) == BZNULL) { + if ((sign & BZ_SIGN1) != 0) { + BzFree(yy); + } + return BZNULL; + } + BnnComplement2(BzToBn(zz), zl); + sign |= BZ_SIGN2; + } else { + zz = z; + } + + if (yl < zl) { + if ((n = BzCopy(zz)) != BZNULL) { + BzSetSign(n, BZ_PLUS); + if ((sign & BZ_SIGN1) == 0) { + for (l = (zl - 1); l >= yl; --l) { + BnnAndDigits(BzToBn(n) + l, BN_ZERO); + } + } + for (l = 0; l < yl; ++l) { + BnnAndDigits(BzToBn(n) + l, *(BzToBn(yy) + l)); + } + } + } else { + if ((n = BzCopy(yy)) != BZNULL) { + BzSetSign(n, BZ_PLUS); + if ((sign & BZ_SIGN2) == 0) { + for (l = (yl - 1); l >= zl; --l) { + BnnAndDigits(BzToBn(n) + l, BN_ZERO); + } + } + for (l = 0; l < zl; ++l) { + BnnAndDigits(BzToBn(n) + l, *(BzToBn(zz) + l)); + } + } + } + + if (n != BZNULL) { + if (BnnIsZero(BzToBn(n), BzNumDigits(n)) == BN_TRUE) { + BzSetSign(n, BZ_ZERO); + } else if (sign == (unsigned int)(BZ_SIGN1 | BZ_SIGN2)) { + BnnComplement2(BzToBn(n), BzNumDigits(n)); + BzSetSign(n, BZ_MINUS); + } + } + + /* + * Free copies. + */ + + if ((sign & BZ_SIGN1) != 0) { + BzFree(yy); + } + + if ((sign & BZ_SIGN2) != 0) { + BzFree(zz); + } + + return n; +} + +/** + * BzOr + * Returns y | z. + * @param [in] y BigZ + * @param [in] z BigZ + * @return BigZ + * @pre y != BZNULL. + * @pre z != BZNULL. + */ +BigZ +BzOr(const BigZ y, const BigZ z) { + BigZ n; + BigZ yy; + BigZ zz; + BigNumLength yl; + BigNumLength zl; + BigNumLength l; + unsigned int sign = 0; + + yl = BzNumDigits(y); + zl = BzNumDigits(z); + + if (BzGetSign(y) == BZ_MINUS) { + if ((yy = BzCopy(y)) == BZNULL) { + return BZNULL; + } + BnnComplement2(BzToBn(yy), yl); + sign |= BZ_SIGN1; + } else { + yy = y; + } + + if (BzGetSign(z) == BZ_MINUS) { + if ((zz = BzCopy(z)) == BZNULL) { + if ((sign & BZ_SIGN1) != 0) { + BzFree(yy); + } + return BZNULL; + } + BnnComplement2(BzToBn(zz), zl); + sign |= BZ_SIGN2; + } else { + zz = z; + } + + if (yl < zl) { + if ((n = BzCopy(zz)) != BZNULL) { + BzSetSign(n, BZ_PLUS); + if ((sign & BZ_SIGN1) != 0) { + for (l = (zl - 1); l >= yl; --l) { + BnnAndDigits(BzToBn(n) + l, BN_ZERO); + } + } + for (l = 0; l < yl; ++l) { + BnnOrDigits(BzToBn(n) + l, *(BzToBn(yy) + l)); + } + } + } else { + if ((n = BzCopy(yy)) != BZNULL) { + BzSetSign(n, BZ_PLUS); + if ((sign & BZ_SIGN2) != 0) { + for (l = (yl - 1); l >= zl; --l) { + BnnAndDigits(BzToBn(n) + l, BN_ZERO); + } + } + for (l = 0; l < zl; ++l) { + BnnOrDigits(BzToBn(n) + l, *(BzToBn(zz) + l)); + } + } + } + + if (n != BZNULL) { + if (BnnIsZero(BzToBn(n), BzNumDigits(n)) == BN_TRUE) { + BzSetSign(n, BZ_ZERO); + } else if (sign != 0) { + BnnComplement2(BzToBn(n), BzNumDigits(n)); + BzSetSign(n, BZ_MINUS); + } + } + + /* + * Free copies. + */ + + if ((sign & BZ_SIGN1) != 0) { + BzFree(yy); + } + + if ((sign & BZ_SIGN2) != 0) { + BzFree(zz); + } + + return n; +} + +/** + * BzXor. + * Returns y ^ z. + * @param [in] y BigZ + * @param [in] z BigZ + * @return BigZ + * @pre y != BZNULL. + * @pre z != BZNULL. + */ +BigZ +BzXor(const BigZ y, const BigZ z) { + /* + * Returns Y ^ Z. + */ + + BigZ n; + BigZ yy; + BigZ zz; + BigNumLength yl; + BigNumLength zl; + BigNumLength l; + int sign = 0; + + yl = BzNumDigits(y); + zl = BzNumDigits(z); + + if (BzGetSign(y) == BZ_MINUS) { + if ((yy = BzCopy(y)) == BZNULL) { + return BZNULL; + } + BnnComplement2(BzToBn(yy), yl); + sign |= BZ_SIGN1; + } else { + yy = y; + } + + if (BzGetSign(z) == BZ_MINUS) { + if ((zz = BzCopy(z)) == BZNULL) { + if ((sign & BZ_SIGN1) != 0) { + BzFree(yy); + } + return BZNULL; + } + BnnComplement2(BzToBn(zz), zl); + sign |= BZ_SIGN2; + } else { + zz = z; + } + + if (yl < zl) { + if ((n = BzCopy(zz)) != BZNULL) { + BzSetSign(n, BZ_PLUS); + if ((sign & BZ_SIGN1) != 0) { + for (l = (zl - 1); l >= yl; --l) { + BnnXorDigits(BzToBn(n) + l, BN_COMPLEMENT); + } + } + for (l = 0; l < yl; ++l) { + BnnXorDigits(BzToBn(n) + l, *(BzToBn(yy) + l)); + } + } + } else { + if ((n = BzCopy(yy)) != BZNULL) { + BzSetSign(n, BZ_PLUS); + if ((sign & BZ_SIGN2) != 0) { + for (l = (yl - 1); l >= zl; --l) { + BnnXorDigits(BzToBn(n) + l, BN_COMPLEMENT); + } + } + for (l = 0; l < zl; ++l) { + BnnXorDigits(BzToBn(n) + l, *(BzToBn(zz) + l)); + } + } + } + + if (n != BZNULL) { + if (BnnIsZero(BzToBn(n), BzNumDigits(n)) == BN_TRUE) { + BzSetSign(n, BZ_ZERO); + } else if (sign == BZ_SIGN1 || sign == BZ_SIGN2) { + BnnComplement2(BzToBn(n), BzNumDigits(n)); + BzSetSign(n, BZ_MINUS); + } + } + + /* + * Free copies. + */ + + if ((sign & BZ_SIGN1) != 0) { + BzFree(yy); + } + + if ((sign & BZ_SIGN2) != 0) { + BzFree(zz); + } + + return n; +} + +/** + * BzTestBit. + * Returns BN_TRUE iff bit is on (i.e. 2**bit is one). It assumes + * that bit is a non-negative integer. + * @param [in] bit BigNumLength + * @param [in] z BigZ + * @return BigNumBool + * @pre z != BZNULL. + */ +BigNumBool +BzTestBit(BigNumLength bit, const BigZ z) { + BigNumLength zl; + BigNumBool r; + + zl = (BigNumLength)(bit / BN_DIGIT_SIZE); + + if (zl >= BzNumDigits(z)) { + return (BzGetSign(z) == BZ_MINUS) ? BN_TRUE : BN_FALSE; + } + + bit = (bit % BN_DIGIT_SIZE); + + if (BzGetSign(z) == BZ_MINUS) { + const BigZ y = BzCopy(z); + BigNumLength yl; + + if (y == BZNULL) { + return BN_FALSE; + } + + yl = BzNumDigits(y); + + BnnComplement2(BzToBn(y), yl); + BzSetSign(y, BZ_PLUS); + r = (BigNumBool)(((*(BzToBn(y) + zl)) & (BN_ONE << bit)) != 0); + BzFree(y); + } else { + r = (BigNumBool)(((*(BzToBn(z) + zl)) & (BN_ONE << bit)) != 0); + } + + return r; +} + +/** + * BzBitCount. + * Returns the number of bits set in z. + * @param [in] z BigZ + * @return BigNumLength + * @pre z != BZNULL. + */ +BigNumLength +BzBitCount(const BigZ z) { + BigNumLength nl = (BigNumLength)0; + BigZ y; + + switch (BzGetSign(z)) { + case BZ_MINUS: + if ((y = BzCopy(z)) == BZNULL) { + return (BigNumLength)0; + } + BnnComplement2(BzToBn(y), BzNumDigits(y)); + BnnComplement(BzToBn(y), BzNumDigits(y)); + nl = BnnNumCount(BzToBn(y), BzNumDigits(y)); + BzFree(y); + break; + case BZ_ZERO: + nl = (BigNumLength)0; + break; + case BZ_PLUS: + nl = BnnNumCount(BzToBn(z), BzNumDigits(z)); + break; + } + + return nl; +} + +/* + * Simple logical equivalence rules. + */ + +/** + * BzNand. + * Returns ~(y & z). + * @param [in] x BigZ + * @param [in] y BigZ + * @return BigZ + * @pre x != BZNULL. + * @pre y != BZNULL. + */ +BigZ +BzNand(const BigZ x, const BigZ y) { + const BigZ tmp = BzAnd(x, y); + BigZ res = BzNot(tmp); + BzFree(tmp); + return res; +} + +/** + * BzNor. + * Returns ~(x | y). + * @param [in] x BigZ + * @param [in] y BigZ + * @return BigZ + * @pre x != BZNULL. + * @pre y != BZNULL. + */ +BigZ +BzNor(const BigZ x, const BigZ y) { + const BigZ tmp = BzOr(x, y); + BigZ res = BzNot(tmp); + BzFree(tmp); + return res; +} + +/** + * BzEqv. + * Returns ~(x ^ y). + * @param [in] x BigZ + * @param [in] y BigZ + * @return BigZ + * @pre x != BZNULL. + * @pre y != BZNULL. + */ +BigZ +BzEqv(const BigZ x, const BigZ y) { + const BigZ tmp = BzXor(x, y); + BigZ res = BzNot(tmp); + BzFree(tmp); + return res; +} + +/** + * BzAndC1. + * Returns ~x & y. + * @param [in] x BigZ + * @param [in] y BigZ + * @return BigZ + * @pre x != BZNULL. + * @pre y != BZNULL. + */ +BigZ +BzAndC1(const BigZ x, const BigZ y) { + const BigZ tmp = BzNot(x); + BigZ res = BzAnd(tmp, y); + BzFree(tmp); + return res; +} + +/** + * BzAndC2. + * Returns x & ~y. + * @param [in] x BigZ + * @param [in] y BigZ + * @return BigZ + * @pre x != BZNULL. + * @pre y != BZNULL. + */ +BigZ +BzAndC2(const BigZ x, const BigZ y) { + const BigZ tmp = BzNot(y); + BigZ res = BzAnd(x, tmp); + BzFree(tmp); + return res; +} + +/** + * BzOrC1. + * Returns ~x | y. + * @param [in] x BigZ + * @param [in] y BigZ + * @return BigZ + * @pre x != BZNULL. + * @pre y != BZNULL. + */ +BigZ +BzOrC1(const BigZ x, const BigZ y) { + const BigZ tmp = BzNot(x); + BigZ res = BzOr(tmp, y); + BzFree(tmp); + return res; +} + +/** + * BzOrC2. + * Returns x | ~y. + * @param [in] x BigZ + * @param [in] y BigZ + * @return BigZ + * @pre x != BZNULL. + * @pre y != BZNULL. + */ +BigZ +BzOrC2(const BigZ x, const BigZ y) { + const BigZ tmp = BzNot(y); + BigZ res = BzOr(x, y); + BzFree(tmp); + return res; +} + +/** + * BzAsh. + * Returns y << n. + * @param [in] y BigZ + * @param [in] n int + * @return BigZ + * @pre y != BZNULL. + */ +BigZ +BzAsh(const BigZ y, int n) { + BigZ z; + + if (n > 0) { + /* + * Create a copy + space for the shift. + */ + + if ((BzNumDigits(y) == (BigNumLength)1) + && (BzGetDigit(y, 0) == (BigNumDigit)1)) { + /* + * Optimize case where y == 1 (power of two). + */ + int bitnb = n + 1; + BigNumDigit digit; + BigNumLength zl; + + /* + * Compute the number of buckets. + */ + + zl = (BigNumLength)bitnb / BN_DIGIT_SIZE; + + if (((BigNumLength)bitnb % BN_DIGIT_SIZE) != 0) { + ++zl; + } + + /* + * Create a zeroed BigZ of zl buckets. + */ + + if ((z = BzCreate(zl)) == BZNULL) { + return z; + } + + /* + * Set the highest bit. + */ + + digit = ((BigNumDigit)1 << ((BzUInt)n % BN_DIGIT_SIZE)); + BzSetDigit(z, zl - 1, digit); + + /* + * Sign is the sign of y. + */ + BzSetSign(z, BzGetSign(y)); + return z; + } else { + int ll; + BigNumLength zl; + int len = (int)BN_DIGIT_SIZE - 1; + + zl = (BigNumLength)n / BN_DIGIT_SIZE; + + if (((BigNumLength)n % BN_DIGIT_SIZE) != 0) { + ++zl; + } + + zl += BzNumDigits(y); + + if ((z = BzCreate(zl)) == BZNULL) { + return z; + } + + BnnAssign(BzToBn(z), BzToBn(y), BzNumDigits(y)); + BzSetSign(z, BzGetSign(y)); + + /* + * Now do the shift by BN_DIGIT_SIZE increment. + */ + + for (ll = n; ll >= (int)BN_DIGIT_SIZE; ll -= len) { + (void)BnnShiftLeft(BzToBn(z), + zl, + (BigNumLength)len); + } + + (void)BnnShiftLeft(BzToBn(z), zl, (BigNumLength)ll); + } + } else if (n < 0) { + BigZ one; + BigZ d; + + if (y == BZNULL) { + return BZNULL; + } else if (BzGetSign(y) == BZ_ZERO) { + return BzCopy(y); + } + + if ((one = BzFromInteger((BzInt)1)) == BZNULL) { + return BZNULL; + } + + if ((d = BzAsh(one, -n)) == BZNULL) { + BzFree(one); + return BZNULL; + } + + z = BzFloor(y, d); + BzFree(d); + BzFree(one); + } else { + /* n == 0, strange but legal. */ + return BzCopy(y); + } + + return z; +} + +/** + * BzSqrt. + * @param [in] z BigZ + * @return BigZ + * @pre z != BZNULL. + */ +BigZ +BzSqrt(const BigZ z) { + BigNumLength n; + BigZ x; + BigZ two; + + if (BzGetSign(z) == BZ_ZERO) { + return BzFromInteger((BzInt)0); + } + + n = BzLength(z); + + if ((n & 1) != 0) { + /* + * n is odd. + */ + n = (n >> 1) + 1; + } else { + /* + * n is even. + */ + n = (n >> 1); + } + + { + const BigZ one = BzFromInteger((BzInt)1); + + x = BzAsh(one, (int)n); + BzFree(one); + } + + two = BzFromInteger((BzInt)2); + + for (;;) { + const BigZ y = BzFloor(z, x); + BigZ v; + + if (BzCompare(x, y) != BZ_GT) { + BzFree(y); + break; + } + + v = BzAdd(x, y); + + BzFree(x); + x = BzFloor(v, two); + BzFree(v); + BzFree(y); + } + + /* + * Free temporary two + */ + + BzFree(two); + + return x; +} + +/** + * BzLcm + * Returns lcm(y, z). + * @param [in] y BigZ + * @param [in] z BigZ + * @return BigZ + * @pre y != BZNULL. + * @pre z != BZNULL. + */ +BigZ +BzLcm(const BigZ y, const BigZ z) { + BigZ a; + BigZ b; + BigZ r; + + r = BZNULL; + + if ((a = BzMultiply(y, z)) != BZNULL) { + if (BzGetSign(a) == BZ_MINUS) { + BzSetSign(a, BZ_PLUS); + } + + if ((b = BzGcd(y, z)) != BZNULL) { + r = BzTruncate(a, b); + BzFree(b); + } + BzFree(a); + } + + return r; +} + +/** + * BzGcd. + * Returns gcd(y, z). + * @param [in] y BigZ + * @param [in] z BigZ + * @return BigZ + * @pre y != BZNULL. + * @pre z != BZNULL. + */ +BigZ +BzGcd(const BigZ y, const BigZ z) { + if (BzGetSign(y) == BZ_ZERO) { + return BzAbs(z); + } else if (BzGetSign(z) == BZ_ZERO) { + return BzAbs(y); + } else { + BigZ yc; + BigZ zc; + + if ((yc = BzAbs(y)) == BZNULL) { + /* a fresh copy failed */ + return BZNULL; + } + + if ((zc = BzAbs(z)) == BZNULL) { + /* a fresh copy failed */ + BzFree(yc); + return BZNULL; + } + + while (BzGetSign(zc) != BZ_ZERO) { + BigZ tmp = BzMod(yc, zc); + BzFree(yc); + + if (tmp == BZNULL) { + yc = BZNULL; + break; + } else { + yc = zc; + zc = tmp; + } + } + + BzFree(zc); + + return yc; + } +} + +static BigNumDigit BzInternalRandom(BzSeed *seed); + +#if defined(HAVE_STDINT_H) +typedef uint32_t BzUInt32; +#else +typedef unsigned int BzUInt32; +#endif + +#define BZ_XORSHIFT + +/** + * BzInternalRandom + * Compute an internal portable pseudo random value. + * @param [in] seed pointer to current BzSeed value. + * @return BigNumDigit in range [0 .. 32767] + */ +static BigNumDigit +BzInternalRandom(BzSeed *seed) { + BzUInt32 s; + + if (seed == (BzSeed *)0) { + /* + * Should not happen. Returns 0 instead of hanging. + */ + return 0; + } + + s = *seed; + +#if !defined(BZ_XORSHIFT) + /* + * http://en.wikipedia.org/wiki/Linear_congruential_generator + * a - 1 is divisible by all prime factors of m. + * (In our case m = 2^32, size of the int, so m has only one prime + * factor = 2) + * a - 1 is a multiple of 4 if m is a multiple of 4. + * (32768 is multiple of 4, and 1103515244 too) + */ + *seed = (BzSeed)(s * 1103515245 + 12345); + return (BigNumDigit)(((BzUInt32)*seed / 65536) % 32768); +#else + /* + * https://en.wikipedia.org/wiki/Xorshift + * Algorithm "xor" from p. 4 of Marsaglia, "Xorshift RNGs" + * We subtract 1 in order to be able to return 0. + */ + if (s == 0) { + s = 1; /* otherwise it will infinitely return 0. */ + } + + s ^= s << 13; + s ^= s >> 17; + s ^= s << 5; + *seed = s; + +#if defined(BZ_TRACE_MIN_MAX_SEED) + { + static BigNumDigit _min = ~0; + static BigNumDigit _max = 0; + + BigNumDigit val = (BigNumDigit)(((BzUInt32)*seed - 1) % 32768); + + if (val > _max) { + (void)printf("New _max = %x\n", val); + _max = val; + } + + if (val < _min) { + (void)printf("New _min = %x\n", val); + _min = val; + } + + return val; + } +#else + return (BigNumDigit)(((BzUInt32)*seed - 1) % 32768); +#endif +#endif +} + +/** + * BzRandom + * @param [in] n BigZ + * @param [in] seed BzSeed + * @return BigZ + * @pre n != BZNULL && seed != NULL + */ +BigZ +BzRandom(const BigZ n, BzSeed *seed) { + BigZ r; + BigZ res; + BigNumLength len; + BigNumLength i; + BigNumLength bytes; + + /* + * Algo: + * - create empty r as big as n + * - call BzInternalRandom() to replace all its bytes (8 bits a a time). + * - return r mod n. + * Assume any bit has an equiprobable [0-1] value. + */ + + if ((r = BzCreate(BzNumDigits(n))) == BZNULL) { + return BZNULL; + } + + len = BzLength(n); + bytes = len / BN_BYTE_SIZE; + + if (len % BN_BYTE_SIZE) { + ++bytes; + } + + for (i = 0; bytes > 0; ++i) { + BigNumLength j; + BigNumDigit *d = BzToBn(r) + i; + + for (j = 0; + (bytes > 0) && (j < (BigNumLength)sizeof(BigNumDigit)); + ++j) { + BigNumDigit chunk = (BzInternalRandom(seed) & 0xff); + *d += (chunk << (j * BN_BYTE_SIZE)); + --bytes; + } + } + + /* + * Call BzMod to ensure result is less than n. + */ + res = BzMod(r, n); + + BzFree(r); + + return res; +} + +/** + * BzPow. + * @param [in] base BigZ + * @param [in] exponent BzUInt + * @return BigZ + * @pre base != BZNULL. + */ +BigZ +BzPow(const BigZ base, BzUInt exponent) { + if (exponent == 0) { + /* + * Any nonzero number raised by the exponent 0 is 1 + */ + return BzFromInteger((BzInt)1); + } else { + BigZ x; + BigZ y; + + if ((x = BzMultiply(base, base)) == BZNULL) { + return BZNULL; + } + + y = BzPow(x, exponent / 2); + + BzFree(x); + + if (y == BZNULL) { + return BZNULL; + } else if ((exponent % (BzUInt)2) != 0) { + x = BzMultiply(y, base); + BzFree(y); + return x; + } else { + return y; + } + } +} + +/** + * BzModExp. + * The operation of modular exponentiation calculates the remainder + * when an integer b (the base) raised to the eth power (the + * exponent), be, is divided by a positive integer m (the modulus). In + * symbols, given base b, exponent e, and modulus m, the modular + * exponentiation c is: c = b**e mod m. From the definition of c, it + * follows that 0 <= c < m. + * + * Implementation uses the Right-to-left binary method + * https://en.wikipedia.org/wiki/Modular_exponentiation + * ~~~{.unparsed} + * function modular_pow(base, exponent, modulus) + * if (modulus == 1) + * return 0 + * result := 1 + * base := base mod modulus + * while exponent > 0 + * if (exponent mod 2 == 1) + * result := (result * base) mod modulus + * exponent := exponent >> 1 + * base := (base * base) mod modulus + * return result + * ~~~ + * @param [in] base BigZ + * @param [in] exponent BigZ + * @param [in] modulus BigZ + * @return BigZ + * @pre base != BZNULL. + * @pre exponent != BZNULL. + */ +BigZ +BzModExp(const BigZ base, const BigZ exponent, const BigZ modulus) { + BigZ result; + BigZ mod; + BigZ expnt; + BigZ b; + int neg; + + if ((result = BzFromInteger((BzInt)1)) == BZNULL) { + return BZNULL; + } + + switch (BzGetSign(exponent)) { + case BZ_ZERO: + /* + * exponent == 0, (base ** 0) == 1, two cases to consider: + */ + if (BzCompare(modulus, result) == BZ_EQ) { + /* + * modulus == 1 => 0 + */ + BnnSetToZero(BzToBn(result), (BigNumLength)1); + BzSetSign(result, BZ_ZERO); + return result; + } else if (BzGetSign(modulus) == BZ_MINUS) { + /* + * modulus < 0 => modulus + 1 + */ + const BigZ tmp = BzAdd(modulus, result); + BzFree(result); + return tmp; + } else { + /* + * modulus > 1 => 1 + */ + return result; + } + case BZ_MINUS: + /* + * Negative exponent is not supported. + */ + BzFree(result); + return BZNULL; + case BZ_PLUS: + default: + if (BzGetSign(modulus) == BZ_PLUS) { + /* + * modulus is positive, don't need to copy it. + */ + neg = 0; + mod = modulus; + } else { + /* + * Make a copy of modulus as positive value. + */ + + if ((mod = BzNegate(modulus)) == BZNULL) { + BzFree(result); + return BZNULL; + } + + neg = 1; + } + + /* + * Copy base as it will be modified. + */ + + if ((b = BzCopy(base)) == BZNULL) { + BzFreeIf(neg, mod); + BzFree(result); + return BZNULL; + } + + /* + * Copy exponent as it will be modified. + */ + + if ((expnt = BzCopy(exponent)) == BZNULL) { + BzFreeIf(neg, mod); + BzFree(b); + BzFree(result); + return BZNULL; + } + + while (BzGetSign(expnt) == BZ_PLUS) { + BigZ tmp; + + if (BzIsOdd(expnt)) { + tmp = BzMultiply(result, b); + BzFree(result); + if (tmp == BZNULL) { + BzFreeIf(neg, mod); + BzFree(expnt); + BzFree(b); + return BZNULL; + } + result = BzMod(tmp, mod); + BzFree(tmp); + if (result == BZNULL) { + BzFreeIf(neg, mod); + BzFree(expnt); + BzFree(b); + return BZNULL; + } + } + /* + * expnt = expnt >> 1; + */ + tmp = BzAsh(expnt, -1); + BzFree(expnt); + expnt = tmp; + tmp = BzMultiply(b, b); + BzFree(b); + if (tmp == BZNULL) { + BzFreeIf(neg, mod); + BzFree(expnt); + BzFree(result); + return BZNULL; + } + b = BzMod(tmp, mod); + BzFree(tmp); + if (b == BZNULL) { + BzFreeIf(neg, mod); + BzFree(expnt); + BzFree(result); + return BZNULL; + } + } + + BzFree(expnt); + BzFree(b); + + if (neg) { + /* + * Modulus is negative. Adjust value. + */ + BigZ tmp; + BzFree(mod); + tmp = BzMod(result, modulus); + BzFree(result); + return tmp; + } else { + return result; + } + } +} + +/** + * BzHash + * The hash value is computed from at most 4 buckets (frist two and last two). + * @param [in] z BigNum + * @return BzUint hash value + */ +BzUInt +BzHash(const BigZ z) { + BzUInt res = 0; + + if (z == BZNULL) { + return res; + } else { + BigNumLength nl; + + switch (BzGetSign(z)) { + case BZ_MINUS : + res = (BzUInt)~0; + break; + case BZ_PLUS : + res = (BzUInt)0; + break; + case BZ_ZERO : + return res; + } + + nl = BzGetSize(z); + + if (nl >= 4) { + /* + * Compute hash value with only the two first and + * two last buckets. + */ + res += BzGetDigit(z, 0) + + BzGetDigit(z, 1) + + BzGetDigit(z, nl - 2) + + BzGetDigit(z, nl - 1); + } else { + /* + * Compute hash value with all buckets. + */ + BigNumLength i; + + for (i = 0; i < nl; ++i) { + res += BzGetDigit(z, i); + } + } + + return res; + } +} diff --git a/modules/bigz/native/bigz.h b/modules/bigz/native/bigz.h new file mode 100644 index 00000000..74729554 --- /dev/null +++ b/modules/bigz/native/bigz.h @@ -0,0 +1,342 @@ +/* + * Simplified BSD License + * + * Copyright (c) 1988-1989, Digital Equipment Corporation & INRIA. + * Copyright (c) 1992-2023, Eligis + * All rights reserved. + * + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions are + * met: + * + * o Redistributions of source code must retain the above copyright + * notice, this list of conditions and the following disclaimer. + * o Redistributions in binary form must reproduce the above copyright + * notice, this list of conditions and the following disclaimer in the + * documentation and/or other materials provided with the distribution. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS + * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT + * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR + * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT + * HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, + * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT + * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, + * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY + * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT + * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE + * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + */ + +/** + * @file bigz.h + * @brief Types and structures for clients of BigZ. + * @version 2.1.0 + * @copyright Digital Equipment Corporation & INRIA, 1988-1989. + * @copyright Eligis, 1992-2023 + * @authors B. Serpette + * @authors J. Vuillemin + * @authors J.C. Hervé + * @authors C. Jullien + * $Revision: 395 $ + * $Date: 2023-02-07 06:53:31 +0000 (Tue, 07 Feb 2023) $ + */ + +#if !defined(__BIGZ_H) +#define __BIGZ_H + +#if !defined(__BIGN_H) +#include "./bign.h" +#endif + +#if defined(__cplusplus) && !defined(CPP_MODULE) +extern "C" { +#endif + +#include +#include + +/** @cond */ +#define BZ_PURE_FUNCTION BN_PURE_FUNCTION +#define BZ_CONST_FUNCTION BN_CONST_FUNCTION +/** @endcond */ + +/** + * BigZ version MM.mm.pp. + */ +#define BZ_VERSION "2.1.0" + +/** + * BigZ sign + */ +typedef enum { + /** Negative */ + BZ_MINUS = -1, + /** Zero */ + BZ_ZERO = 0, + /** Positive */ + BZ_PLUS = 1 +} BzSign; + +/** + * BigZ compare result + */ +typedef enum { + /** Less than value (-1). */ + BZ_LT = BN_LT, + /** Equal than == value (0). */ + BZ_EQ = BN_EQ, + /** Greater than comparison value (1). */ + BZ_GT = BN_GT +} BzCmp; + +/** @cond */ +typedef enum { + BZ_UNTIL_END = 0, + BZ_UNTIL_INVALID = 1, + BZ_UNTIL_SLASH = 2, + BZ_UNTIL_SPACE = 3 +} BzStrFlag; + +/** + * BigZ number header. + */ +typedef struct { + BigNumLength Size; + BzSign Sign; +} BigZHeader; + +/* + * define a dummy positive value to declare a Digits vector. + * BigZ is allocated with required size. Choose a rather large value + * to prevent 'smart' compilers to exchange fields. + * C11/C++11 compilers and above are supposed to support 'dynamic arrays' + */ + +#if (((defined(__STDC_VERSION__) && (__STDC_VERSION__ >= 201112L)) \ + || (defined(__cplusplus) && (__cplusplus >= 201103L))) \ + && !defined(_MSC_VER)) +#define BZ_DUMMY_SIZE 0 +#else +#define BZ_DUMMY_SIZE 32 +#endif + +typedef struct { + BigZHeader Header; + /* + * Digit vector should be the last field to allow allocation + * of the real size (BZ_DUMMY_SIZE is never used). + */ +#if (BZ_DUMMY_SIZE == 0) + BigNumDigit Digits[]; +#else + BigNumDigit Digits[BZ_DUMMY_SIZE]; +#endif +} BigZStruct; + +typedef BigZStruct * __BigZ; +typedef const BigZStruct * __CBigZ; + +/** + * BzToString and related functions add sign only when number is negative ('-'). + */ +#define BZ_DEFAULT_SIGN 0 + +/** + * Force BzToString and related functions to always add a sign [+|-] when + * number is not 0. + */ +#define BZ_FORCE_SIGN 1 + +/* + * macros of bigz.c + */ + +#if !defined(BZ_BIGNUM_TYPE) +#define BZ_BIGNUM_TYPE +typedef __BigZ BigZ; +typedef __CBigZ CBigZ; +#endif + +#if !defined(BZ_CHAR_TYPE) +#define BZ_CHAR_TYPE +typedef char BzChar; +#endif + +#if defined(BZ_EXPERIMENTAL_128BIT) +/* + * No really need to use those types. The most important thing to have is + * bucket size of 128bit as provided by BN_EXPERIMENTAL_128BIT. + */ +#define BZ_INT_TYPE +typedef __int128_t BzInt; +#define BZ_UINT_TYPE +typedef __uint128_t BzUInt; +#endif + +#if !defined(BZ_INT_TYPE) +#define BZ_INT_TYPE +#if (defined(HAVE_STDINT_H) && defined(SIZEOF_VOID_P) && (SIZEOF_VOID_P >= 8)) +typedef int64_t BzInt; +#else +typedef int BzInt; +#endif +#endif + +#if !defined(BZ_UINT_TYPE) +#define BZ_UINT_TYPE +#if (defined(HAVE_STDINT_H) && (SIZEOF_VOID_P >= 8)) +typedef uint64_t BzUInt; +#else +typedef unsigned int BzUInt; +#endif +#endif + +#if !defined(BZ_DBL_MAX) && defined(LDBL_MAX) +#define BZ_DBL_MAX LDBL_MAX +typedef long double BzLDouble; +#endif + +#if !defined(BZ_DBL_MAX) && defined(DBL_MAX) +#define BZ_DBL_MAX DBL_MAX +typedef double BzLDouble; +#endif + +#if !defined(BZ_DBL_MAX) +#define BZ_DBL_MAX ((double)+1.7976931348623158e+308) +typedef double BzLDouble; +#endif + +/* + * Random seed type, by contract it must be an unsigned int. + */ +typedef unsigned int BzSeed; + +#define BZ_OPTIMIZE_PRINT + +#if !defined(BZ_BUCKET_SIZE) +#if (defined(SIZEOF_LONG) && (SIZEOF_LONG == 8)) +#define BZ_BUCKET_SIZE 64 +#else +#define BZ_BUCKET_SIZE 32 +#endif +#endif +/** @endcond */ + +#if !defined(__EXTERNAL_BIGZ_MEMORY) +/** + * User overloadable macro that gets native BigZ implementation from + * a high level object (for example managed in C++ or Lisp). + */ +#define __toBzObj(z) ((__BigZ)z) +/** + * NULL BigZ. + */ +#define BZNULL ((BigZ)0) +/** + * User overloadable macro called to allocate size bytes to store a BigZ. + */ +#define BzAlloc(size) malloc(size) +/** + * User overloadable macro called to free a BigZ allocated by BzAlloc. + */ +#define BzFree(z) free(z) +/** + * User overloadable macro called to allocate size BzChar that represent + * a BzString. + */ +#define BzStringAlloc(size) malloc(size * sizeof(BzChar)) +/** + * User overloadable macro called to free a buffer allocated by BzStringAlloc. + */ +#define BzFreeString(s) free(s) +#endif + +/** @cond */ +#define BzGetSize(z) (__toBzObj(z)->Header.Size) +#define BzGetSign(z) (__toBzObj(z)->Header.Sign) +#define BzGetDigit(z, n) (__toBzObj(z)->Digits[n]) +#define BzToBn(z) (__toBzObj(z)->Digits) +#define BzSetSize(z, s) (__toBzObj(z)->Header.Size = (s)) +#define BzSetSign(z, s) (__toBzObj(z)->Header.Sign = (s)) +#define BzSetDigit(z, n, v) (__toBzObj(z)->Digits[n] = (v)) +/** @endcond */ + +/* + * functions of bigz.c + */ + +extern const char * BzVersion(void) BZ_CONST_FUNCTION; +extern BigZ BzCreate(BigNumLength Size); +extern BigNumLength BzNumDigits(const BigZ z) BZ_PURE_FUNCTION; +extern BigNumLength BzLength(const BigZ z) BZ_PURE_FUNCTION; +extern BigZ BzCopy(const BigZ z); +extern BigZ BzNegate(const BigZ z); +extern BigZ BzAbs(const BigZ z); +extern BzCmp BzCompare(const BigZ y, const BigZ z) BZ_PURE_FUNCTION; +extern BigZ BzAdd(const BigZ y, const BigZ z); +extern BigZ BzSubtract(const BigZ y, const BigZ z); +extern BigZ BzMultiply(const BigZ y, const BigZ z); +extern BigZ BzDivide(const BigZ y, const BigZ z, BigZ *r); +extern BigZ BzDiv(const BigZ y, const BigZ z); +extern BigZ BzTruncate(const BigZ y, const BigZ z); +extern BigZ BzFloor(const BigZ y, const BigZ z); +extern BigZ BzCeiling(const BigZ y, const BigZ z); +extern BigZ BzRound(const BigZ y, const BigZ z); +extern BigZ BzMod(const BigZ y, const BigZ z); +extern BigZ BzRem(const BigZ y, const BigZ z); +extern BigZ BzPow(const BigZ base, BzUInt exponent); +extern BigNumBool BzIsEven(const BigZ y) BZ_PURE_FUNCTION; +extern BigNumBool BzIsOdd(const BigZ y) BZ_PURE_FUNCTION; +extern BzChar * BzToString(const BigZ z, BigNumDigit base, int sign); +extern size_t BzStrLen(const BzChar *s) BN_PURE_FUNCTION; +extern BzChar * BzToStringBuffer(const BigZ z, BigNumDigit base, int sign, /*@null@*/ BzChar * const buf, /*@null@*/ size_t *len); +extern BzChar * BzToStringBufferExt(const BigZ z, BigNumDigit base, int sign, /*@null@*/ BzChar * const buf, /*@null@*/ size_t *len, /*@null@*/ size_t *slen); +extern BigZ BzFromStringLen(const BzChar *s, size_t len, BigNumDigit base, BzStrFlag flag); +//extern BigZ BzFromString(const BzChar *s, BigNumDigit base, BzStrFlag flag); +extern BigZ BzFromString(const BzChar *s, unsigned long base, unsigned int flag); +extern BigZ BzFromInteger(BzInt i); +extern BigZ BzFromUnsignedInteger(BzUInt i); +extern BzInt BzToInteger(const BigZ z) BZ_PURE_FUNCTION; +extern double BzToDouble(const BigZ z) BZ_PURE_FUNCTION; +extern BzLDouble BzToLongDouble(const BigZ z) BZ_PURE_FUNCTION; +extern int BzToIntegerPointer(const BigZ z, BzInt *p); +extern BzUInt BzToUnsignedInteger(const BigZ z) BZ_PURE_FUNCTION; +extern int BzToUnsignedIntegerPointer(const BigZ z, BzUInt *p); +extern BigZ BzFromBigNum(const BigNum n, BigNumLength nl); +extern BigNum BzToBigNum(const BigZ z, BigNumLength *nl); +extern BigNumBool BzTestBit(BigNumLength bit, const BigZ z); +extern BigNumLength BzBitCount(const BigZ z); +extern BigZ BzNot(const BigZ z); +extern BigZ BzAnd(const BigZ y, const BigZ z); +extern BigZ BzOr(const BigZ y, const BigZ z); +extern BigZ BzXor(const BigZ y, const BigZ z); +extern BigZ BzNand(const BigZ x, const BigZ y); +extern BigZ BzNor(const BigZ x, const BigZ y); +extern BigZ BzEqv(const BigZ x, const BigZ y); +extern BigZ BzAndC1(const BigZ x, const BigZ y); +extern BigZ BzAndC2(const BigZ x, const BigZ y); +extern BigZ BzOrC1(const BigZ x, const BigZ y); +extern BigZ BzOrC2(const BigZ x, const BigZ y); +extern BigZ BzAsh(const BigZ y, int n); +extern BigZ BzSqrt(const BigZ z); +extern BigZ BzLcm(const BigZ y, const BigZ z); +extern BigZ BzGcd(const BigZ y, const BigZ z); +extern BigZ BzRandom(const BigZ n, BzSeed *seed); +extern BigZ BzModExp(const BigZ base, const BigZ exponent, const BigZ modulus); +extern BzUInt BzHash(const BigZ z); + +/* +#define BZ_DEBUG +*/ + +#if defined(BZ_DEBUG) +extern void BnDebug(const char *m, const BzChar *bzstr, const BigNum n, BigNumLength nl, BzSign sign); +extern void BzDebug(const char *m, const BigZ y); +#endif + +#if defined(__cplusplus) && !defined(CPP_MODULE) +} +#endif + +#endif /* __BIGZ_H */ diff --git a/modules/bigz/native/bzf.c b/modules/bigz/native/bzf.c new file mode 100644 index 00000000..dc225fd2 --- /dev/null +++ b/modules/bigz/native/bzf.c @@ -0,0 +1,73 @@ +/* + * $Id: bzf.c 392 2023-02-04 06:21:06Z cjullien $ + */ + +/* + * Copyright (c) 1988-1989, Digital Equipment Corporation & INRIA. + * Copyright (c) 1992-2023, Eligis + * All rights reserved. + * + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions are + * met: + * + * o Redistributions of source code must retain the above copyright + * notice, this list of conditions and the following disclaimer. + * o Redistributions in binary form must reproduce the above copyright + * notice, this list of conditions and the following disclaimer in the + * documentation and/or other materials provided with the distribution. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS + * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT + * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR + * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT + * HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, + * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT + * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, + * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY + * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT + * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE + * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + */ + +/** + * @file bzf.c + * @brief Miscellaneous functions built on top of BigZ. + * @version 2.1.0 + * @copyright Digital Equipment Corporation & INRIA, 1988-1989. + * @copyright Eligis, 1992-2023 + * @authors B. Serpette + * @authors J. Vuillemin + * @authors J.C. Hervé + * @authors C. Jullien + * $Revision: 392 $ + * $Date: 2023-02-04 06:21:06 +0000 (Sat, 04 Feb 2023) $ + */ + +#if !defined(__BZF_H) +#include "./bzf.h" +#endif + +BigZ +BzFactorial(BigZ z) { + /* + * Returns Z! + * Assumes Z < Base. + */ + + BigZ f; + BigNumDigit zval; + int fl = 1; + + zval = BnnGetDigit(BzToBn(z)); + f = BzCreate(zval+1); + BnnSetDigit(BzToBn(f), 1); + BzSetSign(f, BzGetSign(z)); + + while (zval-- > 1) { + BnnMultiplyDigit(BzToBn(f), fl+1, BzToBn(f), fl, zval); + fl = BnnNumDigits(BzToBn(f), fl+1); + } + + return (f); +} diff --git a/modules/bigz/native/bzf.h b/modules/bigz/native/bzf.h new file mode 100644 index 00000000..6441a596 --- /dev/null +++ b/modules/bigz/native/bzf.h @@ -0,0 +1,62 @@ +/* + * $Id: bzf.h 392 2023-02-04 06:21:06Z cjullien $ + */ + +/* + * Copyright (c) 1988-1989, Digital Equipment Corporation & INRIA. + * Copyright (c) 1992-2023, Eligis + * All rights reserved. + * + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions are + * met: + * + * o Redistributions of source code must retain the above copyright + * notice, this list of conditions and the following disclaimer. + * o Redistributions in binary form must reproduce the above copyright + * notice, this list of conditions and the following disclaimer in the + * documentation and/or other materials provided with the distribution. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS + * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT + * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR + * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT + * HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, + * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT + * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, + * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY + * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT + * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE + * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + */ + +/** + * @file bzf.h + * @brief Miscellaneous functions built on top of BigZ. + * @version 2.1.0 + * @copyright Digital Equipment Corporation & INRIA, 1988-1989. + * @copyright Eligis, 1992-2023 + * @authors B. Serpette + * @authors J. Vuillemin + * @authors J.C. Hervé + * @authors C. Jullien + * $Revision: 392 $ + * $Date: 2023-02-04 06:21:06 +0000 (Sat, 04 Feb 2023) $ + */ + +#if !defined(__BZF_H) +#define __BZF_H + +#include "./bigz.h" + +#if defined(__cplusplus) && !defined(CPP_MODULE) +extern "C" { +#endif + +extern BigZ BzFactorial(BigZ z); + +#if defined(__cplusplus) && !defined(CPP_MODULE) +} +#endif + +#endif /* __BZF_H */ diff --git a/modules/bigz/structs.wx b/modules/bigz/structs.wx new file mode 100644 index 00000000..a6bd21c9 --- /dev/null +++ b/modules/bigz/structs.wx @@ -0,0 +1,88 @@ + +Namespace bigz + +Const BzZero := New Bz( 0 ) +Const BzOne := New Bz( 1 ) +Const BzTwo := New Bz( 2 ) +Const BzFive := New Bz( 5 ) +Const BzTen := New Bz( 10 ) + +Struct Bz + Field bz:BigZ + + Method New( x:Int ) + Self.bz = BzFromInteger( x ) + End + + Method New( x:BigZ ) + Self.bz = BzCopy( x ) + End + + Method New( x:String ) + Self.bz = BzFromString( x,Cast( 10 ),Cast( 0 ) ) + End + + Operator+:Bz( y:Bz ) + Local t := New Bz( Self.bz ) + t.bz = BzAdd( t.bz,y.bz ) + Return t + End + + Operator-:Bz( y:Bz ) + Local t := New Bz( Self.bz ) + t.bz = BzSubtract( t.bz,y.bz ) + Return t + End + + Operator*:Bz( y:Bz ) + Local t := New Bz( Self.bz ) + t.bz = BzMultiply( t.bz,y.bz ) + Return t + End + + Operator/:Bz( y:Bz ) + Local t := New Bz( Self.bz ) + t.bz = BzDiv( t.bz,y.bz ) + Return t + End + + Operator Mod:Bz( y:Bz ) + Local r := BzZero + r.bz = BzMod( Self.bz,y.bz ) + Return r + End + + Operator+:Bz( y:Int ) + Local t := New Bz( Self.bz ) + t.bz = BzAdd( t.bz,BzFromInteger( y ) ) + Return t + End + + Operator-:Bz( y:Int ) + Local t := New Bz( Self.bz ) + t.bz = BzSubtract( t.bz,BzFromInteger( y ) ) + Return t + End + + Operator*:Bz( y:Int ) + Local t := New Bz( Self.bz ) + t.bz = BzMultiply( t.bz,BzFromInteger( y ) ) + Return t + End + + Operator/:Bz( y:Int ) + Local t := New Bz( Self.bz ) + t.bz = BzDiv( t.bz,BzFromInteger( y ) ) + Return t + End + + Operator Mod:Bz( y:Int ) + Local r := BzZero + r.bz = BzMod( Self.bz,BzFromInteger( y ) ) + Return r + End + + Operator To:String() + Return BzToString( Self.bz,10,0 ) + End +End diff --git a/modules/bigz/tests/00-example.wx b/modules/bigz/tests/00-example.wx new file mode 100644 index 00000000..172fd516 --- /dev/null +++ b/modules/bigz/tests/00-example.wx @@ -0,0 +1,40 @@ + +Namespace myapp + +#Import "" +#Import "" + +Using std.. +Using bigz.. + +Function Fac:String(num:Int) + Local tmp:=BzFromInteger(1) + For Local i := 1 To num + tmp = BzMultiply(tmp, BzFromInteger(i)) + Next + Return BzToString(tmp,10,0) +End + +Function Main() + + Print BzVersion() + Local b1:=BzFromInteger(987654321) + Local b2:=BzFromInteger(123456789) + Local fac:=BzFromInteger(500) + Local a:=BzAdd(b1,b2) + Local s:=BzSubtract(b1,b2) + Local m:=BzMultiply(b1,b2) + Local d:=BzDiv(b1,b2) + Local r:=BzRem(b1,b2) + Local f:=BzFactorial(fac) + + Print "Suma:"+BzToString(a,10,0) + Print "Resta:"+BzToString(s,10,0) + Print "Multiplicacion:"+BzToString(m,10,0) + Print "Division:"+BzToString(d,10,0) + Print "Resto:"+BzToString(r,10,0) + Print "Factorial(500)="+BzToString(f,10,0) + Local buf:=Fac(1000) + Print buf + Print buf.Length +End \ No newline at end of file diff --git a/modules/bigz/tests/01-example.wx b/modules/bigz/tests/01-example.wx new file mode 100644 index 00000000..0dde1b88 --- /dev/null +++ b/modules/bigz/tests/01-example.wx @@ -0,0 +1,33 @@ + +Namespace myapp + +#Import "" +#Import "" + +Using std.. +Using bigz.. + +Function Main() + + Print BzVersion() + + Local z1:=New Bz("9876543210987654321") + Local z2:=New Bz("1234567890123456789") + Local z3:=z1+z2 + Print z1 + Print z2 + Print z1+z2 + Print z1-z2 + Print z1*z2 + Print z1/z2 + Print z1 Mod z2 + z1.bz = BzFromInteger(1) + For Local i:=2 To 20 + z1*=i + Next + Print z1 + Local f:=New Bz(0) + f.bz = BzFactorial(BzFromInteger(20)) + Print f Mod 100000000 + +End \ No newline at end of file