/usr/share/yacas/numbers.rep/GaussianIntegers.ys is in yacas 1.3.1-1.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 | GaussianNorm(z_IsGaussianInteger) <-- Re(z)^2+Im(z)^2;
5 # IsGaussianInteger(x_IsList) <-- False;
// ?????? why is the following rule needed?
// 5 # IsGaussianInteger(ProductPrimesTo257) <-- False;
10 # IsGaussianInteger(x_IsComplex) <-- (IsInteger(Re(x)) And IsInteger(Im(x)));
// to catch IsGaussianInteger(x+2) from Apart
15 # IsGaussianInteger(_x) <-- False;
Function("IsGaussianPrime",{x})
[
if( IsGaussianInteger(x) )[
if( IsZero(Re(x)) )[
( Abs(Im(x)) % 4 = 3 And IsPrime(Abs(Im(x))) );
] else if ( IsZero(Im(x)) ) [
( Abs(Re(x)) % 4 = 3 And IsPrime(Abs(Re(x))) );
] else [
IsPrime(Re(x)^2 + Im(x)^2);
];
] else [
False;
];
];
/*
10 # IsGaussianPrime(p_IsInteger) <-- IsPrime(p) And Mod(p,3)=1;
20 # IsGaussianPrime(p_IsGaussianInteger) <-- IsPrime(GaussianNorm(p));
*/
GaussianMod(z_IsGaussianInteger,w_IsGaussianInteger) <-- z - w * Round(z/w);
10 # GaussianGcd(n_IsGaussianInteger,m_IsGaussianInteger) <--
[
If(N(Abs(m))=0,n, GaussianGcd(m,n - m*Round(n/m) ) );
];
IsGaussianUnit(z_IsGaussianInteger) <-- GaussianNorm(z)=1;
/* GaussianFactorPrime(p): auxiliary function for Gaussian factors.
If p is a rational prime of the form 4n+1, we find a factor of p in the
Gaussian Integers. We compute
a = (2n)!
By Wilson's theorem a^2 is -1 (mod p), it follows that
p| (a+I)(a-I)
in the Gaussian integers. The desired factor is then the Gaussian GCD of a+i
and p. Note: If the result is Complex(a,b), then p=a^2+b^2 */
GaussianFactorPrime(p_IsInteger) <-- [
Local(a,i);
a := 1;
For (i:=2,i<=(p-1)/2,i++) a := Mod(a*i,p);
GaussianGcd(a+I,p);
];
/* AddGaussianFactor: auxiliary function for Gaussian Factors.
L is a lists of factors of the Gaussian integer z and p is a Gaussian prime
that we want to add to the list. We first find the exponent e of p in the
decomposition of z (into Gaussian primes). If it is not zero, we add {p,e}
to the list */
AddGaussianFactor(L_IsList,z_IsGaussianInteger,p_IsGaussianInteger) <--
[
Local(e);
e :=0;
While (IsGaussianInteger(z:= z/p)) e++;
If (e != 0, DestructiveAppend(L,{p,e}));
];
/* GaussianFactors(n) : returns a list of factors of n, in a similar
way to Factors(n).
If n is a rational integer, we factor n in the Gaussian integers, by first
factor it in the rational integers, and after that factoring each of
its integer prime factors. */
10 # GaussianFactors(n_IsInteger) <--
[
Local(ifactors,gfactors,p,alpha);
ifactors := Factors(n);
If (Contains({{-1,1}, {1,1}}, Head(ifactors)), ifactors := Tail(ifactors) );
gfactors := {};
ForEach(p,ifactors)
[
If (p[1]=2, [ DestructiveAppend(gfactors,{1+I,p[2]});
DestructiveAppend(gfactors,{1-I,p[2]}); ]);
If (Mod(p[1],4)=3, DestructiveAppend(gfactors,p));
If (Mod(p[1],4)=1, [ alpha := GaussianFactorPrime(p[1]);
DestructiveAppend(gfactors,{alpha,p[2]});
DestructiveAppend(gfactors,{Conjugate(alpha),p[2]}); ]);
];
gfactors;
];
/* If z is is a Gaussian integer, we find its possible Gassian prime factors,
by factoring its norm */
20 # GaussianFactors(z_IsGaussianInteger) <--
[
Local(n,nfactors,gfactors,p);
gfactors :={};
n := GaussianNorm(z);
nfactors := Factors(n);
If (Contains({{-1,1}, {1,1}}, Head(nfactors)), nfactors := Tail(nfactors) );
ForEach(p,nfactors)
[
If (p[1]=2, [ AddGaussianFactor(gfactors,z,1+I);]);
If (Mod(p[1],4)=3, AddGaussianFactor(gfactors,z,p[1]));
If (Mod(p[1],4)=1, [ Local(alpha);
alpha := GaussianFactorPrime(p[1]);
AddGaussianFactor(gfactors,z,alpha);
AddGaussianFactor(gfactors,z,Conjugate(alpha));
]);
];
gfactors;
];
// Algorithm adapted from: Number Theory: A Programmer's Guide
// Mark Herkommer
// Program 8.7.1c, p 264
// This function needs to be modified to return the factors in
// data structure instead of printing them out
// THIS FUNCTION IS DEPRECATED NOW!
// Use GaussianFactors instead (Pablo)
// I've leave this here so that you can compare the eficiency of one
// function against the other
Function("FactorGaussianInteger",{x}) [
Check( IsGaussianInteger(x), "FactorGaussianInteger: argument must be a Gaussian integer");
Local(re,im,norm,a,b,d,i,j);
re:=Re(x);im:=Im(x);
If(re<0, re:=(-re) );
If(im<0, im:=(-im) );
norm:=re^2+im^2;
if( IsComposite(norm) )[
For(i:=0, i^2 <= norm, i++ )[ // real part
For(j:=0, i^2 + j^2 <= norm, j++)[ // complex part
if( Not( (i = re And j = im) Or
(i = im And j = re) ) )[ // no associates
d:=i^2+j^2;
if( d > 1 )[
a := re * i + im * j;
b := im * i - re * j;
While( (Mod(a,d) = 0) And (Mod(b,d) = 0) ) [
FactorGaussianInteger(Complex(i,j));
re:= a/d;
im:= b/d;
a := re * i + im * j;
b := im * i - re * j;
norm := re^2 + im^2;
];
];
];
];
];
If( re != 1 Or im != 0, Echo(Complex(re,im)) );
] else [
Echo(Complex(re,im));
];
];
|