Context Navigation

source: git/Singular/longrat.cc @ 6f2edc

spielwiese

Last change on this file since 6f2edc was 40edb03, checked in by Olaf Bachmann <obachman@…>, 27 years ago
Fri Apr 25 16:59:31 1997 Olaf Bachmann <obachman@ratchwum.mathematik.uni-kl.de (Olaf Bachmann)> * Various changes to reflect new configure (versions defined in configure.in, changed HAVE_LIBFACTORY into HAVE_FACTORY, data dir is pasted from configure into mod2.h and used from therein feFopen. * Added configure facility, repalced mod2.h by mod2.h.in Makefile by Makefile.in git-svn-id: file:///usr/local/Singular/svn/trunk@189 2c84dea3-7e68-4137-9b89-c4e89433aadc
Property mode set to `100644`
File size: 73.7 KB

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

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

Download in other formats: