Context Navigation

source: git/Singular/longrat.cc @ ea925e

spielwiese

Last change on this file since ea925e was ea925e, checked in by Hans Schönemann <hannes@…>, 27 years ago
* hannes: removed ECHO/ECHOPOST from Makefile changed constant to 1<<28 in longrat.cc git-svn-id: file:///usr/local/Singular/svn/trunk@184 2c84dea3-7e68-4137-9b89-c4e89433aadc
Property mode set to `100644`
File size: 73.6 KB

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

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

Download in other formats: