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

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

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 56 by dashley, Sat Oct 29 01:53:01 2016 UTC revision 71 by dashley, Sat Nov 5 11:07:06 2016 UTC
# Line 1  Line 1 
1  //$Header$  //$Header$
2  //-------------------------------------------------------------------------------------------------  //-------------------------------------------------------------------------------------------------
3  //This file is part of "David T. Ashley's Shared Source Code", a set of shared components  //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.  //integrated into many of David T. Ashley's projects.
5  //-------------------------------------------------------------------------------------------------  //-------------------------------------------------------------------------------------------------
6  //This source code and any program in which it is compiled/used is provided under the MIT License,  //This source code and any program in which it is compiled/used is provided under the MIT License,
7  //reproduced below.  //reproduced below.
8  //-------------------------------------------------------------------------------------------------  //-------------------------------------------------------------------------------------------------
9  //Permission is hereby granted, free of charge, to any person obtaining a copy of  //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  //this software and associated documentation files(the "Software"), to deal in the
11  //Software without restriction, including without limitation the rights to use,  //Software without restriction, including without limitation the rights to use,
12  //copy, modify, merge, publish, distribute, sublicense, and / or sell copies of the  //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,  //Software, and to permit persons to whom the Software is furnished to do so,
14  //subject to the following conditions :  //subject to the following conditions :
15  //  //
16  //The above copyright notice and this permission notice shall be included in all  //The above copyright notice and this permission notice shall be included in all
17  //copies or substantial portions of the Software.  //copies or substantial portions of the Software.
18  //  //
19  //THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR  //THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
20  //IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,  //IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
21  //FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.IN NO EVENT SHALL THE  //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  //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,  //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  //OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
25  //SOFTWARE.  //SOFTWARE.
26  //-------------------------------------------------------------------------------------------------  //-------------------------------------------------------------------------------------------------
27  #define MODULE_GMP_INTS  #define MODULE_GMP_INTS
28    
29  #include <assert.h>  #include <assert.h>
30  #include <ctype.h>  #include <ctype.h>
31  #include <process.h>  #include <process.h>
32  #include <stddef.h>  #include <stddef.h>
33  #include <stdio.h>  #include <stdio.h>
34  #include <string.h>  #include <string.h>
35  #include <time.h>  #include <time.h>
36    
37    
38  /* Only included the guarded allocation header if we are compiling  /* Only included the guarded allocation header if we are compiling
39  ** a DOS console type application.  Other types of applications have  ** a DOS console type application.  Other types of applications have
40  ** other ways of protecting for out of memory.  Including the  ** other ways of protecting for out of memory.  Including the
41  ** header would do no harm in these cases, but do no good, either.  ** header would do no harm in these cases, but do no good, either.
42  */  */
43  #if defined(APP_TYPE_SIMPLE_DOS_CONSOLE)  #if defined(APP_TYPE_SIMPLE_DOS_CONSOLE)
44     #include "ccmalloc.h"     #include "ccmalloc.h"
45  #elif defined(APP_TYPE_IJUSCRIPTER_IJUCONSOLE)  #elif defined(APP_TYPE_IJUSCRIPTER_IJUCONSOLE)
46     #include "tclalloc.h"     #include "tclalloc.h"
47  #else  #else
48     /* Do nothing. */     /* Do nothing. */
49  #endif  #endif
50    
51  #include "bstrfunc.h"  #include "bstrfunc.h"
52  #include "charfunc.h"  #include "charfunc.h"
53  #include "fcmiof.h"  #include "fcmiof.h"
54  #include "gmp_ints.h"  #include "gmp_ints.h"
55  #include "intfunc.h"  #include "intfunc.h"
56    
57    
58  /******************************************************************/  /******************************************************************/
59  /***  CUSTOM ALLOCATION FUNCTIONS   *******************************/  /***  CUSTOM ALLOCATION FUNCTIONS   *******************************/
60  /******************************************************************/  /******************************************************************/
61  /* We need to decide here on how memory not on the stack will be  /* 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  ** allocated (i.e. what will become of the standard functions
63  ** like malloc, free, etc.).  This is dependent on the type of  ** like malloc, free, etc.).  This is dependent on the type of
64  ** application we're making.  The possible types are so far are:  ** application we're making.  The possible types are so far are:
65  **    APP_TYPE_SIMPLE_DOS_CONSOLE :  **    APP_TYPE_SIMPLE_DOS_CONSOLE :
66  **       Simple DOS console application.  **       Simple DOS console application.
67  **    APP_TYPE_IJUSCRIPTER_IJUCONSOLE:  **    APP_TYPE_IJUSCRIPTER_IJUCONSOLE:
68  **       The Tcl tool build.  **       The Tcl tool build.
69  **  **
70  ** The custom allocation functions here are a "portal" or "wrapper"  ** The custom allocation functions here are a "portal" or "wrapper"
71  ** for how the integer and rational number functions should  ** for how the integer and rational number functions should
72  ** get memory.  ** get memory.
73  **  **
74  ** The functions below are standard, except that the GNU MP team  ** The functions below are standard, except that the GNU MP team
75  ** built more generality into what allocation schemes could be  ** built more generality into what allocation schemes could be
76  ** used by including size information in some calls that don't  ** used by including size information in some calls that don't
77  ** normally get it.  That is why there are some extra calls below  ** normally get it.  That is why there are some extra calls below
78  ** where the information is discarded.  Other than that, these are  ** where the information is discarded.  Other than that, these are
79  ** standard allocation calls.  ** standard allocation calls.
80  */  */
81  //07/15/01:  Visual inspection only.  Function deemed too  //07/15/01:  Visual inspection only.  Function deemed too
82  //simple for unit testing.  //simple for unit testing.
83  void *GMP_INTS_malloc( size_t size )  void *GMP_INTS_malloc( size_t size )
84     {     {
85     #if defined(APP_TYPE_SIMPLE_DOS_CONSOLE)     #if defined(APP_TYPE_SIMPLE_DOS_CONSOLE)
86        return(CCMALLOC_malloc(size));        return(CCMALLOC_malloc(size));
87     #elif defined(APP_TYPE_IJUSCRIPTER_IJUCONSOLE)     #elif defined(APP_TYPE_IJUSCRIPTER_IJUCONSOLE)
88        return(TclpAlloc(size));        return(TclpAlloc(size));
89     #else     #else
90        return(malloc(size));        return(malloc(size));
91     #endif     #endif
92     }     }
93    
94    
95  //07/15/01:  Visual inspection only.  Function deemed too  //07/15/01:  Visual inspection only.  Function deemed too
96  //simple for unit testing.  //simple for unit testing.
97  void *GMP_INTS_calloc( size_t num, size_t size )  void *GMP_INTS_calloc( size_t num, size_t size )
98     {     {
99     #if defined(APP_TYPE_SIMPLE_DOS_CONSOLE)     #if defined(APP_TYPE_SIMPLE_DOS_CONSOLE)
100        return(CCMALLOC_calloc(num, size));        return(CCMALLOC_calloc(num, size));
101     #elif defined(APP_TYPE_IJUSCRIPTER_IJUCONSOLE)     #elif defined(APP_TYPE_IJUSCRIPTER_IJUCONSOLE)
102        return(TclpCalloc(num, size));        return(TclpCalloc(num, size));
103     #else     #else
104        return(calloc(num, size));        return(calloc(num, size));
105     #endif     #endif
106     }     }
107    
108    
109  //07/15/01:  Visual inspection only.  Function deemed too  //07/15/01:  Visual inspection only.  Function deemed too
110  //simple for unit testing.  //simple for unit testing.
111  void *GMP_INTS_realloc( void *memblock, size_t size )  void *GMP_INTS_realloc( void *memblock, size_t size )
112     {     {
113     #if defined(APP_TYPE_SIMPLE_DOS_CONSOLE)     #if defined(APP_TYPE_SIMPLE_DOS_CONSOLE)
114        return(CCMALLOC_realloc(memblock, size));        return(CCMALLOC_realloc(memblock, size));
115     #elif defined(APP_TYPE_IJUSCRIPTER_IJUCONSOLE)     #elif defined(APP_TYPE_IJUSCRIPTER_IJUCONSOLE)
116        return(TclpRealloc(memblock, size));        return(TclpRealloc(memblock, size));
117     #else     #else
118        return(realloc(memblock, size));        return(realloc(memblock, size));
119     #endif     #endif
120     }     }
121    
122    
123  //07/15/01:  Visual inspection only.  Function deemed too  //07/15/01:  Visual inspection only.  Function deemed too
124  //simple for unit testing.  //simple for unit testing.
125  void *GMP_INTS_realloc_w_size( void *memblock,  void *GMP_INTS_realloc_w_size( void *memblock,
126                                 size_t old_size,                                 size_t old_size,
127                                 size_t size )                                 size_t size )
128     {     {
129     #if defined(APP_TYPE_SIMPLE_DOS_CONSOLE)     #if defined(APP_TYPE_SIMPLE_DOS_CONSOLE)
130        return(CCMALLOC_realloc(memblock, size));        return(CCMALLOC_realloc(memblock, size));
131     #elif defined(APP_TYPE_IJUSCRIPTER_IJUCONSOLE)     #elif defined(APP_TYPE_IJUSCRIPTER_IJUCONSOLE)
132        return(TclpRealloc(memblock, size));        return(TclpRealloc(memblock, size));
133     #else     #else
134        return(realloc(memblock, size));        return(realloc(memblock, size));
135     #endif     #endif
136     }     }
137    
138    
139  //07/15/01:  Visual inspection only.  Function deemed too  //07/15/01:  Visual inspection only.  Function deemed too
140  //simple for unit testing.  //simple for unit testing.
141  void GMP_INTS_free( void *memblock )  void GMP_INTS_free( void *memblock )
142     {     {
143     #if defined(APP_TYPE_SIMPLE_DOS_CONSOLE)     #if defined(APP_TYPE_SIMPLE_DOS_CONSOLE)
144        CCMALLOC_free(memblock);        CCMALLOC_free(memblock);
145     #elif defined(APP_TYPE_IJUSCRIPTER_IJUCONSOLE)     #elif defined(APP_TYPE_IJUSCRIPTER_IJUCONSOLE)
146        TclpFree(memblock);        TclpFree(memblock);
147     #else     #else
148        free(memblock);        free(memblock);
149     #endif     #endif
150     }     }
151    
152    
153  //07/15/01:  Visual inspection only.  Function deemed too  //07/15/01:  Visual inspection only.  Function deemed too
154  //simple for unit testing.  //simple for unit testing.
155  void GMP_INTS_free_w_size( void *memblock, size_t size )  void GMP_INTS_free_w_size( void *memblock, size_t size )
156     {     {
157     #if defined(APP_TYPE_SIMPLE_DOS_CONSOLE)     #if defined(APP_TYPE_SIMPLE_DOS_CONSOLE)
158        CCMALLOC_free(memblock);        CCMALLOC_free(memblock);
159     #elif defined(APP_TYPE_IJUSCRIPTER_IJUCONSOLE)     #elif defined(APP_TYPE_IJUSCRIPTER_IJUCONSOLE)
160        TclpFree(memblock);        TclpFree(memblock);
161     #else     #else
162        free(memblock);        free(memblock);
163     #endif     #endif
164     }     }
165    
166    
167  /******************************************************************/  /******************************************************************/
168  /***  PORTABILITY CHECK FUNCTIONS   *******************************/  /***  PORTABILITY CHECK FUNCTIONS   *******************************/
169  /******************************************************************/  /******************************************************************/
170  //Because there is the risk that Microsoft Visual C++ might  //Because there is the risk that Microsoft Visual C++ might
171  //change in the future, the following function can be called  //change in the future, the following function can be called
172  //to see if the assumptions about data sizes are valid.  This  //to see if the assumptions about data sizes are valid.  This
173  //function returns TRUE if there is a problem, or FALSE  //function returns TRUE if there is a problem, or FALSE
174  //otherwise.  //otherwise.
175    
176  //07/15/01:  Unit testing complete.  //07/15/01:  Unit testing complete.
177  int GMP_INTS_data_sizes_are_wrong(void)  int GMP_INTS_data_sizes_are_wrong(void)
178     {     {
179     int i;     int i;
180     GMP_INTS_limb_t tv;     GMP_INTS_limb_t tv;
181     _int64 tv64;     _int64 tv64;
182    
183     //Check the number of bit rolls required to get the limb     //Check the number of bit rolls required to get the limb
184     //to go to zero again.  This had better be 32.     //to go to zero again.  This had better be 32.
185     tv = 1;     tv = 1;
186     i = 0;     i = 0;
187     while (tv)     while (tv)
188        {        {
189        tv <<= 1;        tv <<= 1;
190        i++;        i++;
191        }        }
192     if (i != 32)     if (i != 32)
193        return(1);        return(1);
194    
195     //Check that an _int64 is really and truly 64 bits.     //Check that an _int64 is really and truly 64 bits.
196     tv64 = 1;     tv64 = 1;
197     i = 0;     i = 0;
198     while (tv64)     while (tv64)
199        {        {
200        tv64 <<= 1;        tv64 <<= 1;
201        i++;        i++;
202        }        }
203     if (i != 64)     if (i != 64)
204        return(1);        return(1);
205    
206     //Room for additional tests here if needed later.     //Room for additional tests here if needed later.
207    
208     return(0);     return(0);
209     }     }
210    
211    
212  /******************************************************************/  /******************************************************************/
213  /***  ERROR STRING IDENTIFICATION AND PROCESSING FUNCTIONS  *******/  /***  ERROR STRING IDENTIFICATION AND PROCESSING FUNCTIONS  *******/
214  /******************************************************************/  /******************************************************************/
215    
216  int GMP_INTS_identify_nan_string(const char *s)  int GMP_INTS_identify_nan_string(const char *s)
217     {     {
218     assert(s != NULL);     assert(s != NULL);
219    
220     if (!strcmp(s, GMP_INTS_EF_INTOVF_POS_STRING))     if (!strcmp(s, GMP_INTS_EF_INTOVF_POS_STRING))
221        return(0);        return(0);
222     else if (!strcmp(s, GMP_INTS_EF_INTOVF_NEG_STRING))     else if (!strcmp(s, GMP_INTS_EF_INTOVF_NEG_STRING))
223        return(1);        return(1);
224     else if (!strcmp(s, GMP_INTS_EF_INTOVF_TAINT_POS_STRING))     else if (!strcmp(s, GMP_INTS_EF_INTOVF_TAINT_POS_STRING))
225        return(2);        return(2);
226     else if (!strcmp(s, GMP_INTS_EF_INTOVF_TAINT_NEG_STRING))     else if (!strcmp(s, GMP_INTS_EF_INTOVF_TAINT_NEG_STRING))
227        return(3);        return(3);
228     else     else
229        return(-1);        return(-1);
230     }     }
231    
232    
233  const char *GMP_INTS_supply_nan_string(int idx)  const char *GMP_INTS_supply_nan_string(int idx)
234     {     {
235     assert((idx >= 0) && (idx <= 3));     assert((idx >= 0) && (idx <= 3));
236    
237     if (idx==0)     if (idx==0)
238        return(GMP_INTS_EF_INTOVF_POS_STRING);        return(GMP_INTS_EF_INTOVF_POS_STRING);
239     else if (idx==1)     else if (idx==1)
240        return(GMP_INTS_EF_INTOVF_NEG_STRING);        return(GMP_INTS_EF_INTOVF_NEG_STRING);
241     else if (idx==2)     else if (idx==2)
242        return(GMP_INTS_EF_INTOVF_TAINT_POS_STRING);        return(GMP_INTS_EF_INTOVF_TAINT_POS_STRING);
243     else     else
244        return(GMP_INTS_EF_INTOVF_TAINT_NEG_STRING);        return(GMP_INTS_EF_INTOVF_TAINT_NEG_STRING);
245     }     }
246    
247    
248  /******************************************************************/  /******************************************************************/
249  /***  DEBUG PRINTING FUNCTIONS   **********************************/  /***  DEBUG PRINTING FUNCTIONS   **********************************/
250  /******************************************************************/  /******************************************************************/
251  //These functions are for printing out integers and limbs  //These functions are for printing out integers and limbs
252  //and groups of limbs for unit testing and debugging.                          //and groups of limbs for unit testing and debugging.                        
253    
254  //07/15/01:  Exempt from testing because debug/development  //07/15/01:  Exempt from testing because debug/development
255  //function.  //function.
256  void GMP_INTS_print_limb_group(FILE *stream,  void GMP_INTS_print_limb_group(FILE *stream,
257                                 GMP_INTS_limb_srcptr lg,                                 GMP_INTS_limb_srcptr lg,
258                                 GMP_INTS_size_t n,                                 GMP_INTS_size_t n,
259                                 char *desc)                                 char *desc)
260     {     {
261     int i;     int i;
262    
263     assert(stream != NULL);     assert(stream != NULL);
264     assert(n >= 0);     assert(n >= 0);
265     assert(desc != NULL);     assert(desc != NULL);
266    
267     if (!lg)     if (!lg)
268        {        {
269        fprintf(stream, "   %s: NULL\n", desc);        fprintf(stream, "   %s: NULL\n", desc);
270        }        }
271     else     else
272        {        {
273        for (i=n-1; i>=0; i--)        for (i=n-1; i>=0; i--)
274           {           {
275           fprintf(stream, "   %s[%2d]: 0x%8X\n", desc, i, lg[i]);           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,  void GMP_INTS_mpz_print_int(FILE *stream,
282                              const GMP_INTS_mpz_struct *arg,                              const GMP_INTS_mpz_struct *arg,
283                              char *desc)                              char *desc)
284     {     {
285     int i;     int i;
286    
287     assert(stream != NULL);     assert(stream != NULL);
288     assert(arg != NULL);     assert(arg != NULL);
289     assert(desc != NULL);     assert(desc != NULL);
290    
291     fprintf(stream, "Printing integer:\n   %s\n", desc);     fprintf(stream, "Printing integer:\n   %s\n", desc);
292    
293     fprintf(stream, "               flags: %d\n", arg->flags);     fprintf(stream, "               flags: %d\n", arg->flags);
294     fprintf(stream, "   ptr value to body: %p\n", arg);     fprintf(stream, "   ptr value to body: %p\n", arg);
295     fprintf(stream, "            n_allocd: %d\n", arg->n_allocd);     fprintf(stream, "            n_allocd: %d\n", arg->n_allocd);
296     fprintf(stream, "                size: %d\n", arg->size);     fprintf(stream, "                size: %d\n", arg->size);
297     fprintf(stream, "     limbs (ptr val): %p\n", arg->limbs);     fprintf(stream, "     limbs (ptr val): %p\n", arg->limbs);
298    
299     for (i=arg->n_allocd-1; i>=0; i--)     for (i=arg->n_allocd-1; i>=0; i--)
300        {        {
301        fprintf(stream, "   limb[%3d]: %8X\n", i, arg->limbs[i]);        fprintf(stream, "   limb[%3d]: %8X\n", i, arg->limbs[i]);
302        }        }
303     }     }
304    
305    
306  /******************************************************************/  /******************************************************************/
307  /***  LOW-LEVEL MACRO REPLACEMENTS   ******************************/  /***  LOW-LEVEL MACRO REPLACEMENTS   ******************************/
308  /******************************************************************/  /******************************************************************/
309  //The functions in this category are replacements for macros.  //The functions in this category are replacements for macros.
310  //Clarity was gained at the expense of speed.  //Clarity was gained at the expense of speed.
311    
312  int GMP_INTS_mpz_get_flags (const GMP_INTS_mpz_struct *arg)  int GMP_INTS_mpz_get_flags (const GMP_INTS_mpz_struct *arg)
313     {     {
314     assert(arg != NULL);     assert(arg != NULL);
315     assert(arg->n_allocd > 0);     assert(arg->n_allocd > 0);
316    
317     return(arg->flags);     return(arg->flags);
318     }     }
319    
320    
321  //07/15/01:  Visual inspection only.  Function deemed too  //07/15/01:  Visual inspection only.  Function deemed too
322  //simple for unit testing.  //simple for unit testing.
323  GMP_INTS_size_t GMP_INTS_abs_of_size_t(GMP_INTS_size_t arg)  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     //Be sure that the bit pattern does not represent the maximum
326     //negative argument.  Negating this would give the result of     //negative argument.  Negating this would give the result of
327     //zero, which is not what is intended.     //zero, which is not what is intended.
328     assert(arg != 0x80000000);     assert(arg != 0x80000000);
329    
330     if (arg < 0)     if (arg < 0)
331        return(-arg);        return(-arg);
332     else     else
333        return(arg);        return(arg);
334     }     }
335    
336    
337  //07/15/01:  Visual inspection only.  Function deemed too  //07/15/01:  Visual inspection only.  Function deemed too
338  //simple for unit testing.  //simple for unit testing.
339  int GMP_INTS_mpz_sgn(const GMP_INTS_mpz_struct *arg)  int GMP_INTS_mpz_sgn(const GMP_INTS_mpz_struct *arg)
340     {     {
341     assert(arg != NULL);     assert(arg != NULL);
342     assert(arg->n_allocd > 0);     assert(arg->n_allocd > 0);
343    
344     if (arg->size > 0)     if (arg->size > 0)
345        return(1);        return(1);
346     else if (arg->size == 0)     else if (arg->size == 0)
347        return(0);        return(0);
348     else     else
349        return(-1);        return(-1);
350     }     }
351    
352    
353  //07/15/01:  Visual inspection only.  Function deemed too  //07/15/01:  Visual inspection only.  Function deemed too
354  //simple for unit testing.  //simple for unit testing.
355  int GMP_INTS_mpz_is_neg(const GMP_INTS_mpz_struct *arg)  int GMP_INTS_mpz_is_neg(const GMP_INTS_mpz_struct *arg)
356     {     {
357     assert(arg != NULL);     assert(arg != NULL);
358     assert(arg->n_allocd > 0);     assert(arg->n_allocd > 0);
359    
360     if (GMP_INTS_mpz_sgn(arg) == -1)     if (GMP_INTS_mpz_sgn(arg) == -1)
361        return(1);        return(1);
362     else     else
363        return(0);        return(0);
364     }     }
365    
366    
367  //07/15/01:  Visual inspection only.  Function deemed too  //07/15/01:  Visual inspection only.  Function deemed too
368  //simple for unit testing.  //simple for unit testing.
369  int GMP_INTS_mpz_is_zero(const GMP_INTS_mpz_struct *arg)  int GMP_INTS_mpz_is_zero(const GMP_INTS_mpz_struct *arg)
370     {     {
371     assert(arg != NULL);     assert(arg != NULL);
372     assert(arg->n_allocd > 0);     assert(arg->n_allocd > 0);
373    
374     if (GMP_INTS_mpz_sgn(arg) == 0)     if (GMP_INTS_mpz_sgn(arg) == 0)
375        return(1);        return(1);
376     else     else
377        return(0);        return(0);
378     }     }
379    
380    
381  //07/15/01:  Visual inspection only.  Function deemed too  //07/15/01:  Visual inspection only.  Function deemed too
382  //simple for unit testing.  //simple for unit testing.
383  int GMP_INTS_mpz_is_pos(const GMP_INTS_mpz_struct *arg)  int GMP_INTS_mpz_is_pos(const GMP_INTS_mpz_struct *arg)
384     {     {
385     assert(arg != NULL);     assert(arg != NULL);
386     assert(arg->n_allocd > 0);     assert(arg->n_allocd > 0);
387    
388     if (GMP_INTS_mpz_sgn(arg) == 1)     if (GMP_INTS_mpz_sgn(arg) == 1)
389        return(1);        return(1);
390     else     else
391        return(0);        return(0);
392     }     }
393    
394    
395  //07/15/01:  Visual inspection only.  Function deemed too  //07/15/01:  Visual inspection only.  Function deemed too
396  //simple for unit testing.  //simple for unit testing.
397  int GMP_INTS_mpz_is_odd(const GMP_INTS_mpz_struct *arg)  int GMP_INTS_mpz_is_odd(const GMP_INTS_mpz_struct *arg)
398     {     {
399     assert(arg != NULL);     assert(arg != NULL);
400     assert(arg->n_allocd > 0);     assert(arg->n_allocd > 0);
401    
402     if (arg->size == 0)     if (arg->size == 0)
403        return 0;        return 0;
404     else if ((arg->limbs[0] & 0x1) != 0)     else if ((arg->limbs[0] & 0x1) != 0)
405        return 1;        return 1;
406     else     else
407        return 0;        return 0;
408     }     }
409    
410    
411  //07/15/01:  Visual inspection only.  Function deemed too  //07/15/01:  Visual inspection only.  Function deemed too
412  //simple for unit testing.  //simple for unit testing.
413  int GMP_INTS_mpz_is_even(const GMP_INTS_mpz_struct *arg)  int GMP_INTS_mpz_is_even(const GMP_INTS_mpz_struct *arg)
414     {     {
415     assert(arg != NULL);     assert(arg != NULL);
416     assert(arg->n_allocd > 0);     assert(arg->n_allocd > 0);
417    
418     if (GMP_INTS_mpz_is_odd(arg))     if (GMP_INTS_mpz_is_odd(arg))
419        return 0;        return 0;
420     else     else
421        return 1;        return 1;
422     }     }
423    
424    
425  void GMP_INTS_mpz_negate(GMP_INTS_mpz_struct *arg)  void GMP_INTS_mpz_negate(GMP_INTS_mpz_struct *arg)
426     {     {
427     //Eyeball the input parameters.     //Eyeball the input parameters.
428     assert(arg != NULL);     assert(arg != NULL);
429     assert(arg->n_allocd > 0);     assert(arg->n_allocd > 0);
430     assert(arg->limbs != NULL);     assert(arg->limbs != NULL);
431    
432     arg->size = -(arg->size);     arg->size = -(arg->size);
433     }     }
434    
435    
436  //07/15/01:  Visual inspection only.  Function deemed too  //07/15/01:  Visual inspection only.  Function deemed too
437  //simple for unit testing.  //simple for unit testing.
438  void GMP_INTS_mpn_normalize(GMP_INTS_limb_ptr  limb_array,  void GMP_INTS_mpn_normalize(GMP_INTS_limb_ptr  limb_array,
439                              GMP_INTS_size_t   *idx)                              GMP_INTS_size_t   *idx)
440     {     {
441     assert(limb_array != NULL);     assert(limb_array != NULL);
442     assert(idx != NULL);     assert(idx != NULL);
443     assert(idx >= 0);     assert(idx >= 0);
444    
445     while (*idx > 0)     while (*idx > 0)
446        {        {
447        if (limb_array[*idx - 1] != 0)        if (limb_array[*idx - 1] != 0)
448           break;           break;
449        (*idx)--;        (*idx)--;
450        }        }
451     }     }
452    
453    
454  //07/15/01:  Visual inspection only.  Function deemed too  //07/15/01:  Visual inspection only.  Function deemed too
455  //simple for unit testing.  //simple for unit testing.
456  void GMP_INTS_mpn_copy_limbs(GMP_INTS_limb_ptr    dest,  void GMP_INTS_mpn_copy_limbs(GMP_INTS_limb_ptr    dest,
457                               GMP_INTS_limb_srcptr src,                               GMP_INTS_limb_srcptr src,
458                               GMP_INTS_size_t      n)                               GMP_INTS_size_t      n)
459    {    {
460    GMP_INTS_size_t i;    GMP_INTS_size_t i;
461    
462    assert(dest != NULL);    assert(dest != NULL);
463    assert(src != NULL);    assert(src != NULL);
464    assert(n >= 0);    assert(n >= 0);
465    
466    for (i=0; i<n; i++)    for (i=0; i<n; i++)
467      dest[i] = src[i];      dest[i] = src[i];
468    }    }
469    
470    
471  /******************************************************************/  /******************************************************************/
472  /***  LOW-LEVEL ARITHMETIC FUNCTIONS   ****************************/  /***  LOW-LEVEL ARITHMETIC FUNCTIONS   ****************************/
473  /******************************************************************/  /******************************************************************/
474    
475  int GMP_INTS_two_op_flags_map(int flags1, int flags2)  int GMP_INTS_two_op_flags_map(int flags1, int flags2)
476     {     {
477     int cf;     int cf;
478    
479     if (!flags1 && !flags2)     if (!flags1 && !flags2)
480        {        {
481        return(0);        return(0);
482        }        }
483     else     else
484        {        {
485        cf = flags1 | flags2;        cf = flags1 | flags2;
486    
487        if (cf & (GMP_INTS_EF_INTOVF_POS | GMP_INTS_EF_INTOVF_TAINT_POS))        if (cf & (GMP_INTS_EF_INTOVF_POS | GMP_INTS_EF_INTOVF_TAINT_POS))
488           {           {
489           //In either case here, the result will be tainted           //In either case here, the result will be tainted
490           //positive.           //positive.
491           return(GMP_INTS_EF_INTOVF_TAINT_POS);           return(GMP_INTS_EF_INTOVF_TAINT_POS);
492           }           }
493        else if (cf & (GMP_INTS_EF_INTOVF_NEG | GMP_INTS_EF_INTOVF_TAINT_NEG))        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           //In either case here, the result will be tainted
496           //negative.           //negative.
497           return(GMP_INTS_EF_INTOVF_TAINT_NEG);           return(GMP_INTS_EF_INTOVF_TAINT_NEG);
498           }           }
499        else        else
500           {           {
501           //This case is where the flags ints are non-zero, but           //This case is where the flags ints are non-zero, but
502           //no known bits are set.  This is surely some kind of           //no known bits are set.  This is surely some kind of
503           //an internal software error.  In debug mode, want to           //an internal software error.  In debug mode, want to
504           //alert to error.  In actual operation, want to just           //alert to error.  In actual operation, want to just
505           //pretend an ordinary positive taint.           //pretend an ordinary positive taint.
506           assert(0);           assert(0);
507           return(GMP_INTS_EF_INTOVF_TAINT_POS);           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,  GMP_INTS_limb_t GMP_INTS_mpn_add_1 (GMP_INTS_limb_ptr    res_ptr,
514                                      GMP_INTS_limb_srcptr s1_ptr,                                      GMP_INTS_limb_srcptr s1_ptr,
515                                      GMP_INTS_size_t      s1_size,                                      GMP_INTS_size_t      s1_size,
516                                      GMP_INTS_limb_t      s2_limb)                                      GMP_INTS_limb_t      s2_limb)
517     {     {
518     GMP_INTS_limb_t x;     GMP_INTS_limb_t x;
519    
520     assert(res_ptr != NULL);     assert(res_ptr != NULL);
521     assert(s1_ptr != NULL);     assert(s1_ptr != NULL);
522     assert(s1_size > 0);     assert(s1_size > 0);
523    
524     x = *s1_ptr++;     x = *s1_ptr++;
525     s2_limb = x + s2_limb;     s2_limb = x + s2_limb;
526     *res_ptr++ = s2_limb;     *res_ptr++ = s2_limb;
527     //Since limbs are unsigned, the test below tests if there     //Since limbs are unsigned, the test below tests if there
528     //was a carry, i.e. a positive rollover.     //was a carry, i.e. a positive rollover.
529     if (s2_limb < x)     if (s2_limb < x)
530        {        {
531        while (--s1_size != 0)        while (--s1_size != 0)
532           {           {
533           x = *s1_ptr++ + 1;           x = *s1_ptr++ + 1;
534           *res_ptr++ = x;           *res_ptr++ = x;
535           if (x != 0)           if (x != 0)
536              goto fin;              goto fin;
537           }           }
538                
539        return 1;        return 1;
540        }        }
541        
542     fin:     fin:
543     if (res_ptr != s1_ptr)     if (res_ptr != s1_ptr)
544        {        {
545        GMP_INTS_size_t i;        GMP_INTS_size_t i;
546        for (i = 0; i < s1_size - 1; i++)        for (i = 0; i < s1_size - 1; i++)
547           {           {
548           res_ptr[i] = s1_ptr[i];           res_ptr[i] = s1_ptr[i];
549           }           }
550        }        }
551    
552     return 0;     return 0;
553     }     }
554    
555    
556  GMP_INTS_limb_t GMP_INTS_mpn_sub_1(GMP_INTS_limb_ptr    res_ptr,  GMP_INTS_limb_t GMP_INTS_mpn_sub_1(GMP_INTS_limb_ptr    res_ptr,
557                                     GMP_INTS_limb_srcptr s1_ptr,                                     GMP_INTS_limb_srcptr s1_ptr,
558                                     GMP_INTS_size_t      s1_size,                                     GMP_INTS_size_t      s1_size,
559                                     GMP_INTS_limb_t      s2_limb)                                     GMP_INTS_limb_t      s2_limb)
560     {     {
561     GMP_INTS_limb_t x;     GMP_INTS_limb_t x;
562    
563     assert(res_ptr != NULL);     assert(res_ptr != NULL);
564     assert(s1_ptr  != NULL);     assert(s1_ptr  != NULL);
565     assert(s1_size > 0);     assert(s1_size > 0);
566    
567     x = *s1_ptr++;     x = *s1_ptr++;
568     s2_limb = x - s2_limb;     s2_limb = x - s2_limb;
569     *res_ptr++ = s2_limb;     *res_ptr++ = s2_limb;
570     //The test below detects a borrow.     //The test below detects a borrow.
571     if (s2_limb > x)     if (s2_limb > x)
572        {        {
573        while (--s1_size != 0)        while (--s1_size != 0)
574           {           {
575           x = *s1_ptr++;           x = *s1_ptr++;
576           *res_ptr++ = x - 1;           *res_ptr++ = x - 1;
577           if (x != 0)           if (x != 0)
578              goto fin;              goto fin;
579           }           }
580                
581        return 1;        return 1;
582        }        }
583        
584     fin:     fin:
585     if (res_ptr != s1_ptr)     if (res_ptr != s1_ptr)
586        {        {
587        GMP_INTS_size_t i;        GMP_INTS_size_t i;
588        for (i = 0; i < s1_size - 1; i++)        for (i = 0; i < s1_size - 1; i++)
589           {           {
590              res_ptr[i] = s1_ptr[i];              res_ptr[i] = s1_ptr[i];
591           }           }
592        }        }
593        return 0;        return 0;
594     }     }
595    
596    
597  //07/15/01:  Am willing to skip unit-testing on this.  //07/15/01:  Am willing to skip unit-testing on this.
598  //Understand the logic (i.e. passes visual inspection),  //Understand the logic (i.e. passes visual inspection),
599  //and comes from GNU-MP.  Hope any defects here will be  //and comes from GNU-MP.  Hope any defects here will be
600  //caught in testing of GMP_INTS_mpz_mul() and other  //caught in testing of GMP_INTS_mpz_mul() and other
601  //higher-level functions.  //higher-level functions.
602  GMP_INTS_limb_t GMP_INTS_mpn_mul_1 (GMP_INTS_limb_ptr res_ptr,  GMP_INTS_limb_t GMP_INTS_mpn_mul_1 (GMP_INTS_limb_ptr res_ptr,
603                                                                    GMP_INTS_limb_srcptr s1_ptr,                                                                    GMP_INTS_limb_srcptr s1_ptr,
604                                                                    GMP_INTS_size_t s1_size,                                                                    GMP_INTS_size_t s1_size,
605                                                                    GMP_INTS_limb_t s2_limb)                                                                    GMP_INTS_limb_t s2_limb)
606     {     {
607     GMP_INTS_limb_t cy_limb;     GMP_INTS_limb_t cy_limb;
608     GMP_INTS_size_t j;     GMP_INTS_size_t j;
609     GMP_INTS_limb_t prod_high, prod_low;     GMP_INTS_limb_t prod_high, prod_low;
610     unsigned _int64 temp;     unsigned _int64 temp;
611        
612     assert(res_ptr != NULL);     assert(res_ptr != NULL);
613     assert(s1_ptr != NULL);     assert(s1_ptr != NULL);
614     assert(s1_size > 0);     assert(s1_size > 0);
615    
616     /* The loop counter and index J goes from -S1_SIZE to -1.  This way     /* The loop counter and index J goes from -S1_SIZE to -1.  This way
617     the loop becomes faster.  */     the loop becomes faster.  */
618     j = -s1_size;     j = -s1_size;
619        
620     /* Offset the base pointers to compensate for the negative indices.  */     /* Offset the base pointers to compensate for the negative indices.  */
621     s1_ptr -= j;     s1_ptr -= j;
622     res_ptr -= j;     res_ptr -= j;
623        
624     cy_limb = 0;     cy_limb = 0;
625     do     do
626        {        {
627        //The original code here was the following macro:        //The original code here was the following macro:
628        //umul_ppmm (prod_high, prod_low, s1_ptr[j], s2_limb);        //umul_ppmm (prod_high, prod_low, s1_ptr[j], s2_limb);
629        //Will use the 64-bit data type of MSVC++ to achieve        //Will use the 64-bit data type of MSVC++ to achieve
630        //the same effect.        //the same effect.
631        //        //
632        //NOTE AS OF 07/13/01:  I have looked at the assembly-        //NOTE AS OF 07/13/01:  I have looked at the assembly-
633        //language, and the lines below are a real sore spot.        //language, and the lines below are a real sore spot.
634        //The multiply is fairly direct (although there is a        //The multiply is fairly direct (although there is a
635        //function call), but the shift does not behave as        //function call), but the shift does not behave as
636        //expected--there is a function call and a loop to        //expected--there is a function call and a loop to
637        //go through the 32 iterations.  After logical testing,        //go through the 32 iterations.  After logical testing,
638        //may want to clean this out--this would surely        //may want to clean this out--this would surely
639        //result in a speed increase.  This is a sore spot.        //result in a speed increase.  This is a sore spot.
640        temp = ((unsigned _int64)s1_ptr[j]) * ((unsigned _int64)s2_limb);        temp = ((unsigned _int64)s1_ptr[j]) * ((unsigned _int64)s2_limb);
641        prod_low  = (GMP_INTS_limb_t)temp;        prod_low  = (GMP_INTS_limb_t)temp;
642        prod_high = (GMP_INTS_limb_t)(temp >> 32);        prod_high = (GMP_INTS_limb_t)(temp >> 32);
643    
644        prod_low += cy_limb;        prod_low += cy_limb;
645        cy_limb = (prod_low < cy_limb) + prod_high;        cy_limb = (prod_low < cy_limb) + prod_high;
646                
647        res_ptr[j] = prod_low;        res_ptr[j] = prod_low;
648        }        }
649     while (++j != 0);     while (++j != 0);
650                
651     return cy_limb;     return cy_limb;
652     }     }
653    
654    
655  //07/15/01:  Am willing to skip unit-testing on this.  //07/15/01:  Am willing to skip unit-testing on this.
656  //Understand the logic (i.e. passes visual inspection),  //Understand the logic (i.e. passes visual inspection),
657  //and comes from GNU-MP.  Hope any defects here will be  //and comes from GNU-MP.  Hope any defects here will be
658  //caught in testing of GMP_INTS_mpz_add() and other  //caught in testing of GMP_INTS_mpz_add() and other
659  //higher-level functions.  //higher-level functions.
660  GMP_INTS_limb_t GMP_INTS_mpn_add_n(GMP_INTS_limb_ptr    res_ptr,    GMP_INTS_limb_t GMP_INTS_mpn_add_n(GMP_INTS_limb_ptr    res_ptr,  
661                                     GMP_INTS_limb_srcptr s1_ptr,                                     GMP_INTS_limb_srcptr s1_ptr,
662                                     GMP_INTS_limb_srcptr s2_ptr,                                     GMP_INTS_limb_srcptr s2_ptr,
663                                     GMP_INTS_size_t      size)                                     GMP_INTS_size_t      size)
664     {     {
665     GMP_INTS_limb_t x, y, cy;     GMP_INTS_limb_t x, y, cy;
666     GMP_INTS_size_t j;     GMP_INTS_size_t j;
667    
668     assert(res_ptr != NULL);     assert(res_ptr != NULL);
669     assert(s1_ptr  != NULL);     assert(s1_ptr  != NULL);
670     assert(s2_ptr  != NULL);     assert(s2_ptr  != NULL);
671    
672     /* The loop counter and index J goes from -SIZE to -1.  This way     /* The loop counter and index J goes from -SIZE to -1.  This way
673        the loop becomes faster.  */        the loop becomes faster.  */
674     j = -size;     j = -size;
675    
676     /* Offset the base pointers to compensate for the negative indices.  */     /* Offset the base pointers to compensate for the negative indices.  */
677     s1_ptr -= j;     s1_ptr -= j;
678     s2_ptr -= j;     s2_ptr -= j;
679     res_ptr -= j;     res_ptr -= j;
680    
681     cy = 0;     cy = 0;
682     do     do
683        {        {
684        y = s2_ptr[j];        y = s2_ptr[j];
685        x = s1_ptr[j];        x = s1_ptr[j];
686        y += cy;   /* add previous carry to one addend */        y += cy;   /* add previous carry to one addend */
687        cy = (y < cy);  /* get out carry from that addition */        cy = (y < cy);  /* get out carry from that addition */
688        y = x + y;  /* add other addend */        y = x + y;  /* add other addend */
689        cy = (y < x) + cy; /* get out carry from that add, combine */        cy = (y < x) + cy; /* get out carry from that add, combine */
690        res_ptr[j] = y;        res_ptr[j] = y;
691        }        }
692     while (++j != 0);     while (++j != 0);
693    
694     return cy;     return cy;
695     }     }
696    
697    
698  //07/15/01:  Am willing to skip unit-testing on this.  //07/15/01:  Am willing to skip unit-testing on this.
699  //Understand the logic (i.e. passes visual inspection),  //Understand the logic (i.e. passes visual inspection),
700  //and comes from GNU-MP.  Hope any defects here will be  //and comes from GNU-MP.  Hope any defects here will be
701  //caught in testing of GMP_INTS_mpz_mul() and other  //caught in testing of GMP_INTS_mpz_mul() and other
702  //higher-level functions.  //higher-level functions.
703  GMP_INTS_limb_t GMP_INTS_mpn_addmul_1 (GMP_INTS_limb_ptr    res_ptr,  GMP_INTS_limb_t GMP_INTS_mpn_addmul_1 (GMP_INTS_limb_ptr    res_ptr,
704                                         GMP_INTS_limb_srcptr s1_ptr,                                         GMP_INTS_limb_srcptr s1_ptr,
705                                         GMP_INTS_size_t      s1_size,                                         GMP_INTS_size_t      s1_size,
706                                         GMP_INTS_limb_t      s2_limb)                                         GMP_INTS_limb_t      s2_limb)
707     {     {
708     GMP_INTS_limb_t cy_limb;     GMP_INTS_limb_t cy_limb;
709     GMP_INTS_size_t j;     GMP_INTS_size_t j;
710     GMP_INTS_limb_t prod_high, prod_low;     GMP_INTS_limb_t prod_high, prod_low;
711     GMP_INTS_limb_t x;     GMP_INTS_limb_t x;
712     unsigned _int64 temp;     unsigned _int64 temp;
713    
714     //Eyeball the inputs carefully.     //Eyeball the inputs carefully.
715     assert(res_ptr != NULL);     assert(res_ptr != NULL);
716     assert(s1_ptr != NULL);     assert(s1_ptr != NULL);
717     assert(s1_size > 0);     assert(s1_size > 0);
718    
719     /* The loop counter and index J goes from -SIZE to -1.  This way     /* The loop counter and index J goes from -SIZE to -1.  This way
720        the loop becomes faster.  */        the loop becomes faster.  */
721     j = -s1_size;     j = -s1_size;
722    
723     /* Offset the base pointers to compensate for the negative indices.  */     /* Offset the base pointers to compensate for the negative indices.  */
724     res_ptr -= j;     res_ptr -= j;
725     s1_ptr -= j;     s1_ptr -= j;
726    
727     cy_limb = 0;     cy_limb = 0;
728     do     do
729        {        {
730        //The original code here was the following macro:        //The original code here was the following macro:
731        //umul_ppmm (prod_high, prod_low, s1_ptr[j], s2_limb);        //umul_ppmm (prod_high, prod_low, s1_ptr[j], s2_limb);
732        //Will use the 64-bit data type of MSVC++ to achieve        //Will use the 64-bit data type of MSVC++ to achieve
733        //the same effect.        //the same effect.
734        //        //
735        //NOTE AS OF 07/14/01:  I have not looked at the assembly-        //NOTE AS OF 07/14/01:  I have not looked at the assembly-
736        //language, but the assembly-language generated by what        //language, but the assembly-language generated by what
737        //is below is suspected to have performance problems.        //is below is suspected to have performance problems.
738        //May want to come back to this.        //May want to come back to this.
739        temp = ((unsigned _int64)s1_ptr[j]) * ((unsigned _int64)s2_limb);        temp = ((unsigned _int64)s1_ptr[j]) * ((unsigned _int64)s2_limb);
740        prod_low  = (GMP_INTS_limb_t)temp;        prod_low  = (GMP_INTS_limb_t)temp;
741        prod_high = (GMP_INTS_limb_t)(temp >> 32);        prod_high = (GMP_INTS_limb_t)(temp >> 32);
742    
743        prod_low += cy_limb;        prod_low += cy_limb;
744        cy_limb = (prod_low < cy_limb) + prod_high;        cy_limb = (prod_low < cy_limb) + prod_high;
745    
746        x = res_ptr[j];        x = res_ptr[j];
747        prod_low = x + prod_low;        prod_low = x + prod_low;
748        cy_limb += (prod_low < x);        cy_limb += (prod_low < x);
749        res_ptr[j] = prod_low;        res_ptr[j] = prod_low;
750        }        }
751     while (++j != 0);     while (++j != 0);
752    
753     return cy_limb;     return cy_limb;
754     }     }
755    
756    
757  //07/15/01:  Am willing to skip unit-testing on this.  //07/15/01:  Am willing to skip unit-testing on this.
758  //Understand the logic (i.e. passes visual inspection),  //Understand the logic (i.e. passes visual inspection),
759  //and comes from GNU-MP.  //and comes from GNU-MP.
760  GMP_INTS_limb_t GMP_INTS_mpn_add (GMP_INTS_limb_ptr    res_ptr,  GMP_INTS_limb_t GMP_INTS_mpn_add (GMP_INTS_limb_ptr    res_ptr,
761                                    GMP_INTS_limb_srcptr s1_ptr,                                    GMP_INTS_limb_srcptr s1_ptr,
762                                    GMP_INTS_size_t      s1_size,                                    GMP_INTS_size_t      s1_size,
763                                    GMP_INTS_limb_srcptr s2_ptr,                                    GMP_INTS_limb_srcptr s2_ptr,
764                                    GMP_INTS_size_t      s2_size)                                    GMP_INTS_size_t      s2_size)
765     {     {
766     GMP_INTS_limb_t cy_limb = 0;     GMP_INTS_limb_t cy_limb = 0;
767    
768     assert(res_ptr != NULL);     assert(res_ptr != NULL);
769     assert(s1_ptr  != NULL);     assert(s1_ptr  != NULL);
770     assert(s2_ptr  != NULL);     assert(s2_ptr  != NULL);
771    
772     //Numbers apparently must be arranged with sizes so that     //Numbers apparently must be arranged with sizes so that
773     //LIMBS(s1) >= LIMBS(s2).     //LIMBS(s1) >= LIMBS(s2).
774     //Add the parts up to the most significant limb of S2.     //Add the parts up to the most significant limb of S2.
775     if (s2_size != 0)     if (s2_size != 0)
776        cy_limb = GMP_INTS_mpn_add_n (res_ptr,        cy_limb = GMP_INTS_mpn_add_n (res_ptr,
777                                      s1_ptr,                                      s1_ptr,
778                                      s2_ptr,                                      s2_ptr,
779                                      s2_size);                                      s2_size);
780    
781     //Process the carry result, and propagate the carries up through     //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     //the parts of S1 that don't exist in S2, i.e. propagate the
783     //carries upward in S1.     //carries upward in S1.
784     if (s1_size - s2_size != 0)     if (s1_size - s2_size != 0)
785        cy_limb = GMP_INTS_mpn_add_1 (res_ptr + s2_size,        cy_limb = GMP_INTS_mpn_add_1 (res_ptr + s2_size,
786                                      s1_ptr + s2_size,                                      s1_ptr + s2_size,
787                                      s1_size - s2_size,                                      s1_size - s2_size,
788                                      cy_limb);                                      cy_limb);
789     return cy_limb;     return cy_limb;
790     }     }
791    
792    
793  //07/15/01:  Am willing to skip unit-testing on this.  //07/15/01:  Am willing to skip unit-testing on this.
794  //Understand the logic (i.e. passes visual inspection),  //Understand the logic (i.e. passes visual inspection),
795  //and comes from GNU-MP.  //and comes from GNU-MP.
796  GMP_INTS_limb_t GMP_INTS_mpn_sub_n(GMP_INTS_limb_ptr    res_ptr,  GMP_INTS_limb_t GMP_INTS_mpn_sub_n(GMP_INTS_limb_ptr    res_ptr,
797                                     GMP_INTS_limb_srcptr s1_ptr,                                     GMP_INTS_limb_srcptr s1_ptr,
798                                     GMP_INTS_limb_srcptr s2_ptr,                                     GMP_INTS_limb_srcptr s2_ptr,
799                                     GMP_INTS_size_t      size)                                     GMP_INTS_size_t      size)
800     {     {
801     GMP_INTS_limb_t x, y, cy;     GMP_INTS_limb_t x, y, cy;
802     GMP_INTS_size_t j;     GMP_INTS_size_t j;
803    
804     assert(res_ptr != NULL);     assert(res_ptr != NULL);
805     assert(s1_ptr  != NULL);     assert(s1_ptr  != NULL);
806     assert(s2_ptr  != NULL);     assert(s2_ptr  != NULL);
807    
808     /* The loop counter and index J goes from -SIZE to -1.  This way     /* The loop counter and index J goes from -SIZE to -1.  This way
809        the loop becomes faster.  */        the loop becomes faster.  */
810     j = -size;     j = -size;
811    
812     /* Offset the base pointers to compensate for the negative indices.  */     /* Offset the base pointers to compensate for the negative indices.  */
813     s1_ptr -= j;     s1_ptr -= j;
814     s2_ptr -= j;     s2_ptr -= j;
815     res_ptr -= j;     res_ptr -= j;
816    
817     cy = 0;     cy = 0;
818     do     do
819        {        {
820        y = s2_ptr[j];        y = s2_ptr[j];
821        x = s1_ptr[j];        x = s1_ptr[j];
822        y += cy;   /* add previous carry to subtrahend */        y += cy;   /* add previous carry to subtrahend */
823        cy = (y < cy);  /* get out carry from that addition */        cy = (y < cy);  /* get out carry from that addition */
824        y = x - y;  /* main subtract */        y = x - y;  /* main subtract */
825        cy = (y > x) + cy; /* get out carry from the subtract, combine */        cy = (y > x) + cy; /* get out carry from the subtract, combine */
826        res_ptr[j] = y;        res_ptr[j] = y;
827        }        }
828     while (++j != 0);     while (++j != 0);
829    
830     return cy;     return cy;
831     }     }
832    
833    
834  //07/17/01:  Am willing to skip unit-testing on this.  //07/17/01:  Am willing to skip unit-testing on this.
835  //Understand the logic (i.e. passes visual inspection),  //Understand the logic (i.e. passes visual inspection),
836  //and comes from GNU-MP.  //and comes from GNU-MP.
837  GMP_INTS_limb_t GMP_INTS_mpn_sub (GMP_INTS_limb_ptr    res_ptr,  GMP_INTS_limb_t GMP_INTS_mpn_sub (GMP_INTS_limb_ptr    res_ptr,
838                                    GMP_INTS_limb_srcptr s1_ptr,                                    GMP_INTS_limb_srcptr s1_ptr,
839                                    GMP_INTS_size_t      s1_size,                                    GMP_INTS_size_t      s1_size,
840                                    GMP_INTS_limb_srcptr s2_ptr,                                    GMP_INTS_limb_srcptr s2_ptr,
841                                    GMP_INTS_size_t      s2_size)                                    GMP_INTS_size_t      s2_size)
842     {     {
843     GMP_INTS_limb_t cy_limb = 0;     GMP_INTS_limb_t cy_limb = 0;
844    
845     assert(res_ptr != NULL);     assert(res_ptr != NULL);
846     assert(s1_ptr  != NULL);     assert(s1_ptr  != NULL);
847     assert(s2_ptr  != NULL);     assert(s2_ptr  != NULL);
848    
849     if (s2_size != 0)     if (s2_size != 0)
850        cy_limb = GMP_INTS_mpn_sub_n(res_ptr,        cy_limb = GMP_INTS_mpn_sub_n(res_ptr,
851                                     s1_ptr,                                     s1_ptr,
852                                     s2_ptr,                                     s2_ptr,
853                                     s2_size);                                     s2_size);
854    
855     if (s1_size - s2_size != 0)     if (s1_size - s2_size != 0)
856        cy_limb = GMP_INTS_mpn_sub_1(res_ptr + s2_size,        cy_limb = GMP_INTS_mpn_sub_1(res_ptr + s2_size,
857                                     s1_ptr + s2_size,                                     s1_ptr + s2_size,
858                                     s1_size - s2_size,                                     s1_size - s2_size,
859                                     cy_limb);                                     cy_limb);
860     return cy_limb;     return cy_limb;
861     }     }
862    
863    
864  //07/17/01:  Am willing to skip unit-testing on this.  //07/17/01:  Am willing to skip unit-testing on this.
865  //Understand the logic (i.e. passes visual inspection),  //Understand the logic (i.e. passes visual inspection),
866  //and comes from GNU-MP.  //and comes from GNU-MP.
867  GMP_INTS_limb_t GMP_INTS_mpn_lshift(GMP_INTS_limb_ptr    wp,  GMP_INTS_limb_t GMP_INTS_mpn_lshift(GMP_INTS_limb_ptr    wp,
868                                      GMP_INTS_limb_srcptr up,                                      GMP_INTS_limb_srcptr up,
869                                      GMP_INTS_size_t      usize,                                      GMP_INTS_size_t      usize,
870                                      unsigned int         cnt)                                      unsigned int         cnt)
871     {     {
872     GMP_INTS_limb_t high_limb, low_limb;     GMP_INTS_limb_t high_limb, low_limb;
873     unsigned sh_1, sh_2;     unsigned sh_1, sh_2;
874     GMP_INTS_size_t i;     GMP_INTS_size_t i;
875     GMP_INTS_limb_t retval;     GMP_INTS_limb_t retval;
876    
877     assert(wp != NULL);     assert(wp != NULL);
878     assert(up != NULL);     assert(up != NULL);
879     assert(usize > 0);     assert(usize > 0);
880     assert(cnt > 0);     assert(cnt > 0);
881    
882     sh_1 = cnt;     sh_1 = cnt;
883    
884     wp += 1;     wp += 1;
885     sh_2 = GMP_INTS_BITS_PER_LIMB - sh_1;     sh_2 = GMP_INTS_BITS_PER_LIMB - sh_1;
886       //This automatically implies that can't call this function to shift more       //This automatically implies that can't call this function to shift more
887       //than 31 places.       //than 31 places.
888     i = usize - 1;     i = usize - 1;
889     low_limb = up[i];     low_limb = up[i];
890     retval = low_limb >> sh_2;  //Return value is the amount shifted     retval = low_limb >> sh_2;  //Return value is the amount shifted
891                                 //off the top.                                 //off the top.
892     high_limb = low_limb;     high_limb = low_limb;
893     while (--i >= 0)     while (--i >= 0)
894        {        {
895        low_limb = up[i];        low_limb = up[i];
896        wp[i] = (high_limb << sh_1) | (low_limb >> sh_2);        wp[i] = (high_limb << sh_1) | (low_limb >> sh_2);
897        high_limb = low_limb;        high_limb = low_limb;
898        }        }
899     wp[i] = high_limb << sh_1;     wp[i] = high_limb << sh_1;
900    
901     return retval;     return retval;
902     }     }
903    
904    
905  //07/17/01:  Am willing to skip unit-testing on this.  //07/17/01:  Am willing to skip unit-testing on this.
906  //Understand the logic more or less (i.e. passes visual inspection),  //Understand the logic more or less (i.e. passes visual inspection),
907  //and comes from GNU-MP.  //and comes from GNU-MP.
908  /* Shift U (pointed to by UP and USIZE limbs long) CNT bits to the right  /* 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.     and store the USIZE least significant limbs of the result at WP.
910     The bits shifted out to the right are returned.     The bits shifted out to the right are returned.
911    
912     Argument constraints:     Argument constraints:
913     1. 0 < CNT < BITS_PER_MP_LIMB     1. 0 < CNT < BITS_PER_MP_LIMB
914     2. If the result is to be written over the input, WP must be <= UP.     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,  GMP_INTS_limb_t GMP_INTS_mpn_rshift (GMP_INTS_limb_ptr    wp,
917                                       GMP_INTS_limb_srcptr up,                                       GMP_INTS_limb_srcptr up,
918                                       GMP_INTS_size_t      usize,                                       GMP_INTS_size_t      usize,
919                                       unsigned int         cnt)                                       unsigned int         cnt)
920     {     {
921     GMP_INTS_limb_t high_limb, low_limb;     GMP_INTS_limb_t high_limb, low_limb;
922     unsigned sh_1, sh_2;     unsigned sh_1, sh_2;
923     GMP_INTS_size_t i;     GMP_INTS_size_t i;
924     GMP_INTS_limb_t retval;     GMP_INTS_limb_t retval;
925    
926     assert(wp != NULL);     assert(wp != NULL);
927     assert(up != NULL);     assert(up != NULL);
928     assert(usize > 0);     assert(usize > 0);
929     assert(cnt > 0);     assert(cnt > 0);
930    
931     sh_1 = cnt;     sh_1 = cnt;
932    
933     wp -= 1;     wp -= 1;
934     sh_2 = GMP_INTS_BITS_PER_LIMB - sh_1;     sh_2 = GMP_INTS_BITS_PER_LIMB - sh_1;
935     high_limb = up[0];     high_limb = up[0];
936     retval = high_limb << sh_2;     retval = high_limb << sh_2;
937     low_limb = high_limb;     low_limb = high_limb;
938    
939     for (i = 1; i < usize; i++)     for (i = 1; i < usize; i++)
940        {        {
941        high_limb = up[i];        high_limb = up[i];
942        wp[i] = (low_limb >> sh_1) | (high_limb << sh_2);        wp[i] = (low_limb >> sh_1) | (high_limb << sh_2);
943        low_limb = high_limb;        low_limb = high_limb;
944        }        }
945     wp[i] = low_limb >> sh_1;     wp[i] = low_limb >> sh_1;
946    
947     return retval;     return retval;
948     }     }
949    
950    
951  //07/17/01:  Am willing to skip unit-testing on this.  //07/17/01:  Am willing to skip unit-testing on this.
952  //Understand the logic (i.e. passes visual inspection),  //Understand the logic (i.e. passes visual inspection),
953  //and comes from GNU-MP.  //and comes from GNU-MP.
954  int GMP_INTS_mpn_cmp (GMP_INTS_limb_srcptr op1_ptr,  int GMP_INTS_mpn_cmp (GMP_INTS_limb_srcptr op1_ptr,
955                        GMP_INTS_limb_srcptr op2_ptr,                        GMP_INTS_limb_srcptr op2_ptr,
956                        GMP_INTS_size_t      size)                        GMP_INTS_size_t      size)
957     {     {
958     GMP_INTS_size_t i;     GMP_INTS_size_t i;
959     GMP_INTS_limb_t op1_word, op2_word;     GMP_INTS_limb_t op1_word, op2_word;
960    
961     assert(op1_ptr != NULL);     assert(op1_ptr != NULL);
962     assert(op2_ptr != NULL);     assert(op2_ptr != NULL);
963    
964     for (i = size - 1; i >= 0; i--)     for (i = size - 1; i >= 0; i--)
965        {        {
966        op1_word = op1_ptr[i];        op1_word = op1_ptr[i];
967        op2_word = op2_ptr[i];        op2_word = op2_ptr[i];
968        if (op1_word != op2_word)        if (op1_word != op2_word)
969           goto diff;           goto diff;
970        }        }
971     return 0;     return 0;
972    
973     diff:     diff:
974     //This can *not* be simplified to     //This can *not* be simplified to
975     //   op2_word - op2_word     //   op2_word - op2_word
976     //since that expression might give signed overflow.     //since that expression might give signed overflow.
977     return (op1_word > op2_word) ? 1 : -1;     return (op1_word > op2_word) ? 1 : -1;
978     }     }
979    
980    
981  //07/15/01:  Am willing to skip unit-testing on this.  //07/15/01:  Am willing to skip unit-testing on this.
982  //Understand the logic (i.e. passes visual inspection),  //Understand the logic (i.e. passes visual inspection),
983  //and comes from GNU-MP.  Hope any defects here will be  //and comes from GNU-MP.  Hope any defects here will be
984  //caught in testing of GMP_INTS_mpz_mul() and other  //caught in testing of GMP_INTS_mpz_mul() and other
985  //higher-level functions.  //higher-level functions.
986  void GMP_INTS_mpn_mul_basecase (GMP_INTS_limb_ptr prodp,  void GMP_INTS_mpn_mul_basecase (GMP_INTS_limb_ptr prodp,
987                                  GMP_INTS_limb_srcptr up,                                  GMP_INTS_limb_srcptr up,
988                                  GMP_INTS_size_t usize,                                  GMP_INTS_size_t usize,
989                                  GMP_INTS_limb_srcptr vp,                                  GMP_INTS_limb_srcptr vp,
990                                  GMP_INTS_size_t vsize)                                  GMP_INTS_size_t vsize)
991     {     {
992     assert(prodp != NULL);     assert(prodp != NULL);
993     assert(up != NULL);     assert(up != NULL);
994     assert(usize > 0);     assert(usize > 0);
995     assert(vp != NULL);     assert(vp != NULL);
996     assert(vsize > 0);     assert(vsize > 0);
997    
998     /* We first multiply by the low order one or two limbs, as the result can     /* 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        be stored, not added, to PROD.  We also avoid a loop for zeroing this
1000        way.  */        way.  */
1001     prodp[usize] = GMP_INTS_mpn_mul_1 (prodp, up, usize, vp[0]);     prodp[usize] = GMP_INTS_mpn_mul_1 (prodp, up, usize, vp[0]);
1002     prodp++;     prodp++;
1003     vp++;     vp++;
1004     vsize--;     vsize--;
1005    
1006     /* For each iteration in the loop, multiply U with one limb from V, and     /* For each iteration in the loop, multiply U with one limb from V, and
1007        add the result to PROD.  */        add the result to PROD.  */
1008     while (vsize != 0)     while (vsize != 0)
1009        {        {
1010        prodp[usize] = GMP_INTS_mpn_addmul_1 (prodp, up, usize, vp[0]);        prodp[usize] = GMP_INTS_mpn_addmul_1 (prodp, up, usize, vp[0]);
1011        prodp++,        prodp++,
1012        vp++,        vp++,
1013        vsize--;        vsize--;
1014        }        }
1015     }     }
1016    
1017    
1018  //07/15/01:  No unit testing possible--this is a passthrough.  //07/15/01:  No unit testing possible--this is a passthrough.
1019  //In the original GNU MP code, there were several multiplication  //In the original GNU MP code, there were several multiplication
1020  //algorithms, and this function would select one based on the  //algorithms, and this function would select one based on the
1021  //size of the operands and other considerations.  The code has been  //size of the operands and other considerations.  The code has been
1022  //pared so that only simple multiplication is used, which is why  //pared so that only simple multiplication is used, which is why
1023  //this function contains only a single pass-thru function call.  //this function contains only a single pass-thru function call.
1024  void GMP_INTS_mpn_mul_n (GMP_INTS_limb_ptr    p,  void GMP_INTS_mpn_mul_n (GMP_INTS_limb_ptr    p,
1025                           GMP_INTS_limb_srcptr a,                           GMP_INTS_limb_srcptr a,
1026                           GMP_INTS_limb_srcptr b,                           GMP_INTS_limb_srcptr b,
1027                           GMP_INTS_size_t      n)                           GMP_INTS_size_t      n)
1028     {     {
1029     GMP_INTS_mpn_mul_basecase (p, a, n, b, n);     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.  //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,  GMP_INTS_limb_t GMP_INTS_mpn_mul(GMP_INTS_limb_ptr    prodp,
1035                                   GMP_INTS_limb_srcptr up,                                   GMP_INTS_limb_srcptr up,
1036                                   GMP_INTS_size_t      un,                                   GMP_INTS_size_t      un,
1037                                   GMP_INTS_limb_srcptr vp,                                   GMP_INTS_limb_srcptr vp,
1038                                   GMP_INTS_size_t      vn)                                   GMP_INTS_size_t      vn)
1039     {     {
1040     //This is a gutted version of the GNU MP function.  The GNU     //This is a gutted version of the GNU MP function.  The GNU
1041     //MP function considered the case of a square, and also     //MP function considered the case of a square, and also
1042     //better algorithms that pay off with large operands.     //better algorithms that pay off with large operands.
1043     //This gutted version uses only basic multiplication     //This gutted version uses only basic multiplication
1044     //(O(N**2)).     //(O(N**2)).
1045    
1046     //Eyeball the input parameters.     //Eyeball the input parameters.
1047     assert(prodp != NULL);     assert(prodp != NULL);
1048     assert(up != NULL);     assert(up != NULL);
1049     assert(un >= 0);     assert(un >= 0);
1050     assert(vp != NULL);     assert(vp != NULL);
1051     assert(vn >= 0);     assert(vn >= 0);
1052    
1053     /* Basic long multiplication. */     /* Basic long multiplication. */
1054     GMP_INTS_mpn_mul_basecase (prodp, up, un, vp, vn);     GMP_INTS_mpn_mul_basecase (prodp, up, un, vp, vn);
1055    
1056     //Return the most significant limb (which might be zero).       //Return the most significant limb (which might be zero).  
1057     //This is different than     //This is different than
1058     //most other functions, which return the spillover.     //most other functions, which return the spillover.
1059     return prodp[un + vn - 1];     return prodp[un + vn - 1];
1060     }     }
1061    
1062    
1063  /******************************************************************/  /******************************************************************/
1064  /***  LIMB SPACE REALLOCATION FUNCTIONS   *************************/  /***  LIMB SPACE REALLOCATION FUNCTIONS   *************************/
1065  /******************************************************************/  /******************************************************************/
1066    
1067  void *GMP_INTS_mpz_realloc (GMP_INTS_mpz_struct *m,  void *GMP_INTS_mpz_realloc (GMP_INTS_mpz_struct *m,
1068                              GMP_INTS_size_t      new_size)                              GMP_INTS_size_t      new_size)
1069     {     {
1070     /* Never allocate zero space. */     /* Never allocate zero space. */
1071     if (new_size <= 0)     if (new_size <= 0)
1072        new_size = 1;        new_size = 1;
1073    
1074     m->limbs = (GMP_INTS_limb_ptr)     m->limbs = (GMP_INTS_limb_ptr)
1075                 GMP_INTS_realloc_w_size (m->limbs,                 GMP_INTS_realloc_w_size (m->limbs,
1076                                          m->n_allocd * sizeof(GMP_INTS_limb_t),                                          m->n_allocd * sizeof(GMP_INTS_limb_t),
1077                                          new_size * sizeof(GMP_INTS_limb_t));                                          new_size * sizeof(GMP_INTS_limb_t));
1078     m->n_allocd = new_size;     m->n_allocd = new_size;
1079    
1080     return (void *) m->limbs;     return (void *) m->limbs;
1081     }     }
1082    
1083    
1084  /******************************************************************/  /******************************************************************/
1085  /***  PUBLIC INITIALIZATION AND MEMORY MANAGEMENT FUNCTIONS   *****/  /***  PUBLIC INITIALIZATION AND MEMORY MANAGEMENT FUNCTIONS   *****/
1086  /******************************************************************/  /******************************************************************/
1087    
1088  void GMP_INTS_mpz_init (GMP_INTS_mpz_struct *x)  void GMP_INTS_mpz_init (GMP_INTS_mpz_struct *x)
1089     {     {
1090     assert(x != NULL);     assert(x != NULL);
1091    
1092     //The structure (the header block) exists in the     //The structure (the header block) exists in the
1093     //caller's area.  Most likely it is a local variable.     //caller's area.  Most likely it is a local variable.
1094     //This is OK, because it doesn't take up much space.     //This is OK, because it doesn't take up much space.
1095    
1096     //Start off with no errors.     //Start off with no errors.
1097     x->flags = 0;     x->flags = 0;
1098        
1099     //Allocate space for one limb, which is the most     //Allocate space for one limb, which is the most
1100     //basic amount.  This will grow, almost certainly.     //basic amount.  This will grow, almost certainly.
1101     x->limbs = GMP_INTS_malloc(sizeof(GMP_INTS_limb_t));     x->limbs = GMP_INTS_malloc(sizeof(GMP_INTS_limb_t));
1102    
1103     //Indicate that one limb was allocated.     //Indicate that one limb was allocated.
1104     x->n_allocd = 1;     x->n_allocd = 1;
1105    
1106     //Set the size to 0.  This signals a value of zero.     //Set the size to 0.  This signals a value of zero.
1107     x->size = 0;     x->size = 0;
1108     }     }
1109    
1110    
1111  void GMP_INTS_mpz_clear (GMP_INTS_mpz_struct *x)  void GMP_INTS_mpz_clear (GMP_INTS_mpz_struct *x)
1112     {     {
1113     //Be sure the passed pointer is not NULL.     //Be sure the passed pointer is not NULL.
1114     assert(x != NULL);     assert(x != NULL);
1115    
1116     //Be sure that the amount allocated is also above zero.     //Be sure that the amount allocated is also above zero.
1117     //Anything else represents a logical error.     //Anything else represents a logical error.
1118     assert(x->n_allocd > 0);     assert(x->n_allocd > 0);
1119    
1120     //Be sure that the pointer to the allocated limbs     //Be sure that the pointer to the allocated limbs
1121     //is not NULL.  Anything else would be a logical     //is not NULL.  Anything else would be a logical
1122     //error.     //error.
1123     assert(x->limbs != NULL);     assert(x->limbs != NULL);
1124    
1125     //Deallocate the space for the limbs.  The pointer is     //Deallocate the space for the limbs.  The pointer is
1126     //set NULL and the allocated amount set to zero     //set NULL and the allocated amount set to zero
1127     // so in case clear is called again it will be     // so in case clear is called again it will be
1128     //a detectable error.     //a detectable error.
1129     GMP_INTS_free_w_size(x->limbs,     GMP_INTS_free_w_size(x->limbs,
1130                          x->n_allocd * sizeof(GMP_INTS_limb_t));                          x->n_allocd * sizeof(GMP_INTS_limb_t));
1131     x->limbs = NULL;     x->limbs = NULL;
1132     x->n_allocd = 0;     x->n_allocd = 0;
1133     }     }
1134    
1135    
1136  /******************************************************************/  /******************************************************************/
1137  /***  PUBLIC ASSIGNMENT FUNCTIONS   *******************************/  /***  PUBLIC ASSIGNMENT FUNCTIONS   *******************************/
1138  /******************************************************************/  /******************************************************************/
1139    
1140  void GMP_INTS_mpz_copy(      GMP_INTS_mpz_struct *dst,  void GMP_INTS_mpz_copy(      GMP_INTS_mpz_struct *dst,
1141                         const GMP_INTS_mpz_struct *src)                         const GMP_INTS_mpz_struct *src)
1142     {     {
1143     GMP_INTS_size_t i, n;     GMP_INTS_size_t i, n;
1144    
1145     //Eyeball the input parameters.     //Eyeball the input parameters.
1146     assert(dst != NULL);     assert(dst != NULL);
1147     assert(dst->n_allocd > 0);     assert(dst->n_allocd > 0);
1148     assert(dst->limbs != NULL);     assert(dst->limbs != NULL);
1149     assert(src != NULL);     assert(src != NULL);
1150     assert(src->n_allocd > 0);     assert(src->n_allocd > 0);
1151     assert(src->limbs != NULL);     assert(src->limbs != NULL);
1152    
1153     //Source and destination may not be the same.     //Source and destination may not be the same.
1154     assert(src != dst);     assert(src != dst);
1155        
1156     //Figure out the real size of the source.  We need to take the absolute     //Figure out the real size of the source.  We need to take the absolute
1157     //value.     //value.
1158     n = GMP_INTS_abs_of_size_t(src->size);     n = GMP_INTS_abs_of_size_t(src->size);
1159    
1160     //Reallocate the destination to be bigger if necessary.     //Reallocate the destination to be bigger if necessary.
1161     if (dst->n_allocd < n)     if (dst->n_allocd < n)
1162        {        {
1163        GMP_INTS_mpz_realloc (dst, n);        GMP_INTS_mpz_realloc (dst, n);
1164        }        }
1165    
1166     //Copy the non-dynamic fields in the header.     //Copy the non-dynamic fields in the header.
1167     dst->flags     = src->flags;     dst->flags     = src->flags;
1168     dst->size      = src->size;     dst->size      = src->size;
1169    
1170     //Copy the limbs.     //Copy the limbs.
1171     for (i=0; i<n; i++)     for (i=0; i<n; i++)
1172        dst->limbs[i] = src->limbs[i];        dst->limbs[i] = src->limbs[i];
1173     }     }
1174    
1175    
1176  void GMP_INTS_mpz_set_ui (GMP_INTS_mpz_struct *dest,  void GMP_INTS_mpz_set_ui (GMP_INTS_mpz_struct *dest,
1177                            unsigned long int val)                            unsigned long int val)
1178     {     {
1179     assert(dest != NULL);     assert(dest != NULL);
1180        
1181     /* We don't check if the allocation is enough, since the rest of the     /* 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.  */        package ensures it's at least 1, which is what we need here.  */
1183    
1184     dest->flags = 0;     dest->flags = 0;
1185        //A set operation resets any errors.        //A set operation resets any errors.
1186    
1187     if (val > 0)     if (val > 0)
1188        {        {
1189        dest->limbs[0] = val;        dest->limbs[0] = val;
1190        dest->size = 1;        dest->size = 1;
1191        }        }
1192     else     else
1193        {        {
1194        dest->size = 0;        dest->size = 0;
1195        }        }
1196     }     }
1197    
1198    
1199  void GMP_INTS_mpz_set_si (GMP_INTS_mpz_struct *dest,  void GMP_INTS_mpz_set_si (GMP_INTS_mpz_struct *dest,
1200                            signed long int val)                            signed long int val)
1201     {     {
1202     assert(dest != NULL);     assert(dest != NULL);
1203    
1204     /* We don't check if the allocation is enough, since the rest of the     /* 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.  */        package ensures it's at least 1, which is what we need here.  */
1206    
1207     dest->flags = 0;     dest->flags = 0;
1208       //A set operation resets any errors.       //A set operation resets any errors.
1209    
1210     if (val > 0)     if (val > 0)
1211        {        {
1212        dest->limbs[0] = val;        dest->limbs[0] = val;
1213        dest->size = 1;        dest->size = 1;
1214        }        }
1215     else if (val < 0)     else if (val < 0)
1216        {        {
1217        dest->limbs[0] = (unsigned long) -val;        dest->limbs[0] = (unsigned long) -val;
1218        dest->size = -1;        dest->size = -1;
1219        }        }
1220     else     else
1221        {        {
1222        dest->size = 0;        dest->size = 0;
1223        }        }
1224     }     }
1225    
1226    
1227  void GMP_INTS_mpz_set_simple_char_str(GMP_INTS_mpz_struct *z,  void GMP_INTS_mpz_set_simple_char_str(GMP_INTS_mpz_struct *z,
1228                                        const char *s)                                        const char *s)
1229     {     {
1230     int sign=1;     int sign=1;
1231     int digval;     int digval;
1232     GMP_INTS_mpz_struct digvalz, k10;     GMP_INTS_mpz_struct digvalz, k10;
1233    
1234     //Eyeball the arguments.     //Eyeball the arguments.
1235     assert(z != NULL);     assert(z != NULL);
1236     assert(z->n_allocd > 0);     assert(z->n_allocd > 0);
1237     assert(z->limbs != NULL);     assert(z->limbs != NULL);
1238     assert(s != NULL);     assert(s != NULL);
1239    
1240     //Set the arbitrary integer to zero.  This will also kill     //Set the arbitrary integer to zero.  This will also kill
1241     //any error flags.     //any error flags.
1242     GMP_INTS_mpz_set_ui(z, 0);     GMP_INTS_mpz_set_ui(z, 0);
1243    
1244     //Allocate an integer for our private use to hold each digit     //Allocate an integer for our private use to hold each digit
1245     //value.     //value.
1246     GMP_INTS_mpz_init(&digvalz);     GMP_INTS_mpz_init(&digvalz);
1247    
1248     //Allocate the constant 10, which we will use often.     //Allocate the constant 10, which we will use often.
1249     GMP_INTS_mpz_init(&k10);     GMP_INTS_mpz_init(&k10);
1250     GMP_INTS_mpz_set_ui(&k10, 10);     GMP_INTS_mpz_set_ui(&k10, 10);
1251    
1252     //As long as there are are digits and no flags set, keep     //As long as there are are digits and no flags set, keep
1253     //multiplying and adding the value of the digit.  Non-     //multiplying and adding the value of the digit.  Non-
1254     //digits are simply ignored.     //digits are simply ignored.
1255     while (!(z->flags) && (*s))     while (!(z->flags) && (*s))
1256        {        {
1257        if (*s == '-')        if (*s == '-')
1258           {           {
1259           sign = -sign;           sign = -sign;
1260           }           }
1261        else        else
1262           {           {
1263           digval = CHARFUNC_digit_to_val(*s);           digval = CHARFUNC_digit_to_val(*s);
1264           if (digval >= 0)           if (digval >= 0)
1265              {              {
1266              GMP_INTS_mpz_set_ui(&digvalz, digval);              GMP_INTS_mpz_set_ui(&digvalz, digval);
1267              GMP_INTS_mpz_mul(z, z, &k10);              GMP_INTS_mpz_mul(z, z, &k10);
1268              GMP_INTS_mpz_add(z, z, &digvalz);              GMP_INTS_mpz_add(z, z, &digvalz);
1269              }              }
1270           }           }
1271        s++;        s++;
1272        }        }
1273    
1274     //Adjust the final sign of the result.     //Adjust the final sign of the result.
1275     if (sign < 0)     if (sign < 0)
1276        z->size = -(z->size);        z->size = -(z->size);
1277    
1278     //Deallocate our temporary integers.     //Deallocate our temporary integers.
1279     GMP_INTS_mpz_clear(&digvalz);     GMP_INTS_mpz_clear(&digvalz);
1280     GMP_INTS_mpz_clear(&k10);     GMP_INTS_mpz_clear(&k10);
1281     }     }
1282    
1283    
1284  void GMP_INTS_mpz_set_sci_not_num(GMP_INTS_mpz_struct *z,  void GMP_INTS_mpz_set_sci_not_num(GMP_INTS_mpz_struct *z,
1285                                    int *failure,                                    int *failure,
1286                                    const char *s)                                    const char *s)
1287     {     {
1288     int parse_failure;     int parse_failure;
1289        //Return code from the floating point parsing        //Return code from the floating point parsing
1290        //function.        //function.
1291     char mant_sign;     char mant_sign;
1292        //Sign character, if any, from the mantissa,        //Sign character, if any, from the mantissa,
1293        //or N otherwise.        //or N otherwise.
1294     size_t mant_bdp;     size_t mant_bdp;
1295        //The index to the start of the mantissa before        //The index to the start of the mantissa before
1296        //the decimal point.        //the decimal point.
1297     size_t mant_bdp_len;     size_t mant_bdp_len;
1298        //The length of the mantissa before the decimal        //The length of the mantissa before the decimal
1299        //point.  Zero means not defined, i.e. that        //point.  Zero means not defined, i.e. that
1300        //no characters were parsed and interpreted as        //no characters were parsed and interpreted as
1301        //that part of a floating point number.        //that part of a floating point number.
1302     size_t mant_adp;     size_t mant_adp;
1303     size_t mant_adp_len;     size_t mant_adp_len;
1304        //Similar fields for after the decimal point.        //Similar fields for after the decimal point.
1305     char exp_sign;     char exp_sign;
1306        //Sign of the exponent, if any, or N otherwise.        //Sign of the exponent, if any, or N otherwise.
1307     size_t exp;     size_t exp;
1308     size_t exp_len;     size_t exp_len;
1309        //Similar fields as to the mantissa, but for the        //Similar fields as to the mantissa, but for the
1310        //exponent.        //exponent.
1311     size_t si;     size_t si;
1312        //Iteration variable.        //Iteration variable.
1313     int exponent_val;     int exponent_val;
1314        //The value of the exponent.  We can't accept        //The value of the exponent.  We can't accept
1315        //an exponent outside the range of a 24-bit        //an exponent outside the range of a 24-bit
1316        //signed integer.  The 24-bit limit is arbitrary.        //signed integer.  The 24-bit limit is arbitrary.
1317        //For one thing, it gives room to detect overflow        //For one thing, it gives room to detect overflow
1318        //as are adding and multiplying by 10.        //as are adding and multiplying by 10.
1319    
1320     //Eyeball the input parameters.     //Eyeball the input parameters.
1321     assert(z != NULL);     assert(z != NULL);
1322     assert(z->n_allocd > 0);     assert(z->n_allocd > 0);
1323     assert(z->limbs != NULL);     assert(z->limbs != NULL);
1324     assert(failure != NULL);     assert(failure != NULL);
1325     assert(s != NULL);     assert(s != NULL);
1326    
1327     //Start off believing no failure.     //Start off believing no failure.
1328     *failure = 0;     *failure = 0;
1329    
1330     //Set the output to zero.  This is the default case for some     //Set the output to zero.  This is the default case for some
1331     //steps below.     //steps below.
1332     GMP_INTS_mpz_set_ui(z, 0);     GMP_INTS_mpz_set_ui(z, 0);
1333    
1334     //Attempt to parse the number as a general number     //Attempt to parse the number as a general number
1335     //in scientific notation.     //in scientific notation.
1336     BSTRFUNC_parse_gen_sci_not_num(s,     BSTRFUNC_parse_gen_sci_not_num(s,
1337                                    &parse_failure,                                    &parse_failure,
1338                                    &mant_sign,                                    &mant_sign,
1339                                    &mant_bdp,                                    &mant_bdp,
1340                                    &mant_bdp_len,                                    &mant_bdp_len,
1341                                    &mant_adp,                                    &mant_adp,
1342                                    &mant_adp_len,                                    &mant_adp_len,
1343                                    &exp_sign,                                    &exp_sign,
1344                                    &exp,                                    &exp,
1345                                    &exp_len);                                    &exp_len);
1346    
1347     //If it wouldn't parse as a general number, can't go further.     //If it wouldn't parse as a general number, can't go further.
1348     if (parse_failure)     if (parse_failure)
1349        {        {
1350        *failure = 1;        *failure = 1;
1351        return;        return;
1352        }        }
1353     else if (!exp_len && !mant_adp_len)     else if (!exp_len && !mant_adp_len)
1354        {        {
1355        //There was no exponent, and no portion after        //There was no exponent, and no portion after
1356        //the decimal point.  Can just parse as an integer.        //the decimal point.  Can just parse as an integer.
1357        char *temp_buf;        char *temp_buf;
1358    
1359        //Allocate the temporary buffer to be one character longer        //Allocate the temporary buffer to be one character longer
1360        //than the length specified for the parsed mantissa.        //than the length specified for the parsed mantissa.
1361        temp_buf = GMP_INTS_malloc(sizeof(char) * (mant_bdp_len + 1));        temp_buf = GMP_INTS_malloc(sizeof(char) * (mant_bdp_len + 1));
1362    
1363        //Copy from the parsed area into the temporary buffer.        //Copy from the parsed area into the temporary buffer.
1364        for (si=mant_bdp; si<(mant_bdp+mant_bdp_len); si++)        for (si=mant_bdp; si<(mant_bdp+mant_bdp_len); si++)
1365           temp_buf[si-mant_bdp] = s[si];           temp_buf[si-mant_bdp] = s[si];
1366        temp_buf[mant_bdp_len] = 0;        temp_buf[mant_bdp_len] = 0;
1367    
1368        //Set the arbitrary integer to the value of the character        //Set the arbitrary integer to the value of the character
1369        //string.        //string.
1370        GMP_INTS_mpz_set_simple_char_str(z, temp_buf);        GMP_INTS_mpz_set_simple_char_str(z, temp_buf);
1371    
1372        //If the number parsed as negative, invert.        //If the number parsed as negative, invert.
1373        if (mant_sign == '-')        if (mant_sign == '-')
1374           z->size = -z->size;           z->size = -z->size;
1375    
1376        //Deallocate the temporary buffer.        //Deallocate the temporary buffer.
1377        GMP_INTS_free(temp_buf);        GMP_INTS_free(temp_buf);
1378        }        }
1379     else if (!exp_len && mant_adp_len)     else if (!exp_len && mant_adp_len)
1380        {        {
1381        char *temp_buf;        char *temp_buf;
1382    
1383        //In this case, there are digits after the decimal point,        //In this case, there are digits after the decimal point,
1384        //but no exponent specified.  The only way this makes        //but no exponent specified.  The only way this makes
1385        //sense is if all of the digits are zero--otherwise it        //sense is if all of the digits are zero--otherwise it
1386        //cannot be an integer.        //cannot be an integer.
1387        for (si=mant_adp; si<(mant_adp+mant_adp_len); si++)        for (si=mant_adp; si<(mant_adp+mant_adp_len); si++)
1388           {           {
1389           if (s[si] != '0')           if (s[si] != '0')
1390              {              {
1391              *failure = 1;              *failure = 1;
1392              return;              return;
1393              }              }
1394           }           }
1395    
1396        //We're clean.  They are only zeros.  Execute as per        //We're clean.  They are only zeros.  Execute as per
1397        //integer code.        //integer code.
1398    
1399        //Allocate the temporary buffer to be one character longer        //Allocate the temporary buffer to be one character longer
1400        //than the length specified for the parsed mantissa.        //than the length specified for the parsed mantissa.
1401        temp_buf = GMP_INTS_malloc(sizeof(char) * (mant_bdp_len + 1));        temp_buf = GMP_INTS_malloc(sizeof(char) * (mant_bdp_len + 1));
1402    
1403        //Copy from the parsed area into the temporary buffer.        //Copy from the parsed area into the temporary buffer.
1404        for (si=mant_bdp; si<(mant_bdp+mant_bdp_len); si++)        for (si=mant_bdp; si<(mant_bdp+mant_bdp_len); si++)
1405           temp_buf[si-mant_bdp] = s[si];           temp_buf[si-mant_bdp] = s[si];
1406        temp_buf[mant_bdp_len] = 0;        temp_buf[mant_bdp_len] = 0;
1407    
1408        //Set the arbitrary integer to the value of the character        //Set the arbitrary integer to the value of the character
1409        //string.        //string.
1410        GMP_INTS_mpz_set_simple_char_str(z, temp_buf);        GMP_INTS_mpz_set_simple_char_str(z, temp_buf);
1411    
1412        //If the number parsed as negative, invert.        //If the number parsed as negative, invert.
1413        if (mant_sign == '-')        if (mant_sign == '-')
1414           z->size = -z->size;           z->size = -z->size;
1415    
1416        //Deallocate the temporary buffer.        //Deallocate the temporary buffer.
1417        GMP_INTS_free(temp_buf);        GMP_INTS_free(temp_buf);
1418        }        }
1419     else if (exp_len)     else if (exp_len)
1420        {        {
1421        //This is the most difficult case, where an exponent        //This is the most difficult case, where an exponent
1422        //is specified.  There are several complex subcases,        //is specified.  There are several complex subcases,
1423        //such as:        //such as:
1424        //  a)If the exponent is too positive or too negative,        //  a)If the exponent is too positive or too negative,
1425        //    we can't use it.  In general, we won't tackle        //    we can't use it.  In general, we won't tackle
1426        //    an exponent that won't fit in a signed 24-bit        //    an exponent that won't fit in a signed 24-bit
1427        //    integer.  This provides a range of from        //    integer.  This provides a range of from
1428        //    -8,388,608 to +8,388,607.  This dwarfs the        //    -8,388,608 to +8,388,607.  This dwarfs the
1429        //    100,000 or so digit preprocessor limit,        //    100,000 or so digit preprocessor limit,
1430        //    and should be adequate for any practical        //    and should be adequate for any practical
1431        //    application.        //    application.
1432        //  b)If the exponent is zero, we ignore it.        //  b)If the exponent is zero, we ignore it.
1433        //  c)If the exponent is positive, it has to        //  c)If the exponent is positive, it has to
1434        //    be large enough to overcome any        //    be large enough to overcome any
1435        //    digits past the decimal point, otherwise        //    digits past the decimal point, otherwise
1436        //    we don't end up with an integer.        //    we don't end up with an integer.
1437        //  d)If the exponent is negative, there have to        //  d)If the exponent is negative, there have to
1438        //    be enough digits so that an integer remains        //    be enough digits so that an integer remains
1439        //    after the exponent is applied.  This        //    after the exponent is applied.  This
1440        //    generally requires trailing zeros on the        //    generally requires trailing zeros on the
1441        //    part before the decimal point.        //    part before the decimal point.
1442    
1443        //First, tackle the exponent.  Process the        //First, tackle the exponent.  Process the
1444        //exponent into a signed integer.  We have to        //exponent into a signed integer.  We have to
1445        //balk at anything outside of 24 bits.  The        //balk at anything outside of 24 bits.  The
1446        //procedure used automatically handles        //procedure used automatically handles
1447        //leading zeros correctly.        //leading zeros correctly.
1448        exponent_val = 0;        exponent_val = 0;
1449        for (si=exp; si<(exp+exp_len); si++)        for (si=exp; si<(exp+exp_len); si++)
1450           {           {
1451           int val;           int val;
1452    
1453           val = CHARFUNC_digit_to_val(s[si]);           val = CHARFUNC_digit_to_val(s[si]);
1454    
1455           assert(val >= 0 && val <= 9);           assert(val >= 0 && val <= 9);
1456    
1457           exponent_val *= 10;           exponent_val *= 10;
1458           exponent_val += val;           exponent_val += val;
1459    
1460           if (((exp_sign=='-') && (exponent_val>8388608))           if (((exp_sign=='-') && (exponent_val>8388608))
1461              ||              ||
1462              ((exp_sign != '-') && (exponent_val>8388607)))              ((exp_sign != '-') && (exponent_val>8388607)))
1463              {              {
1464              *failure = 1;              *failure = 1;
1465              return;              return;
1466              }              }
1467           }           }
1468    
1469        //If we're here, the exponent has been computed and        //If we're here, the exponent has been computed and
1470        //is within 24 bits.  However, we need to adjust for        //is within 24 bits.  However, we need to adjust for
1471        //the sign.        //the sign.
1472        if (exp_sign == '-')        if (exp_sign == '-')
1473           exponent_val = -exponent_val;           exponent_val = -exponent_val;
1474    
1475        //We need to make accurate assertions about the        //We need to make accurate assertions about the
1476        //portion of the number, if any, after the decimal point.        //portion of the number, if any, after the decimal point.
1477        //This means that we need to effectively discard        //This means that we need to effectively discard
1478        //trailing zeros.  To do this, we do not need to        //trailing zeros.  To do this, we do not need to
1479        //relocate the string, we can just back off the index        //relocate the string, we can just back off the index
1480        //to bypass any trailing zeros.        //to bypass any trailing zeros.
1481        while ((mant_adp_len > 0) && (s[mant_adp + mant_adp_len - 1]=='0'))        while ((mant_adp_len > 0) && (s[mant_adp + mant_adp_len - 1]=='0'))
1482           mant_adp_len--;           mant_adp_len--;
1483    
1484        //We also need to make accurate assertions about the        //We also need to make accurate assertions about the
1485        //portion of the number, if any, before the decimal        //portion of the number, if any, before the decimal
1486        //point.  It is known that the parsing function        //point.  It is known that the parsing function
1487        //isn't tolerant of multiple zeros, but zero is a        //isn't tolerant of multiple zeros, but zero is a
1488        //special case.  Let's advance the pointer to the        //special case.  Let's advance the pointer to the
1489        //part before the decimal point so that zero will        //part before the decimal point so that zero will
1490        //have zero length.        //have zero length.
1491        while ((mant_bdp_len > 0) && (s[mant_bdp]=='0'))        while ((mant_bdp_len > 0) && (s[mant_bdp]=='0'))
1492           {           {
1493           mant_bdp++;           mant_bdp++;
1494           mant_bdp_len--;           mant_bdp_len--;
1495           }           }
1496    
1497        //If we are dealing with zero, who cares about the        //If we are dealing with zero, who cares about the
1498        //exponent?  Just return the value of zero.        //exponent?  Just return the value of zero.
1499        if (!mant_bdp_len && !mant_adp_len)        if (!mant_bdp_len && !mant_adp_len)
1500           {           {
1501           *failure = 0;           *failure = 0;
1502           GMP_INTS_mpz_set_ui(z, 0);           GMP_INTS_mpz_set_ui(z, 0);
1503           return;           return;
1504           }           }
1505                
1506        //Beyond this point, we have something non-zero.        //Beyond this point, we have something non-zero.
1507        //If the exponent is positive, it must be at least        //If the exponent is positive, it must be at least
1508        //as large as the number of digits beyond the        //as large as the number of digits beyond the
1509        //decimal point in order to form an integer.  If the        //decimal point in order to form an integer.  If the
1510        //exponent is zero, there must be no digits after the        //exponent is zero, there must be no digits after the
1511        //decimal point.  If the exponent is negative, there        //decimal point.  If the exponent is negative, there
1512        //must be no digits after the decimal point, and the        //must be no digits after the decimal point, and the
1513        //trailing zeros on the part before the decimal point        //trailing zeros on the part before the decimal point
1514        //must be adequate to handle the right decimal shift.        //must be adequate to handle the right decimal shift.
1515        if (exponent_val == 0)        if (exponent_val == 0)
1516           {           {
1517           if (mant_adp_len)           if (mant_adp_len)
1518              {              {
1519              *failure = 1;              *failure = 1;
1520              return;              return;
1521              }              }
1522           }           }
1523        else if (exponent_val > 0)        else if (exponent_val > 0)
1524           {           {
1525           if ((int)mant_adp_len > exponent_val)           if ((int)mant_adp_len > exponent_val)
1526              {              {
1527              *failure = 1;              *failure = 1;
1528              return;              return;
1529              }              }
1530           }           }
1531        else //exponent_val < 0        else //exponent_val < 0
1532           {           {
1533           if (mant_adp_len)           if (mant_adp_len)
1534              {              {
1535              *failure = 1;              *failure = 1;
1536              return;              return;
1537              }              }
1538           else           else
1539              {              {
1540              //Count the number of trailing zeros on the part              //Count the number of trailing zeros on the part
1541              //before the decimal point.              //before the decimal point.
1542              size_t trailing_zero_count;              size_t trailing_zero_count;
1543              int idx;              int idx;
1544    
1545              trailing_zero_count = 0;              trailing_zero_count = 0;
1546                    
1547              for(idx = mant_bdp + mant_bdp_len - 1;              for(idx = mant_bdp + mant_bdp_len - 1;
1548                  (mant_bdp_len != 0) && (idx >= (int)mant_bdp);                  (mant_bdp_len != 0) && (idx >= (int)mant_bdp);
1549                  idx--)                  idx--)
1550                  {                  {
1551                  if (s[idx] == '0')                  if (s[idx] == '0')
1552                     trailing_zero_count++;                     trailing_zero_count++;
1553                  else                  else
1554                     break;                     break;
1555                  }                  }
1556                            
1557              //Check on the assertion about trailing zeros.              //Check on the assertion about trailing zeros.
1558              if ((int)trailing_zero_count < -exponent_val)              if ((int)trailing_zero_count < -exponent_val)
1559                 {                 {
1560                 *failure = 1;                 *failure = 1;
1561                 return;                 return;
1562                 }                 }
1563              }              }
1564           }           }
1565            
1566           {           {
1567           //Create a string long enough to hold the digits           //Create a string long enough to hold the digits
1568           //before the decimal point plus the ones after and           //before the decimal point plus the ones after and
1569           //convert that to an arbitrary integer.           //convert that to an arbitrary integer.
1570           //Form a power of 10 which is 10 exponentiated to           //Form a power of 10 which is 10 exponentiated to
1571           //the absolute value of the exponent.  If the           //the absolute value of the exponent.  If the
1572           //exponent was positive, multiply by it.  If the           //exponent was positive, multiply by it.  If the
1573           //exponent was negative, divide by it.           //exponent was negative, divide by it.
1574           char *conv_str;           char *conv_str;
1575           size_t sidx;           size_t sidx;
1576           GMP_INTS_mpz_struct power_of_ten, k10, trash;           GMP_INTS_mpz_struct power_of_ten, k10, trash;
1577    
1578           GMP_INTS_mpz_init(&power_of_ten);           GMP_INTS_mpz_init(&power_of_ten);
1579           GMP_INTS_mpz_init(&k10);           GMP_INTS_mpz_init(&k10);
1580           GMP_INTS_mpz_init(&trash);           GMP_INTS_mpz_init(&trash);
1581    
1582           conv_str = GMP_INTS_malloc(sizeof(char) * (mant_bdp_len + mant_adp_len + 1));           conv_str = GMP_INTS_malloc(sizeof(char) * (mant_bdp_len + mant_adp_len + 1));
1583    
1584           sidx=0;           sidx=0;
1585    
1586           for (si=mant_bdp; si<(mant_bdp+mant_bdp_len); si++)           for (si=mant_bdp; si<(mant_bdp+mant_bdp_len); si++)
1587              {              {
1588              conv_str[sidx] = s[si];              conv_str[sidx] = s[si];
1589              sidx++;              sidx++;
1590              }              }
1591           for (si=mant_adp; si<(mant_adp+mant_adp_len); si++)           for (si=mant_adp; si<(mant_adp+mant_adp_len); si++)
1592              {              {
1593              conv_str[sidx] = s[si];              conv_str[sidx] = s[si];
1594              sidx++;              sidx++;
1595              }              }
1596           conv_str[sidx] = 0;           conv_str[sidx] = 0;
1597    
1598           assert(sidx == (mant_bdp_len + mant_adp_len));           assert(sidx == (mant_bdp_len + mant_adp_len));
1599                    
1600           GMP_INTS_mpz_set_simple_char_str(z, conv_str);           GMP_INTS_mpz_set_simple_char_str(z, conv_str);
1601    
1602           GMP_INTS_mpz_set_ui(&k10, 10);           GMP_INTS_mpz_set_ui(&k10, 10);
1603    
1604           if (exponent_val > 0)           if (exponent_val > 0)
1605              GMP_INTS_mpz_pow_ui(&power_of_ten, &k10, exponent_val-mant_adp_len);              GMP_INTS_mpz_pow_ui(&power_of_ten, &k10, exponent_val-mant_adp_len);
1606           else           else
1607              GMP_INTS_mpz_pow_ui(&power_of_ten, &k10, -exponent_val);              GMP_INTS_mpz_pow_ui(&power_of_ten, &k10, -exponent_val);
1608    
1609           if (exponent_val >= 0)           if (exponent_val >= 0)
1610              {              {
1611              GMP_INTS_mpz_mul(z, z, &power_of_ten);              GMP_INTS_mpz_mul(z, z, &power_of_ten);
1612              }              }
1613           else           else
1614              {              {
1615              GMP_INTS_mpz_tdiv_qr (&k10,              GMP_INTS_mpz_tdiv_qr (&k10,
1616                                    &trash,                                    &trash,
1617                                    z,                                    z,
1618                                    &power_of_ten);                                    &power_of_ten);
1619              GMP_INTS_mpz_copy(z, &k10);              GMP_INTS_mpz_copy(z, &k10);
1620              }              }
1621    
1622           //If the argument had a minus sign, invert.           //If the argument had a minus sign, invert.
1623           if (mant_sign == '-')           if (mant_sign == '-')
1624              z->size = -z->size;              z->size = -z->size;
1625    
1626           GMP_INTS_free(conv_str);           GMP_INTS_free(conv_str);
1627    
1628           GMP_INTS_mpz_clear(&trash);           GMP_INTS_mpz_clear(&trash);
1629           GMP_INTS_mpz_clear(&k10);           GMP_INTS_mpz_clear(&k10);
1630           GMP_INTS_mpz_clear(&power_of_ten);           GMP_INTS_mpz_clear(&power_of_ten);
1631    
1632           //Finally, if the arbitrary integer has overflowed, this is           //Finally, if the arbitrary integer has overflowed, this is
1633           //a parse failure.  Must declare as such.           //a parse failure.  Must declare as such.
1634           if (z->flags)           if (z->flags)
1635              *failure = 1;              *failure = 1;
1636           }           }
1637        }        }
1638     else     else
1639        {        {
1640        *failure = 1;        *failure = 1;
1641        return;        return;
1642        }        }
1643     }     }
1644    
1645    
1646  void GMP_INTS_mpz_set_general_int(GMP_INTS_mpz_struct *z,  void GMP_INTS_mpz_set_general_int(GMP_INTS_mpz_struct *z,
1647                                    int *failure,                                    int *failure,
1648                                    const char *s)                                    const char *s)
1649     {     {
1650     //Eyeball the input parameters.     //Eyeball the input parameters.
1651     assert(z != NULL);     assert(z != NULL);
1652     assert(z->n_allocd > 0);     assert(z->n_allocd > 0);
1653     assert(z->limbs != NULL);     assert(z->limbs != NULL);
1654     assert(failure != NULL);     assert(failure != NULL);
1655     assert(s != NULL);     assert(s != NULL);
1656    
1657     //Try to parse it as a simple integer.     //Try to parse it as a simple integer.
1658     if (BSTRFUNC_is_sint_wo_commas(s))     if (BSTRFUNC_is_sint_wo_commas(s))
1659        {        {
1660        GMP_INTS_mpz_set_simple_char_str(z, s);        GMP_INTS_mpz_set_simple_char_str(z, s);
1661        *failure = 0;        *failure = 0;
1662        return;        return;
1663        }        }
1664     //If that didn't work, try to parse it as a simple     //If that didn't work, try to parse it as a simple
1665     //integer with commas.     //integer with commas.
1666     else if (BSTRFUNC_is_sint_w_commas(s))     else if (BSTRFUNC_is_sint_w_commas(s))
1667        {        {
1668        GMP_INTS_mpz_set_simple_char_str(z, s);        GMP_INTS_mpz_set_simple_char_str(z, s);
1669        *failure = 0;        *failure = 0;
1670        return;        return;
1671        }        }
1672     //If neither of those worked, try to parse it as     //If neither of those worked, try to parse it as
1673     //something containing scientific notation.     //something containing scientific notation.
1674     else     else
1675        {        {
1676        GMP_INTS_mpz_set_sci_not_num(z, failure, s);        GMP_INTS_mpz_set_sci_not_num(z, failure, s);
1677    
1678        if (!*failure)        if (!*failure)
1679           {           {
1680           //We were able to parse it that way.           //We were able to parse it that way.
1681           //Everything is set up, just return.           //Everything is set up, just return.
1682           return;           return;
1683           }           }
1684        else        else
1685           {           {
1686           //We're out of options.  All parsing failed.           //We're out of options.  All parsing failed.
1687           GMP_INTS_mpz_set_ui(z, 0);           GMP_INTS_mpz_set_ui(z, 0);
1688           *failure = 1;           *failure = 1;
1689           return;           return;
1690           }           }
1691        }        }
1692     }     }
1693    
1694    
1695  void GMP_INTS_mpz_parse_into_uint32(unsigned *result,  void GMP_INTS_mpz_parse_into_uint32(unsigned *result,
1696                                      int *failure,                                      int *failure,
1697                                      char *s)                                      char *s)
1698     {     {
1699     GMP_INTS_mpz_struct arb_int;     GMP_INTS_mpz_struct arb_int;
1700    
1701     //Eyeball the input parameters.     //Eyeball the input parameters.
1702     assert(result != NULL);     assert(result != NULL);
1703     assert(failure != NULL);     assert(failure != NULL);
1704     assert(s != NULL);     assert(s != NULL);
1705    
1706     //Allocate space for the one arbitrary integer we need.     //Allocate space for the one arbitrary integer we need.
1707     GMP_INTS_mpz_init(&arb_int);     GMP_INTS_mpz_init(&arb_int);
1708    
1709     //Try to parse the string into an arbitrary length integer     //Try to parse the string into an arbitrary length integer
1710     //using all methods known to man.     //using all methods known to man.
1711     GMP_INTS_mpz_set_general_int(&arb_int, failure, s);     GMP_INTS_mpz_set_general_int(&arb_int, failure, s);
1712    
1713     //If the parse failed, we must declare failure and return     //If the parse failed, we must declare failure and return
1714     //0.     //0.
1715     if (*failure)     if (*failure)
1716        {        {
1717        *result = 0;        *result = 0;
1718        *failure = 1;        *failure = 1;
1719        }        }
1720     else     else
1721        {        {
1722        //We might have success, but it might be negative or        //We might have success, but it might be negative or
1723        //too big.        //too big.
1724        if (arb_int.size == 1)        if (arb_int.size == 1)
1725           {           {
1726           *result = arb_int.limbs[0];           *result = arb_int.limbs[0];
1727           *failure = 0;           *failure = 0;
1728           }           }
1729        else if (arb_int.size == 0)        else if (arb_int.size == 0)
1730           {           {
1731           *result =  0;           *result =  0;
1732           *failure = 0;           *failure = 0;
1733           }           }
1734        else        else
1735           {           {
1736           *result = 0;           *result = 0;
1737           *failure = 1;           *failure = 1;
1738           }           }
1739        }        }
1740    
1741     //Deallocate the arbitrary integer.     //Deallocate the arbitrary integer.
1742     GMP_INTS_mpz_clear(&arb_int);     GMP_INTS_mpz_clear(&arb_int);
1743     }     }
1744    
1745    
1746  void GMP_INTS_mpz_swap(GMP_INTS_mpz_struct *a,  void GMP_INTS_mpz_swap(GMP_INTS_mpz_struct *a,
1747                         GMP_INTS_mpz_struct *b)                         GMP_INTS_mpz_struct *b)
1748     {     {
1749     GMP_INTS_mpz_struct temp;     GMP_INTS_mpz_struct temp;
1750    
1751     //Eyeball the input parameters.     //Eyeball the input parameters.
1752     assert(a           != NULL);     assert(a           != NULL);
1753     assert(a->n_allocd >  0);     assert(a->n_allocd >  0);
1754     assert(a->limbs    != NULL);     assert(a->limbs    != NULL);
1755     assert(b           != NULL);     assert(b           != NULL);
1756     assert(b->n_allocd >  0);     assert(b->n_allocd >  0);
1757     assert(b->limbs    != NULL);     assert(b->limbs    != NULL);
1758    
1759     //Make the swap via memory copy.     //Make the swap via memory copy.
1760     memcpy(&temp, a,     sizeof(GMP_INTS_mpz_struct));     memcpy(&temp, a,     sizeof(GMP_INTS_mpz_struct));
1761     memcpy(a,     b,     sizeof(GMP_INTS_mpz_struct));     memcpy(a,     b,     sizeof(GMP_INTS_mpz_struct));
1762     memcpy(b,     &temp, sizeof(GMP_INTS_mpz_struct));     memcpy(b,     &temp, sizeof(GMP_INTS_mpz_struct));
1763     }     }
1764    
1765    
1766  /******************************************************************/  /******************************************************************/
1767  /***  PUBLIC ARITHMETIC FUNCTIONS   *******************************/  /***  PUBLIC ARITHMETIC FUNCTIONS   *******************************/
1768  /******************************************************************/  /******************************************************************/
1769    
1770  //07/15/01:  Unit test and visual inspection passed.  //07/15/01:  Unit test and visual inspection passed.
1771  void GMP_INTS_mpz_add (      GMP_INTS_mpz_struct *w,  void GMP_INTS_mpz_add (      GMP_INTS_mpz_struct *w,
1772                         const GMP_INTS_mpz_struct *u,                         const GMP_INTS_mpz_struct *u,
1773                         const GMP_INTS_mpz_struct *v)                         const GMP_INTS_mpz_struct *v)
1774     {     {
1775     GMP_INTS_limb_srcptr up, vp;     GMP_INTS_limb_srcptr up, vp;
1776     GMP_INTS_limb_ptr wp;     GMP_INTS_limb_ptr wp;
1777     GMP_INTS_size_t usize, vsize, wsize;     GMP_INTS_size_t usize, vsize, wsize;
1778     GMP_INTS_size_t abs_usize;     GMP_INTS_size_t abs_usize;
1779     GMP_INTS_size_t abs_vsize;     GMP_INTS_size_t abs_vsize;
1780    
1781     //Look at the input parameters carefully.     //Look at the input parameters carefully.
1782     assert(w != NULL);     assert(w != NULL);
1783     assert(u != NULL);     assert(u != NULL);
1784     assert(v != NULL);     assert(v != NULL);
1785     assert(w->n_allocd > 0);     assert(w->n_allocd > 0);
1786     assert(u->n_allocd > 0);     assert(u->n_allocd > 0);
1787     assert(v->n_allocd > 0);     assert(v->n_allocd > 0);
1788     assert(w->limbs != NULL);     assert(w->limbs != NULL);
1789     assert(u->limbs != NULL);     assert(u->limbs != NULL);
1790     assert(v->limbs != NULL);     assert(v->limbs != NULL);
1791    
1792     //Handle the case of a tainted result.  If either of the     //Handle the case of a tainted result.  If either of the
1793     //two inputs are either direct overflows or tainted by     //two inputs are either direct overflows or tainted by
1794     //an overflow, mark the result tainted and do not perform     //an overflow, mark the result tainted and do not perform
1795     //any arithmetic operation.     //any arithmetic operation.
1796        {        {
1797        int taint;        int taint;
1798    
1799        taint = GMP_INTS_two_op_flags_map(u->flags, v->flags);        taint = GMP_INTS_two_op_flags_map(u->flags, v->flags);
1800    
1801        w->flags = 0;        w->flags = 0;
1802           //"w" starts off with a clean slate.  Must do this           //"w" starts off with a clean slate.  Must do this
1803           //after taint calculation in case locations of u or v           //after taint calculation in case locations of u or v
1804           //are the same as w.           //are the same as w.
1805        if (taint)        if (taint)
1806           {           {
1807           w->flags = taint;           w->flags = taint;
1808           return;           return;
1809           }           }
1810        }        }
1811    
1812     usize = u->size;     usize = u->size;
1813     vsize = v->size;     vsize = v->size;
1814     abs_usize = GMP_INTS_abs_of_size_t(usize);     abs_usize = GMP_INTS_abs_of_size_t(usize);
1815     abs_vsize = GMP_INTS_abs_of_size_t(vsize);     abs_vsize = GMP_INTS_abs_of_size_t(vsize);
1816    
1817     //Arrange things so that U has at least as many     //Arrange things so that U has at least as many
1818     //limbs as V, i.e. limbs(U) >= limbs(V);     //limbs as V, i.e. limbs(U) >= limbs(V);
1819     if (abs_usize < abs_vsize)     if (abs_usize < abs_vsize)
1820        {        {
1821        const GMP_INTS_mpz_struct *tmp_ptr;        const GMP_INTS_mpz_struct *tmp_ptr;
1822        GMP_INTS_size_t      tmp_size;        GMP_INTS_size_t      tmp_size;
1823    
1824        //Swap U and V.  This does no harm, because we are        //Swap U and V.  This does no harm, because we are
1825        //manipulating only local variables.  This does not        //manipulating only local variables.  This does not
1826        //affect the caller.        //affect the caller.
1827        tmp_ptr = u;        tmp_ptr = u;
1828        u = v;        u = v;
1829        v = tmp_ptr;        v = tmp_ptr;
1830        tmp_size = usize;        tmp_size = usize;
1831        usize = vsize;        usize = vsize;
1832        vsize = tmp_size;        vsize = tmp_size;
1833        tmp_size = abs_usize;        tmp_size = abs_usize;
1834        abs_usize = abs_vsize;        abs_usize = abs_vsize;
1835        abs_vsize = tmp_size;        abs_vsize = tmp_size;
1836        }        }
1837    
1838     /* True: ABS_USIZE >= ABS_VSIZE.  */     /* True: ABS_USIZE >= ABS_VSIZE.  */
1839    
1840     /* If not space for w (and possible carry), increase space.  */     /* If not space for w (and possible carry), increase space.  */
1841     wsize = abs_usize + 1;     wsize = abs_usize + 1;
1842     if (w->n_allocd < wsize)     if (w->n_allocd < wsize)
1843        GMP_INTS_mpz_realloc(w, wsize);        GMP_INTS_mpz_realloc(w, wsize);
1844    
1845     //These pointers must be obtained after realloc.  At this point,     //These pointers must be obtained after realloc.  At this point,
1846     //u or v may be the same as w.     //u or v may be the same as w.
1847     up = u->limbs;     up = u->limbs;
1848     vp = v->limbs;     vp = v->limbs;
1849     wp = w->limbs;     wp = w->limbs;
1850    
1851     if ((usize ^ vsize) < 0)     if ((usize ^ vsize) < 0)
1852        {        {
1853        //U and V have different sign.  Need to compare them to determine        //U and V have different sign.  Need to compare them to determine
1854        //which operand to subtract from which.        //which operand to subtract from which.
1855    
1856        //This test is right since ABS_USIZE >= ABS_VSIZE.        //This test is right since ABS_USIZE >= ABS_VSIZE.
1857        //If the equality case is ruled out, then U has more limbs        //If the equality case is ruled out, then U has more limbs
1858        //than V, which means that it is bigger in magnitude.        //than V, which means that it is bigger in magnitude.
1859        if (abs_usize != abs_vsize)        if (abs_usize != abs_vsize)
1860           {           {
1861           GMP_INTS_mpn_sub (wp, up, abs_usize, vp, abs_vsize);           GMP_INTS_mpn_sub (wp, up, abs_usize, vp, abs_vsize);
1862           wsize = abs_usize;           wsize = abs_usize;
1863    
1864           //Normalize the result.  This was formerly a macro.           //Normalize the result.  This was formerly a macro.
1865           //To normalize in this context means to trim the size           //To normalize in this context means to trim the size
1866           //down to eliminate any leading zero limbs that came           //down to eliminate any leading zero limbs that came
1867           //about because the size of the result of an operation           //about because the size of the result of an operation
1868           //was overestimated.           //was overestimated.
1869           GMP_INTS_mpn_normalize(wp, &wsize);           GMP_INTS_mpn_normalize(wp, &wsize);
1870    
1871           if (usize < 0)           if (usize < 0)
1872              wsize = -wsize;              wsize = -wsize;
1873           }           }
1874        else if (GMP_INTS_mpn_cmp (up, vp, abs_usize) < 0)        else if (GMP_INTS_mpn_cmp (up, vp, abs_usize) < 0)
1875           {           {
1876           GMP_INTS_mpn_sub_n (wp, vp, up, abs_usize);           GMP_INTS_mpn_sub_n (wp, vp, up, abs_usize);
1877           wsize = abs_usize;           wsize = abs_usize;
1878           GMP_INTS_mpn_normalize(wp, &wsize);           GMP_INTS_mpn_normalize(wp, &wsize);
1879           if (usize >= 0)           if (usize >= 0)
1880              wsize = -wsize;              wsize = -wsize;
1881           }           }
1882        else        else
1883           {           {
1884           GMP_INTS_mpn_sub_n (wp, up, vp, abs_usize);           GMP_INTS_mpn_sub_n (wp, up, vp, abs_usize);
1885           wsize = abs_usize;           wsize = abs_usize;
1886           GMP_INTS_mpn_normalize(wp, &wsize);           GMP_INTS_mpn_normalize(wp, &wsize);
1887           if (usize < 0)           if (usize < 0)
1888              wsize = -wsize;              wsize = -wsize;
1889           }           }
1890        }        }
1891     else     else
1892        {        {
1893        /* U and V have same sign.  Add them.  */        /* U and V have same sign.  Add them.  */
1894        GMP_INTS_limb_t cy_limb        GMP_INTS_limb_t cy_limb
1895           = GMP_INTS_mpn_add (wp, up, abs_usize, vp, abs_vsize);           = GMP_INTS_mpn_add (wp, up, abs_usize, vp, abs_vsize);
1896        wp[abs_usize] = cy_limb;        wp[abs_usize] = cy_limb;
1897        wsize = abs_usize + cy_limb;        wsize = abs_usize + cy_limb;
1898        if (usize < 0)        if (usize < 0)
1899           wsize = -wsize;           wsize = -wsize;
1900        }        }
1901    
1902     w->size = wsize;     w->size = wsize;
1903    
1904     //Handle the case of an overflowed result.  If the result     //Handle the case of an overflowed result.  If the result
1905     //of the addition is too big or too small, mark it as     //of the addition is too big or too small, mark it as
1906     //overflowed.     //overflowed.
1907     if (w->size > GMP_INTS_MAXIMUM_LIMBS_PER_INT)     if (w->size > GMP_INTS_MAXIMUM_LIMBS_PER_INT)
1908        {        {
1909        w->flags = GMP_INTS_EF_INTOVF_POS;        w->flags = GMP_INTS_EF_INTOVF_POS;
1910        }        }
1911     else if (w->size < -GMP_INTS_MAXIMUM_LIMBS_PER_INT)     else if (w->size < -GMP_INTS_MAXIMUM_LIMBS_PER_INT)
1912        {        {
1913        w->flags = GMP_INTS_EF_INTOVF_NEG;        w->flags = GMP_INTS_EF_INTOVF_NEG;
1914        }        }
1915     }     }
1916    
1917    
1918  //07/15/01:  Unit testing skipped because of recursive  //07/15/01:  Unit testing skipped because of recursive
1919  //nature.  Visual inspection OK.  //nature.  Visual inspection OK.
1920  void GMP_INTS_mpz_add_ui (      GMP_INTS_mpz_struct *w,  void GMP_INTS_mpz_add_ui (      GMP_INTS_mpz_struct *w,
1921                            const GMP_INTS_mpz_struct *u,                            const GMP_INTS_mpz_struct *u,
1922                            unsigned long int          v)                            unsigned long int          v)
1923     {     {
1924     //The GNU MP version of this is quite efficient, and this     //The GNU MP version of this is quite efficient, and this
1925     //makes sense since it is a common operation.  However,     //makes sense since it is a common operation.  However,
1926     //for simplicity just define this recursively in terms     //for simplicity just define this recursively in terms
1927     //of the ADD function.  This can always be made quicker     //of the ADD function.  This can always be made quicker
1928     //later (by changing back to the GNU MP version).     //later (by changing back to the GNU MP version).
1929     GMP_INTS_mpz_struct temp;     GMP_INTS_mpz_struct temp;
1930    
1931     //Eyeball the inputs carefully.     //Eyeball the inputs carefully.
1932     assert(w != NULL);     assert(w != NULL);
1933     assert(w->n_allocd > 0);     assert(w->n_allocd > 0);
1934     assert(w->limbs != NULL);     assert(w->limbs != NULL);
1935     assert(u != NULL);     assert(u != NULL);
1936     assert(u->n_allocd > 0);     assert(u->n_allocd > 0);
1937     assert(u->limbs != NULL);     assert(u->limbs != NULL);
1938    
1939     //Create a temporary integer.     //Create a temporary integer.
1940     GMP_INTS_mpz_init(&temp);     GMP_INTS_mpz_init(&temp);
1941    
1942     //Set the temporary integer to the value of the input     //Set the temporary integer to the value of the input
1943     //argument.     //argument.
1944     GMP_INTS_mpz_set_ui(&temp, v);     GMP_INTS_mpz_set_ui(&temp, v);
1945    
1946     //Do the actual addition.  This recursive definition     //Do the actual addition.  This recursive definition
1947     //is inherently wasteful, but I'm after clarity, not     //is inherently wasteful, but I'm after clarity, not
1948     //warp speed.     //warp speed.
1949     GMP_INTS_mpz_add(w, u, &temp);     GMP_INTS_mpz_add(w, u, &temp);
1950    
1951     //Destroy the temporary integer (this will reclaim the     //Destroy the temporary integer (this will reclaim the
1952     //memory).     //memory).
1953     GMP_INTS_mpz_clear(&temp);     GMP_INTS_mpz_clear(&temp);
1954     }     }
1955    
1956    
1957  //07/15/01:  Visual inspection passed.  Not unit tested  //07/15/01:  Visual inspection passed.  Not unit tested
1958  //because of symmetry with GMP_INTS_mpz_add().  //because of symmetry with GMP_INTS_mpz_add().
1959  void GMP_INTS_mpz_sub (      GMP_INTS_mpz_struct *w,  void GMP_INTS_mpz_sub (      GMP_INTS_mpz_struct *w,
1960                         const GMP_INTS_mpz_struct *u,                         const GMP_INTS_mpz_struct *u,
1961                         const GMP_INTS_mpz_struct *v)                         const GMP_INTS_mpz_struct *v)
1962     {     {
1963     GMP_INTS_limb_srcptr up, vp;     GMP_INTS_limb_srcptr up, vp;
1964     GMP_INTS_limb_ptr wp;     GMP_INTS_limb_ptr wp;
1965     GMP_INTS_size_t usize, vsize, wsize;     GMP_INTS_size_t usize, vsize, wsize;
1966     GMP_INTS_size_t abs_usize;     GMP_INTS_size_t abs_usize;
1967     GMP_INTS_size_t abs_vsize;     GMP_INTS_size_t abs_vsize;
1968    
1969     //Look at the input parameters carefully.     //Look at the input parameters carefully.
1970     assert(w != NULL);     assert(w != NULL);
1971     assert(u != NULL);     assert(u != NULL);
1972     assert(v != NULL);     assert(v != NULL);
1973     assert(w->n_allocd > 0);     assert(w->n_allocd > 0);
1974     assert(u->n_allocd > 0);     assert(u->n_allocd > 0);
1975     assert(v->n_allocd > 0);     assert(v->n_allocd > 0);
1976     assert(w->limbs != NULL);     assert(w->limbs != NULL);
1977     assert(u->limbs != NULL);     assert(u->limbs != NULL);
1978     assert(v->limbs != NULL);     assert(v->limbs != NULL);
1979    
1980     //Handle the case of a tainted result.  If either of the     //Handle the case of a tainted result.  If either of the
1981     //two inputs are either direct overflows or tainted by     //two inputs are either direct overflows or tainted by
1982     //an overflow, mark the result tainted and do not perform     //an overflow, mark the result tainted and do not perform
1983     //any arithmetic operation.     //any arithmetic operation.
1984        {        {
1985        int taint;        int taint;
1986    
1987        taint = GMP_INTS_two_op_flags_map(u->flags, v->flags);        taint = GMP_INTS_two_op_flags_map(u->flags, v->flags);
1988    
1989        w->flags = 0;        w->flags = 0;
1990           //"w" starts off with a clean slate.  Must do this           //"w" starts off with a clean slate.  Must do this
1991           //after taint calculation in case locations of u or v           //after taint calculation in case locations of u or v
1992           //are the same as w.           //are the same as w.
1993        if (taint)        if (taint)
1994           {           {
1995           w->flags = taint;           w->flags = taint;
1996           return;           return;
1997           }           }
1998        }        }
1999    
2000     usize = u->size;     usize = u->size;
2001     vsize = -(v->size);  /* The "-" makes the difference from mpz_add */     vsize = -(v->size);  /* The "-" makes the difference from mpz_add */
2002     abs_usize = GMP_INTS_abs_of_size_t(usize);     abs_usize = GMP_INTS_abs_of_size_t(usize);
2003     abs_vsize = GMP_INTS_abs_of_size_t(vsize);     abs_vsize = GMP_INTS_abs_of_size_t(vsize);
2004    
2005     if (abs_usize < abs_vsize)     if (abs_usize < abs_vsize)
2006        {        {
2007        const GMP_INTS_mpz_struct *tmp_ptr;        const GMP_INTS_mpz_struct *tmp_ptr;
2008        GMP_INTS_size_t      tmp_size;        GMP_INTS_size_t      tmp_size;
2009    
2010        //Swap U and V.  This does no harm, because we are        //Swap U and V.  This does no harm, because we are
2011        //manipulating only local variables.  This does not        //manipulating only local variables.  This does not
2012        //affect the caller.        //affect the caller.
2013        tmp_ptr = u;        tmp_ptr = u;
2014        u = v;        u = v;
2015        v = tmp_ptr;        v = tmp_ptr;
2016        tmp_size = usize;        tmp_size = usize;
2017        usize = vsize;        usize = vsize;
2018        vsize = tmp_size;        vsize = tmp_size;
2019        tmp_size = abs_usize;        tmp_size = abs_usize;
2020        abs_usize = abs_vsize;        abs_usize = abs_vsize;
2021        abs_vsize = tmp_size;        abs_vsize = tmp_size;
2022        }        }
2023    
2024     /* True: ABS_USIZE >= ABS_VSIZE.  */     /* True: ABS_USIZE >= ABS_VSIZE.  */
2025    
2026     /* If not space for w (and possible carry), increase space.  */     /* If not space for w (and possible carry), increase space.  */
2027     wsize = abs_usize + 1;     wsize = abs_usize + 1;
2028     if (w->n_allocd < wsize)     if (w->n_allocd < wsize)
2029        GMP_INTS_mpz_realloc (w, wsize);        GMP_INTS_mpz_realloc (w, wsize);
2030    
2031     /* These must be after realloc (u or v may be the same as w).  */     /* These must be after realloc (u or v may be the same as w).  */
2032     up = u->limbs;     up = u->limbs;
2033     vp = v->limbs;     vp = v->limbs;
2034     wp = w->limbs;     wp = w->limbs;
2035    
2036     if ((usize ^ vsize) < 0)     if ((usize ^ vsize) < 0)
2037        {        {
2038        //U and V have different sign.  Need to compare them to determine        //U and V have different sign.  Need to compare them to determine
2039        //which operand to subtract from which.        //which operand to subtract from which.
2040    
2041        //This test is right since ABS_USIZE >= ABS_VSIZE.        //This test is right since ABS_USIZE >= ABS_VSIZE.
2042        if (abs_usize != abs_vsize)        if (abs_usize != abs_vsize)
2043           {           {
2044           GMP_INTS_mpn_sub (wp, up, abs_usize, vp, abs_vsize);           GMP_INTS_mpn_sub (wp, up, abs_usize, vp, abs_vsize);
2045           wsize = abs_usize;           wsize = abs_usize;
2046           GMP_INTS_mpn_normalize(wp, &wsize);           GMP_INTS_mpn_normalize(wp, &wsize);
2047           if (usize < 0)           if (usize < 0)
2048              wsize = -wsize;              wsize = -wsize;
2049           }           }
2050        else if (GMP_INTS_mpn_cmp (up, vp, abs_usize) < 0)        else if (GMP_INTS_mpn_cmp (up, vp, abs_usize) < 0)
2051           {           {
2052           GMP_INTS_mpn_sub_n (wp, vp, up, abs_usize);           GMP_INTS_mpn_sub_n (wp, vp, up, abs_usize);
2053           wsize = abs_usize;           wsize = abs_usize;
2054           GMP_INTS_mpn_normalize(wp, &wsize);           GMP_INTS_mpn_normalize(wp, &wsize);
2055           if (usize >= 0)           if (usize >= 0)
2056              wsize = -wsize;              wsize = -wsize;
2057           }           }
2058        else        else
2059           {           {
2060           GMP_INTS_mpn_sub_n (wp, up, vp, abs_usize);           GMP_INTS_mpn_sub_n (wp, up, vp, abs_usize);
2061           wsize = abs_usize;           wsize = abs_usize;
2062           GMP_INTS_mpn_normalize (wp, &wsize);           GMP_INTS_mpn_normalize (wp, &wsize);
2063           if (usize < 0)           if (usize < 0)
2064              wsize = -wsize;              wsize = -wsize;
2065           }           }
2066        }        }
2067     else     else
2068        {        {
2069        /* U and V have same sign.  Add them.  */        /* U and V have same sign.  Add them.  */
2070        GMP_INTS_limb_t cy_limb        GMP_INTS_limb_t cy_limb
2071           = GMP_INTS_mpn_add (wp, up, abs_usize, vp, abs_vsize);           = GMP_INTS_mpn_add (wp, up, abs_usize, vp, abs_vsize);
2072        wp[abs_usize] = cy_limb;        wp[abs_usize] = cy_limb;
2073        wsize = abs_usize + cy_limb;        wsize = abs_usize + cy_limb;
2074        if (usize < 0)        if (usize < 0)
2075           wsize = -wsize;           wsize = -wsize;
2076        }        }
2077    
2078     w->size = wsize;     w->size = wsize;
2079    
2080     //Handle the case of an overflowed result.  If the result     //Handle the case of an overflowed result.  If the result
2081     //of the addition is too big or too small, mark it as     //of the addition is too big or too small, mark it as
2082     //overflowed.     //overflowed.
2083     if (w->size > GMP_INTS_MAXIMUM_LIMBS_PER_INT)     if (w->size > GMP_INTS_MAXIMUM_LIMBS_PER_INT)
2084        {        {
2085        w->flags = GMP_INTS_EF_INTOVF_POS;        w->flags = GMP_INTS_EF_INTOVF_POS;
2086        }        }
2087     else if (w->size < -GMP_INTS_MAXIMUM_LIMBS_PER_INT)     else if (w->size < -GMP_INTS_MAXIMUM_LIMBS_PER_INT)
2088        {        {
2089        w->flags = GMP_INTS_EF_INTOVF_NEG;        w->flags = GMP_INTS_EF_INTOVF_NEG;
2090        }        }
2091     }     }
2092    
2093    
2094  //07/15/01:  Unit testing skipped because of recursive  //07/15/01:  Unit testing skipped because of recursive
2095  //nature.  Visual inspection OK.  //nature.  Visual inspection OK.
2096  void GMP_INTS_mpz_sub_ui (      GMP_INTS_mpz_struct *w,  void GMP_INTS_mpz_sub_ui (      GMP_INTS_mpz_struct *w,
2097                            const GMP_INTS_mpz_struct *u,                            const GMP_INTS_mpz_struct *u,
2098                            unsigned long int          v)                            unsigned long int          v)
2099     {     {
2100     //The GNU MP version of this is quite efficient, and this     //The GNU MP version of this is quite efficient, and this
2101     //makes sense since it is a common operation.  However,     //makes sense since it is a common operation.  However,
2102     //for simplicity just define this recursively in terms     //for simplicity just define this recursively in terms
2103     //of the SUB function.  This can always be made quicker     //of the SUB function.  This can always be made quicker
2104     //later (by changing back to the GNU MP version).     //later (by changing back to the GNU MP version).
2105     GMP_INTS_mpz_struct temp;     GMP_INTS_mpz_struct temp;
2106    
2107     //Eyeball the inputs carefully.     //Eyeball the inputs carefully.
2108     assert(w != NULL);     assert(w != NULL);
2109     assert(w->n_allocd > 0);     assert(w->n_allocd > 0);
2110     assert(w->limbs != NULL);     assert(w->limbs != NULL);
2111     assert(u != NULL);     assert(u != NULL);
2112     assert(u->n_allocd > 0);     assert(u->n_allocd > 0);
2113     assert(u->limbs != NULL);     assert(u->limbs != NULL);
2114    
2115     //Create a temporary integer.     //Create a temporary integer.
2116     GMP_INTS_mpz_init(&temp);     GMP_INTS_mpz_init(&temp);
2117    
2118     //Set the temporary integer to the value of the input     //Set the temporary integer to the value of the input
2119     //argument.     //argument.
2120     GMP_INTS_mpz_set_ui(&temp, v);     GMP_INTS_mpz_set_ui(&temp, v);
2121    
2122     //Do the actual subtraction.  This recursive definition     //Do the actual subtraction.  This recursive definition
2123     //is inherently wasteful, but I'm after clarity, not     //is inherently wasteful, but I'm after clarity, not
2124     //warp speed.     //warp speed.
2125     GMP_INTS_mpz_sub(w, u, &temp);     GMP_INTS_mpz_sub(w, u, &temp);
2126    
2127     //Destroy the temporary integer (this will reclaim the     //Destroy the temporary integer (this will reclaim the
2128     //memory).     //memory).
2129     GMP_INTS_mpz_clear(&temp);     GMP_INTS_mpz_clear(&temp);
2130     }     }
2131    
2132    
2133  void GMP_INTS_mpz_mul (      GMP_INTS_mpz_struct *w,  void GMP_INTS_mpz_mul (      GMP_INTS_mpz_struct *w,
2134                         const GMP_INTS_mpz_struct *u,                         const GMP_INTS_mpz_struct *u,
2135                         const GMP_INTS_mpz_struct *v)                         const GMP_INTS_mpz_struct *v)
2136     {     {
2137     GMP_INTS_size_t    usize = u->size;     GMP_INTS_size_t    usize = u->size;
2138     GMP_INTS_size_t    vsize = v->size;     GMP_INTS_size_t    vsize = v->size;
2139     GMP_INTS_size_t    wsize;     GMP_INTS_size_t    wsize;
2140     GMP_INTS_size_t    sign_product;     GMP_INTS_size_t    sign_product;
2141     GMP_INTS_limb_ptr  up, vp;     GMP_INTS_limb_ptr  up, vp;
2142     GMP_INTS_limb_ptr  wp;     GMP_INTS_limb_ptr  wp;
2143     GMP_INTS_limb_ptr  free_me = NULL;     GMP_INTS_limb_ptr  free_me = NULL;
2144     GMP_INTS_size_t    free_me_size;     GMP_INTS_size_t    free_me_size;
2145     GMP_INTS_limb_t    cy_limb;     GMP_INTS_limb_t    cy_limb;
2146    
2147     //Eyeball the inputs.     //Eyeball the inputs.
2148     assert(w != NULL);     assert(w != NULL);
2149     assert(w->n_allocd > 0);     assert(w->n_allocd > 0);
2150     assert(w->limbs != NULL);     assert(w->limbs != NULL);
2151     assert(u != NULL);     assert(u != NULL);
2152     assert(u->n_allocd > 0);     assert(u->n_allocd > 0);
2153     assert(u->limbs != NULL);     assert(u->limbs != NULL);
2154     assert(v != NULL);     assert(v != NULL);
2155     assert(v->n_allocd > 0);     assert(v->n_allocd > 0);
2156     assert(v->limbs != NULL);     assert(v->limbs != NULL);
2157    
2158     //Handle the case of a tainted result.  If either of the     //Handle the case of a tainted result.  If either of the
2159     //two inputs are either direct overflows or tainted by     //two inputs are either direct overflows or tainted by
2160     //an overflow, mark the result tainted and do not perform     //an overflow, mark the result tainted and do not perform
2161     //any arithmetic operation.     //any arithmetic operation.
2162        {        {
2163        int taint;        int taint;
2164    
2165        taint = GMP_INTS_two_op_flags_map(u->flags, v->flags);        taint = GMP_INTS_two_op_flags_map(u->flags, v->flags);
2166    
2167        w->flags = 0;        w->flags = 0;
2168           //"w" starts off with a clean slate.  Must do this           //"w" starts off with a clean slate.  Must do this
2169           //after taint calculation in case locations of u or v           //after taint calculation in case locations of u or v
2170           //are the same as w.           //are the same as w.
2171        if (taint)        if (taint)
2172           {           {
2173           w->flags = taint;           w->flags = taint;
2174           return;           return;
2175           }           }
2176        }        }
2177    
2178     sign_product = usize ^ vsize;     sign_product = usize ^ vsize;
2179     usize = GMP_INTS_abs_of_size_t(usize);     usize = GMP_INTS_abs_of_size_t(usize);
2180     vsize = GMP_INTS_abs_of_size_t(vsize);     vsize = GMP_INTS_abs_of_size_t(vsize);
2181    
2182     //Handle the case of a certain result overflow (why do the math when     //Handle the case of a certain result overflow (why do the math when
2183     //the result is certain?).  In general, when multiplying two inputs     //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     //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     //either M+N or M+N-1 limbs.  If M+N-1 > MAX_ALLOWED, then can declare
2186     //an early overflow.     //an early overflow.
2187     if ((usize + vsize - 1) > GMP_INTS_MAXIMUM_LIMBS_PER_INT)     if ((usize + vsize - 1) > GMP_INTS_MAXIMUM_LIMBS_PER_INT)
2188        {        {
2189        if (sign_product < 0)        if (sign_product < 0)
2190           w->flags = GMP_INTS_EF_INTOVF_NEG;           w->flags = GMP_INTS_EF_INTOVF_NEG;
2191        else        else
2192           w->flags = GMP_INTS_EF_INTOVF_POS;           w->flags = GMP_INTS_EF_INTOVF_POS;
2193    
2194        return;        return;
2195        }        }
2196    
2197    
2198     if (usize < vsize)     if (usize < vsize)
2199        {        {
2200        //Temporary variables just for the swap.        //Temporary variables just for the swap.
2201        const GMP_INTS_mpz_struct *tmp_ptr;        const GMP_INTS_mpz_struct *tmp_ptr;
2202        GMP_INTS_size_t            tmp_size;        GMP_INTS_size_t            tmp_size;
2203    
2204        //Swap U and V.        //Swap U and V.
2205        tmp_ptr = u;        tmp_ptr = u;
2206        u = v;        u = v;
2207        v = tmp_ptr;        v = tmp_ptr;
2208        tmp_size = usize;        tmp_size = usize;
2209        usize = vsize;        usize = vsize;
2210        vsize = tmp_size;        vsize = tmp_size;
2211        }        }
2212    
2213     //Grab pointers to the arrays of limbs.     //Grab pointers to the arrays of limbs.
2214     up = u->limbs;     up = u->limbs;
2215     vp = v->limbs;     vp = v->limbs;
2216     wp = w->limbs;     wp = w->limbs;
2217    
2218     /* Ensure W has space enough to store the result.  */     /* Ensure W has space enough to store the result.  */
2219     wsize = usize + vsize;     wsize = usize + vsize;
2220     if (w->n_allocd < wsize)     if (w->n_allocd < wsize)
2221        {        {
2222        if (wp == up || wp == vp)        if (wp == up || wp == vp)
2223           {           {
2224           free_me = wp;           free_me = wp;
2225           free_me_size = w->n_allocd;           free_me_size = w->n_allocd;
2226           }           }
2227        else        else
2228           {           {
2229           GMP_INTS_free_w_size (wp, w->n_allocd * sizeof(GMP_INTS_limb_t));           GMP_INTS_free_w_size (wp, w->n_allocd * sizeof(GMP_INTS_limb_t));
2230           }           }
2231    
2232        w->n_allocd = wsize;        w->n_allocd = wsize;
2233        wp = (GMP_INTS_limb_ptr)        wp = (GMP_INTS_limb_ptr)
2234             GMP_INTS_malloc (wsize * sizeof(GMP_INTS_limb_t));             GMP_INTS_malloc (wsize * sizeof(GMP_INTS_limb_t));
2235        w->limbs = wp;        w->limbs = wp;
2236        }        }
2237     else     else
2238        {        {
2239        /* Make U and V not overlap with W.  */        /* Make U and V not overlap with W.  */
2240        if (wp == up)        if (wp == up)
2241           {           {
2242           /* W and U are identical.  Allocate temporary space for U.  */           /* W and U are identical.  Allocate temporary space for U.  */
2243           up = (GMP_INTS_limb_ptr)           up = (GMP_INTS_limb_ptr)
2244                _alloca(usize * sizeof(GMP_INTS_limb_t));                _alloca(usize * sizeof(GMP_INTS_limb_t));
2245           /* Is V identical too?  Keep it identical with U.  */           /* Is V identical too?  Keep it identical with U.  */
2246           if (wp == vp)           if (wp == vp)
2247              vp = up;              vp = up;
2248           /* Copy to the temporary space.  */           /* Copy to the temporary space.  */
2249           GMP_INTS_mpn_copy_limbs(up, wp, usize);           GMP_INTS_mpn_copy_limbs(up, wp, usize);
2250           }           }
2251        else if (wp == vp)        else if (wp == vp)
2252           {           {
2253           /* W and V are identical.  Allocate temporary space for V.  */           /* W and V are identical.  Allocate temporary space for V.  */
2254           vp = (GMP_INTS_limb_ptr)           vp = (GMP_INTS_limb_ptr)
2255                _alloca(vsize * sizeof(GMP_INTS_limb_t));                _alloca(vsize * sizeof(GMP_INTS_limb_t));
2256           /* Copy to the temporary space.  */           /* Copy to the temporary space.  */
2257           GMP_INTS_mpn_copy_limbs(vp, wp, vsize);           GMP_INTS_mpn_copy_limbs(vp, wp, vsize);
2258           }           }
2259        }        }
2260    
2261     if (vsize == 0)     if (vsize == 0)
2262        {        {
2263        wsize = 0;        wsize = 0;
2264        }        }
2265     else     else
2266        {        {
2267        cy_limb = GMP_INTS_mpn_mul (wp, up, usize, vp, vsize);        cy_limb = GMP_INTS_mpn_mul (wp, up, usize, vp, vsize);
2268        wsize = usize + vsize;        wsize = usize + vsize;
2269        wsize -= cy_limb == 0;        wsize -= cy_limb == 0;
2270        }        }
2271    
2272     w->size = sign_product < 0 ? -wsize : wsize;     w->size = sign_product < 0 ? -wsize : wsize;
2273    
2274     if (free_me != NULL)     if (free_me != NULL)
2275        GMP_INTS_free_w_size (free_me, free_me_size * sizeof(GMP_INTS_limb_t));        GMP_INTS_free_w_size (free_me, free_me_size * sizeof(GMP_INTS_limb_t));
2276    
2277     //Final check for overflow.     //Final check for overflow.
2278     if (w->size > GMP_INTS_MAXIMUM_LIMBS_PER_INT)     if (w->size > GMP_INTS_MAXIMUM_LIMBS_PER_INT)
2279        w->flags = GMP_INTS_EF_INTOVF_POS;        w->flags = GMP_INTS_EF_INTOVF_POS;
2280     else if (w->size < -GMP_INTS_MAXIMUM_LIMBS_PER_INT)     else if (w->size < -GMP_INTS_MAXIMUM_LIMBS_PER_INT)
2281        w->flags = GMP_INTS_EF_INTOVF_NEG;        w->flags = GMP_INTS_EF_INTOVF_NEG;
2282     }     }
2283    
2284    
2285  //07/15/01:  Unit testing skipped because of recursive  //07/15/01:  Unit testing skipped because of recursive
2286  //nature.  Visual inspection OK.  //nature.  Visual inspection OK.
2287  void GMP_INTS_mpz_mul_si (      GMP_INTS_mpz_struct *w,  void GMP_INTS_mpz_mul_si (      GMP_INTS_mpz_struct *w,
2288                            const GMP_INTS_mpz_struct *u,                            const GMP_INTS_mpz_struct *u,
2289                                  long int             v)                                  long int             v)
2290     {     {
2291     GMP_INTS_mpz_struct temp;     GMP_INTS_mpz_struct temp;
2292    
2293     //Eyeball the inputs carefully.     //Eyeball the inputs carefully.
2294     assert(w != NULL);     assert(w != NULL);
2295     assert(w->n_allocd > 0);     assert(w->n_allocd > 0);
2296     assert(w->limbs != NULL);     assert(w->limbs != NULL);
2297     assert(u != NULL);     assert(u != NULL);
2298     assert(u->n_allocd > 0);     assert(u->n_allocd > 0);
2299     assert(u->limbs != NULL);     assert(u->limbs != NULL);
2300    
2301     //Create a temporary integer.     //Create a temporary integer.
2302     GMP_INTS_mpz_init(&temp);     GMP_INTS_mpz_init(&temp);
2303    
2304     //Set the temporary integer to the value of the input     //Set the temporary integer to the value of the input
2305     //argument.     //argument.
2306     GMP_INTS_mpz_set_si(&temp, v);     GMP_INTS_mpz_set_si(&temp, v);
2307    
2308     //Do the actual multiplication.  This recursive definition     //Do the actual multiplication.  This recursive definition
2309     //is inherently wasteful, but I'm after clarity, not     //is inherently wasteful, but I'm after clarity, not
2310     //warp speed.     //warp speed.
2311     GMP_INTS_mpz_mul(w, u, &temp);     GMP_INTS_mpz_mul(w, u, &temp);
2312    
2313     //Destroy the temporary integer (this will reclaim the     //Destroy the temporary integer (this will reclaim the
2314     //memory).     //memory).
2315     GMP_INTS_mpz_clear(&temp);     GMP_INTS_mpz_clear(&temp);
2316     }     }
2317    
2318    
2319  //07/15/01:  Unit testing skipped because of recursive  //07/15/01:  Unit testing skipped because of recursive
2320  //nature.  Visual inspection OK.  //nature.  Visual inspection OK.
2321  void GMP_INTS_mpz_mul_ui (      GMP_INTS_mpz_struct *w,  void GMP_INTS_mpz_mul_ui (      GMP_INTS_mpz_struct *w,
2322                            const GMP_INTS_mpz_struct *u,                            const GMP_INTS_mpz_struct *u,
2323                            unsigned long int          v)                            unsigned long int          v)
2324     {     {
2325     GMP_INTS_mpz_struct temp;     GMP_INTS_mpz_struct temp;
2326    
2327     //Eyeball the inputs carefully.     //Eyeball the inputs carefully.
2328     assert(w != NULL);     assert(w != NULL);
2329     assert(w->size >= 0);     assert(w->size >= 0);
2330     assert(w->limbs != NULL);     assert(w->limbs != NULL);
2331     assert(u != NULL);     assert(u != NULL);
2332     assert(u->size >= 0);     assert(u->size >= 0);
2333     assert(u->limbs != NULL);     assert(u->limbs != NULL);
2334    
2335     //Create a temporary integer.     //Create a temporary integer.
2336     GMP_INTS_mpz_init(&temp);     GMP_INTS_mpz_init(&temp);
2337    
2338     //Set the temporary integer to the value of the input     //Set the temporary integer to the value of the input
2339     //argument.     //argument.
2340     GMP_INTS_mpz_set_ui(&temp, v);     GMP_INTS_mpz_set_ui(&temp, v);
2341    
2342     //Do the actual multiplication.  This recursive definition     //Do the actual multiplication.  This recursive definition
2343     //is inherently wasteful, but I'm after clarity, not     //is inherently wasteful, but I'm after clarity, not
2344     //warp speed.     //warp speed.
2345     GMP_INTS_mpz_mul(w, u, &temp);     GMP_INTS_mpz_mul(w, u, &temp);
2346    
2347     //Destroy the temporary integer (this will reclaim the     //Destroy the temporary integer (this will reclaim the
2348     //memory).     //memory).
2349     GMP_INTS_mpz_clear(&temp);     GMP_INTS_mpz_clear(&temp);
2350     }     }
2351    
2352    
2353  void GMP_INTS_mpz_tdiv_qr (      GMP_INTS_mpz_struct *quot,  void GMP_INTS_mpz_tdiv_qr (      GMP_INTS_mpz_struct *quot,
2354                                   GMP_INTS_mpz_struct *rem,                                   GMP_INTS_mpz_struct *rem,
2355                             const GMP_INTS_mpz_struct *num,                             const GMP_INTS_mpz_struct *num,
2356                             const GMP_INTS_mpz_struct *den)                             const GMP_INTS_mpz_struct *den)
2357     {     {
2358     GMP_INTS_size_t   abs_num_size,     GMP_INTS_size_t   abs_num_size,
2359                       abs_den_size,                       abs_den_size,
2360                       quotient_sign,                       quotient_sign,
2361                       remainder_sign,                       remainder_sign,
2362                       numerator_bitsize,                       numerator_bitsize,
2363                       denominator_bitsize,                       denominator_bitsize,
2364                       division_loop_count,                       division_loop_count,
2365                       division_loop_count_mod_32,                       division_loop_count_mod_32,
2366                       division_loop_count_div_32,                       division_loop_count_div_32,
2367                       division_counter,                       division_counter,
2368                       i;                       i;
2369     GMP_INTS_limb_t   temp_limb;     GMP_INTS_limb_t   temp_limb;
2370     GMP_INTS_limb_ptr trial_divisor;     GMP_INTS_limb_ptr trial_divisor;
2371    
2372     //Eyeball the input parameters.     //Eyeball the input parameters.
2373     assert(quot != NULL);     assert(quot != NULL);
2374     assert(quot->n_allocd > 0);     assert(quot->n_allocd > 0);
2375     assert(quot->limbs != NULL);     assert(quot->limbs != NULL);
2376     assert(rem != NULL);     assert(rem != NULL);
2377     assert(rem->n_allocd > 0);     assert(rem->n_allocd > 0);
2378     assert(rem->limbs != NULL);     assert(rem->limbs != NULL);
2379     assert(num != NULL);     assert(num != NULL);
2380     assert(num->n_allocd > 0);     assert(num->n_allocd > 0);
2381     assert(num->limbs != NULL);     assert(num->limbs != NULL);
2382     assert(den != NULL);     assert(den != NULL);
2383     assert(den->n_allocd > 0);     assert(den->n_allocd > 0);
2384     assert(den->limbs != NULL);     assert(den->limbs != NULL);
2385    
2386     //We require for this function that the numerator, denominator, quotient, and     //We require for this function that the numerator, denominator, quotient, and
2387     //remainder all be distinct.     //remainder all be distinct.
2388     assert(quot != rem);     assert(quot != rem);
2389     assert(quot != num);     assert(quot != num);
2390     assert(quot != den);     assert(quot != den);
2391     assert(rem  != num);     assert(rem  != num);
2392     assert(rem  != den);     assert(rem  != den);
2393     assert(num  != den);     assert(num  != den);
2394    
2395     //The GNU code was probably very efficient, but exceeded     //The GNU code was probably very efficient, but exceeded
2396     //my abilities to analyze.  This is the classic     //my abilities to analyze.  This is the classic
2397     //division algorithm.     //division algorithm.
2398    
2399     //First, start off with the quotient and remainder having     //First, start off with the quotient and remainder having
2400     //no error flags set.  These will be set if appropriate.     //no error flags set.  These will be set if appropriate.
2401     quot->flags = 0;     quot->flags = 0;
2402     rem->flags  = 0;     rem->flags  = 0;
2403    
2404     //First, handle tainted inputs.  If the numerator or denominator     //First, handle tainted inputs.  If the numerator or denominator
2405     //are bad or tainted, the quotient and remainder get tainted     //are bad or tainted, the quotient and remainder get tainted
2406     //automatically.     //automatically.
2407        {        {
2408        int taint;        int taint;
2409    
2410        taint = GMP_INTS_two_op_flags_map(num->flags, den->flags);        taint = GMP_INTS_two_op_flags_map(num->flags, den->flags);
2411    
2412        if (taint)        if (taint)
2413           {           {
2414           quot->flags = taint;           quot->flags = taint;
2415           rem->flags  = taint;           rem->flags  = taint;
2416           return;           return;
2417           }           }
2418        }        }
2419    
2420     //The second possible cause for taint is if the divisor is     //The second possible cause for taint is if the divisor is
2421     //zero.  This will get both the value of positive overflow.     //zero.  This will get both the value of positive overflow.
2422     if (den->size == 0)     if (den->size == 0)
2423        {        {
2424        quot->flags = GMP_INTS_EF_INTOVF_POS;        quot->flags = GMP_INTS_EF_INTOVF_POS;
2425        rem->flags  = GMP_INTS_EF_INTOVF_POS;        rem->flags  = GMP_INTS_EF_INTOVF_POS;
2426        return;        return;
2427        }        }
2428    
2429     //Handle the special case of a numerator of zero.  If the numerator     //Handle the special case of a numerator of zero.  If the numerator
2430     //is zero, the quotient and remainder are zero automatically.     //is zero, the quotient and remainder are zero automatically.
2431     if (num->size == 0)     if (num->size == 0)
2432        {        {
2433        GMP_INTS_mpz_set_ui(quot, 0);        GMP_INTS_mpz_set_ui(quot, 0);
2434        GMP_INTS_mpz_set_ui(rem,  0);        GMP_INTS_mpz_set_ui(rem,  0);
2435        return;        return;
2436        }        }
2437    
2438     //Generally, nothing else can go wrong as far as taint.  The     //Generally, nothing else can go wrong as far as taint.  The
2439     //value of the quotient is confined to be no larger than the     //value of the quotient is confined to be no larger than the
2440     //numerator, and the value of the remainder is confined to     //numerator, and the value of the remainder is confined to
2441     //be no larger than denominator-1.  So, generally, if the     //be no larger than denominator-1.  So, generally, if the
2442     //inputs are in size bounds, the outputs will be also.     //inputs are in size bounds, the outputs will be also.
2443    
2444     //Figure out how large in limbs the numerator and denominator actually     //Figure out how large in limbs the numerator and denominator actually
2445     //are.     //are.
2446     abs_num_size = GMP_INTS_abs_of_size_t(num->size);     abs_num_size = GMP_INTS_abs_of_size_t(num->size);
2447     abs_den_size = GMP_INTS_abs_of_size_t(den->size);     abs_den_size = GMP_INTS_abs_of_size_t(den->size);
2448    
2449     //Figure out the sign of things.  We want the following relationship     //Figure out the sign of things.  We want the following relationship
2450     //to be true:     //to be true:
2451     //   num/den = quot + rem/den.     //   num/den = quot + rem/den.
2452     //The way to achieve this is to assign the sign of the quotient in the traditional     //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.     //way, then to assign the remainder to have the same sign as the numerator.
2454     quotient_sign  = num->size ^ den->size;     quotient_sign  = num->size ^ den->size;
2455     remainder_sign = num->size;     remainder_sign = num->size;
2456    
2457     //The remainder starts off with the absolute value of the numerator, and then     //The remainder starts off with the absolute value of the numerator, and then
2458     //we subtract from it as part of the division loop.     //we subtract from it as part of the division loop.
2459     GMP_INTS_mpz_copy(rem, num);     GMP_INTS_mpz_copy(rem, num);
2460        //We know after the copy that the amount of space allocated in the remainder        //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        //MUST be at least as large as the absolute value of the numerator.  So from
2462        //this point forward we use the space.        //this point forward we use the space.
2463     assert(rem->n_allocd >= abs_num_size);     assert(rem->n_allocd >= abs_num_size);
2464    
2465     //Figure out the number of significant bits in the numerator and denominator.     //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.     //This determines the loop count over which we do the shift division loop.
2467     numerator_bitsize = GMP_INTS_BITS_PER_LIMB * abs_num_size;     numerator_bitsize = GMP_INTS_BITS_PER_LIMB * abs_num_size;
2468        
2469     i = abs_num_size - 1;     i = abs_num_size - 1;
2470    
2471     //We need to be extra careful here.  One failure mode is that an integer     //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     //data structure is corrupted and the "size" field reflects limbs
2473     //that are zero.  Need to watch that this kind of failure doesn't     //that are zero.  Need to watch that this kind of failure doesn't
2474     //cause memory access errors.     //cause memory access errors.
2475     assert(num->limbs[i] != 0);     assert(num->limbs[i] != 0);
2476     if (num->limbs[i] == 0)     if (num->limbs[i] == 0)
2477        {        {
2478        quot->flags = GMP_INTS_EF_INTOVF_POS;        quot->flags = GMP_INTS_EF_INTOVF_POS;
2479        rem->flags  = GMP_INTS_EF_INTOVF_POS;        rem->flags  = GMP_INTS_EF_INTOVF_POS;
2480        return;        return;
2481        }        }
2482    
2483     temp_limb = 0x80000000;     temp_limb = 0x80000000;
2484    
2485     while (((num->limbs[i] & temp_limb) == 0) && (temp_limb != 0))     while (((num->limbs[i] & temp_limb) == 0) && (temp_limb != 0))
2486        {        {
2487        numerator_bitsize--;        numerator_bitsize--;
2488        temp_limb >>= 1;        temp_limb >>= 1;
2489        }        }
2490    
2491     denominator_bitsize = GMP_INTS_BITS_PER_LIMB * abs_den_size;     denominator_bitsize = GMP_INTS_BITS_PER_LIMB * abs_den_size;
2492    
2493     i = abs_den_size - 1;     i = abs_den_size - 1;
2494    
2495     //We need to be extra careful here.  One failure mode is that an integer     //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     //data structure is corrupted and the "size" field reflects limbs
2497     //that are zero.  Need to watch that this kind of failure doesn't     //that are zero.  Need to watch that this kind of failure doesn't
2498     //cause memory access errors.     //cause memory access errors.
2499     assert(den->limbs[i] != 0);     assert(den->limbs[i] != 0);
2500     if (den->limbs[i] == 0)     if (den->limbs[i] == 0)
2501        {        {
2502        quot->flags = GMP_INTS_EF_INTOVF_POS;        quot->flags = GMP_INTS_EF_INTOVF_POS;
2503        rem->flags  = GMP_INTS_EF_INTOVF_POS;        rem->flags  = GMP_INTS_EF_INTOVF_POS;
2504        return;        return;
2505        }        }
2506    
2507     temp_limb = 0x80000000;     temp_limb = 0x80000000;
2508    
2509     while (((den->limbs[i] & temp_limb) == 0) && (temp_limb != 0))     while (((den->limbs[i] & temp_limb) == 0) && (temp_limb != 0))
2510        {        {
2511        denominator_bitsize--;        denominator_bitsize--;
2512        temp_limb >>= 1;        temp_limb >>= 1;
2513        }        }
2514    
2515     //The quotient starts off with the value of zero, but we consistently may     //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     //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     //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.     //prepare a block of clear memory as large as the numerator's.
2519     if (quot->n_allocd < abs_num_size)     if (quot->n_allocd < abs_num_size)
2520        {        {
2521        GMP_INTS_mpz_realloc(quot, abs_num_size);  //Make it big enough.        GMP_INTS_mpz_realloc(quot, abs_num_size);  //Make it big enough.
2522        }        }
2523     //Now, zero the memory.     //Now, zero the memory.
2524     for (i=0; i<abs_num_size; i++)     for (i=0; i<abs_num_size; i++)
2525        quot->limbs[i] = 0;        quot->limbs[i] = 0;
2526    
2527     //Determine the division loop count.  This is the difference     //Determine the division loop count.  This is the difference
2528     //in bit sizes between the numerator and denominator.  It is     //in bit sizes between the numerator and denominator.  It is
2529     //possible for this number to be negative, which means that the     //possible for this number to be negative, which means that the
2530     //main division loop will be executed zero times.  This gives the     //main division loop will be executed zero times.  This gives the
2531     //right results.     //right results.
2532     division_loop_count = numerator_bitsize - denominator_bitsize;     division_loop_count = numerator_bitsize - denominator_bitsize;
2533    
2534     //We need to calculate some important numbers from the division loop     //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     //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     //shift the divisor bitwise to line up with the numerator), and we
2537     //also need this number DIV 32 for the limb-wise shift.     //also need this number DIV 32 for the limb-wise shift.
2538     division_loop_count_mod_32 = division_loop_count % 32;     division_loop_count_mod_32 = division_loop_count % 32;
2539     division_loop_count_div_32 = division_loop_count / 32;     division_loop_count_div_32 = division_loop_count / 32;
2540        
2541     //We now need a shift register in which we shift the denominator up     //We now need a shift register in which we shift the denominator up
2542     //for repeated comparisons.  We should dynamically allocate this to     //for repeated comparisons.  We should dynamically allocate this to
2543     //be the same size as the numerator.  Using _alloca() is OK, as one     //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     //of the unit tests is to be sure that _alloca() will handle integer
2545     //of the maximum allowed size.     //of the maximum allowed size.
2546     trial_divisor = _alloca(abs_num_size * sizeof(GMP_INTS_limb_t));     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     //Our trial divisor needs to start off with the divisor shifted up
2549     //so that the most significant bit is aligned with the numerator.     //so that the most significant bit is aligned with the numerator.
2550     for (i = 0; i < abs_num_size; i++)     for (i = 0; i < abs_num_size; i++)
2551        {        {
2552        if ((division_loop_count < 0) || (i < division_loop_count_div_32))        if ((division_loop_count < 0) || (i < division_loop_count_div_32))
2553           {           {
2554           trial_divisor[i] = 0;           trial_divisor[i] = 0;
2555           }           }
2556        else        else
2557           {           {
2558           if ((i-division_loop_count_div_32) < abs_den_size)           if ((i-division_loop_count_div_32) < abs_den_size)
2559              trial_divisor[i] = den->limbs[i - division_loop_count_div_32];              trial_divisor[i] = den->limbs[i - division_loop_count_div_32];
2560           else           else
2561              trial_divisor[i] = 0;              trial_divisor[i] = 0;
2562           }           }
2563        }        }
2564    
2565     //The code above planted the limbs in the right place.  Now need to shift bits     //The code above planted the limbs in the right place.  Now need to shift bits
2566     //upward by the remaining number.     //upward by the remaining number.
2567     if ((division_loop_count > 0) && (division_loop_count_mod_32 > 0))     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.        //There is an existing function we can call to do the left shift.
2570        GMP_INTS_mpn_lshift(trial_divisor,        GMP_INTS_mpn_lshift(trial_divisor,
2571                            trial_divisor,                            trial_divisor,
2572                            abs_num_size,                            abs_num_size,
2573                            division_loop_count_mod_32);                            division_loop_count_mod_32);
2574        }        }
2575    
2576    
2577     //Everything is ready to go.  Now begin the division loop itself.  It is possible     //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     //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     //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.     //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++)     for (division_counter = 0; division_counter < division_loop_count+1; division_counter++)
2582        {        {
2583        //Shift the quotient left one bit.        //Shift the quotient left one bit.
2584        GMP_INTS_mpn_lshift(quot->limbs,        GMP_INTS_mpn_lshift(quot->limbs,
2585                            quot->limbs,                            quot->limbs,
2586                            abs_num_size,                            abs_num_size,
2587                            1);                            1);
2588    
2589        //If the remainder is at least as large as the trial divisor, subtract the trial        //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.        //divisor from the remainder and mask in the quotient.
2591        if (GMP_INTS_mpn_cmp(rem->limbs,        if (GMP_INTS_mpn_cmp(rem->limbs,
2592                             trial_divisor,                             trial_divisor,
2593                             abs_num_size) >= 0)                             abs_num_size) >= 0)
2594           {           {
2595           GMP_INTS_mpn_sub(rem->limbs,           GMP_INTS_mpn_sub(rem->limbs,
2596                            rem->limbs,                            rem->limbs,
2597                            abs_num_size,                            abs_num_size,
2598                            trial_divisor,                            trial_divisor,
2599                            abs_num_size);                            abs_num_size);
2600           quot->limbs[0] |= 1;           quot->limbs[0] |= 1;
2601           }           }
2602    
2603        //Shift the trial divisor right one bit.        //Shift the trial divisor right one bit.
2604        GMP_INTS_mpn_rshift(trial_divisor,        GMP_INTS_mpn_rshift(trial_divisor,
2605                            trial_divisor,                            trial_divisor,
2606                            abs_num_size,                            abs_num_size,
2607                            1);                            1);
2608        }  //End for each iteration of the division loop.        }  //End for each iteration of the division loop.
2609    
2610     //Normalize the quotient and the remainder.  The normalization     //Normalize the quotient and the remainder.  The normalization
2611     //process is to bring the sizes down if we have leading     //process is to bring the sizes down if we have leading
2612     //zeros.     //zeros.
2613     quot->size = abs_num_size;     quot->size = abs_num_size;
2614     GMP_INTS_mpn_normalize(quot->limbs, &(quot->size));     GMP_INTS_mpn_normalize(quot->limbs, &(quot->size));
2615     rem->size  = abs_num_size;     rem->size  = abs_num_size;
2616     GMP_INTS_mpn_normalize(rem->limbs, &(rem->size));     GMP_INTS_mpn_normalize(rem->limbs, &(rem->size));
2617    
2618     //Adjust the signs as required.     //Adjust the signs as required.
2619     if (quotient_sign < 0)     if (quotient_sign < 0)
2620        quot->size = -(quot->size);        quot->size = -(quot->size);
2621     if (remainder_sign < 0)     if (remainder_sign < 0)
2622        rem->size = -(rem->size);        rem->size = -(rem->size);
2623     }     }
2624    
2625    
2626  void GMP_INTS_mpz_fac_ui(GMP_INTS_mpz_struct *result,  void GMP_INTS_mpz_fac_ui(GMP_INTS_mpz_struct *result,
2627                           unsigned long int n)                           unsigned long int n)
2628     {     {
2629     //Just multiply the numbers in ascending order.  The original     //Just multiply the numbers in ascending order.  The original
2630     //GNU library contained a much more elegant algorithm, but     //GNU library contained a much more elegant algorithm, but
2631     //this is more direct.     //this is more direct.
2632    
2633     unsigned long int k;     unsigned long int k;
2634    
2635     GMP_INTS_mpz_set_ui (result, 1L);     GMP_INTS_mpz_set_ui (result, 1L);
2636    
2637     for (k = 2; (k <= n) && !(result->flags); k++)     for (k = 2; (k <= n) && !(result->flags); k++)
2638        GMP_INTS_mpz_mul_ui (result, result, k);        GMP_INTS_mpz_mul_ui (result, result, k);
2639     }     }
2640    
2641    
2642  /******************************************************************/  /******************************************************************/
2643  /***  PUBLIC CONVERSION AND OUTPUT FUNCTIONS   ********************/  /***  PUBLIC CONVERSION AND OUTPUT FUNCTIONS   ********************/
2644  /******************************************************************/  /******************************************************************/
2645  //07/18/01:  Visual inspection OK.  Function returns  //07/18/01:  Visual inspection OK.  Function returns
2646  //reasonable values even out to 100,000 digits--seems OK.  //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)  int GMP_INTS_mpz_size_in_base_10(const GMP_INTS_mpz_struct *arg)
2648     {     {
2649     _int64 n;     _int64 n;
2650    
2651     //Eyeball the input parameter.     //Eyeball the input parameter.
2652     assert(arg != NULL);     assert(arg != NULL);
2653     assert(arg->n_allocd > 0);     assert(arg->n_allocd > 0);
2654     assert(arg->limbs != NULL);     assert(arg->limbs != NULL);
2655    
2656     //Get the number of limbs occupied by the integer.     //Get the number of limbs occupied by the integer.
2657     //Because even the digit zero takes some space,     //Because even the digit zero takes some space,
2658     //don't accept zero for an answer.     //don't accept zero for an answer.
2659     n = GMP_INTS_abs_of_size_t(arg->size);     n = GMP_INTS_abs_of_size_t(arg->size);
2660     if (n==0)     if (n==0)
2661        n = 1;        n = 1;
2662        
2663     //Convert this to the number of bits.  Generously     //Convert this to the number of bits.  Generously
2664     //ignore any unused leading bits.     //ignore any unused leading bits.
2665     n *= 32;     n *= 32;
2666    
2667     //Used a slightly high best rational approximation in F_{65535}     //Used a slightly high best rational approximation in F_{65535}
2668     //to go from the number of bits to the number of     //to go from the number of bits to the number of
2669     //digits.  The division discards, so bump the result     //digits.  The division discards, so bump the result
2670     //up by 1 to compensate for possible truncation.  The number     //up by 1 to compensate for possible truncation.  The number
2671     //we are aproximating is ln(2)/ln(10).     //we are aproximating is ln(2)/ln(10).
2672     n *= 12655;     n *= 12655;
2673     n /= 42039;     n /= 42039;
2674     n++;     n++;
2675    
2676     //Compensate for possible commas in the result.  Again,     //Compensate for possible commas in the result.  Again,
2677     //consider truncation.     //consider truncation.
2678     n *= 4;     n *= 4;
2679     n /= 3;     n /= 3;
2680     n++;     n++;
2681    
2682     //Compensate for the minus sign, the trailing zero,     //Compensate for the minus sign, the trailing zero,
2683     //cosmic rays striking the computer from the martian     //cosmic rays striking the computer from the martian
2684     //listening post camoflaged on the moon, and the     //listening post camoflaged on the moon, and the
2685     //possibility that we might need to put text in the     //possibility that we might need to put text in the
2686     //string if any flag is set.     //string if any flag is set.
2687     n += 100;     n += 100;
2688    
2689     //And that should be a good return value.     //And that should be a good return value.
2690     return((int) n);     return((int) n);
2691     }     }
2692    
2693    
2694  //07/19/01:  Visual inspection and unit test is OK.  //07/19/01:  Visual inspection and unit test is OK.
2695  void GMP_INTS_mpz_to_string(char *out,  void GMP_INTS_mpz_to_string(char *out,
2696                              const GMP_INTS_mpz_struct *in)                              const GMP_INTS_mpz_struct *in)
2697     {     {
2698     //Eyeball the input parameters.     //Eyeball the input parameters.
2699     assert(out != NULL);     assert(out != NULL);
2700     assert(in != NULL);     assert(in != NULL);
2701     assert(in->n_allocd > 0);     assert(in->n_allocd > 0);
2702     assert(in->limbs != NULL);     assert(in->limbs != NULL);
2703    
2704     //If any of the flags are set, stuff in the text.     //If any of the flags are set, stuff in the text.
2705     if (in->flags)     if (in->flags)
2706        {        {
2707        if (in->flags & GMP_INTS_EF_INTOVF_POS)        if (in->flags & GMP_INTS_EF_INTOVF_POS)
2708           {           {
2709           strcpy(out, GMP_INTS_EF_INTOVF_POS_STRING);           strcpy(out, GMP_INTS_EF_INTOVF_POS_STRING);
2710           }           }
2711        else if (in->flags & GMP_INTS_EF_INTOVF_NEG)        else if (in->flags & GMP_INTS_EF_INTOVF_NEG)
2712           {           {
2713           strcpy(out, GMP_INTS_EF_INTOVF_NEG_STRING);           strcpy(out, GMP_INTS_EF_INTOVF_NEG_STRING);
2714           }           }
2715        else if (in->flags & GMP_INTS_EF_INTOVF_TAINT_POS)        else if (in->flags & GMP_INTS_EF_INTOVF_TAINT_POS)
2716           {           {
2717           strcpy(out, GMP_INTS_EF_INTOVF_TAINT_POS_STRING);           strcpy(out, GMP_INTS_EF_INTOVF_TAINT_POS_STRING);
2718           }           }
2719        else if (in->flags & GMP_INTS_EF_INTOVF_TAINT_NEG)        else if (in->flags & GMP_INTS_EF_INTOVF_TAINT_NEG)
2720           {           {
2721           strcpy(out, GMP_INTS_EF_INTOVF_TAINT_NEG_STRING);           strcpy(out, GMP_INTS_EF_INTOVF_TAINT_NEG_STRING);
2722           }           }
2723        else        else
2724           {           {
2725           strcpy(out, "INTERNAL_ERROR");           strcpy(out, "INTERNAL_ERROR");
2726           }           }
2727        }          }  
2728     else     else
2729        {        {
2730        //Ordinary integer conversion.        //Ordinary integer conversion.
2731        GMP_INTS_mpz_struct num, den, quot, rem, k10;        GMP_INTS_mpz_struct num, den, quot, rem, k10;
2732    
2733        //Allocate space for the temporary integers.        //Allocate space for the temporary integers.
2734        GMP_INTS_mpz_init(&num);        GMP_INTS_mpz_init(&num);
2735        GMP_INTS_mpz_init(&den);        GMP_INTS_mpz_init(&den);
2736        GMP_INTS_mpz_init(&quot);        GMP_INTS_mpz_init(&quot);
2737        GMP_INTS_mpz_init(&rem);        GMP_INTS_mpz_init(&rem);
2738        GMP_INTS_mpz_init(&k10);        GMP_INTS_mpz_init(&k10);
2739    
2740        //Assign the constant 10.        //Assign the constant 10.
2741        GMP_INTS_mpz_set_ui(&k10, 10);        GMP_INTS_mpz_set_ui(&k10, 10);
2742    
2743        //If the integer is zero, assign that.        //If the integer is zero, assign that.
2744        if (in->size == 0)        if (in->size == 0)
2745           {           {
2746           strcpy(out, "0");           strcpy(out, "0");
2747           }           }
2748        else        else
2749           {           {
2750           //We have to do a full conversion.  The algorithm           //We have to do a full conversion.  The algorithm
2751           //is division by 10, each time obtaining the least           //is division by 10, each time obtaining the least
2752           //significant digit, until finally the quotient is           //significant digit, until finally the quotient is
2753           //zero.           //zero.
2754           char *ptr;           char *ptr;
2755    
2756           ptr = out;           ptr = out;
2757    
2758           GMP_INTS_mpz_copy(&num, in);           GMP_INTS_mpz_copy(&num, in);
2759           GMP_INTS_mpz_copy(&den, &k10);           GMP_INTS_mpz_copy(&den, &k10);
2760           do           do
2761              {              {
2762  #if 0  #if 0
2763              printf("Values before division:\n");              printf("Values before division:\n");
2764              FCMIOF_hline();              FCMIOF_hline();
2765              GMP_INTS_mpz_print_int(stdout, &num, "Numerator");              GMP_INTS_mpz_print_int(stdout, &num, "Numerator");
2766              FCMIOF_hline();              FCMIOF_hline();
2767              GMP_INTS_mpz_print_int(stdout, &den, "Denominator");              GMP_INTS_mpz_print_int(stdout, &den, "Denominator");
2768              FCMIOF_hline();              FCMIOF_hline();
2769              GMP_INTS_mpz_print_int(stdout, &quot, "Quotient");              GMP_INTS_mpz_print_int(stdout, &quot, "Quotient");
2770              FCMIOF_hline();              FCMIOF_hline();
2771              GMP_INTS_mpz_print_int(stdout, &rem, "Remainder");              GMP_INTS_mpz_print_int(stdout, &rem, "Remainder");
2772              FCMIOF_hline();              FCMIOF_hline();
2773    
2774  if (num.size > 1)  if (num.size > 1)
2775              FCMIOF_hline();              FCMIOF_hline();
2776  #endif  #endif
2777    
2778              GMP_INTS_mpz_tdiv_qr(&quot, &rem, &num, &den);              GMP_INTS_mpz_tdiv_qr(&quot, &rem, &num, &den);
2779  #if 0  #if 0
2780              printf("Values after division:\n");              printf("Values after division:\n");
2781              FCMIOF_hline();              FCMIOF_hline();
2782              GMP_INTS_mpz_print_int(stdout, &num, "Numerator");              GMP_INTS_mpz_print_int(stdout, &num, "Numerator");
2783              FCMIOF_hline();              FCMIOF_hline();
2784              GMP_INTS_mpz_print_int(stdout, &den, "Denominator");              GMP_INTS_mpz_print_int(stdout, &den, "Denominator");
2785              FCMIOF_hline();              FCMIOF_hline();
2786              GMP_INTS_mpz_print_int(stdout, &quot, "Quotient");              GMP_INTS_mpz_print_int(stdout, &quot, "Quotient");
2787              FCMIOF_hline();              FCMIOF_hline();
2788              GMP_INTS_mpz_print_int(stdout, &rem, "Remainder");              GMP_INTS_mpz_print_int(stdout, &rem, "Remainder");
2789              FCMIOF_hline();              FCMIOF_hline();
2790  #endif  #endif
2791    
2792              if (rem.size != 0)              if (rem.size != 0)
2793                 {                 {
2794                 *ptr = '0' + (char)(rem.limbs[0]);                 *ptr = '0' + (char)(rem.limbs[0]);
2795                 }                 }
2796              else              else
2797                 {                 {
2798                 *ptr = '0';                 *ptr = '0';
2799                 }                 }
2800              ptr++;              ptr++;
2801              GMP_INTS_mpz_copy(&num, &quot);              GMP_INTS_mpz_copy(&num, &quot);