%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%% OCTAVE %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % STEP 0: Installing GNU Octave % % Linux: apt install octave liboctave-dev make gfortran % Mac: brew install octave % Windows installer %%%%% Chapter 1 of "Mathematics for Multimedia" %%% Division algorithm -- integers idivide(a,b,"floor"); % returns floor(a/b) rem(a,b); % returns the remainder a%b, not always nonnegative mod(a,b); % returns the remainder a%b in the range {0,1,...,b-1} gcd(a,b) % greatest common divisor lcm(a,b) % least common multiple %%% Primes factor(x) isprime(x) primes(n) % list the primes =< n list_primes(n) % list the first n primes %%% Euler's totient function phi = totient(n) [p,m]=factor(n); % unique primes p with multiplicities m phi = prod(p.^m-p.^(m-1)); % product property of the totient function end %%% Paul Godfrey's version of Euler's totient() function [t] = pgtotient(n) %TOTIENT calculates the totient function (also % called the Euler Phi function) of any % positive integer n. % %Usage: f = pgtotient(n) % % The totient function computes the number of % integers from the set (1...n-1) that are % relatively prime to n. It can be used to % describe the multiplicative structure of % all Galois fields GF(q). % % n can be any size % % Tested under version 5.2.0 % %Ref: "Error Control Systems" by S.B Wicker % %see also Primes, Factor % % Paul Godfrey % pgodfrey@conexant.com % 10-23-2001 % % Comments by MVW [r c]=size(n); % allow matrix inputs n=reshape(n,1,r*c); % reshape into column vector t=zeros(1,r*c); % allocate space for output f=zeros(1,10); % temporary space for prime factors for k=1:r*c; nk=n(k); % kth element in matrix n f=unique(factor(nk)); % unique prime factors t(k)=nk*prod(1-1./f); % real-valued totient end t=reshape(t,r,c); % restore the matrix p=find(n==1); % trivial case: totient of 1 is 1 t(p)=1; t=round(t); % totient is an integer, so round the floats return %a demo of this program is n=(1:50).'; t=mtotient(n); [n t] %%% Next prime function p = nextprime(n) p = n; % initialize while(~isprime(p)) % increment while p is not prime p=p+1; end end %%% Modular inverse (or quasi-inverse) qi*x=d=gcd(x,m) (mod m) function [qi, d] = invmod(x,m) [d, qi, y] = gcd(x,m); % extended Euclidean algorithm end %%% Modular power function pow = powermod(a,k,m) % a^k (mod m) if ( k<0 ) a = invmod(a,m); % ...then use a^(-1) k = -k; % ...and -k, which is positive end if ( a<0 ) a = mod(m-a,m); % ...then use equivalent a>=0 end pow = 1; while( k>0 ) if( rem(k,2)>0 ) pow = mod(a*pow, m); end a = mod(a*a, m); % square of a (mod m) k = floor(k/2); % put next bit of k in the ones place end end %%% RSA encryption and decryption factor(2329) % some big number, check if it is prime isprime(2329) % ...or use a primality test p = nextprime(2329) % 2333 q = nextprime(1776) % 1777 N = p*q % 4145741 modulus, the product of two large primes phi = (p-1)*(q-1) % totient for N=p*q, with p,q, prime phi = totient(N) % 4141632 e = nextprime(floor(sqrt(N))) % encryption exponent should be big gcd(phi,e) % e and phi must be coprime d = invmod(e,phi) % modular inverse: e*d=gcd(e,phi) (mod phi) mod(d*e,phi); % ...so this must be 1 % Alternatively, use the extended Euclidean algorithm: % [gcd,x,y] = gcd(a,b) returns x,y, satisfying x*a+y*b=gcd % This x is the inverse of a (mod b) if gcd(a,b)==1 [g,d,y]=gcd(e,phi); 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 = powermod (cleartext, e, N) % cleartext^e (mod N) % Decrypt: decrypted = powermod(cyphertext, d, N) % cyphertext^d (mod N) %%% Base conversion % positive integers in binary, decimal, and hexadecimal hex2dec("DEADBEEF") % hexadecimal to decimal, string input, case-insensitive bin2dec("1010101") % binary to decimal, string input dec2bin(2017) % decimal to binary, integer input dec2hex(104479) % decimal to hexadecimal, integer input base2dec ("11120", 3) % conversion from base-3 to decimal base2dec ("aeaai", "aei" ) % conversion using symbols for digits 0,1,2 dec2base (123, 3) % conversion from decimal to base-3 dec2base (123, "aei") % conversion using symbols for digits 0,1,2 % Using base conversion to encode words alphabet = "abcdefghijklmnopqrstuvwxyz" base2dec("hello",alphabet) % Real numbers: 1+10 hex digits of pi, with radix point [dec2hex(3) "." dec2hex(round(16^10 * 0.1415926535897932384626433832795))] % IEEE 754 (64-bit and 32-bit floating point format) num2hex(3.1415926535897932384626433832795) % hexadecimal string, 16 characters, 8 bytes num2hex(single(2.718281828459)) % hexadecimal string, 8 characters, 4 bytes dec2bin(hex2dec(num2hex(1/3))) % binary string, IEEE 64-bit float format dec2bin(hex2dec(num2hex(single(1/3)))) % binary string, IEEE 32-bit float format % Special numbers dec2bin(hex2dec(num2hex(2^-1023)),64) % smallest positive normal double dec2bin(hex2dec(num2hex(2^-1074)),64) % smallest positive subnormal double num2hex(0) % all zeros num2hex(-0) % all zeros with negative sign bit num2hex(NaN) % one particular Not-a-Number num2hex(Inf) % special bitstring for positive infinity num2hex(-Inf) % special bitstring for negative infinity %%% Interpolations for a smooth function f = @(x) cos(x); % signal a=0; % Interval of interest is [a,b] b=15; % where a