chiark - git - mdw - catacomb-perl/blob - mp.xs

   1 # ---?---
   2 #
   3 # $Id$
   4 #
   5 # Multiprecision interface
   6 #
   7 # (c) 2000 Straylight/Edgeware
   8 #
   9
  10 #----- Licensing notice -----------------------------------------------------
  11 #
  12 # This file is part of the Perl interface to Catacomb.
  13 #
  14 # Catacomb/Perl is free software; you can redistribute it and/or modify
  15 # it under the terms of the GNU General Public License as published by
  16 # the Free Software Foundation; either version 2 of the License, or
  17 # (at your option) any later version.
  18 #
  19 # Catacomb/Perl is distributed in the hope that it will be useful,
  20 # but WITHOUT ANY WARRANTY; without even the implied warranty of
  21 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  22 # GNU General Public License for more details.
  23 #
  24 # You should have received a copy of the GNU General Public License
  25 # along with Catacomb/Perl; if not, write to the Free Software Foundation,
  26 # Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
  27 #
  28
  29 MODULE = Catacomb PACKAGE = Catacomb::MP PREFIX = mp_
  30
  31 mp *
  32 new(me, sv = 0, radix = 0)
  33         SV *me
  34         SV *sv
  35         int radix
  36         CODE:
  37         RETVAL = sv ? mp_fromsv(sv, "sv", radix, 1) : MP_ZERO;
  38         OUTPUT:
  39         RETVAL
  40
  41 mp *
  42 copy(x)
  43         mp *x
  44         CODE:
  45         RETVAL = MP_COPY(x);
  46         OUTPUT:
  47         RETVAL
  48
  49 mp *
  50 loadb(me, sv)
  51         SV *me
  52         SV *sv
  53         PREINIT:
  54         char *p;
  55         STRLEN len;
  56         CODE:
  57         p = SvPV(sv, len);
  58         RETVAL = mp_loadb(MP_NEW, p, len);
  59         OUTPUT:
  60         RETVAL
  61
  62 mp *
  63 loadl(me, sv)
  64         SV *me
  65         SV *sv
  66         PREINIT:
  67         char *p;
  68         STRLEN len;
  69         CODE:
  70         p = SvPV(sv, len);
  71         RETVAL = mp_loadl(MP_NEW, p, len);
  72         OUTPUT:
  73         RETVAL
  74
  75 mp *
  76 loadb2c(me, sv)
  77         SV *me
  78         SV *sv
  79         PREINIT:
  80         char *p;
  81         STRLEN len;
  82         CODE:
  83         p = SvPV(sv, len);
  84         RETVAL = mp_loadb2c(MP_NEW, p, len);
  85         OUTPUT:
  86         RETVAL
  87
  88 mp *
  89 loadl2c(me, sv)
  90         SV *me
  91         SV *sv
  92         PREINIT:
  93         char *p;
  94         STRLEN len;
  95         CODE:
  96         p = SvPV(sv, len);
  97         RETVAL = mp_loadl2c(MP_NEW, p, len);
  98         OUTPUT:
  99         RETVAL
 100
 101 int
 102 metrics(m)
 103         mp *m
 104         INTERFACE_MACRO:
 105         XSINTERFACE_FUNC
 106         XSINTERFACE_FUNC_SETMP
 107         INTERFACE:
 108         octets bits octets2c
 109
 110 SV *
 111 storeb(m, i = -1)
 112         mp *m
 113         int i
 114         PREINIT:
 115         size_t sz;
 116         CODE:
 117         if (i >= 0)
 118           sz = i;
 119         else {
 120           sz = mp_octets(m);
 121           if (!sz)
 122             sz = 1;
 123         }
 124         RETVAL = NEWSV(0, sz ? sz : 1);
 125         mp_storeb(m, SvPVX(RETVAL), sz);
 126         SvCUR_set(RETVAL, sz);
 127         SvPOK_on(RETVAL);
 128         OUTPUT:
 129         RETVAL
 130
 131 SV *
 132 storel(m, i = -1)
 133         mp *m
 134         int i
 135         PREINIT:
 136         size_t sz;
 137         CODE:
 138         sz = (i < 0) ? mp_octets(m) : i;
 139         RETVAL = NEWSV(0, sz ? sz : 1);
 140         mp_storel(m, SvPVX(RETVAL), sz);
 141         SvCUR_set(RETVAL, sz);
 142         SvPOK_on(RETVAL);
 143         OUTPUT:
 144         RETVAL
 145
 146 SV *
 147 storeb2c(m, i = -1)
 148         mp *m
 149         int i
 150         PREINIT:
 151         size_t sz;
 152         CODE:
 153         sz = (i < 0) ? mp_octets2c(m) : i;
 154         RETVAL = NEWSV(0, sz ? sz : 1);
 155         mp_storeb(m, SvPVX(RETVAL), sz);
 156         SvCUR_set(RETVAL, sz);
 157         SvPOK_on(RETVAL);
 158         OUTPUT:
 159         RETVAL
 160
 161 SV *
 162 storel2c(m, i = -1)
 163         mp *m
 164         int i
 165         PREINIT:
 166         size_t sz;
 167         CODE:
 168         sz = (i < 0) ? mp_octets2c(m) : i;
 169         RETVAL = NEWSV(0, sz ? sz : 1);
 170         mp_storel(m, SvPVX(RETVAL), sz);
 171         SvCUR_set(RETVAL, sz);
 172         SvPOK_on(RETVAL);
 173         OUTPUT:
 174         RETVAL
 175
 176 SV *
 177 tostring(m, radix = 10)
 178         mp *m
 179         int radix
 180         CODE:
 181         RETVAL = NEWSV(0, 0);
 182         mp_writesv(m, RETVAL, radix);
 183         OUTPUT:
 184         RETVAL
 185
 186 void
 187 fromstring(me, s, radix = 10)
 188         SV *me
 189         SV *s
 190         int radix
 191         PREINIT:
 192         mptext_stringctx ms;
 193         STRLEN len;
 194         mp *x;
 195         PPCODE:
 196         ms.buf = SvPV(s, len);
 197         ms.lim = ms.buf + len;
 198         x = mp_read(MP_NEW, radix, &mptext_stringops, &ms);
 199         if (x) {
 200           XPUSHs(RET_MP(x));
 201           if (GIMME_V == G_ARRAY)
 202             XPUSHs(sv_2mortal(newSVpvn(ms.buf, ms.lim - ms.buf)));
 203         }
 204
 205 SV *
 206 toint(m)
 207         mp *m
 208         CODE:
 209         RETVAL = newSViv(mp_toiv(m));
 210         OUTPUT:
 211         RETVAL
 212
 213 SV *
 214 DESTROY(m)
 215         mp *m
 216         CODE:
 217         mp_drop(m);
 218         XSRETURN_UNDEF;
 219
 220 mp *
 221 unop(a)
 222         mp *a
 223         C_ARGS:
 224         MP_NEW, a
 225         INTERFACE_MACRO:
 226         XSINTERFACE_FUNC
 227         XSINTERFACE_FUNC_SETMP
 228         INTERFACE:
 229         not not2c sqr sqrt
 230
 231 mp *
 232 neg(a)
 233         mp *a
 234         CODE:
 235         MP_COPY(a);
 236         RETVAL = mp_split(a);
 237         if (RETVAL->v < RETVAL->vl)
 238           RETVAL->f ^= MP_NEG;
 239         OUTPUT:
 240         RETVAL
 241
 242 mp *
 243 mp_factorial(me, x)
 244         SV *me
 245         IV x
 246         C_ARGS:
 247         x
 248
 249 mp *
 250 binop(a, b)
 251         mp *a
 252         mp *b
 253         C_ARGS:
 254         MP_NEW, a, b
 255         INTERFACE_MACRO:
 256         XSINTERFACE_FUNC
 257         XSINTERFACE_FUNC_SETMP
 258         INTERFACE:
 259         add sub mul and2c or2c nand2c nor2c xor2c and or nand nor xor exp
 260
 261 mp *
 262 shiftop(a, n)
 263         mp *a
 264         int n
 265         C_ARGS:
 266         MP_NEW, a, n
 267         INTERFACE_MACRO:
 268         XSINTERFACE_FUNC
 269         XSINTERFACE_FUNC_SETMP
 270         INTERFACE:
 271         lsl lsr lsl2c lsr2c
 272
 273 bool
 274 testbitop(a, n)
 275         mp *a
 276         unsigned long n
 277         INTERFACE_MACRO:
 278         XSINTERFACE_FUNC
 279         XSINTERFACE_FUNC_SETMP
 280         INTERFACE:
 281         testbit testbit2c
 282
 283 mp *
 284 flipbits(a, n)
 285         mp *a
 286         unsigned long n
 287         C_ARGS:
 288         MP_NEW, a, n
 289         INTERFACE_MACRO:
 290         XSINTERFACE_FUNC
 291         XSINTERFACE_FUNC_SETMP
 292         INTERFACE:
 293         setbit clearbit setbit2c clearbit2c
 294
 295 int
 296 mp_cmp(a, b)
 297         mp *a
 298         mp *b
 299
 300 int
 301 mp_eq(a, b)
 302         mp *a
 303         mp *b
 304
 305 int
 306 jacobi(a, n)
 307         mp *a
 308         mp *n
 309         CODE:
 310         if (!MP_LEN(n) || !(n->v[0] & 1))
 311           croak("n must be odd in Catacomb::MP::jacobi");
 312         RETVAL = mp_jacobi(a, n);
 313         OUTPUT:
 314         RETVAL
 315
 316 mp *
 317 mp_modsqrt(p, x)
 318         mp *p
 319         mp *x
 320         INIT:
 321         if (!MP_POSP(p) || !MP_ODDP(p))
 322           croak("p is not positive and odd");
 323         if (mp_jacobi(x, p) != 1)
 324           croak("x not a quadratic residue mod p");
 325         C_ARGS:
 326         MP_NEW, x, p
 327
 328 void
 329 div(a, b)
 330         mp *a
 331         mp *b
 332         PREINIT:
 333         mp *q = MP_NEW, *r = MP_NEW;
 334         PPCODE:
 335         if (MP_EQ(b, MP_ZERO))
 336           croak("Divide by zero in Catacomb::MP::div");
 337         q = MP_NEW;
 338         switch (GIMME_V) {
 339           case G_ARRAY:
 340             r = MP_NEW;
 341             mp_div(&q, &r, a, b);
 342             EXTEND(SP, 2);
 343             PUSHs(RET_MP(q));
 344             PUSHs(RET_MP(r));
 345             break;
 346           case G_VOID:
 347             break;
 348           default:
 349             mp_div(&q, 0, a, b);
 350             EXTEND(SP, 1);
 351             PUSHs(RET_MP(q));
 352             break;
 353         }
 354
 355 void
 356 gcd(a, b)
 357         mp *a
 358         mp *b
 359         PREINIT:
 360         mp *g = MP_NEW, *x = MP_NEW, *y = MP_NEW;
 361         PPCODE:
 362         switch (GIMME_V) {
 363           case G_ARRAY:
 364             mp_gcd(&g, &x, &y, a, b);
 365             EXTEND(SP, 3);
 366             PUSHs(RET_MP(g));
 367             PUSHs(RET_MP(x));
 368             PUSHs(RET_MP(y));
 369             break;
 370           case G_VOID:
 371             break;
 372           default:
 373             mp_gcd(&g, 0, 0, a, b);
 374             EXTEND(SP, 1);
 375             PUSHs(RET_MP(g));
 376             break;
 377         }
 378
 379 void
 380 odd(m)
 381         mp *m
 382         PREINIT:
 383         mp *t;
 384         size_t s;
 385         PPCODE:
 386         t = mp_odd(MP_NEW, m, &s);
 387         EXTEND(SP, 2);
 388         PUSHs(sv_2mortal(newSViv(s)));
 389         PUSHs(RET_MP(t));
 390
 391 MODULE = Catacomb PACKAGE = Catacomb::MP PREFIX = pgen_
 392
 393 bool
 394 pgen_primep(x, r = &rand_global)
 395         mp *x
 396         grand *r
 397
 398 MODULE = Catacomb PACKAGE = Catacomb::MP PREFIX = pfilt_
 399
 400 int
 401 pfilt_smallfactor(x)
 402         mp *x
 403
 404 MODULE = Catacomb PACKAGE = Catacomb::MP::Mont PREFIX = mpmont_
 405
 406 MP_Mont *
 407 new(me, x)
 408         SV *me
 409         mp *x
 410         CODE:
 411         RETVAL = CREATE(MP_Mont);
 412         if (mpmont_create(RETVAL, x)) {
 413           DESTROY(RETVAL);
 414           RETVAL = 0;
 415         }
 416         OUTPUT:
 417         RETVAL
 418
 419 SV *
 420 DESTROY(mm)
 421         MP_Mont *mm
 422         CODE:
 423         mpmont_destroy(mm);
 424         DESTROY(mm);
 425         XSRETURN_UNDEF;
 426
 427 mp *
 428 mpmont_reduce(mm, x)
 429         MP_Mont *mm
 430         mp *x
 431         C_ARGS:
 432         mm, MP_NEW, x
 433
 434 mp *
 435 mpmont_mul(mm, x, y)
 436         MP_Mont *mm
 437         mp *x
 438         mp *y
 439         C_ARGS:
 440         mm, MP_NEW, x, y
 441
 442 mp *
 443 in(mm, x)
 444         MP_Mont *mm
 445         mp *x
 446         CODE:
 447         RETVAL = MP_NEW;
 448         mp_div(0, &RETVAL, x, mm->m);
 449         RETVAL = mpmont_mul(mm, RETVAL, RETVAL, mm->r2);
 450         OUTPUT:
 451         RETVAL
 452
 453 mp *
 454 mpmont_expr(mm, g, x)
 455         MP_Mont *mm
 456         mp *g
 457         mp *x
 458         C_ARGS:
 459         mm, MP_NEW, g, x
 460
 461 mp *
 462 mpmont_exp(mm, g, x)
 463         MP_Mont *mm
 464         mp *g
 465         mp *x
 466         C_ARGS:
 467         mm, MP_NEW, g, x
 468
 469 mp *
 470 mpmont_mexpr(mm, ...)
 471         MP_Mont *mm
 472         PREINIT:
 473         mp_expfactor *v;
 474         size_t i, j, n;
 475         CODE:
 476         if (items < 3 || !(items & 1)) {
 477           croak("Usage: Catacomb::MP::Mont::mexpr"
 478                 "(mm, g_0, x_0, g_1, x_1, ...");
 479         }
 480         n = (items - 1)/2;
 481         v = xmalloc(n * sizeof(mp_expfactor));
 482         for (i = 1, j = 0; i < items; i += 2, j++) {
 483           v[j].base = mp_fromsv(ST(i), "g_i", 0, 0);
 484           v[j].exp = mp_fromsv(ST(i + 1), "x_i", 0, 0);
 485         }
 486         RETVAL = mpmont_mexpr(mm, MP_NEW, v, n);
 487         xfree(v);
 488         OUTPUT:
 489         RETVAL
 490
 491 mp *
 492 mpmont_mexp(mm, ...)
 493         MP_Mont *mm
 494         PREINIT:
 495         mp_expfactor *v;
 496         size_t i, j, n;
 497         CODE:
 498         if (items < 3 || !(items & 1)) {
 499           croak("Usage: Catacomb::MP::Mont::mexp"
 500                 "(mm, g_0, x_0, g_1, x_1, ...");
 501         }
 502         n = (items - 1)/2;
 503         v = xmalloc(n * sizeof(mp_expfactor));
 504         for (i = 1, j = 0; i < items; i += 2, j++) {
 505           v[j].base = mp_fromsv(ST(i), "g_%lu", 0, 0, (unsigned long)i);
 506           v[j].exp = mp_fromsv(ST(i + 1), "x_%lu", 0, 0, (unsigned long)i);
 507         }
 508         RETVAL = mpmont_mexp(mm, MP_NEW, v, n);
 509         xfree(v);
 510         OUTPUT:
 511         RETVAL
 512
 513 mp *
 514 r(mm)
 515         MP_Mont *mm
 516         CODE:
 517         RETVAL = MP_COPY(mm->r);
 518         OUTPUT:
 519         RETVAL
 520
 521 mp *
 522 r2(mm)
 523         MP_Mont *mm
 524         CODE:
 525         RETVAL = MP_COPY(mm->r2);
 526         OUTPUT:
 527         RETVAL
 528
 529 mp *
 530 m(mm)
 531         MP_Mont *mm
 532         CODE:
 533         RETVAL = MP_COPY(mm->m);
 534         OUTPUT:
 535         RETVAL
 536
 537 MODULE = Catacomb PACKAGE = Catacomb::MP::Barrett PREFIX = mpbarrett_
 538
 539 MP_Barrett *
 540 new(me, x)
 541         SV *me
 542         mp *x
 543         CODE:
 544         RETVAL = CREATE(mpbarrett);
 545         if (mpbarrett_create(RETVAL, x)) {
 546           DESTROY(RETVAL);
 547           RETVAL = 0;
 548         }
 549         OUTPUT:
 550         RETVAL
 551
 552 SV *
 553 DESTROY(mb)
 554         MP_Barrett *mb
 555         CODE:
 556         mpbarrett_destroy(mb);
 557         DESTROY(mb);
 558         XSRETURN_UNDEF;
 559
 560 mp *
 561 mpbarrett_reduce(mb, x)
 562         MP_Barrett *mb
 563         mp *x
 564         C_ARGS:
 565         mb, MP_NEW, x
 566
 567 mp *
 568 mpbarrett_exp(mb, g, x)
 569         MP_Barrett *mb
 570         mp *g
 571         mp *x
 572         C_ARGS:
 573         mb, MP_NEW, g, x
 574
 575 mp *
 576 m(mb)
 577         MP_Barrett *mb
 578         CODE:
 579         RETVAL = MP_COPY(mb->m);
 580         OUTPUT:
 581         RETVAL
 582
 583 MODULE = Catacomb PACKAGE = Catacomb::MP::Reduce PREFIX = mpreduce_
 584
 585 MP_Reduce *
 586 new(me, x)
 587         SV *me
 588         mp *x
 589         CODE:
 590         RETVAL = CREATE(MP_Reduce);
 591         if (mpreduce_create(RETVAL, x)) {
 592           DESTROY(RETVAL);
 593           RETVAL = 0;
 594         }
 595         OUTPUT:
 596         RETVAL
 597
 598 SV *
 599 DESTROY(r)
 600         MP_Reduce *r
 601         CODE:
 602         mpreduce_destroy(r);
 603         DESTROY(r);
 604         XSRETURN_UNDEF;
 605
 606 mp *
 607 reduce(r, x)
 608         MP_Reduce *r
 609         mp *x
 610         CODE:
 611         RETVAL = mpreduce_do(r, MP_NEW, x);
 612         OUTPUT:
 613         RETVAL
 614
 615 mp *
 616 mpreduce_exp(r, x, y)
 617         MP_Reduce *r
 618         mp *x
 619         mp *y
 620         C_ARGS:
 621         r, MP_NEW, x, y
 622
 623 mp *
 624 m(r)
 625         MP_Reduce *r
 626         CODE:
 627         RETVAL = MP_COPY(r->p);
 628         OUTPUT:
 629         RETVAL
 630
 631 MODULE = Catacomb PACKAGE = Catacomb::MP::CRT
 632
 633 MP_CRT *
 634 new(me, ...)
 635         SV *me
 636         PREINIT:
 637         mpcrt_mod *v;
 638         size_t n, i;
 639         CODE:
 640         if (items < 1)
 641           croak("Usage: Catacomb::MP::CRT::new(me, n_0, n_1, ...)");
 642         n = items - 1;
 643         v = xmalloc(n * sizeof(mpcrt_mod));
 644         for (i = 0; i < n; i++) {
 645           v[i].m = mp_copy(mp_fromsv(ST(i + 1), "n_%lu",
 646                            0, 0, (unsigned long)i));
 647           v[i].n = v[i].ni = v[i].nni = 0;
 648         }
 649         RETVAL = CREATE(MP_CRT);
 650         mpcrt_create(RETVAL, v, n, 0);
 651         OUTPUT:
 652         RETVAL
 653
 654 mp *
 655 product(mc)
 656         MP_CRT *mc
 657         CODE:
 658         RETVAL = MP_COPY(mc->mb.m);
 659         OUTPUT:
 660         RETVAL
 661
 662 void
 663 moduli(mc)
 664         MP_CRT *mc
 665         PREINIT:
 666         size_t n, i;
 667         PPCODE:
 668         n = mc->k;
 669         if (GIMME_V == G_SCALAR)
 670           XPUSHs(sv_2mortal(newSViv(n)));
 671         else for (i = 0; i < n; i++)
 672           XPUSHs(RET_MP(MP_COPY(mc->v[i].m)));
 673
 674 SV *
 675 DESTROY(mc)
 676         MP_CRT *mc
 677         CODE:
 678         mpcrt_destroy(mc);
 679         xfree(mc->v);
 680         DESTROY(mc);
 681         XSRETURN_UNDEF;
 682
 683 mp *
 684 solve(mc, ...)
 685         MP_CRT *mc
 686         PREINIT:
 687         mp **v;
 688         size_t n, i;
 689         CODE:
 690         n = mc->k;
 691         if (items - 1 != n)
 692           croak("Wrong number of residues for this CRT context");
 693         v = xmalloc(n * sizeof(mp *));
 694         for (i = 0; i < n; i++)
 695           v[i] = mp_fromsv(ST(i + 1), "r_%lu", 0, 0, (unsigned long)i);
 696         RETVAL = mpcrt_solve(mc, MP_NEW, v);
 697         xfree(v);
 698         OUTPUT:
 699         RETVAL
 700
 701 #----- That's all, folks ----------------------------------------------------