/[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 70 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);
2802              //printf("digit\n");              //printf("digit\n");
2803              }              }
2804           while (!GMP_INTS_mpz_is_zero(&quot));           while (!GMP_INTS_mpz_is_zero(&quot));
2805    
2806           //Finally, if the input was negative, tack on the           //Finally, if the input was negative, tack on the
2807           //minus sign.           //minus sign.
2808           if (GMP_INTS_mpz_is_neg(in))           if (GMP_INTS_mpz_is_neg(in))
2809              {              {
2810              *ptr = '-';              *ptr = '-';
2811              ptr++;              ptr++;
2812              }              }
2813                    
2814           //Finally, tack on the trailing zero terminator.           //Finally, tack on the trailing zero terminator.
2815           *ptr = 0;           *ptr = 0;
2816           ptr++;           ptr++;
2817    
2818           //Reverse the string.           //Reverse the string.
2819           BSTRFUNC_str_reverse(out);           BSTRFUNC_str_reverse(out);
2820           }           }
2821    
2822        //Deallocate the integers.        //Deallocate the integers.
2823        GMP_INTS_mpz_clear(&num);        GMP_INTS_mpz_clear(&num);
2824        GMP_INTS_mpz_clear(&den);        GMP_INTS_mpz_clear(&den);
2825        GMP_INTS_mpz_clear(&quot);        GMP_INTS_mpz_clear(&quot);
2826        GMP_INTS_mpz_clear(&rem);        GMP_INTS_mpz_clear(&rem);
2827        GMP_INTS_mpz_clear(&k10);        GMP_INTS_mpz_clear(&k10);
2828        }        }
2829     }     }
2830    
2831    
2832  void GMP_INTS_mpz_long_int_format_to_stream(FILE *s,  void GMP_INTS_mpz_long_int_format_to_stream(FILE *s,
2833                                              const GMP_INTS_mpz_struct *i,                                              const GMP_INTS_mpz_struct *i,
2834                                              const char *desc)                                              const char *desc)
2835     {     {
2836     int line_len;     int line_len;
2837     int digits_per_line;     int digits_per_line;
2838     char *digits;     char *digits;
2839     int num_digits;     int num_digits;
2840     int nlines;     int nlines;
2841     int cur_line;     int cur_line;
2842     int number_desc_width;     int number_desc_width;
2843    
2844     //Eyeball the inputs, make sure the caller isn't doing     //Eyeball the inputs, make sure the caller isn't doing
2845     //something stupid.     //something stupid.
2846     assert(s != NULL);     assert(s != NULL);
2847     assert(i != NULL);     assert(i != NULL);
2848     assert(i->n_allocd > 0);     assert(i->n_allocd > 0);
2849     assert(i->limbs != NULL);     assert(i->limbs != NULL);
2850     assert(desc != NULL);     assert(desc != NULL);
2851    
2852     //Obtain the line length assumed for formatted output.     //Obtain the line length assumed for formatted output.
2853     line_len = FCMIOF_get_line_len();     line_len = FCMIOF_get_line_len();
2854    
2855     //The description width allowed is 20.     //The description width allowed is 20.
2856     number_desc_width = 20;     number_desc_width = 20;
2857    
2858     /* The number of digits per line that we assume must be a multiple of     /* The number of digits per line that we assume must be a multiple of
2859     ** three.  The formula below was not examined very carefully, but it     ** three.  The formula below was not examined very carefully, but it
2860     ** works fine for a line length of 78.  If line length is changed,     ** works fine for a line length of 78.  If line length is changed,
2861     ** this formula may need to be examined very carefully and rewritten.     ** this formula may need to be examined very carefully and rewritten.
2862     */     */
2863     digits_per_line = INTFUNC_max(3, ((((line_len-42)*3)/4)/3)*3);     digits_per_line = INTFUNC_max(3, ((((line_len-42)*3)/4)/3)*3);
2864     assert(digits_per_line >= 3);     assert(digits_per_line >= 3);
2865    
2866     /* We now need to get a digit string corresponding to this     /* We now need to get a digit string corresponding to this
2867     ** number.  First, need to figure out how much and     ** number.  First, need to figure out how much and
2868     ** allocate the space.     ** allocate the space.
2869     */     */
2870     digits = GMP_INTS_malloc(GMP_INTS_mpz_size_in_base_10(i) * sizeof(char));     digits = GMP_INTS_malloc(GMP_INTS_mpz_size_in_base_10(i) * sizeof(char));
2871     GMP_INTS_mpz_to_string(digits, i);     GMP_INTS_mpz_to_string(digits, i);
2872    
2873     //If the number is negative, delete the leading minus sign.     //If the number is negative, delete the leading minus sign.
2874     //The rest of the display algorithm needs an unsigned     //The rest of the display algorithm needs an unsigned
2875     //series of digits.     //series of digits.
2876     if (*digits == '-')     if (*digits == '-')
2877        {        {
2878        int i = 0;        int i = 0;
2879    
2880        do        do
2881           {           {
2882           digits[i] = digits[i+1];           digits[i] = digits[i+1];
2883           i++;           i++;
2884           }           }
2885        while(digits[i-1]);        while(digits[i-1]);
2886        }        }
2887    
2888     //Figure out how many digits in the string representation.     //Figure out how many digits in the string representation.
2889     num_digits = strlen(digits);     num_digits = strlen(digits);
2890    
2891     /* As the first order of business, figure out how many lines the beast     /* As the first order of business, figure out how many lines the beast
2892     ** will require.     ** will require.
2893     */     */
2894     if (i->flags)     if (i->flags)
2895        {        {
2896        nlines = 1;  /* Only one line required for NAN verbeage. */        nlines = 1;  /* Only one line required for NAN verbeage. */
2897        }        }
2898     else if (GMP_INTS_mpz_is_zero(i))     else if (GMP_INTS_mpz_is_zero(i))
2899        {        {
2900        nlines = 1;  /* The zero value requires one line. */        nlines = 1;  /* The zero value requires one line. */
2901        }        }
2902     else     else
2903        {        {
2904        /* In any other case, have a formula.        /* In any other case, have a formula.
2905        */        */
2906        nlines = 1 + (num_digits - 1) / digits_per_line;        nlines = 1 + (num_digits - 1) / digits_per_line;
2907        }        }
2908    
2909     /* Iterate through each line, spitting out whatever is appropriate. */     /* Iterate through each line, spitting out whatever is appropriate. */
2910     for (cur_line = 0; cur_line < nlines; cur_line++)     for (cur_line = 0; cur_line < nlines; cur_line++)
2911        {        {
2912        int cur_digit_on_line;        int cur_digit_on_line;
2913    
2914        /* If this is the first line, spit out the description, right-aligned.        /* If this is the first line, spit out the description, right-aligned.
2915        ** Otherwise, spit spaces.        ** Otherwise, spit spaces.
2916        */        */
2917        if (!cur_line)        if (!cur_line)
2918           {           {
2919           /* First line. */           /* First line. */
2920           int len;           int len;
2921    
2922           len = strlen(desc);           len = strlen(desc);
2923    
2924           if (len <= number_desc_width)           if (len <= number_desc_width)
2925              {              {
2926              /* Description is shorter or equal, pad on left. */              /* Description is shorter or equal, pad on left. */
2927              FCMIOF_stream_repchar(s, ' ', number_desc_width - len);              FCMIOF_stream_repchar(s, ' ', number_desc_width - len);
2928              fprintf(s, "%s", desc);              fprintf(s, "%s", desc);
2929              }              }
2930           else           else
2931              {              {
2932              /* Description is too long, truncate. */              /* Description is too long, truncate. */
2933              int i;              int i;
2934    
2935              for (i=0; i<number_desc_width; i++)              for (i=0; i<number_desc_width; i++)
2936                 fprintf(s, "%c", desc[i]);                 fprintf(s, "%c", desc[i]);
2937              }              }
2938    
2939           fprintf(s, ": ");           fprintf(s, ": ");
2940    
2941           /* If the number is negative, throw in a minus sign. */           /* If the number is negative, throw in a minus sign. */
2942           if (GMP_INTS_mpz_is_neg(i) && !(i->flags))           if (GMP_INTS_mpz_is_neg(i) && !(i->flags))
2943              {              {
2944              fprintf(s, "- ");              fprintf(s, "- ");
2945              }              }
2946           else           else
2947              {              {
2948              fprintf(s, "  ");              fprintf(s, "  ");
2949              }              }
2950           }           }
2951        else        else
2952           {           {
2953           /* Every line but first line. */           /* Every line but first line. */
2954           FCMIOF_stream_repchar(s, ' ', number_desc_width+4);           FCMIOF_stream_repchar(s, ' ', number_desc_width+4);
2955           }           }
2956    
2957        for(cur_digit_on_line=0; cur_digit_on_line < digits_per_line; cur_digit_on_line++)        for(cur_digit_on_line=0; cur_digit_on_line < digits_per_line; cur_digit_on_line++)
2958           {           {
2959           int idx_into_string;           int idx_into_string;
2960              /* Index into the string which is our digit of interest.              /* Index into the string which is our digit of interest.
2961              */              */
2962    
2963           /* Compute the index.  The equation is based on the ordering           /* Compute the index.  The equation is based on the ordering
2964           ** of presentation, for example,           ** of presentation, for example,
2965           **           **
2966           **            7 6           **            7 6
2967           **          5 4 3           **          5 4 3
2968           **          2 1 0.           **          2 1 0.
2969           **           **
2970           ** With a little thought, the equation should make sense.           ** With a little thought, the equation should make sense.
2971           ** The index won't always be used to index into the string.           ** The index won't always be used to index into the string.
2972           */           */
2973           idx_into_string =           idx_into_string =
2974             ((((nlines-1) - cur_line) * digits_per_line)             ((((nlines-1) - cur_line) * digits_per_line)
2975             +             +
2976             (digits_per_line - 1 - cur_digit_on_line));             (digits_per_line - 1 - cur_digit_on_line));
2977    
2978           /* Print the appropriate digit or a space.  The NAN case and the           /* Print the appropriate digit or a space.  The NAN case and the
2979           ** zero case need to be treated specially.           ** zero case need to be treated specially.
2980           */           */
2981           if (i->flags)           if (i->flags)
2982              {              {
2983              /* Not a number.  Everything is blank, except spell out              /* Not a number.  Everything is blank, except spell out
2984              ** description of condition at the end of the string of              ** description of condition at the end of the string of
2985              ** digits.              ** digits.
2986              */              */
2987              int index_from_right;              int index_from_right;
2988              int virtual_index;              int virtual_index;
2989    
2990              index_from_right = digits_per_line - 1 - cur_digit_on_line;              index_from_right = digits_per_line - 1 - cur_digit_on_line;
2991                 //The index calculated above is calculated so that the                 //The index calculated above is calculated so that the
2992                 //final position on the line has index [0].                 //final position on the line has index [0].
2993              assert(index_from_right >= 0 && index_from_right < digits_per_line);              assert(index_from_right >= 0 && index_from_right < digits_per_line);
2994    
2995              //Now, calculate the "virtual index".  The virtual index              //Now, calculate the "virtual index".  The virtual index
2996              //is the actual number of characters from the right, taking              //is the actual number of characters from the right, taking
2997              //into account commas.              //into account commas.
2998              virtual_index = index_from_right + index_from_right/3;              virtual_index = index_from_right + index_from_right/3;
2999    
3000              if (((index_from_right % 3) == 2) && cur_digit_on_line)              if (((index_from_right % 3) == 2) && cur_digit_on_line)
3001                 {                 {
3002                 //We are one position past a comma.  This means                 //We are one position past a comma.  This means
3003                 //that we might need a "fill" character to go                 //that we might need a "fill" character to go
3004                 //where the comma should have gone.                 //where the comma should have gone.
3005    
3006                 if (virtual_index + 1 < num_digits)                 if (virtual_index + 1 < num_digits)
3007                    {                    {
3008                    //The character we should print exists.                    //The character we should print exists.
3009                    fprintf(s, "%c", digits[num_digits - 2 - virtual_index]);                    fprintf(s, "%c", digits[num_digits - 2 - virtual_index]);
3010                    }                    }
3011                 else                 else
3012                    {                    {
3013                    //The character doesn't exist, because the error                    //The character doesn't exist, because the error
3014                    //string is apparently too short.  Must print a                    //string is apparently too short.  Must print a
3015                    //space, instead.                    //space, instead.
3016                    fprintf(s, " ");                    fprintf(s, " ");
3017                    }                    }
3018                 }                 }
3019    
3020              //We've done the fill character, if the position we're in              //We've done the fill character, if the position we're in
3021              //is one past a comma.  Now, do the ordinary character              //is one past a comma.  Now, do the ordinary character
3022              //corresponding to a digit position.              //corresponding to a digit position.
3023              if (virtual_index < num_digits)              if (virtual_index < num_digits)
3024                 {                 {
3025                 //The character we should print exists.                 //The character we should print exists.
3026                 fprintf(s, "%c", digits[num_digits - 1 - virtual_index]);                 fprintf(s, "%c", digits[num_digits - 1 - virtual_index]);
3027                 }                 }
3028              else              else
3029                 {                 {
3030                 //The character doesn't exist, because the error                 //The character doesn't exist, because the error
3031                 //string is apparently too short.  Must print a                 //string is apparently too short.  Must print a
3032                 //space, instead.                 //space, instead.
3033                 fprintf(s, " ");                 fprintf(s, " ");
3034                 }                 }
3035              }              }
3036           else if (GMP_INTS_mpz_is_zero(i))           else if (GMP_INTS_mpz_is_zero(i))
3037              {              {
3038              /* This is the zero case.  For zero, there is only one line,              /* This is the zero case.  For zero, there is only one line,
3039              ** and every character except the last one is a blank.              ** and every character except the last one is a blank.
3040              */              */
3041              if (cur_digit_on_line == (digits_per_line - 1))              if (cur_digit_on_line == (digits_per_line - 1))
3042                 {                 {
3043                 fprintf(s, "0");                 fprintf(s, "0");
3044                 }                 }
3045              else              else
3046                 {                 {
3047                 fprintf(s, " ");                 fprintf(s, " ");
3048                 }                 }
3049              }              }
3050           else           else
3051              {              {
3052              /* This is a valid number which is not zero.  Need to print              /* This is a valid number which is not zero.  Need to print
3053              ** the digits.              ** the digits.
3054              */              */
3055    
3056              if (idx_into_string < num_digits)              if (idx_into_string < num_digits)
3057                 {                 {
3058                 int actual_index;                 int actual_index;
3059    
3060                 actual_index = num_digits - 1 - idx_into_string;                 actual_index = num_digits - 1 - idx_into_string;
3061                    //This is a string reversal mapping.  The original                    //This is a string reversal mapping.  The original
3062                    //code stored strings least significant digit first,                    //code stored strings least significant digit first,
3063                    //but this code uses most significant digit first.                    //but this code uses most significant digit first.
3064                 assert((actual_index >= 0) && (actual_index < num_digits));                 assert((actual_index >= 0) && (actual_index < num_digits));
3065                 fprintf(s, "%c", digits[actual_index]);                 fprintf(s, "%c", digits[actual_index]);
3066                 }                 }
3067              else              else
3068                 {                 {
3069                 fprintf(s, " ");                 fprintf(s, " ");
3070                 }                 }
3071              } /* End of digit case.              } /* End of digit case.
3072    
3073           /* Now handle the commas.  The rules for commas are straightforward.           /* Now handle the commas.  The rules for commas are straightforward.
3074           **    a)NAN never has a comma.           **    a)NAN never has a comma.
3075           **    b)Zeros never have a comma.           **    b)Zeros never have a comma.
3076           **    c)The final line, last digit never has a comma.           **    c)The final line, last digit never has a comma.
3077           **    d)Everything else in multiples of three ...           **    d)Everything else in multiples of three ...
3078           */           */
3079           if (!(idx_into_string % 3) && (idx_into_string))           if (!(idx_into_string % 3) && (idx_into_string))
3080              {              {
3081              if (i->flags)              if (i->flags)
3082                 {                 {
3083                 //fprintf(s, " ");                 //fprintf(s, " ");
3084                 }                 }
3085              else if (!num_digits)              else if (!num_digits)
3086                 {                 {
3087                 fprintf(s, " ");                 fprintf(s, " ");
3088                 }                 }
3089              else              else
3090                 {                 {
3091                 if (idx_into_string < num_digits)                 if (idx_into_string < num_digits)
3092                    {                    {
3093                    fprintf(s, ",");                    fprintf(s, ",");
3094                    }                    }
3095                 else                 else
3096                    {                    {
3097                    fprintf(s, " ");                    fprintf(s, " ");
3098                    }                    }
3099                 }                 }
3100              }              }
3101           } /* End for each digit on the current line. */           } /* End for each digit on the current line. */
3102                
3103        /* For the first line, print out an informative message        /* For the first line, print out an informative message
3104        ** advising of the number of digits.  For all other lines        ** advising of the number of digits.  For all other lines
3105        ** print nothing.        ** print nothing.
3106        */        */
3107        if (!cur_line && !(i->flags))        if (!cur_line && !(i->flags))
3108           {           {
3109           if (nlines == 1)           if (nlines == 1)
3110              fprintf(s, " ");              fprintf(s, " ");
3111    
3112           if (num_digits <= 1)           if (num_digits <= 1)
3113              {              {
3114              fprintf(s, "  (      1 digit )\n");              fprintf(s, "  (      1 digit )\n");
3115              }              }
3116           else if (num_digits < 1000)           else if (num_digits < 1000)
3117              {              {
3118              fprintf(s, "  (%7d digits)\n", num_digits);              fprintf(s, "  (%7d digits)\n", num_digits);
3119              }              }
3120           else           else
3121              {              {
3122              fprintf(s, "  (%3d,%03d digits)\n", num_digits / 1000, num_digits % 1000);              fprintf(s, "  (%3d,%03d digits)\n", num_digits / 1000, num_digits % 1000);
3123              }              }
3124           }           }
3125        else        else
3126           {           {
3127           fprintf(s, "\n");           fprintf(s, "\n");
3128           }           }
3129        } /* End for each line. */        } /* End for each line. */
3130    
3131     //Deallocate the string space.     //Deallocate the string space.
3132     GMP_INTS_free(digits);     GMP_INTS_free(digits);
3133     }     }
3134    
3135    
3136  void GMP_INTS_mpz_arb_int_raw_to_stream(FILE *s,  void GMP_INTS_mpz_arb_int_raw_to_stream(FILE *s,
3137                                          const GMP_INTS_mpz_struct *i)                                          const GMP_INTS_mpz_struct *i)
3138     {     {
3139     int size_reqd;     int size_reqd;
3140     char *digits;     char *digits;
3141    
3142     //Eyeball the input parameters.     //Eyeball the input parameters.
3143     assert(s != NULL);     assert(s != NULL);
3144     assert(i != NULL);     assert(i != NULL);
3145     assert(i->n_allocd > 0);     assert(i->n_allocd > 0);
3146     assert(i->limbs != NULL);     assert(i->limbs != NULL);
3147    
3148     size_reqd = GMP_INTS_mpz_size_in_base_10(i);     size_reqd = GMP_INTS_mpz_size_in_base_10(i);
3149     digits = GMP_INTS_malloc(size_reqd * sizeof(char));     digits = GMP_INTS_malloc(size_reqd * sizeof(char));
3150     GMP_INTS_mpz_to_string(digits, i);     GMP_INTS_mpz_to_string(digits, i);
3151     fprintf(s, "%s", digits);     fprintf(s, "%s", digits);
3152     GMP_INTS_free(digits);     GMP_INTS_free(digits);
3153     }     }
3154    
3155    
3156  //07/24/01:  Passed visual inspection and unit tests.  //07/24/01:  Passed visual inspection and unit tests.
3157  void GMP_INTS_mpz_pow_ui(         GMP_INTS_mpz_struct *result,  void GMP_INTS_mpz_pow_ui(         GMP_INTS_mpz_struct *result,
3158                           const    GMP_INTS_mpz_struct *base,                           const    GMP_INTS_mpz_struct *base,
3159                           unsigned                      exponent)                           unsigned                      exponent)
3160     {     {
3161     GMP_INTS_mpz_struct temp;     GMP_INTS_mpz_struct temp;
3162        //Temporary location to hold the base raised to        //Temporary location to hold the base raised to
3163        //a binary power (repeated squaring).        //a binary power (repeated squaring).
3164    
3165     //Eyeball the input parameters.     //Eyeball the input parameters.
3166     assert(result != NULL);     assert(result != NULL);
3167     assert(result->n_allocd > 0);     assert(result->n_allocd > 0);
3168     assert(result->limbs != NULL);     assert(result->limbs != NULL);
3169     assert(base != NULL);     assert(base != NULL);
3170     assert(base->n_allocd > 0);     assert(base->n_allocd > 0);
3171     assert(base->limbs != NULL);     assert(base->limbs != NULL);
3172    
3173     //For this function, the base and the result may not     //For this function, the base and the result may not
3174     //be the same object.     //be the same object.
3175     assert(result != base);     assert(result != base);
3176    
3177     //If the base is tained, the output is tainted by association.     //If the base is tained, the output is tainted by association.
3178        {        {
3179        int taint;        int taint;
3180    
3181        taint = GMP_INTS_two_op_flags_map(base->flags, 0);        taint = GMP_INTS_two_op_flags_map(base->flags, 0);
3182    
3183        if (taint)        if (taint)
3184           {           {
3185           result->flags = taint;           result->flags = taint;
3186           return;           return;
3187           }           }
3188        }        }
3189    
3190     //Allocate our temporary variable and set it to the base.     //Allocate our temporary variable and set it to the base.
3191     GMP_INTS_mpz_init(&temp);     GMP_INTS_mpz_init(&temp);
3192     GMP_INTS_mpz_copy(&temp, base);     GMP_INTS_mpz_copy(&temp, base);
3193    
3194     //The result begins with the value of 1.     //The result begins with the value of 1.
3195     GMP_INTS_mpz_set_ui(result, 1);     GMP_INTS_mpz_set_ui(result, 1);
3196    
3197     //Loop through, processing each bit of the exponent.  This is a fairly effective     //Loop through, processing each bit of the exponent.  This is a fairly effective
3198     //algorithm, but not the optimal one (Knuth points this out).     //algorithm, but not the optimal one (Knuth points this out).
3199     while (exponent && !result->flags)     while (exponent && !result->flags)
3200        {        {
3201        if (exponent & 0x1)        if (exponent & 0x1)
3202           {           {
3203           GMP_INTS_mpz_mul(result, result, &temp);           GMP_INTS_mpz_mul(result, result, &temp);
3204           }           }
3205    
3206        //Square the temporary variable.  Because squaring of arb integer        //Square the temporary variable.  Because squaring of arb integer
3207        //may be very expensive, the test against 1 (i.e. last iteration)        //may be very expensive, the test against 1 (i.e. last iteration)
3208        //certainly pays for itself.        //certainly pays for itself.
3209        if (exponent != 1)        if (exponent != 1)
3210           GMP_INTS_mpz_mul(&temp, &temp, &temp);           GMP_INTS_mpz_mul(&temp, &temp, &temp);
3211            
3212        exponent >>= 1;        exponent >>= 1;
3213        }        }
3214    
3215     //Deallocate our temporary variable.     //Deallocate our temporary variable.
3216     GMP_INTS_mpz_clear(&temp);     GMP_INTS_mpz_clear(&temp);
3217     }     }
3218    
3219    
3220  void GMP_INTS_mpz_abs(GMP_INTS_mpz_struct *arg)  void GMP_INTS_mpz_abs(GMP_INTS_mpz_struct *arg)
3221     {     {
3222     //Eyeball the input parameter.     //Eyeball the input parameter.
3223     assert(arg != NULL);     assert(arg != NULL);
3224     assert(arg->n_allocd > 0);     assert(arg->n_allocd > 0);
3225     assert(arg->limbs != NULL);     assert(arg->limbs != NULL);
3226    
3227     //Take the absolute value.     //Take the absolute value.
3228     if (arg->size < 0)     if (arg->size < 0)
3229        arg->size = -arg->size;        arg->size = -arg->size;
3230     }     }
3231    
3232    
3233  //07/29/01:  Visual inspection passed.  Seems to work fine--not explicitly unit-tested  //07/29/01:  Visual inspection passed.  Seems to work fine--not explicitly unit-tested
3234  //directly, but was tested from Tcl.  //directly, but was tested from Tcl.
3235  void GMP_INTS_mpz_gcd(GMP_INTS_mpz_struct       *result,  void GMP_INTS_mpz_gcd(GMP_INTS_mpz_struct       *result,
3236                        const GMP_INTS_mpz_struct *arg1,                        const GMP_INTS_mpz_struct *arg1,
3237                        const GMP_INTS_mpz_struct *arg2)                        const GMP_INTS_mpz_struct *arg2)
3238     {     {
3239     GMP_INTS_mpz_struct u, v, q, r;     GMP_INTS_mpz_struct u, v, q, r;
3240     int loop_count;     int loop_count;
3241    
3242     //Eyeball the inputs carefully.     //Eyeball the inputs carefully.
3243     assert(result != NULL);     assert(result != NULL);
3244     assert(result->n_allocd > 0);     assert(result->n_allocd > 0);
3245     assert(result->limbs != NULL);     assert(result->limbs != NULL);
3246     assert(arg1 != NULL);     assert(arg1 != NULL);
3247     assert(arg1->n_allocd > 0);     assert(arg1->n_allocd > 0);
3248     assert(arg1->limbs != NULL);     assert(arg1->limbs != NULL);
3249     assert(arg2 != NULL);     assert(arg2 != NULL);
3250     assert(arg2->n_allocd > 0);     assert(arg2->n_allocd > 0);
3251     assert(arg2->limbs != NULL);     assert(arg2->limbs != NULL);
3252    
3253     //Args are not allowed to be same object.     //Args are not allowed to be same object.
3254     assert(arg1 != arg2);     assert(arg1 != arg2);
3255    
3256     //If either input is error or taint, taint the output.     //If either input is error or taint, taint the output.
3257        {        {
3258        int taint;        int taint;
3259    
3260        taint = GMP_INTS_two_op_flags_map(arg1->flags, arg2->flags);        taint = GMP_INTS_two_op_flags_map(arg1->flags, arg2->flags);
3261    
3262        result->flags = 0;        result->flags = 0;
3263           //"result" starts off with a clean slate.  Must do this           //"result" starts off with a clean slate.  Must do this
3264           //after taint calculation in case locations of arg1 or arg2           //after taint calculation in case locations of arg1 or arg2
3265           //are the same as result.           //are the same as result.
3266        if (taint)        if (taint)
3267           {           {
3268           result->flags = taint;           result->flags = taint;
3269           return;           return;
3270           }           }
3271        }        }
3272    
3273     //If either input is zero, the result is 1.     //If either input is zero, the result is 1.
3274     if (GMP_INTS_mpz_is_zero(arg1) || GMP_INTS_mpz_is_zero(arg2))     if (GMP_INTS_mpz_is_zero(arg1) || GMP_INTS_mpz_is_zero(arg2))
3275        {        {
3276        GMP_INTS_mpz_set_ui(result, 1);        GMP_INTS_mpz_set_ui(result, 1);
3277        return;        return;
3278        }        }
3279    
3280     //Allocate space for locals.     //Allocate space for locals.
3281     GMP_INTS_mpz_init(&u);     GMP_INTS_mpz_init(&u);
3282     GMP_INTS_mpz_init(&v);     GMP_INTS_mpz_init(&v);
3283     GMP_INTS_mpz_init(&q);     GMP_INTS_mpz_init(&q);
3284     GMP_INTS_mpz_init(&r);     GMP_INTS_mpz_init(&r);
3285    
3286     //We are following Knuth Vol 2, p. 337, the modern Euclidian algorithm.     //We are following Knuth Vol 2, p. 337, the modern Euclidian algorithm.
3287     //Note:  There are faster algorithms for GCD, but because I hacked up the     //Note:  There are faster algorithms for GCD, but because I hacked up the
3288     //GMP multiple-precision library so badly, those aren't included.  This one     //GMP multiple-precision library so badly, those aren't included.  This one
3289     //is logically correct but sub-optimal.  Perhaps at a later time faster     //is logically correct but sub-optimal.  Perhaps at a later time faster
3290     //algorithms will be re-included.     //algorithms will be re-included.
3291     //Copy inputs to u and v.     //Copy inputs to u and v.
3292     GMP_INTS_mpz_copy(&u, arg1);     GMP_INTS_mpz_copy(&u, arg1);
3293     GMP_INTS_mpz_copy(&v, arg2);     GMP_INTS_mpz_copy(&v, arg2);
3294        
3295     //Take the absolute value of each argument.  We know that neither is zero,     //Take the absolute value of each argument.  We know that neither is zero,
3296     //but one or both might be negative.     //but one or both might be negative.
3297     GMP_INTS_mpz_abs(&u);     GMP_INTS_mpz_abs(&u);
3298     GMP_INTS_mpz_abs(&v);     GMP_INTS_mpz_abs(&v);
3299    
3300     //Begin Euclid's algorithm.  There are really three possibilities:     //Begin Euclid's algorithm.  There are really three possibilities:
3301     //  a)We terminate normally.     //  a)We terminate normally.
3302     //  b)Somehow we generate a math error and terminate based on flags.     //  b)Somehow we generate a math error and terminate based on flags.
3303     //  c)Due to some unknown error in the math functions, we go on forever,     //  c)Due to some unknown error in the math functions, we go on forever,
3304     //    and the program locks up.     //    and the program locks up.
3305     GMP_INTS_mpz_set_ui(&r, 1);     GMP_INTS_mpz_set_ui(&r, 1);
3306     loop_count = 0;     loop_count = 0;
3307     while (!GMP_INTS_mpz_is_zero(&r) && !q.flags && !r.flags && (loop_count < 100000))     while (!GMP_INTS_mpz_is_zero(&r) && !q.flags && !r.flags && (loop_count < 100000))
3308        {        {
3309        loop_count++;        loop_count++;
3310    
3311        GMP_INTS_mpz_tdiv_qr(&q, &r, &u, &v);        GMP_INTS_mpz_tdiv_qr(&q, &r, &u, &v);
3312        GMP_INTS_mpz_copy(&u, &v);        GMP_INTS_mpz_copy(&u, &v);
3313        GMP_INTS_mpz_copy(&v, &r);        GMP_INTS_mpz_copy(&v, &r);
3314        }        }
3315    
3316     //Let's hope we didn't get out of the loop based on loop count.     //Let's hope we didn't get out of the loop based on loop count.
3317     assert(loop_count != 100000);     assert(loop_count != 100000);
3318    
3319     //u now contains the answer.     //u now contains the answer.
3320     GMP_INTS_mpz_copy(result, &u);     GMP_INTS_mpz_copy(result, &u);
3321    
3322     //Deallocate space for locals.     //Deallocate space for locals.
3323     GMP_INTS_mpz_clear(&u);     GMP_INTS_mpz_clear(&u);
3324     GMP_INTS_mpz_clear(&v);     GMP_INTS_mpz_clear(&v);
3325     GMP_INTS_mpz_clear(&q);     GMP_INTS_mpz_clear(&q);
3326     GMP_INTS_mpz_clear(&r);     GMP_INTS_mpz_clear(&r);
3327     }     }
3328    
3329    
3330  /******************************************************************/  /******************************************************************/
3331  /***  COMPARISON AND SIZING FUNCTIONS   ***************************/  /***  COMPARISON AND SIZING FUNCTIONS   ***************************/
3332  /******************************************************************/  /******************************************************************/
3333  //07/24/01:  Visual inspection only, due to simplicity.  //07/24/01:  Visual inspection only, due to simplicity.
3334  int GMP_INTS_mpz_fits_uint_p (const GMP_INTS_mpz_struct *src)  int GMP_INTS_mpz_fits_uint_p (const GMP_INTS_mpz_struct *src)
3335     {     {
3336     GMP_INTS_size_t size;     GMP_INTS_size_t size;
3337     GMP_INTS_limb_t mpl;     GMP_INTS_limb_t mpl;
3338    
3339     //Eyeball the input parameter.     //Eyeball the input parameter.
3340     assert(src != NULL);     assert(src != NULL);
3341     assert(src->n_allocd > 0);     assert(src->n_allocd > 0);
3342     assert(src->limbs != NULL);     assert(src->limbs != NULL);
3343    
3344     mpl =  src->limbs[0];     mpl =  src->limbs[0];
3345     size = src->size;     size = src->size;
3346     if (size < 0 || size > 1)     if (size < 0 || size > 1)
3347       return(0);       return(0);
3348    
3349     //The following line came from the original GNU code.     //The following line came from the original GNU code.
3350     //It isn't necessary in our case since limbs and ints are     //It isn't necessary in our case since limbs and ints are
3351     //both 32 bits, but it will do no harm.     //both 32 bits, but it will do no harm.
3352     return (mpl <= (~(unsigned int) 0));     return (mpl <= (~(unsigned int) 0));
3353     }     }
3354    
3355    
3356  unsigned GMP_INTS_mpz_get_limb_zero(const GMP_INTS_mpz_struct *src)  unsigned GMP_INTS_mpz_get_limb_zero(const GMP_INTS_mpz_struct *src)
3357     {     {
3358     //Eyeball the inputs.     //Eyeball the inputs.
3359     assert(src != NULL);     assert(src != NULL);
3360     assert(src->n_allocd > 0);     assert(src->n_allocd > 0);
3361     assert(src->limbs != NULL);     assert(src->limbs != NULL);
3362    
3363     if (!src->size)     if (!src->size)
3364        return(0);        return(0);
3365     else     else
3366        return(src->limbs[0]);        return(src->limbs[0]);
3367     }     }
3368    
3369    
3370  //07/24/01:  Visual inspection only.  Understood the comparisons  //07/24/01:  Visual inspection only.  Understood the comparisons
3371  //and seems like they should work, but ... a little beyond my  //and seems like they should work, but ... a little beyond my
3372  //comfort zone without testing.  Trusting GNU on this one ...  //comfort zone without testing.  Trusting GNU on this one ...
3373  int GMP_INTS_mpz_fits_sint_p (const GMP_INTS_mpz_struct *src)  int GMP_INTS_mpz_fits_sint_p (const GMP_INTS_mpz_struct *src)
3374     {     {
3375     GMP_INTS_size_t size;     GMP_INTS_size_t size;
3376     GMP_INTS_limb_t mpl;     GMP_INTS_limb_t mpl;
3377    
3378     //Eyeball the input parameter.     //Eyeball the input parameter.
3379     assert(src != NULL);     assert(src != NULL);
3380     assert(src->n_allocd > 0);     assert(src->n_allocd > 0);
3381     assert(src->limbs != NULL);     assert(src->limbs != NULL);
3382    
3383     mpl  = src->limbs[0];     mpl  = src->limbs[0];
3384     size = src->size;     size = src->size;
3385     if (size > 0)     if (size > 0)
3386        {        {
3387        if (size > 1)        if (size > 1)
3388                return 0;                return 0;
3389        return mpl < ~((~(unsigned int) 0) >> 1);        return mpl < ~((~(unsigned int) 0) >> 1);
3390        }        }
3391     else     else
3392        {        {
3393        if (size < -1)        if (size < -1)
3394               return 0;               return 0;
3395        return mpl <= ~((~(unsigned int) 0) >> 1);        return mpl <= ~((~(unsigned int) 0) >> 1);
3396        }        }
3397     }     }
3398    
3399    
3400  //07/24/01:  Visual inspection only.  One issue that caught  //07/24/01:  Visual inspection only.  One issue that caught
3401  //my eye is that in one place function returned neg value if  //my eye is that in one place function returned neg value if
3402  //< and pos value if >.  Was within spec, but corrected because  //< and pos value if >.  Was within spec, but corrected because
3403  //it concerned me as I often test against -1 and 1.  Seems  //it concerned me as I often test against -1 and 1.  Seems
3404  //to invite accidents.  //to invite accidents.
3405  int GMP_INTS_mpz_cmp (const GMP_INTS_mpz_struct *u,  int GMP_INTS_mpz_cmp (const GMP_INTS_mpz_struct *u,
3406                        const GMP_INTS_mpz_struct *v)                        const GMP_INTS_mpz_struct *v)
3407     {     {
3408     GMP_INTS_size_t usize = u->size;     GMP_INTS_size_t usize = u->size;
3409     GMP_INTS_size_t vsize = v->size;     GMP_INTS_size_t vsize = v->size;
3410     GMP_INTS_size_t size;     GMP_INTS_size_t size;
3411     GMP_INTS_limb_srcptr up, vp;     GMP_INTS_limb_srcptr up, vp;
3412     int cmp;     int cmp;
3413    
3414     //Eyeball the input parameters.     //Eyeball the input parameters.
3415     assert(u != NULL);     assert(u != NULL);
3416     assert(u->n_allocd > 0);     assert(u->n_allocd > 0);
3417     assert(u->limbs != NULL);     assert(u->limbs != NULL);
3418     assert(v != NULL);     assert(v != NULL);
3419     assert(v->n_allocd > 0);     assert(v->n_allocd > 0);
3420     assert(v->limbs != NULL);     assert(v->limbs != NULL);
3421    
3422     if (usize < vsize)     if (usize < vsize)
3423        return(-1);        return(-1);
3424     else if (usize > vsize)     else if (usize > vsize)
3425        return(1);        return(1);
3426    
3427     if (usize == 0)     if (usize == 0)
3428        return(0);        return(0);
3429    
3430     size = GMP_INTS_abs_of_size_t(usize);     size = GMP_INTS_abs_of_size_t(usize);
3431    
3432     up = u->limbs;     up = u->limbs;
3433     vp = v->limbs;     vp = v->limbs;
3434    
3435     cmp = GMP_INTS_mpn_cmp (up, vp, size);     cmp = GMP_INTS_mpn_cmp (up, vp, size);
3436    
3437     if (cmp == 0)     if (cmp == 0)
3438        return(0);        return(0);
3439    
3440     if ((cmp < 0) == (usize < 0))     if ((cmp < 0) == (usize < 0))
3441        return(1);        return(1);
3442     else     else
3443        return(-1);        return(-1);
3444     }     }
3445    
3446    
3447  //07/24/01:  Not visually inspected.  Relying on  //07/24/01:  Not visually inspected.  Relying on
3448  //GNU ...  //GNU ...
3449  int GMP_INTS_mpz_cmp_ui (const GMP_INTS_mpz_struct *u,  int GMP_INTS_mpz_cmp_ui (const GMP_INTS_mpz_struct *u,
3450                           unsigned long int v_digit)                           unsigned long int v_digit)
3451     {     {
3452     GMP_INTS_size_t usize = u->size;     GMP_INTS_size_t usize = u->size;
3453    
3454     //Eyeball the input parameter.     //Eyeball the input parameter.
3455     assert(u != NULL);     assert(u != NULL);
3456     assert(u->n_allocd > 0);     assert(u->n_allocd > 0);
3457     assert(u->limbs != NULL);     assert(u->limbs != NULL);
3458    
3459     if (usize == 0)     if (usize == 0)
3460        return -(v_digit != 0);        return -(v_digit != 0);
3461    
3462     if (usize == 1)     if (usize == 1)
3463        {        {
3464        GMP_INTS_limb_t u_digit;        GMP_INTS_limb_t u_digit;
3465    
3466        u_digit = u->limbs[0];        u_digit = u->limbs[0];
3467        if (u_digit > v_digit)        if (u_digit > v_digit)
3468                return 1;                return 1;
3469        if (u_digit < v_digit)        if (u_digit < v_digit)
3470                return -1;                return -1;
3471        return 0;        return 0;
3472        }        }
3473    
3474     return (usize > 0) ? 1 : -1;     return (usize > 0) ? 1 : -1;
3475     }     }
3476    
3477    
3478  //07/24/01:  Not visually inspected.  Relying on GNU.  //07/24/01:  Not visually inspected.  Relying on GNU.
3479  int GMP_INTS_mpz_cmp_si (const GMP_INTS_mpz_struct *u,  int GMP_INTS_mpz_cmp_si (const GMP_INTS_mpz_struct *u,
3480                           signed long int v_digit)                           signed long int v_digit)
3481     {     {
3482     GMP_INTS_size_t usize = u->size;     GMP_INTS_size_t usize = u->size;
3483     GMP_INTS_size_t vsize;     GMP_INTS_size_t vsize;
3484     GMP_INTS_limb_t u_digit;     GMP_INTS_limb_t u_digit;
3485    
3486     //Eyeball the input parameter.     //Eyeball the input parameter.
3487     assert(u != NULL);     assert(u != NULL);
3488     assert(u->n_allocd > 0);     assert(u->n_allocd > 0);
3489     assert(u->limbs != NULL);     assert(u->limbs != NULL);
3490    
3491     vsize = 0;     vsize = 0;
3492     if (v_digit > 0)     if (v_digit > 0)
3493        vsize = 1;        vsize = 1;
3494     else if (v_digit < 0)     else if (v_digit < 0)
3495        {        {
3496        vsize = -1;        vsize = -1;
3497        v_digit = -v_digit;        v_digit = -v_digit;
3498        }        }
3499    
3500     if (usize != vsize)     if (usize != vsize)
3501        return usize - vsize;        return usize - vsize;
3502    
3503     if (usize == 0)     if (usize == 0)
3504        return 0;        return 0;
3505    
3506     u_digit = u->limbs[0];     u_digit = u->limbs[0];
3507    
3508     if (u_digit == (GMP_INTS_limb_t) (unsigned long) v_digit)     if (u_digit == (GMP_INTS_limb_t) (unsigned long) v_digit)
3509        return 0;        return 0;
3510    
3511     if (u_digit > (GMP_INTS_limb_t) (unsigned long) v_digit)     if (u_digit > (GMP_INTS_limb_t) (unsigned long) v_digit)
3512        return usize;        return usize;
3513     else     else
3514        return -usize;        return -usize;
3515     }     }
3516    
3517    
3518  //07/24/01:  Not visually inspected.  Counting on GNU.  //07/24/01:  Not visually inspected.  Counting on GNU.
3519  int GMP_INTS_mpz_cmpabs (const GMP_INTS_mpz_struct *u,  int GMP_INTS_mpz_cmpabs (const GMP_INTS_mpz_struct *u,
3520                           const GMP_INTS_mpz_struct *v)                           const GMP_INTS_mpz_struct *v)
3521     {     {
3522     GMP_INTS_size_t usize = u->size;     GMP_INTS_size_t usize = u->size;
3523     GMP_INTS_size_t vsize = v->size;     GMP_INTS_size_t vsize = v->size;
3524     GMP_INTS_limb_srcptr up, vp;     GMP_INTS_limb_srcptr up, vp;
3525     int cmp;     int cmp;
3526    
3527     //Eyeball the input parameters.     //Eyeball the input parameters.
3528     assert(u != NULL);     assert(u != NULL);
3529     assert(u->n_allocd > 0);     assert(u->n_allocd > 0);
3530     assert(u->limbs != NULL);     assert(u->limbs != NULL);
3531     assert(v != NULL);     assert(v != NULL);
3532     assert(v->n_allocd > 0);     assert(v->n_allocd > 0);
3533     assert(v->limbs != NULL);     assert(v->limbs != NULL);
3534    
3535     usize = GMP_INTS_abs_of_size_t(usize);     usize = GMP_INTS_abs_of_size_t(usize);
3536     vsize = GMP_INTS_abs_of_size_t(vsize);     vsize = GMP_INTS_abs_of_size_t(vsize);
3537    
3538     if (usize != vsize)     if (usize != vsize)
3539        return usize - vsize;        return usize - vsize;
3540    
3541     if (usize == 0)     if (usize == 0)
3542        return 0;        return 0;
3543    
3544     up = u->limbs;     up = u->limbs;
3545     vp = v->limbs;     vp = v->limbs;
3546    
3547     cmp = GMP_INTS_mpn_cmp (up, vp, usize);     cmp = GMP_INTS_mpn_cmp (up, vp, usize);
3548    
3549     return cmp;     return cmp;
3550     }     }
3551    
3552  //07/24/01:  Not visually inspected.  Counting on GNU.  //07/24/01:  Not visually inspected.  Counting on GNU.
3553  int GMP_INTS_mpz_cmpabs_ui(const GMP_INTS_mpz_struct *u,  int GMP_INTS_mpz_cmpabs_ui(const GMP_INTS_mpz_struct *u,
3554                             unsigned long int v_digit)                             unsigned long int v_digit)
3555     {     {
3556     GMP_INTS_size_t usize = u->size;     GMP_INTS_size_t usize = u->size;
3557    
3558     //Eyeball the input parameter.     //Eyeball the input parameter.
3559     assert(u != NULL);     assert(u != NULL);
3560     assert(u->n_allocd > 0);     assert(u->n_allocd > 0);
3561     assert(u->limbs != NULL);     assert(u->limbs != NULL);
3562    
3563     if (usize == 0)     if (usize == 0)
3564        return -(v_digit != 0);        return -(v_digit != 0);
3565    
3566     usize = GMP_INTS_abs_of_size_t(usize);     usize = GMP_INTS_abs_of_size_t(usize);
3567    
3568     if (usize == 1)     if (usize == 1)
3569        {        {
3570        GMP_INTS_limb_t u_digit;        GMP_INTS_limb_t u_digit;
3571    
3572        u_digit = u->limbs[0];        u_digit = u->limbs[0];
3573        if (u_digit > v_digit)        if (u_digit > v_digit)
3574                return 1;                return 1;
3575        if (u_digit < v_digit)        if (u_digit < v_digit)
3576                return -1;                return -1;
3577        return 0;        return 0;
3578        }        }
3579    
3580     return 1;     return 1;
3581     }     }
3582    
3583    
3584  /******************************************************************/  /******************************************************************/
3585  /***  VERSION CONTROL IDENTITY FUNCTIONS   ************************/  /***  VERSION CONTROL IDENTITY FUNCTIONS   ************************/
3586  /******************************************************************/  /******************************************************************/
3587    
3588  //07/18/01:  Visual inspection only.  Function deemed too  //07/18/01:  Visual inspection only.  Function deemed too
3589  //simple for unit testing.  //simple for unit testing.
3590  const char *GMP_INTS_cvcinfo(void)  const char *GMP_INTS_cvcinfo(void)
3591     {     {
3592     return("$Header$");     return("$Header$");
3593     }     }
3594    
3595    
3596  //07/18/01:  Visual inspection only.  Function deemed too  //07/18/01:  Visual inspection only.  Function deemed too
3597  //simple for unit testing.  //simple for unit testing.
3598  const char *GMP_INTS_hvcinfo(void)  const char *GMP_INTS_hvcinfo(void)
3599     {     {
3600     return(GMP_INTS_H_VERSION);     return(GMP_INTS_H_VERSION);
3601     }     }
3602    
3603  //End of gmp_ints.c.  //End of gmp_ints.c.

Legend:
Removed from v.70  
changed lines
  Added in v.71

dashley@gmail.com
ViewVC Help
Powered by ViewVC 1.1.25