Context Navigation

source: git/Singular/longrat.cc @ ea5b50

spielwiese

Last change on this file since ea5b50 was ea5b50, checked in by Hans Schönemann <hannes@…>, 27 years ago
* hannes: fix HAVE_LIBGMP1/2 git-svn-id: file:///usr/local/Singular/svn/trunk@248 2c84dea3-7e68-4137-9b89-c4e89433aadc
Property mode set to `100644`
File size: 73.8 KB

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

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

Download in other formats: