ch24aok2.c

稀疏矩阵、链表、图、队列、二叉树、多叉树、排序、遗传算法等的实现

/*
*
*  This code is part of Chapter 24 of "C Unleashed", and is to be     
*  read in conjunction with that chapter. This code deals only with 
*  very large/precise numbers. These may be integers or have fractional
*  parts (decimal places).
*
*  Code by Ian D. K. Kelly
*
*/

/* Ch24AOK2.c
* 
*  Arith - generalized arithmetic routines
* 
*  Copyright (C) 1999  Ian D. K. Kelly, 
*                      idkk Consultancy Ltd.
*                      Macmillan Computer Publishing 
* 
*  This program is free software; you can redistribute it and/or modify 
*  it under the terms of the GNU General Public License as published by 
*  the Free Software Foundation; either version 2 of the License, or 
*  (at your option) any later version. 
* 
*  This program is distributed in the hope that it will be useful, 
*  but WITHOUT ANY WARRANTY; without even the implied warranty of 
*  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the 
*  GNU General Public License for more details. 
* 
*  You should have received a copy of the GNU General Public License 
*  along with this program; if not, write to the Free Software 
*  Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. 
* 
*  Ian Kelly may be contacted at [e-mail]
* 
*/ 
            
#include <stdio.h>
#include <stdlib.h>
#include <limits.h>

/*
*
* These are sample routines for the basic four arithmetic operators
* for extremely large, extremely precise numbers. 
* All very large numbers are assumed to be stored in arrays of INT
* (which may be "int" or "long int", depending on environmental precisions).  
* The first element of each of these arrays is a counter of the total
* number of further elements that exist in that array. The second element
* of each array is a counter of the number of elements within that array
* that are to be considered as being to the right of the decimal point.
* This second value must be no greater than the first value. The value of  
* the number is stored from left to right, from the most significant
* to the least significant. That is, the lowest subscripts have the  
* largest part of the number.                                        
*
* Each routine, after computing its answer, tries to normalize it. In
* the context of this example code, "normalized" means "having the   
* fewest possible number of leading zeros". If there are any places 
* before the point, then these are the ones that are reduced, and NOT
* the places after the point. The trailing zeros after the point are
* also reduced to the minimum.
* In no case is the length of a number zero: it has at least one "value"
* word. Hence the value zero is:
*       2, 0, 0
* the value one would be
*       2, 0, 1
* the value "point three" (to the base of representation) would be
*       2, 1, 3
* and so on. If the base of representation is 10000 (ten thousand - see
* description further on) then "pi" to 36 places (truncated, not rounded)
* might be represented as
*       11, 9, 3, 1415, 9265, 3589, 7932, 3846, 2643, 3832, 7950, 2884
* and the speed of light in centimeters per second could be
*        5, 1, 2, 9979, 2459, 5000
* which is "299,792,459.5". Avogadro's constant (6.022 E+23) would be
*        7, 0, 6022, 0000, 0000, 0000, 0000, 0000 
*and the proton rest mass in kg. (1.672 E-27) would be
*       10, 9, 0000, 0000, 0000, 0000, 0000, 0000, 0001, 6726, 1411
* Hence the rather meaningless number of "Avogadro plus proton rest mass"
* would be
*       16, 9, 6022, 0000, 0000, 0000, 0000, 0000,
*       0000, 0000, 0000, 0000, 0000, 0000, 0001, 6726, 1411    
*/

/*
*
* The number representation can be thought of as using an extremely   
* large base. For our ordinary counting we use base 10 (ten), if we  
* are using octal, then the base is 8, for hexadecimal 16, and so on.
* We set here a base which is large enough to ensure that we can hold
* very large numbers, but small enough to ensure that the addition of
* two numbers does NOT overflow a word. That is, (2 * BASE) < MAX_INT
* The other restriction is that the number base must not be greater 
* than the square root of MAX_INT, as during the divide process we
* have to multiply the remainder by the base, and add on the next
* "digit" (word). This is the tighter of the two restrictions.
*
* We do not HAVE to choose a base which is a power of ten, but it    
* makes the printing of the answers so much simpler if we do.        
* The base chosen, for 32-bit int, is ten thousand, which is a power 
* of ten, and whose square root is 100 (one hundred). This gives exactly
* four decimal digits to each "digit" (word, element) of the stored
* very large numbers.                
*
*/

#if (INT_MAX>100000000) 

/* INT_MAX is large enough to use the definitions described above:    */
#define BASE 10000
#define SQRT_BASE 100
typedef int INT;
#define DIGITS_IN_BASE 4

#else

/*
*
* INT_MAX is not large enough, therefore we have to use "long int" as
* our base element. It might be possible to consider the 64-bit upper
* limit of 9,223,372,036,854,775,807 ( greater than 1E18), hence the
* BASE could be 100000000 (1E8, ten to the eight), and SQRT_BASE 10000 
* (1E4, ten to the four, ten thousand). This depends upon the particular
* machine, and the particular environment. For this sample code, however,
* we stick to the same limits as above: ten thousand, and one hundred.
*
*/
#define BASE 10000
#define SQRT_BASE 100
typedef long int INT;
#define DIGITS_IN_BASE 4

#endif

/*
* The "weak" naming convention addopted in these routines is:
*   each function name contains one Verb (capitalised) that indicates
*        the main action of that routine
*   there is a prefix (not capitalised) to that Verb which indicates:
*        a        operation on a set of arrays
*        casc     cascade: operation on one element, plus an array
*        one      operation on only single items
*        two      operation on one double and one single item
*        internal like an "a" prefix, but no allocations are made
*    variables may have "p", "pp", "i" prefixes indicating "pointer",
*        "pointer to pointer", "integer" and so on.
*
* Most routines whose names begin "a" or "casc" allocate space for their
* answers - except for "aCompare" and "aAbsCompare", which merely compare
* two arrays, and "aNormalize", which MAY allocate an answer array, but also 
* de-allocates the incoming array.
*/

int pairMultiply  ( INT iOne,   INT iTwo,   INT * pAnswer);
int cascMultiply  ( INT * iOne, INT * aTwo, INT ** aAnswer );
int aMultiply     ( INT * aOne, INT * aTwo, INT ** aAnswer);

int oneDivide  ( INT * iNumerator, INT * iDenominator, INT * iAnswerOD );
int pairDivide ( INT * iNumerator, INT * iDenominator, INT * iAnswerPD );
int cascDivide ( INT * aNumerator, INT * iDenominator, 
                 int places,       INT ** aAnswer );
int aDivide ( INT * aNumerator, INT * aDenominator, 
              int places,       INT ** aAnswer );
int aDivNormalize ( INT * pinNumerator, int places, INT ** ppoutNumerator );

int aAdd             ( INT * aOne, INT * aTwo, INT ** aAnswer );
int internalAdd      ( INT * aOne, INT * aTwo, INT ** aAnswer );
int aSubtract        ( INT * aOne, INT * aTwo, INT ** aAnswer );
int internalSubtract ( INT * aOne, INT * aTwo, INT ** aAnswer );

int aCompare    ( INT * aOne,       INT * aTwo );
int aAbsCompare ( INT * aOne,       INT * aTwo );
int aNormalize  ( INT ** aUnNormal, INT ** aNormal );
int aAllocate   ( int iCount,       INT ** aAnswer );

/* Multiply two INTs                                                  */
int pairMultiply ( INT iOne, INT iTwo, INT * pAnswer)
{
    /* This routine does not perform any internal allocations         */

    /*
    * Multiply two INTs giving an array of two INTs. The array is    
    * needed, because the answer (the product) may well have many more
    * digits than the incoming multiplicands.
    * This routine assumes that both input arguments are positive.
    */
    INT iStatus = EXIT_SUCCESS;
    INT iOneTop = 0;
    INT iOneBot = 0;
    INT iTwoTop = 0;
    INT iTwoBot = 0;
    INT iAnsTop = 0;
    INT iAnsBot = 0;
 
    
    /*
    * Split the incoming arguments up by the square root of the base,
    * so that we can never exceed the size of an INT during calculation:
    */
    if (iOne>SQRT_BASE)
    {
        iOneTop = iOne / SQRT_BASE;
        iOneBot = iOne % SQRT_BASE;
    }
    else
    {
        iOneTop = 0;
        iOneBot = iOne;
    }
    
    if (iTwo>SQRT_BASE)
    {
        iTwoTop = iTwo / SQRT_BASE;
        iTwoBot = iTwo % SQRT_BASE;
    }
    else
    {
        iTwoTop = 0;
        iTwoBot = iTwo;
    }
    
    iAnsBot = (iOneBot * iTwoBot) +
        (iOneBot * iTwoTop * SQRT_BASE) +
        (iOneTop * iTwoBot * SQRT_BASE);
    iAnsTop = (iOneTop * iTwoTop);

    /* Is there any carry from the bottom word to the top?            */
    if (iAnsBot>BASE)
    {
        iAnsTop = iAnsTop + (iAnsBot / BASE);
        iAnsBot = iAnsBot % BASE;
    }

    /* Set the caller's answer array:                                 */
    pAnswer[0] = iAnsTop;
    pAnswer[1] = iAnsBot;

    return iStatus;
}

/* Multiply an INT array by one INT                                   */
int cascMultiply  ( INT * iOne, INT * aTwo, INT ** aAnswer )
{
    /*
    * This routine allocates space for the answer, which the caller  
    * must, at some time, release.
    * The single INT is assumed to be positive, and the outgoing
    * sign is taken from the incoming sign on the very large number
    * (the second argument).
    */
    int iStatus = EXIT_SUCCESS;
    INT * pInt;
    INT ** ppInt = &pInt;
    INT w[2];
    int i = 0;
    int j = 0;
    int k = 0;
    int carry = 0;
    int iSign = 1;

    w[0] = 0;
    w[1] = 0;

    /*
    * The length of the product is at most one element longer than   
    * the incoming array. Allocate internal space for calculation:   
    */
    i = abs(*aTwo) + 1;
    k = aAllocate ( i, ppInt );
    /*
    * Remember the sign of the incoming array:
    */
    if (*aTwo<0)
        iSign = -1;
    else
        iSign = 1;
    /*
    * Set the decimal length of the new array to be the same as 
    * the incoming array:
    */
    *(pInt + 1) = *(aTwo + 1);

    /*
    * We multiply from the right to the left, bringing the carry over
    * and calculating the outgoing carry:                            
    */
    for (j=*aTwo; (j>1); j--)
    {
        /* Multiply just this pair of "digits": */
        pairMultiply ( *(aTwo + j), (*iOne), w );

        /* Add the carry into the least significant word: */
        *(pInt + j + 1) = w[1] + carry;

        /*
        * Compute the outgoing carry, doing an integer division, and
        * truncating to the integer part of the quotient:
        */
        carry = w[0] + (*(pInt + j + 1) / BASE);

        /* And the correct the word from which the carry came: */
        *(pInt + j + 1) %= BASE;
    }

    /*
    * The very last carry is used to set the topmost word in the
    * answer. This is the most significant word:
    */
    *(pInt + 2) = carry;

    /* Now we normalize the answer, and hand it back to the caller:   */
    i = aNormalize ( ppInt, aAnswer );

    /* Ensure that the outgoing sign is correct:                      */
    if (iSign<0)
        **aAnswer = - **aAnswer;

    return iStatus;
}

/* Multiply two arrays of INT                                         */
int aMultiply ( INT * aOne, INT * aTwo, INT ** aAnswer)
{
    /*
    * This routine allocates space for the answer, which the caller  
    * must, at some time, release                                    
    */
    int iStatus = EXIT_SUCCESS;
    int i = 0;
    int j = 0;
    int k = 0;
    int m = 0;
    INT * pInt;
    INT ** ppInt = &pInt;
    INT w[2];
    
    /*
    * The function of this routine is to multiply two very large     
    * numbers. Each of them comes in an array of INT. This routine   
    * works out the length required for the answer, allocates    
    * space for that answer, computes the product, and normalizes    
    * the answer, which is then passed back to the caller.           
    */
    
    /* Firstly check whether the incoming pointers are reasonable:    */
    if ((aOne==NULL) || (aTwo==NULL) || (aAnswer==NULL))
    {
        /*
        * There is a fault in the incoming pointers. Indicate this   
        * to the caller:                                             
        */
        iStatus = EXIT_FAILURE;
    }
    else
    {
        /*
        * The incoming pointers seem to be OK. Before we enter the main
        * (and long) body of this routine, we can check for two special
        * cases. These are if either multiplicand is zero or either is
        * one. In these cases we do not have to perform the full operation,
        * but we know the answer already.
        * So is the first operand zero or one?
        */
        if ((abs(aOne[0])==2) && (aOne[1]==0))
        {
            if (aOne[2]==0)
            {
                /*
                * The first multiplicand is zero. Hence the answer is
                * a normalized zero:
                */
                j = aAllocate ( 2, aAnswer );
                if (j!=0)
                    iStatus = EXIT_FAILURE;