home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 10 Tools / 10-Tools.zip / OL.LZH / PROCS.LZH / LARGINT.ICN < prev    next >
Text File  |  1991-07-13  |  4KB  |  174 lines
  1. ############################################################################
  2. #
  3. #    Name:    largint.icn
  4. #
  5. #    Title:    Large integer arithmetic
  6. #
  7. #    Author:    Paul Abrahams and Ralph E. Griswold
  8. #
  9. #    Date:    May 11, 1989
  10. #
  11. ############################################################################
  12. #
  13. #     These procedures perform addition, multiplication, and exponentiation
  14. #  On integers given as strings of numerals:
  15. #
  16. #        add(i,j)      sum of i and j
  17. #
  18. #        mpy(i,j)      product of i and j
  19. #
  20. #        raise(i,j)    i to the power j
  21. #
  22. #  Note:
  23. #
  24. #     The techniques used by add and mpy are different from those used by
  25. #  raise.  These procedures are combined here for organizational reasons.
  26. #  The procedures add and mpy are adapted from the Icon language book.
  27. #  The procedure raise was written by Paul Abrahams.
  28. #
  29. ############################################################################
  30.  
  31. record largint(coeff,nextl)
  32.  
  33. global base, segsize
  34.  
  35. # Add i and j
  36. #
  37. procedure add(i,j)
  38.  
  39.    return lstring(addl(large(i),large(j)))
  40.  
  41. end
  42.  
  43. # Multiply i and j
  44. #
  45. procedure mpy(i,j)
  46.  
  47.    return lstring(mpyl(large(i),large(j)))
  48.  
  49. end
  50.  
  51. # Raise i to power j
  52. #
  53. procedure raise(i,j)
  54.  
  55.      return rstring(ipower(i,binrep(j)))
  56.  
  57. end
  58.  
  59. procedure addl(g1,g2,carry)
  60.    local sum
  61.    /carry := largint(0)    # default carry
  62.    if /g1 & /g2 then return if carry.coeff ~= 0 then carry
  63.    else &null
  64.    if /g1 then return addl(carry,g2)
  65.    if /g2 then return addl(g1,carry)
  66.    sum := g1.coeff + g2.coeff + carry.coeff
  67.    carry := largint(sum / base)
  68.    return largint(sum % base,addl(g1.nextl,g2.nextl,carry))
  69. end
  70.  
  71. procedure large(s)
  72.    initial {
  73.       base := 10000
  74.       segsize := *base - 1
  75.       }
  76.  
  77.    if *s <= segsize then return largint(integer(s))
  78.    else return largint(right(s,segsize),
  79.       large(left(s,*s - segsize)))
  80. end
  81.  
  82. procedure lstring(g)
  83.    local s
  84.  
  85.    if /g.nextl then s := g.coeff
  86.    else s := lstring(g.nextl) || right(g.coeff,segsize,"0")
  87.    s ?:= (tab(upto(~'0') | -1) & tab(0))
  88.    return s
  89. end
  90.  
  91. procedure mpyl(g1,g2)
  92.    local prod
  93.    if /(g1 | g2) then return &null    # zero product
  94.    prod := g1.coeff * g2.coeff
  95.    return largint(prod % base,
  96.       addl(mpyl(largint(g1.coeff),g2.nextl),mpyl(g1.nextl,g2),
  97.       largint(prod / base)))
  98. end
  99.  
  100. # Compute the binary representation of n (as a string)
  101. #
  102. procedure binrep(n)
  103.     local retval
  104.     retval := ""
  105.     while n > 0 do {
  106.         retval := n % 2 || retval
  107.         n /:= 2
  108.         }
  109.     return retval
  110. end
  111.  
  112. # Compute a to the ipower bbits, where bbits is a bit string.
  113. # The result is a list of coefficients for the polynomial a(i)*k^i,
  114. # least significant values first, with k=10000 and zero trailing coefficient
  115. # deleted.
  116. #
  117. procedure ipower(a, bbits)
  118.     local b, m1, retval
  119.     m1 := (if a >= 10000 then [a % 10000, a / 10000] else [a])
  120.     retval := [1]
  121.     every b := !bbits do {
  122.         (retval := product(retval, retval)) | fail
  123.         if b == "1" then
  124.             (retval := product(retval, m1)) | fail
  125.         }
  126.     return retval
  127. end
  128.  
  129. # Compute a*b as a polynomial in the same form as for ipower.
  130. # a and b are also polynomials in this form.
  131. #
  132. procedure product(a,b)
  133.     local i, j, k, retval, x
  134.     if *a + *b > 5001 then
  135.         fail
  136.     retval := list(*a + *b, 0)
  137.     every i := 1 to *a do
  138.         every j := 1 to *b do {
  139.             k := i + j - 1
  140.             retval[k] +:= a[i] * b[j]
  141.             while (x := retval[k]) >= 10000 do {
  142.                 retval[k + 1] +:= x / 10000
  143.                 retval[k] %:= 10000
  144.                 k +:= 1
  145.             }   }
  146.     every i := *retval to 1 by -1 do
  147.         if retval[i] > 0 then
  148.             return retval[1+:i]
  149.     return retval[1+:i]
  150. end
  151.  
  152. procedure rstring(n)
  153.     local ds, i, j, k, result
  154.  
  155.     ds := ""
  156.     every k := *n to 1 by -1 do
  157.         ds ||:= right(n[k], 4, "0")
  158.     ds ?:= (tab(many("0")), tab(0))
  159.     ds := repl("0", 4 - (*ds - 1) % 5) || ds
  160.  
  161.     result := ""
  162.     every i := 1 to *ds by 50 do {
  163.         k := *ds > i + 45 | *ds
  164.         every j := i to k by 5 do {
  165.        ds
  166.            result ||:= ds[j+:5]
  167.        }
  168.         }
  169.    result ? {
  170.       tab(many('0'))
  171.       return tab(0)
  172.       }
  173. end
  174.