GAP 4.8.9 installation with standard packages -- copy to your CoCalc project to get it
#include "typedef.h"1#include "tools.h"2/**************************************************************************\3@---------------------------------------------------------------------------4@---------------------------------------------------------------------------5@ FILE: tools.c6@ containes some general small tools.7@---------------------------------------------------------------------------8@---------------------------------------------------------------------------9@10@11@ Exports the globale variables Zero and One.12@ rational Zero= { 0,1 };13@ rational One = { 1,1 };14@15@ int GGT( int a, int b );16@ calculates the greatest common divisor of a and b.17@ The result is always > 018@19@ int KGV( int a, int b );20@ calculates the least commom multiple of a and b.21@ The result is always > 022@23@ void Normal ( rational *a );24@ divides the numerator and the denominator of a rational number25@ by their greatest common divisor26@27@ void Normal2 ( int *z, *n );28@ divides z and n by their greatest common divisor.29@ Caution: z and n are changed in Normal230@31@ void rat_add ( int *az, int *an, int bz, int bn);32@ cz az bz33@ calculates -- := -- + --34@ cn an bn35@ The result is stored in az and an.36@37@ int *factorize_new( int zahl,int *erg);38@39@ int *factorize ( int zahl );40@ The integer 'zahl' is factorized. If for example41@ zahl = p1^a1 *p2^a2,42@ erg[p1] is set to a1 and erg[p2] is set to a2,43@ all other entries of erg are 0.44@ Caution: the result is a pointer to an integer of length 100.45@ So factorize works only for intergers that46@ have prime factors smaller then 100.47@48@ void gcd_darstell( int a1, int a2, int *v1, int *v2, int *gcd);49@ calculates a presentation of the greatest commom divisor of a1 and a2:50@ gcd = v1*a1 + v2*a251@52@ int p_inv(int a, int p)53@ calculates number ai such that a * ai is kongruent 1 modulo p (if exists)54@55@-------------------------------------------------------------------------56\**************************************************************************/5758/**********************************************************************\59|60| tools.c -- allgemeine kleinere tools:61|62| int GGT( int a, int b );63| int KGV( int a, int b );64| void Normal ( rational *a );65| void Normal2 ( int *z, *n );66| void rat_add ( int *az, int *an, int bz, int bn);67| int *factorize_new( int zahl,int *erg);68| int *factorize ( int zahl );69| void gcd_darstell( int a1, int a2, int *v1, int *v2, int *gcd);70| int p_inv(int a, int p);71| int signum(int a);72|73| exportiert die globale Variable Zero. Eine Moeglichkeit, alles zum74| Absturz zu bringen, ist also z.B. die Anweisung Zero.n = 0 ... :(75|76\**********************************************************************/7778rational Zero= { 0,1 };79rational One = { 1,1 };8081/*{{{}}}*/82/*{{{ GGT*/83int GGT (int _a, int _b)84{85register int a= _a;86register int b= _b;87register int c;8889if ( b == 0 ) {90return abs(a);91} else if ( b == 1 || a == 1 ) {92return 1;93} else {94while ( (c = a%b) != 0 ) {95a = b;96b = c;97}98}99return abs(b);100}101102/*}}} */103/*{{{ KGV*/104/*105|106| int KGV( int a, int b);107|108| berechnet das kgv von a und b. Das Ergebnis ist immer >= 0!109|110*/111int KGV( a, b )112int a, b;113{114int kgv, ggt;115116a = abs(a);117b = abs(b);118ggt= GGT( a, b );119if ( ggt != 0 )120{121if ( a > b )122kgv = (a / ggt) * b; /* die Reihenfolge ist wesentlich! */123/* (a*b)/ggt waere grosser Bloedsinn! */124else125kgv = (b / ggt) * a; /* die Reihenfolge ist wesentlich! */126}127else128kgv = 0;129return kgv;130}131132/*}}} */133/*{{{ Normal*/134135136void Normal (a)137rational *a;138{139register int g;140register int n= a->n;141register int z= a->z;142143if ( n == 0 ) {144printf ("Normal: Error: divide by zero\n");145exit (3);146}147if ( z == 0 ) {148a->n= 1;149} else {150g = GGT ( z, n);151if ( g != 1 ) {152z /= g;153n /= g;154}155if ( n < 0 ) {156a->z = -z;157a->n = -n;158} else {159a->z = z;160a->n = n;161}162}163}164165/*}}} */166/*{{{ Normal2*/167void Normal2 ( _z, _n )168int *_z, *_n;169{170register int z = *_z;171register int n = *_n;172register int g;173174if ( n == 0 ) {175printf ("Normal2: Error: divide by zero\n");176exit (3);177}178if ( z == 0 ) {179*_n = 1;180} else {181g = GGT ( z, n);182if ( g != 1 ) {183z /= g;184n /= g;185}186if ( n < 0) {187*_z = -z;188*_n = -n;189} else {190*_z = z;191*_n = n;192}193}194}195196/*}}} */197/*{{{ rat_add*/198/*199|200| rat_add( &a.z, &a.n, b.z, b.n );201|202| rational a, rational b: addition of a and b, result is stored in a.203|204*/205void rat_add( az, an, bz, bn )206int *az, *an, bz, bn;207{208register int temp_ggt;209210/*211* Normal tests also for divide by zero212*/213Normal2( az, an );214Normal2( &bz, &bn );215temp_ggt = GGT( *an, bn );216*az = (*az) * bn/temp_ggt + bz * (*an)/temp_ggt;217if ( *an > bn ) {218*an = (*an/temp_ggt) * bn;219} else {220*an *= (bn/temp_ggt);221}222Normal2( az, an );223}224225/*}}} */226/*{{{ factorize_new*/227int *factorize_new( zahl, erg)228int zahl;229int *erg;230{231int i;232233if ( erg != NULL )234{235for(i=0; i<100; i++)236erg[i] = 0;237for(i=2; i<100; i++)238{239while((zahl%i) == 0)240{241zahl = zahl/i;242erg[i]++;243}244}245if(zahl != 1)246erg[0] = 1;247}248return( erg );249}250251/*}}} */252/*{{{ factorize*/253int *factorize(zahl)254int zahl;255{256int i;257int *erg;258259erg = (int *) malloc(100 *sizeof(int));260factorize_new( zahl, erg );261return(erg);262}263264/*}}} */265/*{{{ gcd_darstell*/266/*-----------------------------------------------------------------------*\267| gibt darstellung des ggt: gcd = v1*a1 + v2*a2 |268\*-----------------------------------------------------------------------*/269void gcd_darstell(a1, a2, v1, v2, gcd)270int a1, a2;271int *v1, *v2, *gcd;272{273int bn,bn1,bn2,rn,rn1,rn2,q, an, an1, an2;274rn2=a2; rn1=a1; bn2=0; bn1=1;275an2 = 1; an1 = 0;276277if(a1 == 0 && a2 == 0)278{279*gcd = 0; *v1 = 0; *v2 = 0;280}281if(a1 == 0 && a2 != 0)282{283*gcd = a2; *v1 = 0; *v2 = 1;284}285if(a1 != 0 && a2 == 0)286{287*gcd = a1; *v1 = 1; *v2 = 0;288}289if(a1 != 0 && a2 != 0)290{291do292{293rn=rn2%rn1;294q=(rn2-rn)/rn1;295bn=bn2-q*bn1;296an=an2-q*an1;297if(rn != 0)298{299bn2=bn1; bn1=bn; rn2=rn1; rn1=rn; an2=an1; an1 = an;300}301}302while(rn!=0);303if(rn1 < 0)304{305*gcd = -rn1; *v1 = -bn1; *v2 = -an1;306}307if(rn1 > 0)308{309*gcd = rn1; *v1 = bn1; *v2 = an1;310}311}312}313/*}}} */314315/*-----------------------------------------------------------------------*\316| calculates number i such that a * i is kongruent 1 modulo p (if exists)317\*-----------------------------------------------------------------------*/318int p_inv(a, p)319int a,p;320{321int an, an1, an2, rn, rn1, rn2, q;322rn1 = a %p;323if(rn1 == 0)324{325printf("%d is not invertible modulo %d\n", a, p);326exit(3);327}328if(rn1 == 1 || rn1 == -1)329return(rn1);330rn2 = p; an2 = 0; an1 = 1;331do332{333rn=rn2%rn1;334q=(rn2-rn)/rn1;335an=an2-q*an1;336if(rn != 0)337{338rn2=rn1; rn1=rn; an2=an1; an1 = an;339}340}341while(rn!=0);342if(rn1 == -1)343return( -an1);344if(rn1 == 1)345return( an1);346else347{348printf("%d is not invertible modulo %d\n", a, p);349exit(3);350}351}352353354int signum(a)355int a;356{357if(a>0)358return(1);359if(a<0)360return(-1);361return(0);362}363364365