t-matrix22.c
上传用户:qaz666999
上传日期:2022-08-06
资源大小:2570k
文件大小:5k
源码类别:

数学计算

开发平台:

Unix_Linux

  1. /* Tests matrix22_mul.
  2. Copyright 2008 Free Software Foundation, Inc.
  3. This file is part of the GNU MP Library.
  4. The GNU MP Library is free software; you can redistribute it and/or modify
  5. it under the terms of the GNU Lesser General Public License as published by
  6. the Free Software Foundation; either version 3 of the License, or (at your
  7. option) any later version.
  8. The GNU MP Library is distributed in the hope that it will be useful, but
  9. WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  10. or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
  11. License for more details.
  12. You should have received a copy of the GNU Lesser General Public License
  13. along with the GNU MP Library.  If not, see http://www.gnu.org/licenses/.  */
  14. #include <stdio.h>
  15. #include <stdlib.h>
  16. #include "gmp.h"
  17. #include "gmp-impl.h"
  18. #include "tests.h"
  19. struct matrix {
  20.   mp_size_t alloc;
  21.   mp_size_t n;
  22.   mp_ptr e00, e01, e10, e11;
  23. };
  24. static void
  25. matrix_init (struct matrix *M, mp_size_t n)
  26. {
  27.   mp_ptr p = refmpn_malloc_limbs (4*(n+1));
  28.   M->e00 = p; p += n+1;
  29.   M->e01 = p; p += n+1;
  30.   M->e10 = p; p += n+1;
  31.   M->e11 = p;
  32.   M->alloc = n + 1;
  33.   M->n = 0;
  34. }
  35. static void
  36. matrix_clear (struct matrix *M)
  37. {
  38.   refmpn_free_limbs (M->e00);
  39. }
  40. static void
  41. matrix_copy (struct matrix *R, const struct matrix *M)
  42. {
  43.   R->n = M->n;
  44.   MPN_COPY (R->e00, M->e00, M->n);
  45.   MPN_COPY (R->e01, M->e01, M->n);
  46.   MPN_COPY (R->e10, M->e10, M->n);
  47.   MPN_COPY (R->e11, M->e11, M->n);
  48. }
  49. /* Used with same size, so no need for normalization. */
  50. static int
  51. matrix_equal_p (const struct matrix *A, const struct matrix *B)
  52. {
  53.   return (A->n == B->n
  54.   && mpn_cmp (A->e00, B->e00, A->n) == 0
  55.   && mpn_cmp (A->e01, B->e01, A->n) == 0
  56.   && mpn_cmp (A->e10, B->e10, A->n) == 0
  57.   && mpn_cmp (A->e11, B->e11, A->n) == 0);
  58. }
  59. static void
  60. matrix_random(struct matrix *M, mp_size_t n, gmp_randstate_ptr rands)
  61. {
  62.   M->n = n;
  63.   mpn_random (M->e00, n);
  64.   mpn_random (M->e01, n);
  65.   mpn_random (M->e10, n);
  66.   mpn_random (M->e11, n);
  67. }
  68. #define MUL(rp, ap, an, bp, bn) do { 
  69.     if (an > bn)      
  70.       mpn_mul (rp, ap, an, bp, bn);  
  71.     else      
  72.       mpn_mul (rp, bp, bn, ap, an);  
  73.   } while(0)
  74. static void
  75. ref_matrix22_mul (struct matrix *R,
  76.   const struct matrix *A,
  77.   const struct matrix *B, mp_ptr tp)
  78. {
  79.   mp_size_t an, bn, n;
  80.   mp_ptr r00, r01, r10, r11, a00, a01, a10, a11, b00, b01, b10, b11;
  81.   if (A->n >= B->n)
  82.     {
  83.       r00 = R->e00; a00 = A->e00; b00 = B->e00;
  84.       r01 = R->e01; a01 = A->e01; b01 = B->e01;
  85.       r10 = R->e10; a10 = A->e10; b10 = B->e10;
  86.       r11 = R->e11; a11 = A->e11; b11 = B->e11;
  87.       an = A->n, bn = B->n;
  88.     }
  89.   else
  90.     {
  91.       /* Transpose */
  92.       r00 = R->e00; a00 = B->e00; b00 = A->e00;
  93.       r01 = R->e10; a01 = B->e10; b01 = A->e10;
  94.       r10 = R->e01; a10 = B->e01; b10 = A->e01;
  95.       r11 = R->e11; a11 = B->e11; b11 = A->e11;
  96.       an = B->n, bn = A->n;
  97.     }
  98.   n = an + bn;
  99.   R->n = n + 1;
  100.   mpn_mul (r00, a00, an, b00, bn);
  101.   mpn_mul (tp, a01, an, b10, bn);
  102.   r00[n] = mpn_add_n (r00, r00, tp, n);
  103.   mpn_mul (r01, a00, an, b01, bn);
  104.   mpn_mul (tp, a01, an, b11, bn);
  105.   r01[n] = mpn_add_n (r01, r01, tp, n);
  106.   mpn_mul (r10, a10, an, b00, bn);
  107.   mpn_mul (tp, a11, an, b10, bn);
  108.   r10[n] = mpn_add_n (r10, r10, tp, n);
  109.   mpn_mul (r11, a10, an, b01, bn);
  110.   mpn_mul (tp, a11, an, b11, bn);
  111.   r11[n] = mpn_add_n (r11, r11, tp, n);
  112. }
  113. static void
  114. one_test (const struct matrix *A, const struct matrix *B, int i)
  115. {
  116.   struct matrix R;
  117.   struct matrix P;
  118.   mp_ptr tp;
  119.   matrix_init (&R, A->n + B->n + 1);
  120.   matrix_init (&P, A->n + B->n + 1);
  121.   tp = refmpn_malloc_limbs (mpn_matrix22_mul_itch (A->n, B->n));
  122.   ref_matrix22_mul (&R, A, B, tp);
  123.   matrix_copy (&P, A);
  124.   mpn_matrix22_mul (P.e00, P.e01, P.e10, P.e11, A->n,
  125.     B->e00, B->e01, B->e10, B->e11, B->n, tp);
  126.   P.n = A->n + B->n + 1;
  127.   if (!matrix_equal_p (&R, &P))
  128.     {
  129.       fprintf (stderr, "ERROR in test %dn", i);
  130.       gmp_fprintf (stderr, "A = (%Nx, %Nxn      %Nx, %Nx)n"
  131.    "B = (%Nx, %Nxn      %Nx, %Nx)n"
  132.    "R = (%Nx, %Nx (expected)n      %Nx, %Nx)n"
  133.    "P = (%Nx, %Nx (incorrect)n      %Nx, %Nx)n",
  134.    A->e00, A->n, A->e01, A->n, A->e10, A->n, A->e11, A->n,
  135.    B->e00, B->n, B->e01, B->n, B->e10, B->n, B->e11, B->n,
  136.    R.e00, R.n, R.e01, R.n, R.e10, R.n, R.e11, R.n,
  137.    P.e00, P.n, P.e01, P.n, P.e10, P.n, P.e11, P.n);
  138.       abort();
  139.     }
  140.   refmpn_free_limbs (tp);
  141.   matrix_clear (&R);
  142.   matrix_clear (&P);
  143. }
  144. #define MAX_SIZE (2+2*MATRIX22_STRASSEN_THRESHOLD)
  145. int
  146. main (int argc, char **argv)
  147. {
  148.   struct matrix A;
  149.   struct matrix B;
  150.   gmp_randstate_ptr rands;
  151.   mpz_t bs;
  152.   int i;
  153.   tests_start ();
  154.   rands = RANDS;
  155.   matrix_init (&A, MAX_SIZE);
  156.   matrix_init (&B, MAX_SIZE);
  157.   mpz_init (bs);
  158.   for (i = 0; i < 1000; i++)
  159.     {
  160.       mp_size_t an, bn;
  161.       mpz_urandomb (bs, rands, 32);
  162.       an = 1 + mpz_get_ui (bs) % MAX_SIZE;
  163.       mpz_urandomb (bs, rands, 32);
  164.       bn = 1 + mpz_get_ui (bs) % MAX_SIZE;
  165.       matrix_random (&A, an, rands);
  166.       matrix_random (&B, bn, rands);
  167.       one_test (&A, &B, i);
  168.     }
  169.   mpz_clear (bs);
  170.   matrix_clear (&A);
  171.   matrix_clear (&B);
  172.   tests_end ();
  173.   return 0;
  174. }