ch24aok.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 integers, and not numbers or ranges containing a decimal
*  point.
*
*  Code by Ian D. K. Kelly
*
*/

/* Ch24AOK.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 numbers. 
* The very large numbers are assumed here to be integers - whole.    
* All very large numbers are assumed to be stored in arrays of int.  
* The first element of each of these arrays is a counter of the total
* number of further elements that exist in that array. The value of  
* the number being stored is from left to right, from 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".                          
*
*/

/*
*
* 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).
*
* 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).                
*
*/

#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;

#else

/*
*
* INT_MAX is not large enough. 
*
*/
#define BASE 10000
#define SQRT_BASE 100
typedef long int INT;

#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 ** aAnswer );
int aDivide       ( INT * aNumerator,   INT * aDenominator, INT ** aAnswer );
int aDivNormalize ( INT * pinNumerator, 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, becasue 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 );
    if (*aTwo<0)
        iSign = -1;
    else
        iSign = 1;

    /*
    * We multiply from the right to the left, bringing the carry over
    * and calculating the outgoing carry:                            
    */
    for (j=*aTwo; (j>0); 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 + 1) = 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. So now determine the  
        * size of the output array required:                         
        */
        i = abs(aOne[0]) + abs(aTwo[0]) + 1;

        /* Allocate enough space for the answer:                      */
        j = aAllocate(i, ppInt);

        /* Did this allocation work? if not, tell caller:             */
        if (j!=0)
            iStatus = EXIT_FAILURE;
        else
        {
            /*
            * The allocation appeared to work. Now compute the product 
            * of the two incoming, by computing the pairwise product
            * of each pair of incoming INT's, shifting up each time.
            *  
            * Note that multiplication proceeds from right to left:  
            * that is, we use the RIGHTmost digits first - just as we
            * do on pencil and paper.                                
            */
            for (k=abs(aTwo[0]); (k>0); k--)
            {
                for (j=abs(aOne[0]); (j>0); j--)
                {
                    /*
                    * Compute the pairwise product of One[j] by      
                    * Two[k], and place the result in answer[L-j]    
                    * where L = i + 1 - k                            
                    * Note that this may overflow, so we have to add 
                    * the high end of our answer into the next answer
                    * element "to the left" (with lower subscript):  
                    */
                    m = pairMultiply ( aOne[j], aTwo[k], w);
                    pInt[i + 2 - k - j] += w[1];
                    pInt[i + 1 - k - j] += w[0];
                }   /* End of inner "for/j" along first operand       */
            }       /* end of outer "for/k" along second operand      */

            /*
            * Now consider the sign of the output. Currently it is   
            * set to positive, but we need to apply the logic:       
            *              * | +   -                                 
            *             ___________                                
            *              + | +   -                                 
            *              - | -   +                                 
            */
            if (((aOne[0]>0)&&(aTwo[0]<0))
            ||  ((aOne[0]<0)&&(aTwo[0]>0)))
                pInt[0] = - pInt[0];

            /*
            * At this point we have in the local array the possibly  
            * un-normalized product of the input arrays. Normalize:  
            */
            m = aNormalize (ppInt, aAnswer);

            /* If there was a problem in normalization. Tell caller:  */
            if (m!=0)
                iStatus = EXIT_FAILURE;
        }   /* end of "if/else" allocation worked                     */
    }   /* end of "if/else" the incoming pointers were OK             */

    return iStatus;
}

/* Divide one INT by another INT, getting quotient and remainder      */
int oneDivide ( INT * iNumerator, INT * iDenominator, INT * iAnswerOD )
{