1 /* SYMMETRICA nb.c */
2 /* 4.11.91: TPMcD */
3
4 #include <math.h>
5 /* #include <ctype.h> AK 200206 */ /* for isspace AK 160192 */
6
7
8
9
10 #include "def.h"
11 #include "macro.h"
12
myisspace(int i)13 int myisspace (int i)
14 {
15 if (i=='\t') return 1;
16 if (i=='\n') return 1;
17 if (i=='\r') return 1;
18 if (i=='\v') return 1;
19 if (i=='\f') return 1;
20 if (i==' ') return 1;
21 return 0;
22 }
23
24
25 #define INIT_CYCLO(a) \
26 do {\
27 erg += b_ksd_n(CYCLOTOMIC, CALLOCOBJECT(), NULL, a);\
28 } while(0)
29
30
31 #define nb_quores(a,b,c,d) quores(a,b,c,d) /* AK 050393 */
32
33 static INT space_saving = TRUE;
34 static INT basis_type = STD_BASIS;
35
36 /* STATIC VARIABLES RELATING TO THE MAINTENANCE OF CYCLOTOMIC DATA */
37
38 #ifdef CYCLOTRUE
39
40 /* cyclo_table points to an array of CYCLO_DATA with */
41 /* no_cyclos cyclos. cyclo_table_set is a flag which */
42 /* indicates whether the table is present or not. */
43
44 static INT cyclo_table_set = 0L, cyclo_list_set = 0L;
45 static INT zzno_cyclos;
46 static CYCLO_DATA *zzcyclo_table;
47 static OP zzcyclo_list = NULL, zzcyclo_tail = NULL;
48 static INT number_mem;
49
50 #endif /* CYCLOTRUE */
51
52 static INT setup_prime_table();
53 static INT integer_factor_0();
54 static INT integer_factor_1();
55 static INT insert_zero_into_monopoly();
56 static INT integer_coefficients();
57 static INT find_sqrad_data();
58 static INT adjust_sqrad_data();
59 static INT fprint_sqrad();
60 static INT fprint_cyclo();
61 static INT standardise_cyclo_0();
62 static INT make_index_monopoly_cyclo();
63 static INT add_cyclo_cyclo_0();
64 static INT mult_cyclo_cyclo_0();
65 static INT invers_cyclo_norm();
66 static INT SCMPCO();
67
68 # ifdef CYCLOTRUE
69
70 static INT setup_cyclotomic_table();
71 static CYCLO_DATA *cyclo_ptr();
72 static CYCLO_DATA *add_cyclo_data();
73 static INT free_cyclo_list();
74 static INT free_cyclo_table();
75
76 # endif
s_n_s(a)77 OP s_n_s(a) OP a;
78 /* AK 080891 V1.3 */
79 {
80 if (a == NULL)
81 {
82 error("s_n_s:a == NULL");
83 return (OP) NULL;
84 }
85 return (((a)->ob_self).ob_number)->n_self;
86 }
87
c_n_s(a,b)88 INT c_n_s(a,b) OP a,b;
89 /* AK 200891 V1.3 */
90 {
91 ((a->ob_self).ob_number)->n_self = b;
92 return OK;
93 }
94
s_n_d(a)95 OP s_n_d(a) OP a;
96 /* AK 200891 V1.3 */
97 {
98 if (a == NULL)
99 {
100 error("s_n_d:a == NULL");
101 return (OP) NULL;
102 }
103 return (((((a)->ob_self).ob_number)->n_data).o_data);
104 }
105
c_n_d(a,b)106 INT c_n_d(a,b) OP a,b;
107 /* AK 200891 V1.3 */
108 {
109 ((((((a)->ob_self).ob_number)->n_data).o_data) = (b));
110 return OK;
111 }
112
s_n_dci(a)113 OP s_n_dci(a) OP a;
114 /* AK 200891 V1.3 */
115 {
116 return ((((((a)->ob_self).ob_number)->n_data).c_data)->index);
117 }
118
s_n_dcd(a)119 OP s_n_dcd(a) OP a;
120 /* AK 200891 V1.3 */
121 {
122 return ((((((a)->ob_self).ob_number)->n_data).c_data)->deg);
123 }
124
s_n_dcp(a)125 OP s_n_dcp(a)OP a;
126 /* AK 200891 V1.3 */
127 {
128 return ((((((a)->ob_self).ob_number)->n_data).c_data)->poly);
129 }
130
131 /****************** factors.c **********************/
132 /* 26.07.91: TPMcD. */
133 /*************************************************************/
134 /* Determines and returns the number of digits of the */
135 /* integer a. a may be an INTEGER or a LONGINT. */
136
number_of_digits(a)137 INT number_of_digits(a) OP a;
138 /* 04.04.90: TPMcD */
139 /* AK 200891 V1.3 */
140 {
141 INT i = 1L;
142 OP b = CALLOCOBJECT();
143 OP ten = CALLOCOBJECT();
144
145 M_I_I(10L,ten);
146 copy(a,b);
147 if (LT(b,cons_null) == TRUE)
148 addinvers_apply(b);
149
150 while (GE(b,ten) == TRUE)
151 {
152 ganzdiv(b,ten,b);
153 i++;
154 }
155 freeall(b);
156 freeall(ten);
157 return(i);
158 }
159
number_of_bits(a)160 INT number_of_bits(a) OP a;
161 /* 04.04.90: TPMcD */
162 /* AK 200891 V1.3 */
163 /* AK 300102 */
164 /* AK 180902 V2.1 */
165 {
166 INT erg = OK;
167 CTTO(INTEGER,LONGINT,"number_of_bits(1)",a);
168 {
169 INT i = 1L;
170 OP b,ten;
171
172 if (S_O_K(a) == INTEGER) /* AK 280705 */
173 {
174 INT ai = S_I_I(a);
175 while (ai >= 2) { ai /= 2; i++; }
176 return i;
177 }
178
179 b = CALLOCOBJECT();
180
181 COPY(a,b);
182 if (LT(b,cons_null) == TRUE)
183 ADDINVERS_APPLY(b);
184
185 while (GE(b,cons_zwei) == TRUE)
186 {
187 ganzdiv(b,cons_zwei,b);
188 i++;
189 }
190 freeall(b);
191 return(i);
192 }
193 ENDR("number_of_bits");
194 }
195
196
integer_factors_to_integer(l,a)197 INT integer_factors_to_integer(l,a) OP l,a;
198 /* 10.05.90: TPMcD */
199 /* AK 200891 V1.3 */
200 /* 01.10.91: TPMcD */
201 {
202 INT ret = ERROR;
203 #ifdef MONOPOLYTRUE
204 OP b = CALLOCOBJECT();
205 OP ptr;
206
207 if (S_O_K(l) != MONOPOLY)
208 goto exit_label;
209 if (not EMPTYP(a))
210 freeself(a);
211 M_I_I(1L,a);
212 ptr = l;
213 if (EMPTYP(S_PO_S(ptr)))
214 ptr = S_L_N(ptr); /* skip the empty term */
215 while (ptr != NULL)
216 {
217 hoch(S_PO_S(ptr),S_PO_K(ptr),b);
218 mult(a,b,a);
219 ptr = S_L_N(ptr);
220 }
221 ret = OK;
222 exit_label:
223 freeall(b);
224 #else
225 error("integer_factors_to_integer: MONOPOLY not available");
226 #endif
227 return(ret);
228 }
229
230 /*Given the number n, which should be an positive INTEGER or LONGINT */
231 /*or a MONOPOLY representing a factorisation of an integer greater */
232 /*than 1 , the result returns the list of positive integers coprime to n.*/
233
make_coprimes(number,result)234 INT make_coprimes(number,result) OP number, result;
235 /* 01.05.91: TPMcD */ /* AK 200891 V1.3 */ /* 01.10.91: TPMcD */
236 {
237 INT end_flag = 0L, flag= -1L; /* AK 040292 */
238 INT erg = ERROR;
239 OP ptr, ptr_zwei, ptr_drei, num=NULL;
240 #ifdef MONOPOLYTRUE
241 OP new, list=NULL;
242 OP vec,prime,count_eins,count_zwei;
243
244 CTTTO(INTEGER,MONOPOLY,LONGINT,"make_coprimes(1)",number);
245
246 vec = CALLOCOBJECT();
247 prime = CALLOCOBJECT();
248 count_eins = CALLOCOBJECT();
249 count_zwei = CALLOCOBJECT();
250
251 init(LIST,result);
252 if (S_O_K(number) == MONOPOLY)
253 {
254 list = number;
255 num = CALLOCOBJECT();
256 flag = 1L; /* remember to free num */
257 integer_factors_to_integer(list,num);
258 }
259 else
260 {
261 if (((S_O_K(number) != INTEGER) && (S_O_K(number) != LONGINT))
262 || (LT(number,cons_eins) == TRUE))
263 goto exit_label;
264 if (EQ(number,cons_eins) == TRUE)
265 {
266 new = CALLOCOBJECT();
267 M_I_I((INT)1,new);
268 insert(new,result,NULL,NULL);
269 erg = OK;
270 goto exit_label;
271 }
272 num = number;
273 list = CALLOCOBJECT();
274 flag = 0L; /* remember to free list */
275 integer_factor(num,list);
276 }
277 m_i_i((INT)1,count_eins);
278 init(LIST,vec);
279 ptr = vec;
280 while (TRUE)
281 { /* vec is initialised to the list of numbers 1 , . . ., num */
282 S_L_S(ptr) = CALLOCOBJECT();
283 copy(count_eins,S_L_S(ptr));
284 if (LT(count_eins,num) == TRUE)
285 {
286 new = CALLOCOBJECT();
287 S_L_N(ptr) = new;
288 ptr = new;
289 init(LIST,new);
290 INC(count_eins);
291 }
292 else
293 {
294 S_L_N(ptr) = NULL;
295 break;
296 }
297 }
298 ptr = list;
299 while (ptr != NULL)
300 {
301 copy(S_PO_S(ptr),prime);
302 /*
303 copy(cons_eins,count_eins);
304 copy(cons_eins,count_zwei);
305 */
306 M_I_I(1,count_eins);
307 M_I_I(1,count_zwei);
308 ptr_zwei = vec;
309 while (LE(count_eins,num) == TRUE)
310 { /* delete all multiples of prime from vec */
311 if (EQ(count_zwei,prime) == TRUE)
312 {
313 /*
314 copy(cons_eins,count_zwei);
315 */
316 M_I_I(1,count_zwei);
317 if (not EMPTYP(S_L_S(ptr_zwei))) /* AK */
318 FREESELF(S_L_S(ptr_zwei));
319 }
320 else
321 INC_INTEGER(count_zwei);
322 INC_INTEGER(count_eins);
323 ptr_zwei = S_L_N(ptr_zwei);
324 }
325 ptr = S_L_N(ptr);
326 }
327 ptr = result;
328 ptr_drei = result;
329 copy(cons_eins,count_eins);
330 ptr_zwei = vec;
331 while (TRUE)
332 {
333 if (EQ(count_eins,num) == TRUE)
334 end_flag = 1L;
335 if (not EMPTYP(S_L_S(ptr_zwei)))
336 {
337 S_L_S(ptr) = CALLOCOBJECT();
338 copy(count_eins,S_L_S(ptr));
339 if (end_flag)
340 {
341 S_L_N(ptr) = NULL;
342 break;
343 }
344 else
345 {
346 new = CALLOCOBJECT();
347 init(LIST,new);
348 S_L_N(ptr) = new;
349 ptr_drei = ptr;
350 ptr = new;
351 }
352 }
353 if (end_flag)
354 {
355 freeall(ptr);
356 S_L_N(ptr_drei) = NULL;
357 break;
358 }
359 INC(count_eins);
360 ptr_zwei = S_L_N(ptr_zwei);
361 }
362 erg = OK;
363 exit_label:
364 FREEALL4(vec,prime,count_eins,count_zwei);
365 if (flag == 1L) freeall(num);
366 else if (flag != -1L) freeall(list); /* AK 040292 */
367 #endif
368 ENDR("make_coprimes");
369 }
370
euler_phi(a,b)371 INT euler_phi(a,b) OP a,b;
372 /* AK number of numbers coprime to a */
373 /* AK 310191 V1.2 */ /* AK 200891 V1.3 */
374 {
375 OP c = CALLOCOBJECT();
376 INT erg;
377 erg = make_coprimes(a,c);
378 erg += length(c,b);
379 erg += freeall(c);
380 return erg;
381 }
382
prime_power_p(a)383 INT prime_power_p(a) OP a;
384 /* AK 290304 true if power of a prime */
385 /* AK 161204 V3.0 */
386 {
387 INT erg = OK;
388 CTTO(INTEGER,LONGINT,"prime_power_p(1)",a);
389 if (NULLP(a)) return FALSE;
390 if (NEGP(a)) return FALSE;
391
392 {
393 OP b; INT res;
394 b = CALLOCOBJECT();
395 factorize(a,b);
396 if (EQ(S_V_I(b,0),S_V_I(b,S_V_LI(b)-1))) res = TRUE;
397 else res = FALSE;
398 FREEALL(b);
399 return res;
400 }
401
402 ENDR("prime_power_p");
403 }
404
primep(a)405 INT primep(a) OP a;
406 /* AK 220294 true if prime */
407 /* AK 161204 V3.0 */
408 {
409 OP c,d,e;
410 INT erg = TRUE;
411
412 if (EQ(a,cons_zwei))
413 {
414 erg = TRUE;
415 goto p1;
416 }
417 if (negp(a) || NULLP(a) || EVEN(a) )
418 {
419 erg = FALSE;
420 goto p1;
421 }
422 CALLOCOBJECT3(c,d,e);
423
424 ganzsquareroot(a,c);
425 M_I_I(3,d);
426 while (le(d,c))
427 {
428 mod(a,d,e);
429 if (NULLP(e) )
430 { erg = FALSE; break; }
431 ADD_APPLY(cons_zwei,d);
432 }
433 FREEALL3(c,d,e);
434 p1:
435 return erg;
436 ENDR("primep");
437 }
438
439 /* INTEGER SQUARE ROOTS */
440
441 /* If a is a non-negative integer (INTEGER or LONGINT) */
442 /* s is set to the integer part of its square root. */
443 /* In this case, the return value is OK or IMPROPER */
444 /* according as the integer is a perfect square or not. */
445 /* Otherwise, the return value is ERROR. */
446
nb_square_root(a,s)447 static INT nb_square_root(a,s) OP a,s;
448 /* 04.04.90: TPMcD */ /* AK 200891 V1.3 */
449 {
450 INT a_eins,d_eins,e_eins,ret = ERROR;
451 INT erg = OK;
452 OP b,c,d,e,diff;
453 CTTO(INTEGER,LONGINT,"nb_square_root(1)",a);
454
455 b = CALLOCOBJECT();
456 c = CALLOCOBJECT();
457 d = CALLOCOBJECT();
458 e = CALLOCOBJECT();
459 diff = CALLOCOBJECT();
460
461
462
463 if (negp(a)) /* a < 0L */
464 {
465 fprintf(stderr,"Negative number has no real square root\n");
466 goto exit_label;
467 }
468 if (NULLP(a))
469 {
470 m_i_i(0l,s);
471 ret = OK;
472 goto exit_label;
473 }
474 d_eins = number_of_digits(a);
475
476 e_eins = (d_eins + 1L) / 2L;
477
478 M_I_I(10L,d);
479 M_I_I(e_eins-1L,b);
480 hoch(d,b,b);
481 MULT(d,b,c);
482 FREESELF(d); MULT(b,b,d);
483 if (EQ(a,d) == TRUE)
484 {
485 copy(b,s);
486 ret = OK;
487 goto exit_label;
488 }
489 do
490 {
491 FREESELF(d);
492 ADD(b,c,d);
493 if (negp(d))
494 error("square_root : negative integer unexpectedly encountered\n");
495
496 half_apply(d);
497 FREESELF(diff);
498 ADDINVERS(b,diff);
499 ADD_APPLY(c,diff);
500 FREESELF(e); MULT(d,d,e);
501 a_eins = COMP(a,e);
502 if (a_eins < 0L)
503 copy(d,c);
504 else if (a_eins > 0L)
505 copy(d,b);
506 else
507 {
508 copy(d,s);
509 ret = OK;
510 goto exit_label;
511 }
512
513 } while (GE(diff,cons_zwei) == TRUE);
514 copy(b,s);
515 ret = IMPROPER;
516 exit_label:
517 FREEALL5(b,c,d,e,diff);
518 CTTO(INTEGER,LONGINT,"nb_square_root(e2)",s);
519 return(ret);
520 ENDR("nb_square_root");
521 }
522
523 #ifdef LONGINTTRUE
ganzsquareroot_longint(a,b)524 INT ganzsquareroot_longint(a,b) OP a,b;
525 /* AK 040291 */ /* AK 200891 V1.3 */
526 {
527 INT erg ;
528 erg = nb_square_root(a,b);
529 return (erg == IMPROPER ? OK : erg); /* AK 200194 */
530 }
531 #endif /* LONGINTTRUE */
532
ganzsquareroot_integer(a,b)533 INT ganzsquareroot_integer(a,b) OP a,b;
534 /* AK 040291 */ /* AK 200891 V1.3 */
535 {
536 INT erg ;
537 erg = nb_square_root(a,b);
538 return (erg == IMPROPER ? OK : erg); /* AK 200194 */
539 }
540
541 /* INTEGER FACTORISATION */
542
543 /* Routines for prime factorization of integers. */
544 /* prime_table points to an array of INT with the first no_primes primes.*/
545 /* prime_table_set is a flag which indicates whether the table is present*/
546 /* or not. */
547
548 static INT prime_table_set = 0L, no_primes;
549 static INT *prime_table;
550
551 /*Reads the table of primes from the file PRIMES.DAT. The first entry */
552 /*should be no_primes, then the list of primes. Assumes that INT means */
553 /*long int. Returns OK if the table is set; otherwise, returns ERROR. */
554
setup_prime_table()555 static INT setup_prime_table()
556 /* 040490 TPMcD */
557 /* AK 200891 V1.3 */
558 /* 29.10.91 TPMcD */
559 {
560 # ifdef PRIMES_FILE
561
562 FILE *f;
563 if ( (f = fopen("PRIMES.DAT","r")) == NULL ||
564 fscanf(f," %ld",&no_primes) == 0 || no_primes < 1L ||
565 (prime_table = (INT *)SYM_calloc((int)no_primes,sizeof(INT))) == NULL )
566 {
567 no_primes = 0L;
568 return(ERROR);
569 }
570 for (i=0L;i<no_primes;i++)
571 if (fscanf(f," %ld",&prime_table[i]) != 1)
572 {
573 SYM_free(prime_table);
574 no_primes = 0L;
575 return(ERROR);
576 }
577 prime_table_set = 1L;
578 # else
579 no_primes = 15L;
580 if ((prime_table = (INT *)SYM_calloc((int)no_primes,sizeof(INT))) == NULL )
581 {
582 no_primes = 0L; /* The previous version had incompatible */
583 return(ERROR); /* uses for prime_table in the two parts */
584 } /* of the #ifdef construct */
585 prime_table[0] = 2L;
586 prime_table[1] = 3L;
587 prime_table[2] = 5L;
588 prime_table[3] = 7L;
589 prime_table[4] = 11L;
590 prime_table[5] = 13L;
591 prime_table[6] = 17L;
592 prime_table[7] = 19L;
593 prime_table[8] = 23L;
594 prime_table[9] = 29L;
595 prime_table[10] = 31L;
596 prime_table[11] = 37L;
597 prime_table[12] = 41L;
598 prime_table[13] = 43L;
599 prime_table[14] = 47L;
600 prime_table_set = 1L;
601
602 # endif /* PRIMES_FILE */
603
604 return(OK);
605 }
606
607 /* This routine factorizes an integer using the table of primes. */
608 /* a -- the integer to be factored; */
609 /* l -- a list in which the prime factors of a, which are contain- */
610 /* ed in the table, and their exponents are inserted as monomials */
611 /* with the primes as the selfs and the exponents as the koeffs; */
612 /* g -- the remaining factor; */
613 /* m -- the last number tried as a factor. */
614 /* first_prime -- the variable to point to the first prime factor */
615 /* if that is all that is required. For a full factorization, it */
616 /* must be set to NULL. */
617 /* l must be a MONOPOLY. a,l,g,m and first_prime must be different. */
618
619 #ifdef MONOPOLYTRUE
integer_factor_0(a,l,g,m,first_prime)620 static INT integer_factor_0(a,l,g,m,first_prime) OP a,l,g,m, first_prime;
621 /* 04.04.90: TPMcD */
622 /* AK 200891 V1.3 */
623 {
624 INT i,erg = ERROR;
625 OP b = CALLOCOBJECT(), c = CALLOCOBJECT(), d = CALLOCOBJECT(),
626 e = CALLOCOBJECT(), k = CALLOCOBJECT(), f;
627
628 if (S_O_K(a) != INTEGER && S_O_K(a) != LONGINT)
629 goto exit_label;
630 if (S_O_K(l) != MONOPOLY)
631 goto exit_label;
632 if (a == b || a == c || a == l || a == first_prime || b == c || b == l
633 || a == first_prime || c == l || c == first_prime|| l == first_prime)
634 goto exit_label;
635 copy(a,b);
636 if (EQ(b,cons_eins) == TRUE)
637 {
638 f = CALLOCOBJECT();
639 m_sk_mo(cons_eins,cons_eins,f);
640 insert(f,l,add_koeff,NULL);
641 }
642 else if (LT(b,cons_null) == TRUE)
643 {
644 addinvers_apply(b); /* AK 090891 */
645 f = CALLOCOBJECT();
646 M_I_I(-1L,c);
647 m_sk_mo(c,cons_eins,f);
648 insert(f,l,add_koeff,NULL);
649 }
650 if (EINSP(b))
651 {
652 erg = OK;
653 copy(b,g);
654 copy(b,m);
655 goto exit_label;
656 }
657 nb_square_root(b,e); /* e = sqrt(a) */
658
659 CTTO(INTEGER,LONGINT,"integer_factor_0(i1)",e);
660 if (!prime_table_set)
661 if (setup_prime_table() == ERROR) goto exit_label;
662 for (i=0L;i<no_primes;i++)
663 {
664 m_i_i(0,k);
665 /* 29.10.91: TPMcD: type of prime_table changed back to INT[] */
666 m_i_i((INT)prime_table[i],c);
667 if (GT(c,e) == TRUE) /* all primes not greater than */
668 { /* sqrt(b) have been tried. */
669 if (first_prime != NULL)
670 {
671 copy(b,first_prime);
672 erg = OK;
673 goto exit_label;
674 }
675 f = CALLOCOBJECT();
676 m_i_i(1L,d);
677 m_sk_mo(b,d,f);
678 insert(f,l,add_koeff,NULL);
679 m_i_i(1L,b);
680 break;
681 }
682 mod(b,c,d);
683
684 while(NULLP(d))
685 {
686 if (first_prime != NULL)
687 {
688 copy(c,first_prime);
689 erg = OK;
690 goto exit_label;
691 }
692 INC(k);
693 ganzdiv(b,c,b);
694 mod(b,c,d);
695 nb_square_root(b,e); /* e = sqrt(b) */
696 }
697 if (GT(k,cons_null) == TRUE)
698 {
699 f = CALLOCOBJECT();
700 m_sk_mo(c,k,f);
701 insert(f,l,add_koeff,NULL);
702 }
703 if (EINSP(b)) break;
704
705 }
706 erg = OK;
707 copy(b,g);
708 copy(c,m);
709 exit_label:
710 FREEALL(b);
711 FREEALL(c);
712 FREEALL(d);
713 FREEALL(e);
714 FREEALL(k);
715 CTO(MONOPOLY,"integer_factor_0(e2)",l);
716 CTO(INTEGER,"integer_factor_0(e3)",g);
717 ENDR("integer_factor_0");
718 }
719 #endif /* MONOPOLYTRUE */
720
721 /* This routine finds all prime factors between two bounds */
722 /* a -- the integer to be factored; */
723 /* l -- a list in which the prime factors of a, between the given */
724 /* bounds, and their exponents are inserted as monomials */
725 /* with the primes as the selfs and the exponents as the koeffs; */
726 /* if l is not a MONOPOLY, it is initialised to one; */
727 /* b -- the remaining factor; */
728 /* f_eins -- the lower bound on trial factors. */
729 /* If it is even, it is replaced by f_eins+1. */
730 /* f_zwei -- the upper bound on trial factors. */
731 /* first_prime -- the variable to point to the first prime factor */
732 /* if that is all that is required. For a full factorization, it */
733 /* must be set to NULL. */
734 /* l must be a MONOPOLY. a,l,b,f_eins,f_zwei,first_prime must be different. */
735
integer_factor_1(a,f_eins,f_zwei,b,l,first_prime)736 static INT integer_factor_1(a,f_eins,f_zwei,b,l,first_prime) OP a,f_eins,f_zwei,b,l,first_prime;
737 /* 04.04.90: TPMcD */
738 /* AK 200891 V1.3 */
739 {
740 INT flag = 1L;
741 INT erg = OK;
742 INT ret = ERROR, new_factor = 0L, myfirst = 1L;
743 #ifdef MONOPOLYTRUE
744 OP c = CALLOCOBJECT(), f = CALLOCOBJECT(), e = CALLOCOBJECT(),
745 e_eins = CALLOCOBJECT(), q = CALLOCOBJECT(), r = CALLOCOBJECT(),
746 k, g;
747 if (S_O_K(l) != MONOPOLY)
748 init(MONOPOLY,l);
749 copy(a,c);
750 if (LT(c,cons_null) == TRUE)
751 {
752 addinvers(c,c);
753 M_I_I(1L,f);
754 M_I_I(-1L,e);
755 k = CALLOCOBJECT();
756 m_sk_mo(e,f,k);
757 insert(k,l,add_koeff,NULL);
758 }
759 k = CALLOCOBJECT();
760 copy(f_zwei,e);
761 copy(f_eins,f);
762 if (LT(f,cons_eins) == TRUE) /* Make the initial divisor >= 2 */
763 copy(cons_zwei,f);
764 nb_quores(f,cons_zwei,q,r);
765 if (nullp(r))
766 {
767 if (einsp(q)) /* f = 2 */
768 flag = 0L;
769 else /* f is even and greater than 2 */
770 INC(f);
771 }
772 nb_quores(c,f,q,r);
773
774 while (LE(f,e) == TRUE)
775 {
776 while (nullp(r))
777 { /* The value of c entering this loop is divisible
778 by f exactly k times, where k refers to its
779 value exiting the loop. */
780 if (first_prime != NULL)
781 {
782 copy(f,first_prime);
783 ret = OK;
784 goto exit_label;
785 }
786 if (myfirst)
787 {
788 M_I_I(1L,k);
789 new_factor = 1L;
790 myfirst = 0L;
791 }
792 else
793 INC(k);
794 ganzdiv(c,f,c);
795 nb_quores(c,f,q,r);
796 }
797 if (new_factor)
798 { /* make new monomial and insert in the factor list */
799 g = CALLOCOBJECT();
800 m_sk_mo(f,k,g);
801 insert(g,l,add_koeff,NULL);
802 new_factor = 0L;
803 if (EQ(c,cons_eins) == TRUE)
804 break;
805 nb_square_root(c,e_eins);
806 /* reduce the upper limit of the trial factors */
807 if (LT(e_eins,e) == TRUE)
808 copy(e_eins,e);
809 /* the current c is a prime or it has prime factors */
810 /* are less than f_eins, and the factorization is 'complete'. */
811 if (LT(e,f) == TRUE)
812 break;
813 }
814 myfirst = 1L;
815 /* Increase f by 2 and find corresponding q and r */
816 INC(f);
817 if (flag) INC(f);
818 flag = 1L;
819 nb_quores(c,f,q,r);
820 }
821 copy(c,b);
822 ret = OK;
823 exit_label:
824 freeall(c);
825 freeall(q);
826 freeall(r);
827 freeall(f);
828 freeall(k);
829 freeall(e);
830 freeall(e_eins);
831 #else /* MONOPOLYTRUE */
832 error("integer_factor_1: MONOPOLY not available");
833 #endif /* MONOPOLYTRUE */
834 return(ret);
835 }
836
837 /* This is the main integer factorization routine. */
838 /* a -- the integer to be factored; */
839 /* l -- a list in which the prime factors of a and their exponents */
840 /* are inserted as monomials with the primes as the selfs and the */
841 /* exponents as the koeffs. l need not be initialized to a MONOPOLY. */
842
integer_factor(a,l)843 INT integer_factor(a,l) OP a,l;
844 /* 040490: TPMcD */
845 /* AK 200891 V1.3 */
846 /* 04.10.91: TPMcD */
847 {
848 INT erg = ERROR;
849
850 OP b,c,d,e;
851 CTTO(INTEGER,LONGINT,"integer_factor(1)",a);
852 CE2(a,l,integer_factor);
853
854 b = CALLOCOBJECT();
855 c = CALLOCOBJECT();
856 d = CALLOCOBJECT();
857 e = CALLOCOBJECT();
858
859 init(MONOPOLY,l);
860
861 /*First factorize using the list of primes in "PRIMES.DAT" */
862
863 if (integer_factor_0(a,l,b,c,NULL) != OK)
864 {
865 copy(a,b);
866 M_I_I(1L,c);
867 }
868 if (EQ(b,cons_eins) == TRUE) /* Factorization complete. */
869 {
870 erg = OK;
871 goto exit_label;
872 }
873 copy(b,d);
874 if (integer_factor_1(b,c,d,e,l,NULL) == OK)
875 {
876 /* If e > 1 , it is a prime greater than those in l */
877 m_i_i(1L,c);
878 if (gt(e,c) == TRUE)
879 {
880 m_sk_mo(e,c,d);
881 insert(d,l,add_koeff,NULL);
882 m_i_i(1L,e);
883 d = CALLOCOBJECT();
884 }
885 erg = OK;
886 }
887 else
888 printf("\ninteger_factor: factorization error");
889 exit_label:
890 FREEALL4(b,c,d,e);
891 CTO(MONOPOLY,"integer_factor(e2)",l);
892 ENDR("integer_factor");
893 }
894
895 /* This routine finds the smallest prime factor of an integer. */
896 /* a -- the integer; first_prime -- its smallest prime factor. */
897
first_prime_factor(a,first_prime)898 INT first_prime_factor(a,first_prime) OP a,first_prime;
899 /* 04.04.90: TPMcD */
900 /* AK 200891 V1.3 */
901 /* 04.10.91: TPMcD */
902 {
903 INT ret = ERROR;
904 #ifdef MONOPOLYTRUE
905 OP b = CALLOCOBJECT();
906 OP c = CALLOCOBJECT();
907 OP d = CALLOCOBJECT();
908 OP e = CALLOCOBJECT();
909 OP l = CALLOCOBJECT();
910
911 if (S_O_K(a) != INTEGER && S_O_K(a) != LONGINT)
912 goto exit_label;
913 init(MONOPOLY,l);
914 m_i_i(1L,first_prime);
915 copy(a,d);
916 if (LT(d,cons_null) == TRUE)
917 addinvers(d,d);
918 if (einsp(d))
919 {
920 ret = OK;
921 goto exit_label;
922 }
923 if (integer_factor_0(d,l,b,c,first_prime) == OK)
924 if (einsp(first_prime))
925 if (integer_factor_1(d,c,d,e,l,first_prime) != OK
926 || einsp(first_prime))
927 goto exit_label;
928 ret = OK;
929 exit_label:
930 if (ret != OK)
931 printf("\nfirst_prime_factor: factorization error");
932 freeall(b);
933 freeall(c);
934 freeall(d);
935 freeall(e);
936 freeall(l);
937 #else
938 error("integer_factor: MONOPOLY not available");
939 #endif
940 return(ret);
941 }
942
943 /* SQUARE-FREE PARTS */
944
945 /* This routine find the square-free part of the integer, which is */
946 /* given as a prime factors list. */
947 /* la -- a MONOPOLY containing the prime factorization of the integer */
948 /* lb, lc -- return the MONOPOLYs containing the prime factorization */
949 /* of the square-free and square parts, respectively. */
950 /* The parameters la,lb,lc must be distinct. */
951
square_free_part_0(la,lb,lc)952 INT square_free_part_0(la,lb,lc) OP la,lb,lc;
953 /* 14.06.90: TPMcD */
954 /* AK 200891 V1.3 */
955 {
956 INT erg = ERROR;
957 INT flag_b = 1L, flag_c = 1L;
958 OP u , x , y ,
959 z , ptr, w;
960
961 CTO(MONOPOLY,"square_free_part_0(1)",la);
962
963 u = CALLOCOBJECT();
964 x = CALLOCOBJECT();
965 y = CALLOCOBJECT();
966 z = CALLOCOBJECT();
967
968 ptr = la;
969 init(MONOPOLY,lb);
970 init(MONOPOLY,lc);
971 while (ptr != NULL)
972 {
973 copy(S_PO_S(ptr),u); /* the prime */
974 copy(S_PO_K(ptr),x); /* the exponent */
975 if (negp(x))
976 error("square_free_part_0 : unexpected negative exponent");
977 nb_quores(x,cons_zwei,z,y);
978 if (nullp(y)) /* even power */
979 {
980 w = CALLOCOBJECT();
981 m_sk_mo(u,z,w);
982 insert(w,lc,add_koeff,NULL);
983 flag_c = 0L;
984 }
985 else
986 {
987 if (not nullp(z))
988 {
989 w = CALLOCOBJECT();
990 m_sk_mo(u,z,w);
991 insert(w,lc,add_koeff,NULL);
992 flag_c = 0L;
993 }
994 w = CALLOCOBJECT();
995 m_sk_mo(u,cons_eins,w);
996 insert(w,lb,add_koeff,NULL);
997 flag_b = 0L;
998 }
999 ptr = S_L_N(ptr);
1000 }
1001 if (flag_b)
1002 {
1003 w = CALLOCOBJECT();
1004 m_sk_mo(cons_eins,cons_eins,w);
1005 insert(w,lb,add_koeff,NULL);
1006 }
1007 if (flag_c)
1008 {
1009 w = CALLOCOBJECT();
1010 m_sk_mo(cons_eins,cons_eins,w);
1011 insert(w,lc,add_koeff,NULL);
1012 }
1013 erg = OK;
1014
1015 freeall(u);
1016 freeall(x);
1017 freeall(y);
1018 freeall(z);
1019 CTO(MONOPOLY,"square_free_part_0(e2)",lb);
1020 CTO(MONOPOLY,"square_free_part_0(e3)",lc);
1021 ENDR("square_free_part_0");
1022 }
1023
1024 /* This routine find the square-free part of the integer, i.e. the */
1025 /* product of the prime factors (and -1 , if the integer is < 0. */
1026 /* a -- the integer, it is either an INTEGER, LONGINT or a */
1027 /* MONOPOLY containing the prime factorization of an integer. */
1028 /* b -- the square-free part. a = b * c 2 */
1029 /* c -- the square-root of the square part. */
1030 /* la -- returns a MONOPOLY containing the prime factorization of */
1031 /* a if a is not a MONOPOLY and la is not NULL. */
1032 /* lb, lc -- returns the MONOPOLYs containing the prime */
1033 /* factorizations of b and c. */
1034 /* The parameters a,b,c must be distinct. If la,lb,lc are not */
1035 /* NULL they must be distinct also. */
1036
square_free_part(a,b,c,la,lb,lc)1037 INT square_free_part(a,b,c,la,lb,lc) OP a,b,c,la,lb,lc;
1038 /* 14.06.90: TPMcD */
1039 /* AK 200891 V1.3 */
1040 /* 04.10.91: TPMcD */
1041 {
1042 INT erg = ERROR;
1043 OP la_tmp=NULL, lb_tmp=NULL, lc_tmp=NULL;
1044
1045 CTTTO(INTEGER,LONGINT,MONOPOLY,"square_free_part(1)",a);
1046
1047 if (S_O_K(a) == INTEGER || S_O_K(a) == LONGINT)
1048 {
1049 if (la == NULL)
1050 la_tmp = CALLOCOBJECT();
1051 else
1052 la_tmp = la;
1053 init(MONOPOLY,la_tmp);
1054 integer_factor(a,la_tmp);
1055 }
1056 else if (S_O_K(a) == MONOPOLY)
1057 la_tmp = a;
1058 else
1059 goto exit_label;
1060 if (lb == NULL)
1061 lb_tmp = CALLOCOBJECT();
1062 else
1063 lb_tmp = lb;
1064 init(MONOPOLY,lb_tmp);
1065 if (lc == NULL)
1066 lc_tmp = CALLOCOBJECT();
1067 else
1068 lc_tmp = lc;
1069 init(MONOPOLY,lc_tmp);
1070 square_free_part_0(la_tmp,lb_tmp,lc_tmp);
1071 integer_factors_to_integer(lb_tmp,b);
1072 integer_factors_to_integer(lc_tmp,c);
1073 erg = OK;
1074 exit_label:
1075 if (la == NULL && la_tmp != a) freeall(la_tmp);
1076 if (lb == NULL) freeall(lb_tmp);
1077 if (lc == NULL) freeall(lc_tmp);
1078
1079
1080
1081 CTTO(INTEGER,LONGINT,"square_free_part(e2)",b);
1082 CTTO(INTEGER,LONGINT,"square_free_part(e3)",c);
1083 ENDR("square_free_part");
1084 }
1085
1086 /*
1087 The Jacobi Symbol: (a/b) b odd.
1088 a and b are integers. c must point to a location different from a and b.
1089 if a and b have a common factor, c is set to 0 and ERROR is returned.
1090 otherwise, c is set to the the jacobi symbol (a/b). Note that b must
1091 be odd.
1092 */
1093
jacobi(a,b,c)1094 INT jacobi(a,b,c) OP a,b,c;
1095 /* AK 200891 V1.3 */
1096 {
1097 INT erg = OK;
1098 OP x , y , z , w , y_eins , y_zwei ,
1099 z_eins , z_zwei , four , eight ;
1100 INT ret = ERROR;
1101 INT d, f = 0L;
1102
1103 CTTO(INTEGER,LONGINT,"jacobi(1)",a);
1104 CTTO(INTEGER,LONGINT,"jacobi(2)",b);
1105 CE3(a,b,c,jacobi);
1106 M_I_I(0,c);
1107
1108 if (EVEN(b))
1109 {
1110 printf("Jacobi Symbol: Second integer must be odd\n");
1111 goto e2;
1112 }
1113 x =CALLOCOBJECT(), y = CALLOCOBJECT(), z = CALLOCOBJECT(),
1114 w = CALLOCOBJECT(), y_eins = CALLOCOBJECT(), y_zwei = CALLOCOBJECT(),
1115 z_eins = CALLOCOBJECT(), z_zwei = CALLOCOBJECT(),
1116 four = CALLOCOBJECT(), eight = CALLOCOBJECT();
1117 M_I_I(4,four);
1118 M_I_I(8,eight);
1119
1120 COPY(a,x);
1121 COPY(b,y);
1122 if (NEGP(y)) ADDINVERS_APPLY(y);
1123
1124 while (NEQ(y,cons_eins))
1125 {
1126 mod(x,y,z);
1127 if (NULLP(z)) /* The numbers not relatively prime*/
1128 goto exit_label;
1129 if (ODD(z))
1130 {
1131 FREESELF(w);
1132 ADDINVERS(y,w);
1133 ADD_APPLY(w,z);
1134 }
1135 d = 0L;
1136 while (EVEN(z))
1137 {
1138 d = 1L - d;
1139 ganzdiv(z,cons_zwei,z);
1140 if (NULLP(z)) /*The numbers not relatively prime*/
1141 goto exit_label;
1142 }
1143 FREESELF(x);COPY(y,x);
1144 FREESELF(y_eins);COPY(y,y_eins);
1145 DEC(y_eins);
1146 mod(y_eins,eight,y_zwei);
1147 FREESELF(z_eins);COPY(z,z_eins);
1148 DEC(z_eins);
1149 mod(z_eins,eight,z_zwei);
1150 FREESELF(y_eins);MULT(y_zwei,z_zwei,y_eins);
1151 mod(y_eins,eight,y_eins);
1152 if (GE(y_eins,four))
1153 f = 1L - f;
1154 if (d) /* an odd power of two */
1155 {
1156 INC(y);
1157 mod(y,eight,y_eins);
1158 if (GE(y_eins,four))
1159 f = 1L - f;
1160 }
1161 if (NEGP(z))
1162 ADDINVERS_APPLY(z);
1163 FREESELF(y);COPY(z,y);
1164 }
1165 m_i_i(1,c);
1166 if (f)
1167 M_I_I(-1,c);
1168 ret = OK;
1169 exit_label:
1170 FREEALL(x);
1171 FREEALL(y);
1172 FREEALL(z);
1173 FREEALL(w);
1174 FREEALL(y_eins);
1175 FREEALL(y_zwei);
1176 FREEALL(z_eins);
1177 FREEALL(z_zwei);
1178 FREEALL(four);
1179 FREEALL(eight);
1180 e2:
1181 return(ret);
1182 ENDR("jacobi");
1183 }
1184
1185 /*
1186 The Kronecker Symbol: (a/b). a square-free and congruent to 0 or 1 mod 4.
1187 a and b are integers. c must point to a location different from a and b.
1188 if a and b have a common factor, c is set to 0 and ERROR is returned.
1189 otherwise, c is set to the the jacobi symbol (a/b). Note that b must
1190 be odd.
1191 */
1192
kronecker(a,b,c)1193 INT kronecker(a,b,c) OP a,b,c;
1194 /* AK 200891 V1.3 */
1195 {
1196 INT flag = 0L;
1197 INT ret = ERROR;
1198 INT erg = OK;
1199 OP a_eins = CALLOCOBJECT(), a_zwei = CALLOCOBJECT(), b_null = CALLOCOBJECT(),
1200 b_eins = CALLOCOBJECT(), b_zwei = CALLOCOBJECT(), b_drei = CALLOCOBJECT(),
1201 four = CALLOCOBJECT();
1202
1203 if (c == a || c == b || nullp(a) || nullp(b))
1204 goto exit_label;
1205 COPY(a,a_eins);
1206 COPY(b,b_eins);
1207 if (NEGP(b_eins))
1208 ADDINVERS_APPLY(b_eins);
1209 M_I_I(4L,four);
1210 copy(cons_null,c);
1211 mod(a_eins,four,a_zwei);
1212 if (NULLP(a_zwei))
1213 flag = 1L;
1214 else if (!einsp(a_zwei))
1215 goto exit_label;
1216 nb_quores(b_eins,cons_zwei,b_zwei,b_drei);
1217 copy(cons_null,b_null);
1218 if (NULLP(b_drei)) /* b is even. */
1219 {
1220 if (flag) /* a is even also. */
1221 goto exit_label;
1222 do
1223 {
1224 INC(b_null);
1225 copy(b_zwei,b_eins);
1226 nb_quores(b_eins,cons_zwei,b_zwei,b_drei);
1227 } while (nullp(b_drei));
1228 } /* b_eins is the largest odd factor of b. */
1229 nb_quores(b_null,cons_zwei,b_null,b_zwei);
1230 if (NULLP(b_zwei)) /* b is divisible by an odd power of two and a is odd. */
1231 {
1232 m_i_i(8L,b_drei);
1233 mod(a_eins,b_drei,b_zwei);
1234 if (NEQ(b_zwei,cons_eins))
1235 flag = 1L; /* negate the final value */
1236 }
1237 else
1238 flag = 0L;
1239 /* At this point, b_eins is odd */
1240 jacobi(a_eins,b_eins,c);
1241 if (flag)
1242 ADDINVERS_APPLY(c);
1243 ret = OK;
1244 exit_label:
1245 FREEALL(a_eins);
1246 FREEALL(a_zwei);
1247 FREEALL(four);
1248 FREEALL(b_null);
1249 FREEALL(b_eins);
1250 FREEALL(b_zwei);
1251 FREEALL(b_drei);
1252 return(ret);
1253 ENDR("kronecker");
1254 }
1255 /****************** fields_0.c **********************/
1256 /* 26.07.91: TPMcD. */
1257 /*************************************************************/
eq_cyclotomic(a,b)1258 INT eq_cyclotomic(a,b) OP a,b;
1259 /* AK 210202 */
1260 {
1261 INT erg = OK,r;
1262 OP c;
1263 CTO(CYCLOTOMIC,"eq_cyclotomic(1)",a);
1264 if (S_O_K(b) != CYCLOTOMIC) return FALSE;
1265 c = CALLOCOBJECT();
1266 sub(a,b,c);
1267 r = NULLP(c);
1268 FREEALL(c);
1269 return r;
1270 ENDR("eq_cyclotomic");
1271 }
1272
eq_sqrad(a,b)1273 INT eq_sqrad(a,b) OP a,b;
1274 /* AK 180702 */
1275 {
1276 INT erg = OK,r;
1277 OP c;
1278 CTO(SQ_RADICAL,"eq_sqrad(1)",a);
1279 if (S_O_K(b) != SQ_RADICAL) return FALSE;
1280 c = CALLOCOBJECT();
1281 sub(a,b,c);
1282 r = NULLP(c);
1283 FREEALL(c);
1284 return r;
1285 ENDR("eq_sqrad");
1286 }
1287
eq_fieldobject_int(type,a,i)1288 INT eq_fieldobject_int(type,a,i) OBJECTKIND type; OP a; INT i;
1289 /* AK 200891 V1.3 */
1290 {
1291 INT ret = FALSE;
1292 OP b = CALLOCOBJECT();
1293 OP c = CALLOCOBJECT();
1294 INT erg = OK;
1295 M_I_I(-i,b);
1296 switch(S_O_K(a))
1297 {
1298 #ifdef MONOPOLYTRUE
1299 case MONOPOLY:
1300 add_scalar_monopoly(b,a,c);
1301 ret = nullp_monopoly(c);
1302 break;
1303 #endif /* MONOPOLYTRUE */
1304 #ifdef CYCLOTRUE
1305 case CYCLOTOMIC:
1306 add_scalar_cyclo(b,a,c);
1307 ret = nullp_monopoly(S_N_S(c));
1308 break;
1309 #endif /* CYCLOTRUE */
1310 #ifdef SQRADTRUE
1311 case SQ_RADICAL:
1312 add_scalar_sqrad(b,a,c);
1313 ret = nullp_monopoly(S_N_S(c));
1314 break;
1315 #endif /* SQRADTRUE */
1316 default:
1317 error("eq_fieldobject_int: invalid type\n");
1318 }
1319 FREEALL(b);
1320 FREEALL(c);
1321 return(ret);
1322 ENDR("eq_fieldobject_int");
1323 }
1324
1325 #ifdef NUMBERTRUE
1326
callocnumber()1327 static struct number * callocnumber()
1328 /* 22.06.90: TPMcD */
1329 /* AK 200891 V1.3 */
1330 {
1331 struct number *result=(struct number *)
1332 SYM_calloc(1,sizeof(struct number));
1333 if (result == NULL) error("callocnumber:no mem");
1334 number_mem++;
1335 return(result);
1336 }
1337 #endif /* NUMBERTRUE */
1338
m_ksd_n(kind,self,data,result)1339 INT m_ksd_n(kind,self,data,result) OBJECTKIND kind; OP self,data,result;
1340 /* AK 230191 */
1341 /* AK 200891 V1.3 */
1342 {
1343 INT erg = ERROR;
1344 #ifdef NUMBERTRUE
1345 erg = b_ksd_n(kind,CALLOCOBJECT(),CALLOCOBJECT(),result);
1346 if ((S_O_K(self) != MONOPOLY) ||
1347 ((kind == SQ_RADICAL) && (S_O_K(data) != LIST)))
1348 return( error("m_ksd_n: invalid self or data") );
1349 erg += copy(self,S_N_S(result));
1350 erg += copy(data,S_N_D(result));
1351 #endif /* NUMBERTRUE */
1352 return erg;
1353 }
1354
init_sqrad(a)1355 INT init_sqrad(a) OP a;
1356 /* AK 010993 */
1357 {
1358 INT erg = OK;
1359 CTO(EMPTY,"init_sqrad(1)",a);
1360 erg += b_ksd_n(SQ_RADICAL, CALLOCOBJECT(), CALLOCOBJECT(), a);
1361 ENDR("init_sqrad");
1362 }
1363
1364
1365
init_cyclo(a)1366 INT init_cyclo(a) OP a;
1367 /* AK 010993 */
1368 {
1369 INT erg = OK;
1370 CTO(EMPTY,"init_cyclo(1)",a);
1371
1372 erg += b_ksd_n(CYCLOTOMIC, CALLOCOBJECT(), NULL, a);
1373
1374 ENDR("init_cyclo");
1375 }
1376
1377
1378 #ifdef NUMBERTRUE
b_ksd_n(kind,self,data,result)1379 INT b_ksd_n(kind,self,data,result) OBJECTKIND kind; OP self,data,result;
1380 /* 22.06.90: TPMcD */ /* 3.04.91: TPMcD. */
1381 /* AK 200891 V1.3 */
1382 {
1383 OBJECTSELF obself;
1384 if (not EMPTYP(result))
1385 freeself(result);
1386 obself.ob_number = callocnumber();
1387 b_ks_o(kind,obself,result);
1388 if (EMPTYP(self))
1389 init(MONOPOLY,self);
1390 if (kind == SQ_RADICAL && EMPTYP(data))
1391 init(LIST,data);
1392 if ((S_O_K(self) != MONOPOLY) ||
1393 ((kind == SQ_RADICAL) && (S_O_K(data) != LIST)))
1394 return( error("b_ksd_n: invalid self or data") );
1395 S_N_S(result) = self;
1396 S_N_D(result) = data;
1397 return(OK);
1398 }
1399 #endif /* NUMBERTRUE */
1400
objectwrite_number(f,number)1401 INT objectwrite_number(f,number) FILE *f; OP number;
1402 /* AK 200891 V1.3 */
1403 {
1404 #ifdef NUMBERTRUE
1405 fprintf(f, " %" PRIINT "\n" ,(INT)S_O_K(number));
1406 objectwrite(f,S_N_S(number));
1407 switch (S_O_K(number))
1408 {
1409 case SQ_RADICAL:
1410 objectwrite(f,S_N_D(number));
1411 break;
1412 case CYCLOTOMIC:
1413 objectwrite(f,S_N_DCI(number));
1414 break;
1415 default:
1416 error("Invalid number type\n");
1417 }
1418 return(OK);
1419 #else /* NUMBERTRUE */
1420 error("objectwrite_number:NUMBER not available");
1421 return(ERROR);
1422 #endif /* NUMBERTRUE */
1423 }
1424
objectread_number(f,number,type)1425 INT objectread_number(f,number,type) FILE *f; OP number; OBJECTKIND type;
1426 /* AK 200891 V1.3 */
1427 {
1428 #ifdef NUMBERTRUE
1429 init(type,number);
1430 objectread(f,S_N_S(number));
1431 switch (S_O_K(number))
1432 {
1433 case SQ_RADICAL:
1434 objectread(f,S_N_D(number));
1435 break;
1436 case CYCLOTOMIC:
1437 {
1438 OP index = CALLOCOBJECT();
1439 objectread(f,index);
1440 S_N_DC(number) = add_cyclo_data(index);
1441 }
1442 break;
1443 default:
1444 error("Invalid number type\n");
1445 }
1446 return(OK);
1447 #endif /* NUMBERTRUE */
1448 }
1449
fprint_number(f,n)1450 INT fprint_number(f,n) FILE *f; OP n;
1451 /* AK 200891 V1.3 */
1452 {
1453 INT saving;
1454 #ifdef NUMBERTRUE
1455 switch (S_O_K(n))
1456 {
1457 case SQ_RADICAL:
1458 /* Are all coefficients fractions with denominator 2 */
1459 if (S_O_K(S_PO_K(S_N_S(n))) == BRUCH)
1460 {
1461 OP nn = CALLOCOBJECT();
1462 saving = space_saving;
1463 space_saving = FALSE;
1464 mult_scalar_sqrad(cons_zwei,n,nn);
1465 space_saving = saving;
1466 if (integer_coefficients(S_N_S(nn)) == TRUE)
1467 {
1468 fprintf(f," 1/2 (");
1469 fprint_sqrad(f,nn);
1470 fprintf(f,")");
1471 freeall(nn);
1472 zeilenposition += 7L;
1473 return(OK);
1474 }
1475 freeall(nn);
1476 }
1477 fprintf(f," ( ");
1478 fprint_sqrad(f,n);
1479 fprintf(f," )");
1480 zeilenposition += 5L;
1481 break;
1482 case CYCLOTOMIC:
1483 fprintf(f," ( ");
1484 fprint_cyclo(f,n);
1485 fprintf(f," )");
1486 zeilenposition += 5L;
1487 break;
1488 default:
1489 ;
1490 }
1491 #endif /* NUMBERTRUE */
1492 return(OK);
1493 }
1494
freeself_number(n)1495 INT freeself_number(n) OP n;
1496 /* AK 200891 V1.3 */
1497 {
1498 #ifdef NUMBERTRUE
1499 OBJECTSELF d;
1500 INT erg = OK;
1501 d = S_O_S(n);
1502 erg = freeall(S_N_S(n));
1503 if (erg == ERROR)
1504 return error("freeself_number:error in freeall S_N_S");
1505
1506 switch (S_O_K(n))
1507 {
1508 case SQ_RADICAL:
1509 if (not EMPTYP(S_N_D(n)))
1510 freeall(S_N_D(n));
1511 else
1512 error("freeself_number: no data to release");
1513 break;
1514 case CYCLOTOMIC:
1515 break;
1516 default:
1517 ;
1518 }
1519 SYM_free(d.ob_number);
1520 number_mem--;
1521 C_O_K(n,EMPTY);
1522 #endif /* NUMBERTRUE */
1523 return OK;
1524 }
1525
comp_number(a,b)1526 INT comp_number(a,b) OP a,b;
1527 /* 21.07.91 TPMcD: still incomplete */
1528 /* AK 200891 V1.3 */
1529 {
1530 #ifdef NUMBERTRUE
1531 switch (S_O_K(a))
1532 {
1533 case SQ_RADICAL:
1534 comp_sqrad(a,b);
1535 break;
1536 case CYCLOTOMIC:
1537 comp_cyclo(a,b);
1538 break;
1539 default:
1540 return error("comp_number:invalid number type\n");
1541 }
1542 return(OK);
1543 #else /* NUMBERTRUE */
1544 return error("comp_number:NUMBER not available");
1545 #endif /* NUMBERTRUE */
1546 }
1547
copy_number(a,b)1548 INT copy_number(a,b) OP a,b;
1549 /* AK 200891 V1.3 */
1550 /* AK 060498 V2.0 */
1551 {
1552 #ifdef NUMBERTRUE
1553 if (a == b)
1554 error("copy_number: First and second arguments are the same\n");
1555 init(S_O_K(a),b);
1556 if (S_N_S(a) != NULL) copy(S_N_S(a),S_N_S(b));
1557 switch (S_O_K(a))
1558 {
1559 case SQ_RADICAL:
1560 copy(S_N_D(a),S_N_D(b));
1561 break;
1562 case CYCLOTOMIC:
1563 S_N_DC(b) = S_N_DC(a);
1564 break;
1565 default:
1566 return error("copy_number:invalid number type\n");
1567 }
1568 return(OK);
1569 #endif /* NUMBERTRUE */
1570 }
1571
complex_conjugate(a,b)1572 INT complex_conjugate(a,b) OP a,b;
1573 {
1574 OP ptr;
1575 if (a != b)
1576 copy(a,b);
1577 #ifdef NUMBERTRUE
1578 switch (S_O_K(b))
1579 {
1580 case SQ_RADICAL:
1581 ptr = S_N_S(b);
1582 while (ptr != NULL)
1583 {
1584 if (not EMPTYP(S_PO_K(ptr)))
1585 complex_conjugate(S_PO_K(ptr),S_PO_K(ptr));
1586 if (LT(S_PO_S(ptr),cons_null) == TRUE)
1587 addinvers_apply(S_PO_K(ptr));
1588 ptr = S_L_N(ptr);
1589 }
1590 break;
1591 case CYCLOTOMIC:
1592 ptr = S_N_S(b);
1593 while (ptr != NULL)
1594 {
1595 if (not EMPTYP(S_PO_K(ptr)))
1596 {
1597 complex_conjugate(S_PO_K(ptr),S_PO_K(ptr));
1598 addinvers_apply(S_PO_S(ptr));
1599 add_apply(S_N_DCI(b),S_PO_S(ptr));
1600 }
1601 ptr = S_L_N(ptr);
1602 }
1603 break;
1604 default:
1605 break;
1606 }
1607 #endif /* NUMBERTRUE */
1608 return(OK);
1609 }
1610
setup_numbers(basis,saving,filename)1611 INT setup_numbers(basis,saving,filename) INT basis, saving;
1612 char *filename;
1613 /* 29.10.91: TPMcD */
1614 {
1615 #ifdef CYCLOTRUE
1616 number_mem = (INT)0;
1617 #endif /* CYCLOTRUE */
1618 #ifdef MONOPOLYTRUE
1619 setup_prime_table();
1620 #endif /* MONOPOLYTRUE */
1621 #ifdef NUMBERTRUE
1622 basis_type = basis;
1623 space_saving = saving;
1624 setup_cyclotomic_table(filename);
1625 #endif /* NUMBERTRUE */
1626 return(OK);
1627 }
1628
release_numbers()1629 INT release_numbers()
1630 /* 29.10.91: TPMcD */
1631 {
1632 #ifdef MONOPOLYTRUE
1633 if (prime_table_set)
1634 {
1635 SYM_free(prime_table);
1636 prime_table = NULL;
1637 }
1638 #endif
1639 #ifdef NUMBERTRUE
1640 if (cyclo_table_set)
1641 {
1642 free_cyclo_table();
1643 SYM_free(zzcyclo_table);
1644 cyclo_table_set = 0; /* AK 120202 */
1645 }
1646 if (cyclo_list_set)
1647 {
1648 free_cyclo_list();
1649 freeall(zzcyclo_list);
1650 cyclo_list_set = 0; /* AK 120202 */
1651 }
1652 #endif
1653 return(OK);
1654 }
1655
nb_ende()1656 INT nb_ende()
1657 {
1658 #ifdef CYCLOTRUE
1659 if (number_mem != 0L)
1660 fprintf(stderr, "error in number memory %" PRIINT "\n" ,number_mem);
1661 return OK;
1662 #endif
1663 }
1664
reset_basis(basis)1665 INT reset_basis(basis) INT basis;
1666 {
1667 #ifdef NUMBERTRUE
1668 basis_type = basis;
1669 if (basis == NO_REDUCE || basis == POWER_REDUCE)
1670 printf("\nWARNING: not a valid basis\n");
1671 #endif
1672 return(OK);
1673 }
1674
reset_saving(saving)1675 INT reset_saving(saving) INT saving;
1676 {
1677 #ifdef NUMBERTRUE
1678 space_saving = saving;
1679 #endif
1680 return(OK);
1681 }
1682
tex_monom_plus(form,a)1683 INT tex_monom_plus(form,a) INT form; OP a;
1684 /* AK 200891 V1.3 */
1685 {
1686 return tex_monom(a);/* AK return */
1687 }
1688
1689 #ifdef UNDEF
1690 /*Multiplies the entries in two lists pairwise, putting the resulting */
1691 /*objects in a list. Duplicate objects are ignored. */
1692
mult_lists(a,b,c)1693 static INT mult_lists(a,b,c) OP a, b, c;
1694 /* 26.08.90: TPMcD. */ /* AK 200891 V1.3 */ /* 04.10.91: TPMcD */
1695 {
1696 INT flag = 0L;
1697 #ifdef LISTTRUE
1698 OP new, a_ptr, b_ptr, c_tmp;
1699
1700 if (c == a || c == b)
1701 {
1702 flag = 1L;
1703 c_tmp = CALLOCOBJECT();
1704 }
1705 else
1706 c_tmp = c;
1707 init(LIST, c_tmp);
1708 b_ptr = b;
1709 while (b_ptr != NULL)
1710 {
1711 a_ptr = a;
1712 while (a_ptr != NULL)
1713 {
1714 new = CALLOCOBJECT();
1715 mult(S_L_S(a_ptr), S_L_S(b_ptr), new);
1716 insert(new,c_tmp,NULL,NULL);
1717 a_ptr = S_L_N(a_ptr);
1718 };
1719 b_ptr = S_L_N(b_ptr);
1720 }
1721 if (flag)
1722 {
1723 copy(c_tmp,c);
1724 freeall(c_tmp);
1725 }
1726 return(OK);
1727 #else
1728 error("mult_lists: LIST not available");
1729 return(ERROR);
1730 #endif
1731 }
1732 #endif
1733
tidy(a)1734 INT tidy(a) OP a;
1735 /* AK 200891 V1.3 */
1736 /* 04.10.91: TPMcD */
1737 {
1738 INT i,j;
1739 OP ptr;
1740 switch (S_O_K(a))
1741 {
1742 #ifdef BRUCHTRUE
1743 case BRUCH :
1744 tidy(S_B_O(a));
1745 tidy(S_B_U(a));
1746 break;
1747 #endif /* BRUCHTRUE */
1748 case INTEGER :
1749 break;
1750 #ifdef LISTTRUE
1751 case POLYNOM:
1752 case LIST :
1753 ptr = a;
1754 while (ptr != NULL)
1755 {
1756 tidy(S_L_S(ptr));
1757 ptr = S_L_N(ptr);
1758 }
1759 break;
1760 #endif /* LISTTRUE */
1761 #ifdef LONGINTTRUE
1762 case LONGINT :
1763 break;
1764 #endif /* LONGINTTRUE */
1765 #ifdef MATRIXTRUE
1766 case MATRIX :
1767 for (i=0L;i<S_M_LI(a);i++)
1768 for (j=0L;j<S_M_HI(a);j++)
1769 tidy(S_M_IJ(a,i,j));
1770 break;
1771 #endif /* MATRIXTRUE */
1772 #ifdef MONOMTRUE
1773 case MONOM :
1774 tidy(S_MO_S(a));
1775 tidy(S_MO_K(a));
1776 break;
1777 #endif /* MONOMTRUE */
1778 #ifdef NUMBERTRUE
1779 case SQ_RADICAL:
1780 break;
1781 case CYCLOTOMIC:
1782 standardise_cyclo_0(a,basis_type);
1783 break;
1784 #endif /* NUMBERTRUE */
1785 #ifdef VECTORTRUE
1786 case VECTOR :
1787 for (i=0L;i<S_V_LI(a);i++)
1788 tidy(S_V_I(a,i));
1789 break;
1790 #endif /* VECTORTRUE */
1791 default:
1792 return error("tidy: invalid type\n");
1793 }
1794 return(OK);
1795 }
1796 /****************** fields_1.c **********************/
1797 /* 26.07.91: TPMcD. */
1798 /*************************************************************/
1799
1800 /* MONOPOLYS */
1801
1802 #ifdef MONOPOLYTRUE
1803
eval_monopoly(a,b,c)1804 INT eval_monopoly(a,b,c) OP a,b,c;
1805 /*AK 290395 */
1806 {
1807 INT erg = OK;
1808 OP z,d;
1809
1810 CTO(MONOPOLY,"eval_monopoly",a);
1811 CE3(a,b,c,eval_monopoly);
1812
1813 erg += m_i_i((INT)0,c);
1814 z = a;
1815 d = CALLOCOBJECT();
1816 while (z != NULL)
1817 {
1818 CTO(INTEGER,"eval_monopoly:non integer self part",
1819 S_MO_S(S_L_S(z)));
1820 erg += hoch(b,S_MO_S(S_L_S(z)),d);
1821 MULT_APPLY(S_MO_K(S_L_S(z)),d);
1822 erg += add_apply(d,c);
1823 z = S_L_N(z);
1824 }
1825 erg += freeall(d);
1826
1827 ENDR("eval_monopoly");
1828 }
1829
einsp_monopoly(a)1830 INT einsp_monopoly(a) OP a;
1831 {
1832 if (S_L_S(a) == NULL)
1833 return FALSE;
1834 if (S_L_N(a) != NULL)
1835 return FALSE;
1836 if (not nullp(S_MO_S(S_L_S(a))))
1837 return FALSE;
1838 if (not einsp(S_MO_K(S_L_S(a))))
1839 return FALSE;
1840 return TRUE;
1841 }
1842
m_skn_mp(s,k,n,e)1843 INT m_skn_mp(s,k,n,e) OP s,k,n,e;
1844 /* make self koeff next monopoly */
1845 /* AK 040291 */ /* AK 200891 V1.3 */
1846 /* AK 080306 V3.0 */
1847 {
1848 INT erg=OK;
1849 COP("m_skn_mp(4)",e);
1850
1851 if ( n == NULL) {
1852 erg += b_skn_mp(CALLOCOBJECT(),CALLOCOBJECT(),NULL,e);
1853 COPY(s,S_PO_S(e));
1854 COPY(k,S_PO_K(e));
1855 }
1856 else {
1857 erg += b_skn_mp(CALLOCOBJECT(),CALLOCOBJECT(),CALLOCOBJECT(),e);
1858 COPY(s,S_PO_S(e));
1859 COPY(k,S_PO_K(e));
1860 COPY(n,S_PO_N(e));
1861 }
1862 ENDR("m_skn_mp");
1863 }
1864
1865
1866
b_skn_mp(s,k,n,e)1867 INT b_skn_mp(s,k,n,e) OP s,k,n,e;
1868 /* build self koeff next monopoly */
1869 /* AK 040291 */ /* AK 200891 V1.3 */ /* 231091: TPMcD */
1870 /* AK 080306 V3.0 */
1871 {
1872 INT erg=OK;
1873 COP("b_skn_mp(4)",e);
1874
1875 FREESELF(e);
1876 erg += b_sn_l(CALLOCOBJECT(),n,e);
1877 C_O_K(e,MONOPOLY);
1878 erg += b_sk_mo(s,k,S_L_S(e));
1879 ENDR("b_skn_mp");
1880 }
1881
cast_apply_monopoly(a)1882 INT cast_apply_monopoly(a) OP a;
1883 /* tries to convert the parameter in to a MONOPOLY object */
1884 /* AK 200498 V2.0 */
1885 /* AK 271005 V3.1 */
1886 {
1887 INT erg = OK;
1888 EOP("cast_apply_monopoly(1)",a); /* check if empty object */
1889 {
1890 switch (S_O_K(a))
1891 {
1892 case BRUCH:
1893 case INTEGER:
1894 case LONGINT:
1895 case FF:
1896 erg += m_skn_mp(cons_null,a,NULL,a);
1897 break;
1898 case POLYNOM:
1899 NYI("cast_apply_monopoly:POLYNOM->MONOPOLY");
1900 break;
1901 default:
1902 erg += WTO("cast_apply_monopoly:can not convert",a);
1903 break;
1904 }
1905 }
1906 ENDR("cast_apply_monopoly");
1907 }
1908
1909
1910
1911
scan_monopoly(a)1912 INT scan_monopoly(a) OP a;
1913 /* AK 200990 */ /* AK 220191 V1.2 */
1914 /* a becomes a monopoly */ /* AK 200891 V1.3 */
1915 {
1916 OBJECTKIND kt,st;
1917 INT erg = OK;
1918 CTO(EMPTY,"scan_monopoly(1)",a);
1919
1920 erg += printeingabe("type of monopoly self ");
1921 st=scanobjectkind();
1922 erg += printeingabe("type of monopoly koeff ");
1923 kt=scanobjectkind();
1924 erg += SCMPCO(st,kt,a);
1925
1926 ENDR("scan_monopoly");
1927 }
1928 #endif /* MONOPOLYTRUE */
1929
SCMPCO(self_type,coeff_type,result)1930 static INT SCMPCO(self_type,coeff_type,result) /* scan_monopoly_co */
1931 OBJECTKIND self_type,coeff_type;
1932 OP result;
1933 /* 04.06.90: TPMcD. */ /* 1.04.91: TPMcD. */
1934 /* AK 200891 V1.3 */ /* 23.10.91: TPMcD */
1935 {
1936 INT ret = ERROR;
1937 INT i,n;
1938 # ifdef MONOPOLYTRUE
1939 OP x = CALLOCOBJECT(), y = CALLOCOBJECT(), z;
1940 char a[30];
1941
1942 init(MONOPOLY,result);
1943 printeingabe("Length of list: "); /* AK 080891 */
1944 scanf( "%" PRIINT ,&n);
1945 for (i=0L;i<n;i++)
1946 {
1947 sprintf(a, "%" PRIINT "-th monomial (self) " ,i);
1948 printeingabe(a);
1949 scan(self_type,x);
1950 sprintf(a, "%" PRIINT "-th monomial (koeff) " ,i);
1951 printeingabe(a);
1952 scan(coeff_type,y);
1953 if (nullp(y))
1954 continue;
1955 z = CALLOCOBJECT();
1956 m_sk_mo(x,y,z);
1957 insert(z,result,add_koeff,NULL);
1958 }
1959 if (empty_listp(result))
1960 insert_zero_into_monopoly(result);
1961 freeall(x);
1962 freeall(y);
1963 ret = OK;
1964 #else
1965 error("SCMPCO: MONOPOLY not available");
1966 #endif
1967 return(ret);
1968 }
1969
1970 /* Removes those terms from a MONOPOLY with zero coefficients unless */
1971 /* this makes the list empty. In this case, one term with self and */
1972 /* koeff both 0 is left. */
1973
1974 # ifdef MONOPOLYTRUE
remove_zero_terms(a)1975 INT remove_zero_terms(a) OP a;
1976 /* 1.04.91: TPMcD. */
1977 /* AK 200891 V1.3 */
1978 {
1979 OP term, ptr = a, a_tmp;
1980 INT erg = OK;
1981
1982 CTO(MONOPOLY,"remove_zero_terms",a);
1983 /*
1984 a_tmp = CALLOCOBJECT();
1985 erg += init(MONOPOLY,a_tmp);
1986 if (!empty_listp(ptr))
1987 while (ptr != NULL)
1988 {
1989 if (not NULLP(S_PO_K(ptr)))
1990 {
1991 term = CALLOCOBJECT();
1992 COPY(S_L_S(ptr),term);
1993 insert(term,a_tmp,add_koeff,NULL);
1994 }
1995 ptr = S_L_N(ptr);
1996 }
1997 if (empty_listp(a_tmp))
1998 insert_zero_into_monopoly(a_tmp);
1999 SWAP(a_tmp,a);
2000 FREEALL(a_tmp);
2001 */
2002 ptr = a; term = NULL;
2003 again:
2004 if (!empty_listp(ptr))
2005 while (ptr != NULL)
2006 {
2007 if (not NULLP(S_PO_K(ptr)) )
2008 {
2009 term = ptr;
2010 ptr = S_L_N(ptr);
2011 }
2012 else {
2013 a_tmp = ptr;
2014 ptr = S_L_N(ptr);
2015 C_L_N(a_tmp, NULL);
2016 if (term == NULL)
2017 {
2018 /* first element will be removed */
2019 if (ptr == NULL) { M_I_I(0,S_PO_S(a)); goto ende; }
2020 FREESELF(a);
2021 *a = *ptr;
2022 C_O_K(ptr,EMPTY);
2023 FREEALL(ptr);
2024 ptr = a;
2025 goto again;
2026 }
2027 else {
2028 C_L_N(term,ptr);
2029 FREEALL(a_tmp);
2030 }
2031 }
2032 }
2033
2034 if (empty_listp(a))
2035 insert_zero_into_monopoly(a);
2036
2037 ende:
2038 ENDR("remove_zero_terms");
2039 }
2040
2041
2042 /* This is used to convert an empty monopoly into a non-empty one */
2043
2044
insert_zero_into_monopoly(a)2045 static INT insert_zero_into_monopoly(a) OP a;
2046 /* 1.04.91: TPMcD. */ /* AK 200891 V1.3 */
2047 {
2048 OP m = CALLOCOBJECT();
2049 b_sk_mo(CALLOCOBJECT(),CALLOCOBJECT(),m);
2050 M_I_I(0L,S_MO_K(m));
2051 M_I_I(0L,S_MO_S(m));
2052 insert(m,a,add_koeff,NULL);
2053 return(OK);
2054 }
2055 # endif /* MONOPOLYTRUE */
2056
integer_coefficients(a)2057 static INT integer_coefficients(a) OP a;
2058 /* 25.09.91: TPMcD. */
2059 {
2060 OP ptr = a;
2061 # ifdef MONOPOLYTRUE
2062 if (S_O_K(a) != MONOPOLY)
2063 error("integer_coefficients: parameter is not a MONOPOLY");
2064 if (!empty_listp(ptr))
2065 while (ptr != NULL)
2066 {
2067 if (S_O_K(S_PO_K(ptr)) != INTEGER && S_O_K(S_PO_K(ptr)) != LONGINT)
2068 return(FALSE);
2069 ptr = S_L_N(ptr);
2070 }
2071 #else
2072 error("integer_coefficients: MONOPOLY not available");
2073 #endif
2074 return(TRUE);
2075 }
2076
2077 # ifdef MONOPOLYTRUE
add_scalar_monopoly(a,b,c)2078 INT add_scalar_monopoly(a,b,c) OP a,b,c;
2079 /* 30.05.90: TPMcD. */ /* 1.04.91: TPMcD. */
2080 /* AK 200891 V1.3 */
2081 /* 04.10.91: TPMcD */
2082 {
2083 OP ptr;
2084 INT erg = OK;
2085 CTO(MONOPOLY,"add_scalar_monopoly",b);
2086 if (c == a)
2087 erg += ERROR;
2088 if (c != b)
2089 copy(b,c);
2090 ptr = CALLOCOBJECT();
2091 erg += init(MONOPOLY,ptr);
2092 C_L_S(ptr,CALLOCOBJECT());
2093 erg += m_sk_mo(cons_null,a,S_L_S(ptr));
2094 erg += add_apply(ptr,c);
2095 erg += remove_zero_terms(c);
2096 erg += freeall(ptr);
2097 ENDR("add_scalar_monopoly");
2098 }
2099 #endif
2100
mult_apply_scalar_monopoly(a,b)2101 INT mult_apply_scalar_monopoly(a,b) OP a,b;
2102 /* AK 210202 */
2103 {
2104 INT erg = OK;
2105 OP z;
2106 CTO(MONOPOLY,"mult_apply_scalar_monopoly(2)",b);
2107 if (NULLP(a))
2108 erg += init(MONOPOLY,b);
2109 else {
2110 FORALL(z,b,{MULT_APPLY(a,S_MO_K(z)); } );
2111 }
2112 ENDR("mult_apply_scalar_monopoly");
2113 }
2114
mult_scalar_monopoly(a,b,c)2115 INT mult_scalar_monopoly(a,b,c) OP a,b,c;
2116 /* 30.05.90: TPMcD. */ /* 1.04.91: TPMcD. */
2117 /* AK 200891 V1.3 */
2118 /* 04.10.91: TPMcD */
2119 /* AK 200498 V2.0 */
2120 {
2121 INT erg = OK;
2122 CTO(ANYTYPE,"mult_scalar_monopoly(1)",a);
2123 CTO(MONOPOLY,"mult_scalar_monopoly(2)",b);
2124 CTO(EMPTY,"mult_scalar_monopoly(3)",c);
2125 if (NULLP(a)) {
2126 erg += init(MONOPOLY,c);
2127 }
2128 else {
2129 OP z;
2130 COPY(b,c);
2131 FORALL(z,c, { MULT_APPLY(a,S_MO_K(z)); } );
2132 }
2133 CTO(MONOPOLY,"mult_scalar_monopoly(e3)",c);
2134 /*
2135 {
2136 OP d,z;
2137 d = CALLOCOBJECT();
2138 copy(c,d);
2139 FORALL(z,d, { div_apply(S_MO_K(z),a); } );
2140 if (NEQ(d,b)) {
2141 TCTO(ANYTYPE,"mult_scalar_monopoly(b)",b);
2142 TCTO(ANYTYPE,"mult_scalar_monopoly(d)",d);
2143 error("mult_scalar_monopoly:wrong result");
2144 }
2145 FREEALL(d);
2146 }
2147 */
2148 ENDR("mult_scalar_monopoly");
2149 }
2150
2151 # ifdef MONOPOLYTRUE
add_monopoly_monopoly(a,b,c)2152 INT add_monopoly_monopoly(a,b,c) OP a, b, c;
2153 /* 30.05.90: TPMcD. */ /* 1.04.91: TPMcD. */
2154 /* AK 200891 V1.3 */
2155 /* 04.10.91: TPMcD */
2156 {
2157 OP temp_eins, temp_zwei;
2158 INT erg = OK;
2159 CTO(MONOPOLY,"add_monopoly_monopoly(1)",a);
2160 CTO(MONOPOLY,"add_monopoly_monopoly(2)",b);
2161 CTO(EMPTY,"add_monopoly_monopoly(3)",c);
2162
2163 temp_eins = CALLOCOBJECT();
2164 temp_zwei = CALLOCOBJECT();
2165 copy(a,temp_eins);
2166 copy(b,temp_zwei);
2167 init(S_O_K(a), c);
2168 insert(temp_eins,c,add_koeff,NULL);
2169 insert(temp_zwei,c,add_koeff,NULL);
2170 erg += remove_zero_terms(c);
2171 ENDR("add_monopoly_monopoly");
2172 }
2173
2174
mult_monopoly_monopoly(a,b,c)2175 INT mult_monopoly_monopoly(a,b,c) OP a, b, c;
2176 /* 30.05.90: TPMcD. */ /* 1.04.91: TPMcD. */
2177 /* AK 200891 V1.3 */
2178 /* 04.10.91: TPMcD */
2179 {
2180 OP new, a_ptr, b_ptr, c_tmp;
2181 OP temp_eins, temp_zwei;
2182 INT erg = OK;
2183 CTO(MONOPOLY,"mult_monopoly_monopoly(1)",a);
2184 CTO(MONOPOLY,"mult_monopoly_monopoly(2)",b);
2185 CTO(EMPTY,"mult_monopoly_monopoly(3)",c);
2186
2187 c_tmp = c;
2188
2189 erg += init(MONOPOLY, c_tmp);
2190
2191 if (S_L_S(b) == NULL) /* AK 191093 */
2192 goto mmm_ende;
2193 if (S_L_S(a) == NULL) /* AK 220596 */
2194 goto mmm_ende;
2195
2196 b_ptr = b;
2197 if (EMPTYP(S_PO_S(b_ptr))) /* skip the empty term */
2198 b_ptr = S_L_N(b_ptr);
2199 while (b_ptr != NULL)
2200 {
2201 a_ptr = a;
2202 if (EMPTYP(S_PO_S(a_ptr))) /* skip the empty term */
2203 a_ptr = S_L_N(a_ptr);
2204 while (a_ptr != NULL)
2205 {
2206 temp_eins = CALLOCOBJECT();
2207 temp_zwei = CALLOCOBJECT();
2208 ADD(S_PO_S(a_ptr), S_PO_S(b_ptr), temp_eins);
2209 MULT(S_PO_K(a_ptr), S_PO_K(b_ptr), temp_zwei);
2210 new = CALLOCOBJECT();
2211 erg += b_sk_mo(temp_eins,temp_zwei,new);
2212 insert(new,c_tmp,add_koeff,NULL);
2213 a_ptr = S_L_N(a_ptr);
2214 }
2215 b_ptr = S_L_N(b_ptr);
2216 }
2217 erg += remove_zero_terms(c);
2218 mmm_ende:
2219 CTO(MONOPOLY,"mult_monopoly_monopoly(e3)",c);
2220 ENDR("mult_monopoly_monopoly");
2221 }
2222 #endif /* MONOPOLYTRUE */
2223
2224 /* This routine deals with MONOPOLYs which are effectively */
2225 /* LISTs of MONOMs whose coefficients are scalars and selfs */
2226 /* are non-negative integers -- they correspond to polynomials */
2227 /* in one variable. The quotient (qpoly) and remainder (rpoly) */
2228 /* of the division of poly by dpoly are calculated. The */
2229 /* parameters poly,dpoly,qpoly and rpoly must all be different.*/
2230
quores_monopoly_pre300402(poly,dpoly,qpoly,rpoly)2231 INT quores_monopoly_pre300402(poly,dpoly,qpoly,rpoly) OP poly,dpoly,qpoly,rpoly;
2232 /* 30.05.90: TPMcD. */ /* 1.04.91: TPMcD. */
2233 /* AK 200891 V1.3 */
2234 /* 04.10.91: TPMcD */
2235 {
2236 INT ret = ERROR,erg=OK;
2237
2238 OP term , temp, sz, ptr_a, coeff ,
2239 minus , ddeg , dlead ,
2240 rdeg , rlead ;
2241
2242 CTO(MONOPOLY,"quores_monopoly(1)",poly);
2243 CTO(MONOPOLY,"quores_monopoly(2)",dpoly);
2244 SYMCHECK(NULLP(dpoly),"quores_monopoly: divisor is zero");
2245
2246 term = NULL; sz=NULL;
2247 coeff = CALLOCOBJECT();
2248 minus = CALLOCOBJECT(); ddeg = CALLOCOBJECT(); dlead = CALLOCOBJECT();
2249 rdeg = CALLOCOBJECT(); rlead = CALLOCOBJECT();
2250 M_I_I(-1L,minus);
2251
2252 /* Find the leading term of dpoly */
2253 ptr_a = dpoly;
2254 while (ptr_a != NULL)
2255 {
2256 sz = ptr_a;
2257 ptr_a = S_L_N(ptr_a);
2258 }
2259 COPY(S_PO_S(sz),ddeg);
2260 COPY(S_PO_K(sz),dlead);
2261 if (NULLP(ddeg) && NULLP(dlead))
2262 error("quores_monopoly: divisor is zero");
2263 ADDINVERS_APPLY(ddeg);
2264 ADDINVERS_APPLY(dlead);
2265
2266 copy(poly,rpoly);
2267 init(MONOPOLY,qpoly);
2268 ptr_a = rpoly;
2269 while(ptr_a != NULL)
2270 { /* If the remainder is zero, break from the loop */
2271 if (empty_listp(rpoly))
2272 {
2273 insert_zero_into_monopoly(rpoly);
2274 ret = OK;
2275 break;
2276 }
2277 /* Find the leading term of the current rpoly */
2278 while (ptr_a != NULL)
2279 {
2280 sz = ptr_a;
2281 ptr_a = S_L_N(ptr_a);
2282 }
2283 copy(S_PO_S(sz),rdeg);
2284 copy(S_PO_K(sz),rlead);
2285 ADD_APPLY(ddeg,rdeg);
2286 /* The remainder has degree less than the divisor */
2287 if (LT(rdeg,cons_null))
2288 {
2289 ret = OK;
2290 goto exit_label;
2291 }
2292 div(rlead,dlead,coeff);
2293 temp = CALLOCOBJECT();
2294 init(MONOPOLY,temp);
2295 term = CALLOCOBJECT();
2296 m_sk_mo(rdeg,coeff,term);
2297 insert(term,temp,add_koeff,NULL);
2298
2299 term = CALLOCOBJECT();
2300 init(MONOPOLY,term);
2301 C_L_S(term,CALLOCOBJECT());
2302 m_sk_mo(rdeg,cons_eins,S_L_S(term));
2303 MULT_APPLY(coeff,term);
2304
2305 insert(term,qpoly,add_koeff,NULL);
2306 term = CALLOCOBJECT();
2307 mult_monopoly_monopoly(temp,dpoly,term);
2308 insert(term,rpoly,add_koeff,NULL);
2309 freeall(temp);
2310 ptr_a = rpoly;
2311 }
2312 exit_label:
2313 remove_zero_terms(rpoly);
2314 if (empty_listp(qpoly))
2315 insert_zero_into_monopoly(qpoly);
2316 mult_apply_scalar_monopoly(minus,qpoly);
2317 remove_zero_terms(qpoly);
2318 FREEALL(coeff); /* AK 080891 */
2319 FREEALL(minus);
2320 FREEALL(ddeg);
2321 FREEALL(rdeg);
2322 FREEALL(dlead);
2323 FREEALL(rlead);
2324 erg = ret;
2325 {
2326 OP c;
2327 c = CALLOCOBJECT();
2328 MULT(dpoly,qpoly,c);
2329 ADD_APPLY(rpoly,c);
2330 if (NEQ(c,poly)) {
2331 println(poly);
2332 println(dpoly);
2333 println(qpoly);
2334 println(rpoly);
2335 println(c);
2336 error("quores_monopoly:wrong result");
2337 }
2338 FREEALL(c);
2339 }
2340 ENDR("quores_monopoly");
2341 }
2342
2343
2344
quores_monopoly(a,b,c,d)2345 INT quores_monopoly(a,b,c,d) OP a,b,c,d;
2346 /* a/b = c rest d */
2347 /* with monopoly = univariate polynomial */
2348 /* AK 300402 */
2349 /* simple verion, not optimized */
2350 {
2351 INT erg = OK;
2352 OP degree_d,degree_b,leading_coeff_b,leading_coeff_d;
2353 OP h,z,ca;
2354 CTO(MONOPOLY,"quores_monopoly(1)",a);
2355 CTO(MONOPOLY,"quores_monopoly(2)",b);
2356 CTO(EMPTY,"quores_monopoly(3)",c);
2357 CTO(EMPTY,"quores_monopoly(4)",d);
2358 SYMCHECK(c==d,"quores_monopoly:should be different objects");
2359
2360 if (NULLP(b)) error("quores_monopoly(2):null");
2361
2362 if (NULLP(a)) {
2363 erg += copy(a,c);
2364 erg += copy(a,d);
2365 goto endr_ende;
2366 }
2367
2368 init(MONOPOLY,c);
2369 COPY(a,d);
2370
2371 CALLOCOBJECT6(degree_d,degree_b,leading_coeff_b,leading_coeff_d,h,ca);
2372
2373 degree_monopoly(b,degree_b);
2374 ldcf_monopoly(b,leading_coeff_b);
2375
2376 again:
2377 CTO(MONOPOLY,"quores_monopoly(i3)",c);
2378 CTO(MONOPOLY,"quores_monopoly(i4)",d);
2379 FREESELF(degree_d);
2380 degree_monopoly(d,degree_d);
2381 if (LT(degree_d,degree_b)) goto ee; /* finished */
2382
2383 /* else
2384 multiply b by leadcoeff_d/leadcoeff_b
2385 subtract this from d
2386 */
2387 FREESELF(leading_coeff_d);
2388 ldcf_monopoly(d,leading_coeff_d);
2389 m_iindex_monopoly( S_I_I(degree_d) - S_I_I(degree_b),ca);
2390 div(leading_coeff_d,leading_coeff_b,h);
2391 copy(h,S_PO_K(ca));
2392 mult_apply(b,h);
2393 FORALL(z,h,M_I_I(S_I_I(S_MO_S(z)) + S_I_I(degree_d) - S_I_I(degree_b),
2394 S_MO_S(z)));
2395 addinvers_apply(h);
2396 ADD_APPLY(h,d);
2397 ADD_APPLY(ca,c);
2398 goto again;
2399 ee:
2400 FREEALL6(degree_d,degree_b,leading_coeff_b,leading_coeff_d,h,ca);
2401 CTO(MONOPOLY,"quores_monopoly(e3)",c);
2402 CTO(MONOPOLY,"quores_monopoly(e4)",d);
2403 /*
2404 {
2405 OP e;
2406 e = CALLOCOBJECT();
2407 MULT(b,c,e);
2408 ADD_APPLY(d,e);
2409 if (NEQ(e,a)) {
2410 println(a);
2411 println(b);
2412 println(c);
2413 println(d);
2414 println(e);
2415 error("quores_monopoly:wrong result");
2416 }
2417 FREEALL(e);
2418 }
2419 */
2420 ENDR("quores_monopoly");
2421 }
2422
2423 #ifdef MONOPOLYTRUE
add_monopoly(a,b,c)2424 INT add_monopoly(a,b,c) OP a,b,c;
2425 /* AK 200891 V1.3 */
2426 {
2427 INT erg = OK;
2428 OP d;
2429 CTO(MONOPOLY,"add_monopoly(1)",a);
2430 CTO(EMPTY,"add_monopoly(3)",c);
2431
2432 switch(S_O_K(b))
2433 {
2434 case INTEGER:
2435 case FF:
2436 case LONGINT:
2437 erg += add_scalar_monopoly(b,a,c);
2438 goto ee;
2439 case BRUCH:
2440 erg += add_bruch(b,a,c);
2441 goto ee;
2442
2443 case MONOPOLY:
2444 erg += add_monopoly_monopoly(a,b,c);
2445 goto ee;
2446 case LAURENT:
2447 d = CALLOCOBJECT();
2448 erg += t_LAURENT_OBJ(b,d);
2449 erg += add(a,d,c);
2450 erg += freeall(d);
2451 goto ee;
2452 case POLYNOM:
2453 d=CALLOCOBJECT();
2454 erg += t_POLYNOM_MONOPOLY(b,d);
2455 erg += add_monopoly_monopoly(a,d,c);
2456 erg += freeall(d);
2457 goto ee;
2458
2459 default:
2460 WTO("add_monopoly(2)",b);
2461 goto ee;
2462 }
2463 ee:
2464 ENDR("add_monopoly");
2465 }
2466
2467
mult_monopoly_polynom(a,b,c)2468 INT mult_monopoly_polynom(a,b,c) OP a,b,c;
2469 /* CC because of laurent */
2470 {
2471 INT erg = OK;
2472 CTO(POLYNOM,"mult_monopoly_polynom(2)",b);
2473 CTO(MONOPOLY,"mult_monopoly_polynom(1)",a);
2474 CTO(EMPTY,"mult_monopoly_polynom(3)",c);
2475 if(has_one_variable(b)==TRUE) /* CC */
2476 {
2477 OP d;
2478 d=CALLOCOBJECT();
2479 erg += t_POLYNOM_MONOPOLY(b,d);
2480 erg += mult_monopoly_monopoly(a,d,c);
2481 FREEALL(d);
2482 }
2483 else
2484 erg += mult_scalar_polynom(a,b,c);
2485 ENDR("mult_monopoly_polynom");
2486 }
2487
2488
mult_monopoly(a,b,c)2489 INT mult_monopoly(a,b,c) OP a,b,c;
2490 /* 24.07.91: TPMcD. */
2491 /* AK 200891 V1.3 */
2492 {
2493 OP d;
2494 INT erg = OK;
2495 CTO(MONOPOLY,"mult_monopoly(1)",a);
2496 CTO(ANYTYPE,"mult_monopoly(2)",b);
2497 CTO(EMPTY,"mult_monopoly(3)",c);
2498
2499 switch(S_O_K(b))
2500 {
2501 case BRUCH:
2502 if (scalarp(S_B_O(b)) && scalarp(S_B_U(b)))
2503 erg += mult_scalar_monopoly(b,a,c);
2504 else
2505 {
2506 erg += copy(b,c);
2507 erg += mult(a,S_B_O(b),S_B_O(c));
2508 erg += kuerzen(c);
2509 }
2510 goto ee;
2511 case INTEGER:
2512 case FF:
2513 case LONGINT:
2514 erg += mult_scalar_monopoly(b,a,c);
2515 goto ee;
2516 case MONOPOLY:
2517 erg += mult_monopoly_monopoly(a,b,c);
2518 goto ee;
2519 #ifdef MATRIXTRUE
2520 case MATRIX:
2521 erg += mult_scalar_matrix(a,b,c);
2522 goto ee;
2523 #endif /* MATRIXTRUE */
2524 #ifdef MONOMTRUE
2525 case MONOM:
2526 erg += mult_scalar_monom(a,b,c);
2527 break;
2528 #endif /* MONOMTRUE */
2529 #ifdef LAURENTTRUE
2530 case LAURENT:
2531 d = CALLOCOBJECT();
2532 erg += t_LAURENT_OBJ(b,d);
2533 erg += mult(a,d,c);
2534 erg += freeall(d);
2535 break;
2536 #endif /* LAURENTTRUE */
2537 #ifdef POLYTRUE
2538 case POLYNOM:
2539 erg += mult_monopoly_polynom(a,b,c);
2540 break;
2541 #endif /* POLYTRUE */
2542 #ifdef SCHUBERTTRUE
2543 case SCHUBERT: erg += mult_scalar_schubert(a,b,c);
2544 break;
2545 #endif /* SCHUBERTTRUE */
2546 case VECTOR:
2547 erg += mult_scalar_vector(a,b,c);
2548 break;
2549 default:
2550 WTT("mult_monopoly",a,b);
2551 }
2552
2553 ee:
2554 CTO(ANYTYPE,"mult_monopoly(3)",c);
2555 ENDR("mult_monopoly");
2556 }
2557
addinvers_monopoly(a,b)2558 INT addinvers_monopoly(a,b) OP a,b;
2559 /* AK 200891 V1.3 */
2560 /* AK 200498 V2.0 */
2561 /* a and b may be equal */
2562 {
2563 INT erg = OK;
2564 CTO(MONOPOLY,"addinvers_monopoly(1)",a);
2565 CTO(EMPTY,"addinvers_monopoly(2)",b);
2566 erg += mult_scalar_monopoly(cons_negeins,a,b);
2567 ENDR("addinvers_monopoly");
2568 }
2569
2570
add_apply_monopoly(a,b)2571 INT add_apply_monopoly(a,b) OP a,b;
2572 /* AK 200891 V1.3 */
2573 {
2574 INT erg = OK;
2575 OP c;
2576 CTO(MONOPOLY,"add_apply_monopoly(1)",a);
2577 EOP("add_apply_monopoly(2)",b);
2578
2579 c = CALLOCOBJECT();
2580 SWAP(c,b);
2581 erg += add_monopoly(a,c,b);
2582 FREEALL(c);
2583 ENDR("add_apply_monopoly");
2584 }
2585
2586
mult_apply_monopoly(a,b)2587 INT mult_apply_monopoly(a,b) OP a,b;
2588 /* AK 200891 V1.3 */
2589 {
2590 INT erg = OK;
2591 OP c;
2592 CTO(MONOPOLY,"mult_apply_monopoly(1)",a);
2593 EOP("mult_apply_monopoly(2)",b);
2594
2595 c = CALLOCOBJECT();
2596 erg += mult_monopoly(a,b,c);
2597 erg += copy(c,b);
2598 erg += freeall(c);
2599
2600 ENDR("mult_apply_monopoly");
2601 }
2602
addinvers_apply_monopoly(a)2603 INT addinvers_apply_monopoly(a) OP a;
2604 /* AK 200891 V1.3 */
2605 /* AK 290402 */
2606 {
2607 INT erg = OK;
2608 OP z;
2609 CTO(MONOPOLY,"addinvers_apply_monopoly(1)",a);
2610 /* return(addinvers_monopoly(a,a)); */
2611 FORALL(z,a,ADDINVERS_APPLY(S_MO_K(z)));
2612 ENDR("addinvers_apply_monopoly");
2613 }
2614 #endif /* MONOPOLYTRUE */
2615
nullp_monopoly(a)2616 INT nullp_monopoly(a) OP a;
2617 /* 1.04.91: TPMcD. */
2618 /* AK 200891 V1.3 */
2619 {
2620 INT erg = OK;
2621 # ifdef MONOPOLYTRUE
2622 OP ptr = a;
2623 CTO(MONOPOLY,"nullp_monopoly",a);
2624
2625 if (empty_listp(a))
2626 return(TRUE);
2627 if (EMPTYP(S_L_S(ptr)))
2628 ptr = S_L_N(ptr);
2629 if (ptr == NULL || empty_listp(ptr))
2630 return(TRUE);
2631 if (S_L_N(ptr) != NULL)
2632 return(FALSE);
2633 if (EQ(S_PO_S(ptr),cons_null) && EQ(S_PO_K(ptr),cons_null))
2634 return(TRUE);
2635 #endif
2636 return(FALSE);
2637 ENDR("nullp_monopoly");
2638 }
2639
comp_monopoly(a,b)2640 INT comp_monopoly(a,b) OP a,b;
2641 /* AK 200891 V1.3 */
2642 /* 23.10.91: TPMcD */
2643 {
2644 INT erg = OK,ret;
2645 CTO(MONOPOLY,"comp_monopoly(1)",a);
2646 CTO(MONOPOLY,"comp_monopoly(2)",b);
2647 ret = comp_list(a,b);
2648 return ret;
2649 ENDR("comp_monopoly");
2650 }
2651
raise_power_monopoly(a,b)2652 INT raise_power_monopoly(a,b) OP a, b;
2653 /* 1.04.91: TPMcD. */
2654 /* AK 200891 V1.3 */
2655 {
2656 INT erg = OK;
2657 OP ptr;
2658 CTTO(INTEGER,LONGINT,"raise_power_monopoly(1)",a);
2659 CTO(MONOPOLY,"raise_power_monopoly(2)",b);
2660
2661 if (EINSP(a)) {
2662 }
2663 else if (POSP(a))
2664 {
2665 FORALL(ptr,b,{ MULT_APPLY(a,S_MO_S(ptr)); });
2666 }
2667 else {
2668 SYMCHECK(1,"raise_power_monopoly: exponent <= 0");
2669 }
2670 CTO(MONOPOLY,"raise_power_monopoly(e2)",b);
2671 ENDR("raise_power_monopoly");
2672 }
2673
scale_monopoly(a,b)2674 INT scale_monopoly(a,b) OP a, b;
2675 /* MD */
2676 /* AK 200891 V1.3 */
2677 {
2678 OP ptr;
2679 # ifdef MONOPOLYTRUE
2680 OP factor = CALLOCOBJECT(), minus = CALLOCOBJECT();
2681 M_I_I(-1L,minus);
2682 ptr = b;
2683 if (EMPTYP(S_PO_S(ptr)))
2684 ptr = S_L_N(ptr);
2685 if (EQ(a,minus) == TRUE)
2686 while (ptr != NULL)
2687 {
2688 mod(S_PO_S(ptr),cons_zwei,factor);
2689 if (EQ(factor,cons_eins) == TRUE)
2690 mult(minus,S_PO_K(ptr),S_PO_K(ptr));
2691 ptr = S_L_N(ptr);
2692 }
2693 else
2694 while (ptr != NULL)
2695 {
2696 hoch(a,S_PO_S(ptr),factor);
2697 mult(factor,S_PO_K(ptr),S_PO_K(ptr));
2698 ptr = S_L_N(ptr);
2699 }
2700 freeall(factor);
2701 freeall(minus);
2702 #endif
2703 return(OK);
2704 }
2705
objectread_monopoly(f,a)2706 INT objectread_monopoly(f,a) FILE *f; OP a;
2707 /* 4.04.91: TPMcD. */
2708 /* AK 200891 V1.3 */
2709 /* 23.10.91: TPMcD */
2710 {
2711 # ifdef MONOPOLYTRUE
2712 objectread_list(f,a);
2713 C_O_K(a,MONOPOLY);
2714 #endif
2715 return(OK);
2716 }
2717
2718 # ifdef MONOPOLYTRUE
tex_monopoly(a)2719 INT tex_monopoly(a) OP a;
2720 /* 2.04.91: TPMcD. */ /* AK 200891 V1.3 */ /* 04.10.91: TPMcD */
2721 {
2722 INT myfirst = 1L;
2723 OP ptr = a;
2724 fprintf(texout," ");
2725 while (ptr != NULL)
2726 {
2727 if (!negp(S_PO_K(ptr)) && !myfirst)
2728 fprintf(texout," + {");
2729 else
2730 fprintf(texout,"{");
2731 tex(S_PO_K(ptr));
2732 fprintf(texout,"} x {");
2733 tex(S_PO_S(ptr));
2734 fprintf(texout,"}\n");
2735 ptr = S_L_N(ptr);
2736 myfirst = 0L;
2737 texposition += 6;
2738 }
2739 fprintf(texout,"\n");
2740 texposition = 0;
2741 return(OK);
2742 }
2743
2744 /* CONSTRUCTION OF SPECIAL TYPES OF MONOPOLY */
2745
2746 /* Given the number n, which should be an positive INTEGER or LONGINT,*/
2747 /* the result returns the polynomial x**n-1. */
2748
m_iindex_monopoly(ii,mon)2749 INT m_iindex_monopoly(ii,mon) INT ii; OP mon;
2750 {
2751 INT erg = OK;
2752 COP("m_iindex_monopoly",mon);
2753 erg += b_skn_mp(CALLOCOBJECT(),CALLOCOBJECT(),NULL,mon);
2754 erg += m_i_i(ii,S_PO_S(mon));
2755 erg += m_i_i(1L,S_PO_K(mon));
2756 ENDR("m_iindex_monopoly");
2757 }
2758
make_unitary0_monopoly(number,result)2759 INT make_unitary0_monopoly(number,result) OP number, result;
2760 /* 10.06.90: TPMcD. */ /* AK 200891 V1.3 */
2761 {
2762 OP new;
2763 OP e = CALLOCOBJECT(), coeff = CALLOCOBJECT();
2764 init(MONOPOLY,result);
2765 M_I_I(0L,e); /* we use the macro, beacuse e is empty */
2766 M_I_I(-1L,coeff);
2767 new = CALLOCOBJECT();
2768 m_sk_mo(e,coeff,new);
2769 insert(new,result,add_koeff,NULL);
2770 m_i_i(1L,coeff);
2771 new = CALLOCOBJECT();
2772 m_sk_mo(number,coeff,new);
2773 insert(new,result,add_koeff,NULL);
2774 freeall(e);
2775 freeall(coeff);
2776 return(OK);
2777 }
2778 #endif /* MONOPOLYTRUE */
2779
2780 /* Given the number n, which should be an positive INTEGER or LONGINT,*/
2781 /* the result returns the MONOPOLY x**(n-1)+x**(n-2)+...+x+1. */
2782
make_unitary_eins_monopoly(number,result)2783 INT make_unitary_eins_monopoly(number,result) OP number, result;
2784 /* 10.06.90: TPMcD. */ /* AK 200891 V1.3 */
2785 {
2786 OP new, ptr;
2787 INT erg = OK;
2788 # ifdef MONOPOLYTRUE
2789 OP e = CALLOCOBJECT(), coeff = CALLOCOBJECT();
2790 init(MONOPOLY,result);
2791 ptr = result;
2792 M_I_I(0L,e);
2793 M_I_I(1L,coeff);
2794 while (LT(e,number) == TRUE)
2795 {
2796 new = CALLOCOBJECT();
2797 init(MONOPOLY,new);
2798 S_L_S(new) = CALLOCOBJECT();
2799 m_sk_mo(e,coeff,S_L_S(new));
2800 insert(new,result,NULL,NULL);
2801 INC(e);
2802 }
2803 freeall(e);
2804 freeall(coeff);
2805 return(OK);
2806 #else
2807 return(ERROR);
2808 #endif
2809 }
2810
2811 /* Given the number n, which should be an positive INTEGER or LONGINT*/
2812 /* or a MONOPOLY representing a prime factorisation of an integer greater*/
2813 /* than 1 , the result returns the cyclotomic polynomial of index n.*/
2814
make_cyclotomic_monopoly(number,result)2815 INT make_cyclotomic_monopoly(number,result) OP number, result;
2816 /* 100690: TPMcD. */ /* AK 200891 V1.3 */
2817 {
2818 INT ret = ERROR;
2819 OP ptr;
2820 # ifdef MONOPOLYTRUE
2821 OP minus, e, factor, poly_eins, poly_zwei, poly_drei,
2822 list = CALLOCOBJECT();
2823
2824 if (S_O_K(number) == MONOPOLY)
2825 copy(number,list);
2826 else
2827 {
2828 if (((S_O_K(number) != INTEGER) && (S_O_K(number) != LONGINT))
2829 || (LT(number,cons_eins) == TRUE))
2830 goto exit_label;
2831
2832 if (einsp(number))
2833 {
2834 make_unitary0_monopoly(number,result);
2835 ret = OK;
2836 goto exit_label;
2837 }
2838 integer_factor(number,list);
2839 }
2840 ptr = list;
2841 if (EMPTYP(S_PO_S(ptr)))
2842 ptr = S_L_N(ptr);
2843 factor = CALLOCOBJECT();
2844 copy(S_PO_S(ptr),factor);
2845
2846 if (S_L_N(ptr) == NULL) /* list has just one prime factor */
2847 {
2848 make_unitary_eins_monopoly(factor,result);
2849 dec(S_PO_K(ptr));
2850 e = CALLOCOBJECT();
2851 integer_factors_to_integer(ptr,e);
2852 raise_power_monopoly(e,result);
2853 freeall(e);
2854 ret = OK;
2855 freeall(factor); /* AK 060993 */
2856 goto exit_label;
2857 }
2858 else
2859 {
2860 poly_eins = CALLOCOBJECT();
2861 ptr = S_L_N(ptr);
2862 make_cyclotomic_monopoly(ptr,poly_eins);
2863 if (EQ(factor,cons_zwei) == TRUE)
2864 {
2865 minus = CALLOCOBJECT();
2866 M_I_I(-1L,minus);
2867 scale_monopoly(minus,poly_eins);
2868 copy(poly_eins,result);
2869 freeall(minus);
2870 }
2871 else
2872 {
2873 poly_zwei = CALLOCOBJECT();
2874 poly_drei = CALLOCOBJECT();
2875 copy(poly_eins,poly_zwei);
2876 raise_power_monopoly(factor,poly_eins);
2877 quores_monopoly(poly_eins,poly_zwei,result,poly_drei);
2878 freeall(poly_zwei);
2879 freeall(poly_drei);
2880 }
2881 freeall(poly_eins);
2882 ret = OK;
2883 }
2884 ptr = list;
2885 if (EMPTYP(S_PO_S(ptr)))
2886 /* ptr = S_L_N(ptr); */
2887 error("nb internal error 234");
2888 dec(S_PO_K(ptr));
2889 ptr = S_L_N(ptr);
2890 while (ptr != NULL)
2891 {
2892 /* dec(S_PO_K(ptr)); */
2893 m_i_i(0L,S_PO_K(ptr));
2894 ptr = S_L_N(ptr);
2895 }
2896 integer_factors_to_integer(list,factor);
2897 raise_power_monopoly(factor,result);
2898 freeall(factor);
2899 exit_label:
2900 freeall(list);
2901 #endif
2902 return(ret);
2903 }
2904 /****************** fields_2.c **********************/
2905 /* 26.07.91: TPMcD. */
2906 /*************************************************************/
2907
2908 /* SQUARE RADICALS */
2909
2910 #ifdef SQRADTRUE
make_monopoly_sqrad(a,b)2911 INT make_monopoly_sqrad(a,b) OP a,b;
2912 /* 30.05.90: TPMcD. */ /* 3.04.91: TPMcD. */
2913 /* AK 200891 V1.3 */
2914 {
2915 INT flag = 0L;
2916 INT erg = OK;
2917 OP ptr, new, b_tmp, sqfree , sqpart ;
2918
2919 CTO(MONOPOLY,"make_monopoly_sqrad(1)",a);
2920 if (b == a)
2921 {
2922 flag = 1L;
2923 b_tmp = CALLOCOBJECT();
2924 }
2925 else
2926 {
2927 init(SQ_RADICAL,b);
2928 b_tmp = S_N_S(b);
2929 }
2930 sqfree = CALLOCOBJECT();
2931 sqpart = CALLOCOBJECT();
2932 init(MONOPOLY, b_tmp);
2933 ptr = a;
2934 while (ptr != NULL)
2935 {
2936 if (not nullp(S_PO_S(ptr))) /* AK 120891 */
2937 {
2938 square_free_part(S_PO_S(ptr),sqfree,
2939 sqpart,NULL,NULL,NULL);
2940 MULT_APPLY(S_PO_K(ptr),sqpart);
2941 new = CALLOCOBJECT();
2942 m_sk_mo(sqfree,sqpart,new);
2943 insert(new,b_tmp,add_koeff,
2944 comp_monomvector_monomvector);
2945 }
2946 ptr = S_L_N(ptr);
2947 }
2948 if (flag)
2949 {
2950 init(SQ_RADICAL,b);
2951 freeall(S_N_S(b)); /* AK 060993 */
2952 S_N_S(b) = b_tmp;
2953 }
2954 remove_zero_terms(S_N_S(b));
2955 find_sqrad_data(b);
2956 freeall(sqfree);
2957 freeall(sqpart);
2958 ENDR("make_monopoly_sqrad");
2959 }
2960 #endif /* SQRADTRUE */
2961 /* a: the scalar; b: the result. */
2962
2963 #ifdef SQRADTRUE
make_scalar_sqrad(a,b)2964 INT make_scalar_sqrad(a,b) OP a,b;
2965 /* 5.04.91: TPMcD. */ /* AK 200891 V1.3 */ /* 04.10.91: TPMcD */
2966 {
2967 OP c;
2968 INT erg = OK;
2969 EOP("make_scalar_sqrad(1)",a);
2970
2971 c = CALLOCOBJECT();
2972 erg += b_skn_mp(CALLOCOBJECT(),CALLOCOBJECT(),NULL,c);
2973 erg += copy(a,S_PO_K(c));
2974 M_I_I(1L,S_PO_S(c));
2975 erg += make_monopoly_sqrad(c,b);
2976 if (not EMPTYP(S_N_D(b)))
2977 erg += freeself(S_N_D(b));
2978 erg += find_sqrad_data(b);
2979 erg += freeall(c);
2980
2981 ENDR("make_scalar_sqrad");
2982 }
2983 #endif
2984
2985 #ifdef SQRADTRUE
scan_sqrad(a)2986 INT scan_sqrad(a) OP a;
2987 /* AK 300191 V1.2 */ /* a becomes sqrad */
2988 /* 3.04.91: TPMcD. */ /* AK 200891 V1.3 */ /* 04.10.91: TPMcD */
2989 {
2990 OP c = CALLOCOBJECT();
2991 INT erg=OK;
2992 erg += printeingabe("self of sqrad");
2993 erg += scan(MONOPOLY,c);
2994 erg += make_monopoly_sqrad(c,a);
2995 erg += find_sqrad_data(a);
2996 erg += freeall(c);
2997 return erg;
2998 }
2999 #endif /* SQRADTRUE */
3000
3001 #ifdef SQRADTRUE
add_scalar_sqrad(a,b,c)3002 INT add_scalar_sqrad(a,b,c) OP a,b,c;
3003 /* 30.05.90: TPMcD. */ /* AK 200891 V1.3 */ /* 04.10.91: TPMcD */
3004 {
3005 OP ptr;
3006 INT erg = OK;
3007 CTO(SQ_RADICAL,"add_scalar_sqrad(2)",b);
3008 CTO(EMPTY,"add_scalar_sqrad(3)",c);
3009
3010 erg += copy(b,c);
3011 ptr = CALLOCOBJECT();
3012 erg += init(MONOPOLY,ptr);
3013 C_L_S(ptr,CALLOCOBJECT());
3014 erg += m_sk_mo(cons_eins,a,S_L_S(ptr));
3015 ADD_APPLY(ptr,S_N_S(c));
3016 erg += freeall(ptr);
3017 if (space_saving)
3018 convert_sqrad_scalar(c);
3019
3020 ENDR("add_scalar_sqrad");
3021 }
3022
3023
3024
mult_apply_scalar_sqrad(a,b)3025 INT mult_apply_scalar_sqrad(a,b) OP a,b;
3026 /* AK 060498 V2.0 */
3027 {
3028 OP c;
3029 INT erg = OK;
3030 CTO(SQ_RADICAL,"mult_apply_scalar_sqrad(2)",b);
3031 c = CALLOCOBJECT();
3032 SWAP(c,b);
3033 erg += mult_scalar_sqrad(a,c,b);
3034 FREEALL(c);
3035 ENDR("mult_apply_scalar_sqrad");
3036 }
3037
mult_scalar_sqrad(a,b,c)3038 INT mult_scalar_sqrad(a,b,c) OP a, b, c;
3039 /* 30.05.90: TPMcD. */ /* AK 200891 V1.3 */ /* 04.10.91: TPMcD */
3040 {
3041 OP ptr;
3042 INT erg = OK;
3043 CTO(SQ_RADICAL,"mult_scalar_sqrad(2)",b);
3044 CTTO(INTEGER,EMPTY,"mult_scalar_sqrad(3)",c);
3045
3046 if (S_O_K(c)==INTEGER) C_O_K(c,EMPTY);
3047
3048 if (space_saving && nullp(a))
3049 {
3050 COPY(a,c);
3051 goto endr_ende;
3052 }
3053
3054 COPY(b,c);
3055 ptr = S_N_S(c);
3056 if (EMPTYP(S_PO_S(ptr)))
3057 ptr = S_L_N(ptr);
3058 while (ptr != NULL)
3059 {
3060
3061 MULT_APPLY(a,S_PO_K(ptr));
3062 ptr = S_L_N(ptr);
3063 }
3064
3065 remove_zero_terms(S_N_S(c));
3066 if (nullp_sqrad(c))
3067 find_sqrad_data(c);
3068
3069 ENDR("mult_scalar_sqrad");
3070 }
3071
add_sqrad_sqrad(a,b,c)3072 INT add_sqrad_sqrad(a,b,c) OP a, b, c;
3073 /* 23.06.90: TPMcD. */ /* 3.04.91: TPMcD. */
3074 /* AK 200891 V1.3 */ /* 04.10.91: TPMcD */
3075 /* AK 060498 V2.0 */
3076 {
3077 INT flag = 0L;
3078 INT erg = OK;
3079 OP c_tmp, data_a ,data_b;
3080 CTO(SQ_RADICAL,"add_sqrad_sqrad(1)",a);
3081 CTO(SQ_RADICAL,"add_sqrad_sqrad(2)",b);
3082 CTO(EMPTY,"add_sqrad_sqrad(3)",c);
3083
3084 data_a = CALLOCOBJECT();
3085 data_b = CALLOCOBJECT();
3086 find_sqrad_data(a);
3087 find_sqrad_data(b);
3088 copy(S_N_D(a),data_a);
3089 copy(S_N_D(b),data_b);
3090 if (!empty_listp(data_b))
3091 insert(data_b,data_a,NULL,NULL);
3092 else
3093 freeall(data_b);
3094 if (c == a || c == b)
3095 {
3096 flag = 1L;
3097 c_tmp = CALLOCOBJECT();
3098 }
3099 else
3100 {
3101 init(SQ_RADICAL,c);
3102 c_tmp = S_N_S(c);
3103 }
3104 /* erg += init(MONOPOLY,c_tmp); */
3105 FREESELF(c_tmp);
3106 erg += add_monopoly_monopoly(S_N_S(a),S_N_S(b),c_tmp);
3107 if (flag)
3108 {
3109 if (not EMPTYP(S_N_S(c)))
3110 erg += freeall(S_N_S(c));
3111 S_N_S(c) = c_tmp;
3112 }
3113 if (not EMPTYP(S_N_D(c)))
3114 erg += freeall(S_N_D(c));
3115 S_N_D(c) = data_a;
3116 adjust_sqrad_data(c);
3117 ENDR("add_sqrad_sqrad");
3118 }
3119
3120 #endif /* SQRADTRUE */
mult_sqrad_sqrad(a,b,c)3121 INT mult_sqrad_sqrad(a,b,c) OP a, b, c;
3122 /* 30.05.90: TPMcD. */ /* 3.04.91: TPMcD. */
3123 /* AK 200891 V1.3 */ /* 04.10.91: TPMcD */
3124 {
3125 INT flag = 0L;
3126 INT erg = OK;
3127
3128 OP a_ptr, b_ptr, c_tmp, new;
3129 OP temp_eins, temp_zwei, temp_drei, data_a, data_b;
3130 CTO(SQ_RADICAL,"mult_sqrad_sqrad(1)",a);
3131 CTO(SQ_RADICAL,"mult_sqrad_sqrad(2)",b);
3132
3133 find_sqrad_data(a);
3134 find_sqrad_data(b);
3135 data_a = CALLOCOBJECT();
3136 data_b = CALLOCOBJECT();
3137 COPY(S_N_D(a),data_a);
3138 COPY(S_N_D(b),data_b);
3139 if (!empty_listp(data_b))
3140 insert(data_b,data_a,NULL,NULL);
3141 else
3142 freeall(data_b);
3143 if (c == a || c == b)
3144 {
3145 flag = 1L;
3146 c_tmp = CALLOCOBJECT();
3147 }
3148 else
3149 {
3150 init(SQ_RADICAL,c);
3151 c_tmp = S_N_S(c);
3152 }
3153 init(MONOPOLY,c_tmp);
3154 b_ptr = S_N_S(b);
3155 if (EMPTYP(S_PO_S(b_ptr))) /* skip the empty term */
3156 b_ptr = S_L_N(b_ptr);
3157 temp_eins = CALLOCOBJECT();
3158 temp_zwei = CALLOCOBJECT();
3159 temp_drei = CALLOCOBJECT();
3160 while (b_ptr != NULL)
3161 {
3162 if (not nullp(S_PO_S(b_ptr))) /* AK 120891 */
3163 {
3164 a_ptr = S_N_S(a);
3165 while (a_ptr != NULL)
3166 {
3167 if (not nullp(S_PO_S(a_ptr))) /* AK 120891 */
3168 {
3169 ggt(S_PO_S(a_ptr), S_PO_S(b_ptr), temp_eins);
3170 if (NEGP(S_PO_S(a_ptr)) && NEGP( S_PO_S(b_ptr) ) )
3171 ADDINVERS_APPLY(temp_eins); /* AK 259492 */
3172 ganzdiv(S_PO_S(a_ptr), temp_eins, temp_zwei);
3173 ganzdiv(S_PO_S(b_ptr), temp_eins, temp_drei);
3174 MULT_APPLY(S_PO_K(a_ptr),temp_eins);
3175 MULT_APPLY(S_PO_K(b_ptr),temp_eins);
3176 MULT_APPLY(temp_drei,temp_zwei);
3177 new = CALLOCOBJECT();
3178 m_sk_mo(temp_zwei,temp_eins,new);
3179 insert(new,c_tmp,add_koeff,NULL);
3180 }
3181 a_ptr = S_L_N(a_ptr);
3182 }
3183 }
3184 b_ptr = S_L_N(b_ptr);
3185 }
3186 if (empty_listp(c_tmp))
3187 insert_zero_into_monopoly(c_tmp);
3188 if (flag)
3189 {
3190 init(SQ_RADICAL,c);
3191 freeall(S_N_S(c)); /* AK 060993 */
3192 S_N_S(c) = c_tmp;
3193 }
3194 if (not EMPTYP(S_N_D(c)))
3195 freeall(S_N_D(c));
3196 S_N_D(c) = data_a;
3197 adjust_sqrad_data(c);
3198 freeall(temp_eins);
3199 freeall(temp_zwei);
3200 freeall(temp_drei);
3201 return(OK);
3202 ENDR("mult_sqrad_sqrad");
3203 }
3204
3205 #ifdef SQRADTRUE
add_sqrad(a,b,c)3206 INT add_sqrad(a,b,c) OP a,b,c;
3207 /* AK 070891 V1.3 */
3208 /* 04.10.91: TPMcD */
3209 {
3210 INT erg = OK;
3211 CTO(SQ_RADICAL,"add_sqrad(1)",a);
3212 CTO(EMPTY,"add_sqrad(3)",c);
3213
3214 switch(S_O_K(b))
3215 {
3216 case CYCLOTOMIC: /* SQ + CYC = CYC */
3217 erg += add_cyclo(b,a,c);
3218 break;
3219 case INTEGER:
3220 case LONGINT:
3221 case BRUCH:
3222 erg += add_scalar_sqrad(b,a,c);
3223 break;
3224 case POLYNOM:
3225 erg += add_scalar_polynom(a,b,c);
3226 break;
3227 case SQ_RADICAL:
3228 erg += add_sqrad_sqrad(a,b,c);
3229 break;
3230 case EMPTY:
3231 erg += copy_number(a,c);
3232 break;
3233 default:
3234 WTO("add_sqrad(2)",b);
3235 goto ende;
3236 };
3237 if (space_saving)
3238 convert_sqrad_scalar(c);
3239
3240 ende:
3241 ENDR("add_sqrad");
3242 }
3243 #endif /* SQRADDTRUE */
3244
mult_sqrad(a,b,c)3245 INT mult_sqrad(a,b,c) OP a,b,c;
3246 /* AK 200291 V1.2 */ /* 24.07.91: TPMcD. */
3247 /* typ a ist SQ_RADICAL */
3248 /* AK 130891 V1.3 */
3249 /* 04.10.91: TPMcD */
3250 /* AK 060498 V2.0 */
3251 {
3252 INT erg = OK;
3253 CTO(SQ_RADICAL,"mult_sqrad(1)",a);
3254 if (S_O_K(c)!=EMPTY) FREESELF(c);
3255
3256
3257 switch(S_O_K(b))
3258 {
3259 case INTEGER:
3260 case LONGINT:
3261 case CYCLOTOMIC:
3262 case BRUCH:
3263 erg += mult_scalar_sqrad(b,a,c);
3264 goto ende;;
3265
3266 #ifdef MATRIXTRUE
3267 case MATRIX:
3268 erg += mult_scalar_matrix(a,b,c);
3269 break;
3270 #endif /* MATRIXTRUE */
3271
3272 #ifdef MONOMTRUE
3273 case MONOM: /* AK 200291 */
3274 erg += mult_scalar_monom(a,b,c);
3275 break;
3276 #endif /* MONOMTRUE */
3277
3278 case SQ_RADICAL:
3279 erg += mult_sqrad_sqrad(a,b,c);
3280 break;
3281
3282 #ifdef POLYTRUE
3283 case POLYNOM: /* AK 200291 */
3284 erg += mult_scalar_polynom(a,b,c);
3285 break;
3286 #endif /* POLYTRUE */
3287
3288 #ifdef SCHUBERTTRUE
3289 case SCHUBERT: /* AK 200291 */
3290 erg += mult_scalar_schubert(a,b,c);
3291 break;
3292 #endif /* SCHUBERTTRUE */
3293
3294 case VECTOR:
3295 erg += mult_scalar_vector(a,b,c);
3296 break;
3297
3298 default:
3299 WTO("mult_sqrad(2)",b);
3300 break;
3301 }
3302 if (space_saving)
3303 convert_sqrad_scalar(c);
3304 ende:
3305 ENDR("mult_sqrad");
3306 }
3307
3308 #ifdef MONOPOLYTRUE
addinvers_sqrad(a,b)3309 INT addinvers_sqrad(a,b) OP a,b;
3310 /* AK 200891 V1.3 */
3311 {
3312 if (S_O_K(a) != SQ_RADICAL)
3313 return error("addinvers_sqrad:object not sqrad");
3314 find_sqrad_data(a);
3315 mult_scalar_sqrad(cons_negeins,a,b);
3316 return(OK);
3317 }
3318 #endif /* MONOPOLYTRUE */
3319
3320 #ifdef SQRADTRUE
invers_sqrad(a,b)3321 INT invers_sqrad(a,b) OP a,b;
3322 /* 23.06.90: TPMcD. */ /* AK 200891 V1.3 */ /* 23.10.91: TPMcD */
3323 {
3324 INT ret = OK;
3325 INT flag = 0L;
3326 OP b_tmp=NULL, y_tmp, new, ptr, prime, x_tmp = CALLOCOBJECT(),
3327 data_length = CALLOCOBJECT(), mono_length = CALLOCOBJECT();
3328
3329 if (S_O_K(a) != SQ_RADICAL)
3330 {
3331 ret += invers(a,b); /* try the general routine */
3332 goto exit_label;
3333 }
3334 ret += find_sqrad_data(a);
3335 if (nullp_sqrad(a))
3336 error("invers_sqrad: 0 has no inverse\n");
3337 if (b == a)
3338 {
3339 b_tmp = CALLOCOBJECT();
3340 flag = 1L;
3341 }
3342 else
3343 b_tmp = b;
3344 ret += init(SQ_RADICAL,b_tmp);
3345 ret += init(MONOPOLY,S_N_S(b_tmp));
3346 ret += length(S_N_D(a),data_length);
3347 ret += length(S_N_S(a),mono_length);
3348 if (nullp(data_length)) /* No radicals */
3349 {
3350 ptr = S_N_S(a);
3351 ret += invers(S_PO_K(ptr),x_tmp);
3352 new = CALLOCOBJECT();
3353 ret += m_sk_mo(cons_eins,x_tmp,new);
3354 insert_list(new,S_N_S(b_tmp),NULL,NULL);
3355 goto exit_label;
3356 }
3357 else if (einsp(mono_length)) /* One radical */
3358 {
3359 ptr = S_N_S(a);
3360 mult(S_PO_S(ptr),S_PO_K(ptr),x_tmp);
3361 invers(x_tmp,x_tmp);
3362 new = CALLOCOBJECT();
3363 m_sk_mo(S_PO_S(ptr),x_tmp,new);
3364 insert_list(new,S_N_S(b_tmp),NULL,NULL);
3365 ret = OK;
3366 goto exit_label;
3367 }
3368 /* more than one radical */
3369 /* the conjugate is built up in b_tmp */
3370
3371 y_tmp = CALLOCOBJECT();
3372 copy(a,x_tmp);
3373 make_scalar_sqrad(cons_eins,b_tmp);
3374 ptr = S_N_D(a);
3375 while (ptr != NULL)
3376 {
3377 prime = S_L_S(ptr);
3378 if (S_O_K(x_tmp) != SQ_RADICAL)
3379 make_scalar_sqrad(x_tmp,x_tmp);
3380 conj_sqrad(x_tmp,prime,y_tmp);
3381 if (comp_sqrad(x_tmp,y_tmp) != 0L)
3382 {
3383 mult_sqrad_sqrad(x_tmp,y_tmp,x_tmp);
3384 mult_sqrad_sqrad(b_tmp,y_tmp,b_tmp);
3385 }
3386 ptr = S_L_N(ptr);
3387 }
3388 /* at this point x_tmp should be scalar */
3389 if (convert_sqrad_scalar(x_tmp) == ERROR)
3390 {
3391 freeall(y_tmp);
3392 error("invers_sqrad: the norm is not a scalar\n");
3393 goto exit_label;
3394 }
3395 if (negp(x_tmp))
3396 {
3397 ret += mult_apply_scalar_sqrad(cons_negeins,b_tmp);
3398 ret += addinvers_apply(x_tmp);
3399 }
3400
3401 ret += invers(x_tmp,y_tmp);
3402 ret += mult_apply_scalar_sqrad(y_tmp,b_tmp);
3403
3404 ret += freeall(y_tmp);
3405 exit_label:
3406 if (flag)
3407 {
3408 copy(b_tmp,b);
3409 freeall(b_tmp);
3410 }
3411 freeall(x_tmp);
3412 freeall(data_length);
3413 freeall(mono_length);
3414 return(ret);
3415 }
3416
3417
3418
add_apply_sqrad(a,b)3419 INT add_apply_sqrad(a,b) OP a,b;
3420 /* AK 200891 V1.3 */ /* 23.10.91: TPMcD */
3421 {
3422 INT erg=OK;
3423 OP c = CALLOCOBJECT();
3424 erg += add_sqrad(a,b,c);
3425 erg += copy(c,b);
3426 erg += freeall(c);
3427 return erg;
3428 }
3429 #endif /* SQRADTRUE */
3430
mult_apply_sqrad(a,b)3431 INT mult_apply_sqrad(a,b) OP a,b;
3432 /* AK 200891 V1.3 */
3433 /* 23.10.91: TPMcD */
3434 {
3435 INT ret = ERROR;
3436 #ifdef SQRADTRUE
3437 OP c;
3438 c = CALLOCOBJECT();
3439 ret = mult_sqrad(a,b,c);
3440 copy(c,b);
3441 freeall(c);
3442 #endif /* SQRADTRUE */
3443 return(ret);
3444 }
3445
3446 #ifdef SQRADTRUE
addinvers_apply_sqrad(a)3447 INT addinvers_apply_sqrad(a) OP a;
3448 /* AK 200891 V1.3 */
3449 {
3450 INT erg = OK;
3451 OP b;
3452 CTO(SQ_RADICAL,"addinvers_apply_sqrad(1)",a);
3453 b = CALLOCOBJECT();
3454 SWAP(a,b);
3455 erg += addinvers_sqrad(b,a);
3456 FREEALL(b);
3457 ENDR("addinvers_apply_sqrad");
3458 }
3459
3460
3461
nullp_sqrad(a)3462 INT nullp_sqrad(a) OP a;
3463 /* 1.04.91: TPMcD. */
3464 /* AK 200891 V1.3 */
3465 {
3466 INT erg = OK;
3467 OP ptr;
3468 CTO(SQ_RADICAL,"nullp_sqrad(1)",a);
3469 ptr = S_N_S(a);
3470 if (nullp_monopoly(ptr))
3471 return(TRUE);
3472 return(FALSE);
3473 ENDR("nullp_sqrad");
3474 }
3475
einsp_sqrad(a)3476 INT einsp_sqrad(a) OP a;
3477 /* AK 230402 */
3478 {
3479 INT erg = OK;
3480 OP ptr;
3481 CTO(SQ_RADICAL,"einsp_sqrad(1)",a);
3482 ptr = S_N_S(a);
3483 if (ptr == NULL) return FALSE;
3484 if (S_L_S(ptr) == NULL) return FALSE;
3485 if (not EINSP(S_PO_S(ptr))) return FALSE;
3486 if (not EINSP(S_PO_K(ptr))) return FALSE;
3487 if (S_PO_N(ptr) != NULL) return FALSE;
3488 return TRUE;
3489
3490 ENDR("einsp_sqrad");
3491 }
3492 #endif /* SQRADTRUE */
3493
3494 #ifdef CYCLOTRUE
einsp_cyclotomic(a)3495 INT einsp_cyclotomic(a) OP a;
3496 /* AK 240402 */
3497 {
3498 INT erg = OK;
3499 OP ptr;
3500 CTO(CYCLOTOMIC,"einsp_cyclotomic(1)",a);
3501 ptr = S_N_S(a);
3502 if (ptr == NULL) return FALSE;
3503 if (S_L_S(ptr) == NULL) return FALSE;
3504 if (not EINSP(S_PO_S(ptr))) return FALSE;
3505 if (not EINSP(S_PO_K(ptr))) return FALSE;
3506 if (S_PO_N(ptr) != NULL) return FALSE;
3507 return TRUE;
3508
3509 ENDR("einsp_cyclotomic");
3510 }
3511 #endif /*CYCLOTRUE */
3512
comp_sqrad(a,b)3513 INT comp_sqrad(a,b) OP a,b;
3514 /* AK 200891 V1.3 */
3515 {
3516 #ifdef SQRADTRUE
3517 return(comp_list(S_N_S(a),S_N_S(b)));
3518 #endif /* SQRADTRUE */
3519 }
3520
3521 #ifdef SQRADTRUE
fprint_sqrad(f,a)3522 static INT fprint_sqrad(f,a) FILE *f;
3523 OP a;
3524 /* 25.09.91: TPMcD. */
3525 {
3526 INT myfirst = 1L, rational = 0L;
3527 OP ptr;
3528 ptr = S_N_S(a);
3529 zeilenposition += 4L;
3530 if (nullp_sqrad(a))
3531 {
3532 fprintf(f," 0");
3533 return(OK);
3534 }
3535 while (ptr != NULL)
3536 {
3537 if (zeilenposition > 60L)
3538 {
3539 zeilenposition = 0L;
3540 fprintf(f,"\n");
3541 }
3542 if (einsp(S_PO_S(ptr)))
3543 rational = 1L; /* A rational term */
3544 else
3545 rational = 0L;
3546 /* print the coefficient part of a term */
3547 if (!negp(S_PO_K(ptr)) && !myfirst)
3548 fprintf(f," +");
3549 if (negeinsp(S_PO_K(ptr)))
3550 {
3551 if (rational)
3552 fprintf(f," - 1");
3553 else
3554 fprintf(f," -");
3555 }
3556 else if (einsp(S_PO_K(ptr)))
3557 {
3558 if (rational)
3559 fprintf(f," 1");
3560 }
3561 else
3562 {
3563 fprintf(f," ");
3564 fprint(f,S_PO_K(ptr));
3565 }
3566
3567 if (not rational) /* print the radical part of a term */
3568 {
3569 fprintf(f," sqr(");
3570 fprint(f,S_PO_S(ptr));
3571 fprintf(f,")");
3572 zeilenposition += 6L;
3573 }
3574 ptr = S_L_N(ptr);
3575 myfirst = 0L;
3576 }
3577 return(OK);
3578 }
3579 #endif /* SQRADTRUE */
3580
tex_sqrad(a)3581 INT tex_sqrad(a) OP a;
3582 /* 020491: TPMcD. */
3583 /* AK 200891 V1.3 */
3584 /* 041091: TPMcD */
3585 {
3586 INT myfirst = 1L;
3587 # ifdef SQRADTRUE
3588 OP ptr = S_N_S(a);
3589 find_sqrad_data(a);
3590 if (nullp_sqrad(a))
3591 {
3592 fprintf(texout," 0\n");
3593 return(OK);
3594 }
3595 fprintf(texout," ");
3596 while (ptr != NULL)
3597 {
3598 if (!negp(S_PO_K(ptr)) && !myfirst)
3599 fprintf(texout," + {");
3600 else
3601 fprintf(texout,"{");
3602 tex(S_PO_K(ptr));
3603 if (NEQ(S_PO_S(ptr),cons_eins))
3604 {
3605 fprintf(texout,"} \\sqrt{");
3606 tex(S_PO_S(ptr));
3607 }
3608 fprintf(texout,"}");
3609 ptr = S_L_N(ptr);
3610 myfirst = 0L;
3611 }
3612 fprintf(texout," ");
3613 return(OK);
3614 #else
3615 error("tex_sqrad: SQ_RADICAL not available");
3616 return(ERROR);
3617 #endif
3618 }
3619
find_sqrad_data(a)3620 static INT find_sqrad_data(a) OP a;
3621 /* 23.06.90: TPMcD. */ /* AK 200891 V1.3 */ /* 04.10.91: TPMcD */
3622 {
3623 #ifdef SQRADTRUE
3624 OP new, num, next_num, ptr, next_ptr, data_ptr, list_ptr,
3625 list_copy = CALLOCOBJECT(), prime_list = CALLOCOBJECT(),
3626 quo = CALLOCOBJECT(), rem = CALLOCOBJECT();
3627
3628 if (S_N_D(a) == NULL)
3629 S_N_D(a) = CALLOCOBJECT();
3630 data_ptr = S_N_D(a);
3631 /* Assume the data is OK if it is a non-empty LIST */
3632 if (not EMPTYP(data_ptr) && S_O_K(data_ptr) == LIST
3633 && not empty_listp(data_ptr))
3634 {
3635 goto fsd_ende;
3636 }
3637 init(LIST,data_ptr);
3638 copy(S_N_S(a),list_copy);
3639 ptr = list_copy;
3640 num = S_PO_S(ptr);
3641 if (LT(num,cons_null) == TRUE) /* negative radicals */
3642 {
3643 new = CALLOCOBJECT();
3644 M_I_I(-1L,new);
3645 insert_list(new,data_ptr,NULL,NULL);
3646 while (ptr != NULL) /*multiply negative radicals by -1*/
3647 {
3648 num = S_PO_S(ptr);
3649 if (LT(num,cons_null) == TRUE)
3650 addinvers_apply(num);
3651 else
3652 break;
3653 ptr = S_L_N(ptr);
3654 }
3655 }
3656 ptr = list_copy;
3657 while (ptr != NULL)
3658 {
3659 num = S_PO_S(ptr);
3660 if (not einsp(num) && not nullp(num))
3661 {
3662 integer_factor(num,prime_list);
3663 list_ptr = prime_list;
3664 while (list_ptr != NULL)
3665 {
3666 new = CALLOCOBJECT();
3667 copy(S_PO_S(list_ptr),new);/* new is the next prime */
3668 next_ptr = S_L_N(ptr);
3669 while (next_ptr != NULL)
3670 {
3671 next_num = S_PO_S(next_ptr);
3672 if (NEQ(next_num,cons_eins) == TRUE)
3673 {
3674 nb_quores(next_num,new,quo,rem);
3675 if (nullp(rem)) /* AK 120891 */
3676 copy(quo,next_num);
3677 }
3678 next_ptr = S_L_N(next_ptr);
3679 }
3680 insert_list(new,data_ptr,NULL,NULL);
3681 list_ptr = S_L_N(list_ptr);
3682 }
3683 freeself(prime_list);
3684 }
3685 ptr = S_L_N(ptr);
3686 }
3687 fsd_ende:
3688 freeall(list_copy);
3689 freeall(prime_list);
3690 freeall(rem);
3691 freeall(quo);
3692 return(OK);
3693 #else
3694 error("find_sqrad_data: SQ_RADICAL not available");
3695 return(ERROR);
3696 #endif
3697 }
3698
3699 /* a: the sqrad */
3700
adjust_sqrad_data(a)3701 static INT adjust_sqrad_data(a) OP a;
3702 /* 15.04.91: TPMcD. */
3703 /* AK 200891 V1.3 */
3704 {
3705 INT inserted = 1L;
3706 INT erg = OK;
3707
3708 OP new=NULL, quo, rem, ptr, data_ptr, a_copy, prime_list, num_ptr;
3709 if (S_O_K(a) != SQ_RADICAL)
3710 return(ERROR);
3711 if (S_N_D(a) == NULL || EMPTYP(S_N_D(a)))
3712 return(find_sqrad_data(a));
3713 if (empty_listp(S_N_D(a)))
3714 return(OK);
3715 prime_list = CALLOCOBJECT();
3716 init(LIST,prime_list);
3717 a_copy = CALLOCOBJECT();
3718 copy(a,a_copy);
3719 ptr = S_N_S(a_copy);
3720 num_ptr = S_PO_S(ptr);
3721 if (LT(num_ptr,cons_null) == TRUE) /*negative radicals */
3722 {
3723 new = CALLOCOBJECT();
3724 M_I_I(-1L,new);
3725 insert_list(new,prime_list,NULL,NULL);
3726 while (ptr != NULL)/*multiply negative radicals by -1 */
3727 {
3728 num_ptr = S_PO_S(ptr);
3729 if (LT(num_ptr,cons_null) == TRUE)
3730 addinvers_apply(num_ptr);
3731 else
3732 break;
3733 ptr = S_L_N(ptr);
3734 }
3735 }
3736 data_ptr = S_N_D(a);
3737 quo = CALLOCOBJECT();
3738 rem = CALLOCOBJECT();
3739 while (data_ptr != NULL)
3740 {
3741 if (negeinsp(S_L_S(data_ptr))) /* negatives have been taken care of*/
3742 {
3743 data_ptr = S_L_N(data_ptr);
3744 continue;
3745 }
3746 if (inserted)
3747 new = CALLOCOBJECT();
3748 copy(S_L_S(data_ptr),new); /* new is the next prime*/
3749 inserted = 0L;
3750 ptr = S_N_S(a_copy);
3751 while (ptr != NULL)
3752 {
3753 num_ptr = S_PO_S(ptr);
3754 if (einsp(num_ptr) || nullp(num_ptr))
3755 {
3756 ptr = S_L_N(ptr);
3757 continue;
3758 }
3759 nb_quores(num_ptr,new,quo,rem);
3760 if (nullp(rem))
3761 {
3762 copy(quo,num_ptr);
3763 if (not inserted)
3764 {
3765 insert(new,prime_list,NULL,NULL);
3766 inserted = 1L;
3767 }
3768 }
3769 ptr = S_L_N(ptr);
3770 }
3771 data_ptr = S_L_N(data_ptr);
3772 }
3773 if (not inserted)
3774 FREESELF(new);
3775 else
3776 new = CALLOCOBJECT();
3777 make_monopoly_sqrad(S_N_S(a_copy),new); /* reconstitute the sqrad */
3778 if (convert_sqrad_scalar(new) == ERROR)
3779 {
3780 FREESELF(S_N_D(a));
3781 FREEALL(prime_list);
3782 find_sqrad_data(a);
3783 }
3784 else
3785 {
3786 FREEALL(S_N_D(a));
3787 S_N_D(a) = prime_list;
3788 }
3789 FREEALL(quo);
3790 FREEALL(rem);
3791 FREEALL(new);
3792 FREEALL(a_copy);
3793 return(OK);
3794 ENDR("adjust_sqrad_data");
3795 }
3796
3797 /* a: the sqrad; b: the radical; c: the conjugate */
3798
conj_sqrad(a,b,c)3799 INT conj_sqrad(a,b,c) OP a,b,c;
3800 /* AK 200891 V1.3 */ /* 23.10.91: TPMcD */
3801 {
3802 OP la, lb, rem, minus, ptr, new;
3803 #ifdef SQRADTRUE
3804
3805 if (not EMPTYP(c)) /* AK 060993 */
3806 freeself(c);
3807
3808 la = CALLOCOBJECT();
3809 lb = CALLOCOBJECT();
3810 rem = CALLOCOBJECT();
3811 minus = CALLOCOBJECT();
3812 M_I_I(-1L,minus);
3813 init(MONOPOLY,la);
3814 init(MONOPOLY,lb);
3815 ptr = S_N_S(a);
3816 if (EQ(b,minus) == TRUE)
3817 while (ptr != NULL)
3818 {
3819 new = CALLOCOBJECT();
3820 copy(S_L_S(ptr),new);
3821 if (LT(S_MO_S(new),cons_null) == TRUE)
3822 insert_list(new,lb,NULL,NULL);
3823 else
3824 insert_list(new,la,NULL,NULL);
3825 ptr = S_L_N(ptr);
3826 }
3827 else
3828 while (ptr != NULL)
3829 {
3830 new = CALLOCOBJECT();
3831 copy(S_L_S(ptr),new);
3832 mod(S_MO_S(new),b,rem);
3833 if (nullp(rem)) /* AK 120891 */
3834 insert_list(new,lb,NULL,NULL);
3835 else
3836 insert_list(new,la,NULL,NULL);
3837 ptr = S_L_N(ptr);
3838 }
3839 if (empty_listp(lb))
3840 insert_zero_into_monopoly(lb);
3841 mult_apply_scalar_monopoly(minus,lb);
3842 insert(lb,la,NULL,NULL);
3843 if (c == a)
3844 freeall(S_N_S(a));
3845 else
3846 {
3847 init(SQ_RADICAL,c);
3848 copy(S_N_D(a),S_N_D(c));
3849 }
3850 remove_zero_terms(la);
3851 if (S_N_S(c) != NULL)
3852 {
3853 freeall(S_N_S(c)); /* AK 060993 */
3854 }
3855 S_N_S(c) = la;
3856 freeall(rem);
3857 freeall(minus);
3858 return(OK);
3859 #else /* SQRADTRUE */
3860 error("conj_sqrad: SQ_RADICAL not available");
3861 return(ERROR);
3862 #endif /* SQRADTRUE */
3863 }
3864
3865 #ifdef SQRADTRUE
squareroot_integer(a,b)3866 INT squareroot_integer(a,b) OP a,b;
3867 /* AK 040291 V1.2 */ /* AK 200891 V1.3 */
3868 {
3869 INT erg = OK;
3870 OP c;
3871 CTTO(INTEGER,LONGINT,"squareroot_integer(1)",a);
3872 CTO(EMPTY,"squareroot_integer(2)",b);
3873
3874 if (NULLP_INTEGER(a)) {
3875 M_I_I(0,b);
3876 goto ende;
3877 }
3878
3879 c = CALLOCOBJECT();
3880 erg += b_skn_mp(CALLOCOBJECT(),CALLOCOBJECT(),NULL,c);
3881 COPY_INTEGER(a,S_PO_S(c));
3882 M_I_I(1,S_PO_K(c));
3883 erg += make_monopoly_sqrad(c,b);
3884 FREEALL(c);
3885 ende:
3886 ENDR("squareroot_integer");
3887 }
3888
3889
3890
squareroot_longint(a,b)3891 INT squareroot_longint(a,b) OP a,b;
3892 /* AK 040291 V1.2 */ /* AK 200891 V1.3 */
3893 {
3894 return squareroot_integer(a,b);
3895 }
3896 #endif /* SQRADTRUE */
3897
3898 #ifdef SQRADTRUE
squareroot_bruch(a,b)3899 INT squareroot_bruch(a,b) OP a,b;
3900 /* AK 040291 V1.2 */ /* AK 200891 V1.3 */ /* 04.10.91: TPMcD */
3901 {
3902 INT erg=OK;
3903 OP c,d;
3904 CTO(BRUCH,"squareroot_bruch(1)",a);
3905 CTO(EMPTY,"squareroot_bruch(2)",b);
3906
3907
3908 c = CALLOCOBJECT();
3909 d = CALLOCOBJECT();
3910 MULT(S_B_O(a),S_B_U(a),c);
3911 erg += squareroot(c,b);
3912 erg += invers(S_B_U(a),d);
3913 MULT_APPLY(d,b);
3914 FREEALL(c);
3915 FREEALL(d);
3916 CTO(SQ_RADICAL,"squareroot_bruch(e2)",b);
3917 ENDR("squareroot_bruch");
3918 }
3919
3920
3921
convert_sqrad_scalar(a)3922 INT convert_sqrad_scalar(a) OP a;
3923 /* 5.04.91: TPMcD. */ /* AK 200891 V1.3 */
3924 {
3925 INT ret = ERROR;
3926 OP tmp;
3927 if (S_O_K(a) != SQ_RADICAL || S_L_N(S_N_S(a)) != NULL)
3928 return(ret);
3929 tmp = S_PO_S(S_N_S(a));
3930 if (not einsp(tmp) && not nullp(tmp))
3931 return(ret);
3932 ret = OK;
3933 if (nullp(tmp))
3934 {
3935 m_i_i(0L,a);
3936 return(ret);
3937 }
3938 tmp = CALLOCOBJECT();
3939 copy(S_PO_K(S_N_S(a)),tmp);
3940 copy(tmp,a);
3941 freeall(tmp);
3942 return OK;
3943 }
3944 #endif /* SQRADTRUE */
3945
3946 /* Convert the square root of an integer a to a cyclotomic number */
3947
convert_radical_cyclo(a,b)3948 INT convert_radical_cyclo(a,b) OP a,b;
3949 /* AK 200891 V1.3 */ /* 29.10.91: TPMcD */
3950 {
3951 INT myeven, posi, last = 1L,ff=0;
3952 OP new, ptr, mono_ptr;
3953 INT erg = OK;
3954 # ifdef NUMBERTRUE
3955 OP k,m,mpos,x,y,z,w,atmp,four;
3956
3957
3958 CTTO(INTEGER,LONGINT,"convert_radical_cyclo(1)",a);
3959 if (EINSP(a)) ff=1;
3960
3961
3962 k = CALLOCOBJECT();
3963 m = CALLOCOBJECT();
3964 mpos = CALLOCOBJECT();
3965 x = CALLOCOBJECT();
3966 y = CALLOCOBJECT();
3967 z = CALLOCOBJECT();
3968 w = CALLOCOBJECT();
3969 four = CALLOCOBJECT();
3970
3971
3972
3973 if (not negp(a) && nb_square_root(a,k) == OK)
3974 {
3975 make_scalar_cyclo(k,b);
3976 goto exit_label;
3977 }
3978 if (a == b)
3979 {
3980 atmp = CALLOCOBJECT();
3981 copy(a,atmp);
3982 }
3983 else
3984 atmp = a;
3985
3986
3987 FREESELF(b); INIT_CYCLO(b);
3988
3989 mono_ptr = CALLOCOBJECT();
3990 init(MONOPOLY,mono_ptr);
3991 M_I_I(4L,four);
3992 integer_factor_1(atmp,cons_zwei,cons_zwei,m,y,NULL);
3993 /* a = 4k * 2l * m
3994 * l=0 ,m=1(4): return 2k * r(m)
3995 * l=0 ,m=3(4): return 2(k-1) * r(4*m)
3996 * l=1 return 2(k-1) * r(8*m)
3997 */
3998 ptr = y;
3999 if (empty_listp(ptr)) /* a > 0 and a is odd */
4000 {
4001 myeven = 0L;
4002 posi = 1L;
4003 }
4004 else if (EQ(S_PO_S(ptr),cons_zwei)) /* a > 0 and a is even */
4005 {
4006 myeven = 1L;
4007 posi = 1L;
4008 }
4009 else
4010 {
4011 ptr = S_L_N(ptr);
4012 if (ptr == NULL) /* a < 0 and a is odd */
4013 {
4014 myeven = 0L;
4015 posi = 0L;
4016 }
4017 else
4018 {
4019 myeven = 1L;
4020 posi = 0L;
4021 }
4022 }
4023 if (!posi)
4024 addinvers_apply(m);
4025 if (myeven)
4026 {
4027 nb_quores(S_PO_K(ptr),cons_zwei,k,z);
4028 if (NULLP(z)) /* AK 120891 */ /* a = 4k * m */
4029 last = 0L;
4030 hoch(cons_zwei,k,w); /* w = 2k */
4031 }
4032 else
4033 {
4034 copy(cons_eins,w);
4035 last = 0L;
4036 }
4037 if (!last)
4038 {
4039 mod(m,four,z);
4040 if (!EINSP(z))
4041 {
4042 div_apply(w,cons_zwei);
4043 MULT_APPLY(four,m);
4044 }
4045 }
4046 else
4047 {
4048 div_apply(w,cons_zwei); /* w := w/2 */
4049 MULT_APPLY(four,m);
4050 MULT_APPLY(cons_zwei,m);
4051 }
4052 copy(m,mpos);
4053 if (NEGP(mpos))
4054 addinvers_apply(mpos);
4055 make_coprimes(mpos,y);
4056 ptr = y;
4057 while (ptr != NULL)
4058 {
4059 if (kronecker(m,S_L_S(ptr),z) == OK)
4060 {
4061 new = CALLOCOBJECT();
4062 m_sk_mo(S_L_S(ptr),z,new);
4063 insert(new,mono_ptr,add_koeff,NULL);
4064 }
4065 ptr = S_L_N(ptr);
4066 }
4067 remove_zero_terms(mono_ptr);
4068 make_index_monopoly_cyclo(mpos,mono_ptr,b,0L);
4069 MULT_APPLY(w,b);
4070
4071 FREEALL(mono_ptr);
4072 if (a == b)
4073 FREEALL(atmp);
4074 exit_label:
4075 FREEALL(k);
4076 FREEALL(m);
4077 FREEALL(mpos);
4078 FREEALL(x);
4079 FREEALL(y);
4080 FREEALL(z);
4081 FREEALL(w);
4082 FREEALL(four);
4083 erg += standardise_cyclo_0(b,basis_type);
4084 { /* so to get the positive squareroot */
4085 double realpart=0.0;
4086 double fac = 6.2830 / ((double) S_N_DCII(b) );
4087 FORALL(w,S_N_S(b), {
4088 realpart += (double)(S_MO_KI(w)) * cos( fac* (double)S_I_I(S_MO_S(w)) );
4089 });
4090
4091
4092 if (realpart < 0.0 )
4093 if (ff==0) ADDINVERS_APPLY(b);
4094 }
4095 # endif
4096 CTO(ANYTYPE,"convert_radical_cyclo(e2)",b);
4097 ENDR("convert_radical_cyclo");
4098 }
4099
4100
4101 # ifdef NUMBERTRUE
convert_sqrad_cyclo(a,b)4102 INT convert_sqrad_cyclo(a,b) OP a,b;
4103 /* 29.10.91: TPMcD */
4104 {
4105 OP c, ptr;
4106 INT erg = OK;
4107 CTO(SQ_RADICAL,"convert_sqrad_cyclo(1)",a);
4108 CE2(a,b,convert_sqrad_cyclo);
4109
4110 M_I_I(0,b);
4111
4112 c = CALLOCOBJECT();
4113 ptr = S_N_S(a);
4114 while (ptr != NULL)
4115 {
4116 convert_radical_cyclo(S_PO_S(ptr),c);
4117 MULT_APPLY(S_PO_K(ptr),c);
4118 ADD_APPLY(c,b);
4119 ptr = S_L_N(ptr);
4120 }
4121
4122 FREEALL(c);
4123 ENDR("convert_sqrad_cyclo");
4124 }
4125 #endif /* NUMBERTRUE */
4126
4127 /****************** fields_3.c **********************/
4128 /* 26.07.91: TPMcD. */
4129 /*************************************************************/
4130
4131 /* CYCLOTOMIC */
4132
4133 /* a : the index, b : the monopoly, c : the result */
4134
trans_index_monopoly_cyclo(a,b,c)4135 INT trans_index_monopoly_cyclo(a,b,c) OP a,b,c;
4136 /* AK 300791 for compatibility */
4137 /* AK 200891 V1.3 */
4138 {
4139 return make_index_monopoly_cyclo(a,b,c,POWER_REDUCE);
4140 }
4141
make_index_monopoly_cyclo(a,b,c,tid)4142 static INT make_index_monopoly_cyclo(a,b,c,tid) OP a,b,c; INT tid;
4143 /* 30.05.90: TPMcD. */ /* 3.04.91: TPMcD. */
4144 /* AK 200891 V1.3 */
4145 /* 23.10.91: TPMcD */
4146 {
4147 OP c_tmp;
4148 INT flag = 0L;
4149 INT erg = OK;
4150 CYCLO_DATA *c_ptr = NULL;
4151
4152 if (S_O_K(b) != MONOPOLY)
4153 error("make_index_monopoly_cyclo: 2nd parameter is wrong type\n");
4154 if ((c_ptr = add_cyclo_data(a)) == (CYCLO_DATA *) NULL)
4155 error("make_index_monopoly_cyclo: unable to create cyclotomic data\n");
4156 if (c == b)
4157 {
4158 flag = 1L;
4159 c_tmp = CALLOCOBJECT();
4160 }
4161 else
4162 {
4163 FREESELF(c); INIT_CYCLO(c);
4164 c_tmp = S_N_S(c);
4165 }
4166 init(MONOPOLY, c_tmp);
4167 if (empty_listp(c_tmp))
4168 insert_zero_into_monopoly(c_tmp);
4169 copy(b, c_tmp);
4170 if (flag)
4171 {
4172 init(CYCLOTOMIC,c);
4173 S_N_S(c) = c_tmp;
4174 }
4175 S_N_DC(c) = c_ptr;
4176 if (tid != NO_REDUCE)
4177 standardise_cyclo_0(c,tid);
4178 ENDR("make_index_monopoly_cyclo");
4179 }
4180
standardise_cyclo(a)4181 INT standardise_cyclo(a) OP a;
4182 /* 25.10.91: TPMcD */
4183 {
4184 return(standardise_cyclo_0(a,basis_type));
4185 }
4186
standardise_cyclo_0(a,tid)4187 static INT standardise_cyclo_0(a,tid) OP a;
4188 INT tid;
4189 /* 09.09.90: TPMcD. */ /* 4.04.91: TPMcD. */
4190 /* AK 200891 V1.3 */
4191 /* 25.10.91: TPMcD */
4192 {
4193 INT erg=OK, ret = ERROR;
4194 CYCLO_DATA *c_ptr;
4195 OP ptr, new, mono, e, poly_eins, poly_zwei;
4196 /* OP half=NULL; */
4197
4198
4199 if (S_O_K(a) != CYCLOTOMIC || tid == NO_REDUCE)
4200 return(OK);
4201
4202 /*
4203 if (EVEN(S_N_DC(a)->index)) {
4204 half = CALLOCOBJECT();
4205 ganzdiv(S_N_DC(a)->index,cons_zwei,half);
4206 }
4207 */
4208 ptr = S_N_S(a);
4209 c_ptr = S_N_DC(a);
4210 mono = CALLOCOBJECT();
4211 init(MONOPOLY,mono);
4212 e = CALLOCOBJECT();
4213 if ( not empty_listp(ptr))
4214 while (ptr != NULL)
4215 {
4216 erg = mod(S_PO_S(ptr),c_ptr->index,e);
4217 new = CALLOCOBJECT();
4218 m_sk_mo(e,S_PO_K(ptr),new);
4219 /*
4220 if (EVEN(c_ptr->index))
4221 if (GE(e,half)) {
4222 sub_apply(half,S_MO_S(new));
4223 ADDINVERS_APPLY(S_MO_K(new));
4224 }
4225 */
4226 insert(new,mono,add_koeff,NULL);
4227 ptr = S_L_N(ptr);
4228 }
4229 FREEALL(e);
4230
4231 /* if (EVEN(S_N_DC(a)->index)) FREEALL(half); */
4232
4233 if (empty_listp(mono))
4234 insert_zero_into_monopoly(mono);
4235 switch((int)tid)
4236 {
4237 case (int)POWER_REDUCE:
4238 poly_zwei = mono;
4239 break;
4240 case (int)STD_BASIS:
4241 poly_eins = CALLOCOBJECT();
4242 poly_zwei = CALLOCOBJECT();
4243 quores_monopoly(mono,c_ptr->poly,poly_eins,poly_zwei);
4244 FREEALL(mono);
4245 FREEALL(poly_eins);
4246 break;
4247 default:
4248 return error("standardise_cyclo_0: unknown cyclotomic basis");
4249 break;
4250 }
4251 FREEALL(S_N_S(a));
4252 S_N_S(a) = poly_zwei;
4253 ret = OK;
4254 return(ret);
4255 ENDR("standardise_cyclo_0");
4256 }
4257
4258 #ifdef CYCLOTRUE
make_scalar_cyclo(a,b)4259 INT make_scalar_cyclo(a,b) OP a,b;
4260 /* 5.04.91: TPMcD. */
4261 /* AK 200891 V1.3 */
4262 /* 23.10.91: TPMcD */
4263 {
4264 OP c = CALLOCOBJECT();
4265 OP d = CALLOCOBJECT();
4266 M_I_I(1L,c);
4267 b_skn_mp(CALLOCOBJECT(),CALLOCOBJECT(),NULL,d);
4268 copy(a,S_PO_K(d));
4269 M_I_I(0L,S_PO_S(d));
4270 make_index_monopoly_cyclo(c,d,b,NO_REDUCE);
4271 freeall(c);
4272 freeall(d);
4273 return(OK);
4274 }
4275
4276
make_index_power_cyclo(a,c,d)4277 INT make_index_power_cyclo(a,c,d) OP a,c,d;
4278 {
4279 return make_index_coeff_power_cyclo(a,cons_eins,c,d);
4280 }
4281
make_index_coeff_power_cyclo(a,b,c,d)4282 INT make_index_coeff_power_cyclo(a,b,c,d) OP a,b,c,d;
4283 /* 30.05.90: TPMcD. */ /* 17.07.91: TPMcD. */
4284 /* AK 200891 V1.3 */ /* 23.10.91: TPMcD */
4285 {
4286 INT erg = OK;
4287 erg += init(CYCLOTOMIC,d);
4288 erg += b_skn_mp(CALLOCOBJECT(),CALLOCOBJECT(),NULL,S_N_S(d));
4289 erg += mod(c,a,S_PO_S(S_N_S(d)));
4290 erg += copy(b,S_PO_K(S_N_S(d)));
4291
4292 if (S_N_DC(d) != NULL)
4293 error("internal error:MIC1");
4294
4295 S_N_DC(d) = add_cyclo_data(a);
4296
4297 if (space_saving)
4298 convert_cyclo_scalar(d);
4299
4300 ENDR("make_index_coeff_power_cyclo");
4301 }
4302
4303
scan_cyclo(a)4304 INT scan_cyclo(a) OP a;
4305 /* AK 240191 V1.2 */ /* a becomes cyclotomic */
4306 /* AK 200891 V1.3 */ /* 23.10.91: TPMcD */
4307 {
4308 OP b = CALLOCOBJECT();
4309 OP c = CALLOCOBJECT();
4310 INT erg = OK;
4311 erg += printeingabe("degree of cyclotomic field");
4312 erg += scan(INTEGER,b);
4313 erg += printeingabe("self of cyclotomic field");
4314 erg += scan(MONOPOLY,c);
4315 erg += make_index_monopoly_cyclo(b,c,a,basis_type);
4316 erg += freeall(b);
4317 erg += freeall(c);
4318 return erg;
4319 }
4320 #endif /* CYCLOTRUE */
4321
4322 /* a: the scalar, b: the cyclo, c: the result */
4323
add_scalar_cyclo(a,b,c)4324 INT add_scalar_cyclo(a,b,c) OP a,b,c;
4325 /* 30.05.90: TPMcD. */ /* AK 080891 V1.3 */ /* 23.10.91: TPMcD */
4326 {
4327 OP ptr;
4328 INT erg = OK;
4329 #ifdef CYCLOTRUE
4330 if (c == a)
4331 error("First and third arguments equal\n");
4332 if (c != b)
4333 copy(b,c);
4334 ptr = CALLOCOBJECT();
4335 erg += init(MONOPOLY,ptr);
4336 C_L_S(ptr,CALLOCOBJECT());
4337 erg += m_sk_mo(cons_null,a,S_L_S(ptr));
4338 erg += add_apply(ptr,S_N_S(c));
4339 erg += freeall(ptr);
4340 if (space_saving)
4341 convert_cyclo_scalar(c);
4342 #endif /* CYCLOTRUE */
4343 return erg;
4344 }
4345
4346 /* a: the scalar, b: the cyclo, c: the result */
4347
4348 #ifdef CYCLOTRUE
mult_apply_scalar_cyclo(a,b)4349 INT mult_apply_scalar_cyclo(a,b) OP a,b;
4350 /* AK 060498 V2.0 */
4351 {
4352 OP c;
4353 INT erg = OK;
4354 CTO(CYCLOTOMIC,"mult_apply_scalar_cyclo",b);
4355 c = CALLOCOBJECT();
4356 SWAP(c,b);
4357 erg += mult_scalar_cyclo(a,c,b);
4358 erg += freeall(c);
4359 ENDR("mult_apply_scalar_cyclo");
4360 }
4361
mult_scalar_cyclo(a,b,c)4362 INT mult_scalar_cyclo(a,b,c) OP a, b, c;
4363 /* 06.09.90: TPMcD. */ /* AK 200891 V1.3 */ /* 23.10.91: TPMcD */
4364 {
4365 INT erg = OK;
4366 CTO(CYCLOTOMIC,"mult_scalar_cyclo(2)",b);
4367 CTO(EMPTY,"mult_scalar_cyclo(3)",c);
4368
4369 if (NULLP(a)) {
4370 M_I_I(0,c);
4371 }
4372 else {
4373 erg += init(CYCLOTOMIC,c);
4374 FREESELF(S_N_S(c));
4375 erg += mult_scalar_monopoly(a,S_N_S(b),S_N_S(c));
4376 S_N_DC(c) = S_N_DC(b);
4377 if (space_saving)
4378 convert_cyclo_scalar(c);
4379 }
4380 ENDR("mult_scalar_cyclo");
4381 }
4382
4383 #endif
4384 /* a,b: the cyclos, c: the result */
4385
add_cyclo_cyclo(a,b,c)4386 INT add_cyclo_cyclo(a,b,c) OP a,b,c;
4387 /* AK 200891 V1.3 */
4388 /* 23.10.91: TPMcD */
4389 {
4390 return( add_cyclo_cyclo_0(a,b,c,basis_type) );
4391 }
4392
add_cyclo_cyclo_0(a,b,c,tid)4393 static INT add_cyclo_cyclo_0(a,b,c,tid) OP a,b,c; INT tid;
4394 /* 06.09.90: TPMcD. */ /* 5.04.91: TPMcD. */
4395 /* AK 200891 V1.3 */ /* 23.10.91: TPMcD */
4396 {
4397 OP temp_eins, temp_zwei, temp_drei, temp_vier;
4398 INT erg = OK;
4399
4400 if (S_O_K(a) != CYCLOTOMIC || S_O_K(b) != CYCLOTOMIC)
4401 return( error ("add_cyclo_cyclo_0: argument not CYCLOTOMIC") );
4402 temp_eins = CALLOCOBJECT();
4403 temp_zwei = CALLOCOBJECT();
4404 temp_drei = CALLOCOBJECT();
4405 temp_vier = CALLOCOBJECT();
4406 copy(S_N_S(a),temp_eins);
4407 copy(S_N_S(b),temp_zwei);
4408 ggt(S_N_DCI(a),S_N_DCI(b),temp_drei);
4409 ganzdiv(S_N_DCI(a),temp_drei,temp_vier);
4410 raise_power_monopoly(temp_vier,temp_zwei);
4411 ganzdiv(S_N_DCI(b),temp_drei,temp_vier);
4412 raise_power_monopoly(temp_vier,temp_eins);
4413 MULT_APPLY(S_N_DCI(a),temp_vier);
4414 init(CYCLOTOMIC, c);
4415 FREESELF(S_N_S(c));
4416 add_monopoly_monopoly(temp_eins,temp_zwei,S_N_S(c));
4417 S_N_DC(c) = add_cyclo_data(temp_vier);
4418 if (tid != NO_REDUCE)
4419 standardise_cyclo_0(c,tid);
4420 if (space_saving)
4421 convert_cyclo_scalar(c);
4422 FREEALL(temp_eins);
4423 FREEALL(temp_zwei);
4424 FREEALL(temp_drei);
4425 FREEALL(temp_vier);
4426 ENDR("nb.c:add_cyclo_cyclo_0");
4427 }
4428
mult_cyclo_cyclo(a,b,c)4429 INT mult_cyclo_cyclo(a,b,c) OP a,b,c;
4430 /* AK 200891 V1.3 */
4431 {
4432 INT erg = OK;
4433 CTO(CYCLOTOMIC,"mult_cyclo_cyclo(1)",a);
4434 CTO(CYCLOTOMIC,"mult_cyclo_cyclo(2)",b);
4435 CTO(EMPTY,"mult_cyclo_cyclo(3)",c);
4436
4437 erg += mult_cyclo_cyclo_0(a,b,c,basis_type); /* AK return inserted */
4438 CTO(ANYTYPE,"mult_cyclo_cyclo(i3)",c);
4439 erg += standardise_cyclo_0(c,basis_type);
4440 CTO(ANYTYPE,"mult_cyclo_cyclo(e3)",c);
4441 ENDR("mult_cyclo_cyclo");
4442 }
4443
mult_cyclo_cyclo_0(a,b,c,tid)4444 static INT mult_cyclo_cyclo_0(a,b,c,tid) OP a,b,c; INT tid;
4445 /* 06.09.90: TPMcD. */ /* 5.04.91: TPMcD. */
4446 /* AK 200891 V1.3 */
4447 /* 23.10.91: TPMcD */
4448 {
4449 OP temp_eins, temp_zwei, temp_drei, temp_vier;
4450 INT erg = OK;
4451 CTO(CYCLOTOMIC,"mult_cyclo_cyclo_0(1)",a);
4452 CTO(CYCLOTOMIC,"mult_cyclo_cyclo_0(2)",b);
4453
4454 if ( (NULLP(a) || NULLP(b)) && space_saving )
4455 {
4456 m_i_i(0L,c);
4457 return(OK);
4458 }
4459 temp_eins = CALLOCOBJECT();
4460 temp_zwei = CALLOCOBJECT();
4461 temp_drei = CALLOCOBJECT();
4462 temp_vier = CALLOCOBJECT();
4463 COPY(S_N_S(a),temp_eins);
4464 COPY(S_N_S(b),temp_zwei);
4465 ggt(S_N_DCI(a),S_N_DCI(b),temp_drei);
4466 ganzdiv(S_N_DCI(a),temp_drei,temp_vier);
4467 raise_power_monopoly(temp_vier,temp_zwei);
4468 ganzdiv(S_N_DCI(b),temp_drei,temp_vier);
4469 raise_power_monopoly(temp_vier,temp_eins);
4470 MULT_APPLY(S_N_DCI(a),temp_vier);
4471 init(CYCLOTOMIC, c);
4472 FREESELF(S_N_S(c));
4473 mult_monopoly_monopoly(temp_eins,temp_zwei,S_N_S(c));
4474 if ((S_N_DC(c) = add_cyclo_data(temp_vier)) == (CYCLO_DATA *) NULL)
4475 error("Unable to find cyclotomic data\n");
4476 if (tid != NO_REDUCE)
4477 standardise_cyclo_0(c,tid);
4478 if (space_saving)
4479 convert_cyclo_scalar(c);
4480 FREEALL(temp_eins);
4481 FREEALL(temp_zwei);
4482 FREEALL(temp_drei);
4483 FREEALL(temp_vier);
4484 ENDR("mult_cyclo_cyclo_0");
4485 }
4486
add_cyclo(a,b,c)4487 INT add_cyclo(a,b,c) OP a,b,c;
4488 /* AK 070891 V1.3 */
4489 {
4490 INT erg = OK;
4491 #ifdef CYCLOTRUE
4492 switch(S_O_K(b))
4493 {
4494 case INTEGER:
4495 case LONGINT:
4496 case SQ_RADICAL:
4497 case BRUCH:
4498 erg += add_scalar_cyclo(b,a,c);
4499 break;
4500 case CYCLOTOMIC:
4501 erg += add_cyclo_cyclo(a,b,c);
4502 break;
4503 case POLYNOM:
4504 erg += add_scalar_polynom(a,b,c);
4505 break;
4506 default:
4507 printobjectkind(b);
4508 erg += error("add_cyclo: wrong second type\n");
4509 break;
4510 }
4511 convert_cyclo_scalar(c);
4512 #endif
4513 return erg;
4514 }
4515
mult_cyclo(a,b,c)4516 INT mult_cyclo(a,b,c) OP a,b,c;
4517 /* 24.07.91: TPMcD. */
4518 /* AK 200891 V1.3 */
4519 {
4520 INT erg = OK;
4521 CTO(CYCLOTOMIC,"mult_cyclo(1)",a);
4522 CTO(EMPTY,"mult_cyclo(3)",c);
4523
4524 if (NULLP(a)){
4525 M_I_I(0,c);
4526 goto ende;
4527 }
4528 if (NULLP(b)){
4529 M_I_I(0,c);
4530 goto ende;
4531 }
4532 switch(S_O_K(b))
4533 {
4534 case INTEGER:
4535 case SQ_RADICAL: /* AK 220891 */
4536 case LONGINT:
4537 case BRUCH:
4538 erg += mult_scalar_cyclo(b,a,c);
4539 break;
4540 #ifdef MATRIXTRUE
4541 case MATRIX:
4542 erg += mult_scalar_matrix(a,b,c);
4543 break;
4544 #endif /* MATRIXTRUE */
4545 #ifdef POLYTRUE
4546 case SCHUR:
4547 case POW_SYM:
4548 case HOM_SYM:
4549 case ELM_SYM:
4550 case MONOMIAL:
4551 case POLYNOM:
4552 erg += mult_scalar_polynom(a,b,c);
4553 break;
4554 #endif /* POLYTRUE */
4555 #ifdef SCHUBERTTRUE
4556 case SCHUBERT:
4557 erg += mult_scalar_schubert(a,b,c);
4558 break;
4559 #endif /* SCHUBERTTRUE */
4560 case VECTOR:
4561 erg += mult_scalar_vector(a,b,c);
4562 break;
4563 case CYCLOTOMIC:
4564 erg += mult_cyclo_cyclo(a,b,c);
4565 break;
4566 default:
4567 WTO("mult_cyclo(2)",b);
4568 break;
4569 }
4570 convert_cyclo_scalar(c);
4571 ende:
4572 ENDR("mult_cyclo");
4573 }
4574
addinvers_cyclo(a,b)4575 INT addinvers_cyclo(a,b) OP a,b;
4576 /* AK 200891 V1.3 */
4577 {
4578 INT erg = OK;
4579 CTO(CYCLOTOMIC,"addinvers_cyclo(1)",a);
4580 CTO(EMPTY,"addinvers_cyclo(2)",b);
4581 erg += mult_scalar_cyclo(cons_negeins,a,b);
4582 CTO(CYCLOTOMIC,"addinvers_cyclo(e2)",b);
4583 ENDR("addinvers_cyclo");
4584 }
4585
4586 /* a: the cyclo, b: the auto, c: the conjugate */
4587
conj_cyclo(a,b,c)4588 INT conj_cyclo(a,b,c) OP a,b,c;
4589 /* AK 200891 V1.3 */
4590 /* 23.10.91: TPMcD */
4591 {
4592 # ifdef CYCLOTRUE
4593 if (S_O_K(a) != CYCLOTOMIC)
4594 return(ERROR);
4595 if (c != a)
4596 copy(a,c);
4597 raise_power_monopoly(b,S_N_S(c));
4598 standardise_cyclo_0(c,basis_type);
4599 # endif
4600 return(OK);
4601 }
4602
4603 /* a: the cyclo, b: the inverse */
4604
invers_cyclo(a,b)4605 INT invers_cyclo(a,b) OP a,b;
4606 /* AK 200891 V1.3 */
4607 {
4608 INT ret = ERROR;
4609 # ifdef CYCLOTRUE
4610 OP c = CALLOCOBJECT();
4611 ret = invers_cyclo_norm(a,b,c);
4612 freeall(c);
4613 # endif
4614 return(ret);
4615 }
4616
4617 /* a: the cyclo, b: the inverse, c: the norm */
4618
invers_cyclo_norm(a,b,c)4619 static INT invers_cyclo_norm(a,b,c) OP a,b,c;
4620 /* AK 200891 V1.3 */
4621 {
4622 INT flag = 0L, saving = space_saving;
4623 # ifdef CYCLOTRUE
4624 OP ptr, tmp, b_tmp;
4625 if (S_O_K(a) != CYCLOTOMIC)
4626 return(ERROR);
4627 if (nullp_cyclo(a))
4628 return(error("invers_cyclo_norm: cannot invert 0\n"));
4629 if (c == a || c == b)
4630 return(error("invers_cyclo_norm: illegal 3rd parameter\n"));
4631 if (b == a)
4632 {
4633 b_tmp = CALLOCOBJECT();
4634 flag = 1L;
4635 }
4636 else
4637 {
4638 if (not EMPTYP(b)) /* AK */
4639 freeself(b);
4640 b_tmp = b;
4641 }
4642 space_saving = FALSE;
4643 tmp = CALLOCOBJECT();
4644 /*
4645 M_I_I(1L,tmp);
4646 */
4647 make_scalar_cyclo(cons_eins,b_tmp);
4648 ptr = S_N_DC(a)->autos;
4649 ptr = S_L_N(ptr); /* Skip the trivial automorphism */
4650 while (ptr != NULL)
4651 {
4652 conj_cyclo(a,S_L_S(ptr),tmp);
4653 mult_cyclo_cyclo_0(tmp,b_tmp,b_tmp,POWER_REDUCE);
4654 ptr = S_L_N(ptr);
4655 }
4656 mult_cyclo_cyclo_0(a,b_tmp,tmp,basis_type);
4657 if (convert_cyclo_scalar(tmp) == ERROR)
4658 {
4659 freeall(tmp);
4660 if (flag)
4661 freeall(b_tmp);
4662 return(error("invers_cyclo_norm: norm is not scalar\n"));
4663 }
4664 copy(tmp,c);
4665 if (negp(tmp))
4666 {
4667 mult_scalar_sqrad(cons_negeins,b_tmp,b_tmp);
4668 addinvers_apply(tmp);
4669 }
4670 invers(tmp,tmp);
4671 mult_apply_scalar_cyclo(tmp,b_tmp);
4672 if (flag)
4673 {
4674 copy(b_tmp,b);
4675 freeall(b_tmp);
4676 }
4677 freeall(tmp);
4678 space_saving = saving;
4679 # endif
4680 return(OK);
4681 }
4682
add_apply_cyclo(a,b)4683 INT add_apply_cyclo(a,b) OP a,b;
4684 /* AK 200891 V1.3 */
4685 {
4686 INT erg = OK;
4687 OP c;
4688 CTO(CYCLOTOMIC,"add_apply_cyclo(1)",a);
4689 c = CALLOCOBJECT();
4690 SWAP(c,b);
4691 erg += add_cyclo(a,c,b);
4692 FREEALL(c);
4693 ENDR("add_apply_cyclo");
4694 }
4695
mult_apply_cyclo(a,b)4696 INT mult_apply_cyclo(a,b) OP a,b;
4697 /* AK 200891 V1.3 */
4698 {
4699 INT ret;
4700 #ifdef CYCLOTRUE
4701 OP c;
4702 c = CALLOCOBJECT();
4703 ret = mult_cyclo(a,b,c);
4704 copy(c,b);
4705 freeall(c);
4706 #endif
4707 return(ret);
4708 }
4709
addinvers_apply_cyclo(a)4710 INT addinvers_apply_cyclo(a) OP a;
4711 /* AK 200891 V1.3 */
4712 {
4713 OP b;
4714 INT erg = OK;
4715 CTO(CYCLOTOMIC,"addinvers_apply_cyclo(1)",a);
4716 b = CALLOCOBJECT();
4717 SWAP(b,a);
4718 erg += addinvers_cyclo(b,a);
4719 FREEALL(b);
4720 ENDR("addinvers_apply_cyclo");
4721 }
4722
nullp_cyclo(a)4723 INT nullp_cyclo(a) OP a;
4724 /* AK 200891 V1.3 */
4725 {
4726 #ifdef CYCLOTRUE
4727 if (S_O_K(a) != CYCLOTOMIC)
4728 return(ERROR);
4729 if (EMPTYP(S_N_S(a)))
4730 {
4731 error("nullp_cyclo: cyclo with empty self\n");
4732 return(TRUE);
4733 }
4734 return(nullp_monopoly(S_N_S(a)));
4735 #else
4736 return(ERROR);
4737 #endif
4738 }
4739 #ifdef CYCLOTRUE
comp_cyclo(a,b)4740 INT comp_cyclo(a,b) OP a,b;
4741 /* AK 200891 V1.3 */
4742 {
4743 return(comp_list(S_N_S(a),S_N_S(b)));
4744 }
4745 #endif /* CYCLOTRUE */
4746
4747
4748 # ifdef CYCLOTRUE
convert_cyclo_scalar(a)4749 INT convert_cyclo_scalar(a) OP a;
4750 /* 5.04.91: TPMcD. */
4751 /* AK 200891 V1.3 */
4752 {
4753 INT ret = ERROR;
4754 OP tmp;
4755
4756 if (S_O_K(a) != CYCLOTOMIC || S_L_N(S_N_S(a)) != NULL)
4757 goto exit_label;
4758 tmp = S_PO_S(S_N_S(a));
4759 if (not nullp(tmp))
4760 goto exit_label;
4761 tmp = CALLOCOBJECT();
4762 copy(S_PO_K(S_N_S(a)),tmp);
4763 copy(tmp,a);
4764 freeall(tmp);
4765 ret = OK;
4766 exit_label:
4767 return(ret);
4768 }
4769 #endif /* CYCLOTRUE */
4770
4771 #ifdef CYCLOTRUE
fprint_cyclo(f,a)4772 static INT fprint_cyclo(f,a) FILE *f;
4773 OP a;
4774 /* 25.09.91: TPMcD. */
4775 {
4776 INT myfirst = 1L, flag;
4777 OP ptr;
4778
4779 standardise_cyclo_0(a,basis_type);
4780 ptr = S_N_S(a);
4781 zeilenposition += 2L;
4782 if (nullp_cyclo(a))
4783 {
4784 fprintf(f," 0");
4785 return(OK);
4786 }
4787 while (ptr != NULL)
4788 {
4789 flag = 0L;
4790 if (zeilenposition > 60L)
4791 {
4792 zeilenposition = 0L;
4793 fprintf(f,"\n");
4794 }
4795 if (!negp(S_PO_K(ptr)) && !myfirst)
4796 fprintf(f," +");
4797 if (negeinsp(S_PO_K(ptr)))
4798 {
4799 flag = 1L;
4800 fprintf(f," -");
4801 }
4802 else if (!einsp(S_PO_K(ptr)))
4803 {
4804 fprintf(f," ");
4805 fprint(f,S_PO_K(ptr));
4806 }
4807 else
4808 {
4809 fprintf(f," ");
4810 flag = 1L;
4811 }
4812
4813 if (not nullp(S_PO_S(ptr)))
4814 {
4815 fprintf(f," E(");
4816 fprint(f,S_N_DCI(a));
4817 fprintf(f,")");
4818 if (!einsp(S_PO_S(ptr)))
4819 {
4820 fprintf(f,"^");
4821 fprint(f,S_PO_S(ptr));
4822 }
4823 zeilenposition += 5L;
4824 }
4825 else if (flag)
4826 fprintf(f," 1");
4827 ptr = S_L_N(ptr);
4828 myfirst = 0L;
4829 }
4830 return(OK);
4831 }
4832 #endif /* CYCLOTRUE */
4833
4834 #ifdef CYCLOTRUE
tex_cyclo(a)4835 INT tex_cyclo(a) OP a;
4836 /* 4.04.91: TPMcD. */ /* AK 200891 V1.3 */ /* 23.10.91: TPMcD */
4837 {
4838 INT myfirst = 1L;
4839 OP ptr = S_N_S(a);
4840
4841 if (nullp_cyclo(a))
4842 {
4843 fprintf(texout," 0\n");
4844 return(OK);
4845 }
4846 fprintf(texout,"\n");
4847 while (ptr != NULL)
4848 {
4849 if (!negp(S_PO_K(ptr)) && !myfirst)
4850 fprintf(texout," + {");
4851 else
4852 fprintf(texout,"{");
4853 tex(S_PO_K(ptr));
4854
4855 if (not nullp(S_PO_S(ptr)))
4856 {
4857 fprintf(texout,"} \\omega_{");
4858 tex(S_N_DCI(a));
4859 fprintf(texout,"} {");
4860 tex(S_PO_S(ptr));
4861 }
4862 fprintf(texout,"}\n");
4863 ptr = S_L_N(ptr);
4864 myfirst = 0L;
4865 }
4866 fprintf(texout,"\n");
4867 return(OK);
4868 }
4869 #endif /* CYCLOTRUE */
4870
4871 /* ROUTINES RELATING TO THE MAINTENANCE OF CYCLOTOMIC DATA */
4872
4873 # ifdef CYCLOTRUE
4874
4875 /*Reads the table of cyclos from the file CYCLOS.DAT. The first entry */
4876 /*should be no_cyclos, then the list of cyclo_data. Returns OK if the */
4877 /*table is set; otherwise, returns ERROR. */
4878
setup_cyclotomic_table(filename)4879 static INT setup_cyclotomic_table(filename) char *filename;
4880 /* 30.08.90: TPMcD */ /* AK 200891 V1.3 */
4881 {
4882 INT i=0;
4883 FILE *f;
4884 CYCLO_DATA *ptr;
4885 char name[50], *char_ptr;
4886
4887 if (cyclo_table_set || filename == NULL)
4888 return(OK);
4889 if ((f = fopen(filename,"r")) == NULL)
4890 {
4891 printf("\nFile containing cyclo data: ");
4892 char_ptr = name;
4893 while( (*char_ptr = fgetc(stdin)) != '\n')
4894 {
4895 if (myisspace(*char_ptr)) continue;
4896 char_ptr++;
4897 i++;
4898 if (i > (INT)48) break;
4899 }
4900 *char_ptr = /* NULL; AK 290494 */ '\0';
4901 if (strlen(name) == 0)
4902 return(ERROR);
4903 if ((f = fopen(name,"r")) == NULL)
4904 {
4905 printf("Unable to open %s\n",name);
4906 return(ERROR);
4907 }
4908 }
4909 if ( fscanf(f, " %" PRIINT ,&zzno_cyclos) == 0 || zzno_cyclos < 1L ||
4910 (zzcyclo_table
4911 = (CYCLO_DATA *) SYM_calloc((int)zzno_cyclos,sizeof(CYCLO_DATA))
4912 ) == NULL
4913 )
4914 {
4915 zzno_cyclos = 0L;
4916 printf("\nCyclo data table creation error");
4917 return(ERROR);
4918 }
4919 ptr = zzcyclo_table - 1;
4920 for (i=0L;i<zzno_cyclos;i++)
4921 {
4922 ptr++;
4923 ptr->index = CALLOCOBJECT();
4924 objectread(f,ptr->index);
4925 ptr->deg = CALLOCOBJECT();
4926 objectread(f,ptr->deg);
4927 ptr->poly = CALLOCOBJECT();
4928 objectread(f,ptr->poly);
4929 ptr->autos = CALLOCOBJECT();
4930 objectread(f,ptr->autos);
4931 }
4932 cyclo_table_set = 1L;
4933 fclose(f);
4934 return(OK);
4935 }
4936
add_cyclo_data(index)4937 static CYCLO_DATA *add_cyclo_data(index) OP index;
4938 /* AK 200891 V1.3 */
4939 {
4940 CYCLO_DATA *ptr = NULL;
4941 OP ptr_eins, ptr_zwei=NULL;
4942
4943 if ((ptr = cyclo_ptr(index)) != NULL)
4944 return(ptr);
4945 ptr = (CYCLO_DATA *) SYM_calloc(1,sizeof(CYCLO_DATA));
4946 if (ptr == NULL)
4947 return(NULL);
4948 ptr->index = CALLOCOBJECT();
4949 COPY(index,ptr->index);
4950 ptr->poly = CALLOCOBJECT();
4951 if (make_cyclotomic_monopoly(index,ptr->poly) == ERROR)
4952 {
4953 SYM_free(ptr);
4954 return(NULL);
4955 }
4956 ptr_eins = ptr->poly;
4957 while(ptr_eins != NULL)
4958 {
4959 ptr_zwei = ptr_eins;
4960 ptr_eins = S_L_N(ptr_eins);
4961 }
4962 ptr->deg = CALLOCOBJECT();
4963 COPY(S_PO_S(ptr_zwei),ptr->deg);
4964 ptr->autos = CALLOCOBJECT();
4965 make_coprimes(ptr->index,ptr->autos);
4966 ptr_eins = CALLOCOBJECT();
4967 init(LIST,ptr_eins);
4968 /* Some compilers require this cast, others dislike it */
4969 /* (CYCLO_DATA *) S_L_S(ptr_eins) = ptr; */
4970 C_L_S(ptr_eins,ptr);
4971 /* S_L_N(ptr_eins) = NULL; */
4972 C_L_N(ptr_eins,NULL);
4973 if (cyclo_list_set)
4974 S_L_N(zzcyclo_tail) = ptr_eins;
4975 else
4976 {
4977 cyclo_list_set = 1L;
4978 zzcyclo_list = ptr_eins;
4979 }
4980 zzcyclo_tail = ptr_eins;
4981 return(ptr);
4982 }
4983
cyclo_ptr(index)4984 static CYCLO_DATA *cyclo_ptr(index) OP index;
4985 /* AK 200891 V1.3 */
4986 {
4987 CYCLO_DATA *ptr = NULL;
4988 OP list_ptr;
4989 INT i;
4990 if (cyclo_table_set)
4991 {
4992 ptr = zzcyclo_table - 1;
4993 for (i=0L;i<zzno_cyclos;i++)
4994 {
4995 ptr++;
4996 if (EQ(index,ptr->index) == TRUE)
4997 return(ptr);
4998 }
4999 }
5000 if (cyclo_list_set)
5001 {
5002 list_ptr = zzcyclo_list;
5003 while (list_ptr != NULL)
5004 {
5005 ptr = (CYCLO_DATA *) S_L_S(list_ptr);
5006 if (ptr == NULL)
5007 error("cyclo_ptr: null pointer\n");
5008 if (EQ(index,ptr->index) == TRUE)
5009 return(ptr);
5010 list_ptr = S_L_N(list_ptr);
5011 }
5012 }
5013 return(NULL);
5014 }
5015
free_cyclo_list()5016 static INT free_cyclo_list()
5017 /* 29.10.91: TPMcD */
5018 {
5019 OP list_ptr;
5020 OBJECTSELF list_self;
5021 CYCLO_DATA *cp;
5022 list_ptr = zzcyclo_list;
5023 while (list_ptr != NULL)
5024 {
5025 list_self = S_O_S(list_ptr);
5026 cp = (CYCLO_DATA *)S_L_S(list_ptr);
5027 freeall(cp->index);
5028 freeall(cp->deg);
5029 freeall(cp->poly);
5030 freeall(cp->autos);
5031 SYM_free(cp);
5032 C_L_S(list_ptr,NULL); /* Wg speicherverwaltung */
5033 list_ptr = S_L_N(list_ptr);
5034 }
5035 return(OK);
5036 }
5037
print_cyclo_data(ptr)5038 INT print_cyclo_data(ptr) CYCLO_DATA *ptr;
5039 /* AK 200891 V1.3 */
5040 {
5041 printf("Index ");
5042 print(ptr->index);
5043 printf("\tDegree ");
5044 println(ptr->deg);
5045 printf("Polynomial ");
5046 println(ptr->poly);
5047 printf("Automorphism exponents ");
5048 println(ptr->autos);
5049 return OK;
5050 }
5051
free_cyclo_table()5052 static INT free_cyclo_table()
5053 /* AK 310893 */
5054 {
5055 CYCLO_DATA *ptr;
5056 INT i;
5057
5058 if (!cyclo_table_set)
5059 return(ERROR);
5060 ptr = zzcyclo_table;
5061 for (i=0L;i<zzno_cyclos;i++)
5062 {
5063 freeall(ptr->index);
5064 freeall(ptr->deg);
5065 freeall(ptr->poly);
5066 freeall(ptr->autos);
5067 ptr++;
5068 }
5069 return(OK);
5070 }
5071
print_cyclo_table()5072 INT print_cyclo_table()
5073 /* AK 200891 V1.3 */
5074 {
5075 CYCLO_DATA *ptr;
5076 INT i;
5077
5078 if (!cyclo_table_set)
5079 return(ERROR);
5080 printf( "Number of cyclo data on table: %" PRIINT "\n" ,zzno_cyclos);
5081 ptr = zzcyclo_table;
5082 for (i=0L;i<zzno_cyclos;i++)
5083 {
5084 printf( "Table item %" PRIINT ": " ,i);
5085 print_cyclo_data(ptr);
5086 ptr++;
5087 }
5088 return(OK);
5089 }
5090
print_cyclo_list()5091 INT print_cyclo_list()
5092 /* AK 200891 V1.3 */
5093 {
5094 OP list_ptr;
5095 INT i = 0L;
5096
5097 if (!cyclo_list_set)
5098 return(ERROR);
5099 printf("Cyclo data in list:\n");
5100 list_ptr = zzcyclo_list;
5101 while (list_ptr != NULL)
5102 {
5103 printf( "List item %" PRIINT ": " ,i++);
5104 print_cyclo_data((CYCLO_DATA *) S_L_S(list_ptr));
5105 list_ptr = S_L_N(list_ptr);
5106 }
5107 return(OK);
5108 }
5109
save_cyclo_list(filename)5110 INT save_cyclo_list(filename) char *filename;
5111 /* 4.04.91: TPMcD. */ /* AK 200891 V1.3 */
5112 {
5113 CYCLO_DATA *ptr;
5114 OP list_ptr;
5115 INT i = 0L, new = 0L;
5116 FILE *f;
5117 char name[50], *char_ptr;
5118
5119 if (filename == NULL || (f = fopen(filename,"r+")) == NULL)
5120 {
5121 fflush(stdin);
5122 printf("\nFile to receive cyclo data: ");
5123 char_ptr = name;
5124 while( (*char_ptr = fgetc(stdin)) != '\n')
5125 {
5126 if (myisspace(*char_ptr)) continue;
5127 char_ptr++;
5128 i++;
5129 if (i > (INT)48) break;
5130 }
5131 *char_ptr = /* NULL; AK 290494 */ '\0';
5132 if (strlen(name) == 0)
5133 return(ERROR);
5134 if ((f = fopen(name,"r+")) == NULL)
5135 {
5136 if((f = fopen(name,"w")) == NULL)
5137 {
5138 printf("Unable to open %s\n",name);
5139 return(ERROR);
5140 }
5141 else
5142 new = 1L;
5143 }
5144 }
5145 else
5146 strcpy(name,filename);
5147 if (new)
5148 {
5149 fprintf(f," \n\n");
5150 printf("Initialising %s\n",name);
5151 i = 0L;
5152 }
5153 else
5154 {
5155 fseek(f,0L,0);
5156 fscanf(f, "%" SCNINT ,&i);
5157 fseek(f,0L,2);
5158 printf("Cyclo data being appended to file %s.\n",name);
5159 }
5160 list_ptr = zzcyclo_list;
5161 while (list_ptr != NULL)
5162 {
5163 ptr = (CYCLO_DATA *) S_L_S(list_ptr);
5164 fprintf(f,"\n");
5165 objectwrite(f,ptr->index);
5166 objectwrite(f,ptr->deg);
5167 objectwrite(f,ptr->poly);
5168 objectwrite(f,ptr->autos);
5169 list_ptr = S_L_N(list_ptr);
5170 i++;
5171 }
5172 fseek(f,0L,0);
5173 fprintf(f, "%8" PRIINT ,i);
5174 fclose(f);
5175 return(OK);
5176 }
5177
5178 #endif
5179
5180 #ifdef NUMBERTRUE
test_number()5181 INT test_number()
5182 {
5183 OP a = CALLOCOBJECT();
5184 OP b = CALLOCOBJECT();
5185 printeingabe("test_number: squareroot(2L,a)");
5186 squareroot(cons_zwei,a);
5187 println(a);
5188 printeingabe("test_number: squareroot(11L,a)^-1");
5189 m_i_i(19L,b);
5190 squareroot(b,a);
5191 invers(a,b);
5192 println(b);
5193 printeingabe("test_number: euler_phi(311L,a)");
5194 m_i_i(311L,b);
5195 euler_phi(b,a);
5196 println(b);
5197 freeall(a);
5198 freeall(b);
5199 return OK;
5200 }
5201 #endif /* NUMBERTRUE */
5202
t_MONOPOLY_POLYNOM(a,b)5203 INT t_MONOPOLY_POLYNOM(a,b) OP a,b;
5204 /* AK 171194 */
5205 /* AK 170206 V3.0 */
5206 {
5207 INT erg = OK;
5208 CTO(MONOPOLY,"t_MONOPOLY_POLYNOM(1)",a);
5209 CE2(a,b,t_MONOPOLY_POLYNOM);
5210
5211 {
5212 OP c;
5213 init(POLYNOM,b);
5214 if (S_L_S(a) != NULL) {
5215 while (a != NULL)
5216 {
5217 c = CALLOCOBJECT();
5218 m_iindex_iexponent_monom(0,S_I_I(S_PO_S(a)),c);
5219 copy(S_PO_K(a),S_PO_K(c));
5220 insert(c,b,NULL,NULL);
5221 a = S_L_N(a);
5222 }
5223 }
5224 }
5225
5226 ENDR("t_MONOPOLY_POLYNOM");
5227 }
5228
5229
invers_monopoly(lau,res)5230 INT invers_monopoly(lau, res) OP lau,res;
5231 {
5232 INT erg = OK;
5233 CTO(MONOPOLY,"invers_monopoly(1)",lau);
5234 CTO(EMPTY,"invers_monopoly(2)",res);
5235
5236 erg += b_ou_b(CALLOCOBJECT(),CALLOCOBJECT(),res);
5237 M_I_I((INT)1,S_B_O(res));
5238 erg += copy(lau,S_B_U(res));
5239 C_B_I(res,GEKUERZT);
5240
5241 ENDR("invers_monopoly");
5242 }
5243
5244
degree_monopoly(mp,dg)5245 INT degree_monopoly(mp,dg) OP mp,dg;
5246 /*CC 010496*/
5247 /* -1 if null */
5248 {
5249 /* Puts in dg the degree of the MONOPOLY object mp*/
5250 OP z,za=NULL;
5251 INT erg = OK;
5252 CTO(MONOPOLY,"degree_monopoly(1)",mp);
5253 FREESELF(dg);
5254
5255 if(NULLP(mp))
5256 M_I_I(-1L,dg);
5257 else {
5258 z=mp;
5259 while(z !=NULL) { za=z; z=S_L_N(z); }
5260 COPY(S_PO_S(za),dg);
5261 }
5262 ENDR("degree_monopoly");
5263 }
5264
5265 /*
5266 Puts in ld the leading coefficient of the MONOPOLY object mp.
5267 */
5268
ldcf_monopoly(mp,ld)5269 INT ldcf_monopoly(mp,ld) OP mp,ld;
5270 {
5271 INT erg=OK;
5272 OP z,za=NULL;
5273 FREESELF(ld);
5274 if (NULLP(mp)) error("ldcf_monopoly: null monopoly");
5275 else {
5276 z=mp;
5277 while(z !=NULL) { za=z; z=S_L_N(z); }
5278 COPY(S_PO_K(za),ld);
5279 }
5280 ENDR("ldcf_monopoly");
5281 }
5282
5283
5284
has_one_variable(a)5285 INT has_one_variable(a) OP a;
5286 /* AK 310106 */
5287 {
5288 /*
5289 Returns TRUE, if a is an MONOPOLY object, or
5290 is of type POLYNOM with 0 or 1 variable.
5291 */
5292 OP nb;
5293 INT erg =OK;
5294 if(S_O_K(a)==MONOPOLY)
5295 return TRUE;
5296 if(S_O_K(a)==POLYNOM)
5297 {
5298 nb=CALLOCOBJECT();
5299 numberofvariables(a,nb);
5300 if(S_I_I(nb)<=1L)
5301 {
5302 FREEALL(nb);return TRUE;
5303 }
5304 FREEALL(nb);
5305 }
5306 return FALSE;
5307 ENDR("has_one_variable");
5308 }
5309
5310