c_cpp/etc/calc/custom/c_pmodm127.c

   1 /*
   2  * c_pmodm127 - calculate q mod 2^(2^127-1)
   3  *
   4  * Copyright (C) 2004-2007  Landon Curt Noll
   5  *
   6  * Calc is open software; you can redistribute it and/or modify it under
   7  * the terms of the version 2.1 of the GNU Lesser General Public License
   8  * as published by the Free Software Foundation.
   9  *
  10  * Calc is distributed in the hope that it will be useful, but WITHOUT
  11  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  12  * or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General
  13  * Public License for more details.
  14  *
  15  * A copy of version 2.1 of the GNU Lesser General Public License is
  16  * distributed with calc under the filename COPYING-LGPL.  You should have
  17  * received a copy with calc; if not, write to Free Software Foundation, Inc.
  18  * 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
  19  *
  20  * @(#) $Revision: 30.1 $
  21  * @(#) $Id: c_pmodm127.c,v 30.1 2007/03/16 11:10:04 chongo Exp $
  22  * @(#) $Source: /usr/local/src/bin/calc/custom/RCS/c_pmodm127.c,v $
  23  *
  24  * Under source code control:   2004/07/28 22:12:25
  25  * File existed as early as:    2004
  26  *
  27  * Share and enjoy!  :-)        http://www.isthe.com/chongo/tech/comp/calc/
  28  */
  29
  30
  31 #if defined(CUSTOM)
  32
  33 #include <stdio.h>
  34
  35 #include "have_const.h"
  36 #include "value.h"
  37 #include "custom.h"
  38 #include "zmath.h"
  39
  40 #include "have_unused.h"
  41
  42 /* 2^255 */
  43 STATIC HALF h255[] = {
  44 #if BASEB == 32
  45         (HALF)0x00000000, (HALF)0x00000000, (HALF)0x00000000, (HALF)0x00000000,
  46         (HALF)0x00000000, (HALF)0x00000000, (HALF)0x00000000, (HALF)0x80000000
  47 #else /* BASEB == 32 */
  48         (HALF)0x0000, (HALF)0x0000, (HALF)0x0000, (HALF)0x0000,
  49         (HALF)0x0000, (HALF)0x0000, (HALF)0x0000, (HALF)0x0000,
  50         (HALF)0x0000, (HALF)0x0000, (HALF)0x0000, (HALF)0x0000,
  51         (HALF)0x0000, (HALF)0x0000, (HALF)0x0000, (HALF)0x8000
  52 #endif /* BASEB == 32 */
  53 };
  54 ZVALUE p255 = {
  55     h255, 8, 0
  56 };
  57
  58
  59 /* static declarations */
  60 S_FUNC void zmod5_or_zmod(ZVALUE *zp);
  61 STATIC BOOL havelastmod = FALSE;
  62 STATIC ZVALUE lastmod[1];
  63 STATIC ZVALUE lastmodinv[1];
  64
  65
  66 /*
  67  * c_pmodm127 - calculate q mod 2^(2^127-1)
  68  *
  69  * given:
  70  *    count = 1;
  71  *    vals[0]   real number;    (q - potential factor)
  72  *
  73  * returns:
  74  *    result    real number;    (q mod 2^(2^127-1))
  75  */
  76 /*ARGSUSED*/
  77 VALUE
  78 c_pmodm127(char UNUSED *name, int UNUSED count, VALUE **vals)
  79 {
  80         VALUE result;           /* what we will return */
  81         ZVALUE q;               /* test factor */
  82         ZVALUE temp;            /* temp calculation value */
  83         int i;                  /* exponent value */
  84
  85         /*
  86          * arg check
  87          */
  88         result.v_type = V_NULL;
  89         if (vals[0]->v_type != V_NUM) {
  90                 math_error("Non-numeric argument for pmodm127");
  91                 /*NOTREACHED*/
  92         }
  93         if (qisfrac(vals[0]->v_num)) {
  94                 math_error("Non-integer argument for pmodm127");
  95                 /*NOTREACHED*/
  96         }
  97         if (qisneg(vals[0]->v_num) || qiszero(vals[0]->v_num)) {
  98                 math_error("argument for pmodm127 <= 0");
  99                 /*NOTREACHED*/
 100         }
 101
 102         /*
 103          * look at the numerator
 104          */
 105         q = vals[0]->v_num->num;
 106
 107         /*
 108          * setup lastmod with q
 109          */
 110         if (havelastmod && zcmp(q, *lastmod)) {
 111                 zfree(*lastmod);
 112                 zfree(*lastmodinv);
 113                 havelastmod = FALSE;
 114         }
 115         if (!havelastmod) {
 116                 zcopy(q, lastmod);
 117                 zbitvalue(2 * q.len * BASEB, &temp);
 118                 zquo(temp, q, lastmodinv, 0);
 119                 zfree(temp);
 120                 havelastmod = TRUE;
 121         }
 122
 123         /*
 124          * start with 2^255
 125          */
 126         result.v_num = qalloc();
 127         zcopy(p255, &result.v_num->num);
 128         result.v_type = V_NUM;
 129
 130         /*
 131          * compute 2^(2^127-1) mod q by modular exponentation
 132          *
 133          * We implement the following calc code in C:
 134          *
 135          *      (* given q, our test factor, as the arg to this function *)
 136          *      result = 2^255;
 137          *      for (i=8; i < 127; ++i) {
 138          *              result %= q;      (* mod *)
 139          *              result *= result; (* square *)
 140          *              result <<= 1;     (* times 2 *)
 141          *      }
 142          *      result %= q;              (* result is now 2^(2^127-1) % q *)
 143          */
 144         for (i=8; i<127; ++i) {
 145 #if 0
 146                 /* use of zmod is a bit slower than zmod5_or_zmod */
 147                 (void) zmod(result.v_num->num, *lastmod, &temp, 0);
 148                 zfree(result.v_num->num);
 149                 result.v_num->num = temp;
 150 #else
 151                 zmod5_or_zmod(&result.v_num->num);      /* mod */
 152 #endif
 153 #if 0
 154                 /* use of zmul is slightly slower than zsquare */
 155                 zmul(result.v_num->num, result.v_num->num, &temp);  /* square */
 156 #else
 157                 zsquare(result.v_num->num, &temp);      /* square */
 158 #endif
 159                 /*
 160                  * We could manually shift here, but this would o speed
 161                  * up the operation only a very tiny bit at the expense
 162                  * of a bunch of special code.
 163                  */
 164                 zfree(result.v_num->num);
 165                 zshift(temp, 1, &result.v_num->num);    /* times 2 */
 166                 zfree(temp);
 167         }
 168         zmod5_or_zmod(&result.v_num->num);      /* result = 2^(2^127-1) % q */
 169
 170         /*
 171          * cleanup and return result
 172          */
 173         return result;
 174 }
 175
 176
 177 /*
 178  * zmod5_or_zmod - fast integer modulo the modulus computation
 179  *
 180  * "borrowed" from ../zmod.c
 181  *
 182  * Given the address of a positive integer whose word count does not
 183  * exceed twice that of the modulus stored at lastmod, to evaluate and store
 184  * at that address the value of the integer modulo the modulus.
 185  *
 186  * Unlike the static routine in ../zmod.c, we will call the zmod and return
 187  * the result of the zmod5_or_zmod conditions do not apply to the argument
 188  * and saved mod.
 189  */
 190 S_FUNC void
 191 zmod5_or_zmod(ZVALUE *zp)
 192 {
 193         LEN len, modlen, j;
 194         ZVALUE tmp1, tmp2;
 195         ZVALUE z1, z2, z3;
 196         HALF *a, *b;
 197         FULL f;
 198         HALF u;
 199
 200         int subcount = 0;
 201
 202         if (zrel(*zp, *lastmod) < 0)
 203                 return;
 204         modlen = lastmod->len;
 205         len = zp->len;
 206         z1.v = zp->v + modlen - 1;
 207         z1.len = len - modlen + 1;
 208         z1.sign = z2.sign = z3.sign = 0;
 209         if (z1.len > modlen + 1) {
 210                 /* in ../zmod.c we did a math_error("Bad call to zmod5!!!"); */
 211                 /* here we just call the slower zmod and return the result */
 212                 (void) zmod(*zp, *lastmod, &tmp1, 0);
 213                 /* replace zp with tmp1 mod result */
 214                 zfree(*zp);
 215                 *zp = tmp1;
 216                 return;
 217         }
 218         z2.v = lastmodinv->v + modlen + 1 - z1.len;
 219         z2.len = lastmodinv->len - modlen - 1 + z1.len;
 220         zmul(z1, z2, &tmp1);
 221         z3.v = tmp1.v + z1.len;
 222         z3.len = tmp1.len - z1.len;
 223         if (z3.len > 0) {
 224                 zmul(z3, *lastmod, &tmp2);
 225                 j = modlen;
 226                 a = zp->v;
 227                 b = tmp2.v;
 228                 u = 0;
 229                 len = modlen;
 230                 while (j-- > 0) {
 231                         f = (FULL) *a - (FULL) *b++ - (FULL) u;
 232                         *a++ = (HALF) f;
 233                         u = - (HALF) (f >> BASEB);
 234                 }
 235                 if (z1.len > 1) {
 236                         len++;
 237                         if (tmp2.len > modlen)
 238                                 f = (FULL) *a - (FULL) *b - (FULL) u;
 239                         else
 240                                 f = (FULL) *a - (FULL) u;
 241                         *a++ = (HALF) f;
 242                 }
 243                 while (len > 0 && *--a == 0)
 244                         len--;
 245                 zp->len = len;
 246                 zfree(tmp2);
 247         }
 248         zfree(tmp1);
 249         while (len > 0 && zrel(*zp, *lastmod) >= 0) {
 250                 subcount++;
 251                 if (subcount > 2) {
 252                         math_error("Too many subtractions in zmod5_or_zmod");
 253                         /*NOTREACHED*/
 254                 }
 255                 j = modlen;
 256                 a = zp->v;
 257                 b = lastmod->v;
 258                 u = 0;
 259                 while (j-- > 0) {
 260                         f = (FULL) *a - (FULL) *b++ - (FULL) u;
 261                         *a++ = (HALF) f;
 262                         u = - (HALF) (f >> BASEB);
 263                 }
 264                 if (len > modlen) {
 265                         f = (FULL) *a - (FULL) u;
 266                         *a++ = (HALF) f;
 267                 }
 268                 while (len > 0 && *--a == 0)
 269                         len--;
 270                 zp->len = len;
 271         }
 272         if (len == 0)
 273                 zp->len = 1;
 274 }
 275
 276 #endif /* CUSTOM */