modules/gmp.lsp

   1 ;; @module gmp.lsp
   2 ;; @description GNU MP Bignum Library interface
   3 ;; @version 1.5 - new library detection rotuine
   4 ;; @author Lutz Mueller, 2007
   5 ;;
   6 ;; <h3>The GNU MP Bignum Library</h3>
   7 ;; This modules interfaces to libgmp which can be obtained from
   8 ;; @link http://www.swox.com/gmp/ www.swox.com/gmp/ .
   9 ;;
  10 ;; Source code for libgmp.so (UNIX) or lbgmp.dll (Win32) or libgmp.dylib (Mac OS X)
  11 ;; is available at this site.
  12 ;;
  13 ;; This interface module presets the maximum precision to 1024.
  14 ;; The precision can be changed to any other value by changing the
  15 ;; definition of 'MAX_PRECISION' in the source of this module.
  16 ;;
  17 ;; All arguments to the GMP functions in this module must be given as strings,
  18 ;; an error will be thrown otherwise. When starting with a 0 the number is
  19 ;; assumed to be in octal format when starting with 0x in hexadecimal format.
  20 ;;
  21 ;; This file only imports a few functions from the many available in GNU GMP.
  22 ;; See the GMP manual for more functions.
  23 ;;
  24 ;; Note, that since version 8.9.7 newLISP does all integer arithmetik with
  25 ;; 64 bits giving up to 19 digits of precision. For precisions less or equal
  26 ;; 19 digits newLISP's built-in 64-bit integer arithmetik is much faster.
  27 ;;
  28 ;; <h3>Usage</h3>
  29 ;; At the beginning of the programfile include a 'load' statement for the module:
  30 ;; <pre>
  31 ;; (load "/usr/share/newlisp/modules/gmp.lsp")
  32 ;; </pre>
  33 ;;
  34 ;; @example
  35 ;; (GMP:+ "123456789012345678901234567890" "123456789012345678901234567890")
  36 ;; => "246913578024691357802469135780"
  37
  38 ;; <h3>Adding more functions to the library</h3>
  39 ;; When adding functions be aware that inside the GMP context
  40 ;;  +,-,*,/,=,&lt;,&gt;,&lt;=,&gt;=  are overwritten for multiple precision and the
  41 ;; original operators  would have would have to be prefixed with MAIN when used,
  42 ;; inside the 'gmp.lsp' module.
  43
  44 ;; <br/><br/>
  45 ;; <center><h2>Integer arithmetik</h2></center>
  46 ;; @syntax (GMP:+ <arg1> <arg2>)
  47 ;; add two integers in <arg1> and <arg2>
  48 ;; @syntax (GMP:- <arg1> <arg2>)
  49 ;; subtract <arg2> from <arg1>
  50 ;; @syntax (GMP:* <arg1> <arg2>)
  51 ;; multiply <arg1> by <arg2>
  52 ;; @syntax (GMP:/ <arg1> <arg2>)
  53 ;; divide <arg1> by <arg2>, round towards '0' (zero)
  54 ;; @syntax (GMP:% <arg1> <arg2>)
  55 ;; calc rest of division <arg1>/<arg2>
  56 ;; @syntax (GMP:** <arg1> <arg2>)
  57 ;; calc power(<arg1>, <arg2>)
  58 ;; @syntax (GMP:= <arg1> <arg2>)
  59 ;; test for <arg1> equal <arg2>
  60 ;; @syntax (GMP:< <arg1> <arg2>)
  61 ;; test for <arg1> smaller <arg2>
  62 ;; @syntax (GMP:> <arg1> <arg2>)
  63 ;; test for <arg1> bigger <arg2>
  64 ;; @syntax (GMP:<= <arg1> <arg2>)]
  65 ;; test for <arg1> smaller or equal <arg2>
  66 ;; @syntax (GMP:>= <arg1> <arg2>)
  67 ;; test for <arg1> bigger or equal <arg2>
  68 ;;
  69 ;; <center><h2>Bit operations</h2></center>
  70 ;; @syntax (GMP:& <arg1> <arg2>)
  71 ;; bitwise <and> of <arg1>, <arg2>
  72 ;; @syntax (GMP:| <arg1> <arg2>)
  73 ;; bitwise inclusive <or> of <arg1>, <arg2>
  74 ;; @syntax (GMP:^ <arg1> <arg2>)
  75 ;; bitwise exclusive <or> of <arg1>, <arg2>
  76 ;; @syntax (GMP:~ <arg>)
  77 ;; bitwise complement of <arg>
  78 ;;
  79 ;; <center><h2>Number theory</h2></center>
  80 ;; @syntax (GMP:prime? <arg>)
  81 ;; check if <arg> is prime
  82 ;; @syntax (GMP:next-prime <arg>)
  83 ;; calc closes prime greater than <arg>
  84 ;; @syntax (GMP:factor <arg>)
  85 ;; calc a list of prime factors for <arg>
  86 ;; @syntax (GMP:gcd <arg1> <arg2>)
  87 ;; greatest common divisor of <arg1> and <arg2>
  88 ;; @syntax (GMP:bin <arg1> <arg2>)
  89 ;; calc binomial (<arg1> <arg2>)
  90 ;; @syntax (GMP:fac <arg>)
  91 ;; arg! factorial(<arg>)
  92 ;; @syntax (GMP:fib <arg>)
  93 ;; fibonacci(arg)
  94 ;;
  95 ;; <center><h2>Random numbers</h2></center>
  96 ;; @syntax (GMP:seed <arg>)
  97 ;; seed the random generator
  98 ;; @syntax (GMP:rand <arg>)
  99 ;; generate random numbers between 0 and arg - 1
 100 ;;
 101
 102 (constant '+ add '- sub '* mul '/ div)
 103
 104 (context 'GMP)
 105
 106 ; maximum digits, can be set to any value higher if required
 107 ; when choosing a different number the functions GMP:fac and GMP:fib
 108 ; have to be re-calibrated for their maximum value not to overflow
 109 ; the result space
 110 ;
 111 (define MAX_PRECISION 1024)
 112
 113 (define gmp-type-error "String expected in GMP module.")
 114 (define gmp-size-error "Argument too big in GMP module")
 115 (define gmp-divzero-error "Division by zero in GMP module")
 116
 117 (set 'files '(
 118     "/opt/local/lib/libgmp.3.4.2.dylib" ;Mac OSX
 119     "/opt/local/lib/libgmp.3.3.3.dylib" ;Mac OSX
 120     "/opt/local/lib/libgmp.3.dylib" ;Mac OSX
 121     "/opt/local/lib/libgmp.dylib" ;Mac OSX
 122     "/usr/local/lib/libgmp.so" ; Linux, BSDs
 123     "libgmp-3.dll" ; Win32 in path
 124 ))
 125
 126 (set 'library (files (or
 127                                 (find true (map file? files))
 128                                 (begin (println "cannot find GMP library") (exit)))))
 129
 130 ; integer arithmetik
 131 (import library "__gmpz_init")
 132 (import library "__gmpz_add")
 133 (import library "__gmpz_add_ui")
 134 (import library "__gmpz_sub")
 135 (import library "__gmpz_mul")
 136 (import library "__gmpz_tdiv_q")
 137 (import library "__gmpz_tdiv_r")
 138 (import library "__gmpz_cmp")
 139 (import library "__gmpz_cmp_si")
 140 (import library "__gmpz_set_si")
 141 (import library "__gmpz_divisible_p")
 142 (import library "__gmpz_pow_ui")
 143
 144 ; bit operators
 145 (import library "__gmpz_and")
 146 (import library "__gmpz_ior")
 147 (import library "__gmpz_xor")
 148 (import library "__gmpz_com")
 149
 150 ; number theory
 151 (import library "__gmpz_probab_prime_p")
 152 (import library "__gmpz_nextprime")
 153 (import library "__gmpz_gcd")
 154 (import library "__gmpz_bin_ui")
 155 (import library "__gmpz_fac_ui")
 156 (import library "__gmpz_fib_ui")
 157
 158 ; random numbers
 159 (import library "__gmpz_urandomm")
 160 (import library "__gmp_randseed")
 161
 162 ; conversion functions
 163 (import library "__gmpz_set_str")
 164 (import library "__gmpz_get_str")
 165
 166 ; auxiliary functions
 167 (import library "__gmpz_sizeinbase")
 168 (import library "__gmp_randinit_default")
 169 (import library "__gmp_randseed_ui")
 170
 171 ; reserve handles
 172 (define op1 (dup "\000" 12))
 173 (define op2 (dup "\000" 12))
 174 (define rop (dup "\000" 12))
 175
 176 (define randstate (dup "\000" 60))
 177
 178 ; init handles
 179 (__gmpz_init op1)
 180 (__gmpz_init op2)
 181 (__gmpz_init rop)
 182
 183 ; handles to speed up factor
 184 (define mp-n (dup "\000" 12))
 185 (define mp-d (dup "\000" 12))
 186 (define mp-k (dup "\000" 12))
 187 (__gmpz_init mp-n)
 188 (__gmpz_init mp-d)
 189 (__gmpz_init mp-k)
 190
 191
 192 ; init randstate
 193 (__gmp_randinit_default randstate)
 194 (__gmp_randseed_ui randstate (time-of-day))
 195
 196 ; init / reserve result string
 197 (define rops (dup "\000" (MAIN:+ 2 MAX_PRECISION)))
 198
 199 ; add two integers
 200 ;
 201 (define (GMP:+ p1 p2)
 202   (if (or (not (string? p1)) (not (string? p2))) (throw-error gmp-type-error))
 203   (__gmpz_set_str op1 p1 0)
 204   (__gmpz_set_str op2 p2 0)
 205   (__gmpz_add rop op1 op2)
 206   (if (MAIN:> (__gmpz_sizeinbase rop 10) MAX_PRECISION)
 207         (throw-error gmp-size-error))
 208   (get-string (__gmpz_get_str rops 10 rop))
 209 )
 210
 211 ; subtract two integers
 212 ;
 213 (define (GMP:- p1 p2)
 214   (if (or (not (string? p1)) (not (string? p2))) (throw-error gmp-type-error))
 215   (__gmpz_set_str op1 p1 0)
 216   (__gmpz_set_str op2 p2 0)
 217   (__gmpz_sub rop op1 op2)
 218   (if (MAIN:> (__gmpz_sizeinbase rop 10) MAX_PRECISION)
 219         (throw-error gmp-size-error))
 220   (get-string (__gmpz_get_str rops 10 rop))
 221 )
 222
 223 ; multiply two integers
 224 ;
 225 (define (GMP:* p1 p2)
 226   (if (or (not (string? p1)) (not (string? p2))) (throw-error gmp-type-error))
 227   (__gmpz_set_str op1 p1 0)
 228   (__gmpz_set_str op2 p2 0)
 229   (__gmpz_mul rop op1 op2)
 230   (if (MAIN:> (__gmpz_sizeinbase rop 10) MAX_PRECISION)
 231         (throw-error gmp-size-error))
 232   (get-string (__gmpz_get_str rops 10 rop))
 233 )
 234
 235 ; divide two integers
 236 ; return floor value of result
 237 (define (GMP:/ p1 p2)
 238   (if (or (not (string? p1)) (not (string? p2))) (throw-error gmp-type-error))
 239   (__gmpz_set_str op1 p1 0)
 240   (__gmpz_set_str op2 p2 0)
 241   (if (MAIN:= (__gmpz_cmp_si op2 0) 0) (throw-error gmp-divzero-error))
 242   (__gmpz_tdiv_q rop op1 op2)
 243   (get-string (__gmpz_get_str rops 10 rop))
 244 )
 245
 246 ; modulo two integers
 247 ; return rest value of division result
 248 (define (GMP:% p1 p2)
 249   (if (or (not (string? p1)) (not (string? p2))) (throw-error gmp-type-error))
 250   (__gmpz_set_str op1 p1 0)
 251   (__gmpz_set_str op2 p2 0)
 252   (if (MAIN:= (__gmpz_cmp_si op2 0) 0) (throw-error gmp-divzero-error))
 253   (__gmpz_tdiv_r rop op1 op2)
 254   (get-string (__gmpz_get_str rops 10 rop))
 255 )
 256
 257 ; exponentiation function
 258 ; return power(p1 p2)
 259 ;
 260 (define (GMP:** p1 p2 , pexp)
 261   (if (or (not (string? p1)) (not (string? p2))) (throw-error gmp-type-error))
 262   (__gmpz_set_str op1 p1 0)
 263   (set 'pexp (int p2 0))
 264   (__gmpz_pow_ui rop op1 pexp)
 265   (if (MAIN:> (__gmpz_sizeinbase rop 10) MAX_PRECISION)
 266         (throw-error gmp-size-error))
 267   (get-string (__gmpz_get_str rops 10 rop))
 268 )
 269
 270 ; test of two integers are equal, return true if equal
 271 ; otherwise nil
 272 ;
 273 (define (GMP:= p1 p2)
 274   (if (or (not (string? p1)) (not (string? p2))) (throw-error gmp-type-error))
 275   (__gmpz_set_str op1 p1 0)
 276   (__gmpz_set_str op2 p2 0)
 277   (MAIN:= (__gmpz_cmp op1 op2) 0)
 278 )
 279
 280 ; test is p1 is smaller than p2
 281 ; return true or nil
 282 ;
 283 (define (GMP:< p1 p2)
 284   (if (or (not (string? p1)) (not (string? p2))) (throw-error gmp-type-error))
 285   (__gmpz_set_str op1 p1 0)
 286   (__gmpz_set_str op2 p2 0)
 287   (MAIN:< (__gmpz_cmp op1 op2) 0)
 288 )
 289
 290 ; test is p1 is smaller than p2
 291 ; return true or nil
 292 ;
 293 (define (GMP:> p1 p2)
 294   (if (or (not (string? p1)) (not (string? p2))) (throw-error gmp-type-error))
 295   (__gmpz_set_str op1 p1 0)
 296   (__gmpz_set_str op2 p2 0)
 297   (MAIN:> (__gmpz_cmp op1 op2) 0)
 298 )
 299
 300 ; test is p1 is smaller or eaual than p2
 301 ; return true or nil
 302 ;
 303 (define (GMP:<= p1 p2)
 304   (if (or (not (string? p1)) (not (string? p2))) (throw-error gmp-type-error))
 305   (__gmpz_set_str op1 p1 0)
 306   (__gmpz_set_str op2 p2 0)
 307   (MAIN:<= (__gmpz_cmp op1 op2) 0)
 308 )
 309
 310 ; test is p1 is bigger or equal than p2
 311 ; return true or nil
 312 ;
 313 (define (GMP:>= p1 p2)
 314   (if (or (not (string? p1)) (not (string? p2))) (throw-error gmp-type-error))
 315   (__gmpz_set_str op1 p1 0)
 316   (__gmpz_set_str op2 p2 0)
 317   (MAIN:>= (__gmpz_cmp op1 op2) 0)
 318 )
 319
 320 ; bitwise and two integers
 321 ;
 322 (define (GMP:& p1 p2)
 323   (if (or (not (string? p1)) (not (string? p2))) (throw-error gmp-type-error))
 324   (__gmpz_set_str op1 p1 0)
 325   (__gmpz_set_str op2 p2 0)
 326   (__gmpz_and rop op1 op2)
 327   (get-string (__gmpz_get_str rops 10 rop))
 328 )
 329
 330 ; bitwise inclusive or two integers
 331 ;
 332 (define (GMP:| p1 p2)
 333   (if (or (not (string? p1)) (not (string? p2))) (throw-error gmp-type-error))
 334   (__gmpz_set_str op1 p1 0)
 335   (__gmpz_set_str op2 p2 0)
 336   (__gmpz_ior rop op1 op2)
 337   (get-string (__gmpz_get_str rops 10 rop))
 338 )
 339
 340 ; bitwise exclusive or two integers
 341 ;
 342 (define (GMP:^ p1 p2)
 343   (if (or (not (string? p1)) (not (string? p2))) (throw-error gmp-type-error))
 344   (__gmpz_set_str op1 p1 0)
 345   (__gmpz_set_str op2 p2 0)
 346   (__gmpz_xor rop op1 op2)
 347   (get-string (__gmpz_get_str rops 10 rop))
 348 )
 349
 350 ; bitwise complement of one integer
 351 ;
 352 (define (GMP:~ p1)
 353   (if (not (string? p1)) (throw-error gmp-type-error))
 354   (__gmpz_set_str op1 p1 0)
 355   (__gmpz_com rop op1)
 356   (get-string (__gmpz_get_str rops 10 rop))
 357 )
 358
 359 ; check if a prime
 360 ; returns true if a prime
 361 ; returns nil if not a primes or probably not a prime
 362 ;
 363 (define (GMP:prime? p1 , r)
 364   (if (not (string? p1)) (throw-error gmp-type-error))
 365   (__gmpz_set_str op1 p1 0)
 366   (set 'r (__gmpz_probab_prime_p op1 10))
 367   (if
 368     (MAIN:= r 0) nil
 369     (MAIN:= r 1) (MAIN:= (next-prime (- p1 "1")) p1)
 370     (MAIN:= r 2) true)
 371 )
 372
 373 ; get the next prime higher then arg
 374 ;
 375 (define (GMP:next-prime p1)
 376   (if (not (string? p1)) (throw-error gmp-type-error))
 377   (__gmpz_set_str op1 p1 0)
 378   (__gmpz_nextprime rop op1)
 379   (get-string (__gmpz_get_str rops 10 rop))
 380 )
 381
 382 ; factorize n
 383 ;
 384 (define (GMP:factor n)
 385   (if (not (string? n)) (throw-error gmp-type-error))
 386   (set 'factors nil)
 387   (set 'prevfact nil)
 388   (__gmpz_set_str mp-n n 0)
 389   (if (MAIN:> (__gmpz_cmp_si mp-n 2) 0)
 390     (begin
 391       (__gmpz_set_si mp-d 2)
 392       (__gmpz_set_si mp-k 0)
 393
 394       (while (MAIN:!= 0 (__gmpz_divisible_p mp-n mp-d))
 395         (__gmpz_tdiv_q mp-n mp-n mp-d)
 396         (__gmpz_add_ui mp-k mp-k 1)
 397       )
 398
 399       (if (MAIN:> (__gmpz_cmp_si mp-k 0) 0)
 400         (push-factor mp-d mp-k))
 401
 402       (__gmpz_set_si mp-d 3)
 403       (__gmpz_mul op1 mp-d mp-d)
 404       (while (MAIN:<= (__gmpz_cmp op1  mp-n) 0)
 405         (__gmpz_set_si mp-k 0)
 406         (while (MAIN:!= 0 (__gmpz_divisible_p mp-n mp-d))
 407           (__gmpz_tdiv_q mp-n mp-n mp-d)
 408           (__gmpz_add_ui mp-k mp-k 1) )
 409         (if (MAIN:> (__gmpz_cmp_si mp-k 0) 0) (push-factor mp-d mp-k))
 410         (__gmpz_add_ui mp-d mp-d 2)
 411         (__gmpz_mul op1 mp-d mp-d) )
 412     )
 413   )
 414
 415   (if (MAIN:> (__gmpz_cmp_si mp-n 1))
 416     (if prevfact
 417       (begin
 418         (___gmpz_set_si op1 1)
 419         (push-factor mp-n op1))
 420       (begin
 421         (set 'n (get-string (__gmpz_get_str rops 10 mp-n)))
 422         (push n factors -1))))
 423   factors
 424 )
 425
 426 (define (push-factor f k)
 427   (set 'f (get-string (__gmpz_get_str rops 10 f)))
 428   (set 'k (get-string (__gmpz_get_str rops 10 k)))
 429   (if (not prevfact)
 430     (begin
 431       (push f factors -1)
 432       (set 'k (- k "1"))))
 433
 434   (while (> k "0")
 435     (push f factors -1)
 436     (set 'k (- k "1")) )
 437 )
 438
 439
 440 ; get the greates common divisor
 441 ;
 442 (define (GMP:gcd p1 p2)
 443   (if (or (not (string? p1)) (not (string? p2))) (throw-error gmp-type-error))
 444   (__gmpz_set_str op1 p1 0)
 445   (__gmpz_set_str op2 p2 0)
 446   (__gmpz_gcd rop op1 op2)
 447   (get-string (__gmpz_get_str rops 10 rop))
 448 )
 449
 450
 451 ; get binomial of <arg1> <arg2>
 452 ;
 453 (define (GMP:bin n k)
 454   (if (not (string? n)) (throw-error gmp-type-error))
 455   (__gmpz_set_str op1 n 0)
 456   (set 'k (int k 0))
 457   (__gmpz_bin_ui rop op1 k)
 458   (if (MAIN:> (__gmpz_sizeinbase rop 10) MAX_PRECISION)
 459         (throw-error gmp-size-error))
 460   (get-string (__gmpz_get_str rops 10 rop))
 461 )
 462
 463 ; get factorial of arg
 464 ; args may not be bigger than 458
 465 ; for 1024 digits precision in the result
 466 ;
 467 (define (GMP:fac p1 , n)
 468   (if (not (string? p1)) (throw-error gmp-type-error))
 469   (set 'n (int p1))
 470   (__gmpz_fac_ui rop n)
 471   (if (MAIN:> (__gmpz_sizeinbase rop 10) MAX_PRECISION)
 472         (throw-error gmp-size-error))
 473   (get-string (__gmpz_get_str rops 10 rop))
 474 )
 475
 476 ; get fibonacci of arg
 477 ; arg may not be bigger than   for a 1024 digits result
 478 ;
 479 (define (GMP:fib p1 , n)
 480   (if (not (string? p1)) (throw-error gmp-type-error))
 481   (set 'n (int p1))
 482   (__gmpz_fib_ui rop n)
 483   (if (MAIN:> (__gmpz_sizeinbase rop 10) MAX_PRECISION)
 484         (throw-error gmp-size-error))
 485   (get-string (__gmpz_get_str rops 10 rop))
 486 )
 487
 488 ; get a random number between 0 and arg - 1
 489 ;
 490 (define (GMP:rand p1)
 491   (if (not (string? p1)) (throw-error gmp-type-error))
 492   (__gmpz_set_str op1 p1 0)
 493   (__gmpz_urandomm rop randstate op1)
 494   (get-string (__gmpz_get_str rops 10 rop))
 495 )
 496
 497 ; seed the random generator
 498 ;
 499 (define (GMP:seed p1)
 500   (if (not (string? p1)) (throw-error gmp-type-error))
 501   (__gmpz_set_str op1 p1 0)
 502   (__gmp_randseed randstate op1)
 503   true
 504 )
 505
 506 (define (check func <arg1> <arg2> result)
 507    (if (MAIN:= (apply func (list <arg1> <arg2>)) result)
 508         (println "GMP:" (name func) "\t-> Ok")
 509         (println "Problem in GMP:" (name func)))
 510 )
 511
 512 (context 'MAIN)
 513
 514 ; this test assumes that libgmp itself is correct, only the
 515 ; newLISP implementation is tested
 516 ;
 517 (define (test-GMP)
 518   ; INTEGER ARITHMETIK
 519   (GMP:check 'GMP:+ "123" "456" "579")
 520   (GMP:check 'GMP:- "100" "99" "1")
 521   (GMP:check 'GMP:* "123" "456" "56088")
 522   (GMP:check 'GMP:/ "56088" "456" "123")
 523   (GMP:check 'GMP:% "56088" "456" "0")
 524   (GMP:check 'GMP:** "2" "10" "1024")
 525   (GMP:check 'GMP:= "999999" "999999" true)
 526   (GMP:check 'GMP:< "999999" "1000000" true)
 527   (GMP:check 'GMP:> "999999" "1000000" nil)
 528   (GMP:check 'GMP:<= "999999"  "1000000" true)
 529   (GMP:check 'GMP:>= "999999" "1000000" nil)
 530
 531   ; BIT OPERATTIONS
 532   (GMP:check 'GMP:& "0xAAAA" "0xFF00" "43520")
 533   (GMP:check 'GMP:| "0xAAAA" "0x5555" "65535")
 534   (GMP:check 'GMP:^ "0xAAAA" "0xAAAA" "0")
 535   (GMP:check 'GMP:~ "0xFFFF" nil "-65536")
 536
 537   ; NUMBER THEORY
 538   (GMP:check 'GMP:prime? "127" nil true)
 539   (GMP:check 'GMP:next-prime "127" nil "131")
 540   (GMP:check 'GMP:factor "123" nil '("3" "41"))
 541   (GMP:check 'GMP:gcd "20" "8" "4")
 542   (GMP:check 'GMP:bin "10" "2" "45")
 543   (GMP:check 'GMP:fac "5" nil "120")
 544   (GMP:check 'GMP:fib "30" nil "832040")
 545
 546   ; RANDOM NUMBERS
 547   (GMP:check 'GMP:seed "12345" nil true)
 548   (GMP:check 'GMP:rand "1000000" nil "18235")
 549 )
 550
 551 ; eof