Context Navigation

source: git/Singular/longrat.cc @ 43ad8d

fieker-DuValspielwiese

Last change on this file since 43ad8d was 38cfbb, checked in by Hans Schönemann <hannes@…>, 26 years ago
* hannes: some changes for factory 1.3 git-svn-id: file:///usr/local/Singular/svn/trunk@1009 2c84dea3-7e68-4137-9b89-c4e89433aadc
Property mode set to `100644`
File size: 76.0 KB

Line
1	/****************************************
2	* Computer Algebra System SINGULAR *
3	****************************************/
4	/* $Id: longrat.cc,v 1.15 1998-01-12 17:32:47 Singular Exp $ */
5	/*
6	* ABSTRACT: computation with long rational numbers (Hubert Grassmann)
7	*/
8
9	#include <string.h>
10	#include "mod2.h"
11	#include "tok.h"
12	#include "mmemory.h"
13	#include "febase.h"
14	#include "numbers.h"
15	#include "modulop.h"
16	#include "longrat.h"
17
18	#ifndef BYTES_PER_MP_LIMB
19	#ifdef HAVE_LIBGMP2
20	#define BYTES_PER_MP_LIMB sizeof(mp_limb_t)
21	#else
22	#define BYTES_PER_MP_LIMB sizeof(mp_limb)
23	#endif
24	#endif
25
26	static number nlMapP(number from)
27	{
28	number to;
29	to = nlInit(npInt(from));
30	return to;
31	}
32
33	BOOLEAN nlSetMap(int c, char ** par, int nop, number minpol)
34	{
35	if (c == 0)
36	{
37	nMap = nlCopy; /Q -> Q/
38	return TRUE;
39	}
40	if (c>1)
41	{
42	if (par==NULL)
43	{
44	npPrimeM=c;
45	nMap = nlMapP; /* Z/p -> Q */
46	return TRUE;
47	}
48	else
49	return FALSE; /* GF(q) ->Q */
50	}
51	else if (c<0)
52	return FALSE; /* Z/p(a) -> Q*/
53	else if (c==1)
54	return FALSE; /* Q(a) -> Q */
55	return FALSE;
56	}
57
58	/-----------------------------------------------------------------/
59	#ifdef HAVE_GMP
60	/*3
61	* parameter s in number:
62	* 0 (or FALSE): not normalised rational
63	* 1 (or TRUE): normalised rational
64	* 3 : integer with n==NULL
65	*/
66	/*3
67	** 'SR_INT' is the type of those integers small enough to fit into 29 bits.
68	** Therefor the value range of this small integers is: $-2^{28}...2^{28}-1$.
69	**
70	** Small integers are represented by an immediate integer handle, containing
71	** the value instead of pointing to it, which has the following form:
72	**
73	** +-------+-------+-------+-------+- - - -+-------+-------+-------+
74	** \| guard \| sign \| bit \| bit \| \| bit \| tag \| tag \|
75	** \| bit \| bit \| 27 \| 26 \| \| 0 \| 0 \| 1 \|
76	** +-------+-------+-------+-------+- - - -+-------+-------+-------+
77	**
78	** Immediate integers handles carry the tag 'SR_INT', i.e. the last bit is 1.
79	** This distuingishes immediate integers from other handles which point to
80	** structures aligned on 4 byte boundaries and therefor have last bit zero.
81	** (The second bit is reserved as tag to allow extensions of this scheme.)
82	** Using immediates as pointers and dereferencing them gives address errors.
83	**
84	** To aid overflow check the most significant two bits must always be equal,
85	** that is to say that the sign bit of immediate integers has a guard bit.
86	**
87	** The macros 'INT_TO_SR' and 'SR_TO_INT' should be used to convert between
88	** a small integer value and its representation as immediate integer handle.
89	**
90	** Large integers and rationals are represented by z and n
91	** where n may be undefined (if s==3)
92	** NULL represents only deleted values
93	*/
94	#define SR_HDL(A) ((long)(A))
95	/#define SR_INT 1/
96	/#define INT_TO_SR(INT) ((number) (((long)INT << 2) + SR_INT))/
97	#define SR_TO_INT(SR) (((long)SR) >> 2)
98
99	#define MP_SMALL 1
100	//#define mpz_isNeg(A) (mpz_cmp_si(A,(long)0)<0)
101	#ifdef HAVE_LIBGMP1
102	#define mpz_isNeg(A) ((A)->size<0)
103	#define mpz_limb_size(A) ((A)->size)
104	#define mpz_limb_d(A) ((A)->d)
105	#define MPZ_DIV(A,B,C) mpz_div((A),(B),(C))
106	#define MPZ_EXACTDIV(A,B,C) mpz_div((A),(B),(C))
107	#else
108	#define mpz_isNeg(A) ((A)->_mp_size<0)
109	#define mpz_limb_size(A) ((A)->_mp_size)
110	#define mpz_limb_d(A) ((A)->_mp_d)
111	#define MPZ_DIV(A,B,C) mpz_tdiv_q((A),(B),(C))
112	#ifdef HAVE_SMALLGMP
113	#define MPZ_EXACTDIV(A,B,C) mpz_tdiv_q((A),(B),(C))
114	#else
115	#define MPZ_EXACTDIV(A,B,C) mpz_divexact((A),(B),(C))
116	#endif
117	#define nlGmpSimple(A)
118	#endif
119
120
121	#ifdef LDEBUG
122	#define nlTest(a) nlDBTest(a,__FILE__,__LINE__)
123	BOOLEAN nlDBTest(number a, char *f,int l)
124	{
125	if (a==NULL)
126	{
127	Print("!!longrat: NULL in %s:%d\n",f,l);
128	return FALSE;
129	}
130	//if ((int)a==1) Print("!! 0x1 as number ? %s %d\n",f,l);
131	if ((((int)a)&3)==3)
132	{
133	Print(" !!longrat:ptr(3) in %s:%d\n",f,l);
134	return FALSE;
135	}
136	if ((((int)a)&3)==1)
137	{
138	if (((((int)a)<<1)>>1)!=((int)a))
139	{
140	Print(" !!longrat:arith:%x in %s:%d\n",(int)a, f,l);
141	return FALSE;
142	}
143	return TRUE;
144	}
145	#ifdef MDEBUG
146	mmDBTestBlock(a,sizeof(*a),f,l);
147	#endif
148	if (a->debug!=123456)
149	{
150	Print("!!longrat:debug:%d in %s:%d\n",a->debug,f,l);
151	a->debug=123456;
152	return FALSE;
153	}
154	if ((a->s<0)\|\|(a->s>4))
155	{
156	Print("!!longrat:s=%d in %s:%d\n",a->s,f,l);
157	return FALSE;
158	}
159	#ifdef MDEBUG
160	#ifdef HAVE_LIBGMP2
161	mmDBTestBlock(a->z._mp_d,a->z._mp_alloc*BYTES_PER_MP_LIMB,f,l);
162	if (a->z._mp_alloc==0)
163	#else
164	mmDBTestBlock(a->z.d,a->z.alloc*BYTES_PER_MP_LIMB,f,l);
165	if(a->z.alloc==0)
166	#endif
167	Print("!!longrat:z->alloc=0 in %s:%l\n",f,l);
168	#endif
169	if (a->s<2)
170	{
171	#ifdef MDEBUG
172	#ifdef HAVE_LIBGMP2
173	mmDBTestBlock(a->n._mp_d,a->n._mp_alloc*BYTES_PER_MP_LIMB,f,-l);
174	if (a->z._mp_alloc==0)
175	#else
176	mmDBTestBlock(a->n.d,a->n.alloc*BYTES_PER_MP_LIMB,f,-l);
177	if(a->z.alloc==0)
178	#endif
179	Print("!!longrat:n->alloc=0 in %s:%l\n",f,l);
180	#endif
181	if (mpz_cmp_si(&a->n,(long)1)==0)
182	{
183	Print("!!longrat:integer as rational in %s:%d\n",f,l);
184	mpz_clear(&a->n);
185	a->s=3;
186	return FALSE;
187	}
188	else if (mpz_isNeg(&a->n))
189	{
190	Print("!!longrat:div. by negative in %s:%d\n",f,l);
191	mpz_neg(&a->z,&a->z);
192	mpz_neg(&a->n,&a->n);
193	return FALSE;
194	}
195	return TRUE;
196	}
197	if (a->s==2)
198	{
199	Print("!!longrat:s=2 in %s:%d\n",f,l);
200	return FALSE;
201	}
202	if (mpz_size1(&a->z)>MP_SMALL) return TRUE;
203	int ui=(int)mpz_get_si(&a->z);
204	if ((((ui<<3)>>3)==ui)
205	&& (mpz_cmp_si(&a->z,(long)ui)==0))
206	{
207	Print("!!longrat:im int %d in %s:%d\n",ui,f,l);
208	f=NULL;
209	return FALSE;
210	}
211	return TRUE;
212	}
213	#endif
214
215	void nlNew (number * r)
216	{
217	*r=NULL;
218	}
219
220	/*2
221	* z := i
222	*/
223	#ifdef MDEBUG
224	#define nlRInit(A) nlRDBInit(A,__LINE__)
225	static inline number nlRDBInit(int i, int n)
226	#else
227	static inline number nlRInit (int i)
228	#endif
229	{
230	#ifdef MDEBUG
231	number z=(number)mmDBAllocBlock(sizeof(rnumber),"longrat.cc",n);
232	#else
233	number z=(number)Alloc(sizeof(rnumber));
234	#endif
235	#ifdef LDEBUG
236	z->debug=123456;
237	#endif
238	mpz_init_set_si(&z->z,(long)i);
239	z->s = 3;
240	return z;
241	}
242
243	number nlInit (int i)
244	{
245	number n;
246	if ( ((i << 3) >> 3) == i ) n=INT_TO_SR(i);
247	else n=nlRInit(i);
248	#ifdef LDEBUG
249	nlTest(n);
250	#endif
251	return n;
252	}
253
254
255	number nlInit (number u)
256	{
257	if (u->s == 3 && mpz_size1(&u->z)<=MP_SMALL)
258	{
259	int ui=(int)mpz_get_si(&u->z);
260	if ((((ui<<3)>>3)==ui)
261	&& (mpz_cmp_si(&u->z,(long)ui)==0))
262	{
263	mpz_clear(&u->z);
264	Free((ADDRESS)u,sizeof(rnumber));
265	return INT_TO_SR(ui);
266	}
267	}
268	return u;
269	}
270
271	int nlSize(number a)
272	{
273	if (a==INT_TO_SR(0))
274	return 0; /* rational 0*/
275	if (SR_HDL(a) & SR_INT)
276	return 1; /* immidiate int */
277	#ifdef HAVE_LIBGMP2
278	int s=a->z._mp_alloc;
279	#else
280	int s=a->z.alloc;
281	#endif
282	if (a->s<2)
283	{
284	#ifdef HAVE_LIBGMP2
285	s+=a->n._mp_alloc;
286	#else
287	s+=a->n.alloc;
288	#endif
289	}
290	return s;
291	}
292
293	/*
294	* delete leading 0s from a gmp number
295	*/
296	#ifdef HAVE_LIBGMP1
297	static void nlGmpSimple(MP_INT *z)
298	{
299	int k;
300	if (mpz_limb_size(z)<0)
301	{
302	k=-mpz_limb_size(z)-1;
303	while ((k>0) && (mpz_limb_d(z)[k]==0))
304	{ k--;printf("X"); }
305	k++;
306	mpz_limb_size(z)=-k;
307	}
308	else
309	{
310	k=mpz_limb_size(z)-1;
311	while ((k>0) && (mpz_limb_d(z)[k]==0))
312	{ k--;printf("X"); }
313	k++;
314	mpz_limb_size(z)=k;
315	}
316	}
317	#endif
318
319	/*2
320	* convert number to int
321	*/
322	int nlInt(number &i)
323	{
324	#ifdef LDEBUG
325	nlTest(i);
326	#endif
327	nlNormalize(i);
328	if (SR_HDL(i) &SR_INT) return SR_TO_INT(i);
329	nlGmpSimple(&i->z);
330	if ((i->s!=3)\|\|(mpz_size1(&i->z)>MP_SMALL)) return 0;
331	int ul=(int)mpz_get_si(&i->z);
332	if (mpz_cmp_si(&i->z,(long)ul)!=0) return 0;
333	return ul;
334	}
335
336	/*2
337	* delete a
338	*/
339	#ifdef LDEBUG
340	void nlDBDelete (number * a,char *f, int l)
341	#else
342	void nlDelete (number * a)
343	#endif
344	{
345	if (*a!=NULL)
346	{
347	#ifdef LDEBUG
348	nlTest(*a);
349	#endif
350	if ((SR_HDL(*a) & SR_INT)==0)
351	{
352	switch ((*a)->s)
353	{
354	case 0:
355	case 1:
356	mpz_clear(&(*a)->n);
357	case 3:
358	mpz_clear(&(*a)->z);
359	#ifdef LDEBUG
360	(*a)->s=2;
361	(*a)->debug=654321;
362	#endif
363	}
364	#if defined(MDEBUG) && defined(LDEBUG)
365	mmDBFreeBlock((ADDRESS) *a,sizeof(rnumber),f,l);
366	#else
367	Free((ADDRESS) *a,sizeof(rnumber));
368	#endif
369	}
370	*a=NULL;
371	}
372	}
373
374	/*2
375	* copy a to b
376	*/
377	number nlCopy(number a)
378	{
379	number b;
380	if ((SR_HDL(a) & SR_INT)\|\|(a==NULL))
381	{
382	return a;
383	}
384	#ifdef LDEBUG
385	nlTest(a);
386	#endif
387	b=(number)Alloc(sizeof(rnumber));
388	#ifdef LDEBUG
389	b->debug=123456;
390	#endif
391	switch (a->s)
392	{
393	case 0:
394	case 1:
395	nlGmpSimple(&a->n);
396	mpz_init_set(&b->n,&a->n);
397	case 3:
398	nlGmpSimple(&a->z);
399	mpz_init_set(&b->z,&a->z);
400	break;
401	}
402	b->s = a->s;
403	#ifdef LDEBUG
404	nlTest(b);
405	#endif
406	return b;
407	}
408
409	/*2
410	* za:= - za
411	*/
412	number nlNeg (number a)
413	{
414	#ifdef LDEBUG
415	nlTest(a);
416	#endif
417	if(SR_HDL(a) &SR_INT)
418	{
419	int r=SR_TO_INT(a);
420	if (r==(-(1<<28))) a=nlRInit(1<<28);
421	else a=INT_TO_SR(-r);
422	}
423	else
424	{
425	mpz_neg(&a->z,&a->z);
426	if ((a->s==3) && (mpz_size1(&a->z)<=MP_SMALL))
427	{
428	int ui=(int)mpz_get_si(&a->z);
429	if ((((ui<<3)>>3)==ui)
430	&& (mpz_cmp_si(&a->z,(long)ui)==0))
431	{
432	mpz_clear(&a->z);
433	Free((ADDRESS)a,sizeof(rnumber));
434	a=INT_TO_SR(ui);
435	}
436	}
437	}
438	#ifdef LDEBUG
439	nlTest(a);
440	#endif
441	return a;
442	}
443
444	/*
445	* 1/a
446	*/
447	number nlInvers(number a)
448	{
449	#ifdef LDEBUG
450	nlTest(a);
451	#endif
452	number n;
453	if (SR_HDL(a) & SR_INT)
454	{
455	if ((a==INT_TO_SR(1)) \|\| (a==INT_TO_SR(-1)))
456	{
457	return a;
458	}
459	if (nlIsZero(a))
460	{
461	WerrorS("div. 1/0");
462	return INT_TO_SR(0);
463	}
464	n=(number)Alloc(sizeof(rnumber));
465	#ifdef LDEBUG
466	n->debug=123456;
467	#endif
468	n->s=1;
469	if ((int)a>0)
470	{
471	mpz_init_set_si(&n->z,(long)1);
472	mpz_init_set_si(&n->n,(long)SR_TO_INT(a));
473	}
474	else
475	{
476	mpz_init_set_si(&n->z,(long)-1);
477	mpz_init_set_si(&n->n,(long)-SR_TO_INT(a));
478	#ifdef LDEBUG
479	nlTest(n);
480	#endif
481	}
482	#ifdef LDEBUG
483	nlTest(n);
484	#endif
485	return n;
486	}
487	n=(number)Alloc(sizeof(rnumber));
488	#ifdef LDEBUG
489	n->debug=123456;
490	#endif
491	{
492	n->s=a->s;
493	mpz_init_set(&n->n,&a->z);
494	switch (a->s)
495	{
496	case 0:
497	case 1:
498	mpz_init_set(&n->z,&a->n);
499	if (mpz_isNeg(&n->n)) /* && n->s<2*/
500	{
501	mpz_neg(&n->z,&n->z);
502	mpz_neg(&n->n,&n->n);
503	}
504	if (mpz_cmp_si(&n->n,(long)1)==0)
505	{
506	mpz_clear(&n->n);
507	n->s=3;
508	if (mpz_size1(&n->z)<=MP_SMALL)
509	{
510	int ui=(int)mpz_get_si(&n->z);
511	if ((((ui<<3)>>3)==ui)
512	&& (mpz_cmp_si(&n->z,(long)ui)==0))
513	{
514	mpz_clear(&n->z);
515	Free((ADDRESS)n,sizeof(rnumber));
516	n=INT_TO_SR(ui);
517	}
518	}
519	}
520	break;
521	case 3:
522	n->s=1;
523	if (mpz_isNeg(&n->n)) /* && n->s<2*/
524	{
525	mpz_neg(&n->n,&n->n);
526	mpz_init_set_si(&n->z,(long)-1);
527	}
528	else
529	{
530	mpz_init_set_si(&n->z,(long)1);
531	}
532	break;
533	}
534	}
535	#ifdef LDEBUG
536	nlTest(n);
537	#endif
538	return n;
539	}
540
541	/*2
542	* u:= a + b
543	*/
544	number nlAdd (number a, number b)
545	{
546	number u;
547	if (SR_HDL(a) & SR_HDL(b) & SR_INT)
548	{
549	int r=SR_HDL(a)+SR_HDL(b)-1;
550	if ( ((r << 1) >> 1) == r )
551	{
552	return (number)r;
553	}
554	u=(number)Alloc(sizeof(rnumber));
555	u->s=3;
556	mpz_init_set_si(&u->z,(long)SR_TO_INT(r));
557	#ifdef LDEBUG
558	u->debug=123456;
559	nlTest(u);
560	#endif
561	return u;
562	}
563	u=(number)Alloc(sizeof(rnumber));
564	#ifdef LDEBUG
565	u->debug=123456;
566	#endif
567	mpz_init(&u->z);
568	if (SR_HDL(b) & SR_INT)
569	{
570	number x=a;
571	a=b;
572	b=x;
573	}
574	if (SR_HDL(a) & SR_INT)
575	{
576	switch (b->s)
577	{
578	case 0:
579	case 1:/* a:short, b:1 */
580	{
581	MP_INT x;
582	mpz_init(&x);
583	if ((int)a>0)
584	{
585	mpz_mul_ui(&x,&b->n,SR_TO_INT(a));
586	}
587	else
588	{
589	mpz_mul_ui(&x,&b->n,-SR_TO_INT(a));
590	mpz_neg(&x,&x);
591	}
592	mpz_add(&u->z,&b->z,&x);
593	nlGmpSimple(&u->z);
594	mpz_clear(&x);
595	if (mpz_cmp_ui(&u->z,(long)0)==0)
596	{
597	mpz_clear(&u->z);
598	Free((ADDRESS)u,sizeof(rnumber));
599	return INT_TO_SR(0);
600	}
601	if (mpz_cmp(&u->z,&b->n)==0)
602	{
603	mpz_clear(&u->z);
604	Free((ADDRESS)u,sizeof(rnumber));
605	return INT_TO_SR(1);
606	}
607	mpz_init_set(&u->n,&b->n);
608	u->s = 0;
609	break;
610	}
611	case 3:
612	{
613	if ((int)a>0)
614	mpz_add_ui(&u->z,&b->z,SR_TO_INT(a));
615	else
616	mpz_sub_ui(&u->z,&b->z,-SR_TO_INT(a));
617	nlGmpSimple(&u->z);
618	if (mpz_cmp_ui(&u->z,(long)0)==0)
619	{
620	mpz_clear(&u->z);
621	Free((ADDRESS)u,sizeof(rnumber));
622	return INT_TO_SR(0);
623	}
624	//u->s = 3;
625	if (mpz_size1(&u->z)<=MP_SMALL)
626	{
627	int ui=(int)mpz_get_si(&u->z);
628	if ((((ui<<3)>>3)==ui)
629	&& (mpz_cmp_si(&u->z,(long)ui)==0))
630	{
631	mpz_clear(&u->z);
632	Free((ADDRESS)u,sizeof(rnumber));
633	return INT_TO_SR(ui);
634	}
635	}
636	u->s = 3;
637	break;
638	}
639	}
640	}
641	else
642	{
643	switch (a->s)
644	{
645	case 0:
646	case 1:
647	{
648	switch(b->s)
649	{
650	case 0:
651	case 1:
652	{
653	MP_INT x;
654	MP_INT y;
655	mpz_init(&x);
656	mpz_init(&y);
657	mpz_mul(&x,&b->z,&a->n);
658	mpz_mul(&y,&a->z,&b->n);
659	mpz_add(&u->z,&x,&y);
660	nlGmpSimple(&u->z);
661	mpz_clear(&x);
662	mpz_clear(&y);
663	if (mpz_cmp_ui(&u->z,(long)0)==0)
664	{
665	mpz_clear(&u->z);
666	Free((ADDRESS)u,sizeof(rnumber));
667	return INT_TO_SR(0);
668	}
669	mpz_init(&u->n);
670	mpz_mul(&u->n,&a->n,&b->n);
671	nlGmpSimple(&u->n);
672	if (mpz_cmp(&u->z,&u->n)==0)
673	{
674	mpz_clear(&u->z);
675	mpz_clear(&u->n);
676	Free((ADDRESS)u,sizeof(rnumber));
677	return INT_TO_SR(1);
678	}
679	u->s = 0;
680	break;
681	}
682	case 3: /* a:1 b:3 */
683	{
684	MP_INT x;
685	mpz_init(&x);
686	mpz_mul(&x,&b->z,&a->n);
687	mpz_add(&u->z,&a->z,&x);
688	nlGmpSimple(&u->z);
689	mpz_clear(&x);
690	if (mpz_cmp_ui(&u->z,(long)0)==0)
691	{
692	mpz_clear(&u->z);
693	Free((ADDRESS)u,sizeof(rnumber));
694	return INT_TO_SR(0);
695	}
696	if (mpz_cmp(&u->z,&a->n)==0)
697	{
698	mpz_clear(&u->z);
699	Free((ADDRESS)u,sizeof(rnumber));
700	return INT_TO_SR(1);
701	}
702	mpz_init_set(&u->n,&a->n);
703	u->s = 0;
704	break;
705	}
706	} /switch (b->s) /
707	break;
708	}
709	case 3:
710	{
711	switch(b->s)
712	{
713	case 0:
714	case 1:/* a:3, b:1 */
715	{
716	MP_INT x;
717	mpz_init(&x);
718	mpz_mul(&x,&a->z,&b->n);
719	mpz_add(&u->z,&b->z,&x);
720	nlGmpSimple(&u->z);
721	mpz_clear(&x);
722	if (mpz_cmp_ui(&u->z,(long)0)==0)
723	{
724	mpz_clear(&u->z);
725	Free((ADDRESS)u,sizeof(rnumber));
726	return INT_TO_SR(0);
727	}
728	if (mpz_cmp(&u->z,&b->n)==0)
729	{
730	mpz_clear(&u->z);
731	Free((ADDRESS)u,sizeof(rnumber));
732	return INT_TO_SR(1);
733	}
734	mpz_init_set(&u->n,&b->n);
735	u->s = 0;
736	break;
737	}
738	case 3:
739	{
740	mpz_add(&u->z,&a->z,&b->z);
741	nlGmpSimple(&u->z);
742	if (mpz_cmp_ui(&u->z,(long)0)==0)
743	{
744	mpz_clear(&u->z);
745	Free((ADDRESS)u,sizeof(rnumber));
746	return INT_TO_SR(0);
747	}
748	if (mpz_size1(&u->z)<=MP_SMALL)
749	{
750	int ui=(int)mpz_get_si(&u->z);
751	if ((((ui<<3)>>3)==ui)
752	&& (mpz_cmp_si(&u->z,(long)ui)==0))
753	{
754	mpz_clear(&u->z);
755	Free((ADDRESS)u,sizeof(rnumber));
756	return INT_TO_SR(ui);
757	}
758	}
759	u->s = 3;
760	break;
761	}
762	}
763	break;
764	}
765	}
766	}
767	#ifdef LDEBUG
768	nlTest(u);
769	#endif
770	return u;
771	}
772
773	/*2
774	* u:= a - b
775	*/
776	number nlSub (number a, number b)
777	{
778	number u;
779	if (SR_HDL(a) & SR_HDL(b) & SR_INT)
780	{
781	int r=SR_HDL(a)-SR_HDL(b)+1;
782	if ( ((r << 1) >> 1) == r )
783	{
784	return (number)r;
785	}
786	u=(number)Alloc(sizeof(rnumber));
787	u->s=3;
788	mpz_init_set_si(&u->z,(long)SR_TO_INT(r));
789	#ifdef LDEBUG
790	u->debug=123456;
791	nlTest(u);
792	#endif
793	return u;
794	}
795	u=(number)Alloc(sizeof(rnumber));
796	#ifdef LDEBUG
797	u->debug=123456;
798	#endif
799	mpz_init(&u->z);
800	if (SR_HDL(a) & SR_INT)
801	{
802	switch (b->s)
803	{
804	case 0:
805	case 1:/* a:short, b:1 */
806	{
807	MP_INT x;
808	mpz_init(&x);
809	if ((int)a>0)
810	{
811	mpz_mul_ui(&x,&b->n,SR_TO_INT(a));
812	}
813	else
814	{
815	mpz_mul_ui(&x,&b->n,-SR_TO_INT(a));
816	mpz_neg(&x,&x);
817	}
818	mpz_sub(&u->z,&x,&b->z);
819	mpz_clear(&x);
820	if (mpz_cmp_ui(&u->z,(long)0)==0)
821	{
822	mpz_clear(&u->z);
823	Free((ADDRESS)u,sizeof(rnumber));
824	return INT_TO_SR(0);
825	}
826	if (mpz_cmp(&u->z,&b->n)==0)
827	{
828	mpz_clear(&u->z);
829	Free((ADDRESS)u,sizeof(rnumber));
830	return INT_TO_SR(1);
831	}
832	mpz_init_set(&u->n,&b->n);
833	u->s = 0;
834	break;
835	}
836	case 3:
837	{
838	if ((int)a>0)
839	{
840	mpz_sub_ui(&u->z,&b->z,SR_TO_INT(a));
841	mpz_neg(&u->z,&u->z);
842	}
843	else
844	{
845	mpz_add_ui(&u->z,&b->z,-SR_TO_INT(a));
846	mpz_neg(&u->z,&u->z);
847	}
848	nlGmpSimple(&u->z);
849	if (mpz_cmp_ui(&u->z,(long)0)==0)
850	{
851	mpz_clear(&u->z);
852	Free((ADDRESS)u,sizeof(rnumber));
853	return INT_TO_SR(0);
854	}
855	//u->s = 3;
856	if (mpz_size1(&u->z)<=MP_SMALL)
857	{
858	int ui=(int)mpz_get_si(&u->z);
859	if ((((ui<<3)>>3)==ui)
860	&& (mpz_cmp_si(&u->z,(long)ui)==0))
861	{
862	mpz_clear(&u->z);
863	Free((ADDRESS)u,sizeof(rnumber));
864	return INT_TO_SR(ui);
865	}
866	}
867	u->s = 3;
868	break;
869	}
870	}
871	}
872	else if (SR_HDL(b) & SR_INT)
873	{
874	switch (a->s)
875	{
876	case 0:
877	case 1:/* b:short, a:1 */
878	{
879	MP_INT x;
880	mpz_init(&x);
881	if ((int)b>0)
882	{
883	mpz_mul_ui(&x,&a->n,SR_TO_INT(b));
884	}
885	else
886	{
887	mpz_mul_ui(&x,&a->n,-SR_TO_INT(b));
888	mpz_neg(&x,&x);
889	}
890	mpz_sub(&u->z,&a->z,&x);
891	mpz_clear(&x);
892	if (mpz_cmp_ui(&u->z,(long)0)==0)
893	{
894	mpz_clear(&u->z);
895	Free((ADDRESS)u,sizeof(rnumber));
896	return INT_TO_SR(0);
897	}
898	if (mpz_cmp(&u->z,&a->n)==0)
899	{
900	mpz_clear(&u->z);
901	Free((ADDRESS)u,sizeof(rnumber));
902	return INT_TO_SR(1);
903	}
904	mpz_init_set(&u->n,&a->n);
905	u->s = 0;
906	break;
907	}
908	case 3:
909	{
910	if ((int)b>0)
911	{
912	mpz_sub_ui(&u->z,&a->z,SR_TO_INT(b));
913	}
914	else
915	{
916	mpz_add_ui(&u->z,&a->z,-SR_TO_INT(b));
917	}
918	if (mpz_cmp_ui(&u->z,(long)0)==0)
919	{
920	mpz_clear(&u->z);
921	Free((ADDRESS)u,sizeof(rnumber));
922	return INT_TO_SR(0);
923	}
924	//u->s = 3;
925	nlGmpSimple(&u->z);
926	if (mpz_size1(&u->z)<=MP_SMALL)
927	{
928	int ui=(int)mpz_get_si(&u->z);
929	if ((((ui<<3)>>3)==ui)
930	&& (mpz_cmp_si(&u->z,(long)ui)==0))
931	{
932	mpz_clear(&u->z);
933	Free((ADDRESS)u,sizeof(rnumber));
934	return INT_TO_SR(ui);
935	}
936	}
937	u->s = 3;
938	break;
939	}
940	}
941	}
942	else
943	{
944	switch (a->s)
945	{
946	case 0:
947	case 1:
948	{
949	switch(b->s)
950	{
951	case 0:
952	case 1:
953	{
954	MP_INT x;
955	MP_INT y;
956	mpz_init(&x);
957	mpz_init(&y);
958	mpz_mul(&x,&b->z,&a->n);
959	mpz_mul(&y,&a->z,&b->n);
960	mpz_sub(&u->z,&y,&x);
961	mpz_clear(&x);
962	mpz_clear(&y);
963	if (mpz_cmp_ui(&u->z,(long)0)==0)
964	{
965	mpz_clear(&u->z);
966	Free((ADDRESS)u,sizeof(rnumber));
967	return INT_TO_SR(0);
968	}
969	mpz_init(&u->n);
970	mpz_mul(&u->n,&a->n,&b->n);
971	if (mpz_cmp(&u->z,&u->n)==0)
972	{
973	mpz_clear(&u->z);
974	mpz_clear(&u->n);
975	Free((ADDRESS)u,sizeof(rnumber));
976	return INT_TO_SR(1);
977	}
978	u->s = 0;
979	break;
980	}
981	case 3: /* a:1, b:3 */
982	{
983	MP_INT x;
984	mpz_init(&x);
985	mpz_mul(&x,&b->z,&a->n);
986	mpz_sub(&u->z,&a->z,&x);
987	mpz_clear(&x);
988	if (mpz_cmp_ui(&u->z,(long)0)==0)
989	{
990	mpz_clear(&u->z);
991	Free((ADDRESS)u,sizeof(rnumber));
992	return INT_TO_SR(0);
993	}
994	if (mpz_cmp(&u->z,&a->n)==0)
995	{
996	mpz_clear(&u->z);
997	Free((ADDRESS)u,sizeof(rnumber));
998	return INT_TO_SR(1);
999	}
1000	mpz_init_set(&u->n,&a->n);
1001	u->s = 0;
1002	break;
1003	}
1004	}
1005	break;
1006	}
1007	case 3:
1008	{
1009	switch(b->s)
1010	{
1011	case 0:
1012	case 1: /* a:3, b:1 */
1013	{
1014	MP_INT x;
1015	mpz_init(&x);
1016	mpz_mul(&x,&a->z,&b->n);
1017	mpz_sub(&u->z,&x,&b->z);
1018	mpz_clear(&x);
1019	if (mpz_cmp_ui(&u->z,(long)0)==0)
1020	{
1021	mpz_clear(&u->z);
1022	Free((ADDRESS)u,sizeof(rnumber));
1023	return INT_TO_SR(0);
1024	}
1025	if (mpz_cmp(&u->z,&b->n)==0)
1026	{
1027	mpz_clear(&u->z);
1028	Free((ADDRESS)u,sizeof(rnumber));
1029	return INT_TO_SR(1);
1030	}
1031	mpz_init_set(&u->n,&b->n);
1032	u->s = 0;
1033	break;
1034	}
1035	case 3: /* a:3 , b:3 */
1036	{
1037	mpz_sub(&u->z,&a->z,&b->z);
1038	nlGmpSimple(&u->z);
1039	if (mpz_cmp_ui(&u->z,(long)0)==0)
1040	{
1041	mpz_clear(&u->z);
1042	Free((ADDRESS)u,sizeof(rnumber));
1043	return INT_TO_SR(0);
1044	}
1045	//u->s = 3;
1046	if (mpz_size1(&u->z)<=MP_SMALL)
1047	{
1048	int ui=(int)mpz_get_si(&u->z);
1049	if ((((ui<<3)>>3)==ui)
1050	&& (mpz_cmp_si(&u->z,(long)ui)==0))
1051	{
1052	mpz_clear(&u->z);
1053	Free((ADDRESS)u,sizeof(rnumber));
1054	return INT_TO_SR(ui);
1055	}
1056	}
1057	u->s = 3;
1058	break;
1059	}
1060	}
1061	break;
1062	}
1063	}
1064	}
1065	#ifdef LDEBUG
1066	nlTest(u);
1067	#endif
1068	return u;
1069	}
1070
1071	/*2
1072	* u := a * b
1073	*/
1074	number nlMult (number a, number b)
1075	{
1076	number u;
1077
1078	#ifdef LDEBUG
1079	nlTest(a);
1080	nlTest(b);
1081	#endif
1082	if (SR_HDL(a) & SR_HDL(b) & SR_INT)
1083	{
1084	if (b==INT_TO_SR(0)) return INT_TO_SR(0);
1085	int r=(SR_HDL(a)-1)*(SR_HDL(b)>>1);
1086	if ((r/(SR_HDL(b)>>1))==(SR_HDL(a)-1))
1087	{
1088	u=((number) ((r>>1)+SR_INT));
1089	if (((SR_HDL(u)<<1)>>1)==SR_HDL(u)) return (u);
1090	return nlRInit(SR_HDL(u)>>2);
1091	}
1092	u=(number)Alloc(sizeof(rnumber));
1093	#ifdef LDEBUG
1094	u->debug=123456;
1095	#endif
1096	u->s=3;
1097	if ((int)b>0)
1098	{
1099	mpz_init_set_si(&u->z,(long)SR_TO_INT(a));
1100	mpz_mul_ui(&u->z,&u->z,(unsigned long)SR_TO_INT(b));
1101	}
1102	// else if ((int)a>=0)
1103	// {
1104	// mpz_init_set_si(&u->z,(long)SR_TO_INT(b));
1105	// mpz_mul_ui(&u->z,&u->z,(unsigned long)SR_TO_INT(a));
1106	// }
1107	else
1108	{
1109	mpz_init_set_si(&u->z,(long)(-SR_TO_INT(a)));
1110	mpz_mul_ui(&u->z,&u->z,(long)(-SR_TO_INT(b)));
1111	}
1112	#ifdef LDEBUG
1113	nlTest(u);
1114	#endif
1115	return u;
1116	}
1117	else
1118	{
1119	if ((a==INT_TO_SR(0))
1120	\|\|(b==INT_TO_SR(0)))
1121	return INT_TO_SR(0);
1122	u=(number)Alloc(sizeof(rnumber));
1123	#ifdef LDEBUG
1124	u->debug=123456;
1125	#endif
1126	mpz_init(&u->z);
1127	if (SR_HDL(b) & SR_INT)
1128	{
1129	number x=a;
1130	a=b;
1131	b=x;
1132	}
1133	if (SR_HDL(a) & SR_INT)
1134	{
1135	u->s=b->s;
1136	if (u->s==1) u->s=0;
1137	if ((int)a>0)
1138	{
1139	mpz_mul_ui(&u->z,&b->z,(unsigned long)SR_TO_INT(a));
1140	}
1141	else
1142	{
1143	if (a==INT_TO_SR(-1))
1144	{
1145	mpz_set(&u->z,&b->z);
1146	mpz_neg(&u->z,&u->z);
1147	u->s=b->s;
1148	}
1149	else
1150	{
1151	mpz_mul_ui(&u->z,&b->z,(unsigned long)-SR_TO_INT(a));
1152	mpz_neg(&u->z,&u->z);
1153	}
1154	}
1155	nlGmpSimple(&u->z);
1156	if (u->s<2)
1157	{
1158	if (mpz_cmp(&u->z,&b->n)==0)
1159	{
1160	mpz_clear(&u->z);
1161	Free((ADDRESS)u,sizeof(rnumber));
1162	return INT_TO_SR(1);
1163	}
1164	mpz_init_set(&u->n,&b->n);
1165	}
1166	else //u->s==3
1167	{
1168	if (mpz_size1(&u->z)<=MP_SMALL)
1169	{
1170	int ui=(int)mpz_get_si(&u->z);
1171	if ((((ui<<3)>>3)==ui)
1172	&& (mpz_cmp_si(&u->z,(long)ui)==0))
1173	{
1174	mpz_clear(&u->z);
1175	Free((ADDRESS)u,sizeof(rnumber));
1176	return INT_TO_SR(ui);
1177	}
1178	}
1179	}
1180	}
1181	else
1182	{
1183	mpz_mul(&u->z,&a->z,&b->z);
1184	u->s = 0;
1185	if(a->s==3)
1186	{
1187	if(b->s==3)
1188	{
1189	u->s = 3;
1190	}
1191	else
1192	{
1193	if (mpz_cmp(&u->z,&b->n)==0)
1194	{
1195	mpz_clear(&u->z);
1196	Free((ADDRESS)u,sizeof(rnumber));
1197	return INT_TO_SR(1);
1198	}
1199	mpz_init_set(&u->n,&b->n);
1200	}
1201	}
1202	else
1203	{
1204	if(b->s==3)
1205	{
1206	if (mpz_cmp(&u->z,&a->n)==0)
1207	{
1208	mpz_clear(&u->z);
1209	Free((ADDRESS)u,sizeof(rnumber));
1210	return INT_TO_SR(1);
1211	}
1212	mpz_init_set(&u->n,&a->n);
1213	}
1214	else
1215	{
1216	mpz_init(&u->n);
1217	mpz_mul(&u->n,&a->n,&b->n);
1218	if (mpz_cmp(&u->z,&u->n)==0)
1219	{
1220	mpz_clear(&u->z);
1221	mpz_clear(&u->n);
1222	Free((ADDRESS)u,sizeof(rnumber));
1223	return INT_TO_SR(1);
1224	}
1225	}
1226	}
1227	}
1228	}
1229	#ifdef LDEBUG
1230	nlTest(u);
1231	#endif
1232	return u;
1233	}
1234
1235	/*2
1236	* u := a / b in Z, if b \| a (else undefined)
1237	*/
1238	number nlExactDiv(number a, number b)
1239	{
1240	if (b==INT_TO_SR(0))
1241	{
1242	WerrorS("div. by 0");
1243	return INT_TO_SR(0);
1244	}
1245	number u;
1246	if (SR_HDL(a) & SR_HDL(b) & SR_INT)
1247	{
1248	/* the small int -(1<<28) divided by -1 is the large int (1<<28) */
1249	if ((a==INT_TO_SR(-(1<<28)))&&(b==INT_TO_SR(-1)))
1250	{
1251	return nlRInit(1<<28);
1252	}
1253	int aa=SR_TO_INT(a);
1254	int bb=SR_TO_INT(b);
1255	return INT_TO_SR(aa/bb);
1256	}
1257	number bb=NULL;
1258	if (SR_HDL(b) & SR_INT)
1259	{
1260	bb=nlRInit(SR_TO_INT(b));
1261	b=bb;
1262	}
1263	u=(number)Alloc(sizeof(rnumber));
1264	#ifdef LDEBUG
1265	u->debug=123456;
1266	#endif
1267	mpz_init(&u->z);
1268	/* u=a/b */
1269	u->s = 3;
1270	MPZ_EXACTDIV(&u->z,&a->z,&b->z);
1271	if (bb!=NULL)
1272	{
1273	mpz_clear(&bb->z);
1274	Free((ADDRESS)bb,sizeof(rnumber));
1275	}
1276	if (mpz_size1(&u->z)<=MP_SMALL)
1277	{
1278	int ui=(int)mpz_get_si(&u->z);
1279	if ((((ui<<3)>>3)==ui)
1280	&& (mpz_cmp_si(&u->z,(long)ui)==0))
1281	{
1282	mpz_clear(&u->z);
1283	Free((ADDRESS)u,sizeof(rnumber));
1284	u=INT_TO_SR(ui);
1285	}
1286	}
1287	#ifdef LDEBUG
1288	nlTest(u);
1289	#endif
1290	return u;
1291	}
1292
1293	/*2
1294	* u := a / b in Z
1295	*/
1296	number nlIntDiv (number a, number b)
1297	{
1298	if (b==INT_TO_SR(0))
1299	{
1300	WerrorS("div. by 0");
1301	return INT_TO_SR(0);
1302	}
1303	number u;
1304	if (SR_HDL(a) & SR_HDL(b) & SR_INT)
1305	{
1306	/* the small int -(1<<28) divided by -1 is the large int (1<<28) */
1307	if ((a==INT_TO_SR(-(1<<28)))&&(b==INT_TO_SR(-1)))
1308	{
1309	return nlRInit(1<<28);
1310	}
1311	/* consider the signs of a and b:
1312	* + + -> u=a+b-1
1313	* + - -> u=a-b-1
1314	* - + -> u=a-b+1
1315	* - - -> u=a+b+1
1316	*/
1317	int aa=SR_TO_INT(a);
1318	int bb=SR_TO_INT(b);
1319	if (aa>0)
1320	{
1321	if (bb>0)
1322	return INT_TO_SR((aa+bb-1)/bb);
1323	else
1324	return INT_TO_SR((aa-bb-1)/bb);
1325	}
1326	else
1327	{
1328	if (bb>0)
1329	{
1330	return INT_TO_SR((aa-bb+1)/bb);
1331	}
1332	else
1333	{
1334	return INT_TO_SR((aa+bb+1)/bb);
1335	}
1336	}
1337	}
1338	if (SR_HDL(a) & SR_INT)
1339	{
1340	/* the small int -(1<<28) divided by 2^28 is 1 */
1341	if (a==INT_TO_SR(-(1<<28)))
1342	{
1343	if(mpz_cmp_si(&b->z,(long)(1<<28))==0)
1344	{
1345	return INT_TO_SR(-1);
1346	}
1347	}
1348	/* a is a small and b is a large int: -> 0 */
1349	return INT_TO_SR(0);
1350	}
1351	number bb=NULL;
1352	if (SR_HDL(b) & SR_INT)
1353	{
1354	bb=nlRInit(SR_TO_INT(b));
1355	b=bb;
1356	}
1357	u=(number)Alloc(sizeof(rnumber));
1358	#ifdef LDEBUG
1359	u->debug=123456;
1360	#endif
1361	mpz_init_set(&u->z,&a->z);
1362	/* consider the signs of a and b:
1363	* + + -> u=a+b-1
1364	* + - -> u=a-b-1
1365	* - + -> u=a-b+1
1366	* - - -> u=a+b+1
1367	*/
1368	if (mpz_isNeg(&a->z))
1369	{
1370	if (mpz_isNeg(&b->z))
1371	{
1372	mpz_add(&u->z,&u->z,&b->z);
1373	}
1374	else
1375	{
1376	mpz_sub(&u->z,&u->z,&b->z);
1377	}
1378	mpz_add_ui(&u->z,&u->z,1);
1379	}
1380	else
1381	{
1382	if (mpz_isNeg(&b->z))
1383	{
1384	mpz_sub(&u->z,&u->z,&b->z);
1385	}
1386	else
1387	{
1388	mpz_add(&u->z,&u->z,&b->z);
1389	}
1390	mpz_sub_ui(&u->z,&u->z,1);
1391	}
1392	/* u=u/b */
1393	u->s = 3;
1394	MPZ_DIV(&u->z,&u->z,&b->z);
1395	nlGmpSimple(&u->z);
1396	if (bb!=NULL)
1397	{
1398	mpz_clear(&bb->z);
1399	Free((ADDRESS)bb,sizeof(rnumber));
1400	}
1401	if (mpz_size1(&u->z)<=MP_SMALL)
1402	{
1403	int ui=(int)mpz_get_si(&u->z);
1404	if ((((ui<<3)>>3)==ui)
1405	&& (mpz_cmp_si(&u->z,(long)ui)==0))
1406	{
1407	mpz_clear(&u->z);
1408	Free((ADDRESS)u,sizeof(rnumber));
1409	u=INT_TO_SR(ui);
1410	}
1411	}
1412	#ifdef LDEBUG
1413	nlTest(u);
1414	#endif
1415	return u;
1416	}
1417
1418	/*2
1419	* u := a mod b in Z, u>=0
1420	*/
1421	number nlIntMod (number a, number b)
1422	{
1423	if (b==INT_TO_SR(0))
1424	{
1425	WerrorS("div. by 0");
1426	return INT_TO_SR(0);
1427	}
1428	if (a==INT_TO_SR(0))
1429	return INT_TO_SR(0);
1430	number u;
1431	if (SR_HDL(a) & SR_HDL(b) & SR_INT)
1432	{
1433	if ((int)a>0)
1434	{
1435	if ((int)b>0)
1436	return INT_TO_SR(SR_TO_INT(a)%SR_TO_INT(b));
1437	else
1438	return INT_TO_SR(SR_TO_INT(a)%(-SR_TO_INT(b)));
1439	}
1440	else
1441	{
1442	if ((int)b>0)
1443	{
1444	int i=(-SR_TO_INT(a))%SR_TO_INT(b);
1445	if ( i != 0 ) i = (SR_TO_INT(b))-i;
1446	return INT_TO_SR(i);
1447	}
1448	else
1449	{
1450	int i=(-SR_TO_INT(a))%(-SR_TO_INT(b));
1451	if ( i != 0 ) i = (-SR_TO_INT(b))-i;
1452	return INT_TO_SR(i);
1453	}
1454	}
1455	}
1456	if (SR_HDL(a) & SR_INT)
1457	{
1458	/* a is a small and b is a large int: -> a or (a+b) or (a-b) */
1459	if ((int)a<0)
1460	{
1461	if (mpz_isNeg(&b->z))
1462	return nlSub(a,b);
1463	/else/
1464	return nlAdd(a,b);
1465	}
1466	/else/
1467	return a;
1468	}
1469	number bb=NULL;
1470	if (SR_HDL(b) & SR_INT)
1471	{
1472	bb=nlRInit(SR_TO_INT(b));
1473	b=bb;
1474	}
1475	u=(number)Alloc(sizeof(rnumber));
1476	#ifdef LDEBUG
1477	u->debug=123456;
1478	#endif
1479	mpz_init(&u->z);
1480	u->s = 3;
1481	mpz_mod(&u->z,&a->z,&b->z);
1482	if (bb!=NULL)
1483	{
1484	mpz_clear(&bb->z);
1485	Free((ADDRESS)bb,sizeof(rnumber));
1486	}
1487	if (mpz_isNeg(&u->z))
1488	{
1489	if (mpz_isNeg(&b->z))
1490	mpz_sub(&u->z,&u->z,&b->z);
1491	else
1492	mpz_add(&u->z,&u->z,&b->z);
1493	}
1494	nlGmpSimple(&u->z);
1495	if (mpz_size1(&u->z)<=MP_SMALL)
1496	{
1497	int ui=(int)mpz_get_si(&u->z);
1498	if ((((ui<<3)>>3)==ui)
1499	&& (mpz_cmp_si(&u->z,(long)ui)==0))
1500	{
1501	mpz_clear(&u->z);
1502	Free((ADDRESS)u,sizeof(rnumber));
1503	u=INT_TO_SR(ui);
1504	}
1505	}
1506	#ifdef LDEBUG
1507	nlTest(u);
1508	#endif
1509	return u;
1510	}
1511
1512	/*2
1513	* u := a / b
1514	*/
1515	number nlDiv (number a, number b)
1516	{
1517	number u;
1518	if (nlIsZero(b))
1519	{
1520	WerrorS("div. by 0");
1521	return INT_TO_SR(0);
1522	}
1523	u=(number)Alloc(sizeof(rnumber));
1524	u->s=0;
1525	#ifdef LDEBUG
1526	u->debug=123456;
1527	#endif
1528	// ---------- short / short ------------------------------------
1529	if (SR_HDL(a) & SR_HDL(b) & SR_INT)
1530	{
1531	int i=SR_TO_INT(a);
1532	int j=SR_TO_INT(b);
1533	int r=i%j;
1534	if (r==0)
1535	{
1536	Free((ADDRESS)u,sizeof(rnumber));
1537	return INT_TO_SR(i/j);
1538	}
1539	mpz_init_set_si(&u->z,(long)i);
1540	mpz_init_set_si(&u->n,(long)j);
1541	}
1542	else
1543	{
1544	mpz_init(&u->z);
1545	// ---------- short / long ------------------------------------
1546	if (SR_HDL(a) & SR_INT)
1547	{
1548	// short a / (z/n) -> (a*n)/z
1549	if (b->s<2)
1550	{
1551	if ((int)a>=0)
1552	mpz_mul_ui(&u->z,&b->n,SR_TO_INT(a));
1553	else
1554	{
1555	mpz_mul_ui(&u->z,&b->n,-SR_TO_INT(a));
1556	mpz_neg(&u->z,&u->z);
1557	}
1558	}
1559	else
1560	// short a / long z -> a/z
1561	{
1562	mpz_set_si(&u->z,SR_TO_INT(a));
1563	}
1564	if (mpz_cmp(&u->z,&b->z)==0)
1565	{
1566	mpz_clear(&u->z);
1567	Free((ADDRESS)u,sizeof(rnumber));
1568	return INT_TO_SR(1);
1569	}
1570	mpz_init_set(&u->n,&b->z);
1571	}
1572	// ---------- long / short ------------------------------------
1573	else if (SR_HDL(b) & SR_INT)
1574	{
1575	mpz_set(&u->z,&a->z);
1576	// (z/n) / b -> z/(n*b)
1577	if (a->s<2)
1578	{
1579	mpz_init_set(&u->n,&a->n);
1580	if ((int)b>0)
1581	mpz_mul_ui(&u->n,&u->n,SR_TO_INT(b));
1582	else
1583	{
1584	mpz_mul_ui(&u->n,&u->n,-SR_TO_INT(b));
1585	mpz_neg(&u->z,&u->z);
1586	}
1587	}
1588	else
1589	// long z / short b -> z/b
1590	{
1591	//mpz_set(&u->z,&a->z);
1592	mpz_init_set_si(&u->n,SR_TO_INT(b));
1593	}
1594	}
1595	// ---------- long / long ------------------------------------
1596	else
1597	{
1598	//u->s=0;
1599	mpz_set(&u->z,&a->z);
1600	mpz_init_set(&u->n,&b->z);
1601	if (a->s<2) mpz_mul(&u->n,&u->n,&a->n);
1602	if (b->s<2) mpz_mul(&u->z,&u->z,&b->n);
1603	}
1604	}
1605	if (mpz_isNeg(&u->n))
1606	{
1607	mpz_neg(&u->z,&u->z);
1608	mpz_neg(&u->n,&u->n);
1609	}
1610	if (mpz_cmp_si(&u->n,(long)1)==0)
1611	{
1612	mpz_clear(&u->n);
1613	if (mpz_size1(&u->z)<=MP_SMALL)
1614	{
1615	int ui=(int)mpz_get_si(&u->z);
1616	if ((((ui<<3)>>3)==ui)
1617	&& (mpz_cmp_si(&u->z,(long)ui)==0))
1618	{
1619	mpz_clear(&u->z);
1620	Free((ADDRESS)u,sizeof(rnumber));
1621	return INT_TO_SR(ui);
1622	}
1623	}
1624	u->s=3;
1625	}
1626	#ifdef LDEBUG
1627	nlTest(u);
1628	#endif
1629	return u;
1630	}
1631
1632	#if (defined(i386)) && (defined(HAVE_LIBGMP1))
1633	/*
1634	* compute x^exp for x in Z
1635	* there seems to be an bug in mpz_pow_ui for i386
1636	*/
1637	static inline void nlPow (MP_INT * res,MP_INT * x,int exp)
1638	{
1639	if (exp==0)
1640	{
1641	mpz_set_si(res,(long)1);
1642	}
1643	else
1644	{
1645	MP_INT xh;
1646	mpz_init(&xh);
1647	mpz_set(res,x);
1648	exp--;
1649	while (exp!=0)
1650	{
1651	mpz_mul(&xh,res,x);
1652	mpz_set(res,&xh);
1653	exp--;
1654	}
1655	mpz_clear(&xh);
1656	}
1657	}
1658	#endif
1659
1660	/*2
1661	* u:= x ^ exp
1662	*/
1663	void nlPower (number x,int exp,number * u)
1664	{
1665	if (!nlIsZero(x))
1666	{
1667	#ifdef LDEBUG
1668	nlTest(x);
1669	#endif
1670	number aa=NULL;
1671	if (SR_HDL(x) & SR_INT)
1672	{
1673	aa=nlRInit(SR_TO_INT(x));
1674	x=aa;
1675	}
1676	*u=(number)Alloc(sizeof(rnumber));
1677	#ifdef LDEBUG
1678	(*u)->debug=123456;
1679	#endif
1680	mpz_init(&(*u)->z);
1681	(*u)->s = x->s;
1682	#if (!defined(i386)) \|\| (defined(HAVE_LIBGMP2))
1683	mpz_pow_ui(&(*u)->z,&x->z,(unsigned long)exp);
1684	#else
1685	nlPow(&(*u)->z,&x->z,exp);
1686	#endif
1687	if (x->s<2)
1688	{
1689	mpz_init(&(*u)->n);
1690	#if (!defined(i386)) \|\| (defined(HAVE_LIBGMP2))
1691	mpz_pow_ui(&(*u)->n,&x->n,(unsigned long)exp);
1692	#else
1693	nlPow(&(*u)->n,&x->n,exp);
1694	#endif
1695	}
1696	if (aa!=NULL)
1697	{
1698	mpz_clear(&aa->z);
1699	Free((ADDRESS)aa,sizeof(rnumber));
1700	}
1701	if (((u)->s<=2) && (mpz_cmp_si(&(u)->n,(long)1)==0))
1702	{
1703	mpz_clear(&(*u)->n);
1704	(*u)->s=3;
1705	}
1706	if (((u)->s==3) && (mpz_size1(&(u)->z)<=MP_SMALL))
1707	{
1708	int ui=(int)mpz_get_si(&(*u)->z);
1709	if ((((ui<<3)>>3)==ui)
1710	&& (mpz_cmp_si(&(*u)->z,(long)ui)==0))
1711	{
1712	mpz_clear(&(*u)->z);
1713	Free((ADDRESS)*u,sizeof(rnumber));
1714	*u=INT_TO_SR(ui);
1715	}
1716	}
1717	}
1718	else
1719	*u = INT_TO_SR(0);
1720	#ifdef LDEBUG
1721	nlTest(*u);
1722	#endif
1723	}
1724
1725	BOOLEAN nlIsZero (number a)
1726	{
1727	#ifdef LDEBUG
1728	nlTest(a);
1729	#endif
1730	//if (a==INT_TO_SR(0)) return TRUE;
1731	//if (SR_HDL(a) & SR_INT) return FALSE;
1732	//number aa=nlCopy(a);
1733	//nlNormalize(aa);
1734	return (a==INT_TO_SR(0));
1735	}
1736
1737	/*2
1738	* za >= 0 ?
1739	*/
1740	BOOLEAN nlGreaterZero (number a)
1741	{
1742	#ifdef LDEBUG
1743	nlTest(a);
1744	#endif
1745	if (SR_HDL(a) & SR_INT) return SR_HDL(a)>=0;
1746	return (!mpz_isNeg(&a->z));
1747	}
1748
1749	/*2
1750	* a > b ?
1751	*/
1752	BOOLEAN nlGreater (number a, number b)
1753	{
1754	#ifdef LDEBUG
1755	nlTest(a);
1756	nlTest(b);
1757	#endif
1758	number r;
1759	BOOLEAN rr;
1760	r=nlSub(a,b);
1761	rr=(!nlIsZero(r)) && (nlGreaterZero(r));
1762	nlDelete(&r);
1763	return rr;
1764	}
1765
1766	/*2
1767	* a = b ?
1768	*/
1769	BOOLEAN nlEqual (number a, number b)
1770	{
1771	#ifdef LDEBUG
1772	nlTest(a);
1773	nlTest(b);
1774	#endif
1775	// short - short
1776	if (SR_HDL(a) & SR_HDL(b) & SR_INT) return a==b;
1777	// long - short
1778	BOOLEAN bo;
1779	if (SR_HDL(b) & SR_INT)
1780	{
1781	if (a->s!=0) return FALSE;
1782	number n=b; b=a; a=n;
1783	}
1784	// short - long
1785	if (SR_HDL(a) & SR_INT)
1786	{
1787	if (b->s!=0)
1788	return FALSE;
1789	if (((int)a > 0) && (mpz_isNeg(&b->z)))
1790	return FALSE;
1791	if (((int)a < 0) && (!mpz_isNeg(&b->z)))
1792	return FALSE;
1793	MP_INT bb;
1794	mpz_init_set(&bb,&b->n);
1795	if ((int)a<0)
1796	{
1797	mpz_neg(&bb,&bb);
1798	mpz_mul_ui(&bb,&bb,(long)-SR_TO_INT(a));
1799	}
1800	else
1801	{
1802	mpz_mul_ui(&bb,&bb,(long)SR_TO_INT(a));
1803	}
1804	bo=(mpz_cmp(&bb,&b->z)==0);
1805	mpz_clear(&bb);
1806	return bo;
1807	}
1808	// long - long
1809	if (((a->s==1) && (b->s==3))
1810	\|\| ((b->s==1) && (a->s==3)))
1811	return FALSE;
1812	if (mpz_isNeg(&a->z)&&(!mpz_isNeg(&b->z)))
1813	return FALSE;
1814	if (mpz_isNeg(&b->z)&&(!mpz_isNeg(&a->z)))
1815	return FALSE;
1816	nlGmpSimple(&a->z);
1817	nlGmpSimple(&b->z);
1818	if (a->s<2)
1819	nlGmpSimple(&a->n);
1820	if (b->s<2)
1821	nlGmpSimple(&b->n);
1822	MP_INT aa;
1823	MP_INT bb;
1824	mpz_init_set(&aa,&a->z);
1825	mpz_init_set(&bb,&b->z);
1826	if (a->s<2) mpz_mul(&bb,&bb,&a->n);
1827	if (b->s<2) mpz_mul(&aa,&aa,&b->n);
1828	bo=(mpz_cmp(&aa,&bb)==0);
1829	mpz_clear(&aa);
1830	mpz_clear(&bb);
1831	return bo;
1832	}
1833
1834	/*2
1835	* a == 1 ?
1836	*/
1837	BOOLEAN nlIsOne (number a)
1838	{
1839	#ifdef LDEBUG
1840	if (a==NULL) return FALSE;
1841	nlTest(a);
1842	#endif
1843	if (SR_HDL(a) & SR_INT) return (a==INT_TO_SR(1));
1844	return FALSE;
1845	}
1846
1847	/*2
1848	* a == -1 ?
1849	*/
1850	BOOLEAN nlIsMOne (number a)
1851	{
1852	#ifdef LDEBUG
1853	if (a==NULL) return FALSE;
1854	nlTest(a);
1855	#endif
1856	if (SR_HDL(a) & SR_INT) return (a==INT_TO_SR(-1));
1857	return FALSE;
1858	}
1859
1860	/*2
1861	* result =gcd(a,b)
1862	*/
1863	number nlGcd(number a, number b)
1864	{
1865	number result;
1866	#ifdef LDEBUG
1867	nlTest(a);
1868	nlTest(b);
1869	#endif
1870	//nlNormalize(a);
1871	//nlNormalize(b);
1872	if ((SR_HDL(a)==5)\|\|(a==INT_TO_SR(-1))
1873	\|\| (SR_HDL(b)==5)\|\|(b==INT_TO_SR(-1))) return INT_TO_SR(1);
1874	if (SR_HDL(a) & SR_HDL(b) & SR_INT)
1875	{
1876	int i=SR_TO_INT(a);
1877	int j=SR_TO_INT(b);
1878	if((i==0)\|\|(j==0))
1879	return INT_TO_SR(1);
1880	int l;
1881	i=ABS(i);
1882	j=ABS(j);
1883	do
1884	{
1885	l=i%j;
1886	i=j;
1887	j=l;
1888	} while (l!=0);
1889	return INT_TO_SR(i);
1890	}
1891	if (((!(SR_HDL(a) & SR_INT))&&(a->s<2))
1892	\|\| ((!(SR_HDL(b) & SR_INT))&&(b->s<2))) return INT_TO_SR(1);
1893	number aa=NULL;
1894	if (SR_HDL(a) & SR_INT)
1895	{
1896	aa=nlRInit(SR_TO_INT(a));
1897	a=aa;
1898	}
1899	else
1900	if (SR_HDL(b) & SR_INT)
1901	{
1902	aa=nlRInit(SR_TO_INT(b));
1903	b=aa;
1904	}
1905	result=(number)Alloc(sizeof(rnumber));
1906	#ifdef LDEBUG
1907	result->debug=123456;
1908	#endif
1909	mpz_init(&result->z);
1910	mpz_gcd(&result->z,&a->z,&b->z);
1911	nlGmpSimple(&result->z);
1912	result->s = 3;
1913	if (mpz_size1(&result->z)<=MP_SMALL)
1914	{
1915	int ui=(int)mpz_get_si(&result->z);
1916	if ((((ui<<3)>>3)==ui)
1917	&& (mpz_cmp_si(&result->z,(long)ui)==0))
1918	{
1919	mpz_clear(&result->z);
1920	Free((ADDRESS)result,sizeof(rnumber));
1921	result=INT_TO_SR(ui);
1922	}
1923	}
1924	if (aa!=NULL)
1925	{
1926	mpz_clear(&aa->z);
1927	Free((ADDRESS)aa,sizeof(rnumber));
1928	}
1929	#ifdef LDEBUG
1930	nlTest(result);
1931	#endif
1932	return result;
1933	}
1934
1935	/*2
1936	* simplify x
1937	*/
1938	void nlNormalize (number &x)
1939	{
1940	#ifdef LDEBUG
1941	nlTest(x);
1942	#endif
1943	if ((SR_HDL(x) & SR_INT) \|\|(x==NULL))
1944	return;
1945	if (x->s==3)
1946	{
1947	if (mpz_cmp_ui(&x->z,(long)0)==0)
1948	{
1949	nlDelete(&x);
1950	x=nlInit(0);
1951	return;
1952	}
1953	if (mpz_size1(&x->z)<=MP_SMALL)
1954	{
1955	int ui=(int)mpz_get_si(&x->z);
1956	if ((((ui<<3)>>3)==ui)
1957	&& (mpz_cmp_si(&x->z,(long)ui)==0))
1958	{
1959	mpz_clear(&x->z);
1960	Free((ADDRESS)x,sizeof(rnumber));
1961	x=INT_TO_SR(ui);
1962	return;
1963	}
1964	}
1965	}
1966	if (x->s!=1) nlGmpSimple(&x->z);
1967	if (x->s==0)
1968	{
1969	BOOLEAN divided=FALSE;
1970	MP_INT gcd;
1971	mpz_init(&gcd);
1972	mpz_gcd(&gcd,&x->z,&x->n);
1973	x->s=1;
1974	if (mpz_cmp_si(&gcd,(long)1)!=0)
1975	{
1976	MP_INT r;
1977	mpz_init(&r);
1978	MPZ_EXACTDIV(&r,&x->z,&gcd);
1979	mpz_set(&x->z,&r);
1980	MPZ_EXACTDIV(&r,&x->n,&gcd);
1981	mpz_set(&x->n,&r);
1982	mpz_clear(&r);
1983	nlGmpSimple(&x->z);
1984	divided=TRUE;
1985	}
1986	mpz_clear(&gcd);
1987	nlGmpSimple(&x->n);
1988	if (mpz_cmp_si(&x->n,(long)1)==0)
1989	{
1990	mpz_clear(&x->n);
1991	if (mpz_size1(&x->z)<=MP_SMALL)
1992	{
1993	int ui=(int)mpz_get_si(&x->z);
1994	if ((((ui<<3)>>3)==ui)
1995	&& (mpz_cmp_si(&x->z,(long)ui)==0))
1996	{
1997	mpz_clear(&x->z);
1998	#ifdef LDEBUG
1999	x->debug=654324;
2000	#endif
2001	Free((ADDRESS)x,sizeof(rnumber));
2002	x=INT_TO_SR(ui);
2003	return;
2004	}
2005	}
2006	x->s=3;
2007	}
2008	else if (divided)
2009	{
2010	_mpz_realloc(&x->n,mpz_size1(&x->n));
2011	}
2012	if (divided) _mpz_realloc(&x->z,mpz_size1(&x->z));
2013	}
2014	#ifdef LDEBUG
2015	nlTest(x);
2016	#endif
2017	}
2018
2019	/*2
2020	* returns in result->z the lcm(a->z,b->n)
2021	*/
2022	number nlLcm(number a, number b)
2023	{
2024	number result;
2025	#ifdef LDEBUG
2026	nlTest(a);
2027	nlTest(b);
2028	#endif
2029	if ((SR_HDL(b) & SR_INT)
2030	\|\| (b->s==3))
2031	{
2032	// b is 1/(b->n) => b->n is 1 => result is a
2033	return nlCopy(a);
2034	}
2035	number aa=NULL;
2036	if (SR_HDL(a) & SR_INT)
2037	{
2038	aa=nlRInit(SR_TO_INT(a));
2039	a=aa;
2040	}
2041	result=(number)Alloc(sizeof(rnumber));
2042	#ifdef LDEBUG
2043	result->debug=123456;
2044	#endif
2045	result->s=3;
2046	MP_INT gcd;
2047	mpz_init(&gcd);
2048	mpz_init(&result->z);
2049	mpz_gcd(&gcd,&a->z,&b->n);
2050	if (mpz_cmp_si(&gcd,(long)1)!=0)
2051	{
2052	MP_INT bt;
2053	mpz_init_set(&bt,&b->n);
2054	MPZ_EXACTDIV(&bt,&bt,&gcd);
2055	mpz_mul(&result->z,&a->z,&bt);
2056	mpz_clear(&bt);
2057	}
2058	else
2059	mpz_mul(&result->z,&a->z,&b->n);
2060	mpz_clear(&gcd);
2061	if (aa!=NULL)
2062	{
2063	mpz_clear(&aa->z);
2064	Free((ADDRESS)aa,sizeof(rnumber));
2065	}
2066	nlGmpSimple(&result->z);
2067	if (mpz_size1(&result->z)<=MP_SMALL)
2068	{
2069	int ui=(int)mpz_get_si(&result->z);
2070	if ((((ui<<3)>>3)==ui)
2071	&& (mpz_cmp_si(&result->z,(long)ui)==0))
2072	{
2073	mpz_clear(&result->z);
2074	Free((ADDRESS)result,sizeof(rnumber));
2075	return INT_TO_SR(ui);
2076	}
2077	}
2078	#ifdef LDEBUG
2079	nlTest(result);
2080	#endif
2081	return result;
2082	}
2083
2084	int nlModP(number n, int p)
2085	{
2086	if (SR_HDL(n) & SR_INT)
2087	{
2088	int i=SR_TO_INT(n);
2089	if (i<0) return (p-((-i)%p));
2090	return i%p;
2091	}
2092	int iz=mpz_mmod_ui(NULL,&n->z,(unsigned long)p);
2093	if (n->s!=3)
2094	{
2095	int in=mpz_mmod_ui(NULL,&n->n,(unsigned long)p);
2096	return (int)npDiv((number)iz,(number)in);
2097	}
2098	return iz;
2099	}
2100
2101	#else
2102
2103	/0 --------------------------------------------implementation/
2104
2105	typedef unsigned int int32;
2106	#define MAX_INT16 0XFFFF
2107	#define SIGNUM 0X8000
2108	#define BASIS 0X10000
2109	#define MOD_BASIS 0XFFFF
2110	#define SHIFT_BASIS 0X10
2111
2112	/*3
2113	* x%SIGNUM == length(x)
2114	* x/SIGNUM == 1 : x < 0
2115	*/
2116	#define NLLENGTH(a) ((int)(*(a)&0X7FFF))
2117	#define NLNEG(a) (*(a)^SIGNUM)
2118	#define NLSIGN(a) (*(a)>SIGNUM)
2119
2120	/*3
2121	*allocate memory for a long integer
2122	*/
2123	static lint nlAllocLint (int L)
2124	{
2125	lint x;
2126
2127	if (L < SIGNUM)
2128	{
2129	x = (lint)Alloc((L+1)*sizeof(int16));
2130	*x = L;
2131	return x;
2132	}
2133	else
2134	{
2135	WerrorS("number too large");
2136	return NULL;
2137	}
2138	}
2139
2140	/*3
2141	*deallocate Memory of a long integer
2142	*/
2143	static void nlDeleteLint (lint x)
2144	{
2145	if (x)
2146	{
2147	Free((ADDRESS)x,(NLLENGTH(x)+1)*sizeof(int16));
2148	}
2149	}
2150
2151	/*3
2152	*nlShrink the memory used by c
2153	*/
2154	static void nlShrink (lint * c)
2155	{
2156	int j, i;
2157	lint cc = *c;
2158
2159	i = j = NLLENGTH(cc);
2160	while ((i!=0) && (cc[i]==0))
2161	{
2162	i--;
2163	}
2164	if (i != j)
2165	{
2166	cc = nlAllocLint(i);
2167	memcpy(cc+1,(c)+1,isizeof(int16));
2168	nlDeleteLint(*c);
2169	*c = cc;
2170	}
2171	}
2172
2173	/*3
2174	* x == 1 ?
2175	*/
2176	static BOOLEAN nlIsOneLint (lint x)
2177	{
2178	return ((*x == 1) && (x[1] == 1));
2179	}
2180
2181	/*3
2182	*compares absolute values
2183	*abs(x) > abs(y) nlCompAbs(x,y) = 1;
2184	*abs(x) < abs(y) -1;
2185	*abs(x) = abs(y) 0;
2186	*/
2187	static int nlCompAbs (lint x, lint y)
2188	{
2189	int i = NLLENGTH(x);
2190
2191	if (i != NLLENGTH(y))
2192	{
2193	if (i > NLLENGTH(y))
2194	{
2195	return 1;
2196	}
2197	else
2198	{
2199	return -1;
2200	}
2201	}
2202	else
2203	{
2204	do
2205	{
2206	if (x[i] > y[i])
2207	{
2208	return 1;
2209	}
2210	if (x[i] < y[i])
2211	{
2212	return -1;
2213	}
2214	i--;
2215	}
2216	while (i!=0);
2217	return 0;
2218	}
2219	}
2220
2221	/*3
2222	*compares long integers
2223	*(x) > (y) nlComp(x,y) = 1;
2224	*(x) < (y) -1;
2225	*(x) = (y) 0;
2226	*/
2227	static int nlComp (lint x, lint y)
2228	{
2229	int i = NLSIGN(x);
2230
2231	if (i != NLSIGN(y))
2232	{
2233	if (i != 0)
2234	{
2235	return -1;
2236	}
2237	else
2238	{
2239	return 1;
2240	}
2241	}
2242	else if (i != 0)
2243	{
2244	return -nlCompAbs(x, y);
2245	}
2246	else
2247	{
2248	return nlCompAbs(x, y);
2249	}
2250	}
2251
2252	/*3
2253	* copy f to t
2254	*/
2255	static lint nlCopyLint (lint f)
2256	{
2257	int i;
2258	lint t;
2259
2260	if (f!=NULL)
2261	{
2262	i = (NLLENGTH(f)+1)*sizeof(int16);
2263	t = (lint)Alloc(i);
2264	memcpy(t, f, i);
2265	return t;
2266	}
2267	else
2268	{
2269	return NULL;
2270	}
2271	}
2272
2273	/*3
2274	* n := 1
2275	*/
2276	static lint nlInitOne()
2277	{
2278	lint n = (lint)Alloc(2*sizeof(int16));
2279	*n = 1;
2280	n[1] = 1;
2281	return n;
2282	}
2283
2284	/*3
2285	* n := i
2286	*/
2287	static lint nlToLint (int16 i)
2288	{
2289	if (i)
2290	{
2291	lint n = (lint)Alloc(2*sizeof(int16));
2292	*n = 1;
2293	n[1] = i;
2294	return n;
2295	}
2296	else
2297	{
2298	return NULL;
2299	}
2300	}
2301
2302	/*3
2303	* Add longs, signs are equal
2304	*/
2305	static lint nlAddBasic (lint a, lint b)
2306	{
2307	int i, max, min;
2308	int32 w;
2309	lint xx,h;
2310
2311	max = NLLENGTH(a);
2312	min = NLLENGTH(b);
2313	if ((min == 1) && (max == 1))
2314	{
2315	w = (int32) a[1] + (int32) b[1];
2316	if (w < BASIS)
2317	{
2318	xx = (lint)Alloc(2*sizeof(int16));
2319	*xx = 1;
2320	xx[1] = (int16) w;
2321	}
2322	else
2323	{
2324	xx = (lint)Alloc(3 * sizeof(int16));
2325	*xx = 2;
2326	xx[1] = (int16) (MOD_BASIS&w);
2327	xx[2] = 1;
2328	}
2329	if (NLSIGN(a))
2330	{
2331	*xx += SIGNUM;
2332	}
2333	return xx;
2334	}
2335	if(min>max)
2336	{
2337	i = max;
2338	max = min;
2339	min = i;
2340	h = b;
2341	}
2342	else
2343	{
2344	h = a;
2345	}
2346	xx = (lint)Alloc((max+1)*sizeof(int16));
2347	*xx = (int16)max;
2348	w = (int32) a[1] + (int32) b[1];
2349	xx[1] = (int16) (MOD_BASIS&w);
2350	for (i=2; i<=min; i++)
2351	{
2352	w = (w>>SHIFT_BASIS) + (int32) a[i] + (int32) b[i];
2353	xx[i] = (int16) (MOD_BASIS&w);
2354	}
2355	i = min;
2356	loop
2357	{
2358	i++;
2359	if ((w<BASIS) \|\| (i>max))
2360	{
2361	break;
2362	}
2363	w = (w>>SHIFT_BASIS) + (int32) h[i];
2364	xx[i] = (int16) (MOD_BASIS&w);
2365	}
2366	for(; i<=max; i++)
2367	{
2368	xx[i] = h[i];
2369	}
2370	if (w>=BASIS)
2371	{
2372	h = (lint)Alloc((max+2)*sizeof(int16));
2373	memcpy(h+1,xx+1,max*sizeof(int16));
2374	max++;
2375	*h = (int16)max;
2376	h[max] = 1;
2377	nlDeleteLint(xx);
2378	xx = h;
2379	}
2380	if (NLSIGN(a))
2381	{
2382	*xx += SIGNUM;
2383	}
2384	return xx;
2385	}
2386
2387	/*3
2388	* Subtract longs, signs are equal, a > b
2389	*/
2390	static lint nlSubBasic (lint a, lint b)
2391	{
2392	int i,j1,j2;
2393	int16 w;
2394	lint xx;
2395	j1 = NLLENGTH(a);
2396	if (j1 == 1)
2397	{
2398	xx = (lint)Alloc(2*sizeof(int16));
2399	*xx = 1;
2400	xx[1] = a[1]-b[1];
2401	if (NLSIGN(a))
2402	{
2403	*xx += SIGNUM;
2404	}
2405	return xx;
2406	}
2407	j2 = NLLENGTH(b);
2408	xx = (lint)Alloc((j1+1)*sizeof(int16));
2409	*xx = (int16)j1;
2410	if(a[1] < b[1])
2411	{
2412	xx[1] = MAX_INT16-b[1]+a[1]+1;
2413	w = 1;
2414	}
2415	else
2416	{
2417	xx[1] = a[1]-b[1];
2418	w = 0;
2419	}
2420	for (i=2; i<=j2; i++)
2421	{
2422	if(a[i]>b[i])
2423	{
2424	xx[i] = a[i]-b[i]-w;
2425	w = 0;
2426	}
2427	else if(a[i]<b[i])
2428	{
2429	if(w!=0)
2430	{
2431	xx[i] = MAX_INT16 - b[i] + a[i];
2432	}
2433	else
2434	{
2435	xx[i] = MAX_INT16 - b[i] + a[i] + 1;
2436	}
2437	w = 1;
2438	}
2439	else
2440	{
2441	if(w!=0)
2442	{
2443	xx[i] = MAX_INT16;
2444	}
2445	else
2446	{
2447	xx[i] = 0;
2448	}
2449	}
2450	}
2451	i = j2;
2452	while(w!=0)
2453	{
2454	i++;
2455	if(a[i])
2456	{
2457	xx[i] = a[i]-w;
2458	w = 0;
2459	}
2460	else
2461	{
2462	xx[i] = MAX_INT16;
2463	}
2464	}
2465	while(i < j1)
2466	{
2467	i++;
2468	xx[i] = a[i];
2469	}
2470	nlShrink(&xx);
2471	if (NLSIGN(a))
2472	{
2473	*xx += SIGNUM;
2474	}
2475	return xx;
2476	}
2477
2478	/*3
2479	* c:= a + b
2480	*/
2481	static lint nlAddLong (lint a,lint b)
2482	{
2483	if (a==NULL)
2484	return nlCopyLint(b);
2485	else if (b==NULL)
2486	return nlCopyLint(a);
2487	if (NLSIGN(a) == NLSIGN(b))
2488	return nlAddBasic(a,b);
2489	else
2490	{
2491	switch(nlCompAbs(a, b))
2492	{
2493	case -1:
2494	return nlSubBasic(b, a);
2495	case 1:
2496	return nlSubBasic(a, b);
2497	default:
2498	;
2499	}
2500	}
2501	return NULL;
2502	}
2503
2504	/*3
2505	* c:= a - b
2506	*/
2507	static lint nlSubLong (lint a,lint b)
2508	{
2509	lint c;
2510	if(b==NULL)
2511	return nlCopyLint(a);
2512	if(a==NULL)
2513	{
2514	c = nlCopyLint(b);
2515	*c = NLNEG(c);
2516	return c;
2517	}
2518	if (NLSIGN(a) != NLSIGN(b))
2519	return nlAddBasic(a,b);
2520	else
2521	{
2522	switch(nlCompAbs(a, b))
2523	{
2524	case 1:
2525	return nlSubBasic(a, b);
2526	case -1:
2527	c = nlSubBasic(b, a);
2528	*c = NLNEG(c);
2529	return c;
2530	default:
2531	;
2532	}
2533	}
2534	return NULL;
2535	}
2536
2537	/*3
2538	* a[i] = a[i] + b[i]*x, i = 1 to length(b)
2539	*/
2540	static void nlAddMult (lint a, lint b, int32 x)
2541	{
2542	int32 w;
2543	int i, l = NLLENGTH(b);
2544	w = (int32) a[1] + (int32) b[1]*x;
2545	a[1] = (int16) (MOD_BASIS&w);
2546	for (i=2; i<=l; i++)
2547	{
2548	w = (w>>SHIFT_BASIS) + (int32) a[i] + (int32) b[i]*x;
2549	a[i] = (int16) (MOD_BASIS&w);
2550	}
2551	if (w>=BASIS)
2552	a[l+1] = (int16)(w>>SHIFT_BASIS);
2553	}
2554
2555	/*3
2556	* a[i] = b[i]*x, i = 1 to length(b)
2557	*/
2558	static void nlIniMult (lint a, lint b, int32 x)
2559	{
2560	int i, l = NLLENGTH(b);
2561	int32 w;
2562	w = (int32) b[1]*x;
2563	a[1] = (int16) (MOD_BASIS&w);
2564	for (i=2; i<=l; i++)
2565	{
2566	w = (w>>SHIFT_BASIS) + (int32) b[i]*x;
2567	a[i] = (int16) (MOD_BASIS&w);
2568	}
2569	if (w>=BASIS)
2570	a[l+1] = (int16)(w>>SHIFT_BASIS);
2571	}
2572
2573	/*3
2574	* mult long, c:= a * b
2575	*/
2576	static lint nlMultLong (lint a,lint b)
2577	{
2578	int la, lb, i;
2579	int32 h;
2580	lint xx;
2581	if ((a==NULL) \|\| (b==NULL))
2582	return NULL;
2583	la = NLLENGTH(a);
2584	lb = NLLENGTH(b);
2585	if ((la == 1) && (lb == 1))
2586	{
2587	h = (int32) a[1] * (int32) b[1];
2588	if (h < BASIS)
2589	{
2590	xx = (lint)Alloc(2*sizeof(int16));
2591	*xx = 1;
2592	xx[1] = (int16) h;
2593	}
2594	else
2595	{
2596	xx = (lint)Alloc(3 * sizeof(int16));
2597	*xx = 2;
2598	xx[1] = (int16) (MOD_BASIS&h);
2599	xx[2] = (int16) (h>>SHIFT_BASIS);
2600	}
2601	if ((NLSIGN(a)) != (NLSIGN(b)))
2602	*xx += SIGNUM;
2603	return xx;
2604	}
2605	i = la+lb;
2606	xx = (lint)Alloc((i+1)*sizeof(int16));
2607	*xx = (int16)i;
2608	memset(xx+1,0,i*sizeof(int16));
2609	if (la < lb)
2610	{
2611	nlIniMult(xx, b, (int32) a[1]);
2612	for (i=1; i<la; i++)
2613	nlAddMult(xx+i, b, (int32) a[i+1]);
2614	}
2615	else
2616	{
2617	nlIniMult(xx, a, (int32) b[1]);
2618	for (i=1; i<lb; i++)
2619	nlAddMult(xx+i, a, (int32) b[i+1]);
2620	}
2621	nlShrink(&xx);
2622	if ((NLSIGN(a)) != (NLSIGN(b)))
2623	*xx += SIGNUM;
2624	return xx;
2625	}
2626
2627	/*3
2628	* aq := aq div b, b small
2629	*/
2630	static void nlQuotSmall (lint * aq, int16 b)
2631	{
2632	int al,j;
2633	int16 a1,c1;
2634	int32 db,da,dc;
2635	lint a,cc;
2636	db = (int32) b;
2637	a = *aq;
2638	al = NLLENGTH(a);
2639	a1 = a[al];
2640	if (al == 1)
2641	{
2642	a[1] = a1 / b;
2643	return;
2644	}
2645	if (a1 >= b)
2646	{
2647	cc = a;
2648	c1 = a1 / b;
2649	cc[al] = c1;
2650	a1 -= b*c1;
2651	}
2652	else
2653	cc = nlAllocLint(al-1);
2654	loop
2655	{
2656	al--;
2657	if (al==0)
2658	{
2659	if(a != cc)
2660	nlDeleteLint(a);
2661	*aq = cc;
2662	return;
2663	}
2664	if (a1==0)
2665	{
2666	j = al;
2667	loop
2668	{
2669	a1 = a[j];
2670	if (a1 >= b)
2671	{
2672	c1 = a1 / b;
2673	cc[j] = c1;
2674	a1 -= b*c1;
2675	}
2676	else
2677	cc[j] = 0;
2678	j--;
2679	if (j==0)
2680	{
2681	if(a != cc)
2682	nlDeleteLint(a);
2683	*aq = cc;
2684	return;
2685	}
2686	else if (a1!=0)
2687	{
2688	al = j;
2689	break;
2690	}
2691	}
2692	}
2693	da = BASIS*(int32) a1 + (int32) a[al];
2694	dc = da / db;
2695	cc[al] = (int16) dc;
2696	da -= db*dc;
2697	a1 = (int16) da;
2698	}
2699	}
2700
2701	/*3
2702	* a := a - b, b1...bj
2703	*/
2704	static void nlSubSimple (lint a, lint b, int j)
2705	{
2706	int i;
2707	int16 w;
2708	if (a[1] < b[1])
2709	{
2710	a[1] += MAX_INT16-b[1]+1;
2711	w = 1;
2712	}
2713	else
2714	{
2715	a[1] -= b[1];
2716	w = 0;
2717	}
2718	for (i=2; i<=j; i++)
2719	{
2720	if(a[i] > b[i])
2721	{
2722	a[i] -= b[i]+w;
2723	w = 0;
2724	}
2725	else if(a[i] < b[i])
2726	{
2727	if (w)
2728	a[i] += MAX_INT16 - b[i];
2729	else
2730	a[i] += MAX_INT16 - b[i] + 1;
2731	w = 1;
2732	}
2733	else
2734	{
2735	if(w)
2736	a[i] = MAX_INT16;
2737	else
2738	a[i] = 0;
2739	}
2740	}
2741	if(w!=0)
2742	a[j+1]--;
2743	}
2744
2745	/*3
2746	* a := a - b*x, b1...bj
2747	*/
2748	static void nlSubMult(lint a, lint b, int32 x, int l)
2749	{
2750	int i;
2751	int16 k, w;
2752	int32 y;
2753	y = (int32)b[1]*x;
2754	k = (int16)(MOD_BASIS&y);
2755	if (a[1] < k)
2756	{
2757	a[1] += MAX_INT16-k+1;
2758	w = 1;
2759	}
2760	else
2761	{
2762	a[1] -= k;
2763	w = 0;
2764	}
2765	for(i=2; i<=l; i++)
2766	{
2767	y = (y>>SHIFT_BASIS) + (int32)b[i]*x;
2768	k = (int16)(MOD_BASIS&y);
2769	if(a[i] > k)
2770	{
2771	a[i] -= k+w;
2772	w = 0;
2773	}
2774	else if(a[i] < k)
2775	{
2776	if (w!=0)
2777	a[i] += MAX_INT16 - k;
2778	else
2779	a[i] += MAX_INT16 - k + 1;
2780	w = 1;
2781	}
2782	else
2783	{
2784	if(w!=0)
2785	a[i] = MAX_INT16;
2786	else
2787	a[i] = 0;
2788	}
2789	}
2790	w += (int16)(y>>SHIFT_BASIS);
2791	if(w!=0)
2792	a[l+1] -= w;
2793	}
2794
2795	/*3
2796	* compare xi and yi from i=l to 1
2797	*/
2798	static BOOLEAN nlCompHigh (lint x, lint y, int l)
2799	{
2800	int i;
2801	for (i=l; i!=0; i--)
2802	{
2803	if (x[i] > y[i])
2804	return TRUE;
2805	else if (x[i] < y[i])
2806	return FALSE;
2807	}
2808	return TRUE;
2809	}
2810
2811	/*3
2812	* aq := aq div bn, bn long
2813	*/
2814	static void nlQuotLong (lint * aq, lint bn)
2815	{
2816	int bl,al,si,j;
2817	int16 b1,a1,a2;
2818	int32 db,db1,dc,dc1,dc2;
2819	lint an,cc,ansi;
2820	BOOLEAN rel;
2821	bl = NLLENGTH(bn);
2822	b1 = bn[bl];
2823	db = BASIS*(int32)b1 + (int32)bn[bl-1];
2824	if ((bl > 2) && ((b1 < MAX_INT16) \|\| (bn[bl-1] < MAX_INT16)))
2825	db++;
2826	db1 = (int32)(b1+1);
2827	an = *aq;
2828	al = NLLENGTH(an);
2829	a1 = an[al];
2830	a2 = an[al-1];
2831	si = al-bl;
2832	ansi = an+si;
2833	rel = nlCompHigh(ansi, bn, bl);
2834	if (rel)
2835	{
2836	cc = nlAllocLint(si+1);
2837	dc = (BASIS*(int32)a1 + (int32)a2) / db;
2838	if (dc > 1)
2839	{
2840	nlSubMult(ansi, bn, dc, bl);
2841	rel = nlCompHigh(ansi, bn, bl);
2842	}
2843	else
2844	dc = 0;
2845	while (rel)
2846	{
2847	dc++;
2848	nlSubSimple(ansi, bn, bl);
2849	rel = nlCompHigh(ansi, bn, bl);
2850	}
2851	cc[si+1] = (int16)dc;
2852	a1 = an[al];
2853	a2 = an[al-1];
2854	}
2855	else
2856	cc = nlAllocLint(si);
2857	loop
2858	{
2859	al--;
2860	if (al < bl)
2861	{
2862	nlDeleteLint(an);
2863	*aq = cc;
2864	return;
2865	}
2866	si--;
2867	ansi--;
2868	if (a1==0)
2869	{
2870	j = al;
2871	loop
2872	{
2873	a1 = a2;
2874	a2 = an[j-1];
2875	rel = nlCompHigh(ansi, bn, bl);
2876	if (rel)
2877	{
2878	dc = (BASIS*(int32)a1 + (int32)a2) / db;
2879	if (dc > 1)
2880	{
2881	nlSubMult(ansi, bn, dc, bl);
2882	rel = nlCompHigh(ansi, bn, bl);
2883	}
2884	else
2885	dc = 0;
2886	while (rel)
2887	{
2888	dc++;
2889	nlSubSimple(ansi, bn, bl);
2890	rel = nlCompHigh(ansi, bn, bl);
2891	}
2892	cc[si+1] = (int16)dc;
2893	a1 = an[j];
2894	a2 = an[j-1];
2895	}
2896	else
2897	cc[si+1] = 0;
2898	si--;
2899	ansi--;
2900	j--;
2901	if (j < bl)
2902	{
2903	nlDeleteLint(an);
2904	*aq = cc;
2905	return;
2906	}
2907	else if(a1!=0)
2908	{
2909	al = j;
2910	break;
2911	}
2912	}
2913	}
2914	dc2 = dc1 = (BASIS*(int32)a1+(int32)a2) / db1;
2915	loop
2916	{
2917	nlSubMult(ansi, bn, dc1, bl);
2918	a1 = an[al+1];
2919	if (a1!=0)
2920	{
2921	dc1 = (BASIS*(int32)a1+(int32)an[al]) / db1;
2922	dc2 += dc1;
2923	}
2924	else
2925	break;
2926	}
2927	a1 = an[al];
2928	a2 = an[al-1];
2929	rel = nlCompHigh(ansi, bn, bl);
2930	if (rel)
2931	{
2932	dc = (BASIS*(int32)a1 + (int32)a2) / db;
2933	if (dc > 1)
2934	{
2935	nlSubMult(ansi, bn, dc, bl);
2936	rel = nlCompHigh(ansi, bn, bl);
2937	}
2938	else
2939	dc = 0;
2940	while (rel)
2941	{
2942	dc++;
2943	nlSubSimple(ansi, bn, bl);
2944	rel = nlCompHigh(ansi, bn, bl);
2945	}
2946	a1 = an[al];
2947	a2 = an[al-1];
2948	}
2949	else
2950	dc = 0;
2951	cc[si+1] = (int16) (dc2+dc);
2952	}
2953	}
2954
2955	/*3
2956	* aq := aq div b
2957	*/
2958	static void nlQuotient (lint * aq, lint b)
2959	{
2960	lint a;
2961	int bl;
2962	BOOLEAN neg;
2963	a = *aq;
2964	neg = (NLSIGN(a) != NLSIGN(b));
2965	bl = NLLENGTH(b);
2966	if(bl == 1)
2967	nlQuotSmall(&a,b[1]);
2968	else
2969	nlQuotLong(&a,b);
2970	if(neg && !NLSIGN(a))
2971	*a += SIGNUM;
2972	else if (!neg && NLSIGN(a))
2973	*a -= SIGNUM;
2974	*aq = a;
2975	}
2976
2977	/*3
2978	* r := a mod b, b small
2979	* al = length(a)
2980	*/
2981	static void nlRemSmall (lint a, int al, int16 b, int16 * r)
2982	{
2983	int j;
2984	int16 a1,c1;
2985	int32 db,da,dc;
2986	db = (int32)b;
2987	a1 = a[al];
2988	if (al == 1)
2989	{
2990	c1 = a1 / b;
2991	a1 -= b*c1;
2992	*r = a1;
2993	return;
2994	}
2995	if (a1 >= b)
2996	{
2997	c1 = a1 / b;
2998	a1 -= b*c1;
2999	}
3000	loop
3001	{
3002	al--;
3003	if (al==0)
3004	{
3005	*r = a1;
3006	return;
3007	}
3008	if (a1==0)
3009	{
3010	j = al;
3011	loop
3012	{
3013	a1 = a[j];
3014	if (a1 >= b)
3015	{
3016	c1 = a1 / b;
3017	a1 -= b*c1;
3018	}
3019	j--;
3020	if (j==0)
3021	{
3022	*r = a1;
3023	return;
3024	}
3025	else if (a1!=0)
3026	{
3027	al = j;
3028	break;
3029	}
3030	}
3031	}
3032	da = BASIS*(int32)a1 + (int32)a[al];
3033	dc = da / db;
3034	da -= db*dc;
3035	a1 = (int16)da;
3036	}
3037	}
3038
3039	static void nlShrink0(lint c, int j, int * i)
3040	{
3041	while((j!=0) && (c[j]==0)) j--;
3042	*i = j;
3043	}
3044
3045	/*3
3046	* rl := an mod bn, bn long, but <= an
3047	* al = length(an)
3048	*/
3049	static void nlRemLong(lint an, lint bn, int al, int bl, int * rl)
3050	{
3051	int si,j;
3052	int16 b1,a1,a2;
3053	int32 db,db1,dc;
3054	lint ansi;
3055	BOOLEAN rel;
3056	b1 = bn[bl];
3057	db = (int32)bn[bl-1];
3058	db += (int32)b1<<SHIFT_BASIS;
3059	if ((bl > 2) && ((b1 < MAX_INT16) \|\| (bn[bl-1] < MAX_INT16)))
3060	db++;
3061	db1 = (int32)b1+1;
3062	a1 = an[al];
3063	a2 = an[al-1];
3064	si = al-bl;
3065	ansi = an+si;
3066	rel = nlCompHigh(ansi, bn, bl);
3067	if (rel)
3068	{
3069	dc = (BASIS*(int32)a1 + (int32)a2) / db;
3070	if (dc > 1)
3071	{
3072	nlSubMult(ansi, bn, dc, bl);
3073	rel = nlCompHigh(ansi, bn, bl);
3074	}
3075	while (rel)
3076	{
3077	nlSubSimple(ansi, bn, bl);
3078	rel = nlCompHigh(ansi, bn, bl);
3079	}
3080	a1 = an[al];
3081	a2 = an[al-1];
3082	}
3083	loop
3084	{
3085	al--;
3086	if (al < bl)
3087	{
3088	nlShrink0(an,bl,rl);
3089	return;
3090	}
3091	si--;
3092	ansi--;
3093	if(a1==0)
3094	{
3095	j = al;
3096	loop
3097	{
3098	a1 = a2;
3099	a2 = an[j-1];
3100	rel = nlCompHigh(ansi, bn, bl);
3101	if (rel)
3102	{
3103	dc = (BASIS*(int32)a1 + (int32)a2) / db;
3104	if (dc > 1)
3105	{
3106	nlSubMult(ansi, bn, dc, bl);
3107	rel = nlCompHigh(ansi, bn, bl);
3108	}
3109	while (rel)
3110	{
3111	nlSubSimple(ansi, bn, bl);
3112	rel = nlCompHigh(ansi, bn, bl);
3113	}
3114	a1 = an[j];
3115	a2 = an[j-1];
3116	}
3117	si--;
3118	ansi--;
3119	j--;
3120	if (j < bl)
3121	{
3122	nlShrink0(an,bl,rl);
3123	return;
3124	}
3125	else if (a1!=0)
3126	{
3127	al = j;
3128	break;
3129	}
3130	}
3131	}
3132	dc = (BASIS*(int32)a1+(int32)a2) / db1;
3133	loop
3134	{
3135	nlSubMult(ansi, bn, dc, bl);
3136	a1 = an[al+1];
3137	if (a1!=0)
3138	dc = (BASIS*(int32)a1+(int32)an[al]) / db1;
3139	else
3140	break;
3141	}
3142	a1 = an[al];
3143	a2 = an[al-1];
3144	rel = nlCompHigh(ansi, bn, bl);
3145	if (rel)
3146	{
3147	dc = (BASIS*(int32)a1 + (int32)a2) / db;
3148	if (dc > 1)
3149	{
3150	nlSubMult(ansi, bn, dc, bl);
3151	rel = nlCompHigh(ansi, bn, bl);
3152	}
3153	while (rel)
3154	{
3155	nlSubSimple(ansi, bn, bl);
3156	rel = nlCompHigh(ansi, bn, bl);
3157	}
3158	a1 = an[al];
3159	a2 = an[al-1];
3160	}
3161	}
3162	}
3163
3164	/*3
3165	* z = GCD(a, b)
3166	*/
3167	static void nlGcdLong (lint a, lint b, lint * z)
3168	{
3169	lint h, x, y;
3170	int xl,yl,hl;
3171	int16 xc,yc,hc;
3172
3173	x = nlCopyLint(a);
3174	if (NLSIGN(x))
3175	*x %= SIGNUM;
3176	if (nlIsOneLint(x))
3177	{
3178	*z=x;
3179	return;
3180	}
3181	y = nlCopyLint(b);
3182	if (NLSIGN(y))
3183	*y %= SIGNUM;
3184	if (nlIsOneLint(y))
3185	{
3186	*z=y;
3187	nlDeleteLint(x);
3188	return;
3189	}
3190	switch (nlCompAbs(x,y))
3191	{
3192	case -1: h=x;x=y;y=h;h=NULL;
3193	break;
3194	case 0: *z=x;
3195	nlDeleteLint(y);
3196	return;
3197	case 1: break;
3198	}
3199	xl = NLLENGTH(x);
3200	yl = NLLENGTH(y);
3201	/* we have: x > y >=0 as lint, xl>=yl>=0 as int*/
3202	while (yl > 1)
3203	{
3204	nlRemLong(x,y,xl,yl,&hl);
3205	h = x;
3206	x = y;
3207	y = h;
3208	xl = yl;
3209	yl = hl;
3210	}
3211	if(yl==0)
3212	{
3213	/result is x with xl >1/
3214	nlDeleteLint(y);
3215	nlShrink(&x);
3216	*z = x;
3217	return;
3218	}
3219	else
3220	{
3221	yc = y[1];
3222	if (xl>1)
3223	nlRemSmall(x,xl,yc,&hc);
3224	else
3225	{
3226	hc = x[1];
3227	if (hc > yc)
3228	{
3229	xc=yc;
3230	yc=hc;
3231	hc=xc;
3232	}
3233	}
3234	}
3235	nlDeleteLint(x);
3236	nlDeleteLint(y);
3237	/we have now yc >hc/
3238	while(hc!=0)
3239	{
3240	xc = yc;
3241	yc = hc;
3242	hc = xc % yc;
3243	}
3244	*z = nlToLint(yc);
3245	}
3246
3247	/-----------------------------------------------------------------/
3248
3249	/*0
3250	* Begin of operations with rational numbers
3251	*/
3252
3253	/*3
3254	* parameter s in number:
3255	* 0 : integer
3256	* 1 : normalised rational
3257	* -1 : not normalised rational
3258	*/
3259
3260	#define NLISINT(a) !((a)->s)
3261	#define NLISRAT(a) ((a)->s)
3262
3263	#ifdef LDEBUG
3264	#define nlTest(a) nlDBTest(a,__FILE__,__LINE__)
3265	BOOLEAN nlDBTest(number a, char *f,int l)
3266	{
3267	BOOLEAN res;
3268	if (a==NULL)
3269	{
3270	return TRUE;
3271	}
3272	else
3273	{
3274	if (a->z!=NULL)
3275	{
3276	res=mmDBTestBlock(a->z,(NLLENGTH(a->z)+1)*sizeof(int16),f,l);
3277	}
3278	if (NLISRAT(a))
3279	{
3280	res = (res && mmDBTestBlock(a->n,(NLLENGTH(a->n)+1)*sizeof(int16),f,l));
3281	}
3282	return res;
3283	}
3284	}
3285	#endif
3286
3287	/*2
3288	* delete a
3289	*/
3290	#ifdef LDEBUG
3291	void nlDBDelete (number a, char f, int l)
3292	#else
3293	void nlDelete (number *a)
3294	#endif
3295	{
3296	number b=*a;
3297	#ifdef LDEBUG
3298	nlTest(b);
3299	#endif
3300
3301	if (b!=NULL)
3302	{
3303	nlDeleteLint(b->z);
3304	if (NLISRAT(b))
3305	{
3306	nlDeleteLint(b->n);
3307	}
3308	#if defined(LDEBUG) && defined(MDEBUG)
3309	mmDBFreeBlock((ADDRESS)b,sizeof(rnumber),f,l);
3310	#else
3311	Free((ADDRESS)b,sizeof(rnumber));
3312	#endif
3313	*a=NULL;
3314	}
3315	}
3316
3317	/*2
3318	* r = 0
3319	*/
3320	void nlNew (number * r)
3321	{
3322	*r=NULL;
3323	}
3324
3325	/*2
3326	* result->z =gcd(a->z,b->z)
3327	*/
3328	number nlGcd(number a, number b)
3329	{
3330	number result=(number)Alloc(sizeof(rnumber));
3331
3332	result->s = 0;
3333	nlGcdLong(a->z,b->z,&result->z);
3334	#ifdef LDEBUG
3335	nlTest(result);
3336	#endif
3337	return result;
3338	}
3339
3340	/*2
3341	* result->z = lcm(a->z,b->n)
3342	*/
3343	number nlLcm(number a, number b)
3344	{
3345	lint t,bt;
3346	number result=(number)Alloc(sizeof(rnumber));
3347
3348	result->s = 0;
3349	if (NLISINT(b))
3350	{
3351	result->z = nlCopyLint(a->z);
3352	}
3353	else
3354	{
3355	nlGcdLong(a->z,b->n,&t);
3356	if (!nlIsOneLint(t))
3357	{
3358	bt = nlCopyLint(b->n);
3359	nlQuotient(&bt,t);
3360	result->z = nlMultLong(a->z,bt);
3361	nlDeleteLint(bt);
3362	}
3363	else
3364	{
3365	result->z = nlMultLong(a->z,b->n);
3366	}
3367	nlDeleteLint(t);
3368	}
3369	#ifdef LDEBUG
3370	nlTest(result);
3371	#endif
3372	return result;
3373	}
3374
3375	/*2
3376	* a == 0 ?
3377	*/
3378	BOOLEAN nlIsZero (number a)
3379	{
3380	#ifdef LDEBUG
3381	nlTest(a);
3382	#endif
3383	return (a == NULL);
3384	}
3385
3386	/*2
3387	* a = b ?
3388	*/
3389	BOOLEAN nlEqual (number a, number b)
3390	{
3391	number h;
3392
3393	if(a==NULL)
3394	return (b==NULL);
3395	if(b==NULL)
3396	return(a==NULL);
3397	if(NLISINT(a) && NLISINT(b))
3398	{
3399	return !nlComp(a->z,b->z);
3400	}
3401	h = nlSub(a, b);
3402	if (nlIsZero(h))
3403	{
3404	nlDelete(&h);
3405	return TRUE;
3406	}
3407	else
3408	{
3409	nlDelete(&h);
3410	return FALSE;
3411	}
3412	}
3413
3414	/*2
3415	* a == 1 ?
3416	*/
3417	BOOLEAN nlIsOne (number a)
3418	{
3419	#ifdef LDEBUG
3420	nlTest(a);
3421	#endif
3422	if(nlIsZero(a) \|\| NLSIGN(a->z))
3423	{
3424	return FALSE;
3425	}
3426	if(NLISINT(a))
3427	{
3428	return ((a->z[1]==1) && (a->z[0]==1));
3429	}
3430	if (nlCompAbs(a->z, a->n))
3431	{
3432	return FALSE;
3433	}
3434	nlDeleteLint(a->z);
3435	nlDeleteLint(a->n);
3436	a->z = nlInitOne();
3437	a->s = 0;
3438	return TRUE;
3439	}
3440
3441	/*2
3442	* a == -1 ?
3443	*/
3444	BOOLEAN nlIsMOne (number a)
3445	{
3446	#ifdef LDEBUG
3447	nlTest(a);
3448	#endif
3449	if(nlIsZero(a) \|\| !NLSIGN(a->z))
3450	{
3451	return FALSE;
3452	}
3453	if(NLISINT(a))
3454	{
3455	return ((a->z[1]==1) && (a->z[0]==SIGNUM+1));
3456	}
3457	if (nlCompAbs(a->z, a->n))
3458	{
3459	return FALSE;
3460	}
3461	nlDeleteLint(a->z);
3462	nlDeleteLint(a->n);
3463	a->z = nlInitOne();
3464	a->z[0] += SIGNUM;
3465	a->s = 0;
3466	return TRUE;
3467	}
3468
3469	/*2
3470	* z := i
3471	*/
3472	number nlInit (int i)
3473	{
3474	lint nom;
3475	BOOLEAN neg = FALSE;
3476	number z = NULL;
3477
3478	if (i!=0)
3479	{
3480	z = (number)Alloc(sizeof(rnumber));
3481	z->s = 0;
3482	if (i < 0)
3483	{
3484	neg = TRUE;
3485	i = -i;
3486	}
3487	if (i>=BASIS) // init a 4 byte number
3488	{
3489	nom = nlAllocLint(2);
3490	nom[1] = (int16) (MOD_BASIS&i);
3491	nom[2] = (int16) (i>>SHIFT_BASIS);
3492	}
3493	else // init a 2 byte number
3494	{
3495	nom = nlAllocLint(1);
3496	nom[1] = (int16) i;
3497	}
3498	if (neg)
3499	{
3500	*nom += SIGNUM;
3501	}
3502	z->z = nom;
3503	#ifdef LDEBUG
3504	nlTest(z);
3505	#endif
3506	}
3507	return z;
3508	}
3509
3510	/*2
3511	* convert number to int
3512	*/
3513	int nlInt(number &a)
3514	{
3515	lint x;
3516	int res, i;
3517
3518	if (nlIsZero(a))
3519	{
3520	return 0;
3521	}
3522	nlNormalize(a);
3523	if (NLISRAT(a))
3524	{
3525	return 0;
3526	}
3527	x = a->z;
3528	i = NLLENGTH(x);
3529	if (i > 2)
3530	{
3531	return 0;
3532	}
3533	res = x[1];
3534	if (i == 2)
3535	{
3536	if (x[2] < SIGNUM)
3537	{
3538	res += BASIS*x[2];
3539	}
3540	else
3541	{
3542	return 0;
3543	}
3544	}
3545	if (NLSIGN(x))
3546	{
3547	return -res;
3548	}
3549	else
3550	{
3551	return res;
3552	}
3553	}
3554
3555	/*2
3556	* simplify x
3557	*/
3558	void nlNormalize (number &x)
3559	{
3560	lint zh, nh, d;
3561
3562	if ((x==NULL)\|\|(x->s>=0))
3563	{
3564	return;
3565	}
3566	x->s = 1;
3567	zh = x->z;
3568	nh = x->n;
3569	nlGcdLong(zh, nh, &d);
3570	if (!nlIsOneLint(d))
3571	{
3572	nlQuotient(&zh,d);
3573	x->z = zh;
3574	nlQuotient(&nh,d);
3575	if (nlIsOneLint(nh))
3576	{
3577	nlDeleteLint(nh);
3578	x->s = 0;
3579	}
3580	else
3581	{
3582	x->n = nh;
3583	}
3584	}
3585	nlDeleteLint(d);
3586	#ifdef LDEBUG
3587	nlTest(x);
3588	#endif
3589	}
3590
3591	/*2
3592	* za >= 0 ?
3593	*/
3594	BOOLEAN nlGreaterZero (number a)
3595	{
3596	if (nlIsZero(a))
3597	{
3598	return TRUE;
3599	}
3600	return !NLSIGN(a->z);
3601	}
3602
3603	/*2
3604	* a := -a
3605	*/
3606	number nlNeg (number a)
3607	{
3608	if (!nlIsZero(a))
3609	{
3610	a->z[0] = NLNEG(a->z);
3611	}
3612	return a;
3613	}
3614
3615	/*
3616	* 1/a
3617	*/
3618	number nlInvers(number a)
3619	{
3620	number n=(number)Alloc(sizeof(rnumber));
3621
3622	if (nlIsZero(a))
3623	{
3624	WerrorS("div. 1/0");
3625	Free((ADDRESS)n,sizeof(rnumber));
3626	return NULL;
3627	}
3628	nlNormalize(a);
3629	if(NLISRAT(a))
3630	{
3631	n->z = nlCopyLint(a->n);
3632	}
3633	else
3634	{
3635	n->z = nlInitOne();
3636	}
3637	if((a->z[1]!=1) \|\| (NLLENGTH(a->z)!=1))
3638	{
3639	n->n = nlCopyLint(a->z);
3640	n->s = 1;
3641	if(n->n[0]>SIGNUM)
3642	{
3643	n->n[0] -= SIGNUM;
3644	n->z[0] += SIGNUM;
3645	}
3646	}
3647	else
3648	{
3649	if (NLSIGN(a->z))
3650	{
3651	n->z[0] += SIGNUM;
3652	}
3653	n->s = 0;
3654	}
3655	#ifdef LDEBUG
3656	nlTest(n);
3657	#endif
3658	return n;
3659	}
3660
3661	/*2
3662	* copy a to b
3663	*/
3664	number nlCopy(number a)
3665	{
3666	number b=(number)Alloc0(sizeof(rnumber));
3667
3668	if (!nlIsZero(a))
3669	{
3670	b->s = a->s;
3671	b->z = nlCopyLint(a->z);
3672	if (NLISRAT(a))
3673	{
3674	b->n = nlCopyLint(a->n);
3675	}
3676	#ifdef LDEBUG
3677	nlTest(b);
3678	#endif
3679	return b;
3680	}
3681	else
3682	{
3683	Free((ADDRESS)b,sizeof(rnumber));
3684	return NULL;
3685	}
3686	}
3687
3688	/*2
3689	* lu:= la + li
3690	*/
3691	number nlAdd (number la, number li)
3692	{
3693	number lu;
3694	lint x, y;
3695	BOOLEAN a0 = nlIsZero(la), i0 = nlIsZero(li);
3696
3697	if (a0 \|\| i0)
3698	{
3699	if (a0)
3700	{
3701	return nlCopy(li);
3702	}
3703	else
3704	{
3705	return nlCopy(la);
3706	}
3707	}
3708	lu = (number)Alloc(sizeof(rnumber));
3709	if (NLISINT(la))
3710	{
3711	if (NLISINT(li))
3712	{
3713	lu->z = nlAddLong(la->z, li->z);
3714	if (lu->z==NULL)
3715	{
3716	Free((ADDRESS)lu,sizeof(rnumber));
3717	return NULL;
3718	}
3719	lu->s = 0;
3720	#ifdef LDEBUG
3721	nlTest(lu);
3722	#endif
3723	return lu;
3724	}
3725	else
3726	{
3727	x = nlMultLong(la->z, li->n);
3728	lu->z = nlAddLong(x, li->z);
3729	nlDeleteLint(x);
3730	lu->n = nlCopyLint(li->n);
3731	}
3732	}
3733	else
3734	{
3735	if (NLISINT(li))
3736	{
3737	y = nlMultLong(li->z, la->n);
3738	lu->z = nlAddLong(la->z, y);
3739	nlDeleteLint(y);
3740	lu->n = nlCopyLint(la->n);
3741	}
3742	else
3743	{
3744	x = nlMultLong(la->z, li->n);
3745	y = nlMultLong(li->z, la->n);
3746	lu->z = nlAddLong(x, y);
3747	nlDeleteLint(x);
3748	nlDeleteLint(y);
3749	lu->n = nlMultLong(la->n, li->n);
3750	}
3751	}
3752	if (lu->z==NULL)
3753	{
3754	nlDeleteLint(lu->n);
3755	Free((ADDRESS)lu,sizeof(rnumber));
3756	return NULL;
3757	}
3758	lu->s = -1;
3759	#ifdef LDEBUG
3760	nlTest(lu);
3761	#endif
3762	return lu;
3763	}
3764
3765	/*2
3766	* lu:= la - li
3767	*/
3768	number nlSub (number la, number li)
3769	{
3770	number lu;
3771	lint x, y;
3772	BOOLEAN a0 = nlIsZero(la), i0 = nlIsZero(li);
3773
3774	if (a0 \|\| i0)
3775	{
3776	if (a0)
3777	{
3778	lu = nlCopy(li);
3779	lu = nlNeg(lu);
3780	return lu;
3781	}
3782	else
3783	{
3784	return nlCopy(la);
3785	}
3786	}
3787	lu = (number)Alloc(sizeof(rnumber));
3788	if (NLISINT(la))
3789	{
3790	if (NLISINT(li))
3791	{
3792	lu->z = nlSubLong(la->z, li->z);
3793	if (lu->z==NULL)
3794	{
3795	Free((ADDRESS)lu,sizeof(rnumber));
3796	return NULL;
3797	}
3798	lu->s = 0;
3799	#ifdef LDEBUG
3800	nlTest(lu);
3801	#endif
3802	return lu;
3803	}
3804	else
3805	{
3806	x = nlMultLong(la->z, li->n);
3807	lu->z = nlSubLong(x, li->z);
3808	nlDeleteLint(x);
3809	lu->n = nlCopyLint(li->n);
3810	}
3811	}
3812	else
3813	{
3814	if (NLISINT(li))
3815	{
3816	y = nlMultLong(li->z, la->n);
3817	lu->z = nlSubLong(la->z, y);
3818	nlDeleteLint(y);
3819	lu->n = nlCopyLint(la->n);
3820	}
3821	else
3822	{
3823	x = nlMultLong(la->z, li->n);
3824	y = nlMultLong(li->z, la->n);
3825	lu->z = nlSubLong(x, y);
3826	nlDeleteLint(x);
3827	nlDeleteLint(y);
3828	lu->n = nlMultLong(la->n, li->n);
3829	}
3830	}
3831	if (lu->z==NULL)
3832	{
3833	nlDeleteLint(lu->n);
3834	Free((ADDRESS)lu,sizeof(rnumber));
3835	return NULL;
3836	}
3837	lu->s = -1;
3838	#ifdef LDEBUG
3839	nlTest(lu);
3840	#endif
3841	return lu;
3842	}
3843
3844	/*2
3845	* lo := la * li
3846	*/
3847	number nlMult (number la, number li)
3848	{
3849	number lo=(number)Alloc(sizeof(rnumber));
3850
3851	if (nlIsZero(la) \|\| nlIsZero(li))
3852	{
3853	Free((ADDRESS)lo,sizeof(rnumber));
3854	return NULL;
3855	}
3856	lo->z = nlMultLong(la->z,li->z);
3857	if(NLISINT(la))
3858	{
3859	if(NLISINT(li))
3860	{
3861	lo->s = 0;
3862	#ifdef LDEBUG
3863	nlTest(lo);
3864	#endif
3865	return lo;
3866	}
3867	else
3868	{
3869	lo->n = nlCopyLint(li->n);
3870	}
3871	}
3872	else
3873	{
3874	if(NLISINT(li))
3875	{
3876	lo->n = nlCopyLint(la->n);
3877	}
3878	else
3879	{
3880	lo->n = nlMultLong(la->n,li->n);
3881	}
3882	}
3883	lo->s = -1;
3884	#ifdef LDEBUG
3885	nlTest(lo);
3886	#endif
3887	return lo;
3888	}
3889
3890	/*2
3891	* lo := la / li in Z
3892	* with \|li\| > la - (lo * li) >= 0
3893	*/
3894	number nlIntDiv (number la, number li)
3895	{
3896	int c;
3897	number lo=(number)Alloc(sizeof(rnumber));
3898
3899	if (nlIsZero(la) \|\| nlIsZero(li))
3900	{
3901	if (nlIsZero(li))
3902	{
3903	WerrorS("div. by 0");
3904	}
3905	Free((ADDRESS)lo,sizeof(rnumber));
3906	return NULL;
3907	}
3908	lo->s = 0;
3909	c = nlCompAbs(li->z, la->z);
3910	if (c >= 0)
3911	{
3912	if (c == 0)
3913	{
3914	lo->z = nlInitOne();
3915	if (NLSIGN(la->z) != NLSIGN(li->z))
3916	{
3917	lo->z[0] += SIGNUM;
3918	}
3919	}
3920	else if (NLSIGN(la->z))
3921	{
3922	lo->z = nlInitOne();
3923	if (!NLSIGN(li->z))
3924	{
3925	lo->z[0] += SIGNUM;
3926	}
3927	}
3928	else
3929	{
3930	Free((ADDRESS)lo,sizeof(rnumber));
3931	return NULL;
3932	}
3933	#ifdef LDEBUG
3934	nlTest(lo);
3935	#endif
3936	return lo;
3937	}
3938	lo->z = nlCopyLint(la->z);
3939	if (NLLENGTH(la->z)==1)
3940	{
3941	lo->z[1] /= li->z[1];
3942	if (NLSIGN(la->z))
3943	{
3944	int16 h = la->z[1] - (lo->z[1]*li->z[1]);
3945	if (h!=0)
3946	{
3947	lo->z[1]++;
3948	}
3949	}
3950	if (NLSIGN(li->z))
3951	{
3952	lo->z[0] = NLNEG(lo->z);
3953	}
3954	#ifdef LDEBUG
3955	nlTest(lo);
3956	#endif
3957	return lo;
3958	}
3959	if ((li->z[1]!=1) \|\| (NLLENGTH(li->z)>1))
3960	{
3961	nlQuotient(&lo->z, li->z);
3962	if (NLSIGN(la->z))
3963	{
3964	lint o1, t, h = nlMultLong(lo->z, li->z);
3965	if (nlCompAbs(h, la->z))
3966	{
3967	o1 = nlInitOne();
3968	if (NLSIGN(li->z))
3969	{
3970	t = nlAddLong(lo->z,o1);
3971	}
3972	else
3973	{
3974	t = nlSubLong(lo->z,o1);
3975	}
3976	nlDeleteLint(o1);
3977	nlDeleteLint(lo->z);
3978	lo->z = t;
3979	}
3980	nlDeleteLint(h);
3981	}
3982	}
3983	else if (li->z[0] == (SIGNUM+1))
3984	{
3985	lo->z[0] = NLNEG(lo->z);
3986	}
3987	#ifdef LDEBUG
3988	nlTest(lo);
3989	#endif
3990	return lo;
3991	}
3992
3993	/*2
3994	* lo := la mod li == la - ((la / li) * li)
3995	* with \|li\| > lo >= 0
3996	*/
3997	number nlIntMod (number la, number li)
3998	{
3999	int16 xl, yl;
4000	int c;
4001	lint r, t;
4002	number lo=(number)Alloc(sizeof(rnumber));
4003
4004	if (nlIsZero(la) \|\| nlIsZero(li))
4005	{
4006	if (nlIsZero(li))
4007	{
4008	WerrorS("mod. by 0");
4009	}
4010	Free((ADDRESS)lo,sizeof(rnumber));
4011	return NULL;
4012	}
4013	c = nlCompAbs(li->z, la->z);
4014	if (c == 0)
4015	{
4016	Free((ADDRESS)lo,sizeof(rnumber));
4017	return NULL;
4018	}
4019	if (c > 0)
4020	{
4021	if (NLSIGN(la->z))
4022	{
4023	if (NLSIGN(li->z))
4024	{
4025	lo->z = nlSubLong(la->z,li->z);
4026	}
4027	else
4028	{
4029	lo->z = nlAddLong(la->z,li->z);
4030	}
4031	}
4032	else
4033	{
4034	lo->z = nlCopyLint(la->z);
4035	}
4036	#ifdef LDEBUG
4037	nlTest(lo);
4038	#endif
4039	return lo;
4040	}
4041	r = nlCopyLint(la->z);
4042	xl = NLLENGTH(r);
4043	if (xl == 1)
4044	{
4045	yl = r[1] / li->z[1];
4046	r[1] -= yl * li->z[1];
4047	if (r[1] != 0)
4048	{
4049	if (NLSIGN(la->z))
4050	{
4051	r[1] = li->z[1]-r[1];
4052	}
4053	lo->z = r;
4054	#ifdef LDEBUG
4055	nlTest(lo);
4056	#endif
4057	return lo;
4058	}
4059	else
4060	{
4061	nlDeleteLint(r);
4062	Free((ADDRESS)lo,sizeof(rnumber));
4063	return NULL;
4064	}
4065	}
4066	yl = NLLENGTH(li->z);
4067	if (yl > 1)
4068	{
4069	nlRemLong(r,li->z,xl,yl,&c);
4070	xl = c;
4071	if (xl != 0)
4072	{
4073	nlAllocLint(xl);
4074	for (c=xl; c>0; c--)
4075	{
4076	t[c] = r[c];
4077	}
4078	}
4079	}
4080	else if (li->z[1]!=1)
4081	{
4082	nlRemSmall(r,xl,li->z[1],&xl);
4083	if (xl != 0)
4084	{
4085	t = nlAllocLint(1);
4086	t[1] = xl;
4087	}
4088	}
4089	else
4090	{
4091	xl = 0;
4092	}
4093	nlDeleteLint(r);
4094	if (xl == 0)
4095	{
4096	Free((ADDRESS)lo,sizeof(rnumber));
4097	return NULL;
4098	}
4099	if (NLSIGN(la->z))
4100	{
4101	if (NLSIGN(li->z))
4102	{
4103	r = nlSubLong(t,li->z);
4104	}
4105	else
4106	{
4107	r = nlSubLong(li->z,t);
4108	}
4109	nlDeleteLint(t);
4110	t = r;
4111	}
4112	lo->z = t;
4113	#ifdef LDEBUG
4114	nlTest(lo);
4115	#endif
4116	return lo;
4117	}
4118
4119	/*2
4120	* lo := la / li
4121	*/
4122	number nlDiv (number la, number li)
4123	{
4124	number lo=(number)Alloc(sizeof(rnumber));
4125
4126	if (nlIsZero(la) \|\| nlIsZero(li))
4127	{
4128	if (nlIsZero(li))
4129	{
4130	WerrorS("div. by 0");
4131	}
4132	Free((ADDRESS)lo,sizeof(rnumber));
4133	return NULL;
4134	}
4135	if (NLISINT(la))
4136	{
4137	lo->n = nlCopyLint(li->z);
4138	}
4139	else
4140	{
4141	lo->n = nlMultLong(la->n, li->z);
4142	}
4143	if (NLISINT(li))
4144	{
4145	lo->z = nlCopyLint(la->z);
4146	}
4147	else
4148	{
4149	lo->z = nlMultLong(la->z, li->n);
4150	}
4151	if(NLSIGN(lo->n))
4152	{
4153	lo->n[0] -= SIGNUM;
4154	lo->z[0] = NLNEG(lo->z);
4155	}
4156	lo->s = -1;
4157	if(nlIsOneLint(lo->n))
4158	{
4159	nlDeleteLint(lo->n);
4160	lo->s = 0;
4161	}
4162	#ifdef LDEBUG
4163	nlTest(lo);
4164	#endif
4165	return lo;
4166	}
4167
4168	/*3
4169	* Power of long; lu:= x^exp
4170	*/
4171	static void nlPow (lint x,int exp,lint *lu)
4172	{
4173	lint yh, y, w;
4174	y = nlInitOne();
4175	w = nlCopyLint(x);
4176	while (exp)
4177	{
4178	if (exp%2 != 0)
4179	{
4180	yh = nlMultLong(y,w);
4181	nlDeleteLint(y);
4182	y = yh;
4183	}
4184	if (exp > 1)
4185	{
4186	yh = nlMultLong(w,w);
4187	nlDeleteLint(w);
4188	w = yh;
4189	}
4190	exp = exp / 2;
4191	}
4192	*lu = y;
4193	}
4194
4195	/*2
4196	* lu:= x ^ exp
4197	*/
4198	void nlPower (number x,int exp,number * lu)
4199	{
4200	number b=NULL;
4201
4202	if (!nlIsZero(x))
4203	{
4204	b = (number)Alloc(sizeof(rnumber));
4205	if (exp == 0)
4206	{
4207	b->z = nlInitOne();
4208	*lu = b;
4209	#ifdef LDEBUG
4210	nlTest(*lu);
4211	#endif
4212	return;
4213	}
4214	nlNormalize(x);
4215	if (exp > 0)
4216	{
4217	b->s = x->s;
4218	if ((x->z[1]!=1) \|\| (NLLENGTH(x->z)>1))
4219	{
4220	nlPow(x->z,exp,&b->z);
4221	}
4222	else
4223	{
4224	b->z = nlInitOne();
4225	}
4226	if (NLISRAT(x))
4227	{
4228	nlPow(x->n,exp,&b->n);
4229	}
4230	}
4231	else
4232	{
4233	if ((x->z[1]!=1) \|\| (NLLENGTH(x->z)>1))
4234	{
4235	b->s = 1;
4236	nlPow(x->z,exp,&b->n);
4237	if (NLSIGN(b->n))
4238	{
4239	b->n[0] -= SIGNUM;
4240	}
4241	}
4242	else
4243	{
4244	b->s = 0;
4245	}
4246	if (NLISRAT(x))
4247	{
4248	nlPow(x->n,exp,&b->z);
4249	}
4250	else
4251	{
4252	b->z = nlInitOne();
4253	}
4254	}
4255	if (exp%2 != 0) // make sign
4256	{
4257	if(NLSIGN(b->z) != NLSIGN(x->z))
4258	{
4259	b->z[0] = NLNEG(b->z);
4260	}
4261	}
4262	else
4263	{
4264	if(NLSIGN(b->z))
4265	{
4266	b->z[0] -= SIGNUM;
4267	}
4268	}
4269	}
4270	*lu = b;
4271	#ifdef LDEBUG
4272	nlTest(*lu);
4273	#endif
4274	}
4275
4276	/*2
4277	* a > b ?
4278	*/
4279	BOOLEAN nlGreater (number a, number b)
4280	{
4281	number r;
4282	BOOLEAN rr;
4283
4284	r=nlSub(a,b);
4285	rr=(!nlIsZero(r)) && (nlGreaterZero(r));
4286	nlDelete(&r);
4287	return rr;
4288	}
4289
4290	/*3
4291	* a := b * a + r
4292	*/
4293	void nlDivMod(lint a, int16 al, int16 b, int16 * r)
4294	{
4295	int32 da;
4296	const int32 db = (int32)b;
4297	int16 a1 = a[al];
4298
4299	a[al] /= b;
4300	a1 -= b*a[al];
4301	loop
4302	{
4303	al--;
4304	if (al==0)
4305	{
4306	*r = a1;
4307	return;
4308	}
4309	while (a1==0)
4310	{
4311	a1 = a[al];
4312	a[al] /= b;
4313	a1 -= b*a[al];
4314	al--;
4315	if (al==0)
4316	{
4317	*r = a1;
4318	return;
4319	}
4320	}
4321	da = BASIS*(int32) a1 + (int32) a[al];
4322	a[al] = da / db;
4323	da -= db*a[al];
4324	a1 = (int16)da;
4325	}
4326	}
4327
4328	/*2
4329	* returns (n->z mod p) / (n->n mod p)
4330	*/
4331	int nlModP(number n, int p)
4332	{
4333	number a;
4334	int16 iz, in;
4335	int jz, jn, al;
4336
4337	if (nlIsZero(n))
4338	{
4339	return 0;
4340	}
4341	a = nCopy(n);
4342	al = NLLENGTH(a->z);
4343	if (al == 1)
4344	{
4345	iz = (a->z[1])%p;
4346	}
4347	else
4348	{
4349	nlRemSmall(a->z,al,p,&iz);
4350	}
4351	jz = iz;
4352	if (NLSIGN(a->z))
4353	{
4354	jz = p-jz;
4355	}
4356	if (NLISRAT(a))
4357	{
4358	al = *(a->n);
4359	if (al == 1)
4360	{
4361	in = (a->n[1])%p;
4362	}
4363	else
4364	{
4365	nlRemSmall(a->n,al,p,&in);
4366	}
4367	jn = in;
4368	jz = (int)npDiv((number)jz,(number)jn);
4369	}
4370	nlDelete(&a);
4371	return jz;
4372	}
4373
4374	int nlSize(number a)
4375	{
4376	/to be implemtented/ return 0;
4377	}
4378	#endif
4379

Note: See TracBrowser for help on using the repository browser.

Download in other formats: