ch24abad.c

来自「稀疏矩阵、链表、图、队列、二叉树、多叉树、排序、遗传算法等的实现」· C语言 代码 · 共 1,132 行 · 第 1/3 页

/*
*
*  This code is part of Chapter 24 of "C Unleashed", and is to be     
*  read in conjunction with that chapter. It is NOT code of production
*  quality, and should not be so considered.
*
*  Code by Ian D. K. Kelly
*
*/

/* Ch24Abad.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 three of the basic four arithmetic   
* routines for extremely large numbers. This sample deals mostly with
* positive numbers (hence there is no mention here of subtraction).  
* The very large numbers are also assumed 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 in 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 though 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
* There is one other restriction on the size of BASE which you can   
* learn by careful inspection of the code. What is it?               
*
*
* 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).                
*
* As an exercise you might like to look at what is required to change
* this base to one hundred million (ten to the eight), whose square  
* root is ten thousand (ten to the four). You will need to be very   
* cautious about overflow during multiplication and division. There  
* are some hints about this in the implementation of "pairMultiply". 
*
*/

#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. We shall assume, without further testing
* (though we should test!) that "int" is 16 bit, and UINT_MAX is 
* at least 20000. It will be, in fact, at least 65535. Note that this
* means we have to change the definition of all our work variables   
* from being "int" to being "unsigned int". That is why we declare   
* all of them to be of type INT                                      
* By the way, O careful reader, there is another problem here: can   
* you find it? and can you fix it? The solution to this problen may
* be found in the file Ch24AOK.c on the CD.
*
*/
#define BASE 10000
#define SQRT_BASE 100
typedef long int INT;

#endif

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 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.                        
    */
    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;
    }
    
    /*
    * Consider the two-by-two multiplcation:                         
    *                 a   b                                          
    *                 c   d                                          
    *           --------------                                       
    *                da  db                                          
    *            ca  cb                                              
    *           --------------                                       
    *            ca, da+cb db                                        
    *           ==============                                       
    * where we have put spaces between the "sqrt_base" units, and a  
    * comma at "base". This explains the following inscrutable lines 
    */
    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 callers 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                                    
    */
    int iStatus = EXIT_SUCCESS;
    INT * pInt;
    INT ** ppInt = &pInt;
    INT w[2];
    int i = 0;
    int j = 0;
    int k = 0;
    int carry = 0;

    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 = (*aTwo) + 1;
    k = aAllocate ( i, ppInt );

    /*
    * 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 vary 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 );

    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 (m!=0)
            {
                /* There was a problem in normalization. Tell caller: */
                iStatus = EXIT_FAILURE;
            }   /* end of "if/then" there was a normalization problem */
        }   /* 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 )
{
    int iStatus = EXIT_SUCCESS;