/***** MACSYMA *****/
/* Macsyma commands for "Mathematics for Multimedia" */

/**** Chapter 1 ****/

/* Division algorithm -- integers and polynomials */

quotient(a,b);	/*  a/b, not always nonnegative */
remainder(a,b);	/*  a%b, not always nonnegative */
divide(a,b);    /* [quotient(a,b), remainder(a,b)]  */
mod(a,b);   /*  the remainder a%b in the range {0,1,...,b-1} */

gcd(a,b);	/* greatest common divisor */
lcm(a,b)	/* least common multiple */


/* Primes */ 

next_prime(x);		/* Miller-Rabin test for least prime p>x */
factor(x);		/* the prime factors of x */
totient(x);		/* Euler's totient function */

/* Base conversion -- to base 16 */

ibase:10;  obase:16;		/* input base 10, output base 16 */
/* Note: append a decimal point to an integer to force base-10 interpretation */ 
ibase:2; obase:10.;		/* input base 2, output base 10 */


/* IEEE 754 (64-bit and 32-bit floating point format) */
float(2^-1023);			/* smallest positive normal float */
float(2^-1074);			/* smallest positive subnormal float */

/*** Euler's totient and the RSA encryption/decryption algorithm */

/* Find two large primes p,q. */
/* For this demo, 4 digit primes will be used. */ 
factor(2329);                   /* some big number, check if it is prime */
primep(2329);                   /* ...or use a primality test */
p:next_prime(2329);             /* 2333 */
q:next_prime(1776);             /* 1777 */
N:p*q;                          /* 4145741 */
phi:(p-1)*(q-1);                /* totient for N=p*q, with p,q prime */ 
totient(N);                     /* 4141632, slow to compute for big N */
e:next_prime(isqrt(N));         /* encryption exponent should be big */
next_prime(floor(sqrt(N)));     /* ...isqrt() is integer square root */ 
gcd(phi,e);                     /* check that e and phi are coprime */
d:inv_mod(e,phi);               /* modular inverse: e*d=1 (mod phi) */
mod(d*e,phi);                   /* ...so this must be 1   */
/* Alternatively, use the extended Euclidean algorithm: */
/* gcdex(a,b) returns [x,y,gcd] satisfying x*a+y*b=gcd */
/* Thus gcdex(a,b)[1] is the inverse of a (mod b) if gcd(a,b)==1 */ 
gcdex(e,phi);                   
d:gcdex(e,phi)[1];              /* first element is the inverse of e */
mod(d*e,phi);                   /* ...so this must be 1   */
/* Thus e,N is public for encryption; d,phi is private for decryption */
/* Encrypt: */ 
cleartext:123456;               /* message to be sent; no larger than N */
cyphertext: power_mod (cleartext, e, N); /* cleartext^e (mod N) */
/* Decrypt: */
decrypted: power_mod(cyphertext, d, N);  /* cyphertext^d (mod N) */

is(decrypted=cleartext);        /* Are these the same? True */;
power_mod(cyphertext+1, d, N);  /* Transmission error? text is garbled */
power_mod(cyphertext+2, d, N);  /* Transmission error? text is garbled */
power_mod(cyphertext+3, d, N);  /* Transmission error? text is garbled */

/* Finite field computations manual: */
/* http://cs.swan.ac.uk/~csoliver/ok-sat-library/internet_html/doc/doc/Maxima/5.29.0/gf_manual.pdf */