mpir/tests/mpn/t-redc_basecase.c

/* Test mpn_redc_basecase

  Copyright 2009 Jason Moxham

  This file is part of the MPIR Library.

  The MPIR Library is free software; you can redistribute it and/or modify
  it under the terms of the GNU Lesser General Public License as published
  by the Free Software Foundation; either version 2.1 of the License, or (at
  your option) any later version.

  The MPIR Library is distributed in the hope that it will be useful, but
  WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
  License for more details.

  You should have received a copy of the GNU Lesser General Public License
  along with the MPIR Library; see the file COPYING.LIB.  If not, write
  to the Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
  Boston, MA 02110-1301, USA.

*/

#include <stdio.h>
#include <stdlib.h>
#include "mpir.h"
#include "gmp-impl.h"
#include "longlong.h"
#include "tests.h"

/* Set cp[] <- tp[]/R^n mod mp[].  Clobber tp[].
   mp[] is n limbs; tp[] is 2n limbs.  */
void ref_redc_basecase (mp_ptr cp, mp_srcptr mp, mp_size_t n, mp_limb_t Nprim, mp_ptr tp)
{ mp_limb_t cy;  mp_limb_t q;  mp_size_t j;
  ASSERT_MPN (tp, 2*n);
  for (j = 0; j < n; j++)
    { q = (tp[0] * Nprim) & GMP_NUMB_MASK;
      tp[0] = mpn_addmul_1 (tp, mp, n, q);
      tp++;  }
  cy = mpn_add_n (cp, tp, tp - n, n);
  if (cy != 0)  mpn_sub_n (cp, cp, mp, n);
}

int
main (void)
{ 
  gmp_randstate_ptr rands;
  int j,n;
  mp_limb_t cp1[1000], cp2[1000], mp[1000],tp1[1000],tp2[1000],inv;
  tests_start ();
  rands = RANDS;

  for(n=1;n<100;n++)
     {for(j=1;j<100;j++)
         {mpn_random(mp,n);mp[0]|=1;modlimb_invert(inv,mp[0]);inv=-inv;mpn_random(tp1,2*n);
          MPN_COPY(tp2,tp1,2*n);
          ref_redc_basecase(cp1,mp,n,inv,tp1);
          mpn_redc_basecase(cp2,mp,n,inv,tp2);
          if(mpn_cmp(cp1,cp2,n)!=0){printf("mpn_redc_basecase error %d\n",n);abort();}
          if(n!=1 && mpn_cmp(tp1,tp2,2*n)!=0){printf("mpn_redc_basecase possible error\n");abort();}
          /* we dont require the above to be the same but it could be a useful test */      
         }
     }
  for(n=1;n<100;n++)
     {for(j=1;j<100;j++)
         {mpn_random2(mp,n);mp[0]|=1;modlimb_invert(inv,mp[0]);inv=-inv;mpn_random2(tp1,2*n);
          MPN_COPY(tp2,tp1,2*n);
          ref_redc_basecase(cp1,mp,n,inv,tp1);
          mpn_redc_basecase(cp2,mp,n,inv,tp2);
          if(mpn_cmp(cp1,cp2,n)!=0){printf("mpn_redc_basecase error %d\n",n);abort();}
          if(n!=1 && mpn_cmp(tp1,tp2,2*n)!=0){printf("mpn_redc_basecase possible error\n");abort();}
          /* we dont require the above to be the same but it could be a useful test */      
         }
     }



  tests_end ();
  exit (0);
}