ch24aok2.c

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

    /*
    * There was nothing special about the denominator. How about
    * the numerator? We only look for value zero:
    */
    if ((abs(*aNumerator)==2) && (*(aNumerator+1)==0) 
         && (*(aNumerator+2)==0))
    {
        /*
        * The numerator does have the value zero. Allocate a zero
        * value answer, and return to caller:
        */
        j = aAllocate ( 2, aAnswer );
        if (j!=0)
            iStatus = EXIT_FAILURE;
        return iStatus;
    }   /* end of "if/then" numerator has value zero */

    /*
    * If control falls through to here, we have the general, most usual,
    * case, and we will have to perform the bulk of the division
    * operation.
    */

    /*
    * Normalization for division means (a) ensuring that the numerator
    * is longer than the denominator, and (b) the leading "digit" 
    * of the denominator is greater than floor(BASE/2). We can get   
    * this by multiplying both numerator and denominator by "D",     
    * where "D = floor((b-1)/d" and "d" is the leading "digit" of the
    * denominator. "floor" means "round to the next integer below".
    */

    /* Find the length required for the denominator copy, and allocate*/
    iDmax = abs(aDenominator[0]);
    k = aAllocate ( iDmax, ppInter );

    /*
    * Copy from argument to local copy. Note that this is a simple copy
    * which also copies the count of the number of places after the point
    */
    for (i=1;(i<=iDmax);i++)
        *(pInter+i) = *(aDenominator+i);

    /* And normalize the local copy of the denominator:               */
    k = aNormalize ( ppInter, ppInter );

    /*
    * Get the length of the normalized copy, and compute the scaling 
    * factor by which to multiply both the numerator and denominator: 
    */
    iDmax = *pInter;
    littleD = *(pInter + 1);
    bigD = (BASE - 1) / littleD;

    /*
    * Get a scaled (multiplied-up) local copy of the denominator, and
    * release the previous work area.
    */
    cascMultiply ( &bigD, pInter, pplDenominator );
    free ( pInter );

    /*
    * Find out how much space to allocate for a local copy of the 
    * numerator, and get the space, multiplying up at the same time:    
    */
    pplNumerator = &plNumerator;
    plNumerator = NULL;
    cascMultiply ( &bigD, aNumerator, ppInter );

    /*
    * Perform "Division Normalize", which inserts a leading zero, and
    * release the intermediate work area:                            
    */
    k = aDivNormalize ( pInter, places, pplNumerator );
    free ( pInter );

    /*
    * Get the lengths of the normalized numerator, and compute the   
    * length of the quotient:                                        
    * Knuth D1:                                                      
    */
    iNmax = *plNumerator;
    n = iDmax;
    m = iNmax - iDmax;

    /*
    * The length of the quotient was one greater than the difference 
    * in length of the numerator and denominator. Remember that this 
    * routine requires the numerator to be at least as long as the   
    * denominator. The difference in length (and not the "difference 
    * plus one") is required later in the calculation. We have already
    * "lengthened" the numerator, so we need use only the real 
    * difference in length now:                                      
    */

    /* Allocate space for the quotient pre-normalization:             */
    pplQuotient = &plQuotient;
    k = aAllocate ( m, pplQuotient );

    /* Allocate space for holding the intermediate subtraction        */
    pplSubtrahend = &plSubtrahend;
    aAllocate ( n, pplSubtrahend );

    /* Allocate space for holding numerator portion for subtraction:  */
    pplMinuend = &plMinuend;
    aAllocate ( n + 1, pplMinuend );

    /*
    * Division takes place from left to right: we start with the most
    * significant digits - those on the left:                        
    * Knuth D2:                                                      
    */
    for (j=2;(j<=m);j++)
    {
        /* Knuth D3:                                                  */
        wN[0] = *(plNumerator + j); 
        wN[1] = (j==(m+1)) ? 0 : *(plNumerator + j + 1); 
        k = pairDivide ( wN, plDenominator + 1, aW );

        /*
        * Now test for the cases where our putative quotient is too  
        * high. Note that it will sometimes be one too high, and very
        * rarely two too high.                                       
        * If quotient is equal to base, it is too high:              
        */
D3A:    if ((aW[0]==0) && (aW[1]==BASE))
            goto D3B;
        /* If quotient is equal to zero, it cannot be too high:       */
        if ((aW[0]==0) && (aW[1]==0))
            goto D3C;
        /* If there are no more places in the denominator, then there */
        /* is nothing else against which to test, so it's not too high*/
        if (*plDenominator<2)
            goto D3C;

        /* Now we have to compute the inequality of Knuth step D3     */
        pairMultiply (*(plDenominator+2),aW[0],aX);
        pairMultiply (*(plDenominator+2),aW[1],aY);
        aY[0] += aY[1] / BASE; 
        aY[1] %= BASE;
        aX[1] += aY[0] / BASE; 
        aY[0] %= BASE;
        aZ[0] = aW[2];        
        aZ[1] = *(plNumerator + j - m + 3); 
        carry = aZ[1] / BASE; 
        aZ[1] %= BASE;

        if (aX[1]>carry)
            goto D3B;
        if (aY[0]>aZ[0])
            goto D3B;
        if (aY[0]<aZ[0]) 
            goto D3C;
        if (aY[1]>aZ[1]) 
            goto D3B;
        goto D3C;

D3B:    {
            /*
            * The putative quotient is too high. Reduce it, and go   
            * back round to see whether it is still too high:         
            */
            aW[1] -= 1;
            aW[2] += *(plDenominator + 1);
            if (aW[2]<BASE)
                goto D3A;
        }
D3C:    /*
        * This quotient digit is correct.                            
        * Knuth D4:                                                  
        * We now have to subtract from the local copy of the numerator
        * the denominator times the last quotient "digit":        
        */
        if (aW[1]!=0)
        {
            for (i=2; (i<=(n+1)); i++)
            {
                *(plMinuend + i) = *(plNumerator + m - j + i);
                if (i<(n + 1))
                    *(plSubtrahend + i) = *(plDenominator + i);
            }

            /* cascade multiply denominator by quotient digit */
            cascMultiply ( &aW[1], plDenominator, ppInter );

            /* subtract:                                      */
            if ( *pInter > *plMinuend )
                borrow = 1;
            else
            {
                carry = 0;
                k = *pInter;
                l = (*plMinuend) - k;
                for (i=k; (i>=0); i--)
                {
                    if ((i + l)>0)
                    {
                        *(plMinuend + i + l) -= carry; 
                        *(plMinuend + i + l) -= (i>0) ? *(pInter+i) : 0 ;
                        carry = 0;
                        if (*(plMinuend + i + l)<0)
                        {
                            *(plMinuend + i + l) += BASE;
                            carry = 1;
                        }
                        else
                            carry = 0;
                    }   /* end of "if/then" safe to make replacement  */
                }   /* end of "for/i" computing digits to subtract    */

                /*
                * Test the sign of the answer, and if it is negative then
                * set the "borrow" flag. If it is positive or zero then
                * replace a portion of the local copy of the numerator:
                */
                if (*(plMinuend+1)<0)
                    borrow = 1;
                else
                {
                    borrow = 0;
                    /* Replace a portion of our local numerator copy  */
                    for (i=n+1; (i>0); i--)
                    {
                        if ((i + l) <= iNmax)
                           *(plNumerator + i + l) = *(plMinuend + i);
                    }   /* end of "for/i" replace portion of numerator    */
                }   /* end of "if/else" numerator portion was replaced    */
            }   /* end of "if/else" lengths allowed numerator subtraction */

            /* Knuth D5:                                                  */
            *(plQuotient + j) = aW[1];

            /* If the remainder was negative then add back:               */
            if (borrow != 0)
            {
                /* Knuth D6:                                              */
                /*
                * Note that this is a very low probability step. It happens 
                * only in the order of 2/BASE. We have to add back, but
                * ignoring the top end carry:                        
                */
                *(plQuotient + j) = *(plQuotient + j) - 1;
                carry = 0;
                for (i=n+1; (i>1); i--)
                {
                    k = (*plNumerator) - i;
                    *(plNumerator+j+k) += (i==1) ? 0 : 
                                          (*(plDenominator + k) + carry);
                    if (*(plNumerator + j + k) > BASE)
                    {
                        *(plNumerator + j + k) -= BASE;
                        carry = 1;
                    }
                    else
                        carry =0;
                }   /* end of "for/i" adding back in                  */
            }       /* end of "if/then" borrow was non zero           */
        }           /* end of "if/then" quotient digit was non zero   */

        /* Knuth D7:                                                  */

        /* Loop back to get the next quotient "digit"                 */

    }   /* end of "for/j" major loop along digits of numerator        */

    /* Knuth D8:                                                      */
    /* Normalize the result:                                          */
    k = aNormalize ( pplQuotient, aAnswer );

    /* Calculate and apply the sign of the answer:
    *
    *  This is according to the logic: 
    *        /  |  +   - 
    *        _____________
    *        +  |  +   -
    *        -  |  -   +
    *
    */
    if (((*aNumerator<0) && (*aDenominator>0))
    ||  ((*aNumerator>0) && (*aDenominator<0)))
        **aAnswer = - **aAnswer;

    /*
    *  Since this was a non-integral division, there is no real meaning
    *  here to the concept "remainder". The residual error can, however,
    *  be estimated by using routine "cascDivide" to divide plNumerator
    *  by littleD.
    */

    /* Free the allocated intermediate work areas:                    */
    free ( pInter );
    free ( plSubtrahend );
    free ( plMinuend );
    free ( plQuotient );
    free ( plNumerator );
    free ( plDenominator );

    return iStatus;
}


/* Allocate an INT array for the aDivide internal operation which is  */
/* one longer than the incoming array:                                */
int aDivNormalize ( INT * pinNumerator, int places, INT ** ppoutNumerator )
{
    int iStatus = EXIT_SUCCESS;
    int iNmax = 0;
    int i = 0;
    int j = 0;
    int k = 0;
    int iAllocLength = 0;
    INT * pAnswer = NULL;
    INT ** ppAnswer = &pAnswer;

    /* What is the length required for the new allocation to make?    */
    iNmax = abs(*pinNumerator);
    k = iNmax + 1;
    j = ((places + DIGITS_IN_BASE - 1)/DIGITS_IN_BASE)+1;
    iAllocLength = max(j,k);

    /* make the allocation:                                           */
    i = aAllocate ( iAllocLength, ppAnswer );

    /*
    * Indicate the new number of places after the point for the new
    * array. This is the same as the number for the incoming array, 
    * plus any new words that have been added:
    */
    *(pAnswer + 1) = *(pinNumerator + 1) + (j>k) ? (j-k) : 0;

    /* Now copy the original into the new area.                       */
    /* Set up the initial zero:                                       */
    *(pAnswer + 2) = 0;

    /* Now copy the rest of the input array to the output:            */
    if (iAllocLength>=3)
    {
        for (k=3; (k <= iAllocLength); k++)
            *(pAnswer + k) = *(pinNumerator + k - 1);
    }
    /*
    * Note that the remaining words, if any, will already have been
    * set to zero during the allocation process
    */

    /* Point the answer to the newly allocated:                      */
    *ppoutNumerator = pAnswer;

    return iStatus;
}

/* Normalize an INT array                                             */
int aNormalize ( INT ** aUnNormal, INT ** aNormal )
{
    /*
    * Although this routine may allocate space for its answer, if it 
    * does so it frees the space used by its first argument. Hence the   
    * caller does NOT have to consider whether to free the returned  
    * space from this routine.                                       
    */

    int iStatus = EXIT_SUCCESS;
    int i = 0;
    int j = 0;
    int k = 0;
    int n = 0;
    int m = 0;
    int mm = 0;
    int b = 0;
    int iSign = 1;
    INT * plNormal = NULL;
    INT ** pplNormal = &plNormal;
    
    /*
    * This function takes two arguments. The first is a pointer to