home *** CD-ROM | disk | FTP | other *** search
/ Fresh Fish 4 / FreshFish_May-June1994.bin / bbs / gnu / gmp-1.3.2-src.lha / src / amiga / gmp-1.3.2 / tests / tst-mul.c < prev    next >
Encoding:
C/C++ Source or Header  |  1993-05-19  |  6.1 KB  |  262 lines
  1. /* Test mpz_add, mpz_cmp, mpz_cmp_ui, mpz_divmod, mpz_mul.
  2.  
  3. Copyright (C) 1991, 1993 Free Software Foundation, Inc.
  4.  
  5. This file is part of the GNU MP Library.
  6.  
  7. The GNU MP Library is free software; you can redistribute it and/or modify
  8. it under the terms of the GNU General Public License as published by
  9. the Free Software Foundation; either version 2, or (at your option)
  10. any later version.
  11.  
  12. The GNU MP Library is distributed in the hope that it will be useful,
  13. but WITHOUT ANY WARRANTY; without even the implied warranty of
  14. MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  15. GNU General Public License for more details.
  16.  
  17. You should have received a copy of the GNU General Public License
  18. along with the GNU MP Library; see the file COPYING.  If not, write to
  19. the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.  */
  20.  
  21. #include <stdio.h>
  22. #include "gmp.h"
  23. #include "gmp-impl.h"
  24. #include "longlong.h"
  25. #include "urandom.h"
  26.  
  27. void debug_mp ();
  28. mp_size _mpn_mul_classic ();
  29. void mpz_refmul ();
  30.  
  31. #ifndef SIZE
  32. #define SIZE 8
  33. #endif
  34.  
  35. main (argc, argv)
  36.      int argc;
  37.      char **argv;
  38. {
  39.   MP_INT multiplier, multiplicand;
  40.   MP_INT product, ref_product;
  41.   MP_INT quotient, remainder;
  42.   mp_size multiplier_size, multiplicand_size;
  43.   int i;
  44.   int reps = 100000;
  45.  
  46.   if (argc == 2)
  47.      reps = atoi (argv[1]);
  48.  
  49.   mpz_init (&multiplier);
  50.   mpz_init (&multiplicand);
  51.   mpz_init (&product);
  52.   mpz_init (&ref_product);
  53.   mpz_init ("ient);
  54.   mpz_init (&remainder);
  55.  
  56.   for (i = 0; i < reps; i++)
  57.     {
  58.       multiplier_size = urandom () % SIZE - SIZE/2;
  59.       multiplicand_size = urandom () % SIZE - SIZE/2;
  60.  
  61.       mpz_random2 (&multiplier, multiplier_size);
  62.       mpz_random2 (&multiplicand, multiplicand_size);
  63.  
  64.       mpz_mul (&product, &multiplier, &multiplicand);
  65.       mpz_refmul (&ref_product, &multiplier, &multiplicand);
  66.       if (mpz_cmp_ui (&multiplicand, 0) != 0)
  67.     mpz_divmod ("ient, &remainder, &product, &multiplicand);
  68.  
  69.       if (mpz_cmp (&product, &ref_product))
  70.     dump_abort (&multiplier, &multiplicand);
  71.  
  72.       if (mpz_cmp_ui (&multiplicand, 0) != 0)
  73.       if (mpz_cmp_ui (&remainder, 0) || mpz_cmp ("ient, &multiplier))
  74.     dump_abort (&multiplier, &multiplicand);
  75.     }
  76.  
  77.   exit (0);
  78. }
  79.  
  80. void
  81. mpz_refmul (w, u, v)
  82.      MP_INT *w;
  83.      const MP_INT *u;
  84.      const MP_INT *v;
  85. {
  86.   mp_size usize = u->size;
  87.   mp_size vsize = v->size;
  88.   mp_size wsize;
  89.   mp_size sign_product;
  90.   mp_ptr up, vp;
  91.   mp_ptr wp;
  92.   mp_ptr free_me = NULL;
  93.   size_t free_me_size;
  94.  
  95.   sign_product = usize ^ vsize;
  96.   usize = ABS (usize);
  97.   vsize = ABS (vsize);
  98.  
  99.   if (usize < vsize)
  100.     {
  101.       /* Swap U and V.  */
  102.       {const MP_INT *t = u; u = v; v = t;}
  103.       {mp_size t = usize; usize = vsize; vsize = t;}
  104.     }
  105.  
  106.   up = u->d;
  107.   vp = v->d;
  108.   wp = w->d;
  109.  
  110.   /* Ensure W has space enough to store the result.  */
  111.   wsize = usize + vsize;
  112.   if (w->alloc < wsize)
  113.     {
  114.       if (wp == up || wp == vp)
  115.     {
  116.       free_me = wp;
  117.       free_me_size = w->alloc;
  118.     }
  119.       else
  120.     (*_mp_free_func) (wp, w->alloc * BYTES_PER_MP_LIMB);
  121.  
  122.       w->alloc = wsize;
  123.       wp = (mp_ptr) (*_mp_allocate_func) (wsize * BYTES_PER_MP_LIMB);
  124.       w->d = wp;
  125.     }
  126.   else
  127.     {
  128.       /* Make U and V not overlap with W.  */
  129.       if (wp == up)
  130.     {
  131.       /* W and U are identical.  Allocate temporary space for U.  */
  132.       up = (mp_ptr) alloca (usize * BYTES_PER_MP_LIMB);
  133.       /* Is V identical too?  Keep it identical with U.  */
  134.       if (wp == vp)
  135.         vp = up;
  136.       /* Copy to the temporary space.  */
  137.       MPN_COPY (up, wp, usize);
  138.     }
  139.       else if (wp == vp)
  140.     {
  141.       /* W and V are identical.  Allocate temporary space for V.  */
  142.       vp = (mp_ptr) alloca (vsize * BYTES_PER_MP_LIMB);
  143.       /* Copy to the temporary space.  */
  144.       MPN_COPY (vp, wp, vsize);
  145.     }
  146.     }
  147.  
  148.   wsize = _mpn_mul_classic (wp, up, usize, vp, vsize);
  149.   w->size = sign_product < 0 ? -wsize : wsize;
  150.   if (free_me != NULL)
  151.     (*_mp_free_func) (free_me, free_me_size * BYTES_PER_MP_LIMB);
  152.  
  153.   alloca (0);
  154. }
  155.  
  156. mp_size
  157. _mpn_mul_classic (prodp, up, usize, vp, vsize)
  158.      mp_ptr prodp;
  159.      mp_srcptr up;
  160.      mp_size usize;
  161.      mp_srcptr vp;
  162.      mp_size vsize;
  163. {
  164.   mp_size n;
  165.   mp_size prod_size;
  166.   mp_limb cy;
  167.   mp_size i, j;
  168.   mp_limb prod_low, prod_high;
  169.   mp_limb cy_dig;
  170.   mp_limb v_limb, c;
  171.  
  172.   if (vsize == 0)
  173.     return 0;
  174.  
  175.   /* Offset UP and PRODP so that the inner loop can be faster.  */
  176.   up += usize;
  177.   prodp += usize;
  178.  
  179.   /* Multiply by the first limb in V separately, as the result can
  180.      be stored (not added) to PROD.  We also avoid a loop for zeroing.  */
  181.   v_limb = vp[0];
  182.   cy_dig = 0;
  183.   j = -usize;
  184.   do
  185.     {
  186.       umul_ppmm (prod_high, prod_low, up[j], v_limb);
  187.       add_ssaaaa (cy_dig, prodp[j], prod_high, prod_low, 0, cy_dig);
  188.       j++;
  189.     }
  190.   while (j < 0);
  191.  
  192.   prodp[j] = cy_dig;
  193.   prodp++;
  194.  
  195.   /* For each iteration in the outer loop, multiply one limb from
  196.      U with one limb from V, and add it to PROD.  */
  197.   for (i = 1; i < vsize; i++)
  198.     {
  199.       v_limb = vp[i];
  200.       cy_dig = 0;
  201.       j = -usize;
  202.  
  203.       /* Inner loops.  Simulate the carry flag by jumping between
  204.      these loops.  The first is used when there was no carry
  205.      in the previois iteration; the second when there was carry.  */
  206.  
  207.       do
  208.     {
  209.       umul_ppmm (prod_high, prod_low, up[j], v_limb);
  210.       add_ssaaaa (cy_dig, prod_low, prod_high, prod_low, 0, cy_dig);
  211.       c = prodp[j];
  212.       prod_low += c;
  213.       prodp[j] = prod_low;
  214.       if (prod_low < c)
  215.         goto cy_loop;
  216.     ncy_loop:
  217.       j++;
  218.     }
  219.       while (j < 0);
  220.  
  221.       prodp[j] = cy_dig;
  222.       prodp++;
  223.       continue;
  224.  
  225.       do
  226.     {
  227.       umul_ppmm (prod_high, prod_low, up[j], v_limb);
  228.       add_ssaaaa (cy_dig, prod_low, prod_high, prod_low, 0, cy_dig);
  229.       c = prodp[j];
  230.       prod_low += c + 1;
  231.       prodp[j] = prod_low;
  232.       if (prod_low > c)
  233.         goto ncy_loop;
  234.     cy_loop:
  235.       j++;
  236.     }
  237.       while (j < 0);
  238.  
  239.       cy_dig += 1;
  240.       prodp[j] = cy_dig;
  241.       prodp++;
  242.     }
  243.  
  244.   return usize + vsize - (cy_dig == 0);
  245. }
  246.  
  247. dump_abort (multiplier, multiplicand)
  248.      MP_INT *multiplier, *multiplicand;
  249. {
  250.   fprintf (stderr, "ERROR\n");
  251.   fprintf (stderr, "multiplier = "); debug_mp (multiplier, -16);
  252.   fprintf (stderr, "multiplicand  = "); debug_mp (multiplicand, -16);
  253.   abort();
  254. }
  255.  
  256. void
  257. debug_mp (x, base)
  258.      MP_INT *x;
  259. {
  260.   mpz_out_str (stderr, base, x); fputc ('\n', stderr);
  261. }
  262.