RandomSecretQuaternionAlgebra := function();
  A := AssociativeAlgebra(QuaternionAlgebra<Rationals() | -1,1>);
  vecs := [&+[Random(10)*A.i : i in [1..4]] : j in [1..4]];
  Mchange := Matrix(Rationals(),4,4,&cat[Eltseq(vecs[i]) : i in [1..4]]);
  Mchange := Mchange^(-1);
  seq := [<i,j,k,((vecs[i]*vecs[j])*Mchange)[k]> : i,j,k in [1..4]];
  A := AssociativeAlgebra<Rationals(),4 | seq>;
  return A;
end function;

// QUATERNION ALGEBRAS
A<alpha,beta> := QuaternionAlgebra<Rationals() | 3,-5>;
alpha^2;
beta^2;
beta*alpha eq -alpha*beta;
theta := 1+alpha-2*beta+3*alpha*beta;
MinimalPolynomial(theta);
theta^2 - Trace(theta)*theta + Norm(theta);

// RECOGNIZING QUATERNION ALGEBRAS
A := RandomSecretQuaternionAlgebra();
A;
[ A.i*A.j : i,j in [1..4]];
bl, Aquat, phi := IsQuaternionAlgebra(A);
bl;
Aquat;
StandardForm(Aquat);
Aquat.1^2, Aquat.2^2;
phi;

// RECOGNIZING THE MATRIX RING
IsMatrixRing( QuaternionAlgebra<QuadraticField(-3) | -1,-1> );
IsMatrixRing( QuaternionAlgebra<QuadraticField(5) | -1,-1> );

A := QuaternionAlgebra<Rationals() | -1, 1>;
eps := A.3-1;
MinimalPolynomial(eps), Norm(eps);
M2F, phi := MatrixRing(A,eps);
phi;
[<MinimalPolynomial(A.i), MinimalPolynomial(phi(A.i))> : i in [1..3]];

// HILBERT SYMBOL
_<x> := PolynomialRing(Rationals());
K<s> := NumberField(x^2-6);
Z_K := MaximalOrder(K);
S := [x+y*Z_K.2 : x,y in [0..7] | x*y ne 0];
a := Random(S); b := Random(S);
A := QuaternionAlgebra<K | a,b>;
RamifiedPlaces(A);
[<pp, HilbertSymbol(A, pp[1])> : pp in 
   Factorization(ideal<Z_K | 2*Norm(a)*Norm(b)>)];

F<b> := NumberField(x^3-3*x-1);
Foo := InfinitePlaces(F);
Z_F := MaximalOrder(F);
I := ideal<Z_F | 6>;
A := QuaternionAlgebra(I);
FactoredDiscriminant(A);
A := QuaternionAlgebra(ideal<Z_F | 1>, Foo[1..2]);
FactoredDiscriminant(A);

// MAXIMAL ORDERS
A<alpha,beta,alphabeta> := QuaternionAlgebra<F | -3,b>;
O := MaximalOrder(A);
Factorization(Discriminant(O));
PseudoBasis(O);

A := QuaternionAlgebra<QuadraticField(1234830) | -123748, -12349813>;
O := MaximalOrder(A);
Discriminant(O);

// IDEAL CLASSES
pp := Decomposition(Z_F, 17)[1][1];
A := QuaternionAlgebra(pp, Foo);
O := MaximalOrder(A);
time Rideals := RightIdealClasses(O);
#Rideals;
Irand := RandomRightIdeal(O);
IsIsomorphic(Irand, Rideals[1]);
IsIsomorphic(Irand, Rideals[2]);

// UNIT GROUPS
A := QuaternionAlgebra(ideal<Z_F | 2>, Foo);
IsDefinite(A);
O := MaximalOrder(A);
U, h := UnitGroup(O);
U;
#Units(O);
