cat mathset.f
/* mathset.f may 3, 1987 */
/* mathset includes divide(p,h),residue(p,h),sub(p,h),multiply(p,h,r) and
the fundamental functions which these call: copy,pad,add,halve and mod   */
/* version of May 2, 1987, developed while debugging minvert.c           */
divide(p,h) /* Divides p by h. Leaves result in p. */
char *p, *h  ;
{       char *j, *k, *l, *r ;
        j = q + 135 ;
        k = q + 150 ;
        l = q + 165 ;
        r = q + 180 ;
        copy(one,j) ;
        copy(xero,r) ;
        copy(h,k) ;
        copy(p,l) ;
        while(strcmp(k,l) <= 0)
        {       while(strcmp(k,l) <= 0)
                {       add(k,k) ;
                        add(j,j) ;
                }
                halve(k) ;
                halve(j) ;
                add(r,j) ;
                mod(l,k) ;
                copy(h,k) ;
                copy(one,j) ;
        }
        copy(r,p) ;
}

residue(p,h) /* Replaces p with the residue of p mod h */
char *p, *h  ;
{       char *j, *k, *l, *r ;
        j = q + 255 ;
        k = q + 270 ;
        l = q + 285 ;
        r = q + 300 ;
        copy(one,j) ;
        copy(xero,r) ;
        copy(h,k) ;
        copy(p,l) ;
        while(strcmp(k,l) <= 0)
        {       while(strcmp(k,l) <= 0)
                {       add(k,k) ;
                        add(j,j) ;
                }
                halve(k) ;
                halve(j) ;
                add(r,j) ;
                mod(l,k) ;
                copy(h,k) ;
                copy(one,j) ;
        }
        copy(l,p) ;
}












sub(p,h) /* Subtracts number at h from number at p, leaves result at p. */
char *p, *h ;
{       int i, borrow ;
        borrow = 0 ;
        for (i = 12 ; i > -1   ; --i)
        {       *(p + i) = *(p + i) - borrow ;
                if(*(p + i) < *(h + i))
                {       *(p + i) = 10 + *(p + i) - *(h + i) + '0' ;
                        borrow = 1 ;
                }
                else
                {       *(p + i) = *(p + i) - *(h + i) + '0' ;
                        borrow = 0 ;
                }
        }
return(borrow) ;
}

pad(p) /* Expands a numberstring shorter than 13 digits to a 13 digit length */
char *p ;
{       char *c;
        int i,j ;
        c = q + 120 ;
        for (i = 12, j = strlen(p) - 1 ; j > -1   ; --i, --j)
        *(c + i) = *(p + j) ;
        while( i > -1)  *(c + i--) = '0' ;
        copy(c,p) ;
}
copy(p,h) /* Copies the numberstring at p to location h. */
char *p, *h ;
{
        int i ;
        char c ;
        for (i = 0 ; i < 13   ; ++i)
        *(h + i) = *(p + i) ;
}
add(p,h) /* Adds numberstring at h to the one at p.  Result remains in p. */
char *p, *h ;
{       int i, carry ;
        char c ;
        carry = 0 ;
        for (i = 12 ; i > -1   ; --i)
        {       c = ((*(h + i) + *(p + i) - '0' - '0' + carry) % 10 )+ '0' ;
                carry = (*(h + i) + *(p + i) - '0' - '0' + carry) / 10  ;
                *(p + i) = c ;
        }
        return(carry) ;
}
halve(p) /* Halves the number at p. If it was odd, returns integer 1, else 0. */
char *p ;
 {      int i, carry ;
        char c ;
        carry = 0 ;
        for (i = 0 ; i < 13   ; ++i)
        {       c = ((*(p + i) - '0' + carry) / 2 ) + '0' ;
                carry = ((*(p + i) - '0' + carry ) % 2) * 10 ;
                *(p + i) = c ;
        }
        return(carry / 10 ) ;
}








mod(p,h) /* Converts number at p to that number modulo value of number at h. */
char *p, *h ;
{       int i, borrow ;
        while(strcmp(p,h) >= 0)
        {       borrow = 0 ;
                for (i = 12 ; i > -1   ; --i)
                {       *(p + i) = *(p + i) - borrow ;
                        if(*(p + i) < *(h + i))
                        {       *(p + i) = 10 + *(p + i) - *(h + i) + '0' ;
                                borrow = 1 ;
                        }
                        else
                        {       *(p + i) = *(p + i) - *(h + i) + '0' ;
                                borrow = 0 ;
                        }
                }
        }
return(borrow) ;
}





multiply(p,h,m) /* Multiplies p by h modulo m. Leaves result in p. */
char *p, *h, *m ;
{       char *j, *k, *l ;
        j = q + 225 ;
        k = q + 240 ;
        l = q + 255 ;
        j = "0000000000000" ;
        copy(xero,j) ;
        copy(p,k) ;
        copy(h,l) ;
        while(strcmp(xero,l) < 0)
        {       if(halve(l) == 1)
                {       add(j,k) ;
                        mod(j,m) ;
                }
                add(k,k) ;
                mod(k,m) ;
        }
copy(j,p) ;
}

$