/[dtapublic]/projs/trunk/shared_source/c_datd/gmp_ints.c
ViewVC logotype

Annotation of /projs/trunk/shared_source/c_datd/gmp_ints.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 56 - (hide annotations) (download)
Sat Oct 29 01:53:01 2016 UTC (7 years, 8 months ago) by dashley
File MIME type: text/plain
File size: 108843 byte(s)
License and property (keyword) changes.
1 dashley 56 //$Header$
2 dashley 25 //-------------------------------------------------------------------------------------------------
3 dashley 56 //This file is part of "David T. Ashley's Shared Source Code", a set of shared components
4     //integrated into many of David T. Ashley's projects.
5 dashley 25 //-------------------------------------------------------------------------------------------------
6 dashley 56 //This source code and any program in which it is compiled/used is provided under the MIT License,
7     //reproduced below.
8     //-------------------------------------------------------------------------------------------------
9     //Permission is hereby granted, free of charge, to any person obtaining a copy of
10     //this software and associated documentation files(the "Software"), to deal in the
11     //Software without restriction, including without limitation the rights to use,
12     //copy, modify, merge, publish, distribute, sublicense, and / or sell copies of the
13     //Software, and to permit persons to whom the Software is furnished to do so,
14     //subject to the following conditions :
15 dashley 25 //
16 dashley 56 //The above copyright notice and this permission notice shall be included in all
17     //copies or substantial portions of the Software.
18 dashley 25 //
19 dashley 56 //THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
20     //IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
21     //FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.IN NO EVENT SHALL THE
22     //AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
23     //LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
24     //OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
25     //SOFTWARE.
26 dashley 25 //-------------------------------------------------------------------------------------------------
27     #define MODULE_GMP_INTS
28    
29     #include <assert.h>
30     #include <ctype.h>
31     #include <process.h>
32     #include <stddef.h>
33     #include <stdio.h>
34     #include <string.h>
35     #include <time.h>
36    
37    
38     /* Only included the guarded allocation header if we are compiling
39     ** a DOS console type application. Other types of applications have
40     ** other ways of protecting for out of memory. Including the
41     ** header would do no harm in these cases, but do no good, either.
42     */
43     #if defined(APP_TYPE_SIMPLE_DOS_CONSOLE)
44     #include "ccmalloc.h"
45     #elif defined(APP_TYPE_IJUSCRIPTER_IJUCONSOLE)
46     #include "tclalloc.h"
47     #else
48     /* Do nothing. */
49     #endif
50    
51     #include "bstrfunc.h"
52     #include "charfunc.h"
53     #include "fcmiof.h"
54     #include "gmp_ints.h"
55     #include "intfunc.h"
56    
57    
58     /******************************************************************/
59     /*** CUSTOM ALLOCATION FUNCTIONS *******************************/
60     /******************************************************************/
61     /* We need to decide here on how memory not on the stack will be
62     ** allocated (i.e. what will become of the standard functions
63     ** like malloc, free, etc.). This is dependent on the type of
64     ** application we're making. The possible types are so far are:
65     ** APP_TYPE_SIMPLE_DOS_CONSOLE :
66     ** Simple DOS console application.
67     ** APP_TYPE_IJUSCRIPTER_IJUCONSOLE:
68     ** The Tcl tool build.
69     **
70     ** The custom allocation functions here are a "portal" or "wrapper"
71     ** for how the integer and rational number functions should
72     ** get memory.
73     **
74     ** The functions below are standard, except that the GNU MP team
75     ** built more generality into what allocation schemes could be
76     ** used by including size information in some calls that don't
77     ** normally get it. That is why there are some extra calls below
78     ** where the information is discarded. Other than that, these are
79     ** standard allocation calls.
80     */
81     //07/15/01: Visual inspection only. Function deemed too
82     //simple for unit testing.
83     void *GMP_INTS_malloc( size_t size )
84     {
85     #if defined(APP_TYPE_SIMPLE_DOS_CONSOLE)
86     return(CCMALLOC_malloc(size));
87     #elif defined(APP_TYPE_IJUSCRIPTER_IJUCONSOLE)
88     return(TclpAlloc(size));
89     #else
90     return(malloc(size));
91     #endif
92     }
93    
94    
95     //07/15/01: Visual inspection only. Function deemed too
96     //simple for unit testing.
97     void *GMP_INTS_calloc( size_t num, size_t size )
98     {
99     #if defined(APP_TYPE_SIMPLE_DOS_CONSOLE)
100     return(CCMALLOC_calloc(num, size));
101     #elif defined(APP_TYPE_IJUSCRIPTER_IJUCONSOLE)
102     return(TclpCalloc(num, size));
103     #else
104     return(calloc(num, size));
105     #endif
106     }
107    
108    
109     //07/15/01: Visual inspection only. Function deemed too
110     //simple for unit testing.
111     void *GMP_INTS_realloc( void *memblock, size_t size )
112     {
113     #if defined(APP_TYPE_SIMPLE_DOS_CONSOLE)
114     return(CCMALLOC_realloc(memblock, size));
115     #elif defined(APP_TYPE_IJUSCRIPTER_IJUCONSOLE)
116     return(TclpRealloc(memblock, size));
117     #else
118     return(realloc(memblock, size));
119     #endif
120     }
121    
122    
123     //07/15/01: Visual inspection only. Function deemed too
124     //simple for unit testing.
125     void *GMP_INTS_realloc_w_size( void *memblock,
126     size_t old_size,
127     size_t size )
128     {
129     #if defined(APP_TYPE_SIMPLE_DOS_CONSOLE)
130     return(CCMALLOC_realloc(memblock, size));
131     #elif defined(APP_TYPE_IJUSCRIPTER_IJUCONSOLE)
132     return(TclpRealloc(memblock, size));
133     #else
134     return(realloc(memblock, size));
135     #endif
136     }
137    
138    
139     //07/15/01: Visual inspection only. Function deemed too
140     //simple for unit testing.
141     void GMP_INTS_free( void *memblock )
142     {
143     #if defined(APP_TYPE_SIMPLE_DOS_CONSOLE)
144     CCMALLOC_free(memblock);
145     #elif defined(APP_TYPE_IJUSCRIPTER_IJUCONSOLE)
146     TclpFree(memblock);
147     #else
148     free(memblock);
149     #endif
150     }
151    
152    
153     //07/15/01: Visual inspection only. Function deemed too
154     //simple for unit testing.
155     void GMP_INTS_free_w_size( void *memblock, size_t size )
156     {
157     #if defined(APP_TYPE_SIMPLE_DOS_CONSOLE)
158     CCMALLOC_free(memblock);
159     #elif defined(APP_TYPE_IJUSCRIPTER_IJUCONSOLE)
160     TclpFree(memblock);
161     #else
162     free(memblock);
163     #endif
164     }
165    
166    
167     /******************************************************************/
168     /*** PORTABILITY CHECK FUNCTIONS *******************************/
169     /******************************************************************/
170     //Because there is the risk that Microsoft Visual C++ might
171     //change in the future, the following function can be called
172     //to see if the assumptions about data sizes are valid. This
173     //function returns TRUE if there is a problem, or FALSE
174     //otherwise.
175    
176     //07/15/01: Unit testing complete.
177     int GMP_INTS_data_sizes_are_wrong(void)
178     {
179     int i;
180     GMP_INTS_limb_t tv;
181     _int64 tv64;
182    
183     //Check the number of bit rolls required to get the limb
184     //to go to zero again. This had better be 32.
185     tv = 1;
186     i = 0;
187     while (tv)
188     {
189     tv <<= 1;
190     i++;
191     }
192     if (i != 32)
193     return(1);
194    
195     //Check that an _int64 is really and truly 64 bits.
196     tv64 = 1;
197     i = 0;
198     while (tv64)
199     {
200     tv64 <<= 1;
201     i++;
202     }
203     if (i != 64)
204     return(1);
205    
206     //Room for additional tests here if needed later.
207    
208     return(0);
209     }
210    
211    
212     /******************************************************************/
213     /*** ERROR STRING IDENTIFICATION AND PROCESSING FUNCTIONS *******/
214     /******************************************************************/
215    
216     int GMP_INTS_identify_nan_string(const char *s)
217     {
218     assert(s != NULL);
219    
220     if (!strcmp(s, GMP_INTS_EF_INTOVF_POS_STRING))
221     return(0);
222     else if (!strcmp(s, GMP_INTS_EF_INTOVF_NEG_STRING))
223     return(1);
224     else if (!strcmp(s, GMP_INTS_EF_INTOVF_TAINT_POS_STRING))
225     return(2);
226     else if (!strcmp(s, GMP_INTS_EF_INTOVF_TAINT_NEG_STRING))
227     return(3);
228     else
229     return(-1);
230     }
231    
232    
233     const char *GMP_INTS_supply_nan_string(int idx)
234     {
235     assert((idx >= 0) && (idx <= 3));
236    
237     if (idx==0)
238     return(GMP_INTS_EF_INTOVF_POS_STRING);
239     else if (idx==1)
240     return(GMP_INTS_EF_INTOVF_NEG_STRING);
241     else if (idx==2)
242     return(GMP_INTS_EF_INTOVF_TAINT_POS_STRING);
243     else
244     return(GMP_INTS_EF_INTOVF_TAINT_NEG_STRING);
245     }
246    
247    
248     /******************************************************************/
249     /*** DEBUG PRINTING FUNCTIONS **********************************/
250     /******************************************************************/
251     //These functions are for printing out integers and limbs
252     //and groups of limbs for unit testing and debugging.
253    
254     //07/15/01: Exempt from testing because debug/development
255     //function.
256     void GMP_INTS_print_limb_group(FILE *stream,
257     GMP_INTS_limb_srcptr lg,
258     GMP_INTS_size_t n,
259     char *desc)
260     {
261     int i;
262    
263     assert(stream != NULL);
264     assert(n >= 0);
265     assert(desc != NULL);
266    
267     if (!lg)
268     {
269     fprintf(stream, " %s: NULL\n", desc);
270     }
271     else
272     {
273     for (i=n-1; i>=0; i--)
274     {
275     fprintf(stream, " %s[%2d]: 0x%8X\n", desc, i, lg[i]);
276     }
277     }
278     }
279    
280    
281     void GMP_INTS_mpz_print_int(FILE *stream,
282     const GMP_INTS_mpz_struct *arg,
283     char *desc)
284     {
285     int i;
286    
287     assert(stream != NULL);
288     assert(arg != NULL);
289     assert(desc != NULL);
290    
291     fprintf(stream, "Printing integer:\n %s\n", desc);
292    
293     fprintf(stream, " flags: %d\n", arg->flags);
294     fprintf(stream, " ptr value to body: %p\n", arg);
295     fprintf(stream, " n_allocd: %d\n", arg->n_allocd);
296     fprintf(stream, " size: %d\n", arg->size);
297     fprintf(stream, " limbs (ptr val): %p\n", arg->limbs);
298    
299     for (i=arg->n_allocd-1; i>=0; i--)
300     {
301     fprintf(stream, " limb[%3d]: %8X\n", i, arg->limbs[i]);
302     }
303     }
304    
305    
306     /******************************************************************/
307     /*** LOW-LEVEL MACRO REPLACEMENTS ******************************/
308     /******************************************************************/
309     //The functions in this category are replacements for macros.
310     //Clarity was gained at the expense of speed.
311    
312     int GMP_INTS_mpz_get_flags (const GMP_INTS_mpz_struct *arg)
313     {
314     assert(arg != NULL);
315     assert(arg->n_allocd > 0);
316    
317     return(arg->flags);
318     }
319    
320    
321     //07/15/01: Visual inspection only. Function deemed too
322     //simple for unit testing.
323     GMP_INTS_size_t GMP_INTS_abs_of_size_t(GMP_INTS_size_t arg)
324     {
325     //Be sure that the bit pattern does not represent the maximum
326     //negative argument. Negating this would give the result of
327     //zero, which is not what is intended.
328     assert(arg != 0x80000000);
329    
330     if (arg < 0)
331     return(-arg);
332     else
333     return(arg);
334     }
335    
336    
337     //07/15/01: Visual inspection only. Function deemed too
338     //simple for unit testing.
339     int GMP_INTS_mpz_sgn(const GMP_INTS_mpz_struct *arg)
340     {
341     assert(arg != NULL);
342     assert(arg->n_allocd > 0);
343    
344     if (arg->size > 0)
345     return(1);
346     else if (arg->size == 0)
347     return(0);
348     else
349     return(-1);
350     }
351    
352    
353     //07/15/01: Visual inspection only. Function deemed too
354     //simple for unit testing.
355     int GMP_INTS_mpz_is_neg(const GMP_INTS_mpz_struct *arg)
356     {
357     assert(arg != NULL);
358     assert(arg->n_allocd > 0);
359    
360     if (GMP_INTS_mpz_sgn(arg) == -1)
361     return(1);
362     else
363     return(0);
364     }
365    
366    
367     //07/15/01: Visual inspection only. Function deemed too
368     //simple for unit testing.
369     int GMP_INTS_mpz_is_zero(const GMP_INTS_mpz_struct *arg)
370     {
371     assert(arg != NULL);
372     assert(arg->n_allocd > 0);
373    
374     if (GMP_INTS_mpz_sgn(arg) == 0)
375     return(1);
376     else
377     return(0);
378     }
379    
380    
381     //07/15/01: Visual inspection only. Function deemed too
382     //simple for unit testing.
383     int GMP_INTS_mpz_is_pos(const GMP_INTS_mpz_struct *arg)
384     {
385     assert(arg != NULL);
386     assert(arg->n_allocd > 0);
387    
388     if (GMP_INTS_mpz_sgn(arg) == 1)
389     return(1);
390     else
391     return(0);
392     }
393    
394    
395     //07/15/01: Visual inspection only. Function deemed too
396     //simple for unit testing.
397     int GMP_INTS_mpz_is_odd(const GMP_INTS_mpz_struct *arg)
398     {
399     assert(arg != NULL);
400     assert(arg->n_allocd > 0);
401    
402     if (arg->size == 0)
403     return 0;
404     else if ((arg->limbs[0] & 0x1) != 0)
405     return 1;
406     else
407     return 0;
408     }
409    
410    
411     //07/15/01: Visual inspection only. Function deemed too
412     //simple for unit testing.
413     int GMP_INTS_mpz_is_even(const GMP_INTS_mpz_struct *arg)
414     {
415     assert(arg != NULL);
416     assert(arg->n_allocd > 0);
417    
418     if (GMP_INTS_mpz_is_odd(arg))
419     return 0;
420     else
421     return 1;
422     }
423    
424    
425     void GMP_INTS_mpz_negate(GMP_INTS_mpz_struct *arg)
426     {
427     //Eyeball the input parameters.
428     assert(arg != NULL);
429     assert(arg->n_allocd > 0);
430     assert(arg->limbs != NULL);
431    
432     arg->size = -(arg->size);
433     }
434    
435    
436     //07/15/01: Visual inspection only. Function deemed too
437     //simple for unit testing.
438     void GMP_INTS_mpn_normalize(GMP_INTS_limb_ptr limb_array,
439     GMP_INTS_size_t *idx)
440     {
441     assert(limb_array != NULL);
442     assert(idx != NULL);
443     assert(idx >= 0);
444    
445     while (*idx > 0)
446     {
447     if (limb_array[*idx - 1] != 0)
448     break;
449     (*idx)--;
450     }
451     }
452    
453    
454     //07/15/01: Visual inspection only. Function deemed too
455     //simple for unit testing.
456     void GMP_INTS_mpn_copy_limbs(GMP_INTS_limb_ptr dest,
457     GMP_INTS_limb_srcptr src,
458     GMP_INTS_size_t n)
459     {
460     GMP_INTS_size_t i;
461    
462     assert(dest != NULL);
463     assert(src != NULL);
464     assert(n >= 0);
465    
466     for (i=0; i<n; i++)
467     dest[i] = src[i];
468     }
469    
470    
471     /******************************************************************/
472     /*** LOW-LEVEL ARITHMETIC FUNCTIONS ****************************/
473     /******************************************************************/
474    
475     int GMP_INTS_two_op_flags_map(int flags1, int flags2)
476     {
477     int cf;
478    
479     if (!flags1 && !flags2)
480     {
481     return(0);
482     }
483     else
484     {
485     cf = flags1 | flags2;
486    
487     if (cf & (GMP_INTS_EF_INTOVF_POS | GMP_INTS_EF_INTOVF_TAINT_POS))
488     {
489     //In either case here, the result will be tainted
490     //positive.
491     return(GMP_INTS_EF_INTOVF_TAINT_POS);
492     }
493     else if (cf & (GMP_INTS_EF_INTOVF_NEG | GMP_INTS_EF_INTOVF_TAINT_NEG))
494     {
495     //In either case here, the result will be tainted
496     //negative.
497     return(GMP_INTS_EF_INTOVF_TAINT_NEG);
498     }
499     else
500     {
501     //This case is where the flags ints are non-zero, but
502     //no known bits are set. This is surely some kind of
503     //an internal software error. In debug mode, want to
504     //alert to error. In actual operation, want to just
505     //pretend an ordinary positive taint.
506     assert(0);
507     return(GMP_INTS_EF_INTOVF_TAINT_POS);
508     }
509     }
510     }
511    
512    
513     GMP_INTS_limb_t GMP_INTS_mpn_add_1 (GMP_INTS_limb_ptr res_ptr,
514     GMP_INTS_limb_srcptr s1_ptr,
515     GMP_INTS_size_t s1_size,
516     GMP_INTS_limb_t s2_limb)
517     {
518     GMP_INTS_limb_t x;
519    
520     assert(res_ptr != NULL);
521     assert(s1_ptr != NULL);
522     assert(s1_size > 0);
523    
524     x = *s1_ptr++;
525     s2_limb = x + s2_limb;
526     *res_ptr++ = s2_limb;
527     //Since limbs are unsigned, the test below tests if there
528     //was a carry, i.e. a positive rollover.
529     if (s2_limb < x)
530     {
531     while (--s1_size != 0)
532     {
533     x = *s1_ptr++ + 1;
534     *res_ptr++ = x;
535     if (x != 0)
536     goto fin;
537     }
538    
539     return 1;
540     }
541    
542     fin:
543     if (res_ptr != s1_ptr)
544     {
545     GMP_INTS_size_t i;
546     for (i = 0; i < s1_size - 1; i++)
547     {
548     res_ptr[i] = s1_ptr[i];
549     }
550     }
551    
552     return 0;
553     }
554    
555    
556     GMP_INTS_limb_t GMP_INTS_mpn_sub_1(GMP_INTS_limb_ptr res_ptr,
557     GMP_INTS_limb_srcptr s1_ptr,
558     GMP_INTS_size_t s1_size,
559     GMP_INTS_limb_t s2_limb)
560     {
561     GMP_INTS_limb_t x;
562    
563     assert(res_ptr != NULL);
564     assert(s1_ptr != NULL);
565     assert(s1_size > 0);
566    
567     x = *s1_ptr++;
568     s2_limb = x - s2_limb;
569     *res_ptr++ = s2_limb;
570     //The test below detects a borrow.
571     if (s2_limb > x)
572     {
573     while (--s1_size != 0)
574     {
575     x = *s1_ptr++;
576     *res_ptr++ = x - 1;
577     if (x != 0)
578     goto fin;
579     }
580    
581     return 1;
582     }
583    
584     fin:
585     if (res_ptr != s1_ptr)
586     {
587     GMP_INTS_size_t i;
588     for (i = 0; i < s1_size - 1; i++)
589     {
590     res_ptr[i] = s1_ptr[i];
591     }
592     }
593     return 0;
594     }
595    
596    
597     //07/15/01: Am willing to skip unit-testing on this.
598     //Understand the logic (i.e. passes visual inspection),
599     //and comes from GNU-MP. Hope any defects here will be
600     //caught in testing of GMP_INTS_mpz_mul() and other
601     //higher-level functions.
602     GMP_INTS_limb_t GMP_INTS_mpn_mul_1 (GMP_INTS_limb_ptr res_ptr,
603     GMP_INTS_limb_srcptr s1_ptr,
604     GMP_INTS_size_t s1_size,
605     GMP_INTS_limb_t s2_limb)
606     {
607     GMP_INTS_limb_t cy_limb;
608     GMP_INTS_size_t j;
609     GMP_INTS_limb_t prod_high, prod_low;
610     unsigned _int64 temp;
611    
612     assert(res_ptr != NULL);
613     assert(s1_ptr != NULL);
614     assert(s1_size > 0);
615    
616     /* The loop counter and index J goes from -S1_SIZE to -1. This way
617     the loop becomes faster. */
618     j = -s1_size;
619    
620     /* Offset the base pointers to compensate for the negative indices. */
621     s1_ptr -= j;
622     res_ptr -= j;
623    
624     cy_limb = 0;
625     do
626     {
627     //The original code here was the following macro:
628     //umul_ppmm (prod_high, prod_low, s1_ptr[j], s2_limb);
629     //Will use the 64-bit data type of MSVC++ to achieve
630     //the same effect.
631     //
632     //NOTE AS OF 07/13/01: I have looked at the assembly-
633     //language, and the lines below are a real sore spot.
634     //The multiply is fairly direct (although there is a
635     //function call), but the shift does not behave as
636     //expected--there is a function call and a loop to
637     //go through the 32 iterations. After logical testing,
638     //may want to clean this out--this would surely
639     //result in a speed increase. This is a sore spot.
640     temp = ((unsigned _int64)s1_ptr[j]) * ((unsigned _int64)s2_limb);
641     prod_low = (GMP_INTS_limb_t)temp;
642     prod_high = (GMP_INTS_limb_t)(temp >> 32);
643    
644     prod_low += cy_limb;
645     cy_limb = (prod_low < cy_limb) + prod_high;
646    
647     res_ptr[j] = prod_low;
648     }
649     while (++j != 0);
650    
651     return cy_limb;
652     }
653    
654    
655     //07/15/01: Am willing to skip unit-testing on this.
656     //Understand the logic (i.e. passes visual inspection),
657     //and comes from GNU-MP. Hope any defects here will be
658     //caught in testing of GMP_INTS_mpz_add() and other
659     //higher-level functions.
660     GMP_INTS_limb_t GMP_INTS_mpn_add_n(GMP_INTS_limb_ptr res_ptr,
661     GMP_INTS_limb_srcptr s1_ptr,
662     GMP_INTS_limb_srcptr s2_ptr,
663     GMP_INTS_size_t size)
664     {
665     GMP_INTS_limb_t x, y, cy;
666     GMP_INTS_size_t j;
667    
668     assert(res_ptr != NULL);
669     assert(s1_ptr != NULL);
670     assert(s2_ptr != NULL);
671    
672     /* The loop counter and index J goes from -SIZE to -1. This way
673     the loop becomes faster. */
674     j = -size;
675    
676     /* Offset the base pointers to compensate for the negative indices. */
677     s1_ptr -= j;
678     s2_ptr -= j;
679     res_ptr -= j;
680    
681     cy = 0;
682     do
683     {
684     y = s2_ptr[j];
685     x = s1_ptr[j];
686     y += cy; /* add previous carry to one addend */
687     cy = (y < cy); /* get out carry from that addition */
688     y = x + y; /* add other addend */
689     cy = (y < x) + cy; /* get out carry from that add, combine */
690     res_ptr[j] = y;
691     }
692     while (++j != 0);
693    
694     return cy;
695     }
696    
697    
698     //07/15/01: Am willing to skip unit-testing on this.
699     //Understand the logic (i.e. passes visual inspection),
700     //and comes from GNU-MP. Hope any defects here will be
701     //caught in testing of GMP_INTS_mpz_mul() and other
702     //higher-level functions.
703     GMP_INTS_limb_t GMP_INTS_mpn_addmul_1 (GMP_INTS_limb_ptr res_ptr,
704     GMP_INTS_limb_srcptr s1_ptr,
705     GMP_INTS_size_t s1_size,
706     GMP_INTS_limb_t s2_limb)
707     {
708     GMP_INTS_limb_t cy_limb;
709     GMP_INTS_size_t j;
710     GMP_INTS_limb_t prod_high, prod_low;
711     GMP_INTS_limb_t x;
712     unsigned _int64 temp;
713    
714     //Eyeball the inputs carefully.
715     assert(res_ptr != NULL);
716     assert(s1_ptr != NULL);
717     assert(s1_size > 0);
718    
719     /* The loop counter and index J goes from -SIZE to -1. This way
720     the loop becomes faster. */
721     j = -s1_size;
722    
723     /* Offset the base pointers to compensate for the negative indices. */
724     res_ptr -= j;
725     s1_ptr -= j;
726    
727     cy_limb = 0;
728     do
729     {
730     //The original code here was the following macro:
731     //umul_ppmm (prod_high, prod_low, s1_ptr[j], s2_limb);
732     //Will use the 64-bit data type of MSVC++ to achieve
733     //the same effect.
734     //
735     //NOTE AS OF 07/14/01: I have not looked at the assembly-
736     //language, but the assembly-language generated by what
737     //is below is suspected to have performance problems.
738     //May want to come back to this.
739     temp = ((unsigned _int64)s1_ptr[j]) * ((unsigned _int64)s2_limb);
740     prod_low = (GMP_INTS_limb_t)temp;
741     prod_high = (GMP_INTS_limb_t)(temp >> 32);
742    
743     prod_low += cy_limb;
744     cy_limb = (prod_low < cy_limb) + prod_high;
745    
746     x = res_ptr[j];
747     prod_low = x + prod_low;
748     cy_limb += (prod_low < x);
749     res_ptr[j] = prod_low;
750     }
751     while (++j != 0);
752    
753     return cy_limb;
754     }
755    
756    
757     //07/15/01: Am willing to skip unit-testing on this.
758     //Understand the logic (i.e. passes visual inspection),
759     //and comes from GNU-MP.
760     GMP_INTS_limb_t GMP_INTS_mpn_add (GMP_INTS_limb_ptr res_ptr,
761     GMP_INTS_limb_srcptr s1_ptr,
762     GMP_INTS_size_t s1_size,
763     GMP_INTS_limb_srcptr s2_ptr,
764     GMP_INTS_size_t s2_size)
765     {
766     GMP_INTS_limb_t cy_limb = 0;
767    
768     assert(res_ptr != NULL);
769     assert(s1_ptr != NULL);
770     assert(s2_ptr != NULL);
771    
772     //Numbers apparently must be arranged with sizes so that
773     //LIMBS(s1) >= LIMBS(s2).
774     //Add the parts up to the most significant limb of S2.
775     if (s2_size != 0)
776     cy_limb = GMP_INTS_mpn_add_n (res_ptr,
777     s1_ptr,
778     s2_ptr,
779     s2_size);
780    
781     //Process the carry result, and propagate the carries up through
782     //the parts of S1 that don't exist in S2, i.e. propagate the
783     //carries upward in S1.
784     if (s1_size - s2_size != 0)
785     cy_limb = GMP_INTS_mpn_add_1 (res_ptr + s2_size,
786     s1_ptr + s2_size,
787     s1_size - s2_size,
788     cy_limb);
789     return cy_limb;
790     }
791    
792    
793     //07/15/01: Am willing to skip unit-testing on this.
794     //Understand the logic (i.e. passes visual inspection),
795     //and comes from GNU-MP.
796     GMP_INTS_limb_t GMP_INTS_mpn_sub_n(GMP_INTS_limb_ptr res_ptr,
797     GMP_INTS_limb_srcptr s1_ptr,
798     GMP_INTS_limb_srcptr s2_ptr,
799     GMP_INTS_size_t size)
800     {
801     GMP_INTS_limb_t x, y, cy;
802     GMP_INTS_size_t j;
803    
804     assert(res_ptr != NULL);
805     assert(s1_ptr != NULL);
806     assert(s2_ptr != NULL);
807    
808     /* The loop counter and index J goes from -SIZE to -1. This way
809     the loop becomes faster. */
810     j = -size;
811    
812     /* Offset the base pointers to compensate for the negative indices. */
813     s1_ptr -= j;
814     s2_ptr -= j;
815     res_ptr -= j;
816    
817     cy = 0;
818     do
819     {
820     y = s2_ptr[j];
821     x = s1_ptr[j];
822     y += cy; /* add previous carry to subtrahend */
823     cy = (y < cy); /* get out carry from that addition */
824     y = x - y; /* main subtract */
825     cy = (y > x) + cy; /* get out carry from the subtract, combine */
826     res_ptr[j] = y;
827     }
828     while (++j != 0);
829    
830     return cy;
831     }
832    
833    
834     //07/17/01: Am willing to skip unit-testing on this.
835     //Understand the logic (i.e. passes visual inspection),
836     //and comes from GNU-MP.
837     GMP_INTS_limb_t GMP_INTS_mpn_sub (GMP_INTS_limb_ptr res_ptr,
838     GMP_INTS_limb_srcptr s1_ptr,
839     GMP_INTS_size_t s1_size,
840     GMP_INTS_limb_srcptr s2_ptr,
841     GMP_INTS_size_t s2_size)
842     {
843     GMP_INTS_limb_t cy_limb = 0;
844    
845     assert(res_ptr != NULL);
846     assert(s1_ptr != NULL);
847     assert(s2_ptr != NULL);
848    
849     if (s2_size != 0)
850     cy_limb = GMP_INTS_mpn_sub_n(res_ptr,
851     s1_ptr,
852     s2_ptr,
853     s2_size);
854    
855     if (s1_size - s2_size != 0)
856     cy_limb = GMP_INTS_mpn_sub_1(res_ptr + s2_size,
857     s1_ptr + s2_size,
858     s1_size - s2_size,
859     cy_limb);
860     return cy_limb;
861     }
862    
863    
864     //07/17/01: Am willing to skip unit-testing on this.
865     //Understand the logic (i.e. passes visual inspection),
866     //and comes from GNU-MP.
867     GMP_INTS_limb_t GMP_INTS_mpn_lshift(GMP_INTS_limb_ptr wp,
868     GMP_INTS_limb_srcptr up,
869     GMP_INTS_size_t usize,
870     unsigned int cnt)
871     {
872     GMP_INTS_limb_t high_limb, low_limb;
873     unsigned sh_1, sh_2;
874     GMP_INTS_size_t i;
875     GMP_INTS_limb_t retval;
876    
877     assert(wp != NULL);
878     assert(up != NULL);
879     assert(usize > 0);
880     assert(cnt > 0);
881    
882     sh_1 = cnt;
883    
884     wp += 1;
885     sh_2 = GMP_INTS_BITS_PER_LIMB - sh_1;
886     //This automatically implies that can't call this function to shift more
887     //than 31 places.
888     i = usize - 1;
889     low_limb = up[i];
890     retval = low_limb >> sh_2; //Return value is the amount shifted
891     //off the top.
892     high_limb = low_limb;
893     while (--i >= 0)
894     {
895     low_limb = up[i];
896     wp[i] = (high_limb << sh_1) | (low_limb >> sh_2);
897     high_limb = low_limb;
898     }
899     wp[i] = high_limb << sh_1;
900    
901     return retval;
902     }
903    
904    
905     //07/17/01: Am willing to skip unit-testing on this.
906     //Understand the logic more or less (i.e. passes visual inspection),
907     //and comes from GNU-MP.
908     /* Shift U (pointed to by UP and USIZE limbs long) CNT bits to the right
909     and store the USIZE least significant limbs of the result at WP.
910     The bits shifted out to the right are returned.
911    
912     Argument constraints:
913     1. 0 < CNT < BITS_PER_MP_LIMB
914     2. If the result is to be written over the input, WP must be <= UP.
915     */
916     GMP_INTS_limb_t GMP_INTS_mpn_rshift (GMP_INTS_limb_ptr wp,
917     GMP_INTS_limb_srcptr up,
918     GMP_INTS_size_t usize,
919     unsigned int cnt)
920     {
921     GMP_INTS_limb_t high_limb, low_limb;
922     unsigned sh_1, sh_2;
923     GMP_INTS_size_t i;
924     GMP_INTS_limb_t retval;
925    
926     assert(wp != NULL);
927     assert(up != NULL);
928     assert(usize > 0);
929     assert(cnt > 0);
930    
931     sh_1 = cnt;
932    
933     wp -= 1;
934     sh_2 = GMP_INTS_BITS_PER_LIMB - sh_1;
935     high_limb = up[0];
936     retval = high_limb << sh_2;
937     low_limb = high_limb;
938    
939     for (i = 1; i < usize; i++)
940     {
941     high_limb = up[i];
942     wp[i] = (low_limb >> sh_1) | (high_limb << sh_2);
943     low_limb = high_limb;
944     }
945     wp[i] = low_limb >> sh_1;
946    
947     return retval;
948     }
949    
950    
951     //07/17/01: Am willing to skip unit-testing on this.
952     //Understand the logic (i.e. passes visual inspection),
953     //and comes from GNU-MP.
954     int GMP_INTS_mpn_cmp (GMP_INTS_limb_srcptr op1_ptr,
955     GMP_INTS_limb_srcptr op2_ptr,
956     GMP_INTS_size_t size)
957     {
958     GMP_INTS_size_t i;
959     GMP_INTS_limb_t op1_word, op2_word;
960    
961     assert(op1_ptr != NULL);
962     assert(op2_ptr != NULL);
963    
964     for (i = size - 1; i >= 0; i--)
965     {
966     op1_word = op1_ptr[i];
967     op2_word = op2_ptr[i];
968     if (op1_word != op2_word)
969     goto diff;
970     }
971     return 0;
972    
973     diff:
974     //This can *not* be simplified to
975     // op2_word - op2_word
976     //since that expression might give signed overflow.
977     return (op1_word > op2_word) ? 1 : -1;
978     }
979    
980    
981     //07/15/01: Am willing to skip unit-testing on this.
982     //Understand the logic (i.e. passes visual inspection),
983     //and comes from GNU-MP. Hope any defects here will be
984     //caught in testing of GMP_INTS_mpz_mul() and other
985     //higher-level functions.
986     void GMP_INTS_mpn_mul_basecase (GMP_INTS_limb_ptr prodp,
987     GMP_INTS_limb_srcptr up,
988     GMP_INTS_size_t usize,
989     GMP_INTS_limb_srcptr vp,
990     GMP_INTS_size_t vsize)
991     {
992     assert(prodp != NULL);
993     assert(up != NULL);
994     assert(usize > 0);
995     assert(vp != NULL);
996     assert(vsize > 0);
997    
998     /* We first multiply by the low order one or two limbs, as the result can
999     be stored, not added, to PROD. We also avoid a loop for zeroing this
1000     way. */
1001     prodp[usize] = GMP_INTS_mpn_mul_1 (prodp, up, usize, vp[0]);
1002     prodp++;
1003     vp++;
1004     vsize--;
1005    
1006     /* For each iteration in the loop, multiply U with one limb from V, and
1007     add the result to PROD. */
1008     while (vsize != 0)
1009     {
1010     prodp[usize] = GMP_INTS_mpn_addmul_1 (prodp, up, usize, vp[0]);
1011     prodp++,
1012     vp++,
1013     vsize--;
1014     }
1015     }
1016    
1017    
1018     //07/15/01: No unit testing possible--this is a passthrough.
1019     //In the original GNU MP code, there were several multiplication
1020     //algorithms, and this function would select one based on the
1021     //size of the operands and other considerations. The code has been
1022     //pared so that only simple multiplication is used, which is why
1023     //this function contains only a single pass-thru function call.
1024     void GMP_INTS_mpn_mul_n (GMP_INTS_limb_ptr p,
1025     GMP_INTS_limb_srcptr a,
1026     GMP_INTS_limb_srcptr b,
1027     GMP_INTS_size_t n)
1028     {
1029     GMP_INTS_mpn_mul_basecase (p, a, n, b, n);
1030     }
1031    
1032    
1033     //07/16/01: Visual inspection OK. Will not perform unit testing.
1034     GMP_INTS_limb_t GMP_INTS_mpn_mul(GMP_INTS_limb_ptr prodp,
1035     GMP_INTS_limb_srcptr up,
1036     GMP_INTS_size_t un,
1037     GMP_INTS_limb_srcptr vp,
1038     GMP_INTS_size_t vn)
1039     {
1040     //This is a gutted version of the GNU MP function. The GNU
1041     //MP function considered the case of a square, and also
1042     //better algorithms that pay off with large operands.
1043     //This gutted version uses only basic multiplication
1044     //(O(N**2)).
1045    
1046     //Eyeball the input parameters.
1047     assert(prodp != NULL);
1048     assert(up != NULL);
1049     assert(un >= 0);
1050     assert(vp != NULL);
1051     assert(vn >= 0);
1052    
1053     /* Basic long multiplication. */
1054     GMP_INTS_mpn_mul_basecase (prodp, up, un, vp, vn);
1055    
1056     //Return the most significant limb (which might be zero).
1057     //This is different than
1058     //most other functions, which return the spillover.
1059     return prodp[un + vn - 1];
1060     }
1061    
1062    
1063     /******************************************************************/
1064     /*** LIMB SPACE REALLOCATION FUNCTIONS *************************/
1065     /******************************************************************/
1066    
1067     void *GMP_INTS_mpz_realloc (GMP_INTS_mpz_struct *m,
1068     GMP_INTS_size_t new_size)
1069     {
1070     /* Never allocate zero space. */
1071     if (new_size <= 0)
1072     new_size = 1;
1073    
1074     m->limbs = (GMP_INTS_limb_ptr)
1075     GMP_INTS_realloc_w_size (m->limbs,
1076     m->n_allocd * sizeof(GMP_INTS_limb_t),
1077     new_size * sizeof(GMP_INTS_limb_t));
1078     m->n_allocd = new_size;
1079    
1080     return (void *) m->limbs;
1081     }
1082    
1083    
1084     /******************************************************************/
1085     /*** PUBLIC INITIALIZATION AND MEMORY MANAGEMENT FUNCTIONS *****/
1086     /******************************************************************/
1087    
1088     void GMP_INTS_mpz_init (GMP_INTS_mpz_struct *x)
1089     {
1090     assert(x != NULL);
1091    
1092     //The structure (the header block) exists in the
1093     //caller's area. Most likely it is a local variable.
1094     //This is OK, because it doesn't take up much space.
1095    
1096     //Start off with no errors.
1097     x->flags = 0;
1098    
1099     //Allocate space for one limb, which is the most
1100     //basic amount. This will grow, almost certainly.
1101     x->limbs = GMP_INTS_malloc(sizeof(GMP_INTS_limb_t));
1102    
1103     //Indicate that one limb was allocated.
1104     x->n_allocd = 1;
1105    
1106     //Set the size to 0. This signals a value of zero.
1107     x->size = 0;
1108     }
1109    
1110    
1111     void GMP_INTS_mpz_clear (GMP_INTS_mpz_struct *x)
1112     {
1113     //Be sure the passed pointer is not NULL.
1114     assert(x != NULL);
1115    
1116     //Be sure that the amount allocated is also above zero.
1117     //Anything else represents a logical error.
1118     assert(x->n_allocd > 0);
1119    
1120     //Be sure that the pointer to the allocated limbs
1121     //is not NULL. Anything else would be a logical
1122     //error.
1123     assert(x->limbs != NULL);
1124    
1125     //Deallocate the space for the limbs. The pointer is
1126     //set NULL and the allocated amount set to zero
1127     // so in case clear is called again it will be
1128     //a detectable error.
1129     GMP_INTS_free_w_size(x->limbs,
1130     x->n_allocd * sizeof(GMP_INTS_limb_t));
1131     x->limbs = NULL;
1132     x->n_allocd = 0;
1133     }
1134    
1135    
1136     /******************************************************************/
1137     /*** PUBLIC ASSIGNMENT FUNCTIONS *******************************/
1138     /******************************************************************/
1139    
1140     void GMP_INTS_mpz_copy( GMP_INTS_mpz_struct *dst,
1141     const GMP_INTS_mpz_struct *src)
1142     {
1143     GMP_INTS_size_t i, n;
1144    
1145     //Eyeball the input parameters.
1146     assert(dst != NULL);
1147     assert(dst->n_allocd > 0);
1148     assert(dst->limbs != NULL);
1149     assert(src != NULL);
1150     assert(src->n_allocd > 0);
1151     assert(src->limbs != NULL);
1152    
1153     //Source and destination may not be the same.
1154     assert(src != dst);
1155    
1156     //Figure out the real size of the source. We need to take the absolute
1157     //value.
1158     n = GMP_INTS_abs_of_size_t(src->size);
1159    
1160     //Reallocate the destination to be bigger if necessary.
1161     if (dst->n_allocd < n)
1162     {
1163     GMP_INTS_mpz_realloc (dst, n);
1164     }
1165    
1166     //Copy the non-dynamic fields in the header.
1167     dst->flags = src->flags;
1168     dst->size = src->size;
1169    
1170     //Copy the limbs.
1171     for (i=0; i<n; i++)
1172     dst->limbs[i] = src->limbs[i];
1173     }
1174    
1175    
1176     void GMP_INTS_mpz_set_ui (GMP_INTS_mpz_struct *dest,
1177     unsigned long int val)
1178     {
1179     assert(dest != NULL);
1180    
1181     /* We don't check if the allocation is enough, since the rest of the
1182     package ensures it's at least 1, which is what we need here. */
1183    
1184     dest->flags = 0;
1185     //A set operation resets any errors.
1186    
1187     if (val > 0)
1188     {
1189     dest->limbs[0] = val;
1190     dest->size = 1;
1191     }
1192     else
1193     {
1194     dest->size = 0;
1195     }
1196     }
1197    
1198    
1199     void GMP_INTS_mpz_set_si (GMP_INTS_mpz_struct *dest,
1200     signed long int val)
1201     {
1202     assert(dest != NULL);
1203    
1204     /* We don't check if the allocation is enough, since the rest of the
1205     package ensures it's at least 1, which is what we need here. */
1206    
1207     dest->flags = 0;
1208     //A set operation resets any errors.
1209    
1210     if (val > 0)
1211     {
1212     dest->limbs[0] = val;
1213     dest->size = 1;
1214     }
1215     else if (val < 0)
1216     {
1217     dest->limbs[0] = (unsigned long) -val;
1218     dest->size = -1;
1219     }
1220     else
1221     {
1222     dest->size = 0;
1223     }
1224     }
1225    
1226    
1227     void GMP_INTS_mpz_set_simple_char_str(GMP_INTS_mpz_struct *z,
1228     const char *s)
1229     {
1230     int sign=1;
1231     int digval;
1232     GMP_INTS_mpz_struct digvalz, k10;
1233    
1234     //Eyeball the arguments.
1235     assert(z != NULL);
1236     assert(z->n_allocd > 0);
1237     assert(z->limbs != NULL);
1238     assert(s != NULL);
1239    
1240     //Set the arbitrary integer to zero. This will also kill
1241     //any error flags.
1242     GMP_INTS_mpz_set_ui(z, 0);
1243    
1244     //Allocate an integer for our private use to hold each digit
1245     //value.
1246     GMP_INTS_mpz_init(&digvalz);
1247    
1248     //Allocate the constant 10, which we will use often.
1249     GMP_INTS_mpz_init(&k10);
1250     GMP_INTS_mpz_set_ui(&k10, 10);
1251    
1252     //As long as there are are digits and no flags set, keep
1253     //multiplying and adding the value of the digit. Non-
1254     //digits are simply ignored.
1255     while (!(z->flags) && (*s))
1256     {
1257     if (*s == '-')
1258     {
1259     sign = -sign;
1260     }
1261     else
1262     {
1263     digval = CHARFUNC_digit_to_val(*s);
1264     if (digval >= 0)
1265     {
1266     GMP_INTS_mpz_set_ui(&digvalz, digval);
1267     GMP_INTS_mpz_mul(z, z, &k10);
1268     GMP_INTS_mpz_add(z, z, &digvalz);
1269     }
1270     }
1271     s++;
1272     }
1273    
1274     //Adjust the final sign of the result.
1275     if (sign < 0)
1276     z->size = -(z->size);
1277    
1278     //Deallocate our temporary integers.
1279     GMP_INTS_mpz_clear(&digvalz);
1280     GMP_INTS_mpz_clear(&k10);
1281     }
1282    
1283    
1284     void GMP_INTS_mpz_set_sci_not_num(GMP_INTS_mpz_struct *z,
1285     int *failure,
1286     const char *s)
1287     {
1288     int parse_failure;
1289     //Return code from the floating point parsing
1290     //function.
1291     char mant_sign;
1292     //Sign character, if any, from the mantissa,
1293     //or N otherwise.
1294     size_t mant_bdp;
1295     //The index to the start of the mantissa before
1296     //the decimal point.
1297     size_t mant_bdp_len;
1298     //The length of the mantissa before the decimal
1299     //point. Zero means not defined, i.e. that
1300     //no characters were parsed and interpreted as
1301     //that part of a floating point number.
1302     size_t mant_adp;
1303     size_t mant_adp_len;
1304     //Similar fields for after the decimal point.
1305     char exp_sign;
1306     //Sign of the exponent, if any, or N otherwise.
1307     size_t exp;
1308     size_t exp_len;
1309     //Similar fields as to the mantissa, but for the
1310     //exponent.
1311     size_t si;
1312     //Iteration variable.
1313     int exponent_val;
1314     //The value of the exponent. We can't accept
1315     //an exponent outside the range of a 24-bit
1316     //signed integer. The 24-bit limit is arbitrary.
1317     //For one thing, it gives room to detect overflow
1318     //as are adding and multiplying by 10.
1319    
1320     //Eyeball the input parameters.
1321     assert(z != NULL);
1322     assert(z->n_allocd > 0);
1323     assert(z->limbs != NULL);
1324     assert(failure != NULL);
1325     assert(s != NULL);
1326    
1327     //Start off believing no failure.
1328     *failure = 0;
1329    
1330     //Set the output to zero. This is the default case for some
1331     //steps below.
1332     GMP_INTS_mpz_set_ui(z, 0);
1333    
1334     //Attempt to parse the number as a general number
1335     //in scientific notation.
1336     BSTRFUNC_parse_gen_sci_not_num(s,
1337     &parse_failure,
1338     &mant_sign,
1339     &mant_bdp,
1340     &mant_bdp_len,
1341     &mant_adp,
1342     &mant_adp_len,
1343     &exp_sign,
1344     &exp,
1345     &exp_len);
1346    
1347     //If it wouldn't parse as a general number, can't go further.
1348     if (parse_failure)
1349     {
1350     *failure = 1;
1351     return;
1352     }
1353     else if (!exp_len && !mant_adp_len)
1354     {
1355     //There was no exponent, and no portion after
1356     //the decimal point. Can just parse as an integer.
1357     char *temp_buf;
1358    
1359     //Allocate the temporary buffer to be one character longer
1360     //than the length specified for the parsed mantissa.
1361     temp_buf = GMP_INTS_malloc(sizeof(char) * (mant_bdp_len + 1));
1362    
1363     //Copy from the parsed area into the temporary buffer.
1364     for (si=mant_bdp; si<(mant_bdp+mant_bdp_len); si++)
1365     temp_buf[si-mant_bdp] = s[si];
1366     temp_buf[mant_bdp_len] = 0;
1367    
1368     //Set the arbitrary integer to the value of the character
1369     //string.
1370     GMP_INTS_mpz_set_simple_char_str(z, temp_buf);
1371    
1372     //If the number parsed as negative, invert.
1373     if (mant_sign == '-')
1374     z->size = -z->size;
1375    
1376     //Deallocate the temporary buffer.
1377     GMP_INTS_free(temp_buf);
1378     }
1379     else if (!exp_len && mant_adp_len)
1380     {
1381     char *temp_buf;
1382    
1383     //In this case, there are digits after the decimal point,
1384     //but no exponent specified. The only way this makes
1385     //sense is if all of the digits are zero--otherwise it
1386     //cannot be an integer.
1387     for (si=mant_adp; si<(mant_adp+mant_adp_len); si++)
1388     {
1389     if (s[si] != '0')
1390     {
1391     *failure = 1;
1392     return;
1393     }
1394     }
1395    
1396     //We're clean. They are only zeros. Execute as per
1397     //integer code.
1398    
1399     //Allocate the temporary buffer to be one character longer
1400     //than the length specified for the parsed mantissa.
1401     temp_buf = GMP_INTS_malloc(sizeof(char) * (mant_bdp_len + 1));
1402    
1403     //Copy from the parsed area into the temporary buffer.
1404     for (si=mant_bdp; si<(mant_bdp+mant_bdp_len); si++)
1405     temp_buf[si-mant_bdp] = s[si];
1406     temp_buf[mant_bdp_len] = 0;
1407    
1408     //Set the arbitrary integer to the value of the character
1409     //string.
1410     GMP_INTS_mpz_set_simple_char_str(z, temp_buf);
1411    
1412     //If the number parsed as negative, invert.
1413     if (mant_sign == '-')
1414     z->size = -z->size;
1415    
1416     //Deallocate the temporary buffer.
1417     GMP_INTS_free(temp_buf);
1418     }
1419     else if (exp_len)
1420     {
1421     //This is the most difficult case, where an exponent
1422     //is specified. There are several complex subcases,
1423     //such as:
1424     // a)If the exponent is too positive or too negative,
1425     // we can't use it. In general, we won't tackle
1426     // an exponent that won't fit in a signed 24-bit
1427     // integer. This provides a range of from
1428     // -8,388,608 to +8,388,607. This dwarfs the
1429     // 100,000 or so digit preprocessor limit,
1430     // and should be adequate for any practical
1431     // application.
1432     // b)If the exponent is zero, we ignore it.
1433     // c)If the exponent is positive, it has to
1434     // be large enough to overcome any
1435     // digits past the decimal point, otherwise
1436     // we don't end up with an integer.
1437     // d)If the exponent is negative, there have to
1438     // be enough digits so that an integer remains
1439     // after the exponent is applied. This
1440     // generally requires trailing zeros on the
1441     // part before the decimal point.
1442    
1443     //First, tackle the exponent. Process the
1444     //exponent into a signed integer. We have to
1445     //balk at anything outside of 24 bits. The
1446     //procedure used automatically handles
1447     //leading zeros correctly.
1448     exponent_val = 0;
1449     for (si=exp; si<(exp+exp_len); si++)
1450     {
1451     int val;
1452    
1453     val = CHARFUNC_digit_to_val(s[si]);
1454    
1455     assert(val >= 0 && val <= 9);
1456    
1457     exponent_val *= 10;
1458     exponent_val += val;
1459    
1460     if (((exp_sign=='-') && (exponent_val>8388608))
1461     ||
1462     ((exp_sign != '-') && (exponent_val>8388607)))
1463     {
1464     *failure = 1;
1465     return;
1466     }
1467     }
1468    
1469     //If we're here, the exponent has been computed and
1470     //is within 24 bits. However, we need to adjust for
1471     //the sign.
1472     if (exp_sign == '-')
1473     exponent_val = -exponent_val;
1474    
1475     //We need to make accurate assertions about the
1476     //portion of the number, if any, after the decimal point.
1477     //This means that we need to effectively discard
1478     //trailing zeros. To do this, we do not need to
1479     //relocate the string, we can just back off the index
1480     //to bypass any trailing zeros.
1481     while ((mant_adp_len > 0) && (s[mant_adp + mant_adp_len - 1]=='0'))
1482     mant_adp_len--;
1483    
1484     //We also need to make accurate assertions about the
1485     //portion of the number, if any, before the decimal
1486     //point. It is known that the parsing function
1487     //isn't tolerant of multiple zeros, but zero is a
1488     //special case. Let's advance the pointer to the
1489     //part before the decimal point so that zero will
1490     //have zero length.
1491     while ((mant_bdp_len > 0) && (s[mant_bdp]=='0'))
1492     {
1493     mant_bdp++;
1494     mant_bdp_len--;
1495     }
1496    
1497     //If we are dealing with zero, who cares about the
1498     //exponent? Just return the value of zero.
1499     if (!mant_bdp_len && !mant_adp_len)
1500     {
1501     *failure = 0;
1502     GMP_INTS_mpz_set_ui(z, 0);
1503     return;
1504     }
1505    
1506     //Beyond this point, we have something non-zero.
1507     //If the exponent is positive, it must be at least
1508     //as large as the number of digits beyond the
1509     //decimal point in order to form an integer. If the
1510     //exponent is zero, there must be no digits after the
1511     //decimal point. If the exponent is negative, there
1512     //must be no digits after the decimal point, and the
1513     //trailing zeros on the part before the decimal point
1514     //must be adequate to handle the right decimal shift.
1515     if (exponent_val == 0)
1516     {
1517     if (mant_adp_len)
1518     {
1519     *failure = 1;
1520     return;
1521     }
1522     }
1523     else if (exponent_val > 0)
1524     {
1525     if ((int)mant_adp_len > exponent_val)
1526     {
1527     *failure = 1;
1528     return;
1529     }
1530     }
1531     else //exponent_val < 0
1532     {
1533     if (mant_adp_len)
1534     {
1535     *failure = 1;
1536     return;
1537     }
1538     else
1539     {
1540     //Count the number of trailing zeros on the part
1541     //before the decimal point.
1542     size_t trailing_zero_count;
1543     int idx;
1544    
1545     trailing_zero_count = 0;
1546    
1547     for(idx = mant_bdp + mant_bdp_len - 1;
1548     (mant_bdp_len != 0) && (idx >= (int)mant_bdp);
1549     idx--)
1550     {
1551     if (s[idx] == '0')
1552     trailing_zero_count++;
1553     else
1554     break;
1555     }
1556    
1557     //Check on the assertion about trailing zeros.
1558     if ((int)trailing_zero_count < -exponent_val)
1559     {
1560     *failure = 1;
1561     return;
1562     }
1563     }
1564     }
1565    
1566     {
1567     //Create a string long enough to hold the digits
1568     //before the decimal point plus the ones after and
1569     //convert that to an arbitrary integer.
1570     //Form a power of 10 which is 10 exponentiated to
1571     //the absolute value of the exponent. If the
1572     //exponent was positive, multiply by it. If the
1573     //exponent was negative, divide by it.
1574     char *conv_str;
1575     size_t sidx;
1576     GMP_INTS_mpz_struct power_of_ten, k10, trash;
1577    
1578     GMP_INTS_mpz_init(&power_of_ten);
1579     GMP_INTS_mpz_init(&k10);
1580     GMP_INTS_mpz_init(&trash);
1581    
1582     conv_str = GMP_INTS_malloc(sizeof(char) * (mant_bdp_len + mant_adp_len + 1));
1583    
1584     sidx=0;
1585    
1586     for (si=mant_bdp; si<(mant_bdp+mant_bdp_len); si++)
1587     {
1588     conv_str[sidx] = s[si];
1589     sidx++;
1590     }
1591     for (si=mant_adp; si<(mant_adp+mant_adp_len); si++)
1592     {
1593     conv_str[sidx] = s[si];
1594     sidx++;
1595     }
1596     conv_str[sidx] = 0;
1597    
1598     assert(sidx == (mant_bdp_len + mant_adp_len));
1599    
1600     GMP_INTS_mpz_set_simple_char_str(z, conv_str);
1601    
1602     GMP_INTS_mpz_set_ui(&k10, 10);
1603    
1604     if (exponent_val > 0)
1605     GMP_INTS_mpz_pow_ui(&power_of_ten, &k10, exponent_val-mant_adp_len);
1606     else
1607     GMP_INTS_mpz_pow_ui(&power_of_ten, &k10, -exponent_val);
1608    
1609     if (exponent_val >= 0)
1610     {
1611     GMP_INTS_mpz_mul(z, z, &power_of_ten);
1612     }
1613     else
1614     {
1615     GMP_INTS_mpz_tdiv_qr (&k10,
1616     &trash,
1617     z,
1618     &power_of_ten);
1619     GMP_INTS_mpz_copy(z, &k10);
1620     }
1621    
1622     //If the argument had a minus sign, invert.
1623     if (mant_sign == '-')
1624     z->size = -z->size;
1625    
1626     GMP_INTS_free(conv_str);
1627    
1628     GMP_INTS_mpz_clear(&trash);
1629     GMP_INTS_mpz_clear(&k10);
1630     GMP_INTS_mpz_clear(&power_of_ten);
1631    
1632     //Finally, if the arbitrary integer has overflowed, this is
1633     //a parse failure. Must declare as such.
1634     if (z->flags)
1635     *failure = 1;
1636     }
1637     }
1638     else
1639     {
1640     *failure = 1;
1641     return;
1642     }
1643     }
1644    
1645    
1646     void GMP_INTS_mpz_set_general_int(GMP_INTS_mpz_struct *z,
1647     int *failure,
1648     const char *s)
1649     {
1650     //Eyeball the input parameters.
1651     assert(z != NULL);
1652     assert(z->n_allocd > 0);
1653     assert(z->limbs != NULL);
1654     assert(failure != NULL);
1655     assert(s != NULL);
1656    
1657     //Try to parse it as a simple integer.
1658     if (BSTRFUNC_is_sint_wo_commas(s))
1659     {
1660     GMP_INTS_mpz_set_simple_char_str(z, s);
1661     *failure = 0;
1662     return;
1663     }
1664     //If that didn't work, try to parse it as a simple
1665     //integer with commas.
1666     else if (BSTRFUNC_is_sint_w_commas(s))
1667     {
1668     GMP_INTS_mpz_set_simple_char_str(z, s);
1669     *failure = 0;
1670     return;
1671     }
1672     //If neither of those worked, try to parse it as
1673     //something containing scientific notation.
1674     else
1675     {
1676     GMP_INTS_mpz_set_sci_not_num(z, failure, s);
1677    
1678     if (!*failure)
1679     {
1680     //We were able to parse it that way.
1681     //Everything is set up, just return.
1682     return;
1683     }
1684     else
1685     {
1686     //We're out of options. All parsing failed.
1687     GMP_INTS_mpz_set_ui(z, 0);
1688     *failure = 1;
1689     return;
1690     }
1691     }
1692     }
1693    
1694    
1695     void GMP_INTS_mpz_parse_into_uint32(unsigned *result,
1696     int *failure,
1697     char *s)
1698     {
1699     GMP_INTS_mpz_struct arb_int;
1700    
1701     //Eyeball the input parameters.
1702     assert(result != NULL);
1703     assert(failure != NULL);
1704     assert(s != NULL);
1705    
1706     //Allocate space for the one arbitrary integer we need.
1707     GMP_INTS_mpz_init(&arb_int);
1708    
1709     //Try to parse the string into an arbitrary length integer
1710     //using all methods known to man.
1711     GMP_INTS_mpz_set_general_int(&arb_int, failure, s);
1712    
1713     //If the parse failed, we must declare failure and return
1714     //0.
1715     if (*failure)
1716     {
1717     *result = 0;
1718     *failure = 1;
1719     }
1720     else
1721     {
1722     //We might have success, but it might be negative or
1723     //too big.
1724     if (arb_int.size == 1)
1725     {
1726     *result = arb_int.limbs[0];
1727     *failure = 0;
1728     }
1729     else if (arb_int.size == 0)
1730     {
1731     *result = 0;
1732     *failure = 0;
1733     }
1734     else
1735     {
1736     *result = 0;
1737     *failure = 1;
1738     }
1739     }
1740    
1741     //Deallocate the arbitrary integer.
1742     GMP_INTS_mpz_clear(&arb_int);
1743     }
1744    
1745    
1746     void GMP_INTS_mpz_swap(GMP_INTS_mpz_struct *a,
1747     GMP_INTS_mpz_struct *b)
1748     {
1749     GMP_INTS_mpz_struct temp;
1750    
1751     //Eyeball the input parameters.
1752     assert(a != NULL);
1753     assert(a->n_allocd > 0);
1754     assert(a->limbs != NULL);
1755     assert(b != NULL);
1756     assert(b->n_allocd > 0);
1757     assert(b->limbs != NULL);
1758    
1759     //Make the swap via memory copy.
1760     memcpy(&temp, a, sizeof(GMP_INTS_mpz_struct));
1761     memcpy(a, b, sizeof(GMP_INTS_mpz_struct));
1762     memcpy(b, &temp, sizeof(GMP_INTS_mpz_struct));
1763     }
1764    
1765    
1766     /******************************************************************/
1767     /*** PUBLIC ARITHMETIC FUNCTIONS *******************************/
1768     /******************************************************************/
1769    
1770     //07/15/01: Unit test and visual inspection passed.
1771     void GMP_INTS_mpz_add ( GMP_INTS_mpz_struct *w,
1772     const GMP_INTS_mpz_struct *u,
1773     const GMP_INTS_mpz_struct *v)
1774     {
1775     GMP_INTS_limb_srcptr up, vp;
1776     GMP_INTS_limb_ptr wp;
1777     GMP_INTS_size_t usize, vsize, wsize;
1778     GMP_INTS_size_t abs_usize;
1779     GMP_INTS_size_t abs_vsize;
1780    
1781     //Look at the input parameters carefully.
1782     assert(w != NULL);
1783     assert(u != NULL);
1784     assert(v != NULL);
1785     assert(w->n_allocd > 0);
1786     assert(u->n_allocd > 0);
1787     assert(v->n_allocd > 0);
1788     assert(w->limbs != NULL);
1789     assert(u->limbs != NULL);
1790     assert(v->limbs != NULL);
1791    
1792     //Handle the case of a tainted result. If either of the
1793     //two inputs are either direct overflows or tainted by
1794     //an overflow, mark the result tainted and do not perform
1795     //any arithmetic operation.
1796     {
1797     int taint;
1798    
1799     taint = GMP_INTS_two_op_flags_map(u->flags, v->flags);
1800    
1801     w->flags = 0;
1802     //"w" starts off with a clean slate. Must do this
1803     //after taint calculation in case locations of u or v
1804     //are the same as w.
1805     if (taint)
1806     {
1807     w->flags = taint;
1808     return;
1809     }
1810     }
1811    
1812     usize = u->size;
1813     vsize = v->size;
1814     abs_usize = GMP_INTS_abs_of_size_t(usize);
1815     abs_vsize = GMP_INTS_abs_of_size_t(vsize);
1816    
1817     //Arrange things so that U has at least as many
1818     //limbs as V, i.e. limbs(U) >= limbs(V);
1819     if (abs_usize < abs_vsize)
1820     {
1821     const GMP_INTS_mpz_struct *tmp_ptr;
1822     GMP_INTS_size_t tmp_size;
1823    
1824     //Swap U and V. This does no harm, because we are
1825     //manipulating only local variables. This does not
1826     //affect the caller.
1827     tmp_ptr = u;
1828     u = v;
1829     v = tmp_ptr;
1830     tmp_size = usize;
1831     usize = vsize;
1832     vsize = tmp_size;
1833     tmp_size = abs_usize;
1834     abs_usize = abs_vsize;
1835     abs_vsize = tmp_size;
1836     }
1837    
1838     /* True: ABS_USIZE >= ABS_VSIZE. */
1839    
1840     /* If not space for w (and possible carry), increase space. */
1841     wsize = abs_usize + 1;
1842     if (w->n_allocd < wsize)
1843     GMP_INTS_mpz_realloc(w, wsize);
1844    
1845     //These pointers must be obtained after realloc. At this point,
1846     //u or v may be the same as w.
1847     up = u->limbs;
1848     vp = v->limbs;
1849     wp = w->limbs;
1850    
1851     if ((usize ^ vsize) < 0)
1852     {
1853     //U and V have different sign. Need to compare them to determine
1854     //which operand to subtract from which.
1855    
1856     //This test is right since ABS_USIZE >= ABS_VSIZE.
1857     //If the equality case is ruled out, then U has more limbs
1858     //than V, which means that it is bigger in magnitude.
1859     if (abs_usize != abs_vsize)
1860     {
1861     GMP_INTS_mpn_sub (wp, up, abs_usize, vp, abs_vsize);
1862     wsize = abs_usize;
1863    
1864     //Normalize the result. This was formerly a macro.
1865     //To normalize in this context means to trim the size
1866     //down to eliminate any leading zero limbs that came
1867     //about because the size of the result of an operation
1868     //was overestimated.
1869     GMP_INTS_mpn_normalize(wp, &wsize);
1870    
1871     if (usize < 0)
1872     wsize = -wsize;
1873     }
1874     else if (GMP_INTS_mpn_cmp (up, vp, abs_usize) < 0)
1875     {
1876     GMP_INTS_mpn_sub_n (wp, vp, up, abs_usize);
1877     wsize = abs_usize;
1878     GMP_INTS_mpn_normalize(wp, &wsize);
1879     if (usize >= 0)
1880     wsize = -wsize;
1881     }
1882     else
1883     {
1884     GMP_INTS_mpn_sub_n (wp, up, vp, abs_usize);
1885     wsize = abs_usize;
1886     GMP_INTS_mpn_normalize(wp, &wsize);
1887     if (usize < 0)
1888     wsize = -wsize;
1889     }
1890     }
1891     else
1892     {
1893     /* U and V have same sign. Add them. */
1894     GMP_INTS_limb_t cy_limb
1895     = GMP_INTS_mpn_add (wp, up, abs_usize, vp, abs_vsize);
1896     wp[abs_usize] = cy_limb;
1897     wsize = abs_usize + cy_limb;
1898     if (usize < 0)
1899     wsize = -wsize;
1900     }
1901    
1902     w->size = wsize;
1903    
1904     //Handle the case of an overflowed result. If the result
1905     //of the addition is too big or too small, mark it as
1906     //overflowed.
1907     if (w->size > GMP_INTS_MAXIMUM_LIMBS_PER_INT)
1908     {
1909     w->flags = GMP_INTS_EF_INTOVF_POS;
1910     }
1911     else if (w->size < -GMP_INTS_MAXIMUM_LIMBS_PER_INT)
1912     {
1913     w->flags = GMP_INTS_EF_INTOVF_NEG;
1914     }
1915     }
1916    
1917    
1918     //07/15/01: Unit testing skipped because of recursive
1919     //nature. Visual inspection OK.
1920     void GMP_INTS_mpz_add_ui ( GMP_INTS_mpz_struct *w,
1921     const GMP_INTS_mpz_struct *u,
1922     unsigned long int v)
1923     {
1924     //The GNU MP version of this is quite efficient, and this
1925     //makes sense since it is a common operation. However,
1926     //for simplicity just define this recursively in terms
1927     //of the ADD function. This can always be made quicker
1928     //later (by changing back to the GNU MP version).
1929     GMP_INTS_mpz_struct temp;
1930    
1931     //Eyeball the inputs carefully.
1932     assert(w != NULL);
1933     assert(w->n_allocd > 0);
1934     assert(w->limbs != NULL);
1935     assert(u != NULL);
1936     assert(u->n_allocd > 0);
1937     assert(u->limbs != NULL);
1938    
1939     //Create a temporary integer.
1940     GMP_INTS_mpz_init(&temp);
1941    
1942     //Set the temporary integer to the value of the input
1943     //argument.
1944     GMP_INTS_mpz_set_ui(&temp, v);
1945    
1946     //Do the actual addition. This recursive definition
1947     //is inherently wasteful, but I'm after clarity, not
1948     //warp speed.
1949     GMP_INTS_mpz_add(w, u, &temp);
1950    
1951     //Destroy the temporary integer (this will reclaim the
1952     //memory).
1953     GMP_INTS_mpz_clear(&temp);
1954     }
1955    
1956    
1957     //07/15/01: Visual inspection passed. Not unit tested
1958     //because of symmetry with GMP_INTS_mpz_add().
1959     void GMP_INTS_mpz_sub ( GMP_INTS_mpz_struct *w,
1960     const GMP_INTS_mpz_struct *u,
1961     const GMP_INTS_mpz_struct *v)
1962     {
1963     GMP_INTS_limb_srcptr up, vp;
1964     GMP_INTS_limb_ptr wp;
1965     GMP_INTS_size_t usize, vsize, wsize;
1966     GMP_INTS_size_t abs_usize;
1967     GMP_INTS_size_t abs_vsize;
1968    
1969     //Look at the input parameters carefully.
1970     assert(w != NULL);
1971     assert(u != NULL);
1972     assert(v != NULL);
1973     assert(w->n_allocd > 0);
1974     assert(u->n_allocd > 0);
1975     assert(v->n_allocd > 0);
1976     assert(w->limbs != NULL);
1977     assert(u->limbs != NULL);
1978     assert(v->limbs != NULL);
1979    
1980     //Handle the case of a tainted result. If either of the
1981     //two inputs are either direct overflows or tainted by
1982     //an overflow, mark the result tainted and do not perform
1983     //any arithmetic operation.
1984     {
1985     int taint;
1986    
1987     taint = GMP_INTS_two_op_flags_map(u->flags, v->flags);
1988    
1989     w->flags = 0;
1990     //"w" starts off with a clean slate. Must do this
1991     //after taint calculation in case locations of u or v
1992     //are the same as w.
1993     if (taint)
1994     {
1995     w->flags = taint;
1996     return;
1997     }
1998     }
1999    
2000     usize = u->size;
2001     vsize = -(v->size); /* The "-" makes the difference from mpz_add */
2002     abs_usize = GMP_INTS_abs_of_size_t(usize);
2003     abs_vsize = GMP_INTS_abs_of_size_t(vsize);
2004    
2005     if (abs_usize < abs_vsize)
2006     {
2007     const GMP_INTS_mpz_struct *tmp_ptr;
2008     GMP_INTS_size_t tmp_size;
2009    
2010     //Swap U and V. This does no harm, because we are
2011     //manipulating only local variables. This does not
2012     //affect the caller.
2013     tmp_ptr = u;
2014     u = v;
2015     v = tmp_ptr;
2016     tmp_size = usize;
2017     usize = vsize;
2018     vsize = tmp_size;
2019     tmp_size = abs_usize;
2020     abs_usize = abs_vsize;
2021     abs_vsize = tmp_size;
2022     }
2023    
2024     /* True: ABS_USIZE >= ABS_VSIZE. */
2025    
2026     /* If not space for w (and possible carry), increase space. */
2027     wsize = abs_usize + 1;
2028     if (w->n_allocd < wsize)
2029     GMP_INTS_mpz_realloc (w, wsize);
2030    
2031     /* These must be after realloc (u or v may be the same as w). */
2032     up = u->limbs;
2033     vp = v->limbs;
2034     wp = w->limbs;
2035    
2036     if ((usize ^ vsize) < 0)
2037     {
2038     //U and V have different sign. Need to compare them to determine
2039     //which operand to subtract from which.
2040    
2041     //This test is right since ABS_USIZE >= ABS_VSIZE.
2042     if (abs_usize != abs_vsize)
2043     {
2044     GMP_INTS_mpn_sub (wp, up, abs_usize, vp, abs_vsize);
2045     wsize = abs_usize;
2046     GMP_INTS_mpn_normalize(wp, &wsize);
2047     if (usize < 0)
2048     wsize = -wsize;
2049     }
2050     else if (GMP_INTS_mpn_cmp (up, vp, abs_usize) < 0)
2051     {
2052     GMP_INTS_mpn_sub_n (wp, vp, up, abs_usize);
2053     wsize = abs_usize;
2054     GMP_INTS_mpn_normalize(wp, &wsize);
2055     if (usize >= 0)
2056     wsize = -wsize;
2057     }
2058     else
2059     {
2060     GMP_INTS_mpn_sub_n (wp, up, vp, abs_usize);
2061     wsize = abs_usize;
2062     GMP_INTS_mpn_normalize (wp, &wsize);
2063     if (usize < 0)
2064     wsize = -wsize;
2065     }
2066     }
2067     else
2068     {
2069     /* U and V have same sign. Add them. */
2070     GMP_INTS_limb_t cy_limb
2071     = GMP_INTS_mpn_add (wp, up, abs_usize, vp, abs_vsize);
2072     wp[abs_usize] = cy_limb;
2073     wsize = abs_usize + cy_limb;
2074     if (usize < 0)
2075     wsize = -wsize;
2076     }
2077    
2078     w->size = wsize;
2079    
2080     //Handle the case of an overflowed result. If the result
2081     //of the addition is too big or too small, mark it as
2082     //overflowed.
2083     if (w->size > GMP_INTS_MAXIMUM_LIMBS_PER_INT)
2084     {
2085     w->flags = GMP_INTS_EF_INTOVF_POS;
2086     }
2087     else if (w->size < -GMP_INTS_MAXIMUM_LIMBS_PER_INT)
2088     {
2089     w->flags = GMP_INTS_EF_INTOVF_NEG;
2090     }
2091     }
2092    
2093    
2094     //07/15/01: Unit testing skipped because of recursive
2095     //nature. Visual inspection OK.
2096     void GMP_INTS_mpz_sub_ui ( GMP_INTS_mpz_struct *w,
2097     const GMP_INTS_mpz_struct *u,
2098     unsigned long int v)
2099     {
2100     //The GNU MP version of this is quite efficient, and this
2101     //makes sense since it is a common operation. However,
2102     //for simplicity just define this recursively in terms
2103     //of the SUB function. This can always be made quicker
2104     //later (by changing back to the GNU MP version).
2105     GMP_INTS_mpz_struct temp;
2106    
2107     //Eyeball the inputs carefully.
2108     assert(w != NULL);
2109     assert(w->n_allocd > 0);
2110     assert(w->limbs != NULL);
2111     assert(u != NULL);
2112     assert(u->n_allocd > 0);
2113     assert(u->limbs != NULL);
2114    
2115     //Create a temporary integer.
2116     GMP_INTS_mpz_init(&temp);
2117    
2118     //Set the temporary integer to the value of the input
2119     //argument.
2120     GMP_INTS_mpz_set_ui(&temp, v);
2121    
2122     //Do the actual subtraction. This recursive definition
2123     //is inherently wasteful, but I'm after clarity, not
2124     //warp speed.
2125     GMP_INTS_mpz_sub(w, u, &temp);
2126    
2127     //Destroy the temporary integer (this will reclaim the
2128     //memory).
2129     GMP_INTS_mpz_clear(&temp);
2130     }
2131    
2132    
2133     void GMP_INTS_mpz_mul ( GMP_INTS_mpz_struct *w,
2134     const GMP_INTS_mpz_struct *u,
2135     const GMP_INTS_mpz_struct *v)
2136     {
2137     GMP_INTS_size_t usize = u->size;
2138     GMP_INTS_size_t vsize = v->size;
2139     GMP_INTS_size_t wsize;
2140     GMP_INTS_size_t sign_product;
2141     GMP_INTS_limb_ptr up, vp;
2142     GMP_INTS_limb_ptr wp;
2143     GMP_INTS_limb_ptr free_me = NULL;
2144     GMP_INTS_size_t free_me_size;
2145     GMP_INTS_limb_t cy_limb;
2146    
2147     //Eyeball the inputs.
2148     assert(w != NULL);
2149     assert(w->n_allocd > 0);
2150     assert(w->limbs != NULL);
2151     assert(u != NULL);
2152     assert(u->n_allocd > 0);
2153     assert(u->limbs != NULL);
2154     assert(v != NULL);
2155     assert(v->n_allocd > 0);
2156     assert(v->limbs != NULL);
2157    
2158     //Handle the case of a tainted result. If either of the
2159     //two inputs are either direct overflows or tainted by
2160     //an overflow, mark the result tainted and do not perform
2161     //any arithmetic operation.
2162     {
2163     int taint;
2164    
2165     taint = GMP_INTS_two_op_flags_map(u->flags, v->flags);
2166    
2167     w->flags = 0;
2168     //"w" starts off with a clean slate. Must do this
2169     //after taint calculation in case locations of u or v
2170     //are the same as w.
2171     if (taint)
2172     {
2173     w->flags = taint;
2174     return;
2175     }
2176     }
2177    
2178     sign_product = usize ^ vsize;
2179     usize = GMP_INTS_abs_of_size_t(usize);
2180     vsize = GMP_INTS_abs_of_size_t(vsize);
2181    
2182     //Handle the case of a certain result overflow (why do the math when
2183     //the result is certain?). In general, when multiplying two inputs
2184     //whose sizes are M limbs and N limbs, the size of the result will be
2185     //either M+N or M+N-1 limbs. If M+N-1 > MAX_ALLOWED, then can declare
2186     //an early overflow.
2187     if ((usize + vsize - 1) > GMP_INTS_MAXIMUM_LIMBS_PER_INT)
2188     {
2189     if (sign_product < 0)
2190     w->flags = GMP_INTS_EF_INTOVF_NEG;
2191     else
2192     w->flags = GMP_INTS_EF_INTOVF_POS;
2193    
2194     return;
2195     }
2196    
2197    
2198     if (usize < vsize)
2199     {
2200     //Temporary variables just for the swap.
2201     const GMP_INTS_mpz_struct *tmp_ptr;
2202     GMP_INTS_size_t tmp_size;
2203    
2204     //Swap U and V.
2205     tmp_ptr = u;
2206     u = v;
2207     v = tmp_ptr;
2208     tmp_size = usize;
2209     usize = vsize;
2210     vsize = tmp_size;
2211     }
2212    
2213     //Grab pointers to the arrays of limbs.
2214     up = u->limbs;
2215     vp = v->limbs;
2216     wp = w->limbs;
2217    
2218     /* Ensure W has space enough to store the result. */
2219     wsize = usize + vsize;
2220     if (w->n_allocd < wsize)
2221     {
2222     if (wp == up || wp == vp)
2223     {
2224     free_me = wp;
2225     free_me_size = w->n_allocd;
2226     }
2227     else
2228     {
2229     GMP_INTS_free_w_size (wp, w->n_allocd * sizeof(GMP_INTS_limb_t));
2230     }
2231    
2232     w->n_allocd = wsize;
2233     wp = (GMP_INTS_limb_ptr)
2234     GMP_INTS_malloc (wsize * sizeof(GMP_INTS_limb_t));
2235     w->limbs = wp;
2236     }
2237     else
2238     {
2239     /* Make U and V not overlap with W. */
2240     if (wp == up)
2241     {
2242     /* W and U are identical. Allocate temporary space for U. */
2243     up = (GMP_INTS_limb_ptr)
2244     _alloca(usize * sizeof(GMP_INTS_limb_t));
2245     /* Is V identical too? Keep it identical with U. */
2246     if (wp == vp)
2247     vp = up;
2248     /* Copy to the temporary space. */
2249     GMP_INTS_mpn_copy_limbs(up, wp, usize);
2250     }
2251     else if (wp == vp)
2252     {
2253     /* W and V are identical. Allocate temporary space for V. */
2254     vp = (GMP_INTS_limb_ptr)
2255     _alloca(vsize * sizeof(GMP_INTS_limb_t));
2256     /* Copy to the temporary space. */
2257     GMP_INTS_mpn_copy_limbs(vp, wp, vsize);
2258     }
2259     }
2260    
2261     if (vsize == 0)
2262     {
2263     wsize = 0;
2264     }
2265     else
2266     {
2267     cy_limb = GMP_INTS_mpn_mul (wp, up, usize, vp, vsize);
2268     wsize = usize + vsize;
2269     wsize -= cy_limb == 0;
2270     }
2271    
2272     w->size = sign_product < 0 ? -wsize : wsize;
2273    
2274     if (free_me != NULL)
2275     GMP_INTS_free_w_size (free_me, free_me_size * sizeof(GMP_INTS_limb_t));
2276    
2277     //Final check for overflow.
2278     if (w->size > GMP_INTS_MAXIMUM_LIMBS_PER_INT)
2279     w->flags = GMP_INTS_EF_INTOVF_POS;
2280     else if (w->size < -GMP_INTS_MAXIMUM_LIMBS_PER_INT)
2281     w->flags = GMP_INTS_EF_INTOVF_NEG;
2282     }
2283    
2284    
2285     //07/15/01: Unit testing skipped because of recursive
2286     //nature. Visual inspection OK.
2287     void GMP_INTS_mpz_mul_si ( GMP_INTS_mpz_struct *w,
2288     const GMP_INTS_mpz_struct *u,
2289     long int v)
2290     {
2291     GMP_INTS_mpz_struct temp;
2292    
2293     //Eyeball the inputs carefully.
2294     assert(w != NULL);
2295     assert(w->n_allocd > 0);
2296     assert(w->limbs != NULL);
2297     assert(u != NULL);
2298     assert(u->n_allocd > 0);
2299     assert(u->limbs != NULL);
2300    
2301     //Create a temporary integer.
2302     GMP_INTS_mpz_init(&temp);
2303    
2304     //Set the temporary integer to the value of the input
2305     //argument.
2306     GMP_INTS_mpz_set_si(&temp, v);
2307    
2308     //Do the actual multiplication. This recursive definition
2309     //is inherently wasteful, but I'm after clarity, not
2310     //warp speed.
2311     GMP_INTS_mpz_mul(w, u, &temp);
2312    
2313     //Destroy the temporary integer (this will reclaim the
2314     //memory).
2315     GMP_INTS_mpz_clear(&temp);
2316     }
2317    
2318    
2319     //07/15/01: Unit testing skipped because of recursive
2320     //nature. Visual inspection OK.
2321     void GMP_INTS_mpz_mul_ui ( GMP_INTS_mpz_struct *w,
2322     const GMP_INTS_mpz_struct *u,
2323     unsigned long int v)
2324     {
2325     GMP_INTS_mpz_struct temp;
2326    
2327     //Eyeball the inputs carefully.
2328     assert(w != NULL);
2329     assert(w->size >= 0);
2330     assert(w->limbs != NULL);
2331     assert(u != NULL);
2332     assert(u->size >= 0);
2333     assert(u->limbs != NULL);
2334    
2335     //Create a temporary integer.
2336     GMP_INTS_mpz_init(&temp);
2337    
2338     //Set the temporary integer to the value of the input
2339     //argument.
2340     GMP_INTS_mpz_set_ui(&temp, v);
2341    
2342     //Do the actual multiplication. This recursive definition
2343     //is inherently wasteful, but I'm after clarity, not
2344     //warp speed.
2345     GMP_INTS_mpz_mul(w, u, &temp);
2346    
2347     //Destroy the temporary integer (this will reclaim the
2348     //memory).
2349     GMP_INTS_mpz_clear(&temp);
2350     }
2351    
2352    
2353     void GMP_INTS_mpz_tdiv_qr ( GMP_INTS_mpz_struct *quot,
2354     GMP_INTS_mpz_struct *rem,
2355     const GMP_INTS_mpz_struct *num,
2356     const GMP_INTS_mpz_struct *den)
2357     {
2358     GMP_INTS_size_t abs_num_size,
2359     abs_den_size,
2360     quotient_sign,
2361     remainder_sign,
2362     numerator_bitsize,
2363     denominator_bitsize,
2364     division_loop_count,
2365     division_loop_count_mod_32,
2366     division_loop_count_div_32,
2367     division_counter,
2368     i;
2369     GMP_INTS_limb_t temp_limb;
2370     GMP_INTS_limb_ptr trial_divisor;
2371    
2372     //Eyeball the input parameters.
2373     assert(quot != NULL);
2374     assert(quot->n_allocd > 0);
2375     assert(quot->limbs != NULL);
2376     assert(rem != NULL);
2377     assert(rem->n_allocd > 0);
2378     assert(rem->limbs != NULL);
2379     assert(num != NULL);
2380     assert(num->n_allocd > 0);
2381     assert(num->limbs != NULL);
2382     assert(den != NULL);
2383     assert(den->n_allocd > 0);
2384     assert(den->limbs != NULL);
2385    
2386     //We require for this function that the numerator, denominator, quotient, and
2387     //remainder all be distinct.
2388     assert(quot != rem);
2389     assert(quot != num);
2390     assert(quot != den);
2391     assert(rem != num);
2392     assert(rem != den);
2393     assert(num != den);
2394    
2395     //The GNU code was probably very efficient, but exceeded
2396     //my abilities to analyze. This is the classic
2397     //division algorithm.
2398    
2399     //First, start off with the quotient and remainder having
2400     //no error flags set. These will be set if appropriate.
2401     quot->flags = 0;
2402     rem->flags = 0;
2403    
2404     //First, handle tainted inputs. If the numerator or denominator
2405     //are bad or tainted, the quotient and remainder get tainted
2406     //automatically.
2407     {
2408     int taint;
2409    
2410     taint = GMP_INTS_two_op_flags_map(num->flags, den->flags);
2411    
2412     if (taint)
2413     {
2414     quot->flags = taint;
2415     rem->flags = taint;
2416     return;
2417     }
2418     }
2419    
2420     //The second possible cause for taint is if the divisor is
2421     //zero. This will get both the value of positive overflow.
2422     if (den->size == 0)
2423     {
2424     quot->flags = GMP_INTS_EF_INTOVF_POS;
2425     rem->flags = GMP_INTS_EF_INTOVF_POS;
2426     return;
2427     }
2428    
2429     //Handle the special case of a numerator of zero. If the numerator
2430     //is zero, the quotient and remainder are zero automatically.
2431     if (num->size == 0)
2432     {
2433     GMP_INTS_mpz_set_ui(quot, 0);
2434     GMP_INTS_mpz_set_ui(rem, 0);
2435     return;
2436     }
2437    
2438     //Generally, nothing else can go wrong as far as taint. The
2439     //value of the quotient is confined to be no larger than the
2440     //numerator, and the value of the remainder is confined to
2441     //be no larger than denominator-1. So, generally, if the
2442     //inputs are in size bounds, the outputs will be also.
2443    
2444     //Figure out how large in limbs the numerator and denominator actually
2445     //are.
2446     abs_num_size = GMP_INTS_abs_of_size_t(num->size);
2447     abs_den_size = GMP_INTS_abs_of_size_t(den->size);
2448    
2449     //Figure out the sign of things. We want the following relationship
2450     //to be true:
2451     // num/den = quot + rem/den.
2452     //The way to achieve this is to assign the sign of the quotient in the traditional
2453     //way, then to assign the remainder to have the same sign as the numerator.
2454     quotient_sign = num->size ^ den->size;
2455     remainder_sign = num->size;
2456    
2457     //The remainder starts off with the absolute value of the numerator, and then
2458     //we subtract from it as part of the division loop.
2459     GMP_INTS_mpz_copy(rem, num);
2460     //We know after the copy that the amount of space allocated in the remainder
2461     //MUST be at least as large as the absolute value of the numerator. So from
2462     //this point forward we use the space.
2463     assert(rem->n_allocd >= abs_num_size);
2464    
2465     //Figure out the number of significant bits in the numerator and denominator.
2466     //This determines the loop count over which we do the shift division loop.
2467     numerator_bitsize = GMP_INTS_BITS_PER_LIMB * abs_num_size;
2468    
2469     i = abs_num_size - 1;
2470    
2471     //We need to be extra careful here. One failure mode is that an integer
2472     //data structure is corrupted and the "size" field reflects limbs
2473     //that are zero. Need to watch that this kind of failure doesn't
2474     //cause memory access errors.
2475     assert(num->limbs[i] != 0);
2476     if (num->limbs[i] == 0)
2477     {
2478     quot->flags = GMP_INTS_EF_INTOVF_POS;
2479     rem->flags = GMP_INTS_EF_INTOVF_POS;
2480     return;
2481     }
2482    
2483     temp_limb = 0x80000000;
2484    
2485     while (((num->limbs[i] & temp_limb) == 0) && (temp_limb != 0))
2486     {
2487     numerator_bitsize--;
2488     temp_limb >>= 1;
2489     }
2490    
2491     denominator_bitsize = GMP_INTS_BITS_PER_LIMB * abs_den_size;
2492    
2493     i = abs_den_size - 1;
2494    
2495     //We need to be extra careful here. One failure mode is that an integer
2496     //data structure is corrupted and the "size" field reflects limbs
2497     //that are zero. Need to watch that this kind of failure doesn't
2498     //cause memory access errors.
2499     assert(den->limbs[i] != 0);
2500     if (den->limbs[i] == 0)
2501     {
2502     quot->flags = GMP_INTS_EF_INTOVF_POS;
2503     rem->flags = GMP_INTS_EF_INTOVF_POS;
2504     return;
2505     }
2506    
2507     temp_limb = 0x80000000;
2508    
2509     while (((den->limbs[i] & temp_limb) == 0) && (temp_limb != 0))
2510     {
2511     denominator_bitsize--;
2512     temp_limb >>= 1;
2513     }
2514    
2515     //The quotient starts off with the value of zero, but we consistently may
2516     //mask 1 into it and shift left. We need to be sure that we have as much
2517     //shift space there as is in the numerator. For this purpose we need to
2518     //prepare a block of clear memory as large as the numerator's.
2519     if (quot->n_allocd < abs_num_size)
2520     {
2521     GMP_INTS_mpz_realloc(quot, abs_num_size); //Make it big enough.
2522     }
2523     //Now, zero the memory.
2524     for (i=0; i<abs_num_size; i++)
2525     quot->limbs[i] = 0;
2526    
2527     //Determine the division loop count. This is the difference
2528     //in bit sizes between the numerator and denominator. It is
2529     //possible for this number to be negative, which means that the
2530     //main division loop will be executed zero times. This gives the
2531     //right results.
2532     division_loop_count = numerator_bitsize - denominator_bitsize;
2533    
2534     //We need to calculate some important numbers from the division loop
2535     //count. We need to know this number MOD 32 (which tells how far to
2536     //shift the divisor bitwise to line up with the numerator), and we
2537     //also need this number DIV 32 for the limb-wise shift.
2538     division_loop_count_mod_32 = division_loop_count % 32;
2539     division_loop_count_div_32 = division_loop_count / 32;
2540    
2541     //We now need a shift register in which we shift the denominator up
2542     //for repeated comparisons. We should dynamically allocate this to
2543     //be the same size as the numerator. Using _alloca() is OK, as one
2544     //of the unit tests is to be sure that _alloca() will handle integer
2545     //of the maximum allowed size.
2546     trial_divisor = _alloca(abs_num_size * sizeof(GMP_INTS_limb_t));
2547    
2548     //Our trial divisor needs to start off with the divisor shifted up
2549     //so that the most significant bit is aligned with the numerator.
2550     for (i = 0; i < abs_num_size; i++)
2551     {
2552     if ((division_loop_count < 0) || (i < division_loop_count_div_32))
2553     {
2554     trial_divisor[i] = 0;
2555     }
2556     else
2557     {
2558     if ((i-division_loop_count_div_32) < abs_den_size)
2559     trial_divisor[i] = den->limbs[i - division_loop_count_div_32];
2560     else
2561     trial_divisor[i] = 0;
2562     }
2563     }
2564    
2565     //The code above planted the limbs in the right place. Now need to shift bits
2566     //upward by the remaining number.
2567     if ((division_loop_count > 0) && (division_loop_count_mod_32 > 0))
2568     {
2569     //There is an existing function we can call to do the left shift.
2570     GMP_INTS_mpn_lshift(trial_divisor,
2571     trial_divisor,
2572     abs_num_size,
2573     division_loop_count_mod_32);
2574     }
2575    
2576    
2577     //Everything is ready to go. Now begin the division loop itself. It is possible
2578     //for the loop to execute zero times, which will happen if the denominator is longer
2579     //in bits than the numerator. In such cases, we can't execute this loop even once
2580     //because the math assumes that the numerator is at least as long as the denominator.
2581     for (division_counter = 0; division_counter < division_loop_count+1; division_counter++)
2582     {
2583     //Shift the quotient left one bit.
2584     GMP_INTS_mpn_lshift(quot->limbs,
2585     quot->limbs,
2586     abs_num_size,
2587     1);
2588    
2589     //If the remainder is at least as large as the trial divisor, subtract the trial
2590     //divisor from the remainder and mask in the quotient.
2591     if (GMP_INTS_mpn_cmp(rem->limbs,
2592     trial_divisor,
2593     abs_num_size) >= 0)
2594     {
2595     GMP_INTS_mpn_sub(rem->limbs,
2596     rem->limbs,
2597     abs_num_size,
2598     trial_divisor,
2599     abs_num_size);
2600     quot->limbs[0] |= 1;
2601     }
2602    
2603     //Shift the trial divisor right one bit.
2604     GMP_INTS_mpn_rshift(trial_divisor,
2605     trial_divisor,
2606     abs_num_size,
2607     1);
2608     } //End for each iteration of the division loop.
2609    
2610     //Normalize the quotient and the remainder. The normalization
2611     //process is to bring the sizes down if we have leading
2612     //zeros.
2613     quot->size = abs_num_size;
2614     GMP_INTS_mpn_normalize(quot->limbs, &(quot->size));
2615     rem->size = abs_num_size;
2616     GMP_INTS_mpn_normalize(rem->limbs, &(rem->size));
2617    
2618     //Adjust the signs as required.
2619     if (quotient_sign < 0)
2620     quot->size = -(quot->size);
2621     if (remainder_sign < 0)
2622     rem->size = -(rem->size);
2623     }
2624    
2625    
2626     void GMP_INTS_mpz_fac_ui(GMP_INTS_mpz_struct *result,
2627     unsigned long int n)
2628     {
2629     //Just multiply the numbers in ascending order. The original
2630     //GNU library contained a much more elegant algorithm, but
2631     //this is more direct.
2632    
2633     unsigned long int k;
2634    
2635     GMP_INTS_mpz_set_ui (result, 1L);
2636    
2637     for (k = 2; (k <= n) && !(result->flags); k++)
2638     GMP_INTS_mpz_mul_ui (result, result, k);
2639     }
2640    
2641    
2642     /******************************************************************/
2643     /*** PUBLIC CONVERSION AND OUTPUT FUNCTIONS ********************/
2644     /******************************************************************/
2645     //07/18/01: Visual inspection OK. Function returns
2646     //reasonable values even out to 100,000 digits--seems OK.
2647     int GMP_INTS_mpz_size_in_base_10(const GMP_INTS_mpz_struct *arg)
2648     {
2649     _int64 n;
2650    
2651     //Eyeball the input parameter.
2652     assert(arg != NULL);
2653     assert(arg->n_allocd > 0);
2654     assert(arg->limbs != NULL);
2655    
2656     //Get the number of limbs occupied by the integer.
2657     //Because even the digit zero takes some space,
2658     //don't accept zero for an answer.
2659     n = GMP_INTS_abs_of_size_t(arg->size);
2660     if (n==0)
2661     n = 1;
2662    
2663     //Convert this to the number of bits. Generously
2664     //ignore any unused leading bits.
2665     n *= 32;
2666    
2667     //Used a slightly high best rational approximation in F_{65535}
2668     //to go from the number of bits to the number of
2669     //digits. The division discards, so bump the result
2670     //up by 1 to compensate for possible truncation. The number
2671     //we are aproximating is ln(2)/ln(10).
2672     n *= 12655;
2673     n /= 42039;
2674     n++;
2675    
2676     //Compensate for possible commas in the result. Again,
2677     //consider truncation.
2678     n *= 4;
2679     n /= 3;
2680     n++;
2681    
2682     //Compensate for the minus sign, the trailing zero,
2683     //cosmic rays striking the computer from the martian
2684     //listening post camoflaged on the moon, and the
2685     //possibility that we might need to put text in the
2686     //string if any flag is set.
2687     n += 100;
2688    
2689     //And that should be a good return value.
2690     return((int) n);
2691     }
2692    
2693    
2694     //07/19/01: Visual inspection and unit test is OK.
2695     void GMP_INTS_mpz_to_string(char *out,
2696     const GMP_INTS_mpz_struct *in)
2697     {
2698     //Eyeball the input parameters.
2699     assert(out != NULL);
2700     assert(in != NULL);
2701     assert(in->n_allocd > 0);
2702     assert(in->limbs != NULL);
2703    
2704     //If any of the flags are set, stuff in the text.
2705     if (in->flags)
2706     {
2707     if (in->flags & GMP_INTS_EF_INTOVF_POS)
2708     {
2709     strcpy(out, GMP_INTS_EF_INTOVF_POS_STRING);
2710     }
2711     else if (in->flags & GMP_INTS_EF_INTOVF_NEG)
2712     {
2713     strcpy(out, GMP_INTS_EF_INTOVF_NEG_STRING);
2714     }
2715     else if (in->flags & GMP_INTS_EF_INTOVF_TAINT_POS)
2716     {
2717     strcpy(out, GMP_INTS_EF_INTOVF_TAINT_POS_STRING);
2718     }
2719     else if (in->flags & GMP_INTS_EF_INTOVF_TAINT_NEG)
2720     {
2721     strcpy(out, GMP_INTS_EF_INTOVF_TAINT_NEG_STRING);
2722     }
2723     else
2724     {
2725     strcpy(out, "INTERNAL_ERROR");
2726     }
2727     }
2728     else
2729     {
2730     //Ordinary integer conversion.
2731     GMP_INTS_mpz_struct num, den, quot, rem, k10;
2732    
2733     //Allocate space for the temporary integers.
2734     GMP_INTS_mpz_init(&num);
2735     GMP_INTS_mpz_init(&den);
2736     GMP_INTS_mpz_init(&quot);
2737     GMP_INTS_mpz_init(&rem);
2738     GMP_INTS_mpz_init(&k10);
2739    
2740     //Assign the constant 10.
2741     GMP_INTS_mpz_set_ui(&k10, 10);
2742    
2743     //If the integer is zero, assign that.
2744     if (in->size == 0)
2745     {
2746     strcpy(out, "0");
2747     }
2748     else
2749     {
2750     //We have to do a full conversion. The algorithm
2751     //is division by 10, each time obtaining the least
2752     //significant digit, until finally the quotient is
2753     //zero.
2754     char *ptr;
2755    
2756     ptr = out;
2757    
2758     GMP_INTS_mpz_copy(&num, in);
2759     GMP_INTS_mpz_copy(&den, &k10);
2760     do
2761     {
2762     #if 0
2763     printf("Values before division:\n");
2764     FCMIOF_hline();
2765     GMP_INTS_mpz_print_int(stdout, &num, "Numerator");
2766     FCMIOF_hline();
2767     GMP_INTS_mpz_print_int(stdout, &den, "Denominator");
2768     FCMIOF_hline();
2769     GMP_INTS_mpz_print_int(stdout, &quot, "Quotient");
2770     FCMIOF_hline();
2771     GMP_INTS_mpz_print_int(stdout, &rem, "Remainder");
2772     FCMIOF_hline();
2773    
2774     if (num.size > 1)
2775     FCMIOF_hline();
2776     #endif
2777    
2778     GMP_INTS_mpz_tdiv_qr(&quot, &rem, &num, &den);
2779     #if 0
2780     printf("Values after division:\n");
2781     FCMIOF_hline();
2782     GMP_INTS_mpz_print_int(stdout, &num, "Numerator");
2783     FCMIOF_hline();
2784     GMP_INTS_mpz_print_int(stdout, &den, "Denominator");
2785     FCMIOF_hline();
2786     GMP_INTS_mpz_print_int(stdout, &quot, "Quotient");
2787     FCMIOF_hline();
2788     GMP_INTS_mpz_print_int(stdout, &rem, "Remainder");
2789     FCMIOF_hline();
2790     #endif
2791    
2792     if (rem.size != 0)
2793     {
2794     *ptr = '0' + (char)(rem.limbs[0]);
2795     }
2796     else
2797     {
2798     *ptr = '0';
2799     }
2800     ptr++;
2801     GMP_INTS_mpz_copy(&num, &quot);
2802     //printf("digit\n");
2803     }
2804     while (!GMP_INTS_mpz_is_zero(&quot));
2805    
2806     //Finally, if the input was negative, tack on the
2807     //minus sign.
2808     if (GMP_INTS_mpz_is_neg(in))
2809     {
2810     *ptr = '-';
2811     ptr++;
2812     }
2813    
2814     //Finally, tack on the trailing zero terminator.
2815     *ptr = 0;
2816     ptr++;
2817    
2818     //Reverse the string.
2819     BSTRFUNC_str_reverse(out);
2820     }
2821    
2822     //Deallocate the integers.
2823     GMP_INTS_mpz_clear(&num);
2824     GMP_INTS_mpz_clear(&den);
2825     GMP_INTS_mpz_clear(&quot);
2826     GMP_INTS_mpz_clear(&rem);
2827     GMP_INTS_mpz_clear(&k10);
2828     }
2829     }
2830    
2831    
2832     void GMP_INTS_mpz_long_int_format_to_stream(FILE *s,
2833     const GMP_INTS_mpz_struct *i,
2834     const char *desc)
2835     {
2836     int line_len;
2837     int digits_per_line;
2838     char *digits;
2839     int num_digits;
2840     int nlines;
2841     int cur_line;
2842     int number_desc_width;
2843    
2844     //Eyeball the inputs, make sure the caller isn't doing
2845     //something stupid.
2846     assert(s != NULL);
2847     assert(i != NULL);
2848     assert(i->n_allocd > 0);
2849     assert(i->limbs != NULL);
2850     assert(desc != NULL);
2851    
2852     //Obtain the line length assumed for formatted output.
2853     line_len = FCMIOF_get_line_len();
2854    
2855     //The description width allowed is 20.
2856     number_desc_width = 20;
2857    
2858     /* The number of digits per line that we assume must be a multiple of
2859     ** three. The formula below was not examined very carefully, but it
2860     ** works fine for a line length of 78. If line length is changed,
2861     ** this formula may need to be examined very carefully and rewritten.
2862     */
2863     digits_per_line = INTFUNC_max(3, ((((line_len-42)*3)/4)/3)*3);
2864     assert(digits_per_line >= 3);
2865    
2866     /* We now need to get a digit string corresponding to this
2867     ** number. First, need to figure out how much and
2868     ** allocate the space.
2869     */
2870     digits = GMP_INTS_malloc(GMP_INTS_mpz_size_in_base_10(i) * sizeof(char));
2871     GMP_INTS_mpz_to_string(digits, i);
2872    
2873     //If the number is negative, delete the leading minus sign.
2874     //The rest of the display algorithm needs an unsigned
2875     //series of digits.
2876     if (*digits == '-')
2877     {
2878     int i = 0;
2879    
2880     do
2881     {
2882     digits[i] = digits[i+1];
2883     i++;
2884     }
2885     while(digits[i-1]);
2886     }
2887    
2888     //Figure out how many digits in the string representation.
2889     num_digits = strlen(digits);
2890    
2891     /* As the first order of business, figure out how many lines the beast
2892     ** will require.
2893     */
2894     if (i->flags)
2895     {
2896     nlines = 1; /* Only one line required for NAN verbeage. */
2897     }
2898     else if (GMP_INTS_mpz_is_zero(i))
2899     {
2900     nlines = 1; /* The zero value requires one line. */
2901     }
2902     else
2903     {
2904     /* In any other case, have a formula.
2905     */
2906     nlines = 1 + (num_digits - 1) / digits_per_line;
2907     }
2908    
2909     /* Iterate through each line, spitting out whatever is appropriate. */
2910     for (cur_line = 0; cur_line < nlines; cur_line++)
2911     {
2912     int cur_digit_on_line;
2913    
2914     /* If this is the first line, spit out the description, right-aligned.
2915     ** Otherwise, spit spaces.
2916     */
2917     if (!cur_line)
2918     {
2919     /* First line. */
2920     int len;
2921    
2922     len = strlen(desc);
2923    
2924     if (len <= number_desc_width)
2925     {
2926     /* Description is shorter or equal, pad on left. */
2927     FCMIOF_stream_repchar(s, ' ', number_desc_width - len);
2928     fprintf(s, "%s", desc);
2929     }
2930     else
2931     {
2932     /* Description is too long, truncate. */
2933     int i;
2934    
2935     for (i=0; i<number_desc_width; i++)
2936     fprintf(s, "%c", desc[i]);
2937     }
2938    
2939     fprintf(s, ": ");
2940    
2941     /* If the number is negative, throw in a minus sign. */
2942     if (GMP_INTS_mpz_is_neg(i) && !(i->flags))
2943     {
2944     fprintf(s, "- ");
2945     }
2946     else
2947     {
2948     fprintf(s, " ");
2949     }
2950     }
2951     else
2952     {
2953     /* Every line but first line. */
2954     FCMIOF_stream_repchar(s, ' ', number_desc_width+4);
2955     }
2956    
2957     for(cur_digit_on_line=0; cur_digit_on_line < digits_per_line; cur_digit_on_line++)
2958     {
2959     int idx_into_string;
2960     /* Index into the string which is our digit of interest.
2961     */
2962    
2963     /* Compute the index. The equation is based on the ordering
2964     ** of presentation, for example,
2965     **
2966     ** 7 6
2967     ** 5 4 3
2968     ** 2 1 0.
2969     **
2970     ** With a little thought, the equation should make sense.
2971     ** The index won't always be used to index into the string.
2972     */
2973     idx_into_string =
2974     ((((nlines-1) - cur_line) * digits_per_line)
2975     +
2976     (digits_per_line - 1 - cur_digit_on_line));
2977    
2978     /* Print the appropriate digit or a space. The NAN case and the
2979     ** zero case need to be treated specially.
2980     */
2981     if (i->flags)
2982     {
2983     /* Not a number. Everything is blank, except spell out
2984     ** description of condition at the end of the string of
2985     ** digits.
2986     */
2987     int index_from_right;
2988     int virtual_index;
2989    
2990     index_from_right = digits_per_line - 1 - cur_digit_on_line;
2991     //The index calculated above is calculated so that the
2992     //final position on the line has index [0].
2993     assert(index_from_right >= 0 && index_from_right < digits_per_line);
2994    
2995     //Now, calculate the "virtual index". The virtual index
2996     //is the actual number of characters from the right, taking
2997     //into account commas.
2998     virtual_index = index_from_right + index_from_right/3;
2999    
3000     if (((index_from_right % 3) == 2) && cur_digit_on_line)
3001     {
3002     //We are one position past a comma. This means
3003     //that we might need a "fill" character to go
3004     //where the comma should have gone.
3005    
3006     if (virtual_index + 1 < num_digits)
3007     {
3008     //The character we should print exists.
3009     fprintf(s, "%c", digits[num_digits - 2 - virtual_index]);
3010     }
3011     else
3012     {
3013     //The character doesn't exist, because the error
3014     //string is apparently too short. Must print a
3015     //space, instead.
3016     fprintf(s, " ");
3017     }
3018     }
3019    
3020     //We've done the fill character, if the position we're in
3021     //is one past a comma. Now, do the ordinary character
3022     //corresponding to a digit position.
3023     if (virtual_index < num_digits)
3024     {
3025     //The character we should print exists.
3026     fprintf(s, "%c", digits[num_digits - 1 - virtual_index]);
3027     }
3028     else
3029     {
3030     //The character doesn't exist, because the error
3031     //string is apparently too short. Must print a
3032     //space, instead.
3033     fprintf(s, " ");
3034     }
3035     }
3036     else if (GMP_INTS_mpz_is_zero(i))
3037     {
3038     /* This is the zero case. For zero, there is only one line,
3039     ** and every character except the last one is a blank.
3040     */
3041     if (cur_digit_on_line == (digits_per_line - 1))
3042     {
3043     fprintf(s, "0");
3044     }
3045     else
3046     {
3047     fprintf(s, " ");
3048     }
3049     }
3050     else
3051     {
3052     /* This is a valid number which is not zero. Need to print
3053     ** the digits.
3054     */
3055    
3056     if (idx_into_string < num_digits)
3057     {
3058     int actual_index;
3059    
3060     actual_index = num_digits - 1 - idx_into_string;
3061     //This is a string reversal mapping. The original
3062     //code stored strings least significant digit first,
3063     //but this code uses most significant digit first.
3064     assert((actual_index >= 0) && (actual_index < num_digits));
3065     fprintf(s, "%c", digits[actual_index]);
3066     }
3067     else
3068     {
3069     fprintf(s, " ");
3070     }
3071     } /* End of digit case.
3072    
3073     /* Now handle the commas. The rules for commas are straightforward.
3074     ** a)NAN never has a comma.
3075     ** b)Zeros never have a comma.
3076     ** c)The final line, last digit never has a comma.
3077     ** d)Everything else in multiples of three ...
3078     */
3079     if (!(idx_into_string % 3) && (idx_into_string))
3080     {
3081     if (i->flags)
3082     {
3083     //fprintf(s, " ");
3084     }
3085     else if (!num_digits)
3086     {
3087     fprintf(s, " ");
3088     }
3089     else
3090     {
3091     if (idx_into_string < num_digits)
3092     {
3093     fprintf(s, ",");
3094     }
3095     else
3096     {
3097     fprintf(s, " ");
3098     }
3099     }
3100     }
3101     } /* End for each digit on the current line. */
3102    
3103     /* For the first line, print out an informative message
3104     ** advising of the number of digits. For all other lines
3105     ** print nothing.
3106     */
3107     if (!cur_line && !(i->flags))
3108     {
3109     if (nlines == 1)
3110     fprintf(s, " ");
3111    
3112     if (num_digits <= 1)
3113     {
3114     fprintf(s, " ( 1 digit )\n");
3115     }
3116     else if (num_digits < 1000)
3117     {
3118     fprintf(s, " (%7d digits)\n", num_digits);
3119     }
3120     else
3121     {
3122     fprintf(s, " (%3d,%03d digits)\n", num_digits / 1000, num_digits % 1000);
3123     }
3124     }
3125     else
3126     {
3127     fprintf(s, "\n");
3128     }
3129     } /* End for each line. */
3130    
3131     //Deallocate the string space.
3132     GMP_INTS_free(digits);
3133     }
3134    
3135    
3136     void GMP_INTS_mpz_arb_int_raw_to_stream(FILE *s,
3137     const GMP_INTS_mpz_struct *i)
3138     {
3139     int size_reqd;
3140     char *digits;
3141    
3142     //Eyeball the input parameters.
3143     assert(s != NULL);
3144     assert(i != NULL);
3145     assert(i->n_allocd > 0);
3146     assert(i->limbs != NULL);
3147    
3148     size_reqd = GMP_INTS_mpz_size_in_base_10(i);
3149     digits = GMP_INTS_malloc(size_reqd * sizeof(char));
3150     GMP_INTS_mpz_to_string(digits, i);
3151     fprintf(s, "%s", digits);
3152     GMP_INTS_free(digits);
3153     }
3154    
3155    
3156     //07/24/01: Passed visual inspection and unit tests.
3157     void GMP_INTS_mpz_pow_ui( GMP_INTS_mpz_struct *result,
3158     const GMP_INTS_mpz_struct *base,
3159     unsigned exponent)
3160     {
3161     GMP_INTS_mpz_struct temp;
3162     //Temporary location to hold the base raised to
3163     //a binary power (repeated squaring).
3164    
3165     //Eyeball the input parameters.
3166     assert(result != NULL);
3167     assert(result->n_allocd > 0);
3168     assert(result->limbs != NULL);
3169     assert(base != NULL);
3170     assert(base->n_allocd > 0);
3171     assert(base->limbs != NULL);
3172    
3173     //For this function, the base and the result may not
3174     //be the same object.
3175     assert(result != base);
3176    
3177     //If the base is tained, the output is tainted by association.
3178     {
3179     int taint;
3180    
3181     taint = GMP_INTS_two_op_flags_map(base->flags, 0);
3182    
3183     if (taint)
3184     {
3185     result->flags = taint;
3186     return;
3187     }
3188     }
3189    
3190     //Allocate our temporary variable and set it to the base.
3191     GMP_INTS_mpz_init(&temp);
3192     GMP_INTS_mpz_copy(&temp, base);
3193    
3194     //The result begins with the value of 1.
3195     GMP_INTS_mpz_set_ui(result, 1);
3196    
3197     //Loop through, processing each bit of the exponent. This is a fairly effective
3198     //algorithm, but not the optimal one (Knuth points this out).
3199     while (exponent && !result->flags)
3200     {
3201     if (exponent & 0x1)
3202     {
3203     GMP_INTS_mpz_mul(result, result, &temp);
3204     }
3205    
3206     //Square the temporary variable. Because squaring of arb integer
3207     //may be very expensive, the test against 1 (i.e. last iteration)
3208     //certainly pays for itself.
3209     if (exponent != 1)
3210     GMP_INTS_mpz_mul(&temp, &temp, &temp);
3211    
3212     exponent >>= 1;
3213     }
3214    
3215     //Deallocate our temporary variable.
3216     GMP_INTS_mpz_clear(&temp);
3217     }
3218    
3219    
3220     void GMP_INTS_mpz_abs(GMP_INTS_mpz_struct *arg)
3221     {
3222     //Eyeball the input parameter.
3223     assert(arg != NULL);
3224     assert(arg->n_allocd > 0);
3225     assert(arg->limbs != NULL);
3226    
3227     //Take the absolute value.
3228     if (arg->size < 0)
3229     arg->size = -arg->size;
3230     }
3231    
3232    
3233     //07/29/01: Visual inspection passed. Seems to work fine--not explicitly unit-tested
3234     //directly, but was tested from Tcl.
3235     void GMP_INTS_mpz_gcd(GMP_INTS_mpz_struct *result,
3236     const GMP_INTS_mpz_struct *arg1,
3237     const GMP_INTS_mpz_struct *arg2)
3238     {
3239     GMP_INTS_mpz_struct u, v, q, r;
3240     int loop_count;
3241    
3242     //Eyeball the inputs carefully.
3243     assert(result != NULL);
3244     assert(result->n_allocd > 0);
3245     assert(result->limbs != NULL);
3246     assert(arg1 != NULL);
3247     assert(arg1->n_allocd > 0);
3248     assert(arg1->limbs != NULL);
3249     assert(arg2 != NULL);
3250     assert(arg2->n_allocd > 0);
3251     assert(arg2->limbs != NULL);
3252    
3253     //Args are not allowed to be same object.
3254     assert(arg1 != arg2);
3255    
3256     //If either input is error or taint, taint the output.
3257     {
3258     int taint;
3259    
3260     taint = GMP_INTS_two_op_flags_map(arg1->flags, arg2->flags);
3261    
3262     result->flags = 0;
3263     //"result" starts off with a clean slate. Must do this
3264     //after taint calculation in case locations of arg1 or arg2
3265     //are the same as result.
3266     if (taint)
3267     {
3268     result->flags = taint;
3269     return;
3270     }
3271     }
3272    
3273     //If either input is zero, the result is 1.
3274     if (GMP_INTS_mpz_is_zero(arg1) || GMP_INTS_mpz_is_zero(arg2))
3275     {
3276     GMP_INTS_mpz_set_ui(result, 1);
3277     return;
3278     }
3279    
3280     //Allocate space for locals.
3281     GMP_INTS_mpz_init(&u);
3282     GMP_INTS_mpz_init(&v);
3283     GMP_INTS_mpz_init(&q);
3284     GMP_INTS_mpz_init(&r);
3285    
3286     //We are following Knuth Vol 2, p. 337, the modern Euclidian algorithm.
3287     //Note: There are faster algorithms for GCD, but because I hacked up the
3288     //GMP multiple-precision library so badly, those aren't included. This one
3289     //is logically correct but sub-optimal. Perhaps at a later time faster
3290     //algorithms will be re-included.
3291     //Copy inputs to u and v.
3292     GMP_INTS_mpz_copy(&u, arg1);
3293     GMP_INTS_mpz_copy(&v, arg2);
3294    
3295     //Take the absolute value of each argument. We know that neither is zero,
3296     //but one or both might be negative.
3297     GMP_INTS_mpz_abs(&u);
3298     GMP_INTS_mpz_abs(&v);
3299    
3300     //Begin Euclid's algorithm. There are really three possibilities:
3301     // a)We terminate normally.
3302     // b)Somehow we generate a math error and terminate based on flags.
3303     // c)Due to some unknown error in the math functions, we go on forever,
3304     // and the program locks up.
3305     GMP_INTS_mpz_set_ui(&r, 1);
3306     loop_count = 0;
3307     while (!GMP_INTS_mpz_is_zero(&r) && !q.flags && !r.flags && (loop_count < 100000))
3308     {
3309     loop_count++;
3310    
3311     GMP_INTS_mpz_tdiv_qr(&q, &r, &u, &v);
3312     GMP_INTS_mpz_copy(&u, &v);
3313     GMP_INTS_mpz_copy(&v, &r);
3314     }
3315    
3316     //Let's hope we didn't get out of the loop based on loop count.
3317     assert(loop_count != 100000);
3318    
3319     //u now contains the answer.
3320     GMP_INTS_mpz_copy(result, &u);
3321    
3322     //Deallocate space for locals.
3323     GMP_INTS_mpz_clear(&u);
3324     GMP_INTS_mpz_clear(&v);
3325     GMP_INTS_mpz_clear(&q);
3326     GMP_INTS_mpz_clear(&r);
3327     }
3328    
3329    
3330     /******************************************************************/
3331     /*** COMPARISON AND SIZING FUNCTIONS ***************************/
3332     /******************************************************************/
3333     //07/24/01: Visual inspection only, due to simplicity.
3334     int GMP_INTS_mpz_fits_uint_p (const GMP_INTS_mpz_struct *src)
3335     {
3336     GMP_INTS_size_t size;
3337     GMP_INTS_limb_t mpl;
3338    
3339     //Eyeball the input parameter.
3340     assert(src != NULL);
3341     assert(src->n_allocd > 0);
3342     assert(src->limbs != NULL);
3343    
3344     mpl = src->limbs[0];
3345     size = src->size;
3346     if (size < 0 || size > 1)
3347     return(0);
3348    
3349     //The following line came from the original GNU code.
3350     //It isn't necessary in our case since limbs and ints are
3351     //both 32 bits, but it will do no harm.
3352     return (mpl <= (~(unsigned int) 0));
3353     }
3354    
3355    
3356     unsigned GMP_INTS_mpz_get_limb_zero(const GMP_INTS_mpz_struct *src)
3357     {
3358     //Eyeball the inputs.
3359     assert(src != NULL);
3360     assert(src->n_allocd > 0);
3361     assert(src->limbs != NULL);
3362    
3363     if (!src->size)
3364     return(0);
3365     else
3366     return(src->limbs[0]);
3367     }
3368    
3369    
3370     //07/24/01: Visual inspection only. Understood the comparisons
3371     //and seems like they should work, but ... a little beyond my
3372     //comfort zone without testing. Trusting GNU on this one ...
3373     int GMP_INTS_mpz_fits_sint_p (const GMP_INTS_mpz_struct *src)
3374     {
3375     GMP_INTS_size_t size;
3376     GMP_INTS_limb_t mpl;
3377    
3378     //Eyeball the input parameter.
3379     assert(src != NULL);
3380     assert(src->n_allocd > 0);
3381     assert(src->limbs != NULL);
3382    
3383     mpl = src->limbs[0];
3384     size = src->size;
3385     if (size > 0)
3386     {
3387     if (size > 1)
3388     return 0;
3389     return mpl < ~((~(unsigned int) 0) >> 1);
3390     }
3391     else
3392     {
3393     if (size < -1)
3394     return 0;
3395     return mpl <= ~((~(unsigned int) 0) >> 1);
3396     }
3397     }
3398    
3399    
3400     //07/24/01: Visual inspection only. One issue that caught
3401     //my eye is that in one place function returned neg value if
3402     //< and pos value if >. Was within spec, but corrected because
3403     //it concerned me as I often test against -1 and 1. Seems
3404     //to invite accidents.
3405     int GMP_INTS_mpz_cmp (const GMP_INTS_mpz_struct *u,
3406     const GMP_INTS_mpz_struct *v)
3407     {
3408     GMP_INTS_size_t usize = u->size;
3409     GMP_INTS_size_t vsize = v->size;
3410     GMP_INTS_size_t size;
3411     GMP_INTS_limb_srcptr up, vp;
3412     int cmp;
3413    
3414     //Eyeball the input parameters.
3415     assert(u != NULL);
3416     assert(u->n_allocd > 0);
3417     assert(u->limbs != NULL);
3418     assert(v != NULL);
3419     assert(v->n_allocd > 0);
3420     assert(v->limbs != NULL);
3421    
3422     if (usize < vsize)
3423     return(-1);
3424     else if (usize > vsize)
3425     return(1);
3426    
3427     if (usize == 0)
3428     return(0);
3429    
3430     size = GMP_INTS_abs_of_size_t(usize);
3431    
3432     up = u->limbs;
3433     vp = v->limbs;
3434    
3435     cmp = GMP_INTS_mpn_cmp (up, vp, size);
3436    
3437     if (cmp == 0)
3438     return(0);
3439    
3440     if ((cmp < 0) == (usize < 0))
3441     return(1);
3442     else
3443     return(-1);
3444     }
3445    
3446    
3447     //07/24/01: Not visually inspected. Relying on
3448     //GNU ...
3449     int GMP_INTS_mpz_cmp_ui (const GMP_INTS_mpz_struct *u,
3450     unsigned long int v_digit)
3451     {
3452     GMP_INTS_size_t usize = u->size;
3453    
3454     //Eyeball the input parameter.
3455     assert(u != NULL);
3456     assert(u->n_allocd > 0);
3457     assert(u->limbs != NULL);
3458    
3459     if (usize == 0)
3460     return -(v_digit != 0);
3461    
3462     if (usize == 1)
3463     {
3464     GMP_INTS_limb_t u_digit;
3465    
3466     u_digit = u->limbs[0];
3467     if (u_digit > v_digit)
3468     return 1;
3469     if (u_digit < v_digit)
3470     return -1;
3471     return 0;
3472     }
3473    
3474     return (usize > 0) ? 1 : -1;
3475     }
3476    
3477    
3478     //07/24/01: Not visually inspected. Relying on GNU.
3479     int GMP_INTS_mpz_cmp_si (const GMP_INTS_mpz_struct *u,
3480     signed long int v_digit)
3481     {
3482     GMP_INTS_size_t usize = u->size;
3483     GMP_INTS_size_t vsize;
3484     GMP_INTS_limb_t u_digit;
3485    
3486     //Eyeball the input parameter.
3487     assert(u != NULL);
3488     assert(u->n_allocd > 0);
3489     assert(u->limbs != NULL);
3490    
3491     vsize = 0;
3492     if (v_digit > 0)
3493     vsize = 1;
3494     else if (v_digit < 0)
3495     {
3496     vsize = -1;
3497     v_digit = -v_digit;
3498     }
3499    
3500     if (usize != vsize)
3501     return usize - vsize;
3502    
3503     if (usize == 0)
3504     return 0;
3505    
3506     u_digit = u->limbs[0];
3507    
3508     if (u_digit == (GMP_INTS_limb_t) (unsigned long) v_digit)
3509     return 0;
3510    
3511     if (u_digit > (GMP_INTS_limb_t) (unsigned long) v_digit)
3512     return usize;
3513     else
3514     return -usize;
3515     }
3516    
3517    
3518     //07/24/01: Not visually inspected. Counting on GNU.
3519     int GMP_INTS_mpz_cmpabs (const GMP_INTS_mpz_struct *u,
3520     const GMP_INTS_mpz_struct *v)
3521     {
3522     GMP_INTS_size_t usize = u->size;
3523     GMP_INTS_size_t vsize = v->size;
3524     GMP_INTS_limb_srcptr up, vp;
3525     int cmp;
3526    
3527     //Eyeball the input parameters.
3528     assert(u != NULL);
3529     assert(u->n_allocd > 0);
3530     assert(u->limbs != NULL);
3531     assert(v != NULL);
3532     assert(v->n_allocd > 0);
3533     assert(v->limbs != NULL);
3534    
3535     usize = GMP_INTS_abs_of_size_t(usize);
3536     vsize = GMP_INTS_abs_of_size_t(vsize);
3537    
3538     if (usize != vsize)
3539     return usize - vsize;
3540    
3541     if (usize == 0)
3542     return 0;
3543    
3544     up = u->limbs;
3545     vp = v->limbs;
3546    
3547     cmp = GMP_INTS_mpn_cmp (up, vp, usize);
3548    
3549     return cmp;
3550     }
3551    
3552     //07/24/01: Not visually inspected. Counting on GNU.
3553     int GMP_INTS_mpz_cmpabs_ui(const GMP_INTS_mpz_struct *u,
3554     unsigned long int v_digit)
3555     {
3556     GMP_INTS_size_t usize = u->size;
3557    
3558     //Eyeball the input parameter.
3559     assert(u != NULL);
3560     assert(u->n_allocd > 0);
3561     assert(u->limbs != NULL);
3562    
3563     if (usize == 0)
3564     return -(v_digit != 0);
3565    
3566     usize = GMP_INTS_abs_of_size_t(usize);
3567    
3568     if (usize == 1)
3569     {
3570     GMP_INTS_limb_t u_digit;
3571    
3572     u_digit = u->limbs[0];
3573     if (u_digit > v_digit)
3574     return 1;
3575     if (u_digit < v_digit)
3576     return -1;
3577     return 0;
3578     }
3579    
3580     return 1;
3581     }
3582    
3583    
3584     /******************************************************************/
3585     /*** VERSION CONTROL IDENTITY FUNCTIONS ************************/
3586     /******************************************************************/
3587    
3588     //07/18/01: Visual inspection only. Function deemed too
3589     //simple for unit testing.
3590     const char *GMP_INTS_cvcinfo(void)
3591     {
3592 dashley 56 return("$Header$");
3593 dashley 25 }
3594    
3595    
3596     //07/18/01: Visual inspection only. Function deemed too
3597     //simple for unit testing.
3598     const char *GMP_INTS_hvcinfo(void)
3599     {
3600     return(GMP_INTS_H_VERSION);
3601     }
3602    
3603 dashley 56 //End of gmp_ints.c.

Properties

Name Value
svn:keywords Header

dashley@gmail.com
ViewVC Help
Powered by ViewVC 1.1.25