biginteger.c

   1 /* Copyright (c) 2008 Jérémie Laval, Marine Jourdain
   2
   3  * Permission is hereby granted, free of charge, to any person obtaining a copy
   4  * of this software and associated documentation files (the "Software"), to deal
   5  * in the Software without restriction, including without limitation the rights
   6  * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
   7  * copies of the Software, and to permit persons to whom the Software is
   8  * furnished to do so, subject to the following conditions:
   9
  10  * The above copyright notice and this permission notice shall be included in
  11  * all copies or substantial portions of the Software.
  12
  13  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
  14  * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
  15  * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
  16  * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
  17  * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
  18  * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
  19  * THE SOFTWARE.
  20  *
  21  */
  22
  23 #include <stdio.h>
  24 #include <stdlib.h>
  25 #include <math.h>
  26 #include <string.h>
  27
  28 #include "biginteger.h"
  29
  30 // Private Function Impl
  31
  32 #define ATOUI(str) ((unsigned int)strtoul(buffer, NULL, 10))
  33 #define MIN(x, y) (x < y ? x : y)
  34 #define MAX(x, y) (x > y ? x : y)
  35
  36 BigInteger CreateZeroBigInt()
  37 {
  38   BigInteger temp;
  39   temp.Sign = Undefined;
  40   temp.AbsValStart = NULL;
  41   temp.AbsValEnd = NULL;
  42   temp.Count = 0;
  43
  44   return temp;
  45 }
  46
  47 BigUnsignedInteger* InitBUint(unsigned int value)
  48 {
  49   BigUnsignedInteger* bigUint = malloc(sizeof(BigUnsignedInteger));
  50   bigUint->Value = value;
  51
  52   return bigUint;
  53 }
  54
  55 BigUnsignedInteger* CopyBigInteger(BigInteger source, BigInteger* dest, BigUnsignedInteger* start, int count)
  56 {
  57   BigUnsignedInteger* node = (start == NULL) ? source.AbsValStart : start;
  58
  59   do {
  60         PushEndNode(dest, node->Value);
  61         if (--count == 0)
  62           break;
  63   } while ((node = node->Next) != NULL);
  64
  65   return node;
  66 }
  67
  68 void Clean(BigInteger* b)
  69 {
  70   BigUnsignedInteger* node = b->AbsValEnd;
  71   if (node->Value != 0)
  72         return;
  73
  74   while (node->Value == 0) {
  75         BigUnsignedInteger* toFree = node;
  76         node = node->Previous;
  77         free(toFree);
  78   }
  79
  80   // Erase things up
  81   b->AbsValEnd = node;
  82   node->Next = NULL;
  83 }
  84
  85 // base = part1 * 10000^N  + part2
  86 void SplitBigInteger(int N, BigInteger base, BigInteger* part1, BigInteger* part2)
  87 {
  88   // In case N is odd part2 has a 1 more node than part1
  89   int count = (N % 2) == 0 ? N/2 : N/2 + 1;
  90
  91   if (base.Count <= count) {
  92         // Only set part2, part1 == 0
  93         /*part2->AbsValStart = base.AbsValStart;
  94         part2->AbsValEnd = base.AbsValEnd;
  95         part2->Count = base.Count;*/
  96         CopyBigInteger(base, part2, NULL, -1);
  97         part2->Sign = Positive;
  98         return;
  99   }
 100
 101   BigUnsignedInteger* node = CopyBigInteger(base, part2, NULL, count);
 102   part2->Sign = Positive;
 103
 104   CopyBigInteger(base, part1, node->Next, -1);
 105   part1->Sign = Positive;
 106
 107 }
 108
 109 void MultiplyByBase(BigInteger* b, int k)
 110 {
 111   int i = 0;
 112   for (; i < k; i++)
 113         PushFrontNode(b, 0);
 114 }
 115
 116 void PrependNode(BigUnsignedInteger* b1, BigUnsignedInteger* b2)
 117 {
 118   b1->Next = b2;
 119   b2->Previous = b1;
 120 }
 121
 122 // Used to push_front a value
 123 void PushFrontNode(BigInteger* bigInt, unsigned int value)
 124 {
 125   BigUnsignedInteger* bigUint = InitBUint(value);
 126
 127   // Do the magic swapping
 128   if (bigInt->AbsValStart != NULL)
 129         PrependNode(bigUint, bigInt->AbsValStart);
 130   // Update BigInt's head field
 131   bigInt->AbsValStart = bigUint;
 132
 133   // If it's the first node added, correctly set bigInt->AbsValEnd
 134   if (bigInt->AbsValEnd == NULL)
 135         bigInt->AbsValEnd = bigUint;
 136
 137   bigInt->Count++;
 138 }
 139
 140 void PushEndNode(BigInteger* bigInt, unsigned int value)
 141 {
 142   BigUnsignedInteger* buint = InitBUint(value);
 143
 144   // Do the magic swapping
 145   if (bigInt->AbsValEnd != NULL)
 146         PrependNode(bigInt->AbsValEnd, buint);
 147   // Update BigInt's head field
 148   bigInt->AbsValEnd = buint;
 149
 150   // If it's the first node added, correctly set bigInt->AbsValStart
 151   if (bigInt->AbsValStart == NULL)
 152         bigInt->AbsValStart = buint;
 153
 154   bigInt->Count++;
 155 }
 156
 157 // Used to insert a node in the *middle* of the BigInt.
 158 // It Do *not* check if we reproduce Push(Front|End)Node behavior
 159 // Use Push(Front|End)Node if you intend to do what these functions are meant for
 160 void InsertNode(BigInteger* bigInt, BigUnsignedInteger* next, unsigned int value)
 161 {
 162   BigUnsignedInteger* buint = InitBUint(value);
 163   BigUnsignedInteger* temp  = next->Previous;
 164
 165   PrependNode(buint, next);
 166   PrependNode(temp, buint);
 167
 168   bigInt->Count++;
 169 }
 170
 171 // End Private Function Impl
 172
 173
 174
 175 // Public Function Impl
 176
 177 BigInteger newBigInteger(char* desc)
 178 {
 179   int        sLength         = strlen(desc);
 180   BigInteger temp            = CreateZeroBigInt();
 181   int        startIndex = 0;
 182
 183   if (sLength == 0)
 184         return temp;
 185
 186   temp.Sign = Positive;
 187
 188   if (desc[0] == '-' || desc[0] == '+') {
 189         if (desc[0] == '-')
 190           temp.Sign = Negative;
 191         startIndex = 1;
 192   }
 193
 194   char buffer[5];
 195   memset(buffer, '\0', 5);
 196   int i;
 197
 198   for (i = sLength - 4;  i >= startIndex; i -= 4) {
 199         strncpy(buffer, desc + i, 4);
 200         PushEndNode(&temp, ATOUI(buffer));
 201         // Erase the buffer
 202         memset(buffer, '\0', 5);
 203   }
 204   // In case the MSB digits are not packable by four do the fiddling manually
 205   if (sLength % 4 != 0) {
 206         strncpy(buffer, desc + startIndex, i + 4 - startIndex);
 207         PushEndNode(&temp, ATOUI(buffer));
 208   }
 209
 210   return temp;
 211 }
 212
 213 void printBigInteger(BigInteger bint)
 214 {
 215   BigUnsignedInteger* node = bint.AbsValEnd;
 216   if (node == NULL) {
 217         puts("0");
 218         return;
 219   }
 220
 221   if (bint.Sign == Negative)
 222         printf("%c", '-');
 223
 224   do {
 225         if (node != bint.AbsValEnd)
 226           printf("%04u", node->Value);
 227         else
 228           printf("%u", node->Value);
 229   } while((node = node->Previous) != NULL);
 230
 231   printf("\n");
 232 }
 233
 234 Boolean isNull(BigInteger bint)
 235 {
 236   return bint.AbsValEnd == NULL ? true : false;
 237 }
 238
 239 int signBigInt(BigInteger bint)
 240 {
 241   return (int)bint.Sign;
 242 }
 243
 244 Boolean equalsBigInt(BigInteger b1, BigInteger b2)
 245 {
 246   return compareBigInt(b1, b2) == 0 ? true : false;
 247 }
 248
 249 /* b1 > b2 => 1
 250  * b1 < b2 => -1
 251  * b1 == b2 => 0
 252  */
 253 int compareBigInt(BigInteger b1, BigInteger b2)
 254 {
 255   // try to short-circuit
 256   if (b1.Count < b2.Count)
 257         return -1;
 258
 259   if (b1.Count > b2.Count)
 260         return 1;
 261
 262   // Normal processing
 263   BigUnsignedInteger* node1 = b1.AbsValEnd;
 264   BigUnsignedInteger* node2 = b2.AbsValEnd;
 265
 266   // Test for NULL special case
 267   if (node1 == NULL && node2 == NULL)
 268         return 0;
 269
 270   if (b1.Sign != b2.Sign) {
 271         return b1.Sign == Negative ? -1 : 1;
 272   }
 273
 274   do {
 275         if (node1->Value < node2->Value)
 276           return -1;
 277         if (node1->Value > node2->Value)
 278           return 1;
 279   } while ((node1 = node1->Previous) != NULL && (node2 = node2->Previous) != NULL);
 280
 281   return 0;
 282 }
 283
 284 BigInteger sumBigInt(BigInteger b1, BigInteger b2)
 285 {
 286   // Check some sign conditions
 287   // (-b1 + b2) <=> (b2 - (+b2))
 288   if (b1.Sign == Negative && b2.Sign == Positive) {
 289         b1.Sign = Positive;
 290         return diffBigInt(b2, b1);
 291   }
 292   // (b1 + (-b2)) <=> (b1 - (+b2))
 293   if (b1.Sign == Positive && b2.Sign == Negative) {
 294         b2.Sign = Positive;
 295         return diffBigInt(b1, b2);
 296   }
 297
 298   BigInteger result = CreateZeroBigInt();
 299
 300   if (b1.Sign == Negative && b2.Sign == Negative)
 301         result.Sign = Negative;
 302   else
 303         result.Sign = Positive;
 304
 305   BigUnsignedInteger* node1 = b1.AbsValStart;
 306   BigUnsignedInteger* node2 = b2.AbsValStart;
 307
 308   // NULL case
 309   if (node1 == NULL)
 310         return b2;
 311   if (node2 == NULL)
 312         return b1;
 313
 314   unsigned int carry = 0;
 315   Boolean cont = false;
 316
 317   do {
 318         cont = false;
 319
 320         unsigned int newValue = (node1 == NULL ? 0 : node1->Value) + (node2 == NULL ? 0 : node2->Value) + carry;
 321
 322         carry = newValue / BASE;
 323         newValue = newValue % BASE;
 324         PushEndNode(&result, newValue);
 325
 326         if (node1 != NULL) {
 327           cont = (node1 = node1->Next) != NULL;
 328         }
 329         if (node2 != NULL) {
 330           node2 = node2->Next;
 331           cont = node2 != NULL || cont;
 332         }
 333
 334   } while(cont);
 335
 336   if (carry > 0)
 337         PushEndNode(&result, carry);
 338
 339   return result;
 340 }
 341
 342 BigInteger diffBigInt(BigInteger b1, BigInteger b2)
 343 {
 344   // Sign cases
 345   // -b1 - (+b2) <=> -(b1 + b2)
 346   if (b1.Sign == Negative && b2.Sign == Positive) {
 347         b2.Sign = Negative;
 348         return sumBigInt(b1, b2);
 349   }
 350
 351   // b1 - (-b2) <=> b1 + b2
 352   if (b1.Sign == Positive && b2.Sign == Negative) {
 353         b2.Sign = Positive;
 354         return sumBigInt(b1, b2);
 355   }
 356
 357   BigInteger result = CreateZeroBigInt();
 358
 359   BigUnsignedInteger* node1 = b1.AbsValStart;
 360   BigUnsignedInteger* node2 = b2.AbsValStart;
 361
 362   // NULL case
 363   if (node2 == NULL)
 364         return b1;
 365   if (node1 == NULL) {
 366         b2.Sign = b2.Sign == Positive ? Negative : Positive;
 367         return b2;
 368   }
 369
 370   int comp = compareBigInt(b1, b2);
 371   if (comp == 0) {
 372         return result;
 373   } else if (comp > 0) {
 374         result.Sign = Positive;
 375   } else {
 376         // Some magic swapping
 377         result.Sign = Negative;
 378         BigUnsignedInteger* temp = node1;
 379         node1 = node2;
 380         node2 = temp;
 381   }
 382
 383   // Takes back the pointer value as there might have been swapping
 384   // TODO : see if we can exchange directly these pointer in the if rather than b1/b2
 385   /*node1 = b1.AbsValStart;
 386   node2 = b2.AbsValStart;
 387   */
 388
 389   unsigned int carry = 0;
 390   Boolean cont = false;
 391
 392   do {
 393         cont = false;
 394
 395         unsigned int newValue = 0;
 396         if (node1 != NULL && node2 != NULL) {
 397           if (node2->Value > node1->Value) {
 398                 // Do a little trick to get positive coefficient
 399                 // and backport the changement on the upper coefficient with carry
 400                 newValue = (BASE + node1->Value) - node2->Value - carry;
 401                 carry = 1;
 402           } else {
 403                 newValue = node1->Value - node2->Value - carry;
 404                 carry = 0;
 405           }
 406         } else {
 407           // Due to our hypothethis at start of the function (test, swapping and all)
 408           // we know that it's node2 which is NULL, then proceed with node1
 409           newValue = node1->Value - carry;
 410           carry = 0;
 411           // No need for unnecessary 0 in the list
 412           if (newValue == 0)
 413                 break;
 414         }
 415
 416         PushEndNode(&result, newValue);
 417
 418         if (node1 != NULL) {
 419           cont = (node1 = node1->Next) != NULL;
 420         }
 421         if (node2 != NULL) {
 422           node2 = node2->Next;
 423           cont = node2 != NULL || cont;
 424         }
 425
 426   } while(cont);
 427
 428   return result;
 429 }
 430
 431 // Uses Karatsuba algorithm
 432 BigInteger mulBigInt(BigInteger x, BigInteger y)
 433 {
 434   if (x.Sign == 0 || y.Sign == 0)
 435         return CreateZeroBigInt();
 436
 437   // If we have a biginteger with one node no point in further splitting
 438   if (MIN(x.Count, y.Count) == 1) {
 439         BigInteger temp = CreateZeroBigInt();
 440         temp.Sign = Positive;
 441
 442         BigUnsignedInteger* big = (x.Count == 1) ? y.AbsValStart : x.AbsValStart;
 443         unsigned int multiplier = (x.Count == 1) ? x.AbsValStart->Value : y.AbsValStart->Value;
 444
 445         unsigned int carry = 0;
 446
 447         while (big != NULL) {
 448           unsigned int r = big->Value * multiplier + carry;
 449           carry = r / BASE;
 450           PushEndNode(&temp, r % BASE);
 451
 452           big = big->Next;
 453         }
 454         if (carry != 0)
 455           PushEndNode(&temp, carry);
 456
 457         return temp;
 458   }
 459
 460   BigInteger x0 = CreateZeroBigInt();
 461   BigInteger x1 = CreateZeroBigInt();
 462
 463   BigInteger y0 = CreateZeroBigInt();
 464   BigInteger y1 = CreateZeroBigInt();
 465
 466   int N = MAX(x.Count, y.Count);
 467   SplitBigInteger(N, x, &x1, &x0);
 468   SplitBigInteger(N, y, &y1, &y0);
 469
 470   /* a = x1*y1
 471    * b = x0*y0
 472    * c = (x1 - x0)(y1 - y0)
 473    */
 474   BigInteger a = mulBigInt(x1, y1);
 475   BigInteger b = mulBigInt(x0, y0);
 476   BigInteger c = mulBigInt(diffBigInt(x1, x0), diffBigInt(y1, y0));
 477   BigInteger abc = sumBigInt(a, b);
 478   abc = diffBigInt(abc, c);
 479
 480   int k = (N % 2) == 0 ? N/2 : N/2 + 1;
 481   MultiplyByBase(&a, 2*k);
 482   MultiplyByBase(&abc, k);
 483
 484   BigInteger result = sumBigInt(sumBigInt(a, abc), b);
 485   // Remove leading 0s if any
 486   Clean(&result);
 487
 488   return result;
 489 }
 490
 491 // Returns b^n
 492 BigInteger powBigInt(BigInteger b, int n)
 493 {
 494   BigInteger result = b;
 495
 496   while ((n & 0x1) == 0) {
 497         result = mulBigInt(result, result);
 498         n /= 2;
 499   }
 500   if (n != 1) {
 501         int i = n;
 502         for (; i >= 1; i--)
 503           result = mulBigInt(result, b);
 504   }
 505
 506   return result;
 507 }
 508
 509 // End Public Function Impl